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 : if constexpr (Dim == 3) {
255 : lower_n = Index<Dim>{k, i, j};
256 : upper_n = Index<Dim>{k + 1, i, j};
257 : }
258 : }
259 : lower_neighbor_index =
260 : collapsed_index(lower_n, reconstruction_extents);
261 : upper_neighbor_index =
262 : collapsed_index(upper_n, reconstruction_extents);
263 : }
264 :
265 : if (static_cast<int>(DerivOrder) >= 10 or
266 : (static_cast<int>(DerivOrder) < 0 and
267 : min(recons_order[lower_neighbor_index],
268 : recons_order[upper_neighbor_index]) >= 9)) {
269 : correction -=
270 : 5.6689342403628117913e-5 *
271 : (gsl::at(cell_centered_fluxes_for_stencil,
272 : correction_width - 4) +
273 : gsl::at(cell_centered_fluxes_for_stencil,
274 : correction_width + 3) -
275 : 9.8 * (gsl::at(cell_centered_fluxes_for_stencil,
276 : correction_width - 3) +
277 : gsl::at(cell_centered_fluxes_for_stencil,
278 : correction_width + 2)) +
279 : 49.0 * (gsl::at(cell_centered_fluxes_for_stencil,
280 : correction_width - 2) +
281 : gsl::at(cell_centered_fluxes_for_stencil,
282 : correction_width + 1)) -
283 : 245.0 * (gsl::at(cell_centered_fluxes_for_stencil,
284 : correction_width - 1) +
285 : gsl::at(cell_centered_fluxes_for_stencil,
286 : correction_width)) -
287 : 409.6 * second_order_var_correction[storage_index]
288 : [face_storage_index]);
289 : }
290 : if (static_cast<int>(DerivOrder) >= 8 or
291 : (static_cast<int>(DerivOrder) < 0 and
292 : min(recons_order[lower_neighbor_index],
293 : recons_order[upper_neighbor_index]) >= 7)) {
294 : correction +=
295 : 4.7619047619047619047e-4 *
296 : (gsl::at(cell_centered_fluxes_for_stencil,
297 : correction_width - 3) +
298 : gsl::at(cell_centered_fluxes_for_stencil,
299 : correction_width + 2) -
300 : 8.3333333333333333333 *
301 : (gsl::at(cell_centered_fluxes_for_stencil,
302 : correction_width - 2) +
303 : gsl::at(cell_centered_fluxes_for_stencil,
304 : correction_width + 1)) +
305 : 50.0 * (gsl::at(cell_centered_fluxes_for_stencil,
306 : correction_width - 1) +
307 : gsl::at(cell_centered_fluxes_for_stencil,
308 : correction_width)) +
309 : 85.333333333333333333 *
310 : second_order_var_correction[storage_index]
311 : [face_storage_index]);
312 : }
313 : if (static_cast<int>(DerivOrder) >= 6 or
314 : (static_cast<int>(DerivOrder) < 0 and
315 : min(recons_order[lower_neighbor_index],
316 : recons_order[upper_neighbor_index]) >=
317 : (DerivOrder ==
318 : DerivativeOrder::OneHigherThanReconsButFiveToFour
319 : ? 6
320 : : 5))) {
321 : correction -=
322 : 5.5555555555555555555e-3 *
323 : (gsl::at(cell_centered_fluxes_for_stencil,
324 : correction_width - 2) +
325 : gsl::at(cell_centered_fluxes_for_stencil,
326 : correction_width + 1) -
327 : 9.0 * (gsl::at(cell_centered_fluxes_for_stencil,
328 : correction_width - 1) +
329 : gsl::at(cell_centered_fluxes_for_stencil,
330 : correction_width)) -
331 : 16.0 * second_order_var_correction[storage_index]
332 : [face_storage_index]);
333 : }
334 : if (static_cast<int>(DerivOrder) >= 4 or
335 : (static_cast<int>(DerivOrder) < 0 and
336 : min(recons_order[lower_neighbor_index],
337 : recons_order[upper_neighbor_index]) >= 3)) {
338 : correction +=
339 : 0.166666666666666666 *
340 : (gsl::at(cell_centered_fluxes_for_stencil,
341 : correction_width - 1) +
342 : gsl::at(cell_centered_fluxes_for_stencil, correction_width) +
343 : 2.0 * second_order_var_correction[storage_index]
344 : [face_storage_index]);
345 : }
346 :
347 : // Add second-order correction last
348 : correction +=
349 : second_order_var_correction[storage_index][face_storage_index];
350 : }
351 : }
352 : }
353 : }
354 : };
355 :
356 : EXPAND_PACK_LEFT_TO_RIGHT(
357 : impl(EvolvedVarsTags{}, std::integral_constant<size_t, 0>{}));
358 : if constexpr (Dim > 1) {
359 : EXPAND_PACK_LEFT_TO_RIGHT(
360 : impl(EvolvedVarsTags{}, std::integral_constant<size_t, 1>{}));
361 : if constexpr (Dim > 2) {
362 : EXPAND_PACK_LEFT_TO_RIGHT(
363 : impl(EvolvedVarsTags{}, std::integral_constant<size_t, 2>{}));
364 : }
365 : }
366 : }
367 :
368 : template <size_t Dim, typename... EvolvedVarsTags>
369 1 : void cartesian_high_order_fluxes_using_nodes(
370 : const gsl::not_null<
371 : std::array<Variables<tmpl::list<EvolvedVarsTags...>>, Dim>*>
372 : high_order_boundary_corrections_in_logical_direction,
373 :
374 : const std::array<Variables<tmpl::list<EvolvedVarsTags...>>, Dim>&
375 : second_order_boundary_corrections_in_logical_direction,
376 : const Variables<tmpl::list<
377 : ::Tags::Flux<EvolvedVarsTags, tmpl::size_t<Dim>, Frame::Inertial>...>>&
378 : cell_centered_inertial_flux,
379 : const DirectionMap<
380 : Dim, Variables<tmpl::list<::Tags::Flux<
381 : EvolvedVarsTags, tmpl::size_t<Dim>, Frame::Inertial>...>>>&
382 : ghost_cell_inertial_flux,
383 : const Mesh<Dim>& subcell_mesh, const size_t number_of_ghost_cells,
384 : const DerivativeOrder derivative_order,
385 : [[maybe_unused]] const std::array<gsl::span<std::uint8_t>, Dim>&
386 : reconstruction_order = {}) {
387 : switch (derivative_order) {
388 : case DerivativeOrder::OneHigherThanRecons:
389 : cartesian_high_order_fluxes_using_nodes<
390 : DerivativeOrder::OneHigherThanRecons>(
391 : high_order_boundary_corrections_in_logical_direction,
392 : second_order_boundary_corrections_in_logical_direction,
393 : cell_centered_inertial_flux, ghost_cell_inertial_flux, subcell_mesh,
394 : number_of_ghost_cells, reconstruction_order);
395 : break;
396 : case DerivativeOrder::OneHigherThanReconsButFiveToFour:
397 : cartesian_high_order_fluxes_using_nodes<
398 : DerivativeOrder::OneHigherThanReconsButFiveToFour>(
399 : high_order_boundary_corrections_in_logical_direction,
400 : second_order_boundary_corrections_in_logical_direction,
401 : cell_centered_inertial_flux, ghost_cell_inertial_flux, subcell_mesh,
402 : number_of_ghost_cells, reconstruction_order);
403 : break;
404 : case DerivativeOrder::Two:
405 : cartesian_high_order_fluxes_using_nodes<DerivativeOrder::Two>(
406 : high_order_boundary_corrections_in_logical_direction,
407 : second_order_boundary_corrections_in_logical_direction,
408 : cell_centered_inertial_flux, ghost_cell_inertial_flux, subcell_mesh,
409 : number_of_ghost_cells, reconstruction_order);
410 : break;
411 : case DerivativeOrder::Four:
412 : cartesian_high_order_fluxes_using_nodes<DerivativeOrder::Four>(
413 : high_order_boundary_corrections_in_logical_direction,
414 : second_order_boundary_corrections_in_logical_direction,
415 : cell_centered_inertial_flux, ghost_cell_inertial_flux, subcell_mesh,
416 : number_of_ghost_cells, reconstruction_order);
417 : break;
418 : case DerivativeOrder::Six:
419 : cartesian_high_order_fluxes_using_nodes<DerivativeOrder::Six>(
420 : high_order_boundary_corrections_in_logical_direction,
421 : second_order_boundary_corrections_in_logical_direction,
422 : cell_centered_inertial_flux, ghost_cell_inertial_flux, subcell_mesh,
423 : number_of_ghost_cells, reconstruction_order);
424 : break;
425 : case DerivativeOrder::Eight:
426 : cartesian_high_order_fluxes_using_nodes<DerivativeOrder::Eight>(
427 : high_order_boundary_corrections_in_logical_direction,
428 : second_order_boundary_corrections_in_logical_direction,
429 : cell_centered_inertial_flux, ghost_cell_inertial_flux, subcell_mesh,
430 : number_of_ghost_cells, reconstruction_order);
431 : break;
432 : case DerivativeOrder::Ten:
433 : cartesian_high_order_fluxes_using_nodes<DerivativeOrder::Ten>(
434 : high_order_boundary_corrections_in_logical_direction,
435 : second_order_boundary_corrections_in_logical_direction,
436 : cell_centered_inertial_flux, ghost_cell_inertial_flux, subcell_mesh,
437 : number_of_ghost_cells, reconstruction_order);
438 : break;
439 : default:
440 : ERROR("Unsupported correction order " << derivative_order);
441 : };
442 : }
443 : /// @}
444 :
445 : /*!
446 : * \brief Fill the `flux_neighbor_data` with pointers into the
447 : * `all_ghost_data`.
448 : *
449 : * The `all_ghost_data` is stored in the tag
450 : * `evolution::dg::subcell::Tags::GhostDataForReconstruction`, and the
451 : * `ghost_zone_size` should come from the FD reconstructor.
452 : */
453 : template <size_t Dim, typename FluxesTags>
454 1 : void set_cartesian_neighbor_cell_centered_fluxes(
455 : const gsl::not_null<DirectionMap<Dim, Variables<FluxesTags>>*>
456 : flux_neighbor_data,
457 : const DirectionalIdMap<Dim, evolution::dg::subcell::GhostData>&
458 : all_ghost_data,
459 : const Mesh<Dim>& subcell_mesh, const size_t ghost_zone_size,
460 : const size_t number_of_rdmp_values_in_ghost_data) {
461 : for (const auto& [direction_id, ghost_data] : all_ghost_data) {
462 : const size_t neighbor_flux_size =
463 : subcell_mesh.number_of_grid_points() /
464 : subcell_mesh.extents(direction_id.direction().dimension()) *
465 : ghost_zone_size *
466 : Variables<FluxesTags>::number_of_independent_components;
467 : const DataVector& neighbor_data =
468 : ghost_data.neighbor_ghost_data_for_reconstruction();
469 : (*flux_neighbor_data)[direction_id.direction()].set_data_ref(
470 : // NOLINTNEXTLINE(cppcoreguidelines-pro-type-const-cast)
471 : const_cast<double*>(std::next(
472 : neighbor_data.data(),
473 : static_cast<std::ptrdiff_t>(neighbor_data.size() -
474 : number_of_rdmp_values_in_ghost_data -
475 : neighbor_flux_size))),
476 : neighbor_flux_size);
477 : }
478 : }
479 :
480 : /*!
481 : * \brief Computes the high-order Cartesian flux corrections if necessary.
482 : *
483 : * The `cell_centered_fluxes` is stored in the tag
484 : * `evolution::dg::subcell::Tags::CellCenteredFlux`, `fd_derivative_order` is
485 : * from `evolution::dg::subcell::Tags::SubcellOptions`
486 : * (`.finite_difference_derivative_order()`), the `all_ghost_data`
487 : * is stored in the tag
488 : * `evolution::dg::subcell::Tags::GhostDataForReconstruction`, the
489 : * `ghost_zone_size` should come from the FD reconstructor.
490 : *
491 : * By default we assume no RDMP data is in the `ghost_data` buffer. In the
492 : * future we will want to update how we store the data in order to eliminate
493 : * more memory allocations and copies, in which case that value will be
494 : * non-zero.
495 : *
496 : * \note `high_order_corrections` must either not have a value or have all
497 : * elements be of the same size as
498 : * `second_order_boundary_corrections[0].number_of_grid_points()`, where we've
499 : * assumed `second_order_boundary_corrections` is the same in all directions.
500 : */
501 : template <size_t Dim, typename... EvolvedVarsTags,
502 : typename FluxesTags = tmpl::list<::Tags::Flux<
503 : EvolvedVarsTags, tmpl::size_t<Dim>, Frame::Inertial>...>>
504 1 : void cartesian_high_order_flux_corrections(
505 : const gsl::not_null<std::optional<
506 : std::array<Variables<tmpl::list<EvolvedVarsTags...>>, Dim>>*>
507 : high_order_corrections,
508 :
509 : const std::optional<Variables<FluxesTags>>& cell_centered_fluxes,
510 : const std::array<Variables<tmpl::list<EvolvedVarsTags...>>, Dim>&
511 : second_order_boundary_corrections,
512 : const fd::DerivativeOrder& fd_derivative_order,
513 : const DirectionalIdMap<Dim, evolution::dg::subcell::GhostData>&
514 : all_ghost_data,
515 : const Mesh<Dim>& subcell_mesh, const size_t ghost_zone_size,
516 : [[maybe_unused]] const std::array<gsl::span<std::uint8_t>, Dim>&
517 : reconstruction_order = {},
518 : const size_t number_of_rdmp_values_in_ghost_data = 0) {
519 : if (cell_centered_fluxes.has_value()) {
520 : ASSERT(alg::all_of(
521 : second_order_boundary_corrections,
522 : [expected_size = second_order_boundary_corrections[0]
523 : .number_of_grid_points()](const auto& e) {
524 : return e.number_of_grid_points() == expected_size;
525 : }),
526 : "All second-order boundary corrections must be of the same size, "
527 : << second_order_boundary_corrections[0].number_of_grid_points());
528 : if (fd_derivative_order != DerivativeOrder::Two) {
529 : if (not high_order_corrections->has_value()) {
530 : (*high_order_corrections) =
531 : make_array<Dim>(Variables<tmpl::list<EvolvedVarsTags...>>{
532 : second_order_boundary_corrections[0].number_of_grid_points()});
533 : }
534 : ASSERT(
535 : high_order_corrections->has_value() and
536 : alg::all_of(high_order_corrections->value(),
537 : [expected_size =
538 : second_order_boundary_corrections[0]
539 : .number_of_grid_points()](const auto& e) {
540 : return e.number_of_grid_points() == expected_size;
541 : }),
542 : "The high_order_corrections must all have size "
543 : << second_order_boundary_corrections[0].number_of_grid_points());
544 : DirectionMap<Dim, Variables<FluxesTags>> flux_neighbor_data{};
545 : set_cartesian_neighbor_cell_centered_fluxes(
546 : make_not_null(&flux_neighbor_data), all_ghost_data, subcell_mesh,
547 : ghost_zone_size, number_of_rdmp_values_in_ghost_data);
548 :
549 : cartesian_high_order_fluxes_using_nodes(
550 : make_not_null(&(high_order_corrections->value())),
551 : second_order_boundary_corrections, cell_centered_fluxes.value(),
552 : flux_neighbor_data, subcell_mesh, ghost_zone_size,
553 : fd_derivative_order, reconstruction_order);
554 : }
555 : }
556 : }
557 : } // namespace fd
|
