ComputeTimeDerivative.hpp
1 // Distributed under the MIT License.
2 // See LICENSE.txt for details.
3 
4 #pragma once
5 
6 #include <tuple>
7 #include <type_traits>
8 #include <unordered_set>
9 #include <utility>
10 
15 #include "DataStructures/VariablesTag.hpp"
16 #include "Domain/InterfaceHelpers.hpp"
17 #include "Domain/Tags.hpp"
18 #include "Domain/TagsTimeDependent.hpp"
19 #include "NumericalAlgorithms/DiscontinuousGalerkin/Formulation.hpp"
20 #include "NumericalAlgorithms/DiscontinuousGalerkin/MortarHelpers.hpp"
21 #include "NumericalAlgorithms/DiscontinuousGalerkin/Tags.hpp"
25 #include "NumericalAlgorithms/Spectral/Projection.hpp"
26 #include "Utilities/Gsl.hpp"
27 #include "Utilities/TMPL.hpp"
28 
29 /// \cond
30 namespace tuples {
31 template <typename...>
32 class TaggedTuple;
33 } // namespace tuples
34 
35 namespace Parallel {
36 template <typename Metavariables>
37 class GlobalCache;
38 } // namespace Parallel
39 /// \endcond
40 
42 /*!
43  * \brief Computes the time derivative for a DG time step.
44  *
45  * Computes the volume fluxes, the divergence of the fluxes and all additional
46  * interior contributions to the time derivatives (both nonconservative products
47  * and source terms). The internal mortar data is also computed.
48  *
49  * The general first-order hyperbolic evolution equation solved for conservative
50  * systems is:
51  *
52  * \f{align*}{
53  * \frac{\partial u_\alpha}{\partial \hat{t}}
54  * + \partial_{i}
55  * \left(F^i_\alpha - v^i_g u_\alpha\right)
56  * = S_\alpha-u_\alpha\partial_i v^i_g,
57  * \f}
58  *
59  * where \f$F^i_{\alpha}\f$ are the fluxes when the mesh isn't moving,
60  * \f$v^i_g\f$ is the velocity of the mesh, \f$u_{\alpha}\f$ are the evolved
61  * variables, \f$S_{\alpha}\f$ are the source terms, and \f$\hat{t}\f$ is the
62  * time in the logical frame. For evolution equations that do not have any
63  * fluxes and only nonconservative products we evolve:
64  *
65  * \f{align*}{
66  * \frac{\partial u_\alpha}{\partial \hat{t}}
67  * +\left(B^i_{\alpha\beta}-v^i_g \delta_{\alpha\beta}
68  * \right)\partial_{i}u_\beta = S_\alpha.
69  * \f}
70  *
71  * Finally, for equations with both conservative terms and nonconservative
72  * products we use:
73  *
74  * \f{align*}{
75  * \frac{\partial u_\alpha}{\partial \hat{t}}
76  * + \partial_{i}
77  * \left(F^i_\alpha - v^i_g u_\alpha\right)
78  * +B^i_{\alpha\beta}\partial_{i}u_\beta
79  * = S_\alpha-u_\alpha\partial_i v^i_g,
80  * \f}
81  *
82  * where \f$B^i_{\alpha\beta}\f$ is the matrix for the nonconservative products.
83  *
84  * The mesh velocity is added to the flux automatically if the mesh is moving.
85  * That is,
86  *
87  * \f{align*}{
88  * F^i_{\alpha}\to F^i_{\alpha}-v^i_{g} u_{\alpha}
89  * \f}
90  *
91  * The source terms are also altered automatically by adding:
92  *
93  * \f{align*}{
94  * -u_\alpha \partial_i v^i_g,
95  * \f}
96  *
97  * For systems with equations that only contain nonconservative products, the
98  * following mesh velocity is automatically added to the time derivative:
99  *
100  * \f{align*}{
101  * v^i_g \partial_i u_\alpha,
102  * \f}
103  *
104  * Uses:
105  * - System:
106  * - `variables_tag`
107  * - `flux_variables`
108  * - `gradient_variables`
109  * - `compute_volume_time_derivative_terms`
110  *
111  * - DataBox:
112  * - Items in `system::compute_volume_time_derivative_terms::argument_tags`
113  * - `domain::Tags::MeshVelocity<Metavariables::volume_dim>`
114  * - `Metavariables::system::variables_tag`
115  * - `Metavariables::system::flux_variables`
116  * - `Metavariables::system::gradient_variables`
117  * - `domain::Tags::DivMeshVelocity`
118  * - `DirectionsTag`,
119  * - Required interface items for `Metavariables::system::normal_dot_fluxes`
120  *
121  * DataBox changes:
122  * - Adds: nothing
123  * - Removes: nothing
124  * - Modifies:
125  * - db::add_tag_prefix<Tags::Flux, variables_tag,
126  * tmpl::size_t<system::volume_dim>, Frame::Inertial>
127  * - `Tags::dt<system::variable_tags>`
128  * - Tags::Interface<
129  * DirectionsTag, db::add_tag_prefix<Tags::NormalDotFlux, variables_tag>>
130  * - `Tags::Mortars<typename BoundaryScheme::mortar_data_tag, VolumeDim>`
131  */
132 template <typename Metavariables>
134  public:
135  template <typename DbTagsList, typename... InboxTags, typename ArrayIndex,
136  typename ActionList, typename ParallelComponent>
140  const Parallel::GlobalCache<Metavariables>& /*cache*/,
141  const ArrayIndex& /*array_index*/, ActionList /*meta*/,
142  const ParallelComponent* /*meta*/) noexcept; // NOLINT const
143 
144  private:
145  /*
146  * Computes the volume terms for a discontinuous Galerkin scheme.
147  *
148  * The function does the following (in order):
149  *
150  * 1. Compute the partial derivatives of the `System::gradient_variables`.
151  *
152  * The partial derivatives are needed in the nonconservative product terms
153  * of the evolution equations. Any variable whose evolution equation does
154  * not contain a flux must contain a nonconservative product and the
155  * variable must be listed in the `System::gradient_variables` type alias.
156  * The partial derivatives are also needed for adding the moving mesh terms
157  * to the equations that do not have a flux term.
158  *
159  * 2. The volume time derivatives are calculated from
160  * `System::compute_volume_time_derivative_terms`
161  *
162  * The source terms and nonconservative products are contributed directly
163  * to the `dt_vars` arguments passed to the time derivative function, while
164  * the volume fluxes are computed into the `volume_fluxes` arguments. The
165  * divergence of the volume fluxes will be computed and added to the time
166  * derivatives later in the function.
167  *
168  * 3. If the mesh is moving the appropriate mesh velocity terms are added to
169  * the equations.
170  *
171  * For equations with fluxes this means that \f$-v^i_g u_\alpha\f$ is
172  * added to the fluxes and \f$-u_\alpha \partial_i v^i_g\f$ is added
173  * to the time derivatives. For equations without fluxes
174  * \f$v^i\partial_i u_\alpha\f$ is added to the time derivatives.
175  *
176  * 4. Compute flux divergence contribution and add it to the time derivatives.
177  *
178  * Either the weak or strong form can be used. Currently only the strong
179  * form is coded, but adding the weak form is quite easy.
180  *
181  * Note that the computation of the flux divergence and adding that to the
182  * time derivative must be done *after* the mesh velocity is subtracted
183  * from the fluxes.
184  */
185  template <size_t Dim, typename ComputeVolumeTimeDerivatives,
186  typename DbTagsList, typename... VariablesTags,
187  typename... TimeDerivativeArgumentTags,
188  typename... PartialDerivTags, typename... FluxVariablesTags,
189  typename... TemporaryTags>
190  static void volume_terms(
192  gsl::not_null<Variables<tmpl::list<FluxVariablesTags...>>*> volume_fluxes,
193  gsl::not_null<Variables<tmpl::list<PartialDerivTags...>>*> partial_derivs,
194  gsl::not_null<Variables<tmpl::list<TemporaryTags...>>*> temporaries,
195  ::dg::Formulation dg_formulation,
196  tmpl::list<VariablesTags...> /*variables_tags*/,
197  tmpl::list<TimeDerivativeArgumentTags...> /*meta*/) noexcept;
198 
199  template <typename... NormalDotFluxTags, typename... Args>
200  static void apply_flux(
201  gsl::not_null<Variables<tmpl::list<NormalDotFluxTags...>>*> boundary_flux,
202  const Args&... boundary_variables) noexcept {
204  make_not_null(&get<NormalDotFluxTags>(*boundary_flux))...,
205  boundary_variables...);
206  }
207 
208  template <typename DbTagsList>
209  static void boundary_terms_nonconservative_products(
210  gsl::not_null<db::DataBox<DbTagsList>*> box) noexcept;
211 
212  template <size_t VolumeDim, typename BoundaryScheme, typename DbTagsList>
213  static void fill_mortar_data_for_internal_boundaries(
214  gsl::not_null<db::DataBox<DbTagsList>*> box) noexcept;
215 };
216 
217 template <typename Metavariables>
218 template <typename DbTagsList, typename... InboxTags, typename ArrayIndex,
219  typename ActionList, typename ParallelComponent>
224  const Parallel::GlobalCache<Metavariables>& /*cache*/,
225  const ArrayIndex& /*array_index*/, ActionList /*meta*/,
226  const ParallelComponent* const /*meta*/) noexcept { // NOLINT const
227  static constexpr size_t volume_dim = Metavariables::volume_dim;
228  using system = typename Metavariables::system;
229  using variables_tags = typename system::variables_tag::tags_list;
230  using partial_derivative_tags = typename system::gradient_variables;
231  using flux_variables = typename system::flux_variables;
232  using compute_volume_time_derivative_terms =
233  typename system::compute_volume_time_derivative_terms;
234 
235  const Mesh<volume_dim>& mesh = db::get<::domain::Tags::Mesh<volume_dim>>(box);
236  // Allocate the Variables classes needed for the time derivative
237  // computation.
238  //
239  // This is factored out so that we will be able to do ADER-DG/CG where a
240  // spacetime polynomial is constructed by solving implicit equations in time
241  // using a Picard iteration. A high-order initial guess is needed to
242  // efficiently construct the ADER spacetime solution. This initial guess is
243  // obtained using continuous RK methods, and so we will want to reuse
244  // buffers. Thus, the volume_terms function returns by reference rather than
245  // by value.
246  Variables<typename compute_volume_time_derivative_terms::temporary_tags>
247  temporaries{mesh.number_of_grid_points()};
248  Variables<db::wrap_tags_in<::Tags::Flux, flux_variables,
249  tmpl::size_t<volume_dim>, Frame::Inertial>>
250  volume_fluxes{mesh.number_of_grid_points()};
251  Variables<db::wrap_tags_in<::Tags::deriv, partial_derivative_tags,
252  tmpl::size_t<volume_dim>, Frame::Inertial>>
253  partial_derivs{mesh.number_of_grid_points()};
254 
255  volume_terms<volume_dim, compute_volume_time_derivative_terms>(
256  make_not_null(&box), make_not_null(&volume_fluxes),
257  make_not_null(&partial_derivs), make_not_null(&temporaries),
258  Metavariables::dg_formulation, variables_tags{},
259  typename compute_volume_time_derivative_terms::argument_tags{});
260 
261  // The below if-else and fill_mortar_data_for_internal_boundaries are for
262  // compatibility with the current boundary schemes. Once the current
263  // boundary schemes are removed the code will be refactored to reflect that
264  // change.
265  if constexpr (not std::is_same_v<tmpl::list<>, flux_variables>) {
266  using flux_variables_tag = ::Tags::Variables<flux_variables>;
267  using fluxes_tag =
268  db::add_tag_prefix<::Tags::Flux, flux_variables_tag,
269  tmpl::size_t<Metavariables::volume_dim>,
271 
272  // We currently set the fluxes in the DataBox to interface with the
273  // boundary correction communication code.
274  db::mutate<fluxes_tag>(make_not_null(&box),
275  [&volume_fluxes](const auto fluxes_ptr) noexcept {
276  *fluxes_ptr = volume_fluxes;
277  });
278  } else {
279  boundary_terms_nonconservative_products(make_not_null(&box));
280  }
281 
282  // Compute internal boundary quantities
283  fill_mortar_data_for_internal_boundaries<
284  volume_dim, typename Metavariables::boundary_scheme>(make_not_null(&box));
285  return {std::move(box)};
286 }
287 
288 template <typename Metavariables>
289 template <size_t Dim, typename ComputeVolumeTimeDerivatives,
290  typename DbTagsList, typename... VariablesTags,
291  typename... TimeDerivativeArgumentTags, typename... PartialDerivTags,
292  typename... FluxVariablesTags, typename... TemporaryTags>
293 void ComputeTimeDerivative<Metavariables>::volume_terms(
295  const gsl::not_null<Variables<tmpl::list<FluxVariablesTags...>>*>
296  volume_fluxes,
297  const gsl::not_null<Variables<tmpl::list<PartialDerivTags...>>*>
298  partial_derivs,
299  const gsl::not_null<Variables<tmpl::list<TemporaryTags...>>*> temporaries,
300  const ::dg::Formulation dg_formulation,
301 
302  tmpl::list<VariablesTags...> /*variables_tags*/,
303  tmpl::list<TimeDerivativeArgumentTags...> /*meta*/) noexcept {
304  static constexpr bool has_partial_derivs = sizeof...(PartialDerivTags) != 0;
305  static constexpr bool has_fluxes = sizeof...(FluxVariablesTags) != 0;
306  static_assert(
307  has_fluxes or has_partial_derivs,
308  "Must have either fluxes or partial derivatives in a "
309  "DG evolution scheme. This means the evolution system struct (usually in "
310  "Evolution/Systems/YourSystem/System.hpp) being used does not specify "
311  "any flux_variables or gradient_variables. Make sure the type aliases "
312  "are defined, and that at least one of them is a non-empty list of "
313  "tags.");
314 
315  using system = typename Metavariables::system;
316  using variables_tag = typename system::variables_tag;
317  using variables_tags = typename variables_tag::tags_list;
318  using partial_derivative_tags = typename system::gradient_variables;
319  using flux_variables = typename system::flux_variables;
320 
321  const Mesh<Dim>& mesh = db::get<::domain::Tags::Mesh<Dim>>(*box);
322  const auto& logical_to_inertial_inverse_jacobian = db::get<
324  *box);
325  const auto& evolved_vars = db::get<variables_tag>(*box);
326 
327  // Compute d_i u_\alpha for nonconservative products
328  if constexpr (has_partial_derivs) {
329  partial_derivatives<partial_derivative_tags>(
330  partial_derivs, evolved_vars, mesh,
331  logical_to_inertial_inverse_jacobian);
332  }
333 
334  // Compute volume du/dt and fluxes
335 
336  // Compute volume terms that are unrelated to moving meshes
337  db::mutate<db::add_tag_prefix<::Tags::dt, variables_tag>>(
338  box,
339  [&mesh, &partial_derivs, &temporaries, &volume_fluxes](
340  const gsl::not_null<Variables<
342  dt_vars_ptr,
343  const auto&... args) noexcept {
344  // Silence compiler warnings since we genuinely don't always need all
345  // the vars but sometimes do. This warning shows up with empty parameter
346  // packs, which means the packs aren't "used" below in the
347  // ComputeVolumeTimeDerivatives::apply call.
348  (void)partial_derivs;
349  (void)volume_fluxes;
350  (void)temporaries;
351 
352  // For now just zero dt_vars. If this is a performance bottle neck we
353  // can re-evaluate in the future.
354  dt_vars_ptr->initialize(mesh.number_of_grid_points(), 0.0);
355 
357  make_not_null(&get<::Tags::dt<VariablesTags>>(*dt_vars_ptr))...,
358  make_not_null(&get<FluxVariablesTags>(*volume_fluxes))...,
359  make_not_null(&get<TemporaryTags>(*temporaries))...,
360  get<PartialDerivTags>(*partial_derivs)..., args...);
361  },
362  db::get<TimeDerivativeArgumentTags>(*box)...);
363 
364  // Add volume terms for moving meshes
365  if (const auto& mesh_velocity =
367  static_cast<bool>(mesh_velocity)) {
368  db::mutate<db::add_tag_prefix<::Tags::dt, variables_tag>>(
369  box,
370  [&evolved_vars, &mesh_velocity, &volume_fluxes](
371  const gsl::not_null<Variables<db::wrap_tags_in<
372  ::Tags::dt, typename variables_tag::tags_list>>*>
373  dt_vars_ptr,
374  const boost::optional<Scalar<DataVector>>&
375  div_mesh_velocity) noexcept {
376  tmpl::for_each<flux_variables>([&div_mesh_velocity, &dt_vars_ptr,
377  &evolved_vars, &mesh_velocity,
378  &volume_fluxes](auto tag_v) noexcept {
379  // Modify fluxes for moving mesh
380  using var_tag = typename decltype(tag_v)::type;
381  using flux_var_tag =
384  auto& flux_var = get<flux_var_tag>(*volume_fluxes);
385  // Loop over all independent components of flux_var
386  for (size_t flux_var_storage_index = 0;
387  flux_var_storage_index < flux_var.size();
388  ++flux_var_storage_index) {
389  // Get the flux variable's tensor index, e.g. (i,j) for a F^i of
390  // the spatial velocity (or some other spatial tensor).
391  const auto flux_var_tensor_index =
392  flux_var.get_tensor_index(flux_var_storage_index);
393  // Remove the first index from the flux tensor index, gets back
394  // (j)
395  const auto var_tensor_index =
396  all_but_specified_element_of(flux_var_tensor_index, 0);
397  // Set flux_index to (i)
398  const size_t flux_index = gsl::at(flux_var_tensor_index, 0);
399 
400  // We now need to index flux(i,j) -= u(j) * v_g(i)
401  flux_var[flux_var_storage_index] -=
402  get<var_tag>(evolved_vars).get(var_tensor_index) *
403  mesh_velocity->get(flux_index);
404  }
405 
406  // Modify time derivative (i.e. source terms) for moving mesh
407  auto& dt_var = get<::Tags::dt<var_tag>>(*dt_vars_ptr);
408  for (size_t dt_var_storage_index = 0;
409  dt_var_storage_index < dt_var.size(); ++dt_var_storage_index) {
410  // This is S -> S - u d_i v^i_g
411  dt_var[dt_var_storage_index] -=
412  get<var_tag>(evolved_vars)[dt_var_storage_index] *
413  get(*div_mesh_velocity);
414  }
415  });
416  },
417  db::get<::domain::Tags::DivMeshVelocity>(*box));
418 
419  // We add the mesh velocity to all equations that don't have flux terms.
420  // This doesn't need to be equal to the equations that have partial
421  // derivatives. For example, the scalar field evolution equation in
422  // first-order form does not have any partial derivatives but still needs
423  // the velocity term added. This is because the velocity term arises from
424  // transforming the time derivative.
425  using non_flux_tags = tmpl::list_difference<variables_tags, flux_variables>;
426 
427  db::mutate<db::add_tag_prefix<::Tags::dt, variables_tag>>(
428  box, [&mesh_velocity, &partial_derivs](
429  const gsl::not_null<Variables<db::wrap_tags_in<
430  ::Tags::dt, typename variables_tag::tags_list>>*>
431  dt_vars) noexcept {
432  tmpl::for_each<non_flux_tags>(
433  [&dt_vars, &mesh_velocity,
434  &partial_derivs](auto var_tag_v) noexcept {
435  using var_tag = typename decltype(var_tag_v)::type;
436  using dt_var_tag = ::Tags::dt<var_tag>;
437  using deriv_var_tag =
439 
440  const auto& deriv_var = get<deriv_var_tag>(*partial_derivs);
441  auto& dt_var = get<dt_var_tag>(*dt_vars);
442 
443  // Loop over all independent components of the derivative of the
444  // variable.
445  for (size_t deriv_var_storage_index = 0;
446  deriv_var_storage_index < deriv_var.size();
447  ++deriv_var_storage_index) {
448  // We grab the `deriv_tensor_index`, which would be e.g.
449  // `(i, a, b)`, so `(0, 2, 3)`
450  const auto deriv_var_tensor_index =
451  deriv_var.get_tensor_index(deriv_var_storage_index);
452  // Then we drop the derivative index (the first entry) to get
453  // `(a, b)` (or `(2, 3)`)
454  const auto dt_var_tensor_index =
455  all_but_specified_element_of(deriv_var_tensor_index, 0);
456  // Set `deriv_index` to `i` (or `0` in the example)
457  const size_t deriv_index = gsl::at(deriv_var_tensor_index, 0);
458  dt_var.get(dt_var_tensor_index) +=
459  mesh_velocity->get(deriv_index) *
460  deriv_var[deriv_var_storage_index];
461  }
462  });
463  });
464  }
465 
466  // Add the flux divergence term to du_\alpha/dt, which must be done
467  // after the corrections for the moving mesh are made.
468  if constexpr (has_fluxes) {
469  db::mutate<db::add_tag_prefix<::Tags::dt, variables_tag>>(
470  box,
471  [dg_formulation, &logical_to_inertial_inverse_jacobian, &mesh,
472  &volume_fluxes](const gsl::not_null<Variables<db::wrap_tags_in<
473  ::Tags::dt, typename variables_tag::tags_list>>*>
474  dt_vars_ptr) noexcept {
475  if (dg_formulation == ::dg::Formulation::StrongInertial) {
476  const Variables<tmpl::list<::Tags::div<FluxVariablesTags>...>>
477  div_fluxes = divergence(*volume_fluxes, mesh,
478  logical_to_inertial_inverse_jacobian);
479  tmpl::for_each<flux_variables>([&dt_vars_ptr, &div_fluxes](
480  auto var_tag_v) noexcept {
481  using var_tag = typename decltype(var_tag_v)::type;
482  auto& dt_var = get<::Tags::dt<var_tag>>(*dt_vars_ptr);
483  const auto& div_flux_var = get<::Tags::div<
485  div_fluxes);
486  for (size_t storage_index = 0; storage_index < dt_var.size();
487  ++storage_index) {
488  dt_var[storage_index] -= div_flux_var[storage_index];
489  }
490  });
491  } else if (dg_formulation == ::dg::Formulation::WeakInertial) {
492  ERROR("Weak inertial DG formulation not yet implemented");
493  } else {
494  ERROR("Unsupported DG formulation: " << dg_formulation);
495  }
496  });
497  } else {
498  (void)dg_formulation;
499  }
500 }
501 
502 template <typename Metavariables>
503 template <typename DbTagsList>
504 void ComputeTimeDerivative<Metavariables>::
505  boundary_terms_nonconservative_products(
506  const gsl::not_null<db::DataBox<DbTagsList>*> box) noexcept {
507  using system = typename Metavariables::system;
508  using variables_tag = typename system::variables_tag;
509  using DirectionsTag =
511 
512  using interface_normal_dot_fluxes_tag = domain::Tags::Interface<
514 
515  using interface_compute_item_argument_tags = tmpl::transform<
516  typename Metavariables::system::normal_dot_fluxes::argument_tags,
517  tmpl::bind<domain::Tags::Interface, DirectionsTag, tmpl::_1>>;
518 
520  tmpl::list<interface_normal_dot_fluxes_tag>,
521  tmpl::push_front<interface_compute_item_argument_tags, DirectionsTag>>(
522  [](const auto boundary_fluxes_ptr,
524  internal_directions,
525  const auto&... tensors) noexcept {
526  for (const auto& direction : internal_directions) {
527  apply_flux(make_not_null(&boundary_fluxes_ptr->at(direction)),
528  tensors.at(direction)...);
529  }
530  },
531  box);
532 }
533 
534 template <typename Metavariables>
535 template <size_t VolumeDim, typename BoundaryScheme, typename DbTagsList>
536 void ComputeTimeDerivative<Metavariables>::
537  fill_mortar_data_for_internal_boundaries(
538  const gsl::not_null<db::DataBox<DbTagsList>*> box) noexcept {
539  using temporal_id_tag = typename BoundaryScheme::temporal_id_tag;
540  using all_mortar_data_tag =
542  using boundary_data_computer =
543  typename BoundaryScheme::boundary_data_computer;
544 
545  // Collect data on element interfaces
546  auto boundary_data_on_interfaces =
547  interface_apply<domain::Tags::InternalDirections<VolumeDim>,
548  boundary_data_computer>(*box);
549 
550  // Project collected data to all internal mortars and store in DataBox
551  const auto& element = db::get<domain::Tags::Element<VolumeDim>>(*box);
552  const auto& face_meshes =
553  get<domain::Tags::Interface<domain::Tags::InternalDirections<VolumeDim>,
554  domain::Tags::Mesh<VolumeDim - 1>>>(*box);
555  const auto& mortar_meshes =
556  get<::Tags::Mortars<domain::Tags::Mesh<VolumeDim - 1>, VolumeDim>>(*box);
557  const auto& mortar_sizes =
558  get<::Tags::Mortars<::Tags::MortarSize<VolumeDim - 1>, VolumeDim>>(*box);
559  const auto& temporal_id = get<temporal_id_tag>(*box);
560  for (const auto& [direction, neighbors] : element.neighbors()) {
561  const auto& face_mesh = face_meshes.at(direction);
562  for (const auto& neighbor : neighbors) {
563  const auto mortar_id = std::make_pair(direction, neighbor);
564  const auto& mortar_mesh = mortar_meshes.at(mortar_id);
565  const auto& mortar_size = mortar_sizes.at(mortar_id);
566 
567  // Project the data from the face to the mortar.
568  // Where no projection is necessary we `std::move` the data directly to
569  // avoid a copy. We can't move the data or modify it in-place when
570  // projecting, because in that case the face may touch two mortars so we
571  // need to keep the data around.
572  auto boundary_data_on_mortar =
574  ? boundary_data_on_interfaces.at(direction).project_to_mortar(
575  face_mesh, mortar_mesh, mortar_size)
576  : std::move(boundary_data_on_interfaces.at(direction));
577 
578  // Store the boundary data on this side of the mortar
579  db::mutate<all_mortar_data_tag>(
580  box, [&mortar_id, &temporal_id, &boundary_data_on_mortar](
582  all_mortar_data) noexcept {
583  all_mortar_data->at(mortar_id).local_insert(
584  temporal_id, std::move(boundary_data_on_mortar));
585  });
586  }
587  }
588 }
589 } // namespace evolution::dg::Actions
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
DataBoxTag.hpp
utility
Frame::Inertial
Definition: IndexType.hpp:44
Divergence.hpp
get
constexpr Tag::type & get(Variables< TagList > &v) noexcept
Return Tag::type pointing into the contiguous array.
Definition: Variables.hpp:639
Parallel::GlobalCache
Definition: ElementReceiveInterpPoints.hpp:16
unordered_set
Tags::MortarSize
Definition: Tags.hpp:49
Tags.hpp
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:52
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:108
tuple
db::get
const auto & get(const DataBox< TagList > &box) noexcept
Retrieve the item with tag Tag from the DataBox.
Definition: DataBox.hpp:1116
Tags::Flux
Prefix indicating a flux.
Definition: Prefixes.hpp:40
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
ERROR
#define ERROR(m)
prints an error message to the standard error stream and aborts the program.
Definition: Error.hpp:36
db::apply
constexpr auto apply(F &&f, const DataBox< BoxTags > &box, Args &&... args) noexcept
Apply the invokable f with argument Tags TagsList from DataBox box
Definition: DataBox.hpp:1424
Tags::div
Prefix indicating the divergence.
Definition: Divergence.hpp:45
evolution::dg::Actions::ComputeTimeDerivative
Computes the time derivative for a DG time step.
Definition: ComputeTimeDerivative.hpp:133
DataBox.hpp
db::mutate_apply
constexpr auto mutate_apply(F &&f, const gsl::not_null< DataBox< BoxTags > * > box, Args &&... args) noexcept(DataBox_detail::check_mutate_apply_mutate_tags(BoxTags{}, MutateTags{}) and DataBox_detail::check_mutate_apply_argument_tags(BoxTags{}, ArgumentTags{}) and noexcept(DataBox_detail::mutate_apply(f, box, MutateTags{}, ArgumentTags{}, std::forward< Args >(args)...)))
Apply the invokable f mutating items MutateTags and taking as additional arguments ArgumentTags and a...
Definition: DataBox.hpp:1640
domain::Tags::InverseJacobian< Dim, Frame::Logical, Frame::Inertial >
db::item_type
typename DataBox_detail::item_type_impl< TagList, Tag >::type item_type
Get the type that can be written to the Tag. If it is a base tag then a TagList must be passed as a s...
Definition: DataBoxTag.hpp:247
evolution::dg::Actions
Actions for using the discontinuous Galerkin to evolve hyperbolic partial differential equations.
Definition: ComputeTimeDerivative.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:127
tuples::TaggedTuple
An associative container that is indexed by structs.
Definition: TaggedTuple.hpp:271
Variables.hpp
Mesh
Holds the number of grid points, basis, and quadrature in each direction of the computational grid.
Definition: Mesh.hpp:47
Tags::dt
Prefix indicating a time derivative.
Definition: Prefixes.hpp:29
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:75
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
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 list of Tags by wrapping each tag in TagList using the Wrapper.
Definition: PrefixHelpers.hpp:30
Tags::deriv
Prefix indicating spatial derivatives.
Definition: PartialDerivatives.hpp:53
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
Mesh.hpp
PartialDerivatives.hpp
dg::Formulation
Formulation
The DG formulation to use.
Definition: Formulation.hpp:26
dg::mortar_mesh
Mesh< Dim > mortar_mesh(const Mesh< Dim > &face_mesh1, const Mesh< Dim > &face_mesh2) noexcept
Definition: MortarHelpers.cpp:19
Prefixes.hpp
db::DataBox
Definition: InterpolationTargetWedgeSectionTorus.hpp:24
Parallel
Contains functions that forward to Charm++ parallel functions.
Definition: ElementReceiveInterpPoints.hpp:14
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
gsl::not_null
Require a pointer to not be a nullptr
Definition: Gsl.hpp:183
Tags::Mortars
Definition: Tags.hpp:38