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/DataVector.hpp"
      10             : #include "DataStructures/Variables.hpp"
      11             : #include "Evolution/DgSubcell/ReconstructionMethod.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 reconstruct_impl(gsl::span<double> dg_u,
      31             :                       gsl::span<const double> subcell_u_times_projected_det_jac,
      32             :                       const Mesh<Dim>& dg_mesh,
      33             :                       const Index<Dim>& subcell_extents,
      34             :                       ReconstructionMethod reconstruction_method,
      35             :                       Spectral::Parity parity);
      36             : 
      37             : /*!
      38             :  * \brief Reconstruct `subcell_u` onto the DG grid, sorting even- and
      39             :  * odd-parity components into separate batches when a ZernikeB1 basis is
      40             :  * present.
      41             :  *
      42             :  * For non-ZernikeB1 meshes this falls through to `reconstruct_impl` with
      43             :  * `Parity::Uninitialized`. The `TagList` must be the full Variables tag list
      44             :  * so that `Spectral::compute_parity_list` can determine per-component parity.
      45             :  */
      46             : template <typename TagList, size_t Dim>
      47             : void reconstruct_impl_with_tag_list(
      48             :     gsl::span<double> dg_u, gsl::span<const double> subcell_u,
      49             :     const Mesh<Dim>& dg_mesh, const Index<Dim>& subcell_extents,
      50             :     const ReconstructionMethod reconstruction_method) {
      51             :   if (dg_mesh.basis(0) == Spectral::Basis::ZernikeB1) {
      52             :     ASSERT(Variables<TagList>::number_of_independent_components *
      53             :                    dg_mesh.number_of_grid_points() ==
      54             :                dg_u.size(),
      55             :            "Passed TagList does not have the same components, "
      56             :                << Variables<TagList>::number_of_independent_components
      57             :                << ", as dg_u holds, "
      58             :                << dg_u.size() / dg_mesh.number_of_grid_points());
      59             :     constexpr auto parity_info = Spectral::compute_parity_list<TagList>();
      60             :     constexpr auto parity_list = std::get<0>(parity_info);
      61             :     constexpr size_t num_even = std::get<1>(parity_info);
      62             :     constexpr size_t num_odd = std::get<2>(parity_info);
      63             : 
      64             :     const size_t num_dg_pts = dg_mesh.number_of_grid_points();
      65             :     const size_t num_subcell_pts = subcell_extents.product();
      66             : 
      67             :     // NOLINTNEXTLINE(modernize-avoid-c-arrays)
      68             :     auto buffer = cpp20::make_unique_for_overwrite<double[]>(
      69             :         (num_even + num_odd) * (num_subcell_pts + num_dg_pts));
      70             :     DataVector even_subcell_input{&buffer[0], num_even * num_subcell_pts};
      71             :     DataVector odd_subcell_input{&buffer[num_even * num_subcell_pts],
      72             :                                  num_odd * num_subcell_pts};
      73             :     DataVector even_dg_output{&buffer[(num_even + num_odd) * num_subcell_pts],
      74             :                               num_even * num_dg_pts};
      75             :     DataVector odd_dg_output{
      76             :         &buffer[(num_even + num_odd) * num_subcell_pts + num_even * num_dg_pts],
      77             :         num_odd * num_dg_pts};
      78             : 
      79             :     // Sort subcell input into even/odd parity buffers
      80             :     const double* p_in = subcell_u.data();
      81             :     double* p_even_in = even_subcell_input.data();
      82             :     double* p_odd_in = odd_subcell_input.data();
      83             :     bool is_even = true;
      84             :     for (const size_t seg_size : parity_list) {
      85             :       if (seg_size == 0) {
      86             :         if (is_even) {
      87             :           is_even = false;
      88             :           continue;
      89             :         } else {
      90             :           break;
      91             :         }
      92             :       }
      93             :       if (is_even) {
      94             :         std::copy(p_in, p_in + seg_size * num_subcell_pts,  // NOLINT
      95             :                   p_even_in);
      96             :         p_even_in += seg_size * num_subcell_pts;  // NOLINT
      97             :       } else {
      98             :         std::copy(p_in, p_in + seg_size * num_subcell_pts,  // NOLINT
      99             :                   p_odd_in);
     100             :         p_odd_in += seg_size * num_subcell_pts;  // NOLINT
     101             :       }
     102             :       p_in += seg_size * num_subcell_pts;  // NOLINT
     103             :       is_even = not is_even;
     104             :     }
     105             : 
     106             :     // Reconstruct each parity batch with the appropriate reconstruction
     107             :     // matrix
     108             :     if constexpr (num_even > 0) {
     109             :       reconstruct_impl(
     110             :           gsl::span<double>{even_dg_output.data(), even_dg_output.size()},
     111             :           gsl::span<const double>{even_subcell_input.data(),
     112             :                                   even_subcell_input.size()},
     113             :           dg_mesh, subcell_extents, reconstruction_method,
     114             :           Spectral::Parity::Even);
     115             :     }
     116             :     if constexpr (num_odd > 0) {
     117             :       reconstruct_impl(
     118             :           gsl::span<double>{odd_dg_output.data(), odd_dg_output.size()},
     119             :           gsl::span<const double>{odd_subcell_input.data(),
     120             :                                   odd_subcell_input.size()},
     121             :           dg_mesh, subcell_extents, reconstruction_method,
     122             :           Spectral::Parity::Odd);
     123             :     }
     124             : 
     125             :     // Reassemble output in original component order
     126             :     double* p_out = dg_u.data();
     127             :     const double* p_even_out = even_dg_output.data();
     128             :     const double* p_odd_out = odd_dg_output.data();
     129             :     is_even = true;
     130             :     for (const size_t seg_size : parity_list) {
     131             :       if (seg_size == 0) {
     132             :         if (is_even) {
     133             :           is_even = false;
     134             :           continue;
     135             :         } else {
     136             :           break;
     137             :         }
     138             :       }
     139             :       if (is_even) {
     140             :         // NOLINTNEXTLINE
     141             :         std::copy(p_even_out, p_even_out + seg_size * num_dg_pts, p_out);
     142             :         p_even_out += seg_size * num_dg_pts;  // NOLINT
     143             :       } else {
     144             :         // NOLINTNEXTLINE
     145             :         std::copy(p_odd_out, p_odd_out + seg_size * num_dg_pts, p_out);
     146             :         p_odd_out += seg_size * num_dg_pts;  // NOLINT
     147             :       }
     148             :       p_out += seg_size * num_dg_pts;  // NOLINT
     149             :       is_even = not is_even;
     150             :     }
     151             :     return;
     152             :   }
     153             :   reconstruct_impl(dg_u, subcell_u, dg_mesh, subcell_extents,
     154             :                    reconstruction_method, Spectral::Parity::Uninitialized);
     155             : }
     156             : }  // namespace detail
     157             : 
     158             : /// @{
     159             : /*!
     160             :  * \ingroup DgSubcellGroup
     161             :  * \brief reconstruct the variable `subcell_u_times_projected_det_jac` onto the
     162             :  * DG grid `dg_mesh`.
     163             :  *
     164             :  * In general we wish that the reconstruction operator is the pseudo-inverse of
     165             :  * the projection operator. On curved meshes this means we either need to
     166             :  * compute a (time-dependent) reconstruction and projection matrix on each DG
     167             :  * element, or we expand the determinant of the Jacobian on the basis, accepting
     168             :  * the aliasing errors from that. We accept the aliasing errors in favor of the
     169             :  * significantly reduced computational overhead. This means that the projection
     170             :  * and reconstruction operators are only inverses of each other if both operate
     171             :  * on \f$u J\f$ where \f$u\f$ is the variable being projected and \f$J\f$ is the
     172             :  * determinant of the Jacobian. That is, the matrices are guaranteed to satisfy
     173             :  * \f$\mathcal{R}(\mathcal{P}(u J))=u J\f$. If the mesh is regular Cartesian,
     174             :  * then this isn't an issue. Furthermore, if we reconstruct
     175             :  * \f$uJ/\mathcal{P}(J)\f$ we again recover the exact DG solution. Doing the
     176             :  * latter has the advantage that, in general, we are ideally projecting to the
     177             :  * subcells much more often than reconstructing from them (a statement that we
     178             :  * would rather use DG more than the subcells).
     179             :  *
     180             :  * When the DG mesh uses a ZernikeB1 basis the Variables overloads deduce
     181             :  * per-component parity from the tag list automatically. The raw `DataVector`
     182             :  * overloads accepting a `Spectral::Parity` are for single-component data where
     183             :  * the caller already knows the parity. Only `DimByDim` reconstruction is
     184             :  * supported for ZernikeB1 meshes.
     185             :  */
     186             : template <size_t Dim>
     187           1 : DataVector reconstruct(
     188             :     const DataVector& subcell_u_times_projected_det_jac,
     189             :     const Mesh<Dim>& dg_mesh, const Index<Dim>& subcell_extents,
     190             :     ReconstructionMethod reconstruction_method,
     191             :     Spectral::Parity parity = Spectral::Parity::Uninitialized);
     192             : 
     193             : template <size_t Dim>
     194           1 : void reconstruct(gsl::not_null<DataVector*> dg_u,
     195             :                  const DataVector& subcell_u_times_projected_det_jac,
     196             :                  const Mesh<Dim>& dg_mesh, const Index<Dim>& subcell_extents,
     197             :                  ReconstructionMethod reconstruction_method,
     198             :                  Spectral::Parity parity = Spectral::Parity::Uninitialized);
     199             : 
     200             : template <typename SubcellTagList, typename DgTagList, size_t Dim>
     201           1 : void reconstruct(const gsl::not_null<Variables<DgTagList>*> dg_u,
     202             :                  const Variables<SubcellTagList>& subcell_u,
     203             :                  const Mesh<Dim>& dg_mesh, const Index<Dim>& subcell_extents,
     204             :                  const ReconstructionMethod reconstruction_method) {
     205             :   ASSERT(subcell_u.number_of_grid_points() == subcell_extents.product(),
     206             :          "Incorrect subcell size of u: " << subcell_u.number_of_grid_points()
     207             :                                          << " but should be "
     208             :                                          << subcell_extents.product());
     209             :   if (UNLIKELY(dg_u->number_of_grid_points() !=
     210             :                dg_mesh.number_of_grid_points())) {
     211             :     dg_u->initialize(dg_mesh.number_of_grid_points(), 0.0);
     212             :   }
     213             :   detail::reconstruct_impl_with_tag_list<DgTagList>(
     214             :       gsl::span<double>{dg_u->data(), dg_u->size()},
     215             :       gsl::span<const double>{subcell_u.data(), subcell_u.size()}, dg_mesh,
     216             :       subcell_extents, reconstruction_method);
     217             : }
     218             : 
     219             : template <typename TagList, size_t Dim>
     220           1 : Variables<TagList> reconstruct(
     221             :     const Variables<TagList>& subcell_u, const Mesh<Dim>& dg_mesh,
     222             :     const Index<Dim>& subcell_extents,
     223             :     const ReconstructionMethod reconstruction_method) {
     224             :   Variables<TagList> dg_u(dg_mesh.number_of_grid_points());
     225             :   reconstruct(make_not_null(&dg_u), subcell_u, dg_mesh, subcell_extents,
     226             :               reconstruction_method);
     227             :   return dg_u;
     228             : }
     229             : /// @}
     230             : }  // namespace evolution::dg::subcell::fd