Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : #pragma once
5 :
6 : #include <array>
7 : #include <cstddef>
8 : #include <iterator>
9 : #include <limits>
10 : #include <type_traits>
11 : #include <utility>
12 :
13 : #include "DataStructures/DataBox/Prefixes.hpp"
14 : #include "DataStructures/DataVector.hpp"
15 : #include "DataStructures/Index.hpp"
16 : #include "DataStructures/Tensor/Tensor.hpp"
17 : #include "DataStructures/Variables.hpp"
18 : #include "Domain/Structure/Direction.hpp"
19 : #include "Domain/Structure/DirectionMap.hpp"
20 : #include "Domain/Structure/DirectionalIdMap.hpp"
21 : #include "Domain/Structure/ElementId.hpp"
22 : #include "Evolution/DgSubcell/GhostData.hpp"
23 : #include "NumericalAlgorithms/FiniteDifference/DerivativeOrder.hpp"
24 : #include "NumericalAlgorithms/Spectral/Mesh.hpp"
25 : #include "Utilities/Algorithm.hpp"
26 : #include "Utilities/ErrorHandling/Error.hpp"
27 : #include "Utilities/Gsl.hpp"
28 : #include "Utilities/OptionalHelpers.hpp"
29 : #include "Utilities/TMPL.hpp"
30 :
31 : namespace fd {
32 : /// @{
33 : /*!
34 : * \brief Computes a high-order boundary correction $G$ at the FD interface.
35 : *
36 : * The correction to the second-order boundary correction is given by
37 : *
38 : * \f{align*}{
39 : * G=G^{(2)}-G^{(4)}+G^{(6)}-G^{(8)}+G^{(10)},
40 : * \f}
41 : *
42 : * where
43 : *
44 : *\f{align*}{
45 : * G^{(4)}_{j+1/2}&=\frac{1}{6}\left(G_j -2 G^{(2)} +
46 : * G_{j+1}\right), \\
47 : * G^{(6)}_{j+1/2}&=\frac{1}{180}\left(G_{j-1} - 9 G_j + 16 G^{(2)}
48 : * -9 G_{j+1} + G_{j+2}\right), \\
49 : * G^{(8)}_{j+1/2}&=\frac{1}{2100}\left(G_{j-2} - \frac{25}{3} G_{j-1}
50 : * + 50 G_j - \frac{256}{3} G^{(2)} + 50 G_{j+1}
51 : * - \frac{25}{3} G_{j+2} +G_{j+3}\right), \\
52 : * G^{(10)}_{j+1/2}&=\frac{1}{17640}
53 : * \left(G_{j-3} - \frac{49}{5} G_{j-2}
54 : * + 49 G_{j-1} - 245 G_j + \frac{2048}{5} G^{(2)}\right.
55 : * \nonumber \\
56 : * &\left.- 245 G_{j+1}+ 49 G_{j+2} - \frac{49}{5} G_{j+3}
57 : * + G_{j+4}\right),
58 : * \f}
59 : *
60 : * where
61 : *
62 : * \f{align*}{
63 : * G_{j} &= F^i_j n_i^{j+1/2}, \\
64 : * G_{j\pm1} &= F^i_{j\pm1} n_i^{j+1/2}, \\
65 : * G_{j\pm2} &= F^i_{j\pm2} n_i^{j+1/2}, \\
66 : * G_{j\pm3} &= F^i_{j\pm3} n_i^{j+1/2}, \\
67 : * G_{j\pm4} &= F^i_{j\pm4} n_i^{j+1/2}.
68 : * \f}
69 : *
70 : * This is a generalization of the correction presented in \cite CHEN2016604.
71 : *
72 : * This high-order flux can be fed into a flux limiter, e.g. to guarantee
73 : * positivity.
74 : *
75 : * \note This implementation should be profiled and optimized.
76 : *
77 : * \warning This documentation is for the general case. In the restricted
78 : * Cartesian case we use the cell-centered flux as opposed to `G^{(4)}`, which
79 : * differs by a minus sign. This amounts to a minus sign change in front of the
80 : * $G^{(k)}$ terms in computing $G$ for $k>2$, and also a sign change in front
81 : * of $G^{(2)}$ in all $G^{(k)}$ for $k>2$.
82 : */
83 : template <DerivativeOrder DerivOrder, size_t Dim, typename... EvolvedVarsTags>
84 1 : void cartesian_high_order_fluxes_using_nodes(
85 : const gsl::not_null<
86 : std::array<Variables<tmpl::list<EvolvedVarsTags...>>, Dim>*>
87 : high_order_boundary_corrections_in_logical_direction,
88 :
89 : const std::array<Variables<tmpl::list<EvolvedVarsTags...>>, Dim>&
90 : second_order_boundary_corrections_in_logical_direction,
91 : const Variables<tmpl::list<
92 : ::Tags::Flux<EvolvedVarsTags, tmpl::size_t<Dim>, Frame::Inertial>...>>&
93 : cell_centered_inertial_flux,
94 : const DirectionMap<
95 : Dim, Variables<tmpl::list<::Tags::Flux<
96 : EvolvedVarsTags, tmpl::size_t<Dim>, Frame::Inertial>...>>>&
97 : ghost_cell_inertial_flux,
98 : const Mesh<Dim>& subcell_mesh, const size_t number_of_ghost_cells,
99 : [[maybe_unused]] const std::array<gsl::span<std::uint8_t>, Dim>&
100 : reconstruction_order = {}) {
101 : using std::min;
102 : constexpr int max_correction_order = 10;
103 : static_assert(static_cast<int>(DerivOrder) <= max_correction_order);
104 : constexpr size_t stencil_size = static_cast<int>(DerivOrder) < 0
105 : ? 8
106 : : (static_cast<size_t>(DerivOrder) - 2);
107 : const size_t correction_width =
108 : min(static_cast<size_t>(DerivOrder) / 2 - 1,
109 : min(number_of_ghost_cells, stencil_size / 2));
110 : ASSERT(correction_width <= number_of_ghost_cells,
111 : "The width of the derivative correction ("
112 : << correction_width
113 : << ") must be less than or equal to the number of ghost cells "
114 : << number_of_ghost_cells);
115 : ASSERT(alg::all_of(reconstruction_order,
116 : [](const auto& t) { return not t.empty(); }) or
117 : static_cast<int>(DerivOrder) > 0,
118 : "For adaptive derivative orders the reconstruction_order must be set");
119 : for (size_t dim = 0; dim < Dim; ++dim) {
120 : gsl::at(*high_order_boundary_corrections_in_logical_direction, dim)
121 : .initialize(
122 : gsl::at(second_order_boundary_corrections_in_logical_direction, dim)
123 : .number_of_grid_points());
124 : }
125 :
126 : // Reconstruction order is always first-varying fastest since we don't
127 : // transpose that back to {x,y,z} ordering.
128 : Index<Dim> reconstruction_extents = subcell_mesh.extents();
129 : reconstruction_extents[0] += 2;
130 :
131 : const auto impl = [&cell_centered_inertial_flux, &ghost_cell_inertial_flux,
132 : &high_order_boundary_corrections_in_logical_direction,
133 : number_of_ghost_cells,
134 : &second_order_boundary_corrections_in_logical_direction,
135 : &subcell_mesh, &correction_width, &reconstruction_order,
136 : &reconstruction_extents](auto tag_v, auto dim_v) {
137 : (void)reconstruction_extents;
138 : using tag = decltype(tag_v);
139 : constexpr size_t dim = decltype(dim_v)::value;
140 :
141 : auto& high_order_var_correction =
142 : get<tag>((*high_order_boundary_corrections_in_logical_direction)[dim]);
143 : const auto& second_order_var_correction =
144 : get<tag>(second_order_boundary_corrections_in_logical_direction[dim]);
145 : const auto& recons_order = reconstruction_order[dim];
146 : const auto& cell_centered_flux =
147 : get<::Tags::Flux<tag, tmpl::size_t<Dim>, Frame::Inertial>>(
148 : cell_centered_inertial_flux);
149 : const auto& lower_neighbor_cell_centered_flux =
150 : get<::Tags::Flux<tag, tmpl::size_t<Dim>, Frame::Inertial>>(
151 : ghost_cell_inertial_flux.at(Direction<Dim>{dim, Side::Lower}));
152 : const auto& upper_neighbor_cell_centered_flux =
153 : get<::Tags::Flux<tag, tmpl::size_t<Dim>, Frame::Inertial>>(
154 : ghost_cell_inertial_flux.at(Direction<Dim>{dim, Side::Upper}));
155 : using FluxTensor = std::decay_t<decltype(cell_centered_flux)>;
156 : const auto& subcell_extents = subcell_mesh.extents();
157 : auto subcell_face_extents = subcell_extents;
158 : ++subcell_face_extents[dim];
159 : auto neighbor_extents = subcell_extents;
160 : neighbor_extents[dim] = number_of_ghost_cells;
161 : const size_t number_of_components = second_order_var_correction.size();
162 : for (size_t storage_index = 0; storage_index < number_of_components;
163 : ++storage_index) {
164 : const auto flux_multi_index = prepend(
165 : second_order_var_correction.get_tensor_index(storage_index), dim);
166 : const size_t flux_storage_index =
167 : FluxTensor::get_storage_index(flux_multi_index);
168 : // Loop over each face
169 : for (size_t k = 0; k < (Dim == 3 ? subcell_face_extents[2] : 1); ++k) {
170 : for (size_t j = 0; j < (Dim >= 2 ? subcell_face_extents[1] : 1); ++j) {
171 : for (size_t i = 0; i < subcell_face_extents[0]; ++i) {
172 : const Index<Dim> face_index = [i, j, k]() -> Index<Dim> {
173 : if constexpr (Dim == 3) {
174 : return Index<Dim>{i, j, k};
175 : } else if constexpr (Dim == 2) {
176 : (void)k;
177 : return Index<Dim>{i, j};
178 : } else {
179 : (void)k, (void)j;
180 : return Index<Dim>{i};
181 : }
182 : }();
183 : const size_t face_storage_index =
184 : collapsed_index(face_index, subcell_face_extents);
185 : Index<Dim> neighbor_index{};
186 : for (size_t l = 0; l < Dim; ++l) {
187 : if (l != dim) {
188 : neighbor_index[l] = face_index[l];
189 : }
190 : }
191 :
192 : double& correction =
193 : high_order_var_correction[storage_index][face_storage_index] =
194 : 0.0;
195 :
196 : std::array<double, stencil_size> cell_centered_fluxes_for_stencil{};
197 : // fill if we have to retrieve from lower neighbor
198 : size_t stencil_index = 0;
199 : for (int grid_index = static_cast<int>(face_index[dim]) -
200 : static_cast<int>(correction_width);
201 : grid_index < static_cast<int>(face_index[dim]) +
202 : static_cast<int>(correction_width);
203 : ++grid_index, ++stencil_index) {
204 : if (grid_index < 0) {
205 : neighbor_index[dim] = static_cast<size_t>(
206 : static_cast<int>(number_of_ghost_cells) + grid_index);
207 : gsl::at(cell_centered_fluxes_for_stencil, stencil_index) =
208 : lower_neighbor_cell_centered_flux[flux_storage_index]
209 : [collapsed_index(
210 : neighbor_index,
211 : neighbor_extents)];
212 : } else if (grid_index >= static_cast<int>(subcell_extents[dim])) {
213 : neighbor_index[dim] = static_cast<size_t>(
214 : grid_index - static_cast<int>(subcell_extents[dim]));
215 : gsl::at(cell_centered_fluxes_for_stencil, stencil_index) =
216 : upper_neighbor_cell_centered_flux[flux_storage_index]
217 : [collapsed_index(
218 : neighbor_index,
219 : neighbor_extents)];
220 : } else {
221 : Index<Dim> volume_index = face_index;
222 : volume_index[dim] = static_cast<size_t>(grid_index);
223 : gsl::at(cell_centered_fluxes_for_stencil, stencil_index) =
224 : cell_centered_flux[flux_storage_index][collapsed_index(
225 : volume_index, subcell_extents)];
226 : }
227 : }
228 :
229 : size_t lower_neighbor_index = std::numeric_limits<size_t>::max();
230 : size_t upper_neighbor_index = std::numeric_limits<size_t>::max();
231 : if constexpr (static_cast<int>(DerivOrder) < 0) {
232 : Index<Dim> lower_n{};
233 : Index<Dim> upper_n{};
234 : if constexpr (dim == 0) {
235 : if constexpr (Dim == 1) {
236 : lower_n = Index<Dim>{i};
237 : upper_n = Index<Dim>{i + 1};
238 : } else if constexpr (Dim == 2) {
239 : lower_n = Index<Dim>{i, j};
240 : upper_n = Index<Dim>{i + 1, j};
241 : } else if constexpr (Dim == 3) {
242 : lower_n = Index<Dim>{i, j, k};
243 : upper_n = Index<Dim>{i + 1, j, k};
244 : }
245 : } else if constexpr (dim == 1) {
246 : if constexpr (Dim == 2) {
247 : lower_n = Index<Dim>{j, i};
248 : upper_n = Index<Dim>{j + 1, i};
249 : } else if constexpr (Dim == 3) {
250 : lower_n = Index<Dim>{j, k, i};
251 : upper_n = Index<Dim>{j + 1, k, i};
252 : }
253 : } else if constexpr (dim == 2) {
254 : lower_n = Index<Dim>{k, i, j};
255 : upper_n = Index<Dim>{k + 1, i, j};
256 : }
257 : lower_neighbor_index =
258 : collapsed_index(lower_n, reconstruction_extents);
259 : upper_neighbor_index =
260 : collapsed_index(upper_n, reconstruction_extents);
261 : }
262 :
263 : if (static_cast<int>(DerivOrder) >= 10 or
264 : (static_cast<int>(DerivOrder) < 0 and
265 : min(recons_order[lower_neighbor_index],
266 : recons_order[upper_neighbor_index]) >= 9)) {
267 : correction -=
268 : 5.6689342403628117913e-5 *
269 : (gsl::at(cell_centered_fluxes_for_stencil,
270 : correction_width - 4) +
271 : gsl::at(cell_centered_fluxes_for_stencil,
272 : correction_width + 3) -
273 : 9.8 * (gsl::at(cell_centered_fluxes_for_stencil,
274 : correction_width - 3) +
275 : gsl::at(cell_centered_fluxes_for_stencil,
276 : correction_width + 2)) +
277 : 49.0 * (gsl::at(cell_centered_fluxes_for_stencil,
278 : correction_width - 2) +
279 : gsl::at(cell_centered_fluxes_for_stencil,
280 : correction_width + 1)) -
281 : 245.0 * (gsl::at(cell_centered_fluxes_for_stencil,
282 : correction_width - 1) +
283 : gsl::at(cell_centered_fluxes_for_stencil,
284 : correction_width)) -
285 : 409.6 * second_order_var_correction[storage_index]
286 : [face_storage_index]);
287 : }
288 : if (static_cast<int>(DerivOrder) >= 8 or
289 : (static_cast<int>(DerivOrder) < 0 and
290 : min(recons_order[lower_neighbor_index],
291 : recons_order[upper_neighbor_index]) >= 7)) {
292 : correction +=
293 : 4.7619047619047619047e-4 *
294 : (gsl::at(cell_centered_fluxes_for_stencil,
295 : correction_width - 3) +
296 : gsl::at(cell_centered_fluxes_for_stencil,
297 : correction_width + 2) -
298 : 8.3333333333333333333 *
299 : (gsl::at(cell_centered_fluxes_for_stencil,
300 : correction_width - 2) +
301 : gsl::at(cell_centered_fluxes_for_stencil,
302 : correction_width + 1)) +
303 : 50.0 * (gsl::at(cell_centered_fluxes_for_stencil,
304 : correction_width - 1) +
305 : gsl::at(cell_centered_fluxes_for_stencil,
306 : correction_width)) +
307 : 85.333333333333333333 *
308 : second_order_var_correction[storage_index]
309 : [face_storage_index]);
310 : }
311 : if (static_cast<int>(DerivOrder) >= 6 or
312 : (static_cast<int>(DerivOrder) < 0 and
313 : min(recons_order[lower_neighbor_index],
314 : recons_order[upper_neighbor_index]) >=
315 : (DerivOrder ==
316 : DerivativeOrder::OneHigherThanReconsButFiveToFour
317 : ? 6
318 : : 5))) {
319 : correction -=
320 : 5.5555555555555555555e-3 *
321 : (gsl::at(cell_centered_fluxes_for_stencil,
322 : correction_width - 2) +
323 : gsl::at(cell_centered_fluxes_for_stencil,
324 : correction_width + 1) -
325 : 9.0 * (gsl::at(cell_centered_fluxes_for_stencil,
326 : correction_width - 1) +
327 : gsl::at(cell_centered_fluxes_for_stencil,
328 : correction_width)) -
329 : 16.0 * second_order_var_correction[storage_index]
330 : [face_storage_index]);
331 : }
332 : if (static_cast<int>(DerivOrder) >= 4 or
333 : (static_cast<int>(DerivOrder) < 0 and
334 : min(recons_order[lower_neighbor_index],
335 : recons_order[upper_neighbor_index]) >= 3)) {
336 : correction +=
337 : 0.166666666666666666 *
338 : (gsl::at(cell_centered_fluxes_for_stencil,
339 : correction_width - 1) +
340 : gsl::at(cell_centered_fluxes_for_stencil, correction_width) +
341 : 2.0 * second_order_var_correction[storage_index]
342 : [face_storage_index]);
343 : }
344 :
345 : // Add second-order correction last
346 : correction +=
347 : second_order_var_correction[storage_index][face_storage_index];
348 : }
349 : }
350 : }
351 : }
352 : };
353 :
354 : EXPAND_PACK_LEFT_TO_RIGHT(
355 : impl(EvolvedVarsTags{}, std::integral_constant<size_t, 0>{}));
356 : if constexpr (Dim > 1) {
357 : EXPAND_PACK_LEFT_TO_RIGHT(
358 : impl(EvolvedVarsTags{}, std::integral_constant<size_t, 1>{}));
359 : if constexpr (Dim > 2) {
360 : EXPAND_PACK_LEFT_TO_RIGHT(
361 : impl(EvolvedVarsTags{}, std::integral_constant<size_t, 2>{}));
362 : }
363 : }
364 : }
365 :
366 : template <size_t Dim, typename... EvolvedVarsTags>
367 1 : void cartesian_high_order_fluxes_using_nodes(
368 : const gsl::not_null<
369 : std::array<Variables<tmpl::list<EvolvedVarsTags...>>, Dim>*>
370 : high_order_boundary_corrections_in_logical_direction,
371 :
372 : const std::array<Variables<tmpl::list<EvolvedVarsTags...>>, Dim>&
373 : second_order_boundary_corrections_in_logical_direction,
374 : const Variables<tmpl::list<
375 : ::Tags::Flux<EvolvedVarsTags, tmpl::size_t<Dim>, Frame::Inertial>...>>&
376 : cell_centered_inertial_flux,
377 : const DirectionMap<
378 : Dim, Variables<tmpl::list<::Tags::Flux<
379 : EvolvedVarsTags, tmpl::size_t<Dim>, Frame::Inertial>...>>>&
380 : ghost_cell_inertial_flux,
381 : const Mesh<Dim>& subcell_mesh, const size_t number_of_ghost_cells,
382 : const DerivativeOrder derivative_order,
383 : [[maybe_unused]] const std::array<gsl::span<std::uint8_t>, Dim>&
384 : reconstruction_order = {}) {
385 : switch (derivative_order) {
386 : case DerivativeOrder::OneHigherThanRecons:
387 : cartesian_high_order_fluxes_using_nodes<
388 : DerivativeOrder::OneHigherThanRecons>(
389 : high_order_boundary_corrections_in_logical_direction,
390 : second_order_boundary_corrections_in_logical_direction,
391 : cell_centered_inertial_flux, ghost_cell_inertial_flux, subcell_mesh,
392 : number_of_ghost_cells, reconstruction_order);
393 : break;
394 : case DerivativeOrder::OneHigherThanReconsButFiveToFour:
395 : cartesian_high_order_fluxes_using_nodes<
396 : DerivativeOrder::OneHigherThanReconsButFiveToFour>(
397 : high_order_boundary_corrections_in_logical_direction,
398 : second_order_boundary_corrections_in_logical_direction,
399 : cell_centered_inertial_flux, ghost_cell_inertial_flux, subcell_mesh,
400 : number_of_ghost_cells, reconstruction_order);
401 : break;
402 : case DerivativeOrder::Two:
403 : cartesian_high_order_fluxes_using_nodes<DerivativeOrder::Two>(
404 : high_order_boundary_corrections_in_logical_direction,
405 : second_order_boundary_corrections_in_logical_direction,
406 : cell_centered_inertial_flux, ghost_cell_inertial_flux, subcell_mesh,
407 : number_of_ghost_cells, reconstruction_order);
408 : break;
409 : case DerivativeOrder::Four:
410 : cartesian_high_order_fluxes_using_nodes<DerivativeOrder::Four>(
411 : high_order_boundary_corrections_in_logical_direction,
412 : second_order_boundary_corrections_in_logical_direction,
413 : cell_centered_inertial_flux, ghost_cell_inertial_flux, subcell_mesh,
414 : number_of_ghost_cells, reconstruction_order);
415 : break;
416 : case DerivativeOrder::Six:
417 : cartesian_high_order_fluxes_using_nodes<DerivativeOrder::Six>(
418 : high_order_boundary_corrections_in_logical_direction,
419 : second_order_boundary_corrections_in_logical_direction,
420 : cell_centered_inertial_flux, ghost_cell_inertial_flux, subcell_mesh,
421 : number_of_ghost_cells, reconstruction_order);
422 : break;
423 : case DerivativeOrder::Eight:
424 : cartesian_high_order_fluxes_using_nodes<DerivativeOrder::Eight>(
425 : high_order_boundary_corrections_in_logical_direction,
426 : second_order_boundary_corrections_in_logical_direction,
427 : cell_centered_inertial_flux, ghost_cell_inertial_flux, subcell_mesh,
428 : number_of_ghost_cells, reconstruction_order);
429 : break;
430 : case DerivativeOrder::Ten:
431 : cartesian_high_order_fluxes_using_nodes<DerivativeOrder::Ten>(
432 : high_order_boundary_corrections_in_logical_direction,
433 : second_order_boundary_corrections_in_logical_direction,
434 : cell_centered_inertial_flux, ghost_cell_inertial_flux, subcell_mesh,
435 : number_of_ghost_cells, reconstruction_order);
436 : break;
437 : default:
438 : ERROR("Unsupported correction order " << derivative_order);
439 : };
440 : }
441 : /// @}
442 :
443 : /*!
444 : * \brief Fill the `flux_neighbor_data` with pointers into the
445 : * `all_ghost_data`.
446 : *
447 : * The `all_ghost_data` is stored in the tag
448 : * `evolution::dg::subcell::Tags::GhostDataForReconstruction`, and the
449 : * `ghost_zone_size` should come from the FD reconstructor.
450 : */
451 : template <size_t Dim, typename FluxesTags>
452 1 : void set_cartesian_neighbor_cell_centered_fluxes(
453 : const gsl::not_null<DirectionMap<Dim, Variables<FluxesTags>>*>
454 : flux_neighbor_data,
455 : const DirectionalIdMap<Dim, evolution::dg::subcell::GhostData>&
456 : all_ghost_data,
457 : const Mesh<Dim>& subcell_mesh, const size_t ghost_zone_size,
458 : const size_t number_of_rdmp_values_in_ghost_data) {
459 : for (const auto& [direction_id, ghost_data] : all_ghost_data) {
460 : const size_t neighbor_flux_size =
461 : subcell_mesh.number_of_grid_points() /
462 : subcell_mesh.extents(direction_id.direction.dimension()) *
463 : ghost_zone_size *
464 : Variables<FluxesTags>::number_of_independent_components;
465 : const DataVector& neighbor_data =
466 : ghost_data.neighbor_ghost_data_for_reconstruction();
467 : (*flux_neighbor_data)[direction_id.direction].set_data_ref(
468 : // NOLINTNEXTLINE(cppcoreguidelines-pro-type-const-cast)
469 : const_cast<double*>(std::next(
470 : neighbor_data.data(),
471 : static_cast<std::ptrdiff_t>(neighbor_data.size() -
472 : number_of_rdmp_values_in_ghost_data -
473 : neighbor_flux_size))),
474 : neighbor_flux_size);
475 : }
476 : }
477 :
478 : /*!
479 : * \brief Computes the high-order Cartesian flux corrections if necessary.
480 : *
481 : * The `cell_centered_fluxes` is stored in the tag
482 : * `evolution::dg::subcell::Tags::CellCenteredFlux`, `fd_derivative_order` is
483 : * from `evolution::dg::subcell::Tags::SubcellOptions`
484 : * (`.finite_difference_derivative_order()`), the `all_ghost_data`
485 : * is stored in the tag
486 : * `evolution::dg::subcell::Tags::GhostDataForReconstruction`, the
487 : * `ghost_zone_size` should come from the FD reconstructor.
488 : *
489 : * By default we assume no RDMP data is in the `ghost_data` buffer. In the
490 : * future we will want to update how we store the data in order to eliminate
491 : * more memory allocations and copies, in which case that value will be
492 : * non-zero.
493 : *
494 : * \note `high_order_corrections` must either not have a value or have all
495 : * elements be of the same size as
496 : * `second_order_boundary_corrections[0].number_of_grid_points()`, where we've
497 : * assumed `second_order_boundary_corrections` is the same in all directions.
498 : */
499 : template <size_t Dim, typename... EvolvedVarsTags,
500 : typename FluxesTags = tmpl::list<::Tags::Flux<
501 : EvolvedVarsTags, tmpl::size_t<Dim>, Frame::Inertial>...>>
502 1 : void cartesian_high_order_flux_corrections(
503 : const gsl::not_null<std::optional<
504 : std::array<Variables<tmpl::list<EvolvedVarsTags...>>, Dim>>*>
505 : high_order_corrections,
506 :
507 : const std::optional<Variables<FluxesTags>>& cell_centered_fluxes,
508 : const std::array<Variables<tmpl::list<EvolvedVarsTags...>>, Dim>&
509 : second_order_boundary_corrections,
510 : const fd::DerivativeOrder& fd_derivative_order,
511 : const DirectionalIdMap<Dim, evolution::dg::subcell::GhostData>&
512 : all_ghost_data,
513 : const Mesh<Dim>& subcell_mesh, const size_t ghost_zone_size,
514 : [[maybe_unused]] const std::array<gsl::span<std::uint8_t>, Dim>&
515 : reconstruction_order = {},
516 : const size_t number_of_rdmp_values_in_ghost_data = 0) {
517 : if (cell_centered_fluxes.has_value()) {
518 : ASSERT(alg::all_of(
519 : second_order_boundary_corrections,
520 : [expected_size = second_order_boundary_corrections[0]
521 : .number_of_grid_points()](const auto& e) {
522 : return e.number_of_grid_points() == expected_size;
523 : }),
524 : "All second-order boundary corrections must be of the same size, "
525 : << second_order_boundary_corrections[0].number_of_grid_points());
526 : if (fd_derivative_order != DerivativeOrder::Two) {
527 : if (not high_order_corrections->has_value()) {
528 : (*high_order_corrections) =
529 : make_array<Dim>(Variables<tmpl::list<EvolvedVarsTags...>>{
530 : second_order_boundary_corrections[0].number_of_grid_points()});
531 : }
532 : ASSERT(
533 : high_order_corrections->has_value() and
534 : alg::all_of(high_order_corrections->value(),
535 : [expected_size =
536 : second_order_boundary_corrections[0]
537 : .number_of_grid_points()](const auto& e) {
538 : return e.number_of_grid_points() == expected_size;
539 : }),
540 : "The high_order_corrections must all have size "
541 : << second_order_boundary_corrections[0].number_of_grid_points());
542 : DirectionMap<Dim, Variables<FluxesTags>> flux_neighbor_data{};
543 : set_cartesian_neighbor_cell_centered_fluxes(
544 : make_not_null(&flux_neighbor_data), all_ghost_data, subcell_mesh,
545 : ghost_zone_size, number_of_rdmp_values_in_ghost_data);
546 :
547 : cartesian_high_order_fluxes_using_nodes(
548 : make_not_null(&(high_order_corrections->value())),
549 : second_order_boundary_corrections, cell_centered_fluxes.value(),
550 : flux_neighbor_data, subcell_mesh, ghost_zone_size,
551 : fd_derivative_order, reconstruction_order);
552 : }
553 : }
554 : }
555 : } // namespace fd
|