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