PrimReconstructor.hpp
1 // Distributed under the MIT License.
2 // See LICENSE.txt for details.
3 
4 #pragma once
5 
7 
8 #include <array>
9 #include <cstddef>
10 #include <utility>
11 
12 #include "DataStructures/DataBox/PrefixHelpers.hpp"
14 #include "DataStructures/DataBox/TagName.hpp"
15 #include "DataStructures/DataVector.hpp"
16 #include "DataStructures/Index.hpp"
22 #include "Domain/Structure/DirectionMap.hpp"
25 #include "Domain/Structure/MaxNumberOfNeighbors.hpp"
27 #include "Evolution/DgSubcell/NeighborData.hpp"
28 #include "Evolution/DgSubcell/SliceData.hpp"
29 #include "Evolution/DiscontinuousGalerkin/Actions/NormalCovectorAndMagnitude.hpp"
30 #include "Evolution/Systems/GrMhd/ValenciaDivClean/ConservativeFromPrimitive.hpp"
31 #include "Evolution/Systems/GrMhd/ValenciaDivClean/FiniteDifference/Reconstructor.hpp"
32 #include "Evolution/Systems/GrMhd/ValenciaDivClean/System.hpp"
33 #include "Evolution/Systems/GrMhd/ValenciaDivClean/Tags.hpp"
36 #include "PointwiseFunctions/Hydro/EquationsOfState/EquationOfState.hpp"
37 #include "PointwiseFunctions/Hydro/EquationsOfState/IdealFluid.hpp"
38 #include "PointwiseFunctions/Hydro/EquationsOfState/PolytropicFluid.hpp"
39 #include "PointwiseFunctions/Hydro/SpecificEnthalpy.hpp"
40 #include "PointwiseFunctions/Hydro/Tags.hpp"
41 #include "Utilities/Gsl.hpp"
42 #include "Utilities/TMPL.hpp"
43 
44 namespace TestHelpers::grmhd::ValenciaDivClean::fd {
45 namespace detail {
46 template <typename F>
50  boost::hash<std::pair<Direction<3>, ElementId<3>>>>
51 compute_neighbor_data(
52  const Mesh<3>& subcell_mesh,
53  const tnsr::I<DataVector, 3, Frame::Logical>& volume_logical_coords,
54  const DirectionMap<3, Neighbors<3>>& neighbors,
55  const size_t ghost_zone_size, const F& compute_variables_of_neighbor_data) {
59  boost::hash<std::pair<Direction<3>, ElementId<3>>>>
60  neighbor_data{};
61  for (const auto& [direction, neighbors_in_direction] : neighbors) {
62  REQUIRE(neighbors_in_direction.size() == 1);
63  const ElementId<3>& neighbor_id = *neighbors_in_direction.begin();
64  auto neighbor_logical_coords = volume_logical_coords;
65  neighbor_logical_coords.get(direction.dimension()) +=
66  direction.sign() * 2.0;
67  const auto neighbor_vars_for_reconstruction =
68  compute_variables_of_neighbor_data(neighbor_logical_coords);
69 
70  DirectionMap<3, bool> directions_to_slice{};
71  directions_to_slice[direction.opposite()] = true;
72  const auto sliced_data = evolution::dg::subcell::detail::slice_data_impl(
73  gsl::make_span(neighbor_vars_for_reconstruction.data(),
74  neighbor_vars_for_reconstruction.size()),
75  subcell_mesh.extents(), ghost_zone_size, directions_to_slice);
76  REQUIRE(sliced_data.size() == 1);
77  REQUIRE(sliced_data.contains(direction.opposite()));
78  neighbor_data[std::pair{direction, neighbor_id}].data_for_reconstruction =
79  sliced_data.at(direction.opposite());
80  }
81  return neighbor_data;
82 }
83 
84 template <size_t ThermodynamicDim, typename Reconstructor>
85 void test_prim_reconstructor_impl(
86  const size_t points_per_dimension,
87  const Reconstructor& derived_reconstructor,
89  // 1. Create linear prims to reconstruct
90  // 2. send through reconstruction
91  // 3. check prims and cons were computed correctly
92  namespace mhd = ::grmhd::ValenciaDivClean;
93  const mhd::fd::Reconstructor& reconstructor = derived_reconstructor;
94  static_assert(
95  tmpl::list_contains_v<typename mhd::fd::Reconstructor::creatable_classes,
96  Reconstructor>);
97 
99  using Pressure = hydro::Tags::Pressure<DataVector>;
103  using VelocityW =
105  using SpecificInternalEnergy =
107  using SpecificEnthalpy = hydro::Tags::SpecificEnthalpy<DataVector>;
108  using LorentzFactor = hydro::Tags::LorentzFactor<DataVector>;
109 
110  using prims_tags = hydro::grmhd_tags<DataVector>;
111  using cons_tags = typename mhd::System::variables_tag::tags_list;
114  using prim_tags_for_reconstruction =
115  tmpl::list<Rho, Pressure, VelocityW, MagField, Phi>;
116 
117  const Mesh<3> subcell_mesh{points_per_dimension,
118  Spectral::Basis::FiniteDifference,
119  Spectral::Quadrature::CellCentered};
120  auto logical_coords = logical_coordinates(subcell_mesh);
121  // Make the logical coordinates different in each direction
122  for (size_t i = 1; i < 3; ++i) {
123  logical_coords.get(i) += 4.0 * i;
124  }
125 
126  DirectionMap<3, Neighbors<3>> neighbors{};
127  for (size_t i = 0; i < 2 * 3; ++i) {
128  neighbors[gsl::at(Direction<3>::all_directions(), i)] =
129  Neighbors<3>{{ElementId<3>{i + 1, {}}}, {}};
130  }
131  const Element<3> element{ElementId<3>{0, {}}, neighbors};
132  const auto compute_solution = [](const auto& coords) {
133  Variables<prim_tags_for_reconstruction> vars{get<0>(coords).size(), 0.0};
134  for (size_t i = 0; i < 3; ++i) {
135  get(get<Rho>(vars)) += coords.get(i);
136  get(get<Pressure>(vars)) += coords.get(i);
137  get(get<Phi>(vars)) += coords.get(i);
138  for (size_t j = 0; j < 3; ++j) {
139  get<VelocityW>(vars).get(j) += coords.get(i);
140  get<MagField>(vars).get(j) += coords.get(i);
141  }
142  }
143  get(get<Rho>(vars)) += 2.0;
144  get(get<Pressure>(vars)) += 30.0;
145  get(get<Phi>(vars)) += 50.0;
146  for (size_t j = 0; j < 3; ++j) {
147  get<VelocityW>(vars).get(j) += 1.0e-2 * (j + 2.0) + 10.0;
148  get<MagField>(vars).get(j) += 1.0e-2 * (j + 2.0) + 60.0;
149  }
150  return vars;
151  };
152 
156  boost::hash<std::pair<Direction<3>, ElementId<3>>>>
157  neighbor_data = compute_neighbor_data(
158  subcell_mesh, logical_coords, element.neighbors(),
159  reconstructor.ghost_zone_size(), compute_solution);
160 
161  const size_t reconstructed_num_pts =
162  (subcell_mesh.extents(0) + 1) *
163  subcell_mesh.extents().slice_away(0).product();
164 
165  using dg_package_data_argument_tags = tmpl::append<
166  cons_tags,
168  prims_tags,
170  flux_tags,
171  tmpl::remove_duplicates<tmpl::push_back<tmpl::list<
177  evolution::dg::Actions::detail::NormalVector<3>>>>>;
178  tnsr::ii<DataVector, 3, Frame::Inertial> lower_face_spatial_metric{
179  reconstructed_num_pts, 0.0};
180  tnsr::ii<DataVector, 3, Frame::Inertial> upper_face_spatial_metric{
181  reconstructed_num_pts, 0.0};
182  for (size_t i = 0; i < 3; ++i) {
183  lower_face_spatial_metric.get(i, i) = 1.0 + 0.01 * i;
184  upper_face_spatial_metric.get(i, i) = 1.0 - 0.01 * i;
185  }
186  const Scalar<DataVector> lower_face_sqrt_det_spatial_metric{
187  sqrt(get(determinant(lower_face_spatial_metric)))};
188  const Scalar<DataVector> upper_face_sqrt_det_spatial_metric{
189  sqrt(get(determinant(upper_face_spatial_metric)))};
190 
192  make_array<3>(
193  Variables<dg_package_data_argument_tags>(reconstructed_num_pts));
195  make_array<3>(
196  Variables<dg_package_data_argument_tags>(reconstructed_num_pts));
197  for (size_t i = 0; i < 3; ++i) {
198  get<gr::Tags::SqrtDetSpatialMetric<DataVector>>(
199  gsl::at(vars_on_lower_face, i)) = lower_face_sqrt_det_spatial_metric;
200  get<gr::Tags::SqrtDetSpatialMetric<DataVector>>(
201  gsl::at(vars_on_upper_face, i)) = upper_face_sqrt_det_spatial_metric;
202 
203  get<gr::Tags::SpatialMetric<3>>(gsl::at(vars_on_lower_face, i)) =
204  lower_face_spatial_metric;
205  get<gr::Tags::SpatialMetric<3>>(gsl::at(vars_on_upper_face, i)) =
206  upper_face_spatial_metric;
207  }
208 
209  Variables<prims_tags> volume_prims{subcell_mesh.number_of_grid_points()};
210  {
211  const auto volume_prims_for_recons = compute_solution(logical_coords);
212  tmpl::for_each<tmpl::list<Rho, Pressure, MagField, Phi>>(
213  [&volume_prims, &volume_prims_for_recons](auto tag_v) {
214  using tag = tmpl::type_from<decltype(tag_v)>;
215  get<tag>(volume_prims) = get<tag>(volume_prims_for_recons);
216  });
217  // Metric is identity at cell center
218  get(get<LorentzFactor>(volume_prims)) =
219  sqrt(1.0 + get(dot_product(get<VelocityW>(volume_prims_for_recons),
220  get<VelocityW>(volume_prims_for_recons))));
221  for (size_t i = 0; i < 3; ++i) {
222  get<Velocity>(volume_prims).get(i) =
223  get<VelocityW>(volume_prims_for_recons).get(i) /
224  get(get<LorentzFactor>(volume_prims));
225  }
226  }
227 
228  // Now we have everything to call the reconstruction
229  dynamic_cast<const Reconstructor&>(reconstructor)
230  .reconstruct(make_not_null(&vars_on_lower_face),
231  make_not_null(&vars_on_upper_face), volume_prims, eos,
232  element, neighbor_data, subcell_mesh);
233 
234  for (size_t dim = 0; dim < 3; ++dim) {
235  CAPTURE(dim);
236  const auto basis = make_array<3>(Spectral::Basis::FiniteDifference);
237  auto quadrature = make_array<3>(Spectral::Quadrature::CellCentered);
238  auto extents = make_array<3>(points_per_dimension);
239  gsl::at(extents, dim) = points_per_dimension + 1;
240  gsl::at(quadrature, dim) = Spectral::Quadrature::FaceCentered;
241  const Mesh<3> face_centered_mesh{extents, basis, quadrature};
242  auto logical_coords_face_centered = logical_coordinates(face_centered_mesh);
243  for (size_t i = 1; i < 3; ++i) {
244  logical_coords_face_centered.get(i) =
245  logical_coords_face_centered.get(i) + 4.0 * i;
246  }
247  Variables<dg_package_data_argument_tags> expected_lower_face_values{
248  face_centered_mesh.number_of_grid_points()};
249  expected_lower_face_values.assign_subset(
250  compute_solution(logical_coords_face_centered));
251  if constexpr (ThermodynamicDim == 2) {
252  get<SpecificInternalEnergy>(expected_lower_face_values) =
253  eos.specific_internal_energy_from_density_and_pressure(
254  get<Rho>(expected_lower_face_values),
255  get<Pressure>(expected_lower_face_values));
256  } else {
257  get<SpecificInternalEnergy>(expected_lower_face_values) =
258  eos.specific_internal_energy_from_density(
259  get<Rho>(expected_lower_face_values));
260  }
261  get<SpecificEnthalpy>(expected_lower_face_values) =
263  get<Rho>(expected_lower_face_values),
264  get<SpecificInternalEnergy>(expected_lower_face_values),
265  get<Pressure>(expected_lower_face_values));
266 
267  Variables<dg_package_data_argument_tags> expected_upper_face_values =
268  expected_lower_face_values;
269  get(get<LorentzFactor>(expected_lower_face_values)) =
270  sqrt(1.0 + get(dot_product(get<VelocityW>(expected_lower_face_values),
271  get<VelocityW>(expected_lower_face_values),
272  lower_face_spatial_metric)));
273  get(get<LorentzFactor>(expected_upper_face_values)) =
274  sqrt(1.0 + get(dot_product(get<VelocityW>(expected_upper_face_values),
275  get<VelocityW>(expected_upper_face_values),
276  upper_face_spatial_metric)));
277  for (size_t i = 0; i < 3; ++i) {
278  get<Velocity>(expected_lower_face_values).get(i) =
279  get<VelocityW>(expected_lower_face_values).get(i) /
280  get(get<LorentzFactor>(expected_lower_face_values));
281  get<Velocity>(expected_upper_face_values).get(i) =
282  get<VelocityW>(expected_upper_face_values).get(i) /
283  get(get<LorentzFactor>(expected_upper_face_values));
284  }
285 
286  get<gr::Tags::SqrtDetSpatialMetric<DataVector>>(
287  expected_lower_face_values) = lower_face_sqrt_det_spatial_metric;
288  get<gr::Tags::SpatialMetric<3>>(expected_lower_face_values) =
289  lower_face_spatial_metric;
290  get<gr::Tags::SqrtDetSpatialMetric<DataVector>>(
291  expected_upper_face_values) = upper_face_sqrt_det_spatial_metric;
292  get<gr::Tags::SpatialMetric<3>>(expected_upper_face_values) =
293  upper_face_spatial_metric;
294 
295  mhd::ConservativeFromPrimitive::apply(
296  make_not_null(&get<mhd::Tags::TildeD>(expected_lower_face_values)),
297  make_not_null(&get<mhd::Tags::TildeTau>(expected_lower_face_values)),
298  make_not_null(&get<mhd::Tags::TildeS<>>(expected_lower_face_values)),
299  make_not_null(&get<mhd::Tags::TildeB<>>(expected_lower_face_values)),
300  make_not_null(&get<mhd::Tags::TildePhi>(expected_lower_face_values)),
301  get<Rho>(expected_lower_face_values),
302  get<SpecificInternalEnergy>(expected_lower_face_values),
303  get<SpecificEnthalpy>(expected_lower_face_values),
304  get<Pressure>(expected_lower_face_values),
305  get<Velocity>(expected_lower_face_values),
306  get<LorentzFactor>(expected_lower_face_values),
307  get<MagField>(expected_lower_face_values),
309  expected_lower_face_values),
310  get<gr::Tags::SpatialMetric<3>>(expected_lower_face_values),
311  get<Phi>(expected_lower_face_values));
312  mhd::ConservativeFromPrimitive::apply(
313  make_not_null(&get<mhd::Tags::TildeD>(expected_upper_face_values)),
314  make_not_null(&get<mhd::Tags::TildeTau>(expected_upper_face_values)),
315  make_not_null(&get<mhd::Tags::TildeS<>>(expected_upper_face_values)),
316  make_not_null(&get<mhd::Tags::TildeB<>>(expected_upper_face_values)),
317  make_not_null(&get<mhd::Tags::TildePhi>(expected_upper_face_values)),
318  get<Rho>(expected_upper_face_values),
319  get<SpecificInternalEnergy>(expected_upper_face_values),
320  get<SpecificEnthalpy>(expected_upper_face_values),
321  get<Pressure>(expected_upper_face_values),
322  get<Velocity>(expected_upper_face_values),
323  get<LorentzFactor>(expected_upper_face_values),
324  get<MagField>(expected_upper_face_values),
326  expected_upper_face_values),
327  get<gr::Tags::SpatialMetric<3>>(expected_upper_face_values),
328  get<Phi>(expected_upper_face_values));
329 
330  tmpl::for_each<tmpl::append<cons_tags, prims_tags>>(
331  [dim, &expected_lower_face_values, &expected_upper_face_values,
332  &vars_on_lower_face, &vars_on_upper_face](auto tag_to_check_v) {
333  using tag_to_check = tmpl::type_from<decltype(tag_to_check_v)>;
334  CAPTURE(db::tag_name<tag_to_check>());
336  get<tag_to_check>(gsl::at(vars_on_lower_face, dim)),
337  get<tag_to_check>(expected_lower_face_values));
339  get<tag_to_check>(gsl::at(vars_on_upper_face, dim)),
340  get<tag_to_check>(expected_upper_face_values));
341  });
342  }
343 }
344 } // namespace detail
345 
346 template <typename Reconstructor>
347 void test_prim_reconstructor(const size_t points_per_dimension,
348  const Reconstructor& derived_reconstructor) {
349  detail::test_prim_reconstructor_impl(points_per_dimension,
350  derived_reconstructor,
352  detail::test_prim_reconstructor_impl(
353  points_per_dimension, derived_reconstructor,
355 }
356 } // namespace TestHelpers::grmhd::ValenciaDivClean::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
hydro::Tags::Pressure
The fluid pressure .
Definition: Tags.hpp:119
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
hydro::Tags::SpecificInternalEnergy
The specific internal energy .
Definition: Tags.hpp:173
gr::Tags::SpatialMetric
Definition: Tags.hpp:26
sqrt
auto sqrt(const TensorExpression< T, typename T::type, tmpl::list<>, tmpl::list<>, tmpl::list<>> &t)
Returns the tensor expression representing the square root of a tensor expression that evaluates to a...
Definition: SquareRoot.hpp:83
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
hydro::Tags::DivergenceCleaningField
The divergence-cleaning field .
Definition: Tags.hpp:48
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
determinant
void determinant(const gsl::not_null< Scalar< T > * > det_tensor, const Tensor< T, Symm, index_list< Index0, Index1 >> &tensor) noexcept
Computes the determinant of a rank-2 Tensor tensor.
Definition: Determinant.hpp:139
Spectral.hpp
Direction
Definition: Direction.hpp:23
Element
Definition: Element.hpp:29
dot_product
void dot_product(const gsl::not_null< Scalar< DataType > * > dot_product, const Tensor< DataType, Symmetry< 1 >, index_list< Index >> &vector_a, const Tensor< DataType, Symmetry< 1 >, index_list< Index >> &vector_b) noexcept
Compute the Euclidean dot product of two vectors or one forms.
Definition: DotProduct.hpp:25
domain::push_back
CoordinateMap< SourceFrame, TargetFrame, Maps..., NewMap > push_back(CoordinateMap< SourceFrame, TargetFrame, Maps... > old_map, NewMap new_map) noexcept
Creates a CoordinateMap by appending the new map to the end of the old maps.
ElementId
An ElementId uniquely labels an Element.
Definition: ElementId.hpp:51
ElementId.hpp
DotProduct.hpp
hydro::Tags::LorentzFactor
The Lorentz factor , where is the spatial velocity of the fluid.
Definition: Tags.hpp:64
cstddef
Neighbors.hpp
EquationsOfState::IdealFluid< true >
hydro::Tags::SpecificEnthalpy
The relativistic specific enthalpy .
Definition: Tags.hpp:167
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< true >
hydro::relativistic_specific_enthalpy
void relativistic_specific_enthalpy(gsl::not_null< Scalar< DataType > * > result, const Scalar< DataType > &rest_mass_density, const Scalar< DataType > &specific_internal_energy, const Scalar< DataType > &pressure) noexcept
Computes the relativistic specific enthalpy as: where is the specific internal energy,...
gr::Tags::Shift
Definition: Tags.hpp:48
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
hydro::Tags::LorentzFactorTimesSpatialVelocity
The Lorentz factor times the spatial velocity .
Definition: Tags.hpp:185
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.
Scalar
Tensor< T, Symmetry<>, index_list<> > Scalar
Definition: TypeAliases.hpp:21
Gsl.hpp
hydro::Tags::SpatialVelocity
The spatial velocity of the fluid, where . Here is the spatial part of the 4-velocity of the fluid,...
Definition: Tags.hpp:142
hydro::Tags::MagneticField
The magnetic field measured by an Eulerian observer, where is the normal to the spatial hypersurfac...
Definition: Tags.hpp:80
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
Determinant.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
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
gr::Tags::Lapse
Definition: Tags.hpp:52
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.
gr::Tags::InverseSpatialMetric
Inverse of the spatial metric.
Definition: Tags.hpp:33
Mesh.hpp
hydro::Tags::RestMassDensity
The rest-mass density .
Definition: Tags.hpp:125
gr::Tags::SqrtDetSpatialMetric
Definition: Tags.hpp:44