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 : const std::array<gsl::span<std::uint8_t>, Dim>& reconstruction_order = {},
100 : const bool aligned_coordinates = true,
101 : const std::array<tnsr::i<DataVector, Dim, Frame::Inertial>, Dim>& normal =
102 : {},
103 : const std::array<
104 : DirectionMap<Dim, tnsr::i<DataVector, Dim, Frame::Inertial>>, Dim>&
105 : ghost_cell_normal = {}) {
106 : using std::min;
107 : static_assert(Dim > 0 and Dim <= 3, "Dim must be 1, 2, or 3");
108 : constexpr int max_correction_order = 10;
109 : static_assert(static_cast<int>(DerivOrder) <= max_correction_order);
110 : constexpr size_t stencil_size = static_cast<int>(DerivOrder) < 0
111 : ? 8
112 : : (static_cast<size_t>(DerivOrder) - 2);
113 : const size_t correction_width =
114 : min(static_cast<size_t>(DerivOrder) / 2 - 1,
115 : min(number_of_ghost_cells, stencil_size / 2));
116 : ASSERT(correction_width <= number_of_ghost_cells,
117 : "The width of the derivative correction ("
118 : << correction_width
119 : << ") must be less than or equal to the number of ghost cells "
120 : << number_of_ghost_cells);
121 : ASSERT(alg::all_of(reconstruction_order,
122 : [](const auto& t) { return not t.empty(); }) or
123 : static_cast<int>(DerivOrder) > 0,
124 : "For adaptive derivative orders the reconstruction_order must be set");
125 : ASSERT(normal[0][0].size() != 0 or aligned_coordinates,
126 : "Normal vectors must be specified when coordinate systems are not "
127 : "aligned");
128 : for (size_t dim = 0; dim < Dim; ++dim) {
129 : gsl::at(*high_order_boundary_corrections_in_logical_direction, dim)
130 : .initialize(
131 : gsl::at(second_order_boundary_corrections_in_logical_direction, dim)
132 : .number_of_grid_points());
133 : }
134 :
135 : // Reconstruction order is always first-varying fastest since we don't
136 : // transpose that back to {x,y,z} ordering.
137 : Index<Dim> reconstruction_extents = subcell_mesh.extents();
138 : reconstruction_extents[0] += 2;
139 :
140 : const auto impl = [&cell_centered_inertial_flux, &ghost_cell_inertial_flux,
141 : &high_order_boundary_corrections_in_logical_direction,
142 : number_of_ghost_cells,
143 : &second_order_boundary_corrections_in_logical_direction,
144 : &subcell_mesh, &correction_width, &reconstruction_order,
145 : &reconstruction_extents, &normal,
146 : &ghost_cell_normal]<typename tag, size_t dim,
147 : bool AlignedCoordinates>(
148 : tmpl::type_<tag> /*meta*/,
149 : std::integral_constant<size_t, dim> /*meta*/,
150 : std::integral_constant<bool,
151 : AlignedCoordinates> /*meta*/) {
152 : (void)reconstruction_extents;
153 :
154 : auto& high_order_var_correction =
155 : get<tag>((*high_order_boundary_corrections_in_logical_direction)[dim]);
156 : const auto& second_order_var_correction =
157 : get<tag>(second_order_boundary_corrections_in_logical_direction[dim]);
158 : const auto& recons_order = reconstruction_order[dim];
159 : const auto& cell_centered_flux =
160 : get<::Tags::Flux<tag, tmpl::size_t<Dim>, Frame::Inertial>>(
161 : cell_centered_inertial_flux);
162 : const auto& lower_neighbor_cell_centered_flux =
163 : get<::Tags::Flux<tag, tmpl::size_t<Dim>, Frame::Inertial>>(
164 : ghost_cell_inertial_flux.at(Direction<Dim>{dim, Side::Lower}));
165 : const auto& upper_neighbor_cell_centered_flux =
166 : get<::Tags::Flux<tag, tmpl::size_t<Dim>, Frame::Inertial>>(
167 : ghost_cell_inertial_flux.at(Direction<Dim>{dim, Side::Upper}));
168 : using FluxTensor = std::decay_t<decltype(cell_centered_flux)>;
169 : const auto& subcell_extents = subcell_mesh.extents();
170 : auto subcell_face_extents = subcell_extents;
171 : ++subcell_face_extents[dim];
172 : auto neighbor_extents = subcell_extents;
173 : neighbor_extents[dim] = number_of_ghost_cells;
174 : const size_t number_of_components = second_order_var_correction.size();
175 : for (size_t storage_index = 0; storage_index < number_of_components;
176 : ++storage_index) {
177 : std::array<size_t, Dim> flux_storage_indices;
178 : for (size_t n = 0; n < Dim; ++n) {
179 : const auto temp_multi_index = prepend(
180 : second_order_var_correction.get_tensor_index(storage_index), n);
181 : flux_storage_indices[n] =
182 : FluxTensor::get_storage_index(temp_multi_index);
183 : }
184 : for (size_t k = 0; k < (Dim == 3 ? subcell_face_extents[2] : 1); ++k) {
185 : for (size_t j = 0; j < (Dim >= 2 ? subcell_face_extents[1] : 1); ++j) {
186 : for (size_t i = 0; i < subcell_face_extents[0]; ++i) {
187 : const Index<Dim> face_index = [i, j, k]() -> Index<Dim> {
188 : if constexpr (Dim == 3) {
189 : return Index<Dim>{i, j, k};
190 : } else if constexpr (Dim == 2) {
191 : (void)k;
192 : return Index<Dim>{i, j};
193 : } else {
194 : (void)k, (void)j;
195 : return Index<Dim>{i};
196 : }
197 : }();
198 : const size_t face_storage_index =
199 : collapsed_index(face_index, subcell_face_extents);
200 : Index<Dim> neighbor_index{};
201 : for (size_t l = 0; l < Dim; ++l) {
202 : if (l != dim) {
203 : neighbor_index[l] = face_index[l];
204 : }
205 : }
206 :
207 : double& correction =
208 : high_order_var_correction[storage_index][face_storage_index] =
209 : 0.0;
210 :
211 : std::array<double, stencil_size> cell_centered_fluxes_for_stencil{};
212 : // fill if we have to retrieve from lower neighbor; compute a
213 : // dot-product with normal components when coordinates are not
214 : // aligned
215 : size_t stencil_index = 0;
216 : for (int grid_index = static_cast<int>(face_index[dim]) -
217 : static_cast<int>(correction_width);
218 : grid_index < static_cast<int>(face_index[dim]) +
219 : static_cast<int>(correction_width);
220 : ++grid_index, ++stencil_index) {
221 : if (grid_index < 0) {
222 : neighbor_index[dim] =
223 : number_of_ghost_cells + static_cast<size_t>(grid_index);
224 : const size_t idx_flat =
225 : collapsed_index(neighbor_index, neighbor_extents);
226 : double normal_flux = 0.0;
227 : for (size_t n = 0; n < Dim; ++n) {
228 : const double normal_in_direction =
229 : AlignedCoordinates
230 : ? (n == dim ? 1.0 : 0.0)
231 : : ghost_cell_normal[dim]
232 : .at(Direction<Dim>{dim, Side::Lower})
233 : .get(n)[idx_flat];
234 : normal_flux +=
235 : lower_neighbor_cell_centered_flux[flux_storage_indices[n]]
236 : [idx_flat] *
237 : normal_in_direction;
238 : }
239 : gsl::at(cell_centered_fluxes_for_stencil, stencil_index) =
240 : normal_flux;
241 : } else if (grid_index >= static_cast<int>(subcell_extents[dim])) {
242 : neighbor_index[dim] = static_cast<size_t>(
243 : grid_index - static_cast<int>(subcell_extents[dim]));
244 : const size_t idx_flat =
245 : collapsed_index(neighbor_index, neighbor_extents);
246 : double normal_flux = 0.0;
247 : for (size_t n = 0; n < Dim; ++n) {
248 : const double normal_in_direction =
249 : AlignedCoordinates
250 : ? (n == dim ? 1.0 : 0.0)
251 : : ghost_cell_normal[dim]
252 : .at(Direction<Dim>{dim, Side::Upper})
253 : .get(n)[idx_flat];
254 : normal_flux +=
255 : upper_neighbor_cell_centered_flux[flux_storage_indices[n]]
256 : [idx_flat] *
257 : normal_in_direction;
258 : }
259 : gsl::at(cell_centered_fluxes_for_stencil, stencil_index) =
260 : normal_flux;
261 : } else {
262 : Index<Dim> volume_index = face_index;
263 : volume_index[dim] = static_cast<size_t>(grid_index);
264 : const size_t idx_flat =
265 : collapsed_index(volume_index, subcell_extents);
266 : double normal_flux = 0.0;
267 : for (size_t n = 0; n < Dim; ++n) {
268 : const double normal_in_direction =
269 : AlignedCoordinates ? (n == dim ? 1.0 : 0.0)
270 : : normal[dim].get(n)[idx_flat];
271 : normal_flux +=
272 : cell_centered_flux[flux_storage_indices[n]][idx_flat] *
273 : normal_in_direction;
274 : }
275 : gsl::at(cell_centered_fluxes_for_stencil, stencil_index) =
276 : normal_flux;
277 : }
278 : }
279 :
280 : size_t lower_neighbor_index = std::numeric_limits<size_t>::max();
281 : size_t upper_neighbor_index = std::numeric_limits<size_t>::max();
282 : if constexpr (static_cast<int>(DerivOrder) < 0) {
283 : Index<Dim> lower_n{};
284 : Index<Dim> upper_n{};
285 : if constexpr (dim == 0) {
286 : if constexpr (Dim == 1) {
287 : lower_n = Index<Dim>{i};
288 : upper_n = Index<Dim>{i + 1};
289 : } else if constexpr (Dim == 2) {
290 : lower_n = Index<Dim>{i, j};
291 : upper_n = Index<Dim>{i + 1, j};
292 : } else if constexpr (Dim == 3) {
293 : lower_n = Index<Dim>{i, j, k};
294 : upper_n = Index<Dim>{i + 1, j, k};
295 : }
296 : } else if constexpr (dim == 1) {
297 : if constexpr (Dim == 2) {
298 : lower_n = Index<Dim>{j, i};
299 : upper_n = Index<Dim>{j + 1, i};
300 : } else if constexpr (Dim == 3) {
301 : lower_n = Index<Dim>{j, k, i};
302 : upper_n = Index<Dim>{j + 1, k, i};
303 : }
304 : } else if constexpr (dim == 2) {
305 : if constexpr (Dim == 3) {
306 : lower_n = Index<Dim>{k, i, j};
307 : upper_n = Index<Dim>{k + 1, i, j};
308 : }
309 : }
310 : lower_neighbor_index =
311 : collapsed_index(lower_n, reconstruction_extents);
312 : upper_neighbor_index =
313 : collapsed_index(upper_n, reconstruction_extents);
314 : }
315 :
316 : if (static_cast<int>(DerivOrder) >= 10 or
317 : (static_cast<int>(DerivOrder) < 0 and
318 : min(recons_order[lower_neighbor_index],
319 : recons_order[upper_neighbor_index]) >= 9)) {
320 : correction -=
321 : 5.6689342403628117913e-5 *
322 : (gsl::at(cell_centered_fluxes_for_stencil,
323 : correction_width - 4) +
324 : gsl::at(cell_centered_fluxes_for_stencil,
325 : correction_width + 3) -
326 : 9.8 * (gsl::at(cell_centered_fluxes_for_stencil,
327 : correction_width - 3) +
328 : gsl::at(cell_centered_fluxes_for_stencil,
329 : correction_width + 2)) +
330 : 49.0 * (gsl::at(cell_centered_fluxes_for_stencil,
331 : correction_width - 2) +
332 : gsl::at(cell_centered_fluxes_for_stencil,
333 : correction_width + 1)) -
334 : 245.0 * (gsl::at(cell_centered_fluxes_for_stencil,
335 : correction_width - 1) +
336 : gsl::at(cell_centered_fluxes_for_stencil,
337 : correction_width)) -
338 : 409.6 * second_order_var_correction[storage_index]
339 : [face_storage_index]);
340 : }
341 : if (static_cast<int>(DerivOrder) >= 8 or
342 : (static_cast<int>(DerivOrder) < 0 and
343 : min(recons_order[lower_neighbor_index],
344 : recons_order[upper_neighbor_index]) >= 7)) {
345 : correction +=
346 : 4.7619047619047619047e-4 *
347 : (gsl::at(cell_centered_fluxes_for_stencil,
348 : correction_width - 3) +
349 : gsl::at(cell_centered_fluxes_for_stencil,
350 : correction_width + 2) -
351 : 8.3333333333333333333 *
352 : (gsl::at(cell_centered_fluxes_for_stencil,
353 : correction_width - 2) +
354 : gsl::at(cell_centered_fluxes_for_stencil,
355 : correction_width + 1)) +
356 : 50.0 * (gsl::at(cell_centered_fluxes_for_stencil,
357 : correction_width - 1) +
358 : gsl::at(cell_centered_fluxes_for_stencil,
359 : correction_width)) +
360 : 85.333333333333333333 *
361 : second_order_var_correction[storage_index]
362 : [face_storage_index]);
363 : }
364 : if (static_cast<int>(DerivOrder) >= 6 or
365 : (static_cast<int>(DerivOrder) < 0 and
366 : min(recons_order[lower_neighbor_index],
367 : recons_order[upper_neighbor_index]) >=
368 : (DerivOrder ==
369 : DerivativeOrder::OneHigherThanReconsButFiveToFour
370 : ? 6
371 : : 5))) {
372 : correction -=
373 : 5.5555555555555555555e-3 *
374 : (gsl::at(cell_centered_fluxes_for_stencil,
375 : correction_width - 2) +
376 : gsl::at(cell_centered_fluxes_for_stencil,
377 : correction_width + 1) -
378 : 9.0 * (gsl::at(cell_centered_fluxes_for_stencil,
379 : correction_width - 1) +
380 : gsl::at(cell_centered_fluxes_for_stencil,
381 : correction_width)) -
382 : 16.0 * second_order_var_correction[storage_index]
383 : [face_storage_index]);
384 : }
385 : if (static_cast<int>(DerivOrder) >= 4 or
386 : (static_cast<int>(DerivOrder) < 0 and
387 : min(recons_order[lower_neighbor_index],
388 : recons_order[upper_neighbor_index]) >= 3)) {
389 : correction +=
390 : 0.166666666666666666 *
391 : (gsl::at(cell_centered_fluxes_for_stencil,
392 : correction_width - 1) +
393 : gsl::at(cell_centered_fluxes_for_stencil, correction_width) +
394 : 2.0 * second_order_var_correction[storage_index]
395 : [face_storage_index]);
396 : }
397 :
398 : // Add second-order correction last
399 : correction +=
400 : second_order_var_correction[storage_index][face_storage_index];
401 : }
402 : }
403 : }
404 : }
405 : };
406 :
407 : if (aligned_coordinates) {
408 : EXPAND_PACK_LEFT_TO_RIGHT(impl(tmpl::type_<EvolvedVarsTags>{},
409 : std::integral_constant<size_t, 0>{},
410 : std::true_type{}));
411 : if constexpr (Dim > 1) {
412 : EXPAND_PACK_LEFT_TO_RIGHT(impl(tmpl::type_<EvolvedVarsTags>{},
413 : std::integral_constant<size_t, 1>{},
414 : std::true_type{}));
415 : if constexpr (Dim > 2) {
416 : EXPAND_PACK_LEFT_TO_RIGHT(impl(tmpl::type_<EvolvedVarsTags>{},
417 : std::integral_constant<size_t, 2>{},
418 : std::true_type{}));
419 : }
420 : }
421 : } else {
422 : EXPAND_PACK_LEFT_TO_RIGHT(impl(tmpl::type_<EvolvedVarsTags>{},
423 : std::integral_constant<size_t, 0>{},
424 : std::false_type{}));
425 : if constexpr (Dim > 1) {
426 : EXPAND_PACK_LEFT_TO_RIGHT(impl(tmpl::type_<EvolvedVarsTags>{},
427 : std::integral_constant<size_t, 1>{},
428 : std::false_type{}));
429 : if constexpr (Dim > 2) {
430 : EXPAND_PACK_LEFT_TO_RIGHT(impl(tmpl::type_<EvolvedVarsTags>{},
431 : std::integral_constant<size_t, 2>{},
432 : std::false_type{}));
433 : }
434 : }
435 : }
436 : }
437 :
438 : template <size_t Dim, typename... EvolvedVarsTags>
439 1 : void cartesian_high_order_fluxes_using_nodes(
440 : const gsl::not_null<
441 : std::array<Variables<tmpl::list<EvolvedVarsTags...>>, Dim>*>
442 : high_order_boundary_corrections_in_logical_direction,
443 :
444 : const std::array<Variables<tmpl::list<EvolvedVarsTags...>>, Dim>&
445 : second_order_boundary_corrections_in_logical_direction,
446 : const Variables<tmpl::list<
447 : ::Tags::Flux<EvolvedVarsTags, tmpl::size_t<Dim>, Frame::Inertial>...>>&
448 : cell_centered_inertial_flux,
449 : const DirectionMap<
450 : Dim, Variables<tmpl::list<::Tags::Flux<
451 : EvolvedVarsTags, tmpl::size_t<Dim>, Frame::Inertial>...>>>&
452 : ghost_cell_inertial_flux,
453 : const Mesh<Dim>& subcell_mesh, const size_t number_of_ghost_cells,
454 : const DerivativeOrder derivative_order,
455 : const std::array<gsl::span<std::uint8_t>, Dim>& reconstruction_order = {},
456 : const bool aligned_coordinates = true,
457 : const std::array<tnsr::i<DataVector, Dim, Frame::Inertial>, Dim>& normal =
458 : {},
459 : const std::array<
460 : DirectionMap<Dim, tnsr::i<DataVector, Dim, Frame::Inertial>>, Dim>&
461 : ghost_cell_normal = {}) {
462 : switch (derivative_order) {
463 : case DerivativeOrder::OneHigherThanRecons:
464 : cartesian_high_order_fluxes_using_nodes<
465 : DerivativeOrder::OneHigherThanRecons>(
466 : high_order_boundary_corrections_in_logical_direction,
467 : second_order_boundary_corrections_in_logical_direction,
468 : cell_centered_inertial_flux, ghost_cell_inertial_flux, subcell_mesh,
469 : number_of_ghost_cells, reconstruction_order, aligned_coordinates,
470 : normal, ghost_cell_normal);
471 : break;
472 : case DerivativeOrder::OneHigherThanReconsButFiveToFour:
473 : cartesian_high_order_fluxes_using_nodes<
474 : DerivativeOrder::OneHigherThanReconsButFiveToFour>(
475 : high_order_boundary_corrections_in_logical_direction,
476 : second_order_boundary_corrections_in_logical_direction,
477 : cell_centered_inertial_flux, ghost_cell_inertial_flux, subcell_mesh,
478 : number_of_ghost_cells, reconstruction_order, aligned_coordinates,
479 : normal, ghost_cell_normal);
480 : break;
481 : case DerivativeOrder::Two:
482 : cartesian_high_order_fluxes_using_nodes<DerivativeOrder::Two>(
483 : high_order_boundary_corrections_in_logical_direction,
484 : second_order_boundary_corrections_in_logical_direction,
485 : cell_centered_inertial_flux, ghost_cell_inertial_flux, subcell_mesh,
486 : number_of_ghost_cells, reconstruction_order, aligned_coordinates,
487 : normal, ghost_cell_normal);
488 : break;
489 : case DerivativeOrder::Four:
490 : cartesian_high_order_fluxes_using_nodes<DerivativeOrder::Four>(
491 : high_order_boundary_corrections_in_logical_direction,
492 : second_order_boundary_corrections_in_logical_direction,
493 : cell_centered_inertial_flux, ghost_cell_inertial_flux, subcell_mesh,
494 : number_of_ghost_cells, reconstruction_order, aligned_coordinates,
495 : normal, ghost_cell_normal);
496 : break;
497 : case DerivativeOrder::Six:
498 : cartesian_high_order_fluxes_using_nodes<DerivativeOrder::Six>(
499 : high_order_boundary_corrections_in_logical_direction,
500 : second_order_boundary_corrections_in_logical_direction,
501 : cell_centered_inertial_flux, ghost_cell_inertial_flux, subcell_mesh,
502 : number_of_ghost_cells, reconstruction_order, aligned_coordinates,
503 : normal, ghost_cell_normal);
504 : break;
505 : case DerivativeOrder::Eight:
506 : cartesian_high_order_fluxes_using_nodes<DerivativeOrder::Eight>(
507 : high_order_boundary_corrections_in_logical_direction,
508 : second_order_boundary_corrections_in_logical_direction,
509 : cell_centered_inertial_flux, ghost_cell_inertial_flux, subcell_mesh,
510 : number_of_ghost_cells, reconstruction_order, aligned_coordinates,
511 : normal, ghost_cell_normal);
512 : break;
513 : case DerivativeOrder::Ten:
514 : cartesian_high_order_fluxes_using_nodes<DerivativeOrder::Ten>(
515 : high_order_boundary_corrections_in_logical_direction,
516 : second_order_boundary_corrections_in_logical_direction,
517 : cell_centered_inertial_flux, ghost_cell_inertial_flux, subcell_mesh,
518 : number_of_ghost_cells, reconstruction_order, aligned_coordinates,
519 : normal, ghost_cell_normal);
520 : break;
521 : default:
522 : ERROR("Unsupported correction order " << derivative_order);
523 : };
524 : }
525 : /// @}
526 :
527 : /*!
528 : * \brief Fill the `flux_neighbor_data` with pointers into the
529 : * `all_ghost_data`.
530 : *
531 : * The `all_ghost_data` is stored in the tag
532 : * `evolution::dg::subcell::Tags::GhostDataForReconstruction`, and the
533 : * `ghost_zone_size` should come from the FD reconstructor.
534 : */
535 : template <size_t Dim, typename FluxesTags>
536 1 : void set_cartesian_neighbor_cell_centered_fluxes(
537 : const gsl::not_null<DirectionMap<Dim, Variables<FluxesTags>>*>
538 : flux_neighbor_data,
539 : const DirectionalIdMap<Dim, evolution::dg::subcell::GhostData>&
540 : all_ghost_data,
541 : const Mesh<Dim>& subcell_mesh, const size_t ghost_zone_size,
542 : const size_t number_of_rdmp_values_in_ghost_data) {
543 : for (const auto& [direction_id, ghost_data] : all_ghost_data) {
544 : const size_t neighbor_flux_size =
545 : subcell_mesh.number_of_grid_points() /
546 : subcell_mesh.extents(direction_id.direction().dimension()) *
547 : ghost_zone_size *
548 : Variables<FluxesTags>::number_of_independent_components;
549 : const DataVector& neighbor_data =
550 : ghost_data.neighbor_ghost_data_for_reconstruction();
551 : (*flux_neighbor_data)[direction_id.direction()].set_data_ref(
552 : // NOLINTNEXTLINE(cppcoreguidelines-pro-type-const-cast)
553 : const_cast<double*>(std::next(
554 : neighbor_data.data(),
555 : static_cast<std::ptrdiff_t>(neighbor_data.size() -
556 : number_of_rdmp_values_in_ghost_data -
557 : neighbor_flux_size))),
558 : neighbor_flux_size);
559 : }
560 : }
561 :
562 : /*!
563 : * \brief Computes the high-order Cartesian flux corrections if necessary.
564 : *
565 : * The `cell_centered_fluxes` is stored in the tag
566 : * `evolution::dg::subcell::Tags::CellCenteredFlux`, `fd_derivative_order` is
567 : * from `evolution::dg::subcell::Tags::SubcellOptions`
568 : * (`.finite_difference_derivative_order()`), the `all_ghost_data`
569 : * is stored in the tag
570 : * `evolution::dg::subcell::Tags::GhostDataForReconstruction`, the
571 : * `ghost_zone_size` should come from the FD reconstructor.
572 : *
573 : * By default we assume no RDMP data is in the `ghost_data` buffer. In the
574 : * future we will want to update how we store the data in order to eliminate
575 : * more memory allocations and copies, in which case that value will be
576 : * non-zero.
577 : *
578 : * \note `high_order_corrections` must either not have a value or have all
579 : * elements be of the same size as
580 : * `second_order_boundary_corrections[0].number_of_grid_points()`, where we've
581 : * assumed `second_order_boundary_corrections` is the same in all directions.
582 : */
583 : template <size_t Dim, typename... EvolvedVarsTags,
584 : typename FluxesTags = tmpl::list<::Tags::Flux<
585 : EvolvedVarsTags, tmpl::size_t<Dim>, Frame::Inertial>...>>
586 1 : void cartesian_high_order_flux_corrections(
587 : const gsl::not_null<std::optional<
588 : std::array<Variables<tmpl::list<EvolvedVarsTags...>>, Dim>>*>
589 : high_order_corrections,
590 :
591 : const std::optional<Variables<FluxesTags>>& cell_centered_fluxes,
592 : const std::array<Variables<tmpl::list<EvolvedVarsTags...>>, Dim>&
593 : second_order_boundary_corrections,
594 : const fd::DerivativeOrder& fd_derivative_order,
595 : const DirectionalIdMap<Dim, evolution::dg::subcell::GhostData>&
596 : all_ghost_data,
597 : const Mesh<Dim>& subcell_mesh, const size_t ghost_zone_size,
598 : const std::array<gsl::span<std::uint8_t>, Dim>& reconstruction_order = {},
599 : const bool aligned_coordinates = true,
600 : const std::array<tnsr::i<DataVector, Dim, Frame::Inertial>, Dim>& normal =
601 : {},
602 : const std::array<
603 : DirectionMap<Dim, tnsr::i<DataVector, Dim, Frame::Inertial>>, Dim>&
604 : ghost_cell_normal = {},
605 : const size_t number_of_rdmp_values_in_ghost_data = 0) {
606 : if (cell_centered_fluxes.has_value()) {
607 : ASSERT(alg::all_of(
608 : second_order_boundary_corrections,
609 : [expected_size = second_order_boundary_corrections[0]
610 : .number_of_grid_points()](const auto& e) {
611 : return e.number_of_grid_points() == expected_size;
612 : }),
613 : "All second-order boundary corrections must be of the same size, "
614 : << second_order_boundary_corrections[0].number_of_grid_points());
615 : if (fd_derivative_order != DerivativeOrder::Two) {
616 : if (not high_order_corrections->has_value()) {
617 : (*high_order_corrections) =
618 : make_array<Dim>(Variables<tmpl::list<EvolvedVarsTags...>>{
619 : second_order_boundary_corrections[0].number_of_grid_points()});
620 : }
621 : ASSERT(
622 : high_order_corrections->has_value() and
623 : alg::all_of(high_order_corrections->value(),
624 : [expected_size =
625 : second_order_boundary_corrections[0]
626 : .number_of_grid_points()](const auto& e) {
627 : return e.number_of_grid_points() == expected_size;
628 : }),
629 : "The high_order_corrections must all have size "
630 : << second_order_boundary_corrections[0].number_of_grid_points());
631 : DirectionMap<Dim, Variables<FluxesTags>> flux_neighbor_data{};
632 : set_cartesian_neighbor_cell_centered_fluxes(
633 : make_not_null(&flux_neighbor_data), all_ghost_data, subcell_mesh,
634 : ghost_zone_size, number_of_rdmp_values_in_ghost_data);
635 :
636 : cartesian_high_order_fluxes_using_nodes(
637 : make_not_null(&(high_order_corrections->value())),
638 : second_order_boundary_corrections, cell_centered_fluxes.value(),
639 : flux_neighbor_data, subcell_mesh, ghost_zone_size,
640 : fd_derivative_order, reconstruction_order, aligned_coordinates,
641 : normal, ghost_cell_normal);
642 : }
643 : }
644 : }
645 : } // namespace fd
|
