PrimReconstructor.hpp
1 // Distributed under the MIT License.
2 // See LICENSE.txt for details.
3 
4 #pragma once
5 
7 
8 #include <cstddef>
9 #include <utility>
10 
11 #include "DataStructures/DataBox/PrefixHelpers.hpp"
13 #include "DataStructures/DataBox/TagName.hpp"
14 #include "DataStructures/DataVector.hpp"
15 #include "DataStructures/Index.hpp"
19 #include "Domain/Structure/DirectionMap.hpp"
22 #include "Domain/Structure/MaxNumberOfNeighbors.hpp"
24 #include "Evolution/DgSubcell/NeighborData.hpp"
25 #include "Evolution/DgSubcell/SliceData.hpp"
26 #include "Evolution/Systems/NewtonianEuler/ConservativeFromPrimitive.hpp"
27 #include "Evolution/Systems/NewtonianEuler/FiniteDifference/Reconstructor.hpp"
28 #include "Evolution/Systems/NewtonianEuler/Tags.hpp"
31 #include "PointwiseFunctions/Hydro/EquationsOfState/EquationOfState.hpp"
32 #include "PointwiseFunctions/Hydro/EquationsOfState/IdealFluid.hpp"
33 #include "PointwiseFunctions/Hydro/EquationsOfState/PolytropicFluid.hpp"
34 #include "Utilities/Gsl.hpp"
35 #include "Utilities/TMPL.hpp"
36 
37 namespace TestHelpers::NewtonianEuler::fd {
38 template <size_t Dim, typename F>
42  boost::hash<std::pair<Direction<Dim>, ElementId<Dim>>>>
43 compute_neighbor_data(
44  const Mesh<Dim>& subcell_mesh,
45  const tnsr::I<DataVector, Dim, Frame::Logical>& volume_logical_coords,
46  const DirectionMap<Dim, Neighbors<Dim>>& neighbors,
47  const size_t ghost_zone_size, const F& compute_variables_of_neighbor_data) {
51  boost::hash<std::pair<Direction<Dim>, ElementId<Dim>>>>
52  neighbor_data{};
53  for (const auto& [direction, neighbors_in_direction] : neighbors) {
54  REQUIRE(neighbors_in_direction.size() == 1);
55  const ElementId<Dim>& neighbor_id = *neighbors_in_direction.begin();
56  auto neighbor_logical_coords = volume_logical_coords;
57  neighbor_logical_coords.get(direction.dimension()) +=
58  direction.sign() * 2.0;
59  const auto neighbor_vars_for_reconstruction =
60  compute_variables_of_neighbor_data(neighbor_logical_coords);
61 
62  DirectionMap<Dim, bool> directions_to_slice{};
63  directions_to_slice[direction.opposite()] = true;
64  const auto sliced_data = evolution::dg::subcell::detail::slice_data_impl(
65  gsl::make_span(neighbor_vars_for_reconstruction.data(),
66  neighbor_vars_for_reconstruction.size()),
67  subcell_mesh.extents(), ghost_zone_size, directions_to_slice);
68  REQUIRE(sliced_data.size() == 1);
69  REQUIRE(sliced_data.contains(direction.opposite()));
70  neighbor_data[std::pair{direction, neighbor_id}].data_for_reconstruction =
71  sliced_data.at(direction.opposite());
72  }
73  return neighbor_data;
74 }
75 
76 namespace detail {
77 template <size_t Dim, size_t ThermodynamicDim, typename Reconstructor>
78 void test_prim_reconstructor_impl(
79  const size_t points_per_dimension,
80  const Reconstructor& derived_reconstructor,
82  // 1. Create linear prims to reconstruct
83  // 2. send through reconstruction
84  // 3. check prims and cons were computed correctly
85  namespace ne = ::NewtonianEuler;
86  const ne::fd::Reconstructor<Dim>& reconstructor = derived_reconstructor;
87  static_assert(tmpl::list_contains_v<
88  typename ne::fd::Reconstructor<Dim>::creatable_classes,
89  Reconstructor>);
90 
91  using MassDensityCons = ne::Tags::MassDensityCons;
92  using EnergyDensity = ne::Tags::EnergyDensity;
93  using MomentumDensity = ne::Tags::MomentumDensity<Dim>;
94 
95  // Primitive vars tags
96  using MassDensity = ne::Tags::MassDensity<DataVector>;
97  using Velocity = ne::Tags::Velocity<DataVector, Dim>;
98  using SpecificInternalEnergy = ne::Tags::SpecificInternalEnergy<DataVector>;
99  using Pressure = ne::Tags::Pressure<DataVector>;
100 
101  using prims_tags =
102  tmpl::list<MassDensity, Velocity, SpecificInternalEnergy, Pressure>;
103  using cons_tags = tmpl::list<MassDensityCons, MomentumDensity, EnergyDensity>;
106  using prim_tags_for_reconstruction =
107  tmpl::list<MassDensity, Velocity, Pressure>;
108 
109  const Mesh<Dim> subcell_mesh{points_per_dimension,
110  Spectral::Basis::FiniteDifference,
111  Spectral::Quadrature::CellCentered};
112  auto logical_coords = logical_coordinates(subcell_mesh);
113  // Make the logical coordinates different in each direction
114  for (size_t i = 1; i < Dim; ++i) {
115  logical_coords.get(i) += 4.0 * i;
116  }
117 
119  for (size_t i = 0; i < 2 * Dim; ++i) {
120  neighbors[gsl::at(Direction<Dim>::all_directions(), i)] =
121  Neighbors<Dim>{{ElementId<Dim>{i + 1, {}}}, {}};
122  }
123  const Element<Dim> element{ElementId<Dim>{0, {}}, neighbors};
124  const auto compute_solution = [](const auto& coords) {
125  Variables<prim_tags_for_reconstruction> vars{get<0>(coords).size(), 0.0};
126  for (size_t i = 0; i < Dim; ++i) {
127  get(get<MassDensity>(vars)) += coords.get(i);
128  get(get<Pressure>(vars)) += coords.get(i);
129  for (size_t j = 0;j <Dim;++j) {
130  get<Velocity>(vars).get(j) += coords.get(i);
131  }
132  }
133  get(get<MassDensity>(vars)) += 2.0;
134  get(get<Pressure>(vars)) += 30.0;
135  for (size_t j = 0; j < Dim; ++j) {
136  get<Velocity>(vars).get(j) += 1.0e-2 * (j + 2.0) + 10.0;
137  }
138  return vars;
139  };
140 
144  boost::hash<std::pair<Direction<Dim>, ElementId<Dim>>>>
145  neighbor_data = compute_neighbor_data(
146  subcell_mesh, logical_coords, element.neighbors(),
147  reconstructor.ghost_zone_size(), compute_solution);
148 
149  const size_t reconstructed_num_pts =
150  (subcell_mesh.extents(0) + 1) *
151  subcell_mesh.extents().slice_away(0).product();
152 
153  using dg_package_data_argument_tags =
154  tmpl::append<cons_tags, prims_tags, flux_tags>;
155  std::array<Variables<dg_package_data_argument_tags>, Dim> vars_on_lower_face =
156  make_array<Dim>(
157  Variables<dg_package_data_argument_tags>(reconstructed_num_pts));
158  std::array<Variables<dg_package_data_argument_tags>, Dim> vars_on_upper_face =
159  make_array<Dim>(
160  Variables<dg_package_data_argument_tags>(reconstructed_num_pts));
161 
162  Variables<prims_tags> volume_prims{subcell_mesh.number_of_grid_points()};
163  volume_prims.assign_subset(compute_solution(logical_coords));
164 
165  // Now we have everything to call the reconstruction
166  dynamic_cast<const Reconstructor&>(reconstructor).reconstruct(
167  make_not_null(&vars_on_lower_face), make_not_null(&vars_on_upper_face),
168  volume_prims, eos, element, neighbor_data, subcell_mesh);
169 
170  for (size_t dim = 0; dim < Dim; ++dim) {
171  CAPTURE(dim);
172  const auto basis = make_array<Dim>(Spectral::Basis::FiniteDifference);
173  auto quadrature = make_array<Dim>(Spectral::Quadrature::CellCentered);
174  auto extents = make_array<Dim>(points_per_dimension);
175  gsl::at(extents, dim) = points_per_dimension + 1;
176  gsl::at(quadrature, dim) = Spectral::Quadrature::FaceCentered;
177  const Mesh<Dim> face_centered_mesh{extents, basis, quadrature};
178  auto logical_coords_face_centered = logical_coordinates(face_centered_mesh);
179  for (size_t i = 1; i < Dim; ++i) {
180  logical_coords_face_centered.get(i) =
181  logical_coords_face_centered.get(i) + 4.0 * i;
182  }
183  Variables<dg_package_data_argument_tags> expected_face_values{
184  face_centered_mesh.number_of_grid_points()};
185  expected_face_values.assign_subset(
186  compute_solution(logical_coords_face_centered));
187  if constexpr (ThermodynamicDim == 2) {
188  get<SpecificInternalEnergy>(expected_face_values) =
189  eos.specific_internal_energy_from_density_and_pressure(
190  get<MassDensity>(expected_face_values),
191  get<Pressure>(expected_face_values));
192  } else {
193  get<SpecificInternalEnergy>(expected_face_values) =
194  eos.specific_internal_energy_from_density(
195  get<MassDensity>(expected_face_values));
196  }
197  ne::ConservativeFromPrimitive<Dim>::apply(
198  make_not_null(&get<MassDensityCons>(expected_face_values)),
199  make_not_null(&get<MomentumDensity>(expected_face_values)),
200  make_not_null(&get<EnergyDensity>(expected_face_values)),
201  get<MassDensity>(expected_face_values),
202  get<Velocity>(expected_face_values),
203  get<SpecificInternalEnergy>(expected_face_values));
204 
205  tmpl::for_each<tmpl::append<cons_tags, prims_tags>>(
206  [dim, &expected_face_values, &vars_on_lower_face,
207  &vars_on_upper_face](auto tag_to_check_v) {
208  using tag_to_check = tmpl::type_from<decltype(tag_to_check_v)>;
209  CAPTURE(db::tag_name<tag_to_check>());
211  get<tag_to_check>(gsl::at(vars_on_lower_face, dim)),
212  get<tag_to_check>(expected_face_values));
214  get<tag_to_check>(gsl::at(vars_on_upper_face, dim)),
215  get<tag_to_check>(expected_face_values));
216  });
217  }
218 }
219 } // namespace detail
220 
221 template <size_t Dim, typename Reconstructor>
222 void test_prim_reconstructor(const size_t points_per_dimension,
223  const Reconstructor& derived_reconstructor) {
224  detail::test_prim_reconstructor_impl<Dim>(
225  points_per_dimension, derived_reconstructor,
227  detail::test_prim_reconstructor_impl<Dim>(
228  points_per_dimension, derived_reconstructor,
230 }
231 } // namespace TestHelpers::NewtonianEuler::fd
gsl::at
constexpr T & at(std::array< T, N > &arr, Size index)
Retrieve a entry from a container, with checks in Debug mode that the index being retrieved is valid.
Definition: Gsl.hpp:125
maximum_number_of_neighbors
constexpr size_t maximum_number_of_neighbors(const size_t dim)
Definition: MaxNumberOfNeighbors.hpp:15
utility
Frame::Inertial
Definition: IndexType.hpp:44
EquationsOfState::EquationOfState
Base class for equations of state depending on whether or not the system is relativistic,...
Definition: EquationOfState.hpp:63
std::pair
TestingFramework.hpp
db::get
const auto & get(const DataBox< TagList > &box) noexcept
Retrieve the item with tag Tag from the DataBox.
Definition: DataBox.hpp:791
CHECK_ITERABLE_APPROX
#define CHECK_ITERABLE_APPROX(a, b)
A wrapper around Catch's CHECK macro that checks approximate equality of entries in iterable containe...
Definition: TestingFramework.hpp:122
Spectral.hpp
Direction
Definition: Direction.hpp:23
Element
Definition: Element.hpp:29
ElementId< Dim >
ElementId.hpp
cstddef
Neighbors.hpp
EquationsOfState::IdealFluid< false >
std::array
LogicalCoordinates.hpp
Mesh::number_of_grid_points
size_t number_of_grid_points() const noexcept
The total number of grid points in all dimensions.
Definition: Mesh.hpp:165
Index.hpp
EquationsOfState::PolytropicFluid< false >
Neighbors
Definition: Neighbors.hpp:28
Element.hpp
Mesh
Holds the number of grid points, basis, and quadrature in each direction of the computational grid.
Definition: Mesh.hpp:49
evolution::dg::subcell::NeighborData
Holds neighbor data needed for the TCI and reconstruction.
Definition: NeighborData.hpp:24
DirectionMap
Definition: DirectionMap.hpp:15
gsl::make_span
constexpr span< ElementType > make_span(ElementType *ptr, typename span< ElementType >::index_type count)
Definition: Gsl.hpp:801
logical_coordinates
void logical_coordinates(gsl::not_null< tnsr::I< DataVector, VolumeDim, Frame::Logical > * > logical_coords, const Mesh< VolumeDim > &mesh) noexcept
Compute the logical coordinates in an Element.
Gsl.hpp
Tensor.hpp
db::wrap_tags_in
tmpl::transform< TagList, tmpl::bind< Wrapper, tmpl::_1, tmpl::pin< Args >... > > wrap_tags_in
Create a new tmpl::list of tags by wrapping each tag in TagList in Wrapper<_, Args....
Definition: PrefixHelpers.hpp:30
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
Direction.hpp
FixedHashMap
A hash table with a compile-time specified maximum size and ability to efficiently handle perfect has...
Definition: FixedHashMap.hpp:82
Mesh::extents
const Index< Dim > & extents() const noexcept
The number of grid points in each dimension of the grid.
Definition: Mesh.hpp:148
Prefixes.hpp
TMPL.hpp
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.
Mesh.hpp