ComputeTimeDerivative.hpp
1 // Distributed under the MIT License.
2 // See LICENSE.txt for details.
3 
4 #pragma once
5 
6 #include <optional>
7 #include <tuple>
8 #include <type_traits>
9 #include <unordered_set>
10 #include <utility>
11 #include <vector>
12 
14 #include "DataStructures/DataBox/PrefixHelpers.hpp"
16 #include "DataStructures/DataVector.hpp"
19 #include "DataStructures/VariablesTag.hpp"
20 #include "Domain/CoordinateMaps/Tags.hpp"
21 #include "Domain/InterfaceHelpers.hpp"
23 #include "Domain/Structure/OrientationMapHelpers.hpp"
24 #include "Domain/Tags.hpp"
25 #include "Domain/TagsTimeDependent.hpp"
26 #include "Evolution/BoundaryCorrectionTags.hpp"
27 #include "Evolution/DiscontinuousGalerkin/Actions/BoundaryConditionsImpl.hpp"
28 #include "Evolution/DiscontinuousGalerkin/Actions/ComputeTimeDerivativeHelpers.hpp"
29 #include "Evolution/DiscontinuousGalerkin/Actions/InternalMortarDataImpl.hpp"
30 #include "Evolution/DiscontinuousGalerkin/Actions/NormalCovectorAndMagnitude.hpp"
31 #include "Evolution/DiscontinuousGalerkin/Actions/PackageDataImpl.hpp"
32 #include "Evolution/DiscontinuousGalerkin/Actions/VolumeTermsImpl.hpp"
33 #include "Evolution/DiscontinuousGalerkin/InboxTags.hpp"
34 #include "Evolution/DiscontinuousGalerkin/MortarData.hpp"
35 #include "Evolution/DiscontinuousGalerkin/MortarTags.hpp"
36 #include "Evolution/DiscontinuousGalerkin/NormalVectorTags.hpp"
37 #include "Evolution/DiscontinuousGalerkin/UsingSubcell.hpp"
38 #include "NumericalAlgorithms/DiscontinuousGalerkin/Formulation.hpp"
39 #include "NumericalAlgorithms/DiscontinuousGalerkin/MortarHelpers.hpp"
40 #include "NumericalAlgorithms/DiscontinuousGalerkin/Tags.hpp"
41 #include "NumericalAlgorithms/DiscontinuousGalerkin/Tags/Formulation.hpp"
44 #include "NumericalAlgorithms/Spectral/Projection.hpp"
45 #include "Parallel/GlobalCache.hpp"
46 #include "Time/Actions/SelfStartActions.hpp"
47 #include "Time/Tags.hpp"
48 #include "Time/TakeStep.hpp"
49 #include "Utilities/Algorithm.hpp"
50 #include "Utilities/Gsl.hpp"
51 #include "Utilities/TMPL.hpp"
52 
53 /// \cond
54 namespace evolution::dg::subcell {
55 // We use a forward declaration instead of including a header file to avoid
56 // coupling to the DG-subcell libraries for executables that don't use subcell.
57 template <typename Metavariables, typename DbTagsList>
59 prepare_neighbor_data(gsl::not_null<db::DataBox<DbTagsList>*> box) noexcept;
60 } // namespace evolution::dg::subcell
61 namespace tuples {
62 template <typename...>
63 class TaggedTuple;
64 } // namespace tuples
65 /// \endcond
66 
67 namespace evolution::dg::Actions {
68 /*!
69  * \brief Computes the time derivative for a DG time step.
70  *
71  * Computes the volume fluxes, the divergence of the fluxes and all additional
72  * interior contributions to the time derivatives (both nonconservative products
73  * and source terms). The internal mortar data is also computed.
74  *
75  * The general first-order hyperbolic evolution equation solved for conservative
76  * systems is:
77  *
78  * \f{align*}{
79  * \frac{\partial u_\alpha}{\partial \hat{t}}
80  * + \partial_{i}
81  * \left(F^i_\alpha - v^i_g u_\alpha\right)
82  * = S_\alpha-u_\alpha\partial_i v^i_g,
83  * \f}
84  *
85  * where \f$F^i_{\alpha}\f$ are the fluxes when the mesh isn't moving,
86  * \f$v^i_g\f$ is the velocity of the mesh, \f$u_{\alpha}\f$ are the evolved
87  * variables, \f$S_{\alpha}\f$ are the source terms, and \f$\hat{t}\f$ is the
88  * time in the logical frame. For evolution equations that do not have any
89  * fluxes and only nonconservative products we evolve:
90  *
91  * \f{align*}{
92  * \frac{\partial u_\alpha}{\partial \hat{t}}
93  * +\left(B^i_{\alpha\beta}-v^i_g \delta_{\alpha\beta}
94  * \right)\partial_{i}u_\beta = S_\alpha.
95  * \f}
96  *
97  * Finally, for equations with both conservative terms and nonconservative
98  * products we use:
99  *
100  * \f{align*}{
101  * \frac{\partial u_\alpha}{\partial \hat{t}}
102  * + \partial_{i}
103  * \left(F^i_\alpha - v^i_g u_\alpha\right)
104  * +B^i_{\alpha\beta}\partial_{i}u_\beta
105  * = S_\alpha-u_\alpha\partial_i v^i_g,
106  * \f}
107  *
108  * where \f$B^i_{\alpha\beta}\f$ is the matrix for the nonconservative products.
109  *
110  * ### Volume Terms
111  *
112  * The mesh velocity is added to the flux automatically if the mesh is moving.
113  * That is,
114  *
115  * \f{align*}{
116  * F^i_{\alpha}\to F^i_{\alpha}-v^i_{g} u_{\alpha}
117  * \f}
118  *
119  * The source terms are also altered automatically by adding:
120  *
121  * \f{align*}{
122  * -u_\alpha \partial_i v^i_g,
123  * \f}
124  *
125  * For systems with equations that only contain nonconservative products, the
126  * following mesh velocity is automatically added to the time derivative:
127  *
128  * \f{align*}{
129  * v^i_g \partial_i u_\alpha,
130  * \f}
131  *
132  * \note The term is always added in the `Frame::Inertial` frame, and the plus
133  * sign arises because we add it to the time derivative.
134  *
135  * Here are examples of the `TimeDerivative` struct used to compute the volume
136  * time derivative. This struct is what the type alias
137  * `System::compute_volume_time_derivative` points to. The time derivatives are
138  * as `gsl::not_null` first, then the temporary tags as `gsl::not_null`,
139  * followed by the `argument_tags`. These type aliases are given by
140  *
141  * \snippet ComputeTimeDerivativeImpl.tpp dt_ta
142  *
143  * for the examples. For a conservative system without primitives the `apply`
144  * function would look like
145  *
146  * \snippet ComputeTimeDerivativeImpl.tpp dt_con
147  *
148  * For a nonconservative system it would be
149  *
150  * \snippet ComputeTimeDerivativeImpl.tpp dt_nc
151  *
152  * And finally, for a mixed conservative-nonconservative system with primitive
153  * variables
154  *
155  * \snippet ComputeTimeDerivativeImpl.tpp dt_mp
156  *
157  * Uses:
158  * - System:
159  * - `variables_tag`
160  * - `flux_variables`
161  * - `gradient_variables`
162  * - `compute_volume_time_derivative_terms`
163  *
164  * - DataBox:
165  * - Items in `system::compute_volume_time_derivative_terms::argument_tags`
166  * - `domain::Tags::MeshVelocity<Metavariables::volume_dim>`
167  * - `Metavariables::system::variables_tag`
168  * - `Metavariables::system::flux_variables`
169  * - `Metavariables::system::gradient_variables`
170  * - `domain::Tags::DivMeshVelocity`
171  * - `DirectionsTag`,
172  * - Required interface items for `Metavariables::system::normal_dot_fluxes`
173  *
174  * DataBox changes:
175  * - Adds: nothing
176  * - Removes: nothing
177  * - Modifies:
178  * - db::add_tag_prefix<Tags::Flux, variables_tag,
179  * tmpl::size_t<system::volume_dim>, Frame::Inertial>
180  * - `Tags::dt<system::variable_tags>`
181  * - Tags::Interface<
182  * DirectionsTag, db::add_tag_prefix<Tags::NormalDotFlux, variables_tag>>
183  * - `Tags::Mortars<typename BoundaryScheme::mortar_data_tag, VolumeDim>`
184  *
185  * ### Internal Boundary Terms
186  *
187  * Internal boundary terms are computed from the
188  * `System::boundary_correction_base` type alias. This type alias must point to
189  * a base class with `creatable_classes`. Each concrete boundary correction must
190  * specify:
191  *
192  * - type alias template `dg_package_field_tags`. These are what will be
193  * returned by `gsl::not_null` from the `dg_package_data` member function.
194  *
195  * - type alias template `dg_package_temporary_tags`. These are temporary tags
196  * that are projected to the face and then passed to the `dg_package_data`
197  * function.
198  *
199  * - type alias template `dg_package_primitive_tags`. These are the primitive
200  * variables (if any) that are projected to the face and then passed to
201  * `dg_package_data`.
202  *
203  * - type alias template `dg_package_volume_tags`. These are tags that are not
204  * projected to the interface and are retrieved directly from the `DataBox`.
205  * The equation of state for hydrodynamics systems is an example of what would
206  * be a "volume tag".
207  *
208  * A `static constexpr bool need_normal_vector` must be specified. If `true`
209  * then the normal vector is computed from the normal covector. This is
210  * currently not implemented.
211  *
212  * The `dg_package_data` function takes as arguments `gsl::not_null` of the
213  * `dg_package_data_field_tags`, then the projected evolved variables, the
214  * projected fluxes, the projected temporaries, the projected primitives, the
215  * unit normal covector, mesh velocity, normal dotted into the mesh velocity,
216  * the `volume_tags`, and finally the `dg::Formulation`. The `dg_package_data`
217  * function must compute all ingredients for the boundary correction, including
218  * mesh-velocity-corrected characteristic speeds. However, the projected fluxes
219  * passed in are \f$F^i - u v^i_g\f$ (the mesh velocity term is already
220  * included). The `dg_package_data` function must also return a `double` that is
221  * the maximum absolute characteristic speed over the entire face. This will be
222  * used for checking that the time step doesn't violate the CFL condition.
223  *
224  * Here is an example of the type aliases and `bool`:
225  *
226  * \snippet ComputeTimeDerivativeImpl.tpp bt_ta
227  *
228  * The normal vector requirement is:
229  *
230  * \snippet ComputeTimeDerivativeImpl.tpp bt_nnv
231  *
232  * For a conservative system with primitive variables and using the `TimeStepId`
233  * as a volume tag the `dg_package_data` function looks like:
234  *
235  * \snippet ComputeTimeDerivativeImpl.tpp bt_cp
236  *
237  * For a mixed conservative-nonconservative system with primitive variables and
238  * using the `TimeStepId` as a volume tag the `dg_package_data` function looks
239  * like:
240  *
241  * \snippet ComputeTimeDerivativeImpl.tpp bt_mp
242  *
243  * Uses:
244  * - System:
245  * - `boundary_correction`
246  * - `variables_tag`
247  * - `flux_variables`
248  * - `gradients_tags`
249  * - `compute_volume_time_derivative`
250  * - `has_primitive_and_conservative_vars`
251  * - `primitive_variables_tag` if system has primitive variables
252  *
253  * - DataBox:
254  * - `domain::Tags::Element<Dim>`
255  * - `domain::Tags::Mesh<Dim>`
256  * - `evolution::dg::Tags::MortarMesh<Dim>`
257  * - `evolution::dg::Tags::MortarSize<Dim>`
258  * - `evolution::dg::Tags::MortarData<Dim>`
259  * - `Tags::TimeStepId`
260  * - \code{.cpp}
261  * domain::Tags::Interface<domain::Tags::InternalDirections<Dim>,
262  * domain::Tags::Mesh<Dim - 1>>
263  * \endcode
264  * - \code{.cpp}
265  * domain::Tags::Interface<
266  * domain::Tags::InternalDirections<Dim>,
267  * ::Tags::Normalized<
268  * domain::Tags::UnnormalizedFaceNormal<Dim, Frame::Inertial>>>
269  * \endcode
270  * - \code{.cpp}
271  * domain::Tags::Interface<
272  * domain::Tags::InternalDirections<Dim>,
273  * domain::Tags::MeshVelocity<Dim, Frame::Inertial>>
274  * \endcode
275  * - `Metavariables::system::variables_tag`
276  * - `Metavariables::system::flux_variables`
277  * - `Metavariables::system::primitive_tags` if exists
278  * - `system::boundary_correction_base::dg_package_data_volume_tags`
279  *
280  * DataBox changes:
281  * - Adds: nothing
282  * - Removes: nothing
283  * - Modifies:
284  * - `evolution::dg::Tags::MortarData<Dim>`
285  */
286 template <typename Metavariables>
288  using inbox_tags =
290  Metavariables::volume_dim>>;
291  using const_global_cache_tags = tmpl::append<tmpl::list<
294 
295  template <typename DbTagsList, typename... InboxTags, typename ArrayIndex,
296  typename ActionList, typename ParallelComponent>
297  static std::tuple<db::DataBox<DbTagsList>&&> apply(
298  db::DataBox<DbTagsList>& box,
301  const ArrayIndex& /*array_index*/, ActionList /*meta*/,
302  const ParallelComponent* /*meta*/) noexcept; // NOLINT const
303 
304  private:
305  // The below functions will be removed once we've finished writing the new
306  // method of handling boundaries.
307  template <typename... NormalDotFluxTags, typename... Args>
308  static void apply_flux(
309  gsl::not_null<Variables<tmpl::list<NormalDotFluxTags...>>*> boundary_flux,
310  const Args&... boundary_variables) noexcept {
311  Metavariables::system::normal_dot_fluxes::apply(
312  make_not_null(&get<NormalDotFluxTags>(*boundary_flux))...,
313  boundary_variables...);
314  }
315 
316  template <typename DbTagsList>
317  static void boundary_terms_nonconservative_products(
318  gsl::not_null<db::DataBox<DbTagsList>*> box) noexcept;
319 
320  template <size_t VolumeDim, typename BoundaryScheme, typename DbTagsList>
321  static void fill_mortar_data_for_internal_boundaries(
322  gsl::not_null<db::DataBox<DbTagsList>*> box) noexcept;
323 
324  template <typename ParallelComponent, typename DbTagsList>
325  static void send_data_for_fluxes(
327  gsl::not_null<db::DataBox<DbTagsList>*> box) noexcept;
328 };
329 
330 template <typename Metavariables>
331 template <typename DbTagsList, typename... InboxTags, typename ArrayIndex,
332  typename ActionList, typename ParallelComponent>
335  db::DataBox<DbTagsList>& box,
338  const ArrayIndex& /*array_index*/, ActionList /*meta*/,
339  const ParallelComponent* const /*meta*/) noexcept { // NOLINT const
340  static constexpr size_t volume_dim = Metavariables::volume_dim;
341  using system = typename Metavariables::system;
342  using variables_tag = typename system::variables_tag;
343  using dt_variables_tag = db::add_tag_prefix<::Tags::dt, variables_tag>;
344  using partial_derivative_tags = typename system::gradient_variables;
345  using flux_variables = typename system::flux_variables;
346  using compute_volume_time_derivative_terms =
347  typename system::compute_volume_time_derivative_terms;
348 
349  const Mesh<volume_dim>& mesh = db::get<::domain::Tags::Mesh<volume_dim>>(box);
350  const ::dg::Formulation dg_formulation =
351  db::get<::dg::Tags::Formulation>(box);
352  ASSERT(alg::all_of(mesh.basis(),
353  [&mesh](const Spectral::Basis current_basis) noexcept {
354  return current_basis == mesh.basis(0);
355  }),
356  "An isotropic basis must be used in the evolution code. While "
357  "theoretically this restriction could be lifted, the simplification "
358  "it offers are quite substantial. Relaxing this assumption is likely "
359  "to require quite a bit of careful code refactoring and debugging.");
361  mesh.quadrature(),
362  [&mesh](const Spectral::Quadrature current_quadrature) noexcept {
363  return current_quadrature == mesh.quadrature(0);
364  }),
365  "An isotropic quadrature must be used in the evolution code. While "
366  "theoretically this restriction could be lifted, the simplification "
367  "it offers are quite substantial. Relaxing this assumption is likely "
368  "to require quite a bit of careful code refactoring and debugging.");
369 
370  // Allocate the Variables classes needed for the time derivative
371  // computation.
372  //
373  // This is factored out so that we will be able to do ADER-DG/CG where a
374  // spacetime polynomial is constructed by solving implicit equations in time
375  // using a Picard iteration. A high-order initial guess is needed to
376  // efficiently construct the ADER spacetime solution. This initial guess is
377  // obtained using continuous RK methods, and so we will want to reuse
378  // buffers. Thus, the volume_terms function returns by reference rather than
379  // by value.
380  Variables<typename compute_volume_time_derivative_terms::temporary_tags>
381  temporaries{mesh.number_of_grid_points()};
382  Variables<db::wrap_tags_in<::Tags::Flux, flux_variables,
383  tmpl::size_t<volume_dim>, Frame::Inertial>>
384  volume_fluxes{mesh.number_of_grid_points()};
385  Variables<db::wrap_tags_in<::Tags::deriv, partial_derivative_tags,
386  tmpl::size_t<volume_dim>, Frame::Inertial>>
387  partial_derivs{mesh.number_of_grid_points()};
388 
389  const Scalar<DataVector>* det_inverse_jacobian = nullptr;
390  if constexpr (tmpl::size<flux_variables>::value != 0) {
391  if (dg_formulation == ::dg::Formulation::WeakInertial) {
392  det_inverse_jacobian = &db::get<
394  }
395  }
397  tmpl::list<dt_variables_tag>,
398  typename compute_volume_time_derivative_terms::argument_tags>(
399  [&dg_formulation, &det_inverse_jacobian,
400  &div_mesh_velocity = db::get<::domain::Tags::DivMeshVelocity>(box),
401  &evolved_variables = db::get<variables_tag>(box),
402  &inertial_coordinates =
403  db::get<domain::Tags::Coordinates<volume_dim, Frame::Inertial>>(box),
404  &logical_to_inertial_inv_jacobian =
406  Frame::Inertial>>(box),
407  &mesh,
408  &mesh_velocity = db::get<::domain::Tags::MeshVelocity<volume_dim>>(box),
409  &partial_derivs, &temporaries, &volume_fluxes](
410  const gsl::not_null<Variables<
412  dt_vars_ptr,
413  const auto&... time_derivative_args) noexcept {
414  detail::volume_terms<compute_volume_time_derivative_terms>(
415  dt_vars_ptr, make_not_null(&volume_fluxes),
416  make_not_null(&partial_derivs), make_not_null(&temporaries),
417  evolved_variables, dg_formulation, mesh, inertial_coordinates,
418  logical_to_inertial_inv_jacobian, det_inverse_jacobian,
419  mesh_velocity, div_mesh_velocity, time_derivative_args...);
420  },
421  make_not_null(&box));
422 
423  const auto& boundary_correction =
424  db::get<evolution::Tags::BoundaryCorrection<system>>(box);
425  using derived_boundary_corrections =
426  typename std::decay_t<decltype(boundary_correction)>::creatable_classes;
427 
428  const Variables<detail::get_primitive_vars_tags_from_system<system>>*
429  primitive_vars{nullptr};
430  if constexpr (system::has_primitive_and_conservative_vars) {
431  primitive_vars = &db::get<typename system::primitive_variables_tag>(box);
432  }
433 
434  static_assert(
435  tmpl::all<derived_boundary_corrections, std::is_final<tmpl::_1>>::value,
436  "All createable classes for boundary corrections must be marked "
437  "final.");
438  tmpl::for_each<derived_boundary_corrections>(
439  [&boundary_correction, &box, &partial_derivs, &primitive_vars,
440  &temporaries, &volume_fluxes](auto derived_correction_v) noexcept {
441  using DerivedCorrection =
442  tmpl::type_from<decltype(derived_correction_v)>;
443  if (typeid(boundary_correction) == typeid(DerivedCorrection)) {
444  // Compute internal boundary quantities on the mortar for sides
445  // of the element that have neighbors, i.e. they are not an
446  // external side.
447  // Note: this call mutates:
448  // - evolution::dg::Tags::NormalCovectorAndMagnitude<Dim>,
449  // - evolution::dg::Tags::MortarData<Dim>
450  detail::internal_mortar_data<system, volume_dim>(
451  make_not_null(&box),
452  dynamic_cast<const DerivedCorrection&>(boundary_correction),
453  db::get<variables_tag>(box), volume_fluxes, temporaries,
454  primitive_vars,
455  typename DerivedCorrection::dg_package_data_volume_tags{});
456 
457  detail::apply_boundary_conditions_on_all_external_faces<system,
458  volume_dim>(
459  make_not_null(&box),
460  dynamic_cast<const DerivedCorrection&>(boundary_correction),
461  temporaries, volume_fluxes, partial_derivs, primitive_vars);
462  }
463  });
464 
465  if constexpr (Metavariables::local_time_stepping) {
466  take_step(make_not_null(&box), cache);
467  }
468 
469  send_data_for_fluxes<ParallelComponent>(make_not_null(&cache),
470  make_not_null(&box));
471  return {std::move(box)};
472 }
473 
474 template <typename Metavariables>
475 template <typename ParallelComponent, typename DbTagsList>
476 void ComputeTimeDerivative<Metavariables>::send_data_for_fluxes(
478  const gsl::not_null<db::DataBox<DbTagsList>*> box) noexcept {
479  constexpr size_t volume_dim = Metavariables::volume_dim;
480  auto& receiver_proxy =
481  Parallel::get_parallel_component<ParallelComponent>(*cache);
482  const auto& element = db::get<domain::Tags::Element<volume_dim>>(*box);
483 
484  const auto& time_step_id = db::get<::Tags::TimeStepId>(*box);
485  const auto& all_mortar_data =
486  db::get<evolution::dg::Tags::MortarData<volume_dim>>(*box);
487  const auto& mortar_meshes =
488  get<evolution::dg::Tags::MortarMesh<volume_dim>>(*box);
489 
491  all_neighbor_data_for_reconstruction = std::nullopt;
492  if constexpr (using_subcell_v<Metavariables>) {
493  all_neighbor_data_for_reconstruction =
494  evolution::dg::subcell::prepare_neighbor_data<Metavariables>(box);
495  }
496 
497  for (const auto& [direction, neighbors] : element.neighbors()) {
498  const auto& orientation = neighbors.orientation();
499  const auto direction_from_neighbor = orientation(direction.opposite());
500 
501  std::vector<double> ghost_and_subcell_data{};
502  if constexpr (using_subcell_v<Metavariables>) {
503  ASSERT(all_neighbor_data_for_reconstruction.has_value(),
504  "Trying to do DG-subcell but the ghost and subcell data for the "
505  "neighbor has not been set.");
506  ghost_and_subcell_data =
507  std::move(all_neighbor_data_for_reconstruction.value()[direction]);
508  }
509 
510  for (const auto& neighbor : neighbors) {
511  const std::pair mortar_id{direction, neighbor};
512 
513  std::pair<Mesh<volume_dim - 1>, std::vector<double>>
514  neighbor_boundary_data_on_mortar{};
515  ASSERT(time_step_id == all_mortar_data.at(mortar_id).time_step_id(),
516  "The current time step id of the volume is "
517  << time_step_id
518  << "but the time step id on the mortar with mortar id "
519  << mortar_id << " is "
520  << all_mortar_data.at(mortar_id).time_step_id());
521 
522  // Reorient the data to the neighbor orientation if necessary
523  if (LIKELY(orientation.is_aligned())) {
524  neighbor_boundary_data_on_mortar =
525  *all_mortar_data.at(mortar_id).local_mortar_data();
526  } else {
527  const auto& slice_extents = mortar_meshes.at(mortar_id).extents();
528  neighbor_boundary_data_on_mortar.first =
529  all_mortar_data.at(mortar_id).local_mortar_data()->first;
530  neighbor_boundary_data_on_mortar.second = orient_variables_on_slice(
531  all_mortar_data.at(mortar_id).local_mortar_data()->second,
532  slice_extents, direction.dimension(), orientation);
533  }
534 
535  const TimeStepId& next_time_step_id = [&box]() noexcept {
536  if (Metavariables::local_time_stepping) {
537  return db::get<::Tags::Next<::Tags::TimeStepId>>(*box);
538  } else {
539  return db::get<::Tags::TimeStepId>(*box);
540  }
541  }();
542 
545  data{neighbor_boundary_data_on_mortar.first,
546  ghost_and_subcell_data,
547  {std::move(neighbor_boundary_data_on_mortar.second)},
548  next_time_step_id};
549 
550  // Send mortar data (the `std::tuple` named `data`) to neighbor
553  volume_dim>>(
554  receiver_proxy[neighbor], time_step_id,
555  std::make_pair(std::pair{direction_from_neighbor, element.id()},
556  data));
557  }
558  }
559 
560  if constexpr (Metavariables::local_time_stepping) {
561  using variables_tag = typename Metavariables::system::variables_tag;
562  using dt_variables_tag = db::add_tag_prefix<::Tags::dt, variables_tag>;
563  // We assume isotropic quadrature, i.e. the quadrature is the same in
564  // all directions.
565  const bool using_gauss_points =
566  db::get<domain::Tags::Mesh<volume_dim>>(*box).quadrature() ==
567  make_array<volume_dim>(Spectral::Quadrature::Gauss);
568 
569  const Scalar<DataVector> volume_det_inv_jacobian{};
570  if (using_gauss_points) {
571  // NOLINTNEXTLINE
572  const_cast<DataVector&>(get(volume_det_inv_jacobian))
573  .set_data_ref(make_not_null(&const_cast<DataVector&>( // NOLINT
575  Frame::Logical, Frame::Inertial>>(*box)))));
576  }
577 
578  // Add face normal and Jacobian determinants to the local mortar data. We
579  // only need the Jacobians if we are using Gauss points. Then copy over
580  // into the boundary history, since that's what the LTS steppers use.
581  //
582  // The boundary history coupling computation (which computes the _lifted_
583  // boundary correction) returns a Variables<dt<EvolvedVars>> instead of
584  // using the `NormalDotNumericalFlux` prefix tag. This is because the
585  // returned quantity is more a `dt` quantity than a
586  // `NormalDotNormalDotFlux` since it's been lifted to the volume.
588  const auto integration_order =
589  db::get<::Tags::HistoryEvolvedVariables<>>(*box).integration_order();
590  db::mutate<evolution::dg::Tags::MortarData<volume_dim>,
592  volume_dim, typename dt_variables_tag::type>>(
593  box,
594  [&element, integration_order, &time_step_id, using_gauss_points,
595  &volume_det_inv_jacobian](
596  const gsl::not_null<
598  boost::hash<Key>>*>
599  mortar_data,
600  const gsl::not_null<
601  std::unordered_map<Key,
605  typename dt_variables_tag::type>,
606  boost::hash<Key>>*>
607  boundary_data_history,
608  const Mesh<volume_dim>& volume_mesh,
609  const DirectionMap<
610  volume_dim,
611  std::optional<Variables<tmpl::list<
614  normal_covector_and_magnitude) noexcept {
615  Scalar<DataVector> volume_det_jacobian{};
616  Scalar<DataVector> face_det_jacobian{};
617  if (using_gauss_points) {
618  get(volume_det_jacobian) = 1.0 / get(volume_det_inv_jacobian);
619  }
620  for (const auto& [direction, neighbors_in_direction] :
621  element.neighbors()) {
622  // We can perform projections once for all neighbors in the
623  // direction because we care about the _face_ mesh, not the mortar
624  // mesh.
625  ASSERT(
626  normal_covector_and_magnitude.at(direction).has_value(),
627  "The normal covector and magnitude have not been computed.");
628  const Scalar<DataVector>& face_normal_magnitude =
629  get<evolution::dg::Tags::MagnitudeOfNormal>(
630  *normal_covector_and_magnitude.at(direction));
631  if (using_gauss_points) {
632  const Matrix identity{};
633  auto interpolation_matrices =
634  make_array<Metavariables::volume_dim>(std::cref(identity));
635  const std::pair<Matrix, Matrix>& matrices =
637  volume_mesh.slice_through(direction.dimension()));
638  gsl::at(interpolation_matrices, direction.dimension()) =
639  direction.side() == Side::Upper ? matrices.second
640  : matrices.first;
641  if (get(face_det_jacobian).size() !=
642  get(face_normal_magnitude).size()) {
643  get(face_det_jacobian) =
644  DataVector{get(face_normal_magnitude).size()};
645  }
646  apply_matrices(make_not_null(&get(face_det_jacobian)),
647  interpolation_matrices, get(volume_det_jacobian),
648  volume_mesh.extents());
649  }
650 
651  for (const auto& neighbor : neighbors_in_direction) {
652  const std::pair mortar_id{direction, neighbor};
653  if (using_gauss_points) {
654  mortar_data->at(mortar_id).insert_local_geometric_quantities(
655  volume_det_inv_jacobian, face_det_jacobian,
656  face_normal_magnitude);
657  } else {
658  mortar_data->at(mortar_id).insert_local_face_normal_magnitude(
659  face_normal_magnitude);
660  }
661  ASSERT(boundary_data_history->count(mortar_id) != 0,
662  "Could not insert the mortar data for "
663  << mortar_id
664  << " because the unordered map has not been "
665  "initialized "
666  "to have the mortar id.");
667  boundary_data_history->at(mortar_id).local_insert(
668  time_step_id, std::move(mortar_data->at(mortar_id)));
669  boundary_data_history->at(mortar_id).integration_order(
670  integration_order);
671  mortar_data->at(mortar_id) =
672  MortarData<Metavariables::volume_dim>{};
673  }
674  }
675  },
676  db::get<domain::Tags::Mesh<volume_dim>>(*box),
677  db::get<evolution::dg::Tags::NormalCovectorAndMagnitude<volume_dim>>(
678  *box));
679  }
680 }
681 } // namespace evolution::dg::Actions
Matrix
A dynamically sized matrix of doubles with column-major storage.
Definition: Matrix.hpp:19
apply_matrices
void apply_matrices(const gsl::not_null< Variables< VariableTags > * > result, const std::array< MatrixType, Dim > &matrices, const Variables< VariableTags > &u, const Index< Dim > &extents) noexcept
Multiply by matrices in each dimension.
Definition: ApplyMatrices.hpp:67
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
utility
Frame::Inertial
Definition: IndexType.hpp:44
Spectral::Basis
Basis
The choice of basis functions for computing collocation points and weights.
Definition: Spectral.hpp:70
Parallel::GlobalCache
Definition: ElementReceiveInterpPoints.hpp:15
unordered_set
std::pair
GlobalCache.hpp
Tags.hpp
vector
std::system
T system(T... args)
std::size
T size(T... args)
db::add_tag_prefix
typename detail::add_tag_prefix_impl< Prefix, Tag, Args... >::type add_tag_prefix
Definition: PrefixHelpers.hpp:51
tuple
db::get
const auto & get(const DataBox< TagList > &box) noexcept
Retrieve the item with tag Tag from the DataBox.
Definition: DataBox.hpp:791
orient_variables_on_slice
Variables< TagsList > orient_variables_on_slice(const Variables< TagsList > &variables_on_slice, const Index< VolumeDim - 1 > &slice_extents, const size_t sliced_dim, const OrientationMap< VolumeDim > &orientation_of_neighbor) noexcept
Definition: OrientationMapHelpers.hpp:83
evolution::dg::subcell::fd::mesh
Mesh< Dim > mesh(const Mesh< Dim > &dg_mesh) noexcept
Computes the cell-centered finite-difference mesh from the DG mesh, using grid points per dimension,...
Tags::Flux
Prefix indicating a flux.
Definition: Prefixes.hpp:40
evolution::dg::subcell::prepare_neighbor_data
DirectionMap< Metavariables::volume_dim, std::vector< double > > prepare_neighbor_data(const gsl::not_null< db::DataBox< DbTagsList > * > box) noexcept
Add local data for our and our neighbor's relaxed discrete maximum principle troubled-cell indicator,...
Definition: PrepareNeighborData.hpp:52
evolution::Tags::BoundaryCorrection
The boundary correction used for coupling together neighboring cells or applying boundary conditions.
Definition: BoundaryCorrectionTags.hpp:37
ElementId
An ElementId uniquely labels an Element.
Definition: ElementId.hpp:51
evolution::dg::MortarData
Data on the mortar used to compute the boundary correction for the DG scheme.
Definition: MortarData.hpp:50
evolution::dg::Actions::ComputeTimeDerivative
Computes the time derivative for a DG time step.
Definition: ComputeTimeDerivative.hpp:287
DataBox.hpp
domain::Tags::InverseJacobian
The inverse Jacobian from the source frame to the target frame.
Definition: Tags.hpp:165
identity
tnsr::Ij< DataType, Dim, Frame::NoFrame > identity(const DataType &used_for_type) noexcept
returns the Identity matrix
Definition: Identity.hpp:14
evolution::dg::Actions
Actions for using the discontinuous Galerkin to evolve hyperbolic partial differential equations.
Definition: ApplyBoundaryCorrections.hpp:67
Spectral::boundary_interpolation_matrices
const std::pair< Matrix, Matrix > & boundary_interpolation_matrices(const Mesh< 1 > &mesh) noexcept
Matrices that interpolate to the lower and upper boundaries of the element.
tuples::TaggedTuple
An associative container that is indexed by structs.
Definition: TaggedTuple.hpp:271
evolution::dg::Tags::MagnitudeOfNormal
The magnitude of the unnormalized normal covector to the interface.
Definition: NormalVectorTags.hpp:18
DataVector
Stores a collection of function values.
Definition: DataVector.hpp:46
Spectral::Quadrature
Quadrature
The choice of quadrature method to compute integration weights.
Definition: Spectral.hpp:90
TimeStepId
Definition: TimeStepId.hpp:25
std::decay_t
ASSERT
#define ASSERT(a, m)
Assert that an expression should be true.
Definition: Assert.hpp:49
Variables.hpp
Parallel::receive_data
void receive_data(Proxy &&proxy, typename ReceiveTag::temporal_id temporal_id, ReceiveDataType &&receive_data, const bool enable_if_disabled=false) noexcept
Send the data args... to the algorithm running on proxy, and tag the message with the identifier temp...
Definition: Invoke.hpp:32
Mesh
Holds the number of grid points, basis, and quadrature in each direction of the computational grid.
Definition: Mesh.hpp:49
DirectionMap
Definition: DirectionMap.hpp:15
ActionTesting::cache
Parallel::GlobalCache< Metavariables > & cache(MockRuntimeSystem< Metavariables > &runner, const ArrayIndex &array_index) noexcept
Returns the GlobalCache of Component with index array_index.
Definition: MockRuntimeSystemFreeFunctions.hpp:382
evolution::dg::subcell
Implementation of a generic finite volume/conservative finite difference subcell limiter.
Definition: Actions.hpp:6
Scalar
Tensor< T, Symmetry<>, index_list<> > Scalar
Definition: TypeAliases.hpp:21
std::is_final
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
Tensor.hpp
TimeSteppers::BoundaryHistory
Definition: BoundaryHistory.hpp:31
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
LIKELY
#define LIKELY(x)
Definition: Gsl.hpp:67
evolution::dg::Tags::BoundaryCorrectionAndGhostCellsInbox
The inbox tag for boundary correction communication and DG-subcell ghost zone cells.
Definition: InboxTags.hpp:94
Tags::deriv
Prefix indicating spatial derivatives.
Definition: PartialDerivatives.hpp:52
Frame::Logical
Definition: IndexType.hpp:42
optional
evolution::dg::Tags::NormalCovector
The normal covector to the interface.
Definition: NormalVectorTags.hpp:24
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::Tags::MortarDataHistory
History of the data on mortars, indexed by (Direction, ElementId) pairs, and used by the linear multi...
Definition: MortarTags.hpp:43
Direction.hpp
PartialDerivatives.hpp
domain::Tags::DetInvJacobian
The determinant of the inverse Jacobian from the source frame to the target frame.
Definition: Tags.hpp:235
dg::Formulation
Formulation
The DG formulation to use.
Definition: Formulation.hpp:34
Prefixes.hpp
alg::all_of
decltype(auto) all_of(const Container &c, UnaryPredicate &&unary_predicate)
Convenience wrapper around std::all_of.
Definition: Algorithm.hpp:181
std::unordered_map
std::data
T data(T... args)
type_traits
TMPL.hpp
Mesh.hpp
gsl::not_null
Require a pointer to not be a nullptr
Definition: ReadSpecPiecewisePolynomial.hpp:13