Line data    Source code

       1           0 : // Distributed under the MIT License.
       2             : // See LICENSE.txt for details.
       3             : 
       4             : #pragma once
       5             : 
       6             : #include <algorithm>
       7             : #include <cstddef>
       8             : 
       9             : #include "DataStructures/DataBox/PrefixHelpers.hpp"
      10             : #include "DataStructures/DataVector.hpp"
      11             : #include "DataStructures/Variables.hpp"
      12             : #include "NumericalAlgorithms/Spectral/Basis.hpp"
      13             : #include "NumericalAlgorithms/Spectral/Parity.hpp"
      14             : #include "NumericalAlgorithms/Spectral/ParityFromSymmetry.hpp"
      15             : #include "Utilities/ErrorHandling/Assert.hpp"
      16             : #include "Utilities/Gsl.hpp"
      17             : #include "Utilities/MemoryHelpers.hpp"
      18             : #include "Utilities/TMPL.hpp"
      19             : 
      20             : /// \cond
      21             : template <size_t>
      22             : class Index;
      23             : template <size_t>
      24             : class Mesh;
      25             : /// \endcond
      26             : 
      27             : namespace evolution::dg::subcell::fd {
      28             : namespace detail {
      29             : template <size_t Dim>
      30             : void project_impl(gsl::span<double> subcell_u, gsl::span<const double> dg_u,
      31             :                   const Mesh<Dim>& dg_mesh, const Index<Dim>& subcell_extents,
      32             :                   Spectral::Parity parity);
      33             : template <size_t Dim>
      34             : void project_to_faces_impl(gsl::span<double> subcell_u,
      35             :                            gsl::span<const double> dg_u,
      36             :                            const Mesh<Dim>& dg_mesh,
      37             :                            const Index<Dim>& subcell_extents,
      38             :                            const size_t& face_direction,
      39             :                            Spectral::Parity parity);
      40             : 
      41             : /*!
      42             :  * \brief Project `dg_u` onto the subcell grid, sorting even- and odd-parity
      43             :  * components into separate batches when a ZernikeB1 basis is present.
      44             :  *
      45             :  * For non-ZernikeB1 meshes this falls through to the default implementation.
      46             :  * The `TagList` must be the full Variables tag list.
      47             :  */
      48             : template <typename TagList, size_t Dim>
      49             : void project_impl_with_tag_list(gsl::span<double> subcell_u,
      50             :                                 gsl::span<const double> dg_u,
      51             :                                 const Mesh<Dim>& dg_mesh,
      52             :                                 const Index<Dim>& subcell_extents) {
      53             :   if (dg_mesh.basis(0) == Spectral::Basis::ZernikeB1) {
      54             :     ASSERT(Variables<TagList>::number_of_independent_components *
      55             :                    dg_mesh.number_of_grid_points() ==
      56             :                dg_u.size(),
      57             :            "Passed TagList does not have the same components, "
      58             :                << Variables<TagList>::number_of_independent_components
      59             :                << ", as dg_u holds, "
      60             :                << dg_u.size() / dg_mesh.number_of_grid_points());
      61             :     constexpr auto parity_info = Spectral::compute_parity_list<TagList>();
      62             :     constexpr auto parity_list = std::get<0>(parity_info);
      63             :     constexpr size_t num_even = std::get<1>(parity_info);
      64             :     constexpr size_t num_odd = std::get<2>(parity_info);
      65             : 
      66             :     const size_t num_dg_pts = dg_mesh.number_of_grid_points();
      67             :     const size_t num_subcell_pts = subcell_extents.product();
      68             : 
      69             :     // NOLINTNEXTLINE(modernize-avoid-c-arrays)
      70             :     auto buffer = cpp20::make_unique_for_overwrite<double[]>(
      71             :         (num_even + num_odd) * (num_dg_pts + num_subcell_pts));
      72             :     DataVector even_input{&buffer[0], num_even * num_dg_pts};
      73             :     DataVector odd_input(&buffer[num_even * num_dg_pts], num_odd * num_dg_pts);
      74             :     DataVector even_output(&buffer[(num_even + num_odd) * num_dg_pts],
      75             :                            num_even * num_subcell_pts);
      76             :     DataVector odd_output(
      77             :         &buffer[(num_even + num_odd) * num_dg_pts + num_even * num_subcell_pts],
      78             :         num_odd * num_subcell_pts);
      79             : 
      80             :     // Sort input components into even/odd parity buffers
      81             :     const double* p_in = dg_u.data();
      82             :     double* p_even_in = even_input.data();
      83             :     double* p_odd_in = odd_input.data();
      84             :     bool is_even = true;
      85             :     for (const size_t seg_size : parity_list) {
      86             :       if (seg_size == 0) {
      87             :         if (is_even) {
      88             :           is_even = false;
      89             :           continue;
      90             :         } else {
      91             :           break;
      92             :         }
      93             :       }
      94             :       if (is_even) {
      95             :         std::copy(p_in, p_in + seg_size * num_dg_pts, p_even_in);  // NOLINT
      96             :         p_even_in += seg_size * num_dg_pts;                        // NOLINT
      97             :       } else {
      98             :         std::copy(p_in, p_in + seg_size * num_dg_pts, p_odd_in);  // NOLINT
      99             :         p_odd_in += seg_size * num_dg_pts;                        // NOLINT
     100             :       }
     101             :       p_in += seg_size * num_dg_pts;  // NOLINT
     102             :       is_even = not is_even;
     103             :     }
     104             : 
     105             :     // Project each parity batch with the appropriate projection matrix
     106             :     if constexpr (num_even > 0) {
     107             :       project_impl(
     108             :           gsl::span<double>{even_output.data(), even_output.size()},
     109             :           gsl::span<const double>{even_input.data(), even_input.size()},
     110             :           dg_mesh, subcell_extents, Spectral::Parity::Even);
     111             :     }
     112             :     if constexpr (num_odd > 0) {
     113             :       project_impl(gsl::span<double>{odd_output.data(), odd_output.size()},
     114             :                    gsl::span<const double>{odd_input.data(), odd_input.size()},
     115             :                    dg_mesh, subcell_extents, Spectral::Parity::Odd);
     116             :     }
     117             : 
     118             :     // Reassemble output in original component order
     119             :     double* p_out = subcell_u.data();
     120             :     const double* p_even_out = even_output.data();
     121             :     const double* p_odd_out = odd_output.data();
     122             :     is_even = true;
     123             :     for (const size_t seg_size : parity_list) {
     124             :       if (seg_size == 0) {
     125             :         if (is_even) {
     126             :           is_even = false;
     127             :           continue;
     128             :         } else {
     129             :           break;
     130             :         }
     131             :       }
     132             :       if (is_even) {
     133             :         // NOLINTNEXTLINE
     134             :         std::copy(p_even_out, p_even_out + seg_size * num_subcell_pts, p_out);
     135             :         p_even_out += seg_size * num_subcell_pts;  // NOLINT
     136             :       } else {
     137             :         // NOLINTNEXTLINE
     138             :         std::copy(p_odd_out, p_odd_out + seg_size * num_subcell_pts, p_out);
     139             :         p_odd_out += seg_size * num_subcell_pts;  // NOLINT
     140             :       }
     141             :       p_out += seg_size * num_subcell_pts;  // NOLINT
     142             :       is_even = not is_even;
     143             :     }
     144             :     return;
     145             :   }
     146             :   project_impl(subcell_u, dg_u, dg_mesh, subcell_extents,
     147             :                Spectral::Parity::Uninitialized);
     148             : }
     149             : }  // namespace detail
     150             : 
     151             : /// @{
     152             : /*!
     153             :  * \ingroup DgSubcellGroup
     154             :  * \brief Project the variable `dg_u` onto the subcell grid with extents
     155             :  * `subcell_extents`.
     156             :  *
     157             :  * When the DG mesh uses a ZernikeB1 basis the Variables overloads deduce
     158             :  * per-component parity from the tag list automatically.  The raw `DataVector`
     159             :  * overloads accepting a `Spectral::Parity` are for single-component data
     160             :  * where the caller already knows the parity; the overloads accepting a
     161             :  * `tmpl::list<TagList>` meta-parameter project a multi-component DataVector
     162             :  * whose tensor-parity structure is encoded in `TagList`.
     163             :  *
     164             :  * \note In the return-by-`gsl::not_null` with `Variables` interface, the
     165             :  * `SubcellTagList` and the `DgTagList` must be the same when all tag prefixes
     166             :  * are removed.
     167             :  */
     168             : template <size_t Dim>
     169           1 : DataVector project(const DataVector& dg_u, const Mesh<Dim>& dg_mesh,
     170             :                    const Index<Dim>& subcell_extents,
     171             :                    Spectral::Parity parity = Spectral::Parity::Uninitialized);
     172             : 
     173             : template <size_t Dim>
     174           1 : void project(gsl::not_null<DataVector*> subcell_u, const DataVector& dg_u,
     175             :              const Mesh<Dim>& dg_mesh, const Index<Dim>& subcell_extents,
     176             :              Spectral::Parity parity = Spectral::Parity::Uninitialized);
     177             : 
     178             : template <typename TagList, size_t Dim>
     179           1 : void project(const gsl::not_null<DataVector*> subcell_u, const DataVector& dg_u,
     180             :              const Mesh<Dim>& dg_mesh, const Index<Dim>& subcell_extents,
     181             :              TagList /*meta*/) {
     182             :   ASSERT(dg_u.size() % dg_mesh.number_of_grid_points() == 0,
     183             :          "The vector dg_u must have size that is a multiple of the number of "
     184             :          "grid points "
     185             :              << dg_mesh.number_of_grid_points() << " but got " << dg_u.size());
     186             :   subcell_u->destructive_resize(subcell_extents.product() * dg_u.size() /
     187             :                                 dg_mesh.number_of_grid_points());
     188             :   detail::project_impl_with_tag_list<TagList>(
     189             :       gsl::span<double>{subcell_u->data(), subcell_u->size()},
     190             :       gsl::span<const double>{dg_u.data(), dg_u.size()}, dg_mesh,
     191             :       subcell_extents);
     192             : }
     193             : 
     194             : template <typename TagList, size_t Dim>
     195           1 : DataVector project(const DataVector& dg_u, const Mesh<Dim>& dg_mesh,
     196             :                    const Index<Dim>& subcell_extents, TagList /*meta*/) {
     197             :   ASSERT(dg_u.size() % dg_mesh.number_of_grid_points() == 0,
     198             :          "The vector dg_u must have size that is a multiple of the number of "
     199             :          "grid points "
     200             :              << dg_mesh.number_of_grid_points() << " but got " << dg_u.size());
     201             :   DataVector subcell_u(subcell_extents.product() * dg_u.size() /
     202             :                        dg_mesh.number_of_grid_points());
     203             :   project(make_not_null(&subcell_u), dg_u, dg_mesh, subcell_extents, TagList{});
     204             :   return subcell_u;
     205             : }
     206             : 
     207             : template <typename SubcellTagList, typename DgTagList, size_t Dim>
     208           1 : void project(const gsl::not_null<Variables<SubcellTagList>*> subcell_u,
     209             :              const Variables<DgTagList>& dg_u, const Mesh<Dim>& dg_mesh,
     210             :              const Index<Dim>& subcell_extents) {
     211             :   static_assert(
     212             :       std::is_same_v<
     213             :           tmpl::transform<SubcellTagList,
     214             :                           tmpl::bind<db::remove_all_prefixes, tmpl::_1>>,
     215             :           tmpl::transform<DgTagList,
     216             :                           tmpl::bind<db::remove_all_prefixes, tmpl::_1>>>,
     217             :       "DG and subcell tag lists must be the same once prefix tags "
     218             :       "are removed.");
     219             :   ASSERT(dg_u.number_of_grid_points() == dg_mesh.number_of_grid_points(),
     220             :          "dg_u has incorrect size " << dg_u.number_of_grid_points()
     221             :                                     << " since the mesh is size "
     222             :                                     << dg_mesh.number_of_grid_points());
     223             :   if (UNLIKELY(subcell_u->number_of_grid_points() !=
     224             :                subcell_extents.product())) {
     225             :     subcell_u->initialize(subcell_extents.product());
     226             :   }
     227             :   detail::project_impl_with_tag_list<DgTagList>(
     228             :       gsl::span<double>{subcell_u->data(), subcell_u->size()},
     229             :       gsl::span<const double>{dg_u.data(), dg_u.size()}, dg_mesh,
     230             :       subcell_extents);
     231             : }
     232             : 
     233             : template <typename TagList, size_t Dim>
     234           1 : Variables<TagList> project(const Variables<TagList>& dg_u,
     235             :                            const Mesh<Dim>& dg_mesh,
     236             :                            const Index<Dim>& subcell_extents) {
     237             :   Variables<TagList> subcell_u(subcell_extents.product());
     238             :   project(make_not_null(&subcell_u), dg_u, dg_mesh, subcell_extents);
     239             :   return subcell_u;
     240             : }
     241             : 
     242             : template <size_t Dim>
     243           1 : DataVector project_to_faces(const DataVector& dg_u, const Mesh<Dim>& dg_mesh,
     244             :                             const Index<Dim>& subcell_extents,
     245             :                             const size_t& face_direction,
     246             :                             Spectral::Parity parity);
     247             : 
     248             : template <size_t Dim>
     249           1 : void project_to_faces(gsl::not_null<DataVector*> subcell_u,
     250             :                       const DataVector& dg_u, const Mesh<Dim>& dg_mesh,
     251             :                       const Index<Dim>& subcell_extents,
     252             :                       const size_t& face_direction, Spectral::Parity parity);
     253             : 
     254             : template <typename SubcellTagList, typename DgTagList, size_t Dim>
     255           1 : void project_to_faces(const gsl::not_null<Variables<SubcellTagList>*> subcell_u,
     256             :                       const Variables<DgTagList>& dg_u,
     257             :                       const Mesh<Dim>& dg_mesh,
     258             :                       const Index<Dim>& subcell_extents,
     259             :                       const size_t& face_direction, Spectral::Parity parity) {
     260             :   static_assert(
     261             :       std::is_same_v<
     262             :           tmpl::transform<SubcellTagList,
     263             :                           tmpl::bind<db::remove_all_prefixes, tmpl::_1>>,
     264             :           tmpl::transform<DgTagList,
     265             :                           tmpl::bind<db::remove_all_prefixes, tmpl::_1>>>,
     266             :       "DG and subcell tag lists must be the same once prefix tags "
     267             :       "are removed.");
     268             :   ASSERT(dg_u.number_of_grid_points() == dg_mesh.number_of_grid_points(),
     269             :          "dg_u has incorrect size " << dg_u.number_of_grid_points()
     270             :                                     << " since the mesh is size "
     271             :                                     << dg_mesh.number_of_grid_points());
     272             :   if (UNLIKELY(subcell_u->number_of_grid_points() !=
     273             :                subcell_extents.product())) {
     274             :     subcell_u->initialize(subcell_extents.product());
     275             :   }
     276             :   detail::project_to_faces_impl(
     277             :       gsl::span<double>{subcell_u->data(), subcell_u->size()},
     278             :       gsl::span<const double>{dg_u.data(), dg_u.size()}, dg_mesh,
     279             :       subcell_extents, face_direction, parity);
     280             : }
     281             : 
     282             : template <typename TagList, size_t Dim>
     283           1 : Variables<TagList> project_to_faces(const Variables<TagList>& dg_u,
     284             :                                     const Mesh<Dim>& dg_mesh,
     285             :                                     const Index<Dim>& subcell_extents,
     286             :                                     const size_t& face_direction,
     287             :                                     Spectral::Parity parity) {
     288             :   Variables<TagList> subcell_u(subcell_extents.product());
     289             :   project_to_faces(make_not_null(&subcell_u), dg_u, dg_mesh, subcell_extents,
     290             :                    face_direction, parity);
     291             :   return subcell_u;
     292             : }
     293             : /// @}
     294             : }  // namespace evolution::dg::subcell::fd