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 "NumericalAlgorithms/DiscontinuousGalerkin/Formulation.hpp"
38 #include "NumericalAlgorithms/DiscontinuousGalerkin/MortarHelpers.hpp"
39 #include "NumericalAlgorithms/DiscontinuousGalerkin/Tags.hpp"
40 #include "NumericalAlgorithms/DiscontinuousGalerkin/Tags/Formulation.hpp"
43 #include "NumericalAlgorithms/Spectral/Projection.hpp"
44 #include "Parallel/GlobalCache.hpp"
45 #include "ParallelAlgorithms/DiscontinuousGalerkin/FluxCommunication.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 tuples {
55 template <typename...>
56 class TaggedTuple;
57 } // namespace tuples
58 /// \endcond
59 
60 namespace evolution::dg::Actions {
61 /*!
62  * \brief Computes the time derivative for a DG time step.
63  *
64  * Computes the volume fluxes, the divergence of the fluxes and all additional
65  * interior contributions to the time derivatives (both nonconservative products
66  * and source terms). The internal mortar data is also computed.
67  *
68  * The general first-order hyperbolic evolution equation solved for conservative
69  * systems is:
70  *
71  * \f{align*}{
72  * \frac{\partial u_\alpha}{\partial \hat{t}}
73  * + \partial_{i}
74  * \left(F^i_\alpha - v^i_g u_\alpha\right)
75  * = S_\alpha-u_\alpha\partial_i v^i_g,
76  * \f}
77  *
78  * where \f$F^i_{\alpha}\f$ are the fluxes when the mesh isn't moving,
79  * \f$v^i_g\f$ is the velocity of the mesh, \f$u_{\alpha}\f$ are the evolved
80  * variables, \f$S_{\alpha}\f$ are the source terms, and \f$\hat{t}\f$ is the
81  * time in the logical frame. For evolution equations that do not have any
82  * fluxes and only nonconservative products we evolve:
83  *
84  * \f{align*}{
85  * \frac{\partial u_\alpha}{\partial \hat{t}}
86  * +\left(B^i_{\alpha\beta}-v^i_g \delta_{\alpha\beta}
87  * \right)\partial_{i}u_\beta = S_\alpha.
88  * \f}
89  *
90  * Finally, for equations with both conservative terms and nonconservative
91  * products we use:
92  *
93  * \f{align*}{
94  * \frac{\partial u_\alpha}{\partial \hat{t}}
95  * + \partial_{i}
96  * \left(F^i_\alpha - v^i_g u_\alpha\right)
97  * +B^i_{\alpha\beta}\partial_{i}u_\beta
98  * = S_\alpha-u_\alpha\partial_i v^i_g,
99  * \f}
100  *
101  * where \f$B^i_{\alpha\beta}\f$ is the matrix for the nonconservative products.
102  *
103  * ### Volume Terms
104  *
105  * The mesh velocity is added to the flux automatically if the mesh is moving.
106  * That is,
107  *
108  * \f{align*}{
109  * F^i_{\alpha}\to F^i_{\alpha}-v^i_{g} u_{\alpha}
110  * \f}
111  *
112  * The source terms are also altered automatically by adding:
113  *
114  * \f{align*}{
115  * -u_\alpha \partial_i v^i_g,
116  * \f}
117  *
118  * For systems with equations that only contain nonconservative products, the
119  * following mesh velocity is automatically added to the time derivative:
120  *
121  * \f{align*}{
122  * v^i_g \partial_i u_\alpha,
123  * \f}
124  *
125  * \note The term is always added in the `Frame::Inertial` frame, and the plus
126  * sign arises because we add it to the time derivative.
127  *
128  * Here are examples of the `TimeDerivative` struct used to compute the volume
129  * time derivative. This struct is what the type alias
130  * `System::compute_volume_time_derivative` points to. The time derivatives are
131  * as `gsl::not_null` first, then the temporary tags as `gsl::not_null`,
132  * followed by the `argument_tags`. These type aliases are given by
133  *
134  * \snippet ComputeTimeDerivativeImpl.tpp dt_ta
135  *
136  * for the examples. For a conservative system without primitives the `apply`
137  * function would look like
138  *
139  * \snippet ComputeTimeDerivativeImpl.tpp dt_con
140  *
141  * For a nonconservative system it would be
142  *
143  * \snippet ComputeTimeDerivativeImpl.tpp dt_nc
144  *
145  * And finally, for a mixed conservative-nonconservative system with primitive
146  * variables
147  *
148  * \snippet ComputeTimeDerivativeImpl.tpp dt_mp
149  *
150  * Uses:
151  * - System:
152  * - `variables_tag`
153  * - `flux_variables`
154  * - `gradient_variables`
155  * - `compute_volume_time_derivative_terms`
156  *
157  * - DataBox:
158  * - Items in `system::compute_volume_time_derivative_terms::argument_tags`
159  * - `domain::Tags::MeshVelocity<Metavariables::volume_dim>`
160  * - `Metavariables::system::variables_tag`
161  * - `Metavariables::system::flux_variables`
162  * - `Metavariables::system::gradient_variables`
163  * - `domain::Tags::DivMeshVelocity`
164  * - `DirectionsTag`,
165  * - Required interface items for `Metavariables::system::normal_dot_fluxes`
166  *
167  * DataBox changes:
168  * - Adds: nothing
169  * - Removes: nothing
170  * - Modifies:
171  * - db::add_tag_prefix<Tags::Flux, variables_tag,
172  * tmpl::size_t<system::volume_dim>, Frame::Inertial>
173  * - `Tags::dt<system::variable_tags>`
174  * - Tags::Interface<
175  * DirectionsTag, db::add_tag_prefix<Tags::NormalDotFlux, variables_tag>>
176  * - `Tags::Mortars<typename BoundaryScheme::mortar_data_tag, VolumeDim>`
177  *
178  * ### Internal Boundary Terms
179  *
180  * Internal boundary terms are computed from the
181  * `System::boundary_correction_base` type alias. This type alias must point to
182  * a base class with `creatable_classes`. Each concrete boundary correction must
183  * specify:
184  *
185  * - type alias template `dg_package_field_tags`. These are what will be
186  * returned by `gsl::not_null` from the `dg_package_data` member function.
187  *
188  * - type alias template `dg_package_temporary_tags`. These are temporary tags
189  * that are projected to the face and then passed to the `dg_package_data`
190  * function.
191  *
192  * - type alias template `dg_package_primitive_tags`. These are the primitive
193  * variables (if any) that are projected to the face and then passed to
194  * `dg_package_data`.
195  *
196  * - type alias template `dg_package_volume_tags`. These are tags that are not
197  * projected to the interface and are retrieved directly from the `DataBox`.
198  * The equation of state for hydrodynamics systems is an example of what would
199  * be a "volume tag".
200  *
201  * A `static constexpr bool need_normal_vector` must be specified. If `true`
202  * then the normal vector is computed from the normal covector. This is
203  * currently not implemented.
204  *
205  * The `dg_package_data` function takes as arguments `gsl::not_null` of the
206  * `dg_package_data_field_tags`, then the projected evolved variables, the
207  * projected fluxes, the projected temporaries, the projected primitives, the
208  * unit normal covector, mesh velocity, normal dotted into the mesh velocity,
209  * the `volume_tags`, and finally the `dg::Formulation`. The `dg_package_data`
210  * function must compute all ingredients for the boundary correction, including
211  * mesh-velocity-corrected characteristic speeds. However, the projected fluxes
212  * passed in are \f$F^i - u v^i_g\f$ (the mesh velocity term is already
213  * included). The `dg_package_data` function must also return a `double` that is
214  * the maximum absolute characteristic speed over the entire face. This will be
215  * used for checking that the time step doesn't violate the CFL condition.
216  *
217  * Here is an example of the type aliases and `bool`:
218  *
219  * \snippet ComputeTimeDerivativeImpl.tpp bt_ta
220  *
221  * The normal vector requirement is:
222  *
223  * \snippet ComputeTimeDerivativeImpl.tpp bt_nnv
224  *
225  * For a conservative system with primitive variables and using the `TimeStepId`
226  * as a volume tag the `dg_package_data` function looks like:
227  *
228  * \snippet ComputeTimeDerivativeImpl.tpp bt_cp
229  *
230  * For a mixed conservative-nonconservative system with primitive variables and
231  * using the `TimeStepId` as a volume tag the `dg_package_data` function looks
232  * like:
233  *
234  * \snippet ComputeTimeDerivativeImpl.tpp bt_mp
235  *
236  * Uses:
237  * - System:
238  * - `boundary_correction`
239  * - `variables_tag`
240  * - `flux_variables`
241  * - `gradients_tags`
242  * - `compute_volume_time_derivative`
243  * - `has_primitive_and_conservative_vars`
244  * - `primitive_variables_tag` if system has primitive variables
245  *
246  * - DataBox:
247  * - `domain::Tags::Element<Dim>`
248  * - `domain::Tags::Mesh<Dim>`
249  * - `evolution::dg::Tags::MortarMesh<Dim>`
250  * - `evolution::dg::Tags::MortarSize<Dim>`
251  * - `evolution::dg::Tags::MortarData<Dim>`
252  * - `Tags::TimeStepId`
253  * - \code{.cpp}
254  * domain::Tags::Interface<domain::Tags::InternalDirections<Dim>,
255  * domain::Tags::Mesh<Dim - 1>>
256  * \endcode
257  * - \code{.cpp}
258  * domain::Tags::Interface<
259  * domain::Tags::InternalDirections<Dim>,
260  * ::Tags::Normalized<
261  * domain::Tags::UnnormalizedFaceNormal<Dim, Frame::Inertial>>>
262  * \endcode
263  * - \code{.cpp}
264  * domain::Tags::Interface<
265  * domain::Tags::InternalDirections<Dim>,
266  * domain::Tags::MeshVelocity<Dim, Frame::Inertial>>
267  * \endcode
268  * - `Metavariables::system::variables_tag`
269  * - `Metavariables::system::flux_variables`
270  * - `Metavariables::system::primitive_tags` if exists
271  * - `system::boundary_correction_base::dg_package_data_volume_tags`
272  *
273  * DataBox changes:
274  * - Adds: nothing
275  * - Removes: nothing
276  * - Modifies:
277  * - `evolution::dg::Tags::MortarData<Dim>`
278  */
279 template <typename Metavariables>
281  private:
282  template <typename LocalMetaVars, typename = std::void_t<>>
283  struct GetFluxInboxTag {
284  using type = tmpl::list<>;
285  };
286 
287  template <typename LocalMetaVars>
288  struct GetFluxInboxTag<LocalMetaVars,
290  using type = tmpl::list<
292  };
293 
294  public:
295  using inbox_tags = tmpl::append<
296  typename GetFluxInboxTag<Metavariables>::type,
298  Metavariables::volume_dim>>>;
299  using const_global_cache_tags = tmpl::append<
300  tmpl::conditional_t<detail::has_boundary_correction_base_v<
301  typename Metavariables::system>,
302  tmpl::list<::dg::Tags::Formulation,
304  typename Metavariables::system>>,
305  tmpl::list<::dg::Tags::Formulation>>,
306  tmpl::conditional_t<
307  Metavariables::local_time_stepping,
308  tmpl::list<
311  tmpl::list<>>>;
312 
313  template <typename DbTagsList, typename... InboxTags, typename ArrayIndex,
314  typename ActionList, typename ParallelComponent>
315  static std::tuple<db::DataBox<DbTagsList>&&> apply(
316  db::DataBox<DbTagsList>& box,
319  const ArrayIndex& /*array_index*/, ActionList /*meta*/,
320  const ParallelComponent* /*meta*/) noexcept; // NOLINT const
321 
322  private:
323  // The below functions will be removed once we've finished writing the new
324  // method of handling boundaries.
325  template <typename... NormalDotFluxTags, typename... Args>
326  static void apply_flux(
327  gsl::not_null<Variables<tmpl::list<NormalDotFluxTags...>>*> boundary_flux,
328  const Args&... boundary_variables) noexcept {
329  Metavariables::system::normal_dot_fluxes::apply(
330  make_not_null(&get<NormalDotFluxTags>(*boundary_flux))...,
331  boundary_variables...);
332  }
333 
334  template <typename DbTagsList>
335  static void boundary_terms_nonconservative_products(
336  gsl::not_null<db::DataBox<DbTagsList>*> box) noexcept;
337 
338  template <size_t VolumeDim, typename BoundaryScheme, typename DbTagsList>
339  static void fill_mortar_data_for_internal_boundaries(
340  gsl::not_null<db::DataBox<DbTagsList>*> box) noexcept;
341 
342  template <typename ParallelComponent, typename DbTagsList>
343  static void send_data_for_fluxes(
345  gsl::not_null<db::DataBox<DbTagsList>*> box) noexcept;
346 };
347 
348 template <typename Metavariables>
349 template <typename DbTagsList, typename... InboxTags, typename ArrayIndex,
350  typename ActionList, typename ParallelComponent>
353  db::DataBox<DbTagsList>& box,
356  const ArrayIndex& /*array_index*/, ActionList /*meta*/,
357  const ParallelComponent* const /*meta*/) noexcept { // NOLINT const
358  static constexpr size_t volume_dim = Metavariables::volume_dim;
359  using system = typename Metavariables::system;
360  using variables_tag = typename system::variables_tag;
361  using dt_variables_tag = db::add_tag_prefix<::Tags::dt, variables_tag>;
362  using partial_derivative_tags = typename system::gradient_variables;
363  using flux_variables = typename system::flux_variables;
364  using compute_volume_time_derivative_terms =
365  typename system::compute_volume_time_derivative_terms;
366 
367  const Mesh<volume_dim>& mesh = db::get<::domain::Tags::Mesh<volume_dim>>(box);
368  const ::dg::Formulation dg_formulation =
369  db::get<::dg::Tags::Formulation>(box);
370  ASSERT(alg::all_of(mesh.basis(),
371  [&mesh](const Spectral::Basis current_basis) noexcept {
372  return current_basis == mesh.basis(0);
373  }),
374  "An isotropic basis must be used in the evolution code. While "
375  "theoretically this restriction could be lifted, the simplification "
376  "it offers are quite substantial. Relaxing this assumption is likely "
377  "to require quite a bit of careful code refactoring and debugging.");
379  mesh.quadrature(),
380  [&mesh](const Spectral::Quadrature current_quadrature) noexcept {
381  return current_quadrature == mesh.quadrature(0);
382  }),
383  "An isotropic quadrature must be used in the evolution code. While "
384  "theoretically this restriction could be lifted, the simplification "
385  "it offers are quite substantial. Relaxing this assumption is likely "
386  "to require quite a bit of careful code refactoring and debugging.");
387 
388  // Allocate the Variables classes needed for the time derivative
389  // computation.
390  //
391  // This is factored out so that we will be able to do ADER-DG/CG where a
392  // spacetime polynomial is constructed by solving implicit equations in time
393  // using a Picard iteration. A high-order initial guess is needed to
394  // efficiently construct the ADER spacetime solution. This initial guess is
395  // obtained using continuous RK methods, and so we will want to reuse
396  // buffers. Thus, the volume_terms function returns by reference rather than
397  // by value.
398  Variables<typename compute_volume_time_derivative_terms::temporary_tags>
399  temporaries{mesh.number_of_grid_points()};
400  Variables<db::wrap_tags_in<::Tags::Flux, flux_variables,
401  tmpl::size_t<volume_dim>, Frame::Inertial>>
402  volume_fluxes{mesh.number_of_grid_points()};
403  Variables<db::wrap_tags_in<::Tags::deriv, partial_derivative_tags,
404  tmpl::size_t<volume_dim>, Frame::Inertial>>
405  partial_derivs{mesh.number_of_grid_points()};
406 
407  const Scalar<DataVector>* det_inverse_jacobian = nullptr;
408  if constexpr (tmpl::size<flux_variables>::value != 0) {
409  if (dg_formulation == ::dg::Formulation::WeakInertial) {
410  det_inverse_jacobian = &db::get<
412  }
413  }
415  tmpl::list<dt_variables_tag>,
416  typename compute_volume_time_derivative_terms::argument_tags>(
417  [&dg_formulation, &det_inverse_jacobian,
418  &div_mesh_velocity = db::get<::domain::Tags::DivMeshVelocity>(box),
419  &evolved_variables = db::get<variables_tag>(box),
420  &inertial_coordinates =
421  db::get<domain::Tags::Coordinates<volume_dim, Frame::Inertial>>(box),
422  &logical_to_inertial_inv_jacobian =
424  Frame::Inertial>>(box),
425  &mesh,
426  &mesh_velocity = db::get<::domain::Tags::MeshVelocity<volume_dim>>(box),
427  &partial_derivs, &temporaries, &volume_fluxes](
428  const gsl::not_null<Variables<
430  dt_vars_ptr,
431  const auto&... time_derivative_args) noexcept {
432  detail::volume_terms<compute_volume_time_derivative_terms>(
433  dt_vars_ptr, make_not_null(&volume_fluxes),
434  make_not_null(&partial_derivs), make_not_null(&temporaries),
435  evolved_variables, dg_formulation, mesh, inertial_coordinates,
436  logical_to_inertial_inv_jacobian, det_inverse_jacobian,
437  mesh_velocity, div_mesh_velocity, time_derivative_args...);
438  },
439  make_not_null(&box));
440 
441  if constexpr (detail::has_boundary_correction_base_v<system>) {
442  const auto& boundary_correction =
443  db::get<evolution::Tags::BoundaryCorrection<system>>(box);
444  using derived_boundary_corrections =
445  typename std::decay_t<decltype(boundary_correction)>::creatable_classes;
446 
447  const Variables<detail::get_primitive_vars_tags_from_system<system>>*
448  primitive_vars{nullptr};
449  if constexpr (system::has_primitive_and_conservative_vars) {
450  primitive_vars = &db::get<typename system::primitive_variables_tag>(box);
451  }
452 
453  static_assert(
454  tmpl::all<derived_boundary_corrections, std::is_final<tmpl::_1>>::value,
455  "All createable classes for boundary corrections must be marked "
456  "final.");
457  tmpl::for_each<derived_boundary_corrections>(
458  [&boundary_correction, &box, &partial_derivs, &primitive_vars,
459  &temporaries, &volume_fluxes](auto derived_correction_v) noexcept {
460  using DerivedCorrection =
461  tmpl::type_from<decltype(derived_correction_v)>;
462  if (typeid(boundary_correction) == typeid(DerivedCorrection)) {
463  // Compute internal boundary quantities on the mortar for sides
464  // of the element that have neighbors, i.e. they are not an
465  // external side.
466  // Note: this call mutates:
467  // - evolution::dg::Tags::NormalCovectorAndMagnitude<Dim>,
468  // - evolution::dg::Tags::MortarData<Dim>
469  detail::internal_mortar_data<system, volume_dim>(
470  make_not_null(&box),
471  dynamic_cast<const DerivedCorrection&>(boundary_correction),
472  db::get<variables_tag>(box), volume_fluxes, temporaries,
473  primitive_vars,
474  typename DerivedCorrection::dg_package_data_volume_tags{});
475 
476  if constexpr (detail::has_boundary_conditions_base_v<system>) {
477  detail::apply_boundary_conditions_on_all_external_faces<
478  system, volume_dim>(
479  make_not_null(&box),
480  dynamic_cast<const DerivedCorrection&>(boundary_correction),
481  temporaries, volume_fluxes, partial_derivs, primitive_vars);
482  } else {
483  (void)partial_derivs;
484  }
485  }
486  });
487  } else {
488  // The below if-else and fill_mortar_data_for_internal_boundaries are for
489  // compatibility with the current boundary schemes. Once the current
490  // boundary schemes are removed the code will be refactored to reflect that
491  // change.
492  if constexpr (not std::is_same_v<tmpl::list<>, flux_variables>) {
493  using flux_variables_tag = ::Tags::Variables<flux_variables>;
494  using fluxes_tag =
495  db::add_tag_prefix<::Tags::Flux, flux_variables_tag,
496  tmpl::size_t<Metavariables::volume_dim>,
498 
499  // We currently set the fluxes in the DataBox to interface with the
500  // boundary correction communication code.
501  db::mutate<fluxes_tag>(make_not_null(&box),
502  [&volume_fluxes](const auto fluxes_ptr) noexcept {
503  *fluxes_ptr = volume_fluxes;
504  });
505  } else {
506  boundary_terms_nonconservative_products(make_not_null(&box));
507  }
508 
509  // Compute internal boundary quantities
510  fill_mortar_data_for_internal_boundaries<
511  volume_dim, typename Metavariables::boundary_scheme>(
512  make_not_null(&box));
513  }
514 
515  if constexpr (Metavariables::local_time_stepping) {
516  take_step<typename Metavariables::step_choosers>(make_not_null(&box),
517  cache);
518  }
519 
520  send_data_for_fluxes<ParallelComponent>(make_not_null(&cache),
521  make_not_null(&box));
522  return {std::move(box)};
523 }
524 
525 template <typename Metavariables>
526 template <typename DbTagsList>
527 void ComputeTimeDerivative<Metavariables>::
528  boundary_terms_nonconservative_products(
529  const gsl::not_null<db::DataBox<DbTagsList>*> box) noexcept {
530  using system = typename Metavariables::system;
531  using variables_tag = typename system::variables_tag;
532  using DirectionsTag =
534 
535  using interface_normal_dot_fluxes_tag = domain::Tags::Interface<
537 
538  using interface_compute_item_argument_tags = tmpl::transform<
539  typename Metavariables::system::normal_dot_fluxes::argument_tags,
540  tmpl::bind<domain::Tags::Interface, DirectionsTag, tmpl::_1>>;
541 
543  tmpl::list<interface_normal_dot_fluxes_tag>,
544  tmpl::push_front<interface_compute_item_argument_tags, DirectionsTag>>(
545  [](const auto boundary_fluxes_ptr,
547  internal_directions,
548  const auto&... tensors) noexcept {
549  for (const auto& direction : internal_directions) {
550  apply_flux(make_not_null(&boundary_fluxes_ptr->at(direction)),
551  tensors.at(direction)...);
552  }
553  },
554  box);
555 }
556 
557 template <typename Metavariables>
558 template <size_t VolumeDim, typename BoundaryScheme, typename DbTagsList>
559 void ComputeTimeDerivative<Metavariables>::
560  fill_mortar_data_for_internal_boundaries(
561  const gsl::not_null<db::DataBox<DbTagsList>*> box) noexcept {
562  using temporal_id_tag = typename BoundaryScheme::temporal_id_tag;
563  using all_mortar_data_tag =
565  using boundary_data_computer =
566  typename BoundaryScheme::boundary_data_computer;
567 
568  // Collect data on element interfaces
569  auto boundary_data_on_interfaces =
570  interface_apply<domain::Tags::InternalDirections<VolumeDim>,
571  boundary_data_computer>(*box);
572 
573  // Project collected data to all internal mortars and store in DataBox
574  const auto& element = db::get<domain::Tags::Element<VolumeDim>>(*box);
575  const auto& face_meshes =
576  get<domain::Tags::Interface<domain::Tags::InternalDirections<VolumeDim>,
577  domain::Tags::Mesh<VolumeDim - 1>>>(*box);
578  const auto& mortar_meshes =
579  get<::Tags::Mortars<domain::Tags::Mesh<VolumeDim - 1>, VolumeDim>>(*box);
580  const auto& mortar_sizes =
581  get<::Tags::Mortars<::Tags::MortarSize<VolumeDim - 1>, VolumeDim>>(*box);
582  const auto& temporal_id = get<temporal_id_tag>(*box);
583  const auto integration_order =
584  db::get<::Tags::HistoryEvolvedVariables<>>(*box).integration_order();
585  for (const auto& [direction, neighbors] : element.neighbors()) {
586  const auto& face_mesh = face_meshes.at(direction);
587  for (const auto& neighbor : neighbors) {
588  const auto mortar_id = std::make_pair(direction, neighbor);
589  const auto& mortar_mesh = mortar_meshes.at(mortar_id);
590  const auto& mortar_size = mortar_sizes.at(mortar_id);
591 
592  // Project the data from the face to the mortar.
593  // Where no projection is necessary we `std::move` the data directly to
594  // avoid a copy. We can't move the data or modify it in-place when
595  // projecting, because in that case the face may touch more than one
596  // mortar so we need to keep the data around.
597  auto boundary_data_on_mortar =
599  ? boundary_data_on_interfaces.at(direction).project_to_mortar(
600  face_mesh, mortar_mesh, mortar_size)
601  : std::move(boundary_data_on_interfaces.at(direction));
602 
603  // Store the boundary data on this side of the mortar
604  db::mutate<all_mortar_data_tag>(
605  box, [&mortar_id, &temporal_id, &boundary_data_on_mortar,
606  &integration_order](
608  all_mortar_data) noexcept {
609  auto& mortar_data = all_mortar_data->at(mortar_id);
610  mortar_data.integration_order(integration_order);
611  mortar_data.local_insert(temporal_id,
612  std::move(boundary_data_on_mortar));
613  });
614  }
615  }
616 }
617 
618 template <typename Metavariables>
619 template <typename ParallelComponent, typename DbTagsList>
620 void ComputeTimeDerivative<Metavariables>::send_data_for_fluxes(
622  const gsl::not_null<db::DataBox<DbTagsList>*> box) noexcept {
623  using system = typename Metavariables::system;
624  constexpr size_t volume_dim = Metavariables::volume_dim;
625  auto& receiver_proxy =
626  Parallel::get_parallel_component<ParallelComponent>(*cache);
627  const auto& element = db::get<domain::Tags::Element<volume_dim>>(*box);
628 
629  if constexpr (detail::has_boundary_correction_base_v<system>) {
630  const auto& time_step_id = db::get<::Tags::TimeStepId>(*box);
631  const auto& all_mortar_data =
632  db::get<evolution::dg::Tags::MortarData<volume_dim>>(*box);
633  const auto& mortar_meshes =
634  get<evolution::dg::Tags::MortarMesh<volume_dim>>(*box);
635 
636  for (const auto& [direction, neighbors] : element.neighbors()) {
637  const auto& orientation = neighbors.orientation();
638  const auto direction_from_neighbor = orientation(direction.opposite());
639 
640  for (const auto& neighbor : neighbors) {
641  const std::pair mortar_id{direction, neighbor};
642 
643  std::pair<Mesh<volume_dim - 1>, std::vector<double>>
644  neighbor_boundary_data_on_mortar{};
645  ASSERT(time_step_id == all_mortar_data.at(mortar_id).time_step_id(),
646  "The current time step id of the volume is "
647  << time_step_id
648  << "but the time step id on the mortar with mortar id "
649  << mortar_id << " is "
650  << all_mortar_data.at(mortar_id).time_step_id());
651 
652  // Reorient the data to the neighbor orientation if necessary
653  if (LIKELY(orientation.is_aligned())) {
654  neighbor_boundary_data_on_mortar =
655  *all_mortar_data.at(mortar_id).local_mortar_data();
656  } else {
657  const auto& slice_extents = mortar_meshes.at(mortar_id).extents();
658  neighbor_boundary_data_on_mortar.first =
659  all_mortar_data.at(mortar_id).local_mortar_data()->first;
660  neighbor_boundary_data_on_mortar.second = orient_variables_on_slice(
661  all_mortar_data.at(mortar_id).local_mortar_data()->second,
662  slice_extents, direction.dimension(), orientation);
663  }
664 
665  const TimeStepId& next_time_step_id = [&box]() noexcept {
666  if (Metavariables::local_time_stepping) {
667  return db::get<::Tags::Next<::Tags::TimeStepId>>(*box);
668  } else {
669  return db::get<::Tags::TimeStepId>(*box);
670  }
671  }();
672 
675  data{neighbor_boundary_data_on_mortar.first,
676  {},
677  {std::move(neighbor_boundary_data_on_mortar.second)},
678  next_time_step_id};
679 
680  // Send mortar data (the `std::tuple` named `data`) to neighbor
683  volume_dim>>(
684  receiver_proxy[neighbor], time_step_id,
685  std::make_pair(std::pair{direction_from_neighbor, element.id()},
686  data));
687  }
688  }
689 
690  if constexpr (Metavariables::local_time_stepping) {
691  using variables_tag = typename Metavariables::system::variables_tag;
692  using dt_variables_tag = db::add_tag_prefix<::Tags::dt, variables_tag>;
693  // We assume isotropic quadrature, i.e. the quadrature is the same in
694  // all directions.
695  const bool using_gauss_points =
696  db::get<domain::Tags::Mesh<volume_dim>>(*box).quadrature() ==
697  make_array<volume_dim>(Spectral::Quadrature::Gauss);
698 
699  const Scalar<DataVector> volume_det_inv_jacobian{};
700  if (using_gauss_points) {
701  // NOLINTNEXTLINE
702  const_cast<DataVector&>(get(volume_det_inv_jacobian))
703  .set_data_ref(make_not_null(&const_cast<DataVector&>( // NOLINT
705  Frame::Logical, Frame::Inertial>>(*box)))));
706  }
707 
708  // Add face normal and Jacobian determinants to the local mortar data. We
709  // only need the Jacobians if we are using Gauss points. Then copy over
710  // into the boundary history, since that's what the LTS steppers use.
711  //
712  // The boundary history coupling computation (which computes the _lifted_
713  // boundary correction) returns a Variables<dt<EvolvedVars>> instead of
714  // using the `NormalDotNumericalFlux` prefix tag. This is because the
715  // returned quantity is more a `dt` quantity than a
716  // `NormalDotNormalDotFlux` since it's been lifted to the volume.
718  const auto integration_order =
719  db::get<::Tags::HistoryEvolvedVariables<>>(*box).integration_order();
720  db::mutate<evolution::dg::Tags::MortarData<volume_dim>,
722  volume_dim, typename dt_variables_tag::type>>(
723  box,
724  [&element, integration_order, &time_step_id, using_gauss_points,
725  &volume_det_inv_jacobian](
726  const gsl::not_null<
728  boost::hash<Key>>*>
729  mortar_data,
730  const gsl::not_null<
731  std::unordered_map<Key,
735  typename dt_variables_tag::type>,
736  boost::hash<Key>>*>
737  boundary_data_history,
738  const Mesh<volume_dim>& volume_mesh,
739  const DirectionMap<
740  volume_dim,
741  std::optional<Variables<tmpl::list<
744  normal_covector_and_magnitude) noexcept {
745  Scalar<DataVector> volume_det_jacobian{};
746  Scalar<DataVector> face_det_jacobian{};
747  if (using_gauss_points) {
748  get(volume_det_jacobian) = 1.0 / get(volume_det_inv_jacobian);
749  }
750  for (const auto& [direction, neighbors_in_direction] :
751  element.neighbors()) {
752  // We can perform projections once for all neighbors in the
753  // direction because we care about the _face_ mesh, not the mortar
754  // mesh.
755  ASSERT(
756  normal_covector_and_magnitude.at(direction).has_value(),
757  "The normal covector and magnitude have not been computed.");
758  const Scalar<DataVector>& face_normal_magnitude =
759  get<evolution::dg::Tags::MagnitudeOfNormal>(
760  *normal_covector_and_magnitude.at(direction));
761  if (using_gauss_points) {
762  const Matrix identity{};
763  auto interpolation_matrices =
764  make_array<Metavariables::volume_dim>(std::cref(identity));
765  const std::pair<Matrix, Matrix>& matrices =
767  volume_mesh.slice_through(direction.dimension()));
768  gsl::at(interpolation_matrices, direction.dimension()) =
769  direction.side() == Side::Upper ? matrices.second
770  : matrices.first;
771  if (get(face_det_jacobian).size() !=
772  get(face_normal_magnitude).size()) {
773  get(face_det_jacobian) =
774  DataVector{get(face_normal_magnitude).size()};
775  }
776  apply_matrices(make_not_null(&get(face_det_jacobian)),
777  interpolation_matrices, get(volume_det_jacobian),
778  volume_mesh.extents());
779  }
780 
781  for (const auto& neighbor : neighbors_in_direction) {
782  const std::pair mortar_id{direction, neighbor};
783  if (using_gauss_points) {
784  mortar_data->at(mortar_id).insert_local_geometric_quantities(
785  volume_det_inv_jacobian, face_det_jacobian,
786  face_normal_magnitude);
787  } else {
788  mortar_data->at(mortar_id).insert_local_face_normal_magnitude(
789  face_normal_magnitude);
790  }
791  ASSERT(boundary_data_history->count(mortar_id) != 0,
792  "Could not insert the mortar data for "
793  << mortar_id
794  << " because the unordered map has not been "
795  "initialized "
796  "to have the mortar id.");
797  boundary_data_history->at(mortar_id).local_insert(
798  time_step_id, std::move(mortar_data->at(mortar_id)));
799  boundary_data_history->at(mortar_id).integration_order(
800  integration_order);
801  mortar_data->at(mortar_id) =
802  MortarData<Metavariables::volume_dim>{};
803  }
804  }
805  },
806  db::get<domain::Tags::Mesh<volume_dim>>(*box),
807  db::get<evolution::dg::Tags::NormalCovectorAndMagnitude<volume_dim>>(
808  *box));
809  }
810  } else {
811  using BoundaryScheme = typename Metavariables::boundary_scheme;
812  using temporal_id_tag = typename BoundaryScheme::temporal_id_tag;
813  using receive_temporal_id_tag =
814  typename BoundaryScheme::receive_temporal_id_tag;
815  using fluxes_inbox_tag = ::dg::FluxesInboxTag<BoundaryScheme>;
816  using all_mortar_data_tag =
818 
819  const auto& all_mortar_data = get<all_mortar_data_tag>(*box);
820  const auto& temporal_id = db::get<temporal_id_tag>(*box);
821  const auto& receive_temporal_id = db::get<receive_temporal_id_tag>(*box);
822  const auto& mortar_meshes = db::get<
823  ::Tags::Mortars<domain::Tags::Mesh<volume_dim - 1>, volume_dim>>(*box);
824 
825  // Iterate over neighbors
826  for (const auto& direction_and_neighbors : element.neighbors()) {
827  const auto& direction = direction_and_neighbors.first;
828  const size_t dimension = direction.dimension();
829  const auto& neighbors_in_direction = direction_and_neighbors.second;
830  const auto& orientation = neighbors_in_direction.orientation();
831  const auto direction_from_neighbor = orientation(direction.opposite());
832 
833  for (const auto& neighbor : neighbors_in_direction) {
834  const ::dg::MortarId<volume_dim> mortar_id{direction, neighbor};
835 
836  // Make a copy of the local boundary data on the mortar to send to the
837  // neighbor
838  ASSERT(all_mortar_data.find(mortar_id) != all_mortar_data.end(),
839  "Mortar data on mortar "
840  << mortar_id
841  << " not available for sending. Did you forget to collect "
842  "the data on mortars?");
843  auto neighbor_boundary_data_on_mortar =
844  all_mortar_data.at(mortar_id).local_data(temporal_id);
845 
846  // Reorient the data to the neighbor orientation if necessary
847  if (not orientation.is_aligned()) {
848  neighbor_boundary_data_on_mortar.orient_on_slice(
849  mortar_meshes.at(mortar_id).extents(), dimension, orientation);
850  }
851 
852  // Send mortar data (named 'neighbor_boundary_data_on_mortar') to
853  // neighbor
854  Parallel::receive_data<fluxes_inbox_tag>(
855  receiver_proxy[neighbor], temporal_id,
856  std::make_pair(
857  ::dg::MortarId<volume_dim>{direction_from_neighbor,
858  element.id()},
859  std::make_pair(receive_temporal_id,
860  std::move(neighbor_boundary_data_on_mortar))));
861  }
862  }
863  }
864 }
865 } // 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
Spectral::needs_projection
bool needs_projection(const Mesh< Dim > &mesh1, const Mesh< Dim > &mesh2, const std::array< ChildSize, Dim > &child_sizes) noexcept
Determine whether data needs to be projected between a child mesh and its parent mesh....
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
get
constexpr Tag::type & get(Variables< TagList > &v) noexcept
Return Tag::type pointing into the contiguous array.
Definition: Variables.hpp:660
Parallel::GlobalCache
Definition: ElementReceiveInterpPoints.hpp:15
unordered_set
std::pair
GlobalCache.hpp
Tags::StepController
Tag for a StepController.
Definition: Tags.hpp:242
Tags::MortarSize
Definition: Tags.hpp:48
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
std::void_t
Tags::Variables
Definition: VariablesTag.hpp:21
domain::Tags::Mesh
The computational grid of the Element in the DataBox.
Definition: Tags.hpp:107
tuple
db::get
const auto & get(const DataBox< TagList > &box) noexcept
Retrieve the item with tag Tag from the DataBox.
Definition: DataBox.hpp:786
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::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:50
evolution::dg::MortarData
Data on the mortar used to compute the boundary correction for the DG scheme.
Definition: MortarData.hpp:50
dg::FluxesInboxTag
The inbox tag for flux communication.
Definition: FluxCommunication.hpp:35
evolution::dg::Actions::ComputeTimeDerivative
Computes the time derivative for a DG time step.
Definition: ComputeTimeDerivative.hpp:280
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:43
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
Tags::StepChoosers
Tag for a vector of StepChoosers.
Definition: Tags.hpp:229
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:48
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:380
Scalar
Tensor< T, Symmetry<>, index_list<> > Scalar
Definition: TypeAliases.hpp:21
std::is_final
Gsl.hpp
domain::Tags::Interface
Tag which is either a SimpleTag for quantities on an interface, base tag to a compute item which acts...
Definition: Tags.hpp:366
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:1193
dg::mortar_mesh
Mesh< Dim > mortar_mesh(const Mesh< Dim > &face_mesh1, const Mesh< Dim > &face_mesh2) noexcept
Tags::TimeStepper
Tag for a TimeStepper of type StepperType.
Definition: Tags.hpp:206
Tensor.hpp
TimeSteppers::BoundaryHistory
Definition: BoundaryHistory.hpp:31
domain::Tags::InternalDirections
Definition: Tags.hpp:270
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
dg::mortar_size
MortarSize< Dim - 1 > mortar_size(const ElementId< Dim > &self, const ElementId< Dim > &neighbor, size_t dimension, const OrientationMap< Dim > &orientation) noexcept
TMPL.hpp
Mesh.hpp
gsl::not_null
Require a pointer to not be a nullptr
Definition: ReadSpecThirdOrderPiecewisePolynomial.hpp:13
Tags::Mortars
Definition: Tags.hpp:37