TimeDerivative.hpp
1 // Distributed under the MIT License.
2 // See LICENSE.txt for details.
3 
4 #pragma once
5 
6 #include <array>
7 #include <cstddef>
8 #include <optional>
9 #include <type_traits>
10 
12 #include "DataStructures/DataBox/PrefixHelpers.hpp"
14 #include "DataStructures/DataVector.hpp"
18 #include "Domain/Tags.hpp"
19 #include "Evolution/BoundaryCorrectionTags.hpp"
20 #include "Evolution/DgSubcell/CartesianFluxDivergence.hpp"
21 #include "Evolution/DgSubcell/ComputeBoundaryTerms.hpp"
22 #include "Evolution/DgSubcell/CorrectPackagedData.hpp"
23 #include "Evolution/DgSubcell/Tags/Coordinates.hpp"
24 #include "Evolution/DgSubcell/Tags/Mesh.hpp"
25 #include "Evolution/DgSubcell/Tags/OnSubcellFaces.hpp"
26 #include "Evolution/DiscontinuousGalerkin/Actions/NormalCovectorAndMagnitude.hpp"
27 #include "Evolution/DiscontinuousGalerkin/Actions/PackageDataImpl.hpp"
28 #include "Evolution/DiscontinuousGalerkin/MortarTags.hpp"
29 #include "Evolution/Systems/GrMhd/ValenciaDivClean/BoundaryCorrections/BoundaryCorrection.hpp"
30 #include "Evolution/Systems/GrMhd/ValenciaDivClean/FiniteDifference/Reconstructor.hpp"
31 #include "Evolution/Systems/GrMhd/ValenciaDivClean/FiniteDifference/Tag.hpp"
32 #include "Evolution/Systems/GrMhd/ValenciaDivClean/Fluxes.hpp"
33 #include "Evolution/Systems/GrMhd/ValenciaDivClean/Sources.hpp"
34 #include "Evolution/Systems/GrMhd/ValenciaDivClean/Subcell/ComputeFluxes.hpp"
35 #include "Evolution/Systems/GrMhd/ValenciaDivClean/System.hpp"
37 #include "PointwiseFunctions/Hydro/Tags.hpp"
39 #include "Utilities/FakeVirtual.hpp"
40 #include "Utilities/Gsl.hpp"
41 #include "Utilities/TMPL.hpp"
42 
44 /*!
45  * \brief Compute the time derivative on the subcell grid using FD
46  * reconstruction.
47  *
48  * The code makes the following unchecked assumptions:
49  * - Assumes Cartesian coordinates with a diagonal Jacobian matrix
50  * from the logical to the inertial frame
51  * - Assumes the mesh is not moving (grid and inertial frame are the same)
52  */
54  template <typename DbTagsList>
55  static void apply(
56  const gsl::not_null<db::DataBox<DbTagsList>*> box,
57  const InverseJacobian<DataVector, 3, Frame::Logical, Frame::Grid>&
58  cell_centered_logical_to_grid_inv_jacobian,
59  const Scalar<DataVector>& /*cell_centered_det_inv_jacobian*/) noexcept {
60  using evolved_vars_tag = typename System::variables_tag;
61  using evolved_vars_tags = typename evolved_vars_tag::tags_list;
62  using prim_tags = typename System::primitive_variables_tag::tags_list;
63  using recons_prim_tags = tmpl::push_back<
64  prim_tags,
66  using fluxes_tags = db::wrap_tags_in<::Tags::Flux, evolved_vars_tags,
67  tmpl::size_t<3>, Frame::Inertial>;
68 
69  // The copy of Mesh is intentional to avoid a GCC-7 internal compiler error.
70  const Mesh<3> subcell_mesh =
71  db::get<evolution::dg::subcell::Tags::Mesh<3>>(*box);
72  ASSERT(
73  subcell_mesh == Mesh<3>(subcell_mesh.extents(0), subcell_mesh.basis(0),
74  subcell_mesh.quadrature(0)),
75  "The subcell/FD mesh must be isotropic for the FD time derivative but "
76  "got "
77  << subcell_mesh);
78  const size_t num_pts = subcell_mesh.number_of_grid_points();
79  const size_t reconstructed_num_pts =
80  (subcell_mesh.extents(0) + 1) *
81  subcell_mesh.extents().slice_away(0).product();
82 
83  const tnsr::I<DataVector, 3, Frame::Logical>& cell_centered_logical_coords =
84  db::get<evolution::dg::subcell::Tags::Coordinates<3, Frame::Logical>>(
85  *box);
86  std::array<double, 3> one_over_delta_xi{};
87  for (size_t i = 0; i < 3; ++i) {
88  // Note: assumes isotropic extents
89  gsl::at(one_over_delta_xi, i) =
90  1.0 / (get<0>(cell_centered_logical_coords)[1] -
91  get<0>(cell_centered_logical_coords)[0]);
92  }
93 
95  db::get<grmhd::ValenciaDivClean::fd::Tags::Reconstructor>(*box);
96 
97  const Element<3>& element = db::get<domain::Tags::Element<3>>(*box);
98  ASSERT(element.external_boundaries().size() == 0,
99  "Can't have external boundaries right now with subcell. ElementID "
100  << element.id());
101 
102  // Now package the data and compute the correction
103  const auto& boundary_correction =
104  db::get<evolution::Tags::BoundaryCorrection<System>>(*box);
105  using derived_boundary_corrections =
106  typename std::decay_t<decltype(boundary_correction)>::creatable_classes;
107  std::array<Variables<evolved_vars_tags>, 3> boundary_corrections{};
108  tmpl::for_each<
109  derived_boundary_corrections>([&](auto derived_correction_v) noexcept {
110  using DerivedCorrection = tmpl::type_from<decltype(derived_correction_v)>;
111  if (typeid(boundary_correction) == typeid(DerivedCorrection)) {
112  using dg_package_data_temporary_tags =
113  typename DerivedCorrection::dg_package_data_temporary_tags;
114  using dg_package_data_argument_tags = tmpl::append<
115  evolved_vars_tags, recons_prim_tags, fluxes_tags,
116  tmpl::remove_duplicates<tmpl::push_back<
117  dg_package_data_temporary_tags, gr::Tags::SpatialMetric<3>,
120  evolution::dg::Actions::detail::NormalVector<3>>>>;
121  // Computed prims and cons on face via reconstruction
122  auto vars_on_lower_face = make_array<3>(
123  Variables<dg_package_data_argument_tags>(reconstructed_num_pts));
124  auto vars_on_upper_face = make_array<3>(
125  Variables<dg_package_data_argument_tags>(reconstructed_num_pts));
126  // Copy over the face values of the metric quantities.
127  using spacetime_vars_to_copy = tmpl::list<
133  tmpl::for_each<spacetime_vars_to_copy>(
134  [&vars_on_lower_face, &vars_on_upper_face,
135  &spacetime_vars_on_face =
137  typename System::flux_spacetime_variables_tag, 3>>(*box)](
138  auto tag_v) noexcept {
139  using tag = tmpl::type_from<decltype(tag_v)>;
140  for (size_t d = 0; d < 3; ++d) {
141  get<tag>(gsl::at(vars_on_lower_face, d)) =
142  get<tag>(gsl::at(spacetime_vars_on_face, d));
143  get<tag>(gsl::at(vars_on_upper_face, d)) =
144  get<tag>(gsl::at(spacetime_vars_on_face, d));
145  }
146  });
147 
148  // Reconstruct data to the face
149  call_with_dynamic_type<void, typename grmhd::ValenciaDivClean::fd::
150  Reconstructor::creatable_classes>(
151  &recons, [&box, &vars_on_lower_face,
152  &vars_on_upper_face](const auto& reconstructor) noexcept {
153  db::apply<typename std::decay_t<decltype(
154  *reconstructor)>::reconstruction_argument_tags>(
155  [&vars_on_lower_face, &vars_on_upper_face,
156  &reconstructor](const auto&... args) noexcept {
157  reconstructor->reconstruct(
158  make_not_null(&vars_on_lower_face),
159  make_not_null(&vars_on_upper_face), args...);
160  },
161  *box);
162  });
163 
164  using dg_package_field_tags =
165  typename DerivedCorrection::dg_package_field_tags;
166  // Allocated outside for loop to reduce allocations
167  Variables<dg_package_field_tags> upper_packaged_data{
168  reconstructed_num_pts};
169  Variables<dg_package_field_tags> lower_packaged_data{
170  reconstructed_num_pts};
171 
172  // Compute fluxes on faces
173  for (size_t i = 0; i < 3; ++i) {
174  auto& vars_upper_face = gsl::at(vars_on_upper_face, i);
175  auto& vars_lower_face = gsl::at(vars_on_lower_face, i);
177  make_not_null(&vars_upper_face));
179  make_not_null(&vars_lower_face));
180 
181  // Normal vectors in curved spacetime normalized by inverse
182  // spatial metric. Since we assume a Cartesian grid, this is
183  // relatively easy. Note that we use the sign convention on
184  // the normal vectors to be compatible with DG.
185  //
186  // Note that these normal vectors are on all faces inside the DG
187  // element since there are a bunch of subcells. We don't use the
188  // NormalCovectorAndMagnitude tag in the DataBox right now to avoid
189  // conflicts with the DG solver. We can explore in the future if it's
190  // possible to reuse that allocation.
191  const Scalar<DataVector> normalization{sqrt(
193  DataVector>>(vars_upper_face)
194  .get(i, i))};
195 
196  tnsr::i<DataVector, 3, Frame::Inertial> lower_outward_conormal{
197  reconstructed_num_pts, 0.0};
198  lower_outward_conormal.get(i) = 1.0 / get(normalization);
199 
200  tnsr::i<DataVector, 3, Frame::Inertial> upper_outward_conormal{
201  reconstructed_num_pts, 0.0};
202  upper_outward_conormal.get(i) = -lower_outward_conormal.get(i);
203  // Note: we probably should compute the normal vector in addition to
204  // the co-vector. Not a huge issue since we'll get an FPE right now if
205  // it's used by a Riemann solver.
206 
207  // Compute the packaged data
208  using dg_package_data_projected_tags = tmpl::append<
209  evolved_vars_tags, fluxes_tags, dg_package_data_temporary_tags,
210  typename DerivedCorrection::dg_package_data_primitive_tags>;
211  evolution::dg::Actions::detail::dg_package_data<System>(
212  make_not_null(&upper_packaged_data),
213  dynamic_cast<const DerivedCorrection&>(boundary_correction),
214  vars_upper_face, upper_outward_conormal, {std::nullopt}, *box,
215  typename DerivedCorrection::dg_package_data_volume_tags{},
216  dg_package_data_projected_tags{});
217 
218  evolution::dg::Actions::detail::dg_package_data<System>(
219  make_not_null(&lower_packaged_data),
220  dynamic_cast<const DerivedCorrection&>(boundary_correction),
221  vars_lower_face, lower_outward_conormal, {std::nullopt}, *box,
222  typename DerivedCorrection::dg_package_data_volume_tags{},
223  dg_package_data_projected_tags{});
224 
225  // Now need to check if any of our neighbors are doing DG,
226  // because if so then we need to use whatever boundary data
227  // they sent instead of what we computed locally.
228  //
229  // Note: We could check this beforehand to avoid the extra
230  // work of reconstruction and flux computations at the
231  // boundaries.
232  evolution::dg::subcell::correct_package_data<true>(
233  make_not_null(&lower_packaged_data),
234  make_not_null(&upper_packaged_data), i, element, subcell_mesh,
236 
237  // Compute the corrections on the faces. We only need to
238  // compute this once because we can just flip the normal
239  // vectors then
240  gsl::at(boundary_corrections, i).initialize(reconstructed_num_pts);
241  evolution::dg::subcell::compute_boundary_terms(
242  make_not_null(&gsl::at(boundary_corrections, i)),
243  dynamic_cast<const DerivedCorrection&>(boundary_correction),
244  upper_packaged_data, lower_packaged_data);
245  // We need to multiply by the normal vector normalization
246  gsl::at(boundary_corrections, i) *= get(normalization);
247  }
248  }
249  });
250 
251  // Now compute the actual time derivatives.
252  using variables_tag = typename System::variables_tag;
253  using dt_variables_tag = db::add_tag_prefix<::Tags::dt, variables_tag>;
255  tmpl::list<dt_variables_tag>,
256  tmpl::list<
271  tmpl::size_t<3>, Frame::Inertial>,
275  tmpl::size_t<3>, Frame::Inertial>,
279  [&cell_centered_logical_to_grid_inv_jacobian, &num_pts,
280  &boundary_corrections, &subcell_mesh, &one_over_delta_xi](
281  const auto dt_vars_ptr, const auto&... source_args) noexcept {
282  dt_vars_ptr->initialize(num_pts, 0.0);
288 
289  auto& dt_tilde_d = get<::Tags::dt<TildeD>>(*dt_vars_ptr);
290  auto& dt_tilde_tau = get<::Tags::dt<TildeTau>>(*dt_vars_ptr);
291  auto& dt_tilde_s = get<::Tags::dt<TildeS>>(*dt_vars_ptr);
292  auto& dt_tilde_b = get<::Tags::dt<TildeB>>(*dt_vars_ptr);
293  auto& dt_tilde_phi = get<::Tags::dt<TildePhi>>(*dt_vars_ptr);
294 
295  grmhd::ValenciaDivClean::ComputeSources::apply(
296  make_not_null(&dt_tilde_tau), make_not_null(&dt_tilde_s),
297  make_not_null(&dt_tilde_b), make_not_null(&dt_tilde_phi),
298  source_args...);
299 
300  for (size_t dim = 0; dim < 3; ++dim) {
301  Scalar<DataVector>& tilde_d_density_correction =
302  get<TildeD>(gsl::at(boundary_corrections, dim));
303  Scalar<DataVector>& tilde_tau_density_correction =
304  get<TildeTau>(gsl::at(boundary_corrections, dim));
305  tnsr::i<DataVector, 3, Frame::Inertial>&
306  tilde_s_density_correction =
307  get<TildeS>(gsl::at(boundary_corrections, dim));
308  tnsr::I<DataVector, 3, Frame::Inertial>&
309  tilde_b_density_correction =
310  get<TildeB>(gsl::at(boundary_corrections, dim));
311  Scalar<DataVector>& tilde_phi_density_correction =
312  get<TildePhi>(gsl::at(boundary_corrections, dim));
314  make_not_null(&get(dt_tilde_d)),
315  gsl::at(one_over_delta_xi, dim),
316  cell_centered_logical_to_grid_inv_jacobian.get(dim, dim),
317  get(tilde_d_density_correction), subcell_mesh.extents(), dim);
319  make_not_null(&get(dt_tilde_tau)),
320  gsl::at(one_over_delta_xi, dim),
321  cell_centered_logical_to_grid_inv_jacobian.get(dim, dim),
322  get(tilde_tau_density_correction), subcell_mesh.extents(), dim);
323  for (size_t d = 0; d < 3; ++d) {
325  make_not_null(&dt_tilde_s.get(d)),
326  gsl::at(one_over_delta_xi, dim),
327  cell_centered_logical_to_grid_inv_jacobian.get(dim, dim),
328  tilde_s_density_correction.get(d), subcell_mesh.extents(),
329  dim);
331  make_not_null(&dt_tilde_b.get(d)),
332  gsl::at(one_over_delta_xi, dim),
333  cell_centered_logical_to_grid_inv_jacobian.get(dim, dim),
334  tilde_b_density_correction.get(d), subcell_mesh.extents(),
335  dim);
336  }
338  make_not_null(&get(dt_tilde_phi)),
339  gsl::at(one_over_delta_xi, dim),
340  cell_centered_logical_to_grid_inv_jacobian.get(dim, dim),
341  get(tilde_phi_density_correction), subcell_mesh.extents(), dim);
342  }
343  },
344  box);
345  }
346 };
347 } // namespace grmhd::ValenciaDivClean::subcell
db::apply
constexpr auto apply(F &&f, const DataBox< BoxTags > &box, Args &&... args) noexcept
Apply the invokable f with argument Tags TagsList from DataBox box
Definition: DataBox.hpp:1148
grmhd::ValenciaDivClean::subcell::compute_fluxes
void compute_fluxes(const gsl::not_null< Variables< TagsList > * > vars) noexcept
Helper function that calls ComputeFluxes by retrieving the return and argument tags from vars.
Definition: ComputeFluxes.hpp:27
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
grmhd::ValenciaDivClean::subcell::TimeDerivative
Compute the time derivative on the subcell grid using FD reconstruction.
Definition: TimeDerivative.hpp:53
hydro::Tags::Pressure
The fluid pressure .
Definition: Tags.hpp:119
Frame::Inertial
Definition: IndexType.hpp:44
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
grmhd::ValenciaDivClean::subcell
Code required by the DG-subcell/FD hybrid solver.
Definition: ComputeFluxes.hpp:11
Tags.hpp
call_with_dynamic_type
Result call_with_dynamic_type(Base *const obj, Callable &&f) noexcept
Call a functor with the derived type of a base class pointer.
Definition: FakeVirtual.hpp:103
Mesh::basis
const std::array< Spectral::Basis, Dim > & basis() const noexcept
The basis chosen in each dimension of the grid.
Definition: Mesh.hpp:182
db::add_tag_prefix
typename detail::add_tag_prefix_impl< Prefix, Tag, Args... >::type add_tag_prefix
Definition: PrefixHelpers.hpp:51
Tags::Variables
Definition: VariablesTag.hpp:21
db::get
const auto & get(const DataBox< TagList > &box) noexcept
Retrieve the item with tag Tag from the DataBox.
Definition: DataBox.hpp:791
Tags::Flux
Prefix indicating a flux.
Definition: Prefixes.hpp:40
Element
Definition: Element.hpp:29
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.
DataBox.hpp
hydro::Tags::LorentzFactor
The Lorentz factor , where is the spatial velocity of the fluid.
Definition: Tags.hpp:64
cstddef
Assert.hpp
grmhd::ValenciaDivClean::Tags::TildeS
The densitized momentum density .
Definition: Tags.hpp:42
hydro::Tags::SpecificEnthalpy
The relativistic specific enthalpy .
Definition: Tags.hpp:167
array
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
grmhd::ValenciaDivClean::fd::Reconstructor
The base class from which all reconstruction schemes must inherit.
Definition: Reconstructor.hpp:20
evolution::dg::subcell::Tags::OnSubcellFaces
Mark a tag as the being on the subcell faces.
Definition: OnSubcellFaces.hpp:20
grmhd::ValenciaDivClean::Tags::ConstraintDampingParameter
The constraint damping parameter for divergence cleaning.
Definition: Tags.hpp:88
DataVector
Stores a collection of function values.
Definition: DataVector.hpp:46
evolution::dg::subcell::add_cartesian_flux_divergence
void add_cartesian_flux_divergence(gsl::not_null< DataVector * > dt_var, double one_over_delta, const DataVector &inv_jacobian, const DataVector &boundary_correction, const Index< 1 > &subcell_extents, size_t dimension) noexcept
Compute and add the 2nd-order flux divergence on a Cartesian mesh to the cell-centered time derivativ...
Element::id
const ElementId< VolumeDim > & id() const noexcept
A unique ID for the Element.
Definition: Element.hpp:63
gr::Tags::Shift
Definition: Tags.hpp:48
std::decay_t
grmhd::ValenciaDivClean::Tags::TildeB
The densitized magnetic field .
Definition: Tags.hpp:49
ASSERT
#define ASSERT(a, m)
Assert that an expression should be true.
Definition: Assert.hpp:49
Element::external_boundaries
const std::unordered_set< Direction< VolumeDim > > & external_boundaries() const noexcept
The directions of the faces of the Element that are external boundaries.
Definition: Element.hpp:51
Variables.hpp
Element.hpp
Mesh
Holds the number of grid points, basis, and quadrature in each direction of the computational grid.
Definition: Mesh.hpp:49
hydro::Tags::LorentzFactorTimesSpatialVelocity
The Lorentz factor times the spatial velocity .
Definition: Tags.hpp:185
Scalar
Tensor< T, Symmetry<>, index_list<> > Scalar
Definition: TypeAliases.hpp:21
Gsl.hpp
db::mutate_apply
constexpr decltype(auto) mutate_apply(F &&f, const gsl::not_null< DataBox< BoxTags > * > box, Args &&... args) noexcept
Apply the invokable f mutating items MutateTags and taking as additional arguments ArgumentTags and a...
Definition: DataBox.hpp:1257
evolution::dg::Tags::MortarData
Data on mortars, indexed by (Direction, ElementId) pairs.
Definition: MortarTags.hpp:27
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
Tags::deriv
Prefix indicating spatial derivatives.
Definition: PartialDerivatives.hpp:52
optional
grmhd::ValenciaDivClean::Tags::TildeD
The densitized rest-mass density .
Definition: Tags.hpp:31
Mesh::quadrature
const std::array< Spectral::Quadrature, Dim > & quadrature() const noexcept
The quadrature chosen in each dimension of the grid.
Definition: Mesh.hpp:196
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
grmhd::ValenciaDivClean::Tags::TildePhi
The densitized divergence-cleaning field .
Definition: Tags.hpp:55
gr::Tags::ExtrinsicCurvature
Definition: Tags.hpp:116
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
type_traits
gr::Tags::Lapse
Definition: Tags.hpp:52
TMPL.hpp
gr::Tags::InverseSpatialMetric
Inverse of the spatial metric.
Definition: Tags.hpp:33
Mesh.hpp
gsl::not_null
Require a pointer to not be a nullptr
Definition: ReadSpecPiecewisePolynomial.hpp:13
hydro::Tags::RestMassDensity
The rest-mass density .
Definition: Tags.hpp:125
grmhd::ValenciaDivClean::Tags::TildeTau
The densitized energy density .
Definition: Tags.hpp:36
gr::Tags::SqrtDetSpatialMetric
Definition: Tags.hpp:44