Line data    Source code

       1           0 : // Distributed under the MIT License.
       2             : // See LICENSE.txt for details.
       3             : 
       4             : #pragma once
       5             : 
       6             : #include <cstddef>
       7             : 
       8             : #include "DataStructures/Variables.hpp"
       9             : #include "Evolution/DgSubcell/ReconstructionMethod.hpp"
      10             : #include "Utilities/ErrorHandling/Assert.hpp"
      11             : #include "Utilities/Gsl.hpp"
      12             : 
      13             : /// \cond
      14             : class DataVector;
      15             : template <size_t>
      16             : class Index;
      17             : template <size_t>
      18             : class Mesh;
      19             : /// \endcond
      20             : 
      21             : namespace evolution::dg::subcell::fd {
      22             : namespace detail {
      23             : template <size_t Dim>
      24             : void reconstruct_impl(gsl::span<double> dg_u,
      25             :                       gsl::span<const double> subcell_u_times_projected_det_jac,
      26             :                       const Mesh<Dim>& dg_mesh,
      27             :                       const Index<Dim>& subcell_extents,
      28             :                       ReconstructionMethod reconstruction_method);
      29             : }  // namespace detail
      30             : 
      31             : /// @{
      32             : /*!
      33             :  * \ingroup DgSubcellGroup
      34             :  * \brief reconstruct the variable `subcell_u_times_projected_det_jac` onto the
      35             :  * DG grid `dg_mesh`.
      36             :  *
      37             :  * In general we wish that the reconstruction operator is the pseudo-inverse of
      38             :  * the projection operator. On curved meshes this means we either need to
      39             :  * compute a (time-dependent) reconstruction and projection matrix on each DG
      40             :  * element, or we expand the determinant of the Jacobian on the basis, accepting
      41             :  * the aliasing errors from that. We accept the aliasing errors in favor of the
      42             :  * significantly reduced computational overhead. This means that the projection
      43             :  * and reconstruction operators are only inverses of each other if both operate
      44             :  * on \f$u J\f$ where \f$u\f$ is the variable being projected and \f$J\f$ is the
      45             :  * determinant of the Jacobian. That is, the matrices are guaranteed to satisfy
      46             :  * \f$\mathcal{R}(\mathcal{P}(u J))=u J\f$. If the mesh is regular Cartesian,
      47             :  * then this isn't an issue. Furthermore, if we reconstruct
      48             :  * \f$uJ/\mathcal{P}(J)\f$ we again recover the exact DG solution. Doing the
      49             :  * latter has the advantage that, in general, we are ideally projecting to the
      50             :  * subcells much more often than reconstructing from them (a statement that we
      51             :  * would rather use DG more than the subcells).
      52             :  */
      53             : template <size_t Dim>
      54           1 : DataVector reconstruct(const DataVector& subcell_u_times_projected_det_jac,
      55             :                        const Mesh<Dim>& dg_mesh,
      56             :                        const Index<Dim>& subcell_extents,
      57             :                        ReconstructionMethod reconstruction_method);
      58             : 
      59             : template <size_t Dim>
      60           1 : void reconstruct(gsl::not_null<DataVector*> dg_u,
      61             :                  const DataVector& subcell_u_times_projected_det_jac,
      62             :                  const Mesh<Dim>& dg_mesh, const Index<Dim>& subcell_extents,
      63             :                  ReconstructionMethod reconstruction_method);
      64             : 
      65             : template <typename SubcellTagList, typename DgTagList, size_t Dim>
      66           1 : void reconstruct(const gsl::not_null<Variables<DgTagList>*> dg_u,
      67             :                  const Variables<SubcellTagList>& subcell_u,
      68             :                  const Mesh<Dim>& dg_mesh, const Index<Dim>& subcell_extents,
      69             :                  const ReconstructionMethod reconstruction_method) {
      70             :   ASSERT(subcell_u.number_of_grid_points() == subcell_extents.product(),
      71             :          "Incorrect subcell size of u: " << subcell_u.number_of_grid_points()
      72             :                                          << " but should be "
      73             :                                          << subcell_extents.product());
      74             :   if (UNLIKELY(dg_u->number_of_grid_points() !=
      75             :                dg_mesh.number_of_grid_points())) {
      76             :     dg_u->initialize(dg_mesh.number_of_grid_points(), 0.0);
      77             :   }
      78             :   detail::reconstruct_impl(
      79             :       gsl::span<double>{dg_u->data(), dg_u->size()},
      80             :       gsl::span<const double>{subcell_u.data(), subcell_u.size()}, dg_mesh,
      81             :       subcell_extents, reconstruction_method);
      82             : }
      83             : 
      84             : template <typename TagList, size_t Dim>
      85           1 : Variables<TagList> reconstruct(
      86             :     const Variables<TagList>& subcell_u, const Mesh<Dim>& dg_mesh,
      87             :     const Index<Dim>& subcell_extents,
      88             :     const ReconstructionMethod reconstruction_method) {
      89             :   Variables<TagList> dg_u(dg_mesh.number_of_grid_points());
      90             :   reconstruct(make_not_null(&dg_u), subcell_u, dg_mesh, subcell_extents,
      91             :               reconstruction_method);
      92             :   return dg_u;
      93             : }
      94             : /// @}
      95             : }  // namespace evolution::dg::subcell::fd