SpECTRE Documentation Coverage Report
Current view: top level - NumericalAlgorithms/FiniteDifference - HighOrderFluxCorrection.hpp Hit Total Coverage
Commit: 9a905b0737f373631c1b8e8389b8f26e67fa5313 Lines: 4 5 80.0 %
Date: 2024-03-28 09:03:18
Legend: Lines: hit not hit

          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

Generated by: LCOV version 1.14