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/Tags.hpp"
24 #include "Domain/TagsTimeDependent.hpp"
25 #include "Evolution/BoundaryCorrectionTags.hpp"
26 #include "Evolution/DiscontinuousGalerkin/Actions/NormalCovectorAndMagnitude.hpp"
27 #include "Evolution/DiscontinuousGalerkin/InboxTags.hpp"
28 #include "Evolution/DiscontinuousGalerkin/MortarData.hpp"
29 #include "Evolution/DiscontinuousGalerkin/MortarTags.hpp"
30 #include "Evolution/DiscontinuousGalerkin/NormalVectorTags.hpp"
31 #include "Evolution/DiscontinuousGalerkin/ProjectToBoundary.hpp"
32 #include "NumericalAlgorithms/DiscontinuousGalerkin/Formulation.hpp"
33 #include "NumericalAlgorithms/DiscontinuousGalerkin/MetricIdentityJacobian.hpp"
34 #include "NumericalAlgorithms/DiscontinuousGalerkin/MortarHelpers.hpp"
35 #include "NumericalAlgorithms/DiscontinuousGalerkin/Tags.hpp"
36 #include "NumericalAlgorithms/DiscontinuousGalerkin/Tags/Formulation.hpp"
39 #include "NumericalAlgorithms/LinearOperators/WeakDivergence.hpp"
41 #include "NumericalAlgorithms/Spectral/Projection.hpp"
42 #include "Parallel/GlobalCache.hpp"
43 #include "ParallelAlgorithms/DiscontinuousGalerkin/FluxCommunication.hpp"
44 #include "Utilities/Gsl.hpp"
45 #include "Utilities/TMPL.hpp"
46 #include "Utilities/TypeTraits/CreateHasTypeAlias.hpp"
47 
48 /// \cond
49 namespace tuples {
50 template <typename...>
51 class TaggedTuple;
52 } // namespace tuples
53 /// \endcond
54 
55 namespace evolution::dg::Actions {
56 namespace detail {
57 CREATE_HAS_TYPE_ALIAS(boundary_correction)
58 CREATE_HAS_TYPE_ALIAS_V(boundary_correction)
59 
60 template <bool HasInverseSpatialMetricTag = false>
61 struct inverse_spatial_metric_tag_impl {
62  template <typename System>
63  using f = tmpl::list<>;
64 };
65 
66 template <>
67 struct inverse_spatial_metric_tag_impl<true> {
68  template <typename System>
69  using f = tmpl::list<typename System::inverse_spatial_metric_tag>;
70 };
71 
72 template <typename System>
73 using inverse_spatial_metric_tag = typename inverse_spatial_metric_tag_impl<
74  has_inverse_spatial_metric_tag_v<System>>::template f<System>;
75 
76 template <bool HasPrimitiveVars = false>
77 struct get_primitive_vars {
78  template <typename BoundaryCorrection>
79  using f = tmpl::list<>;
80 };
81 
82 template <>
83 struct get_primitive_vars<true> {
84  template <typename BoundaryCorrection>
85  using f = typename BoundaryCorrection::dg_package_data_primitive_tags;
86 };
87 
88 // Helper function to get parameter packs so we can forward `Tensor`s instead
89 // of `Variables` to the boundary corrections. Returns the maximum absolute
90 // char speed on the face, which can be used for setting or monitoring the CFL
91 // without having to compute the speeds for each dimension in the volume.
92 // Whether using only the face speeds is accurate enough to guarantee
93 // stability is yet to be determined. However, if the CFL condition is
94 // violated on the boundaries we are definitely in trouble, so it can at least
95 // be used as a cheap diagnostic.
96 template <typename System, typename BoundaryCorrection,
97  typename... PackagedFieldTags, typename... ProjectedFieldTags,
98  typename... ProjectedFieldTagsForCorrection, size_t Dim,
99  typename DbTagsList, typename... VolumeTags>
100 double dg_package_data(
101  const gsl::not_null<Variables<tmpl::list<PackagedFieldTags...>>*>
102  packaged_data,
103  const BoundaryCorrection& boundary_correction,
104  const Variables<tmpl::list<ProjectedFieldTags...>>& projected_fields,
105  const tnsr::i<DataVector, Dim, Frame::Inertial>& unit_normal_covector,
106  const std::optional<tnsr::I<DataVector, Dim, Frame::Inertial>>&
107  mesh_velocity,
108  const db::DataBox<DbTagsList>& box, tmpl::list<VolumeTags...> /*meta*/,
109  tmpl::list<ProjectedFieldTagsForCorrection...> /*meta*/) noexcept {
110  std::optional<Scalar<DataVector>> normal_dot_mesh_velocity{};
111  if (mesh_velocity.has_value()) {
112  normal_dot_mesh_velocity =
113  dot_product(*mesh_velocity, unit_normal_covector);
114  }
115 
116  if constexpr (evolution::dg::Actions::detail::
117  has_inverse_spatial_metric_tag_v<System>) {
118  return boundary_correction.dg_package_data(
119  make_not_null(&get<PackagedFieldTags>(*packaged_data))...,
120  get<ProjectedFieldTagsForCorrection>(projected_fields)...,
121  unit_normal_covector,
122  get<evolution::dg::Actions::detail::NormalVector<Dim>>(
123  projected_fields),
124  mesh_velocity, normal_dot_mesh_velocity, db::get<VolumeTags>(box)...);
125  } else {
126  return boundary_correction.dg_package_data(
127  make_not_null(&get<PackagedFieldTags>(*packaged_data))...,
128  get<ProjectedFieldTagsForCorrection>(projected_fields)...,
129  unit_normal_covector, mesh_velocity, normal_dot_mesh_velocity,
130  db::get<VolumeTags>(box)...);
131  }
132 }
133 } // namespace detail
134 
135 /*!
136  * \brief Computes the time derivative for a DG time step.
137  *
138  * Computes the volume fluxes, the divergence of the fluxes and all additional
139  * interior contributions to the time derivatives (both nonconservative products
140  * and source terms). The internal mortar data is also computed.
141  *
142  * The general first-order hyperbolic evolution equation solved for conservative
143  * systems is:
144  *
145  * \f{align*}{
146  * \frac{\partial u_\alpha}{\partial \hat{t}}
147  * + \partial_{i}
148  * \left(F^i_\alpha - v^i_g u_\alpha\right)
149  * = S_\alpha-u_\alpha\partial_i v^i_g,
150  * \f}
151  *
152  * where \f$F^i_{\alpha}\f$ are the fluxes when the mesh isn't moving,
153  * \f$v^i_g\f$ is the velocity of the mesh, \f$u_{\alpha}\f$ are the evolved
154  * variables, \f$S_{\alpha}\f$ are the source terms, and \f$\hat{t}\f$ is the
155  * time in the logical frame. For evolution equations that do not have any
156  * fluxes and only nonconservative products we evolve:
157  *
158  * \f{align*}{
159  * \frac{\partial u_\alpha}{\partial \hat{t}}
160  * +\left(B^i_{\alpha\beta}-v^i_g \delta_{\alpha\beta}
161  * \right)\partial_{i}u_\beta = S_\alpha.
162  * \f}
163  *
164  * Finally, for equations with both conservative terms and nonconservative
165  * products we use:
166  *
167  * \f{align*}{
168  * \frac{\partial u_\alpha}{\partial \hat{t}}
169  * + \partial_{i}
170  * \left(F^i_\alpha - v^i_g u_\alpha\right)
171  * +B^i_{\alpha\beta}\partial_{i}u_\beta
172  * = S_\alpha-u_\alpha\partial_i v^i_g,
173  * \f}
174  *
175  * where \f$B^i_{\alpha\beta}\f$ is the matrix for the nonconservative products.
176  *
177  * ### Volume Terms
178  *
179  * The mesh velocity is added to the flux automatically if the mesh is moving.
180  * That is,
181  *
182  * \f{align*}{
183  * F^i_{\alpha}\to F^i_{\alpha}-v^i_{g} u_{\alpha}
184  * \f}
185  *
186  * The source terms are also altered automatically by adding:
187  *
188  * \f{align*}{
189  * -u_\alpha \partial_i v^i_g,
190  * \f}
191  *
192  * For systems with equations that only contain nonconservative products, the
193  * following mesh velocity is automatically added to the time derivative:
194  *
195  * \f{align*}{
196  * v^i_g \partial_i u_\alpha,
197  * \f}
198  *
199  * \note The term is always added in the `Frame::Inertial` frame, and the plus
200  * sign arises because we add it to the time derivative.
201  *
202  * Here are examples of the `TimeDerivative` struct used to compute the volume
203  * time derivative. This struct is what the type alias
204  * `System::compute_volume_time_derivative` points to. The time derivatives are
205  * as `gsl::not_null` first, then the temporary tags as `gsl::not_null`,
206  * followed by the `argument_tags`. These type aliases are given by
207  *
208  * \snippet ComputeTimeDerivativeImpl.tpp dt_ta
209  *
210  * for the examples. For a conservative system without primitives the `apply`
211  * function would look like
212  *
213  * \snippet ComputeTimeDerivativeImpl.tpp dt_con
214  *
215  * For a nonconservative system it would be
216  *
217  * \snippet ComputeTimeDerivativeImpl.tpp dt_nc
218  *
219  * And finally, for a mixed conservative-nonconservative system with primitive
220  * variables
221  *
222  * \snippet ComputeTimeDerivativeImpl.tpp dt_mp
223  *
224  * Uses:
225  * - System:
226  * - `variables_tag`
227  * - `flux_variables`
228  * - `gradient_variables`
229  * - `compute_volume_time_derivative_terms`
230  *
231  * - DataBox:
232  * - Items in `system::compute_volume_time_derivative_terms::argument_tags`
233  * - `domain::Tags::MeshVelocity<Metavariables::volume_dim>`
234  * - `Metavariables::system::variables_tag`
235  * - `Metavariables::system::flux_variables`
236  * - `Metavariables::system::gradient_variables`
237  * - `domain::Tags::DivMeshVelocity`
238  * - `DirectionsTag`,
239  * - Required interface items for `Metavariables::system::normal_dot_fluxes`
240  *
241  * DataBox changes:
242  * - Adds: nothing
243  * - Removes: nothing
244  * - Modifies:
245  * - db::add_tag_prefix<Tags::Flux, variables_tag,
246  * tmpl::size_t<system::volume_dim>, Frame::Inertial>
247  * - `Tags::dt<system::variable_tags>`
248  * - Tags::Interface<
249  * DirectionsTag, db::add_tag_prefix<Tags::NormalDotFlux, variables_tag>>
250  * - `Tags::Mortars<typename BoundaryScheme::mortar_data_tag, VolumeDim>`
251  *
252  * ### Internal Boundary Terms
253  *
254  * Internal boundary terms are computed from the `System::boundary_correction`
255  * type alias. This type alias must point to a base class with
256  * `creatable_classes`. Each concrete boundary correction must specify:
257  *
258  * - type alias template `dg_package_field_tags`. These are what will be
259  * returned by `gsl::not_null` from the `dg_package_data` member function.
260  *
261  * - type alias template `dg_package_temporary_tags`. These are temporary tags
262  * that are projected to the face and then passed to the `dg_package_data`
263  * function.
264  *
265  * - type alias template `dg_package_primitive_tags`. These are the primitive
266  * variables (if any) that are projected to the face and then passed to
267  * `dg_package_data`.
268  *
269  * - type alias template `dg_package_volume_tags`. These are tags that are not
270  * projected to the interface and are retrieved directly from the `DataBox`.
271  * The equation of state for hydrodynamics systems is an example of what would
272  * be a "volume tag".
273  *
274  * A `static constexpr bool need_normal_vector` must be specified. If `true`
275  * then the normal vector is computed from the normal covector. This is
276  * currently not implemented.
277  *
278  * The `dg_package_data` function takes as arguments `gsl::not_null` of the
279  * `dg_package_data_field_tags`, then the projected evolved variables, the
280  * projected fluxes, the projected temporaries, the projected primitives, the
281  * unit normal covector, mesh velocity, normal dotted into the mesh velocity,
282  * the `volume_tags`, and finally the `dg::Formulation`. The `dg_package_data`
283  * function must compute all ingredients for the boundary correction, including
284  * mesh-velocity-corrected characteristic speeds. However, the projected fluxes
285  * passed in are \f$F^i - u v^i_g\f$ (the mesh velocity term is already
286  * included). The `dg_package_data` function must also return a `double` that is
287  * the maximum absolute characteristic speed over the entire face. This will be
288  * used for checking that the time step doesn't violate the CFL condition.
289  *
290  * Here is an example of the type aliases and `bool`:
291  *
292  * \snippet ComputeTimeDerivativeImpl.tpp bt_ta
293  *
294  * The normal vector requirement is:
295  *
296  * \snippet ComputeTimeDerivativeImpl.tpp bt_nnv
297  *
298  * For a conservative system with primitive variables and using the `TimeStepId`
299  * as a volume tag the `dg_package_data` function looks like:
300  *
301  * \snippet ComputeTimeDerivativeImpl.tpp bt_cp
302  *
303  * For a mixed conservative-nonconservative system with primitive variables and
304  * using the `TimeStepId` as a volume tag the `dg_package_data` function looks
305  * like:
306  *
307  * \snippet ComputeTimeDerivativeImpl.tpp bt_mp
308  *
309  * Uses:
310  * - System:
311  * - `boundary_correction`
312  * - `variables_tag`
313  * - `flux_variables`
314  * - `gradients_tags`
315  * - `compute_volume_time_derivative`
316  * - `has_primitive_and_conservative_vars`
317  * - `primitive_variables_tag` if system has primitive variables
318  *
319  * - DataBox:
320  * - `domain::Tags::Element<Dim>`
321  * - `domain::Tags::Mesh<Dim>`
322  * - `evolution::dg::Tags::MortarMesh<Dim>`
323  * - `evolution::dg::Tags::MortarSize<Dim>`
324  * - `evolution::dg::Tags::MortarData<Dim>`
325  * - `Tags::TimeStepId`
326  * - \code{.cpp}
327  * domain::Tags::Interface<domain::Tags::InternalDirections<Dim>,
328  * domain::Tags::Mesh<Dim - 1>>
329  * \endcode
330  * - \code{.cpp}
331  * domain::Tags::Interface<
332  * domain::Tags::InternalDirections<Dim>,
333  * ::Tags::Normalized<
334  * domain::Tags::UnnormalizedFaceNormal<Dim, Frame::Inertial>>>
335  * \endcode
336  * - \code{.cpp}
337  * domain::Tags::Interface<
338  * domain::Tags::InternalDirections<Dim>,
339  * domain::Tags::MeshVelocity<Dim, Frame::Inertial>>
340  * \endcode
341  * - `Metavariables::system::variables_tag`
342  * - `Metavariables::system::flux_variables`
343  * - `Metavariables::system::primitive_tags` if exists
344  * - `Metavariables::system::boundary_correction::dg_package_data_volume_tags`
345  *
346  * DataBox changes:
347  * - Adds: nothing
348  * - Removes: nothing
349  * - Modifies:
350  * - `evolution::dg::Tags::MortarData<Dim>`
351  */
352 template <typename Metavariables>
354  public:
355  using inbox_tags =
356  tmpl::list<::dg::FluxesInboxTag<typename Metavariables::boundary_scheme>,
358  Metavariables::volume_dim>>;
359  using const_global_cache_tags = tmpl::conditional_t<
360  detail::has_boundary_correction_v<typename Metavariables::system>,
362  typename Metavariables::system>>,
363  tmpl::list<::dg::Tags::Formulation>>;
364 
365  template <typename DbTagsList, typename... InboxTags, typename ArrayIndex,
366  typename ActionList, typename ParallelComponent>
367  static std::tuple<db::DataBox<DbTagsList>&&> apply(
368  db::DataBox<DbTagsList>& box,
371  const ArrayIndex& /*array_index*/, ActionList /*meta*/,
372  const ParallelComponent* /*meta*/) noexcept; // NOLINT const
373 
374  private:
375  /*
376  * Computes the volume terms for a discontinuous Galerkin scheme.
377  *
378  * The function does the following (in order):
379  *
380  * 1. Compute the partial derivatives of the `System::gradient_variables`.
381  *
382  * The partial derivatives are needed in the nonconservative product terms
383  * of the evolution equations. Any variable whose evolution equation does
384  * not contain a flux must contain a nonconservative product and the
385  * variable must be listed in the `System::gradient_variables` type alias.
386  * The partial derivatives are also needed for adding the moving mesh terms
387  * to the equations that do not have a flux term.
388  *
389  * 2. The volume time derivatives are calculated from
390  * `System::compute_volume_time_derivative_terms`
391  *
392  * The source terms and nonconservative products are contributed directly
393  * to the `dt_vars` arguments passed to the time derivative function, while
394  * the volume fluxes are computed into the `volume_fluxes` arguments. The
395  * divergence of the volume fluxes will be computed and added to the time
396  * derivatives later in the function.
397  *
398  * 3. If the mesh is moving the appropriate mesh velocity terms are added to
399  * the equations.
400  *
401  * For equations with fluxes this means that \f$-v^i_g u_\alpha\f$ is
402  * added to the fluxes and \f$-u_\alpha \partial_i v^i_g\f$ is added
403  * to the time derivatives. For equations without fluxes
404  * \f$v^i\partial_i u_\alpha\f$ is added to the time derivatives.
405  *
406  * 4. Compute flux divergence contribution and add it to the time derivatives.
407  *
408  * Either the weak or strong form can be used. Currently only the strong
409  * form is coded, but adding the weak form is quite easy.
410  *
411  * Note that the computation of the flux divergence and adding that to the
412  * time derivative must be done *after* the mesh velocity is subtracted
413  * from the fluxes.
414  */
415  template <size_t Dim, typename ComputeVolumeTimeDerivatives,
416  typename DbTagsList, typename... VariablesTags,
417  typename... TimeDerivativeArgumentTags,
418  typename... PartialDerivTags, typename... FluxVariablesTags,
419  typename... TemporaryTags>
420  static void volume_terms(
421  gsl::not_null<db::DataBox<DbTagsList>*> box,
422  gsl::not_null<Variables<tmpl::list<FluxVariablesTags...>>*> volume_fluxes,
423  gsl::not_null<Variables<tmpl::list<PartialDerivTags...>>*> partial_derivs,
424  gsl::not_null<Variables<tmpl::list<TemporaryTags...>>*> temporaries,
425  ::dg::Formulation dg_formulation,
426  tmpl::list<VariablesTags...> /*variables_tags*/,
427  tmpl::list<TimeDerivativeArgumentTags...> /*meta*/) noexcept;
428 
429  template <size_t Dim, typename BoundaryCorrection, typename TemporaryTags,
430  typename DbTagsList, typename VariablesTags,
431  typename... PackageDataVolumeTags>
432  static void compute_internal_mortar_data(
433  gsl::not_null<db::DataBox<DbTagsList>*> box,
434  const BoundaryCorrection& boundary_correction,
435  const Variables<VariablesTags>& volume_evolved_vars,
436  const Variables<db::wrap_tags_in<
437  ::Tags::Flux, typename Metavariables::system::flux_variables,
438  tmpl::size_t<Dim>, Frame::Inertial>>& volume_fluxes,
439  const Variables<TemporaryTags>& volume_temporaries) noexcept;
440 
441  // The below functions will be removed once we've finished writing the new
442  // method of handling boundaries.
443  template <typename... NormalDotFluxTags, typename... Args>
444  static void apply_flux(
445  gsl::not_null<Variables<tmpl::list<NormalDotFluxTags...>>*> boundary_flux,
446  const Args&... boundary_variables) noexcept {
447  Metavariables::system::normal_dot_fluxes::apply(
448  make_not_null(&get<NormalDotFluxTags>(*boundary_flux))...,
449  boundary_variables...);
450  }
451 
452  template <typename DbTagsList>
453  static void boundary_terms_nonconservative_products(
454  gsl::not_null<db::DataBox<DbTagsList>*> box) noexcept;
455 
456  template <size_t VolumeDim, typename BoundaryScheme, typename DbTagsList>
457  static void fill_mortar_data_for_internal_boundaries(
458  gsl::not_null<db::DataBox<DbTagsList>*> box) noexcept;
459 
460  template <typename ParallelComponent, typename DbTagsList>
461  static void send_data_for_fluxes(
463  const db::DataBox<DbTagsList>& box) noexcept;
464 };
465 
466 template <typename Metavariables>
467 template <typename DbTagsList, typename... InboxTags, typename ArrayIndex,
468  typename ActionList, typename ParallelComponent>
471  db::DataBox<DbTagsList>& box,
474  const ArrayIndex& /*array_index*/, ActionList /*meta*/,
475  const ParallelComponent* const /*meta*/) noexcept { // NOLINT const
476  static constexpr size_t volume_dim = Metavariables::volume_dim;
477  using system = typename Metavariables::system;
478  using variables_tag = typename system::variables_tag;
479  using variables_tags = typename variables_tag::tags_list;
480  using partial_derivative_tags = typename system::gradient_variables;
481  using flux_variables = typename system::flux_variables;
482  using compute_volume_time_derivative_terms =
483  typename system::compute_volume_time_derivative_terms;
484 
485  const Mesh<volume_dim>& mesh = db::get<::domain::Tags::Mesh<volume_dim>>(box);
486  // Allocate the Variables classes needed for the time derivative
487  // computation.
488  //
489  // This is factored out so that we will be able to do ADER-DG/CG where a
490  // spacetime polynomial is constructed by solving implicit equations in time
491  // using a Picard iteration. A high-order initial guess is needed to
492  // efficiently construct the ADER spacetime solution. This initial guess is
493  // obtained using continuous RK methods, and so we will want to reuse
494  // buffers. Thus, the volume_terms function returns by reference rather than
495  // by value.
496  Variables<typename compute_volume_time_derivative_terms::temporary_tags>
497  temporaries{mesh.number_of_grid_points()};
498  Variables<db::wrap_tags_in<::Tags::Flux, flux_variables,
499  tmpl::size_t<volume_dim>, Frame::Inertial>>
500  volume_fluxes{mesh.number_of_grid_points()};
501  Variables<db::wrap_tags_in<::Tags::deriv, partial_derivative_tags,
502  tmpl::size_t<volume_dim>, Frame::Inertial>>
503  partial_derivs{mesh.number_of_grid_points()};
504 
505  volume_terms<volume_dim, compute_volume_time_derivative_terms>(
506  make_not_null(&box), make_not_null(&volume_fluxes),
507  make_not_null(&partial_derivs), make_not_null(&temporaries),
508  db::get<::dg::Tags::Formulation>(box), variables_tags{},
509  typename compute_volume_time_derivative_terms::argument_tags{});
510 
511  // The below if-else and fill_mortar_data_for_internal_boundaries are for
512  // compatibility with the current boundary schemes. Once the current
513  // boundary schemes are removed the code will be refactored to reflect that
514  // change.
515  if constexpr (not std::is_same_v<tmpl::list<>, flux_variables>) {
516  using flux_variables_tag = ::Tags::Variables<flux_variables>;
517  using fluxes_tag =
518  db::add_tag_prefix<::Tags::Flux, flux_variables_tag,
519  tmpl::size_t<Metavariables::volume_dim>,
521 
522  // We currently set the fluxes in the DataBox to interface with the
523  // boundary correction communication code.
524  db::mutate<fluxes_tag>(make_not_null(&box),
525  [&volume_fluxes](const auto fluxes_ptr) noexcept {
526  *fluxes_ptr = volume_fluxes;
527  });
528  } else {
529  boundary_terms_nonconservative_products(make_not_null(&box));
530  }
531 
532  if constexpr (detail::has_boundary_correction_v<system>) {
533  const auto& boundary_correction =
534  db::get<evolution::Tags::BoundaryCorrection<system>>(box);
535  using derived_boundary_corrections =
536  typename std::decay_t<decltype(boundary_correction)>::creatable_classes;
537 
538  static_assert(
539  tmpl::all<derived_boundary_corrections, std::is_final<tmpl::_1>>::value,
540  "All createable classes for boundary corrections must be marked "
541  "final.");
542  tmpl::for_each<derived_boundary_corrections>(
543  [&boundary_correction, &box, &temporaries,
544  &volume_fluxes](auto derived_correction_v) noexcept {
545  using DerivedCorrection =
546  tmpl::type_from<decltype(derived_correction_v)>;
547  if (typeid(boundary_correction) == typeid(DerivedCorrection)) {
548  // Compute internal boundary quantities on the mortar for sides
549  // of the element that have neighbors, i.e. they are not an
550  // external side.
551  compute_internal_mortar_data<volume_dim>(
552  make_not_null(&box),
553  dynamic_cast<const DerivedCorrection&>(boundary_correction),
554  db::get<variables_tag>(box), volume_fluxes, temporaries);
555  }
556  });
557  } else {
558  // Compute internal boundary quantities
559  fill_mortar_data_for_internal_boundaries<
560  volume_dim, typename Metavariables::boundary_scheme>(
561  make_not_null(&box));
562  }
563 
564  send_data_for_fluxes<ParallelComponent>(make_not_null(&cache), box);
565  return {std::move(box)};
566 }
567 
568 template <typename Metavariables>
569 template <size_t Dim, typename ComputeVolumeTimeDerivatives,
570  typename DbTagsList, typename... VariablesTags,
571  typename... TimeDerivativeArgumentTags, typename... PartialDerivTags,
572  typename... FluxVariablesTags, typename... TemporaryTags>
573 void ComputeTimeDerivative<Metavariables>::volume_terms(
574  const gsl::not_null<db::DataBox<DbTagsList>*> box,
575  const gsl::not_null<Variables<tmpl::list<FluxVariablesTags...>>*>
576  volume_fluxes,
577  const gsl::not_null<Variables<tmpl::list<PartialDerivTags...>>*>
578  partial_derivs,
579  const gsl::not_null<Variables<tmpl::list<TemporaryTags...>>*> temporaries,
580  const ::dg::Formulation dg_formulation,
581 
582  tmpl::list<VariablesTags...> /*variables_tags*/,
583  tmpl::list<TimeDerivativeArgumentTags...> /*meta*/) noexcept {
584  static constexpr bool has_partial_derivs = sizeof...(PartialDerivTags) != 0;
585  static constexpr bool has_fluxes = sizeof...(FluxVariablesTags) != 0;
586  static_assert(
587  has_fluxes or has_partial_derivs,
588  "Must have either fluxes or partial derivatives in a "
589  "DG evolution scheme. This means the evolution system struct (usually in "
590  "Evolution/Systems/YourSystem/System.hpp) being used does not specify "
591  "any flux_variables or gradient_variables. Make sure the type aliases "
592  "are defined, and that at least one of them is a non-empty list of "
593  "tags.");
594 
595  using system = typename Metavariables::system;
596  using variables_tag = typename system::variables_tag;
597  using variables_tags = typename variables_tag::tags_list;
598  using partial_derivative_tags = typename system::gradient_variables;
599  using flux_variables = typename system::flux_variables;
600 
601  const Mesh<Dim>& mesh = db::get<::domain::Tags::Mesh<Dim>>(*box);
602  const auto& logical_to_inertial_inverse_jacobian = db::get<
604  *box);
605  const auto& evolved_vars = db::get<variables_tag>(*box);
606 
607  // Compute d_i u_\alpha for nonconservative products
608  if constexpr (has_partial_derivs) {
609  partial_derivatives<partial_derivative_tags>(
610  partial_derivs, evolved_vars, mesh,
611  logical_to_inertial_inverse_jacobian);
612  }
613 
614  // Compute volume du/dt and fluxes
615 
616  // Compute volume terms that are unrelated to moving meshes
617  db::mutate<db::add_tag_prefix<::Tags::dt, variables_tag>>(
618  box,
619  [&mesh, &partial_derivs, &temporaries, &volume_fluxes](
620  const gsl::not_null<Variables<
622  dt_vars_ptr,
623  const auto&... args) noexcept {
624  // Silence compiler warnings since we genuinely don't always need all
625  // the vars but sometimes do. This warning shows up with empty parameter
626  // packs, which means the packs aren't "used" below in the
627  // ComputeVolumeTimeDerivatives::apply call.
628  (void)partial_derivs;
629  (void)volume_fluxes;
630  (void)temporaries;
631 
632  // For now just zero dt_vars. If this is a performance bottle neck we
633  // can re-evaluate in the future.
634  dt_vars_ptr->initialize(mesh.number_of_grid_points(), 0.0);
635 
636  ComputeVolumeTimeDerivatives::apply(
637  make_not_null(&get<::Tags::dt<VariablesTags>>(*dt_vars_ptr))...,
638  make_not_null(&get<FluxVariablesTags>(*volume_fluxes))...,
639  make_not_null(&get<TemporaryTags>(*temporaries))...,
640  get<PartialDerivTags>(*partial_derivs)..., args...);
641  },
642  db::get<TimeDerivativeArgumentTags>(*box)...);
643 
644  // Add volume terms for moving meshes
645  if (const auto& mesh_velocity =
647  mesh_velocity.has_value()) {
648  db::mutate<db::add_tag_prefix<::Tags::dt, variables_tag>>(
649  box,
650  [&evolved_vars, &mesh_velocity, &volume_fluxes](
651  const gsl::not_null<Variables<db::wrap_tags_in<
652  ::Tags::dt, typename variables_tag::tags_list>>*>
653  dt_vars_ptr,
655  div_mesh_velocity) noexcept {
656  tmpl::for_each<flux_variables>([&div_mesh_velocity, &dt_vars_ptr,
657  &evolved_vars, &mesh_velocity,
658  &volume_fluxes](auto tag_v) noexcept {
659  // Modify fluxes for moving mesh
660  using var_tag = typename decltype(tag_v)::type;
661  using flux_var_tag =
664  auto& flux_var = get<flux_var_tag>(*volume_fluxes);
665  // Loop over all independent components of flux_var
666  for (size_t flux_var_storage_index = 0;
667  flux_var_storage_index < flux_var.size();
668  ++flux_var_storage_index) {
669  // Get the flux variable's tensor index, e.g. (i,j) for a F^i of
670  // the spatial velocity (or some other spatial tensor).
671  const auto flux_var_tensor_index =
672  flux_var.get_tensor_index(flux_var_storage_index);
673  // Remove the first index from the flux tensor index, gets back
674  // (j)
675  const auto var_tensor_index =
676  all_but_specified_element_of(flux_var_tensor_index, 0);
677  // Set flux_index to (i)
678  const size_t flux_index = gsl::at(flux_var_tensor_index, 0);
679 
680  // We now need to index flux(i,j) -= u(j) * v_g(i)
681  flux_var[flux_var_storage_index] -=
682  get<var_tag>(evolved_vars).get(var_tensor_index) *
683  mesh_velocity->get(flux_index);
684  }
685 
686  // Modify time derivative (i.e. source terms) for moving mesh
687  auto& dt_var = get<::Tags::dt<var_tag>>(*dt_vars_ptr);
688  for (size_t dt_var_storage_index = 0;
689  dt_var_storage_index < dt_var.size(); ++dt_var_storage_index) {
690  // This is S -> S - u d_i v^i_g
691  dt_var[dt_var_storage_index] -=
692  get<var_tag>(evolved_vars)[dt_var_storage_index] *
693  get(*div_mesh_velocity);
694  }
695  });
696  },
697  db::get<::domain::Tags::DivMeshVelocity>(*box));
698 
699  // We add the mesh velocity to all equations that don't have flux terms.
700  // This doesn't need to be equal to the equations that have partial
701  // derivatives. For example, the scalar field evolution equation in
702  // first-order form does not have any partial derivatives but still needs
703  // the velocity term added. This is because the velocity term arises from
704  // transforming the time derivative.
705  using non_flux_tags = tmpl::list_difference<variables_tags, flux_variables>;
706 
707  db::mutate<db::add_tag_prefix<::Tags::dt, variables_tag>>(
708  box, [&mesh_velocity, &partial_derivs](
709  const gsl::not_null<Variables<db::wrap_tags_in<
710  ::Tags::dt, typename variables_tag::tags_list>>*>
711  dt_vars) noexcept {
712  tmpl::for_each<non_flux_tags>(
713  [&dt_vars, &mesh_velocity,
714  &partial_derivs](auto var_tag_v) noexcept {
715  using var_tag = typename decltype(var_tag_v)::type;
716  using dt_var_tag = ::Tags::dt<var_tag>;
717  using deriv_var_tag =
719 
720  const auto& deriv_var = get<deriv_var_tag>(*partial_derivs);
721  auto& dt_var = get<dt_var_tag>(*dt_vars);
722 
723  // Loop over all independent components of the derivative of the
724  // variable.
725  for (size_t deriv_var_storage_index = 0;
726  deriv_var_storage_index < deriv_var.size();
727  ++deriv_var_storage_index) {
728  // We grab the `deriv_tensor_index`, which would be e.g.
729  // `(i, a, b)`, so `(0, 2, 3)`
730  const auto deriv_var_tensor_index =
731  deriv_var.get_tensor_index(deriv_var_storage_index);
732  // Then we drop the derivative index (the first entry) to get
733  // `(a, b)` (or `(2, 3)`)
734  const auto dt_var_tensor_index =
735  all_but_specified_element_of(deriv_var_tensor_index, 0);
736  // Set `deriv_index` to `i` (or `0` in the example)
737  const size_t deriv_index = gsl::at(deriv_var_tensor_index, 0);
738  dt_var.get(dt_var_tensor_index) +=
739  mesh_velocity->get(deriv_index) *
740  deriv_var[deriv_var_storage_index];
741  }
742  });
743  });
744  }
745 
746  // Add the flux divergence term to du_\alpha/dt, which must be done
747  // after the corrections for the moving mesh are made.
748  if constexpr (has_fluxes) {
749  // Compute the divergence outside of the `db::mutate` call so that we don't
750  // have to pass in the inertial coordinates and the Jacobian when using the
751  // strong form. This is a minor and trivial "optimization".
752  Variables<tmpl::list<::Tags::div<FluxVariablesTags>...>> div_fluxes{
753  mesh.number_of_grid_points()};
754  if (dg_formulation == ::dg::Formulation::StrongInertial) {
755  divergence(make_not_null(&div_fluxes), *volume_fluxes, mesh,
756  logical_to_inertial_inverse_jacobian);
757  } else if (dg_formulation == ::dg::Formulation::WeakInertial) {
758  // We should ideally not recompute the
759  // det_jac_times_inverse_jacobian for non-moving meshes.
760  if constexpr (Dim == 1) {
761  weak_divergence(make_not_null(&div_fluxes), *volume_fluxes, mesh, {});
762  } else {
763  const auto& inertial_coords =
764  db::get<domain::Tags::Coordinates<Dim, Frame::Inertial>>(*box);
765  // The Jacobian should be computed as a compute tag
766  const auto jacobian =
769  Frame::Inertial>>(*box))
770  .second;
771  InverseJacobian<DataVector, Dim, Frame::Logical, Frame::Inertial>
772  det_jac_times_inverse_jacobian{};
774  make_not_null(&det_jac_times_inverse_jacobian), mesh,
775  inertial_coords, jacobian);
776  weak_divergence(make_not_null(&div_fluxes), *volume_fluxes, mesh,
777  det_jac_times_inverse_jacobian);
778  }
779  div_fluxes *=
780  get(db::get<
782  *box));
783  } else {
784  ERROR("Unsupported DG formulation: " << dg_formulation);
785  }
786 
787  db::mutate<db::add_tag_prefix<::Tags::dt, variables_tag>>(
788  box, [dg_formulation,
789  &div_fluxes](const gsl::not_null<Variables<db::wrap_tags_in<
790  ::Tags::dt, typename variables_tag::tags_list>>*>
791  dt_vars_ptr) noexcept {
792  tmpl::for_each<flux_variables>(
793  [&dg_formulation, &dt_vars_ptr,
794  &div_fluxes](auto var_tag_v) noexcept {
795  using var_tag = typename decltype(var_tag_v)::type;
796  auto& dt_var = get<::Tags::dt<var_tag>>(*dt_vars_ptr);
797  const auto& div_flux = get<::Tags::div<
799  div_fluxes);
800  if (dg_formulation == ::dg::Formulation::StrongInertial) {
801  for (size_t storage_index = 0; storage_index < dt_var.size();
802  ++storage_index) {
803  dt_var[storage_index] -= div_flux[storage_index];
804  }
805  } else {
806  for (size_t storage_index = 0; storage_index < dt_var.size();
807  ++storage_index) {
808  dt_var[storage_index] += div_flux[storage_index];
809  }
810  }
811  });
812  });
813  } else {
814  (void)dg_formulation;
815  }
816 }
817 
818 template <typename Metavariables>
819 template <size_t Dim, typename BoundaryCorrection, typename TemporaryTags,
820  typename DbTagsList, typename VariablesTags,
821  typename... PackageDataVolumeTags>
822 void ComputeTimeDerivative<Metavariables>::compute_internal_mortar_data(
823  const gsl::not_null<db::DataBox<DbTagsList>*> box,
824  const BoundaryCorrection& boundary_correction,
825  const Variables<VariablesTags>& volume_evolved_vars,
826  const Variables<db::wrap_tags_in<
827  ::Tags::Flux, typename Metavariables::system::flux_variables,
828  tmpl::size_t<Dim>, Frame::Inertial>>& volume_fluxes,
829  const Variables<TemporaryTags>& volume_temporaries) noexcept {
830  using system = typename Metavariables::system;
831  using variables_tags = typename system::variables_tag::tags_list;
832  using flux_variables = typename system::flux_variables;
833  using fluxes_tags = db::wrap_tags_in<::Tags::Flux, flux_variables,
834  tmpl::size_t<Metavariables::volume_dim>,
836  using temporary_tags_for_face =
837  typename BoundaryCorrection::dg_package_data_temporary_tags;
838  using primitive_tags_for_face = typename detail::get_primitive_vars<
839  system::has_primitive_and_conservative_vars>::
840  template f<BoundaryCorrection>;
841  using mortar_tags_list = typename BoundaryCorrection::dg_package_field_tags;
842  static_assert(
843  std::is_same_v<VariablesTags, variables_tags>,
844  "The evolved variables passed in should match the evolved variables of "
845  "the system.");
846 
847  const Element<Dim>& element = db::get<domain::Tags::Element<Dim>>(*box);
848  const Mesh<Dim>& volume_mesh = db::get<domain::Tags::Mesh<Dim>>(*box);
849  const auto& face_meshes =
850  get<domain::Tags::Interface<domain::Tags::InternalDirections<Dim>,
851  domain::Tags::Mesh<Dim - 1>>>(*box);
852  const auto& mortar_meshes = db::get<Tags::MortarMesh<Dim>>(*box);
853  const auto& mortar_sizes = db::get<Tags::MortarSize<Dim>>(*box);
854  const TimeStepId& temporal_id = db::get<::Tags::TimeStepId>(*box);
856  moving_mesh_map =
858  Frame::Inertial>>(
859  *box);
860  const std::unordered_map<Direction<Dim>, tnsr::i<DataVector, Dim>>&
861  unnormalized_normal_covectors = db::get<domain::Tags::Interface<
864 
867  face_mesh_velocities = db::get<domain::Tags::Interface<
870 
871  using dg_package_data_projected_tags =
872  tmpl::append<variables_tags, fluxes_tags, temporary_tags_for_face,
873  primitive_tags_for_face>;
874  Variables<tmpl::remove_duplicates<tmpl::push_back<
875  tmpl::append<dg_package_data_projected_tags,
876  detail::inverse_spatial_metric_tag<system>>,
877  detail::OneOverNormalVectorMagnitude, detail::NormalVector<Dim>>>>
878  fields_on_face{};
879  for (const auto& [direction, neighbors_in_direction] : element.neighbors()) {
880  (void)neighbors_in_direction; // unused variable
881  // In order to reduce memory allocations we handle both the upper and
882  // lower neighbors in each direction together since computing the
883  // contributions to the faces is guaranteed to require the same number of
884  // grid points.
885  if (direction.side() == Side::Upper and
886  element.neighbors().count(direction.opposite()) != 0) {
887  continue;
888  }
889 
890  const auto internal_mortars =
891  [&box, &boundary_correction, &element, &face_mesh_velocities,
892  &fields_on_face, &mortar_meshes, &mortar_sizes, &moving_mesh_map,
893  &temporal_id, &unnormalized_normal_covectors, &volume_evolved_vars,
894  &volume_fluxes, &volume_temporaries,
895  &volume_mesh](const Mesh<Dim - 1>& face_mesh,
896  const Direction<Dim>& local_direction) noexcept {
897  // We may not need to bring the volume fluxes or temporaries to the
898  // boundary since that depends on the specific boundary correction we
899  // are using. Silence compilers warnings about them being unused.
900  (void)volume_fluxes;
901  (void)volume_temporaries;
902  // We only need the box for volume tags and primitive tags. Silence
903  // unused variable compiler warnings in the case where we don't have
904  // both of those.
905  (void)box;
906 
907  // This helper does the following:
908  //
909  // 1. Use a helper function to get data onto the faces. Done either by
910  // slicing (Gauss-Lobatto points) or interpolation (Gauss points).
911  // This is done using the `project_contiguous_data_to_boundary` and
912  // `project_tensors_to_boundary` functions.
913  //
914  // 2. Invoke the boundary correction to get the packaged data. Note
915  // that this is done on the *face* and NOT the mortar.
916  //
917  // 3. Project the packaged data onto the DG mortars (these might need
918  // re-projection onto subcell mortars later).
919 
920  // Perform step 1
922  volume_evolved_vars, volume_mesh,
923  local_direction);
924  if constexpr (tmpl::size<fluxes_tags>::value != 0) {
926  volume_fluxes, volume_mesh,
927  local_direction);
928  }
929  if constexpr (tmpl::size<tmpl::append<
930  temporary_tags_for_face,
931  detail::inverse_spatial_metric_tag<system>>>::
932  value != 0) {
934  tmpl::append<temporary_tags_for_face,
935  detail::inverse_spatial_metric_tag<system>>>(
936  make_not_null(&fields_on_face), volume_temporaries, volume_mesh,
937  local_direction);
938  }
939  if constexpr (system::has_primitive_and_conservative_vars and
940  tmpl::size<primitive_tags_for_face>::value != 0) {
941  project_tensors_to_boundary<primitive_tags_for_face>(
942  make_not_null(&fields_on_face),
943  db::get<typename system::primitive_variables_tag>(*box),
944  volume_mesh, local_direction);
945  }
946 
947  // Normalize the normal vectors. In the case of a time-independent
948  // mesh with no inverse spatial metric tag (i.e. flat background) we
949  // could cache the normal vectors as simple tags in the DataBox as an
950  // optimization. This will require making the
951  // Normalized<UnnormalizedFaceNormal> a simple tag instead of a
952  // compute tag.
953  db::mutate<evolution::dg::Tags::InternalFace::
954  NormalCovectorAndMagnitude<Dim>>(
955  box, [&fields_on_face, &local_direction, &moving_mesh_map,
956  &unnormalized_normal_covectors](
957  const auto normal_covector_and_magnitude_ptr) noexcept {
958  detail::unit_normal_vector_and_covector_and_magnitude<system>(
959  normal_covector_and_magnitude_ptr,
960  make_not_null(&fields_on_face), local_direction,
961  unnormalized_normal_covectors, moving_mesh_map);
962  });
963 
964  // Perform step 2
966  NormalCovectorAndMagnitude<Dim>>(*box)
967  .at(local_direction)
968  .has_value(),
969  "The magnitude of the normal vector and the unit normal "
970  "covector have not been computed, even though they should "
971  "have been. Direction: "
972  << local_direction);
973 
974  Variables<mortar_tags_list> packaged_data{
975  face_mesh.number_of_grid_points()};
976  // The DataBox is passed in for retrieving the `volume_tags`
977  const double max_abs_char_speed_on_face =
978  detail::dg_package_data<system>(
979  make_not_null(&packaged_data), boundary_correction,
980  fields_on_face,
983  NormalCovectorAndMagnitude<Dim>>(*box)
984  .at(local_direction)),
985  face_mesh_velocities.at(local_direction), *box,
986  typename BoundaryCorrection::dg_package_data_volume_tags{},
987  dg_package_data_projected_tags{});
988  (void)max_abs_char_speed_on_face;
989 
990  // Perform step 3
991  const auto& neighbors_in_local_direction =
992  element.neighbors().at(local_direction);
993  for (const auto& neighbor : neighbors_in_local_direction) {
994  const auto mortar_id = std::make_pair(local_direction, neighbor);
995  const auto& mortar_mesh = mortar_meshes.at(mortar_id);
996  const auto& mortar_size = mortar_sizes.at(mortar_id);
997 
998  // Project the data from the face to the mortar.
999  // Where no projection is necessary we `std::move` the data
1000  // directly to avoid a copy. We can't move the data or modify it
1001  // in-place when projecting, because in that case the face may
1002  // touch two mortars so we need to keep the data around.
1003  auto boundary_data_on_mortar =
1005  // NOLINTNEXTLINE(bugprone-use-after-move)
1006  ? ::dg::project_to_mortar(packaged_data, face_mesh,
1008  : std::move(packaged_data);
1009 
1010  // Store the boundary data on this side of the mortar in a way
1011  // that is agnostic to the type of boundary correction used. This
1012  // currently requires an additional allocation that could be
1013  // eliminated either by:
1014  //
1015  // 1. Having non-owning Variables
1016  //
1017  // 2. Allow stealing the allocation out of a Variables (and
1018  // inserting an allocation).
1019  std::vector<double> type_erased_boundary_data_on_mortar{
1020  boundary_data_on_mortar.data(),
1021  boundary_data_on_mortar.data() +
1022  boundary_data_on_mortar.size()};
1023  db::mutate<Tags::MortarData<Dim>>(
1024  box, [&face_mesh, &mortar_id, &temporal_id,
1025  &type_erased_boundary_data_on_mortar](
1026  const auto mortar_data_ptr) noexcept {
1027  mortar_data_ptr->at(mortar_id).insert_local_mortar_data(
1028  temporal_id, face_mesh,
1029  std::move(type_erased_boundary_data_on_mortar));
1030  });
1031  }
1032  };
1033 
1034  if (fields_on_face.number_of_grid_points() !=
1035  face_meshes.at(direction).number_of_grid_points()) {
1036  fields_on_face.initialize(
1037  face_meshes.at(direction).number_of_grid_points());
1038  }
1039  internal_mortars(face_meshes.at(direction), direction);
1040 
1041  if (element.neighbors().count(direction.opposite()) != 0) {
1042  if (fields_on_face.number_of_grid_points() !=
1043  face_meshes.at(direction.opposite()).number_of_grid_points()) {
1044  fields_on_face.initialize(
1045  face_meshes.at(direction.opposite()).number_of_grid_points());
1046  }
1047  internal_mortars(face_meshes.at(direction.opposite()),
1048  direction.opposite());
1049  }
1050  }
1051 }
1052 
1053 template <typename Metavariables>
1054 template <typename DbTagsList>
1055 void ComputeTimeDerivative<Metavariables>::
1056  boundary_terms_nonconservative_products(
1057  const gsl::not_null<db::DataBox<DbTagsList>*> box) noexcept {
1058  using system = typename Metavariables::system;
1059  using variables_tag = typename system::variables_tag;
1060  using DirectionsTag =
1062 
1063  using interface_normal_dot_fluxes_tag = domain::Tags::Interface<
1065 
1066  using interface_compute_item_argument_tags = tmpl::transform<
1067  typename Metavariables::system::normal_dot_fluxes::argument_tags,
1068  tmpl::bind<domain::Tags::Interface, DirectionsTag, tmpl::_1>>;
1069 
1071  tmpl::list<interface_normal_dot_fluxes_tag>,
1072  tmpl::push_front<interface_compute_item_argument_tags, DirectionsTag>>(
1073  [](const auto boundary_fluxes_ptr,
1075  internal_directions,
1076  const auto&... tensors) noexcept {
1077  for (const auto& direction : internal_directions) {
1078  apply_flux(make_not_null(&boundary_fluxes_ptr->at(direction)),
1079  tensors.at(direction)...);
1080  }
1081  },
1082  box);
1083 }
1084 
1085 template <typename Metavariables>
1086 template <size_t VolumeDim, typename BoundaryScheme, typename DbTagsList>
1087 void ComputeTimeDerivative<Metavariables>::
1088  fill_mortar_data_for_internal_boundaries(
1089  const gsl::not_null<db::DataBox<DbTagsList>*> box) noexcept {
1090  using temporal_id_tag = typename BoundaryScheme::temporal_id_tag;
1091  using all_mortar_data_tag =
1093  using boundary_data_computer =
1094  typename BoundaryScheme::boundary_data_computer;
1095 
1096  // Collect data on element interfaces
1097  auto boundary_data_on_interfaces =
1098  interface_apply<domain::Tags::InternalDirections<VolumeDim>,
1099  boundary_data_computer>(*box);
1100 
1101  // Project collected data to all internal mortars and store in DataBox
1102  const auto& element = db::get<domain::Tags::Element<VolumeDim>>(*box);
1103  const auto& face_meshes =
1104  get<domain::Tags::Interface<domain::Tags::InternalDirections<VolumeDim>,
1105  domain::Tags::Mesh<VolumeDim - 1>>>(*box);
1106  const auto& mortar_meshes =
1107  get<::Tags::Mortars<domain::Tags::Mesh<VolumeDim - 1>, VolumeDim>>(*box);
1108  const auto& mortar_sizes =
1109  get<::Tags::Mortars<::Tags::MortarSize<VolumeDim - 1>, VolumeDim>>(*box);
1110  const auto& temporal_id = get<temporal_id_tag>(*box);
1111  for (const auto& [direction, neighbors] : element.neighbors()) {
1112  const auto& face_mesh = face_meshes.at(direction);
1113  for (const auto& neighbor : neighbors) {
1114  const auto mortar_id = std::make_pair(direction, neighbor);
1115  const auto& mortar_mesh = mortar_meshes.at(mortar_id);
1116  const auto& mortar_size = mortar_sizes.at(mortar_id);
1117 
1118  // Project the data from the face to the mortar.
1119  // Where no projection is necessary we `std::move` the data directly to
1120  // avoid a copy. We can't move the data or modify it in-place when
1121  // projecting, because in that case the face may touch more than one
1122  // mortar so we need to keep the data around.
1123  auto boundary_data_on_mortar =
1125  ? boundary_data_on_interfaces.at(direction).project_to_mortar(
1126  face_mesh, mortar_mesh, mortar_size)
1127  : std::move(boundary_data_on_interfaces.at(direction));
1128 
1129  // Store the boundary data on this side of the mortar
1130  db::mutate<all_mortar_data_tag>(
1131  box, [&mortar_id, &temporal_id, &boundary_data_on_mortar ](
1133  all_mortar_data) noexcept {
1134  all_mortar_data->at(mortar_id).local_insert(
1135  temporal_id, std::move(boundary_data_on_mortar));
1136  });
1137  }
1138  }
1139 }
1140 
1141 template <typename Metavariables>
1142 template <typename ParallelComponent, typename DbTagsList>
1143 void ComputeTimeDerivative<Metavariables>::send_data_for_fluxes(
1145  const db::DataBox<DbTagsList>& box) noexcept {
1146  using system = typename Metavariables::system;
1147  constexpr size_t volume_dim = Metavariables::volume_dim;
1148  auto& receiver_proxy =
1149  Parallel::get_parallel_component<ParallelComponent>(*cache);
1150  const auto& element = db::get<domain::Tags::Element<volume_dim>>(box);
1151 
1152  if constexpr (detail::has_boundary_correction_v<system>) {
1153  const auto& time_step_id = db::get<::Tags::TimeStepId>(box);
1154  const auto& all_mortar_data =
1155  db::get<evolution::dg::Tags::MortarData<volume_dim>>(box);
1156 
1157  for (const auto& [direction, neighbors] : element.neighbors()) {
1158  const auto& orientation = neighbors.orientation();
1159  const auto direction_from_neighbor = orientation(direction.opposite());
1160 
1161  for (const auto& neighbor : neighbors) {
1162  const std::pair mortar_id{direction, neighbor};
1163 
1164  std::pair<Mesh<volume_dim - 1>, std::vector<double>>
1165  neighbor_boundary_data_on_mortar =
1166  *all_mortar_data.at(mortar_id).local_mortar_data();
1167  ASSERT(time_step_id == all_mortar_data.at(mortar_id).time_step_id(),
1168  "The current time step id of the volume is "
1169  << time_step_id
1170  << "but the time step id on the mortar with mortar id "
1171  << mortar_id << " is "
1172  << all_mortar_data.at(mortar_id).time_step_id());
1173 
1174  // Reorient the data to the neighbor orientation if necessary
1175  if (not orientation.is_aligned()) {
1176  ERROR("Currently don't support unaligned meshes");
1177  // neighbor_boundary_data_on_mortar.orient_on_slice(
1178  // mortar_meshes.at(mortar_id).extents(), direction.dimension(),
1179  // orientation);
1180  }
1181 
1184  data{neighbor_boundary_data_on_mortar.first,
1185  {},
1186  {std::move(neighbor_boundary_data_on_mortar.second)},
1187  db::get<tmpl::conditional_t<Metavariables::local_time_stepping,
1189  ::Tags::TimeStepId>>(box)};
1190 
1191  // Send mortar data (the `std::tuple` named `data`) to neighbor
1194  volume_dim>>(
1195  receiver_proxy[neighbor], time_step_id,
1196  std::make_pair(std::pair{direction_from_neighbor, element.id()},
1197  data));
1198  }
1199  }
1200  } else {
1201  using BoundaryScheme = typename Metavariables::boundary_scheme;
1202  using temporal_id_tag = typename BoundaryScheme::temporal_id_tag;
1203  using receive_temporal_id_tag =
1204  typename BoundaryScheme::receive_temporal_id_tag;
1205  using fluxes_inbox_tag = ::dg::FluxesInboxTag<BoundaryScheme>;
1206  using all_mortar_data_tag =
1208 
1209  const auto& all_mortar_data = get<all_mortar_data_tag>(box);
1210  const auto& temporal_id = db::get<temporal_id_tag>(box);
1211  const auto& receive_temporal_id = db::get<receive_temporal_id_tag>(box);
1212  const auto& mortar_meshes = db::get<
1213  ::Tags::Mortars<domain::Tags::Mesh<volume_dim - 1>, volume_dim>>(box);
1214 
1215  // Iterate over neighbors
1216  for (const auto& direction_and_neighbors : element.neighbors()) {
1217  const auto& direction = direction_and_neighbors.first;
1218  const size_t dimension = direction.dimension();
1219  const auto& neighbors_in_direction = direction_and_neighbors.second;
1220  const auto& orientation = neighbors_in_direction.orientation();
1221  const auto direction_from_neighbor = orientation(direction.opposite());
1222 
1223  for (const auto& neighbor : neighbors_in_direction) {
1224  const ::dg::MortarId<volume_dim> mortar_id{direction, neighbor};
1225 
1226  // Make a copy of the local boundary data on the mortar to send to the
1227  // neighbor
1228  ASSERT(all_mortar_data.find(mortar_id) != all_mortar_data.end(),
1229  "Mortar data on mortar "
1230  << mortar_id
1231  << " not available for sending. Did you forget to collect "
1232  "the data on mortars?");
1233  auto neighbor_boundary_data_on_mortar =
1234  all_mortar_data.at(mortar_id).local_data(temporal_id);
1235 
1236  // Reorient the data to the neighbor orientation if necessary
1237  if (not orientation.is_aligned()) {
1238  neighbor_boundary_data_on_mortar.orient_on_slice(
1239  mortar_meshes.at(mortar_id).extents(), dimension, orientation);
1240  }
1241 
1242  // Send mortar data (named 'neighbor_boundary_data_on_mortar') to
1243  // neighbor
1244  Parallel::receive_data<fluxes_inbox_tag>(
1245  receiver_proxy[neighbor], temporal_id,
1246  std::make_pair(
1247  ::dg::MortarId<volume_dim>{direction_from_neighbor,
1248  element.id()},
1249  std::make_pair(receive_temporal_id,
1250  std::move(neighbor_boundary_data_on_mortar))));
1251  }
1252  }
1253  }
1254 }
1255 } // namespace evolution::dg::Actions
domain::Tags::UnnormalizedFaceNormal
Definition: FaceNormal.hpp:112
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
all_but_specified_element_of
std::array< T, Dim - 1 > all_but_specified_element_of(const std::array< T, Dim > &a, const size_t element_to_remove) noexcept
Construct an array from an existing array omitting one element.
Definition: StdArrayHelpers.hpp:98
db::mutate_apply
constexpr void 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:1417
utility
Frame::Inertial
Definition: IndexType.hpp:44
domain::CoordinateMapBase
Abstract base class for CoordinateMap.
Definition: CoordinateMap.hpp:45
Divergence.hpp
Tags::TimeStepId
Tag for TimeStepId for the algorithm state.
Definition: Tags.hpp:33
get
constexpr Tag::type & get(Variables< TagList > &v) noexcept
Return Tag::type pointing into the contiguous array.
Definition: Variables.hpp:638
Parallel::GlobalCache
Definition: ElementReceiveInterpPoints.hpp:15
dg::project_to_mortar
Variables< Tags > project_to_mortar(const Variables< Tags > &vars, const Mesh< Dim > &face_mesh, const Mesh< Dim > &mortar_mesh, const MortarSize< Dim > &mortar_size) noexcept
Definition: MortarHelpers.hpp:85
Frame::Grid
Definition: IndexType.hpp:43
unordered_set
std::pair
GlobalCache.hpp
Tags::MortarSize
Definition: Tags.hpp:48
Tags.hpp
vector
std::system
T system(T... args)
db::add_tag_prefix
typename detail::add_tag_prefix_impl< Prefix, Tag, Args... >::type add_tag_prefix
Definition: PrefixHelpers.hpp:51
evolution::dg::Tags::InternalFace::NormalCovector
The normal covector to the interface.
Definition: NormalVectorTags.hpp:25
Tags::Variables
Definition: VariablesTag.hpp:21
domain::Tags::MeshVelocity
The mesh velocity.
Definition: TagsTimeDependent.hpp:138
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:969
Tags::Next
Prefix indicating the value a quantity will take on the next iteration of the algorithm.
Definition: Prefixes.hpp:118
Tags::Flux
Prefix indicating a flux.
Definition: Prefixes.hpp:40
Direction
Definition: Direction.hpp:23
db::mutate
void mutate(const gsl::not_null< DataBox< TagList > * > box, Invokable &&invokable, Args &&... args) noexcept
Allows changing the state of one or more non-computed elements in the DataBox.
Definition: DataBox.hpp:859
divergence
Scalar< DataVector > divergence(const tnsr::I< DataVector, Dim, DerivativeFrame > &input, const Mesh< Dim > &mesh, const InverseJacobian< DataVector, Dim, Frame::Logical, DerivativeFrame > &inverse_jacobian) noexcept
Compute the divergence of the vector input
Definition: Divergence.cpp:26
Element
Definition: Element.hpp:29
dot_product
void dot_product(const gsl::not_null< Scalar< DataType > * > dot_product, const Tensor< DataType, Symmetry< 1 >, index_list< Index >> &vector_a, const Tensor< DataType, Symmetry< 1 >, index_list< Index >> &vector_b) noexcept
Compute the Euclidean dot product of two vectors or one forms.
Definition: DotProduct.hpp:24
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.
evolution::Tags::BoundaryCorrection
The boundary correction used for coupling together neighboring cells or applying boundary conditions.
Definition: BoundaryCorrectionTags.hpp:37
ERROR
#define ERROR(m)
prints an error message to the standard error stream and aborts the program.
Definition: Error.hpp:36
CREATE_HAS_TYPE_ALIAS
#define CREATE_HAS_TYPE_ALIAS(ALIAS_NAME)
Generate a type trait to check if a class has a type alias with a particular name,...
Definition: CreateHasTypeAlias.hpp:27
determinant_and_inverse
void determinant_and_inverse(const gsl::not_null< Scalar< T > * > det, const gsl::not_null< Tensor< T, Symm, tmpl::list< change_index_up_lo< Index1 >, change_index_up_lo< Index0 >>> * > inv, const Tensor< T, Symm, tmpl::list< Index0, Index1 >> &tensor) noexcept
Computes the determinant and inverse of a rank-2 Tensor.
Definition: DeterminantAndInverse.hpp:371
Tags::div
Prefix indicating the divergence.
Definition: Divergence.hpp:44
weak_divergence
void weak_divergence(const gsl::not_null< Variables< tmpl::list< Tags::div< FluxTags >... >> * > divergence_of_fluxes, const Variables< tmpl::list< FluxTags... >> &fluxes, const Mesh< Dim > &mesh, const InverseJacobian< DataVector, Dim, Frame::Logical, Frame::Inertial > &det_jac_times_inverse_jacobian) noexcept
Compute the weak form divergence of fluxes.
Definition: WeakDivergence.hpp:65
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:353
DataBox.hpp
domain::Tags::InverseJacobian< Dim, Frame::Logical, Frame::Inertial >
evolution::dg::Tags::InternalFace
Tags used on the interior faces of the elements
Definition: NormalVectorTags.hpp:17
evolution::dg::Actions
Actions for using the discontinuous Galerkin to evolve hyperbolic partial differential equations.
Definition: ApplyBoundaryCorrections.hpp:41
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:156
tuples::TaggedTuple
An associative container that is indexed by structs.
Definition: TaggedTuple.hpp:271
dg::metric_identity_det_jac_times_inv_jac
void metric_identity_det_jac_times_inv_jac(gsl::not_null< InverseJacobian< DataVector, Dim, Frame::Logical, Frame::Inertial > * > det_jac_times_inverse_jacobian, const Mesh< Dim > &mesh, const tnsr::I< DataVector, Dim, Frame::Inertial > &inertial_coords, const Jacobian< DataVector, Dim, Frame::Logical, Frame::Inertial > &jacobian) noexcept
Compute the Jacobian determinant times the inverse Jacobian so that the result is divergence-free.
TimeStepId
Definition: TimeStepId.hpp:25
domain::CoordinateMaps::Tags::CoordinateMap
Definition: Tags.hpp:30
std::decay_t
ASSERT
#define ASSERT(a, m)
Assert that an expression should be true.
Definition: Assert.hpp:49
evolution::dg::project_tensors_to_boundary
void project_tensors_to_boundary(const gsl::not_null< Variables< FaceVarsTagsList > * > face_fields, const Variables< VolumeVarsTagsList > &volume_fields, const Mesh< Dim > &volume_mesh, const Direction< Dim > &direction) noexcept
Projects a subset of the tensors in the volume_fields onto the face.
Definition: ProjectToBoundary.hpp:162
Variables.hpp
Mesh
Holds the number of grid points, basis, and quadrature in each direction of the computational grid.
Definition: Mesh.hpp:48
Tags::dt
Prefix indicating a time derivative.
Definition: Prefixes.hpp:29
std::is_final
dg::needs_projection
bool needs_projection(const Mesh< Dim > &face_mesh, const Mesh< Dim > &mortar_mesh, const MortarSize< Dim > &mortar_size) noexcept
Definition: MortarHelpers.hpp:74
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:49
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:365
evolution::dg::project_contiguous_data_to_boundary
void project_contiguous_data_to_boundary(const gsl::not_null< Variables< FaceVarsTagsList > * > face_fields, const Variables< VolumeVarsTagsList > &volume_fields, const Mesh< Dim > &volume_mesh, const Direction< Dim > &direction) noexcept
Projects a Variables of volume data to a contiguous subset of a boundary Variables
Definition: ProjectToBoundary.hpp:53
Tensor.hpp
domain::Tags::InternalDirections
Definition: Tags.hpp:269
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
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
make_not_null
gsl::not_null< T * > make_not_null(T *ptr) noexcept
Construct a not_null from a pointer. Often this will be done as an implicit conversion,...
Definition: Gsl.hpp:880
Direction.hpp
Mesh.hpp
PartialDerivatives.hpp
domain::Tags::DetInvJacobian
The determinant of the inverse Jacobian from the source frame to the target frame.
Definition: Tags.hpp:234
dg::Formulation
Formulation
The DG formulation to use.
Definition: Formulation.hpp:34
dg::mortar_mesh
Mesh< Dim > mortar_mesh(const Mesh< Dim > &face_mesh1, const Mesh< Dim > &face_mesh2) noexcept
Definition: MortarHelpers.cpp:19
Prefixes.hpp
std::unordered_map
std::data
T data(T... args)
type_traits
dg::mortar_size
std::array< Spectral::MortarSize, Dim - 1 > mortar_size(const ElementId< Dim > &self, const ElementId< Dim > &neighbor, const size_t dimension, const OrientationMap< Dim > &orientation) noexcept
Definition: MortarHelpers.cpp:45
TMPL.hpp
Element::neighbors
const Neighbors_t & neighbors() const noexcept
Information about the neighboring Elements.
Definition: Element.hpp:66
gsl::not_null
Require a pointer to not be a nullptr
Definition: ReadSpecThirdOrderPiecewisePolynomial.hpp:13
Tags::Mortars
Definition: Tags.hpp:37