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