Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : #pragma once
5 :
6 : #include <limits>
7 : #include <optional>
8 : #include <tuple>
9 : #include <type_traits>
10 : #include <unordered_set>
11 : #include <utility>
12 :
13 : #include "DataStructures/DataBox/DataBox.hpp"
14 : #include "DataStructures/DataBox/PrefixHelpers.hpp"
15 : #include "DataStructures/DataBox/Prefixes.hpp"
16 : #include "DataStructures/DataVector.hpp"
17 : #include "DataStructures/Tensor/Tensor.hpp"
18 : #include "DataStructures/Variables.hpp"
19 : #include "DataStructures/VariablesTag.hpp"
20 : #include "Domain/CoordinateMaps/Tags.hpp"
21 : #include "Domain/Creators/Tags/ExternalBoundaryConditions.hpp"
22 : #include "Domain/InterfaceHelpers.hpp"
23 : #include "Domain/Structure/Direction.hpp"
24 : #include "Domain/Structure/DirectionMap.hpp"
25 : #include "Domain/Structure/OrientationMapHelpers.hpp"
26 : #include "Domain/Tags.hpp"
27 : #include "Domain/TagsTimeDependent.hpp"
28 : #include "Evolution/BoundaryCorrection.hpp"
29 : #include "Evolution/BoundaryCorrectionTags.hpp"
30 : #include "Evolution/DiscontinuousGalerkin/Actions/BoundaryConditionsImpl.hpp"
31 : #include "Evolution/DiscontinuousGalerkin/Actions/ComputeTimeDerivativeHelpers.hpp"
32 : #include "Evolution/DiscontinuousGalerkin/Actions/InternalMortarDataImpl.hpp"
33 : #include "Evolution/DiscontinuousGalerkin/Actions/NormalCovectorAndMagnitude.hpp"
34 : #include "Evolution/DiscontinuousGalerkin/Actions/PackageDataImpl.hpp"
35 : #include "Evolution/DiscontinuousGalerkin/Actions/VolumeTermsImpl.hpp"
36 : #include "Evolution/DiscontinuousGalerkin/BoundaryData.hpp"
37 : #include "Evolution/DiscontinuousGalerkin/InboxTags.hpp"
38 : #include "Evolution/DiscontinuousGalerkin/InterfaceDataPolicy.hpp"
39 : #include "Evolution/DiscontinuousGalerkin/InterpolatedBoundaryData.hpp"
40 : #include "Evolution/DiscontinuousGalerkin/MortarData.hpp"
41 : #include "Evolution/DiscontinuousGalerkin/MortarDataHolder.hpp"
42 : #include "Evolution/DiscontinuousGalerkin/MortarInfo.hpp"
43 : #include "Evolution/DiscontinuousGalerkin/MortarTags.hpp"
44 : #include "Evolution/DiscontinuousGalerkin/NormalVectorTags.hpp"
45 : #include "Evolution/DiscontinuousGalerkin/TimeSteppingPolicy.hpp"
46 : #include "Evolution/DiscontinuousGalerkin/UsingSubcell.hpp"
47 : #include "NumericalAlgorithms/DiscontinuousGalerkin/Formulation.hpp"
48 : #include "NumericalAlgorithms/DiscontinuousGalerkin/MortarHelpers.hpp"
49 : #include "NumericalAlgorithms/DiscontinuousGalerkin/Tags.hpp"
50 : #include "NumericalAlgorithms/DiscontinuousGalerkin/Tags/Formulation.hpp"
51 : #include "NumericalAlgorithms/LinearOperators/Divergence.hpp"
52 : #include "NumericalAlgorithms/LinearOperators/PartialDerivatives.hpp"
53 : #include "NumericalAlgorithms/Spectral/Basis.hpp"
54 : #include "NumericalAlgorithms/Spectral/BoundaryInterpolationMatrices.hpp"
55 : #include "NumericalAlgorithms/Spectral/Mesh.hpp"
56 : #include "NumericalAlgorithms/Spectral/Quadrature.hpp"
57 : #include "Parallel/AlgorithmExecution.hpp"
58 : #include "Parallel/ArrayCollection/IsDgElementCollection.hpp"
59 : #include "Parallel/ArrayCollection/SendDataToElement.hpp"
60 : #include "Parallel/GlobalCache.hpp"
61 : #include "Parallel/Invoke.hpp"
62 : #include "Time/BoundaryHistory.hpp"
63 : #include "Time/ChangeStepSize.hpp"
64 : #include "Utilities/Algorithm.hpp"
65 : #include "Utilities/Gsl.hpp"
66 : #include "Utilities/TMPL.hpp"
67 :
68 : /// \cond
69 : namespace Tags {
70 : template <typename Tag>
71 : struct HistoryEvolvedVariables;
72 : struct TimeStepId;
73 : } // namespace Tags
74 : namespace evolution::dg::Tags {
75 : template <size_t Dim>
76 : struct MortarInfo;
77 : } // namespace evolution::dg::Tags
78 :
79 : namespace evolution::dg::subcell {
80 : // We use a forward declaration instead of including a header file to avoid
81 : // coupling to the DG-subcell libraries for executables that don't use subcell.
82 : template <typename Metavariables, typename DbTagsList, size_t Dim>
83 : void prepare_neighbor_data(
84 : gsl::not_null<DirectionMap<Dim, DataVector>*>
85 : all_neighbor_data_for_reconstruction,
86 : gsl::not_null<Mesh<Dim>*> ghost_data_mesh,
87 : gsl::not_null<db::DataBox<DbTagsList>*> box,
88 : [[maybe_unused]] const Variables<db::wrap_tags_in<
89 : ::Tags::Flux, typename Metavariables::system::flux_variables,
90 : tmpl::size_t<Dim>, Frame::Inertial>>& volume_fluxes);
91 : template <typename DbTagsList>
92 : int get_tci_decision(const db::DataBox<DbTagsList>& box);
93 : } // namespace evolution::dg::subcell
94 : namespace tuples {
95 : template <typename...>
96 : class TaggedTuple;
97 : } // namespace tuples
98 : /// \endcond
99 :
100 : namespace evolution::dg::Actions {
101 : namespace detail {
102 : template <typename T>
103 : struct get_dg_package_temporary_tags {
104 : using type = typename T::dg_package_data_temporary_tags;
105 : };
106 : template <typename T>
107 : struct get_dg_package_field_tags {
108 : using type = typename T::dg_package_field_tags;
109 : };
110 : template <typename System, typename T>
111 : struct get_primitive_tags_for_face {
112 : using type = typename get_primitive_vars<
113 : System::has_primitive_and_conservative_vars>::template f<T>;
114 : };
115 : } // namespace detail
116 :
117 : /*!
118 : * \brief Computes the time derivative for a DG time step.
119 : *
120 : * Computes the volume fluxes, the divergence of the fluxes and all additional
121 : * interior contributions to the time derivatives (both nonconservative products
122 : * and source terms). The internal mortar data is also computed.
123 : *
124 : * The general first-order hyperbolic evolution equation solved for conservative
125 : * systems is:
126 : *
127 : * \f{align*}{
128 : * \frac{\partial u_\alpha}{\partial \hat{t}}
129 : * + \partial_{i}
130 : * \left(F^i_\alpha - v^i_g u_\alpha\right)
131 : * = S_\alpha-u_\alpha\partial_i v^i_g,
132 : * \f}
133 : *
134 : * where \f$F^i_{\alpha}\f$ are the fluxes when the mesh isn't moving,
135 : * \f$v^i_g\f$ is the velocity of the mesh, \f$u_{\alpha}\f$ are the evolved
136 : * variables, \f$S_{\alpha}\f$ are the source terms, \f$\hat{t}\f$ is the
137 : * time in the logical frame, \f$t\f$ is the time in the inertial frame, hatted
138 : * indices correspond to logical frame quantites, and unhatted indices to
139 : * inertial frame quantities (e.g. \f$\partial_i\f$ is the derivative with
140 : * respect to the inertial coordinates). For evolution equations that do not
141 : * have any fluxes and only nonconservative products we evolve:
142 : *
143 : * \f{align*}{
144 : * \frac{\partial u_\alpha}{\partial \hat{t}}
145 : * +\left(B^i_{\alpha\beta}-v^i_g \delta_{\alpha\beta}
146 : * \right)\partial_{i}u_\beta = S_\alpha.
147 : * \f}
148 : *
149 : * Finally, for equations with both conservative terms and nonconservative
150 : * products we use:
151 : *
152 : * \f{align*}{
153 : * \frac{\partial u_\alpha}{\partial \hat{t}}
154 : * + \partial_{i}
155 : * \left(F^i_\alpha - v^i_g u_\alpha\right)
156 : * +B^i_{\alpha\beta}\partial_{i}u_\beta
157 : * = S_\alpha-u_\alpha\partial_i v^i_g,
158 : * \f}
159 : *
160 : * where \f$B^i_{\alpha\beta}\f$ is the matrix for the nonconservative products.
161 : *
162 : * ### Volume Terms
163 : *
164 : * The mesh velocity is added to the flux automatically if the mesh is moving.
165 : * That is,
166 : *
167 : * \f{align*}{
168 : * F^i_{\alpha}\to F^i_{\alpha}-v^i_{g} u_{\alpha}
169 : * \f}
170 : *
171 : * The source terms are also altered automatically by adding:
172 : *
173 : * \f{align*}{
174 : * -u_\alpha \partial_i v^i_g,
175 : * \f}
176 : *
177 : * For systems with equations that only contain nonconservative products, the
178 : * following mesh velocity is automatically added to the time derivative:
179 : *
180 : * \f{align*}{
181 : * v^i_g \partial_i u_\alpha,
182 : * \f}
183 : *
184 : * \note The term is always added in the `Frame::Inertial` frame, and the plus
185 : * sign arises because we add it to the time derivative.
186 : *
187 : * \warning The mesh velocity terms are added to the time derivatives before
188 : * invoking the boundary conditions. This means that the time derivatives passed
189 : * to the boundary conditions are with respect to \f$\hat{t}\f$, not \f$t\f$.
190 : * This is especially important in the TimeDerivative/Bjorhus boundary
191 : * conditions.
192 : *
193 : * Here are examples of the `TimeDerivative` struct used to compute the volume
194 : * time derivative. This struct is what the type alias
195 : * `System::compute_volume_time_derivative` points to. The time derivatives are
196 : * as `gsl::not_null` first, then the temporary tags as `gsl::not_null`,
197 : * followed by the `argument_tags`. These type aliases are given by
198 : *
199 : * \snippet ComputeTimeDerivativeImpl.tpp dt_ta
200 : *
201 : * for the examples. For a conservative system without primitives the `apply`
202 : * function would look like
203 : *
204 : * \snippet ComputeTimeDerivativeImpl.tpp dt_con
205 : *
206 : * For a nonconservative system it would be
207 : *
208 : * \snippet ComputeTimeDerivativeImpl.tpp dt_nc
209 : *
210 : * And finally, for a mixed conservative-nonconservative system with primitive
211 : * variables
212 : *
213 : * \snippet ComputeTimeDerivativeImpl.tpp dt_mp
214 : *
215 : * In addition to each variable being passed individually, if the time
216 : * derivative struct inherits from `evolution::PassVariables`, then the time
217 : * derivatives, fluxes, and temporaries are passed as
218 : * `gsl::not_null<Variables<...>>`. This is useful for systems where
219 : * additional quantities are sometimes evolved, and just generally nice for
220 : * keeping the number of arguments reasonable. Below are the above examples
221 : * but with `Variables` being passed.
222 : *
223 : * \snippet ComputeTimeDerivativeImpl.tpp dt_con_variables
224 : *
225 : * \snippet ComputeTimeDerivativeImpl.tpp dt_nc_variables
226 : *
227 : * \snippet ComputeTimeDerivativeImpl.tpp dt_mp_variables
228 : *
229 : * Uses:
230 : * - System:
231 : * - `variables_tag`
232 : * - `flux_variables`
233 : * - `gradient_variables`
234 : * - `compute_volume_time_derivative_terms`
235 : *
236 : * - DataBox:
237 : * - Items in `system::compute_volume_time_derivative_terms::argument_tags`
238 : * - `domain::Tags::MeshVelocity<Metavariables::volume_dim>`
239 : * - `Metavariables::system::variables_tag`
240 : * - `Metavariables::system::flux_variables`
241 : * - `Metavariables::system::gradient_variables`
242 : * - `domain::Tags::DivMeshVelocity`
243 : * - `DirectionsTag`,
244 : * - Required interface items for `Metavariables::system::normal_dot_fluxes`
245 : *
246 : * DataBox changes:
247 : * - Adds: nothing
248 : * - Removes: nothing
249 : * - Modifies:
250 : * - db::add_tag_prefix<Tags::Flux, variables_tag,
251 : * tmpl::size_t<system::volume_dim>, Frame::Inertial>
252 : * - `Tags::dt<system::variable_tags>`
253 : * - Tags::Interface<
254 : * DirectionsTag, db::add_tag_prefix<Tags::NormalDotFlux, variables_tag>>
255 : * - `Tags::Mortars<typename BoundaryScheme::mortar_data_tag, VolumeDim>`
256 : *
257 : * ### Internal Boundary Terms
258 : *
259 : * Internal boundary terms must be derived from
260 : * `evolution::BoundaryCorrection`. Each concrete boundary correction
261 : * must specify:
262 : *
263 : * - type alias `dg_package_field_tags`. These are what will be returned by
264 : * `gsl::not_null` from the `dg_package_data` member function.
265 : *
266 : * - type alias `dg_package_data_temporary_tags`. These are temporary tags
267 : * that are projected to the face and then passed to the `dg_package_data`
268 : * function.
269 : *
270 : * - type alias `dg_package_data_primitive_tags`. These are the primitive
271 : * variables (if any) that are projected to the face and then passed to
272 : * `dg_package_data`.
273 : *
274 : * - type alias `dg_package_data_volume_tags`. These are tags that are not
275 : * projected to the interface and are retrieved directly from the `DataBox`.
276 : * The equation of state for hydrodynamics systems is an example of what
277 : * would be a "volume tag".
278 : *
279 : * A `static constexpr bool need_normal_vector` must be specified. If `true`
280 : * then the normal vector is computed from the normal covector. This is
281 : * currently not implemented.
282 : *
283 : * The `dg_package_data` function takes as arguments `gsl::not_null` of the
284 : * `dg_package_field_tags`, then the projected evolved variables, the
285 : * projected fluxes, the projected temporaries, the projected primitives, the
286 : * unit normal covector, mesh velocity, normal dotted into the mesh velocity,
287 : * the `volume_tags`, and finally the `dg::Formulation`. The `dg_package_data`
288 : * function must compute all ingredients for the boundary correction, including
289 : * mesh-velocity-corrected characteristic speeds. However, the projected fluxes
290 : * passed in are \f$F^i - u v^i_g\f$ (the mesh velocity term is already
291 : * included). The `dg_package_data` function must also return a `double` that is
292 : * the maximum absolute characteristic speed over the entire face. This will be
293 : * used for checking that the time step doesn't violate the CFL condition.
294 : *
295 : * Here is an example of the type aliases and `bool`:
296 : *
297 : * \snippet ComputeTimeDerivativeImpl.tpp bt_ta
298 : *
299 : * The normal vector requirement is:
300 : *
301 : * \snippet ComputeTimeDerivativeImpl.tpp bt_nnv
302 : *
303 : * For a conservative system with primitive variables and using the `TimeStepId`
304 : * as a volume tag the `dg_package_data` function looks like:
305 : *
306 : * \snippet ComputeTimeDerivativeImpl.tpp bt_cp
307 : *
308 : * For a mixed conservative-nonconservative system with primitive variables and
309 : * using the `TimeStepId` as a volume tag the `dg_package_data` function looks
310 : * like:
311 : *
312 : * \snippet ComputeTimeDerivativeImpl.tpp bt_mp
313 : *
314 : * Uses:
315 : * - System:
316 : * - `boundary_correction`
317 : * - `variables_tag`
318 : * - `flux_variables`
319 : * - `gradients_tags`
320 : * - `compute_volume_time_derivative`
321 : * - `has_primitive_and_conservative_vars`
322 : * - `primitive_variables_tag` if system has primitive variables
323 : *
324 : * - DataBox:
325 : * - `domain::Tags::Element<Dim>`
326 : * - `domain::Tags::Mesh<Dim>`
327 : * - `evolution::dg::Tags::MortarMesh<Dim>`
328 : * - `evolution::dg::Tags::MortarData<Dim>`
329 : * - `Tags::TimeStepId`
330 : * - \code{.cpp}
331 : * domain::Tags::Interface<domain::Tags::InternalDirections<Dim>,
332 : * domain::Tags::Mesh<Dim - 1>>
333 : * \endcode
334 : * - \code{.cpp}
335 : * domain::Tags::Interface<
336 : * domain::Tags::InternalDirections<Dim>,
337 : * ::Tags::Normalized<
338 : * domain::Tags::UnnormalizedFaceNormal<Dim, Frame::Inertial>>>
339 : * \endcode
340 : * - \code{.cpp}
341 : * domain::Tags::Interface<
342 : * domain::Tags::InternalDirections<Dim>,
343 : * domain::Tags::MeshVelocity<Dim, Frame::Inertial>>
344 : * \endcode
345 : * - `Metavariables::system::variables_tag`
346 : * - `Metavariables::system::flux_variables`
347 : * - `Metavariables::system::primitive_tags` if exists
348 : * - boundary correction `dg_package_data_volume_tags`
349 : *
350 : * DataBox changes:
351 : * - Adds: nothing
352 : * - Removes: nothing
353 : * - Modifies:
354 : * - `evolution::dg::Tags::MortarData<Dim>`
355 : */
356 : template <size_t Dim, typename EvolutionSystem, typename DgStepChoosers,
357 : bool LocalTimeStepping, bool UseNodegroupDgElements>
358 1 : struct ComputeTimeDerivative {
359 0 : using inbox_tags =
360 : tmpl::list<evolution::dg::Tags::BoundaryCorrectionAndGhostCellsInbox<
361 : Dim, UseNodegroupDgElements>>;
362 0 : using const_global_cache_tags = tmpl::append<
363 : tmpl::list<::dg::Tags::Formulation, evolution::Tags::BoundaryCorrection,
364 : domain::Tags::ExternalBoundaryConditions<Dim>>,
365 : tmpl::conditional_t<
366 : LocalTimeStepping,
367 : typename ChangeStepSize<DgStepChoosers>::const_global_cache_tags,
368 : tmpl::list<>>>;
369 :
370 : template <typename DbTagsList, typename... InboxTags, typename ArrayIndex,
371 : typename ActionList, typename ParallelComponent,
372 : typename Metavariables>
373 0 : static Parallel::iterable_action_return_t apply(
374 : db::DataBox<DbTagsList>& box,
375 : tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
376 : Parallel::GlobalCache<Metavariables>& cache,
377 : const ArrayIndex& /*array_index*/, ActionList /*meta*/,
378 : const ParallelComponent* /*meta*/); // NOLINT const
379 :
380 : private:
381 : template <typename ParallelComponent, typename DbTagsList,
382 : typename Metavariables>
383 0 : static void send_data_for_fluxes(
384 : gsl::not_null<Parallel::GlobalCache<Metavariables>*> cache,
385 : gsl::not_null<db::DataBox<DbTagsList>*> box,
386 : [[maybe_unused]] const Variables<db::wrap_tags_in<
387 : ::Tags::Flux, typename EvolutionSystem::flux_variables,
388 : tmpl::size_t<Dim>, Frame::Inertial>>& volume_fluxes);
389 : };
390 :
391 : template <size_t Dim, typename EvolutionSystem, typename DgStepChoosers,
392 : bool LocalTimeStepping, bool UseNodegroupDgElements>
393 : template <typename DbTagsList, typename... InboxTags, typename ArrayIndex,
394 : typename ActionList, typename ParallelComponent,
395 : typename Metavariables>
396 : Parallel::iterable_action_return_t
397 : ComputeTimeDerivative<Dim, EvolutionSystem, DgStepChoosers, LocalTimeStepping,
398 : UseNodegroupDgElements>::
399 : apply(db::DataBox<DbTagsList>& box,
400 : tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
401 : Parallel::GlobalCache<Metavariables>& cache,
402 : const ArrayIndex& /*array_index*/, ActionList /*meta*/,
403 : const ParallelComponent* const /*meta*/) { // NOLINT const
404 : static_assert(UseNodegroupDgElements ==
405 : Parallel::is_dg_element_collection_v<ParallelComponent>,
406 : "The action ComputeTimeDerivative is told by the "
407 : "template parameter UseNodegroupDgElements that it is being "
408 : "used with a DgElementCollection, but the ParallelComponent "
409 : "is not a DgElementCollection. You need to change the "
410 : "template parameter on the ComputeTimeDerivative action "
411 : "in your action list.");
412 :
413 : using variables_tag = typename EvolutionSystem::variables_tag;
414 : using dt_variables_tag = db::add_tag_prefix<::Tags::dt, variables_tag>;
415 : using partial_derivative_tags = typename EvolutionSystem::gradient_variables;
416 : using flux_variables = typename EvolutionSystem::flux_variables;
417 : using compute_volume_time_derivative_terms =
418 : typename EvolutionSystem::compute_volume_time_derivative_terms;
419 :
420 : const Mesh<Dim>& mesh = db::get<::domain::Tags::Mesh<Dim>>(box);
421 : const Element<Dim>& element = db::get<domain::Tags::Element<Dim>>(box);
422 : const ::dg::Formulation dg_formulation =
423 : db::get<::dg::Tags::Formulation>(box);
424 : ASSERT(alg::all_of(mesh.basis(),
425 : [&mesh](const Spectral::Basis current_basis) {
426 : return current_basis == mesh.basis(0);
427 : }) or
428 : element.topologies() != domain::topologies::hypercube<Dim>,
429 : "An isotropic basis must be used in the evolution code. While "
430 : "theoretically this restriction could be lifted, the simplification "
431 : "it offers are quite substantial. Relaxing this assumption is likely "
432 : "to require quite a bit of careful code refactoring and debugging.");
433 : ASSERT(alg::all_of(mesh.quadrature(),
434 : [&mesh](const Spectral::Quadrature current_quadrature) {
435 : return current_quadrature == mesh.quadrature(0);
436 : }) or
437 : element.topologies() != domain::topologies::hypercube<Dim>,
438 : "An isotropic quadrature must be used in the evolution code. While "
439 : "theoretically this restriction could be lifted, the simplification "
440 : "it offers are quite substantial. Relaxing this assumption is likely "
441 : "to require quite a bit of careful code refactoring and debugging.");
442 :
443 : const auto& boundary_correction =
444 : db::get<evolution::Tags::BoundaryCorrection>(box);
445 : using derived_boundary_corrections =
446 : tmpl::at<typename Metavariables::factory_creation::factory_classes,
447 : evolution::BoundaryCorrection>;
448 :
449 : // To avoid a second allocation in internal_mortar_data, we allocate the
450 : // variables needed to construct the fields on the faces here along with
451 : // everything else. This requires us to know all the tags necessary to apply
452 : // boundary corrections. However, since we pick boundary corrections at
453 : // runtime, we just gather all possible tags from all possible boundary
454 : // corrections and lump them into the allocation. This may result in a
455 : // larger-than-necessary allocation, but it won't be that much larger.
456 : using all_dg_package_temporary_tags =
457 : tmpl::transform<derived_boundary_corrections,
458 : detail::get_dg_package_temporary_tags<tmpl::_1>>;
459 : using all_primitive_tags_for_face =
460 : tmpl::transform<derived_boundary_corrections,
461 : detail::get_primitive_tags_for_face<
462 : tmpl::pin<EvolutionSystem>, tmpl::_1>>;
463 : using fluxes_tags = db::wrap_tags_in<::Tags::Flux, flux_variables,
464 : tmpl::size_t<Dim>, Frame::Inertial>;
465 : using dg_package_data_projected_tags =
466 : tmpl::list<typename variables_tag::tags_list, fluxes_tags,
467 : all_dg_package_temporary_tags, all_primitive_tags_for_face>;
468 : using all_face_temporary_tags =
469 : tmpl::remove_duplicates<tmpl::flatten<tmpl::push_back<
470 : tmpl::list<dg_package_data_projected_tags,
471 : detail::inverse_spatial_metric_tag<EvolutionSystem>>,
472 : detail::OneOverNormalVectorMagnitude, detail::NormalVector<Dim>>>>;
473 : // To avoid additional allocations in internal_mortar_data, we provide a
474 : // buffer used to compute the packaged data before it has to be projected to
475 : // the mortar. We get all mortar tags for similar reasons as described above
476 : using all_mortar_tags = tmpl::remove_duplicates<tmpl::flatten<
477 : tmpl::transform<derived_boundary_corrections,
478 : detail::get_dg_package_field_tags<tmpl::_1>>>>;
479 :
480 : // We also don't use the number of volume mesh grid points. We instead use the
481 : // max number of grid points from each face. That way, our allocation will be
482 : // large enough to hold any face and we can reuse the allocation for each face
483 : // without having to resize it.
484 : size_t num_face_temporary_grid_points = 0;
485 : {
486 : for (const auto& [direction, neighbors_in_direction] :
487 : element.neighbors()) {
488 : (void)neighbors_in_direction;
489 : const auto face_mesh = mesh.slice_away(direction.dimension());
490 : num_face_temporary_grid_points = std::max(
491 : num_face_temporary_grid_points, face_mesh.number_of_grid_points());
492 : }
493 : }
494 :
495 : // Allocate the Variables classes needed for the time derivative
496 : // computation.
497 : //
498 : // This is factored out so that we will be able to do ADER-DG/CG where a
499 : // spacetime polynomial is constructed by solving implicit equations in time
500 : // using a Picard iteration. A high-order initial guess is needed to
501 : // efficiently construct the ADER spacetime solution. This initial guess is
502 : // obtained using continuous RK methods, and so we will want to reuse
503 : // buffers. Thus, the volume_terms function returns by reference rather than
504 : // by value.
505 : using VarsTemporaries =
506 : Variables<typename compute_volume_time_derivative_terms::temporary_tags>;
507 : using VarsFluxes =
508 : Variables<db::wrap_tags_in<::Tags::Flux, flux_variables,
509 : tmpl::size_t<Dim>, Frame::Inertial>>;
510 : using VarsPartialDerivatives =
511 : Variables<db::wrap_tags_in<::Tags::deriv, partial_derivative_tags,
512 : tmpl::size_t<Dim>, Frame::Inertial>>;
513 : using VarsDivFluxes = Variables<db::wrap_tags_in<
514 : ::Tags::div, db::wrap_tags_in<::Tags::Flux, flux_variables,
515 : tmpl::size_t<Dim>, Frame::Inertial>>>;
516 : using VarsFaceTemporaries = Variables<all_face_temporary_tags>;
517 : using DgPackagedDataVarsOnFace = Variables<all_mortar_tags>;
518 : const size_t number_of_grid_points = mesh.number_of_grid_points();
519 : const size_t buffer_size =
520 : (VarsTemporaries::number_of_independent_components +
521 : VarsFluxes::number_of_independent_components +
522 : VarsPartialDerivatives::number_of_independent_components +
523 : VarsDivFluxes::number_of_independent_components) *
524 : number_of_grid_points +
525 : // Different number of grid points. See explanation above where
526 : // num_face_temporary_grid_points is defined
527 : (VarsFaceTemporaries::number_of_independent_components +
528 : DgPackagedDataVarsOnFace::number_of_independent_components) *
529 : num_face_temporary_grid_points;
530 : auto buffer = cpp20::make_unique_for_overwrite<double[]>(buffer_size);
531 : #ifdef SPECTRE_NAN_INIT
532 : std::fill(&buffer[0], &buffer[buffer_size],
533 : std::numeric_limits<double>::signaling_NaN());
534 : #endif
535 : VarsTemporaries temporaries{
536 : &buffer[0], VarsTemporaries::number_of_independent_components *
537 : number_of_grid_points};
538 : VarsFluxes volume_fluxes{
539 : &buffer[VarsTemporaries::number_of_independent_components *
540 : number_of_grid_points],
541 : VarsFluxes::number_of_independent_components * number_of_grid_points};
542 : VarsPartialDerivatives partial_derivs{
543 : &buffer[(VarsTemporaries::number_of_independent_components +
544 : VarsFluxes::number_of_independent_components) *
545 : number_of_grid_points],
546 : VarsPartialDerivatives::number_of_independent_components *
547 : number_of_grid_points};
548 : VarsDivFluxes div_fluxes{
549 : &buffer[(VarsTemporaries::number_of_independent_components +
550 : VarsFluxes::number_of_independent_components +
551 : VarsPartialDerivatives::number_of_independent_components) *
552 : number_of_grid_points],
553 : VarsDivFluxes::number_of_independent_components * number_of_grid_points};
554 : // Lighter weight data structure than a Variables to avoid passing even more
555 : // templates to internal_mortar_data.
556 : gsl::span<double> face_temporaries = gsl::make_span<double>(
557 : &buffer[(VarsTemporaries::number_of_independent_components +
558 : VarsFluxes::number_of_independent_components +
559 : VarsPartialDerivatives::number_of_independent_components +
560 : VarsDivFluxes::number_of_independent_components) *
561 : number_of_grid_points],
562 : // Different number of grid points. See explanation above where
563 : // num_face_temporary_grid_points is defined
564 : VarsFaceTemporaries::number_of_independent_components *
565 : num_face_temporary_grid_points);
566 : gsl::span<double> packaged_data_buffer = gsl::make_span<double>(
567 : &buffer[(VarsTemporaries::number_of_independent_components +
568 : VarsFluxes::number_of_independent_components +
569 : VarsPartialDerivatives::number_of_independent_components +
570 : VarsDivFluxes::number_of_independent_components) *
571 : number_of_grid_points +
572 : VarsFaceTemporaries::number_of_independent_components *
573 : num_face_temporary_grid_points],
574 : // Different number of grid points. See explanation above where
575 : // num_face_temporary_grid_points is defined
576 : DgPackagedDataVarsOnFace::number_of_independent_components *
577 : num_face_temporary_grid_points);
578 :
579 : const Scalar<DataVector>* det_inverse_jacobian = nullptr;
580 : if constexpr (tmpl::size<flux_variables>::value != 0) {
581 : if (dg_formulation == ::dg::Formulation::WeakInertial) {
582 : det_inverse_jacobian = &db::get<
583 : domain::Tags::DetInvJacobian<Frame::ElementLogical, Frame::Inertial>>(
584 : box);
585 : }
586 : }
587 : db::mutate_apply<
588 : tmpl::list<dt_variables_tag>,
589 : typename compute_volume_time_derivative_terms::argument_tags>(
590 : [&dg_formulation, &div_fluxes, &det_inverse_jacobian,
591 : &div_mesh_velocity = db::get<::domain::Tags::DivMeshVelocity>(box),
592 : &evolved_variables = db::get<variables_tag>(box),
593 : &inertial_coordinates =
594 : db::get<domain::Tags::Coordinates<Dim, Frame::Inertial>>(box),
595 : &logical_to_inertial_inv_jacobian =
596 : db::get<::domain::Tags::InverseJacobian<Dim, Frame::ElementLogical,
597 : Frame::Inertial>>(box),
598 : &mesh, &mesh_velocity = db::get<::domain::Tags::MeshVelocity<Dim>>(box),
599 : &partial_derivs, &temporaries, &volume_fluxes](
600 : const gsl::not_null<Variables<
601 : db::wrap_tags_in<::Tags::dt, typename variables_tag::tags_list>>*>
602 : dt_vars_ptr,
603 : const auto&... time_derivative_args) {
604 : detail::volume_terms<compute_volume_time_derivative_terms>(
605 : dt_vars_ptr, make_not_null(&volume_fluxes),
606 : make_not_null(&partial_derivs), make_not_null(&temporaries),
607 : make_not_null(&div_fluxes), evolved_variables, dg_formulation, mesh,
608 : inertial_coordinates, logical_to_inertial_inv_jacobian,
609 : det_inverse_jacobian, mesh_velocity, div_mesh_velocity,
610 : time_derivative_args...);
611 : },
612 : make_not_null(&box));
613 :
614 : const Variables<detail::get_primitive_vars_tags_from_system<EvolutionSystem>>*
615 : primitive_vars{nullptr};
616 : if constexpr (EvolutionSystem::has_primitive_and_conservative_vars) {
617 : primitive_vars =
618 : &db::get<typename EvolutionSystem::primitive_variables_tag>(box);
619 : }
620 :
621 : static_assert(
622 : tmpl::all<derived_boundary_corrections, std::is_final<tmpl::_1>>::value,
623 : "All createable classes for boundary corrections must be marked "
624 : "final.");
625 : tmpl::for_each<derived_boundary_corrections>(
626 : [&boundary_correction, &box, &partial_derivs, &primitive_vars,
627 : &temporaries, &volume_fluxes, &packaged_data_buffer,
628 : &face_temporaries](auto derived_correction_v) {
629 : using DerivedCorrection =
630 : tmpl::type_from<decltype(derived_correction_v)>;
631 : if (typeid(boundary_correction) == typeid(DerivedCorrection)) {
632 : // Compute internal boundary quantities on the mortar for sides
633 : // of the element that have neighbors, i.e. they are not an
634 : // external side.
635 : // Note: this call mutates:
636 : // - evolution::dg::Tags::NormalCovectorAndMagnitude<Dim>,
637 : // - evolution::dg::Tags::MortarData<Dim>
638 : detail::internal_mortar_data<EvolutionSystem, Dim>(
639 : make_not_null(&box), make_not_null(&face_temporaries),
640 : make_not_null(&packaged_data_buffer),
641 : dynamic_cast<const DerivedCorrection&>(boundary_correction),
642 : db::get<variables_tag>(box), volume_fluxes, temporaries,
643 : primitive_vars,
644 : typename DerivedCorrection::dg_package_data_volume_tags{});
645 :
646 : detail::apply_boundary_conditions_on_all_external_faces<
647 : EvolutionSystem, Dim>(
648 : make_not_null(&box),
649 : dynamic_cast<const DerivedCorrection&>(boundary_correction),
650 : temporaries, volume_fluxes, partial_derivs, primitive_vars);
651 : }
652 : });
653 :
654 : if constexpr (LocalTimeStepping) {
655 : db::mutate_apply<ChangeStepSize<DgStepChoosers>>(make_not_null(&box));
656 : }
657 :
658 : send_data_for_fluxes<ParallelComponent>(make_not_null(&cache),
659 : make_not_null(&box), volume_fluxes);
660 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
661 : }
662 :
663 : template <size_t Dim, typename EvolutionSystem, typename DgStepChoosers,
664 : bool LocalTimeStepping, bool UseNodegroupDgElements>
665 : template <typename ParallelComponent, typename DbTagsList,
666 : typename Metavariables>
667 : void ComputeTimeDerivative<Dim, EvolutionSystem, DgStepChoosers,
668 0 : LocalTimeStepping, UseNodegroupDgElements>::
669 : send_data_for_fluxes(
670 : const gsl::not_null<Parallel::GlobalCache<Metavariables>*> cache,
671 : const gsl::not_null<db::DataBox<DbTagsList>*> box,
672 : [[maybe_unused]] const Variables<db::wrap_tags_in<
673 : ::Tags::Flux, typename EvolutionSystem::flux_variables,
674 : tmpl::size_t<Dim>, Frame::Inertial>>& volume_fluxes) {
675 : using variables_tag = typename EvolutionSystem::variables_tag;
676 :
677 : auto& receiver_proxy =
678 : Parallel::get_parallel_component<ParallelComponent>(*cache);
679 : const auto& element = db::get<domain::Tags::Element<Dim>>(*box);
680 :
681 : const auto& time_step_id = db::get<::Tags::TimeStepId>(*box);
682 : const auto integration_order =
683 : db::get<::Tags::HistoryEvolvedVariables<variables_tag>>(*box)
684 : .integration_order();
685 : const auto& all_mortar_data =
686 : db::get<evolution::dg::Tags::MortarData<Dim>>(*box);
687 : const auto& mortar_meshes = get<evolution::dg::Tags::MortarMesh<Dim>>(*box);
688 : const auto& mortar_info = get<evolution::dg::Tags::MortarInfo<Dim>>(*box);
689 :
690 : std::optional<DirectionMap<Dim, DataVector>>
691 : all_neighbor_data_for_reconstruction = std::nullopt;
692 : int tci_decision = 0;
693 : const Mesh<Dim>& volume_mesh = db::get<domain::Tags::Mesh<Dim>>(*box);
694 : std::optional<Mesh<Dim>> ghost_data_mesh = std::nullopt;
695 : if constexpr (using_subcell_v<Metavariables>) {
696 : if (not all_neighbor_data_for_reconstruction.has_value()) {
697 : all_neighbor_data_for_reconstruction = DirectionMap<Dim, DataVector>{};
698 : }
699 :
700 : evolution::dg::subcell::prepare_neighbor_data<Metavariables>(
701 : make_not_null(&all_neighbor_data_for_reconstruction.value()),
702 : make_not_null(&ghost_data_mesh), box, volume_fluxes);
703 : tci_decision = evolution::dg::subcell::get_tci_decision(*box);
704 : }
705 :
706 : for (const auto& [direction, neighbors] : element.neighbors()) {
707 : std::optional<DataVector> ghost_and_subcell_data = std::nullopt;
708 : if constexpr (using_subcell_v<Metavariables>) {
709 : ASSERT(all_neighbor_data_for_reconstruction.has_value(),
710 : "Trying to do DG-subcell but the ghost and subcell data for the "
711 : "neighbor has not been set.");
712 : ghost_and_subcell_data =
713 : std::move(all_neighbor_data_for_reconstruction.value()[direction]);
714 : }
715 :
716 : const size_t total_neighbors = neighbors.size();
717 : // If there are multiple non-conforming neighbors, we only create a single
718 : // mortar labeled by the host ElementId. This is done because the data
719 : // from all neighbors will be combined onto a single mortar as it makes no
720 : // sense to have multiple mortars between non-conforming Elements.
721 : const bool has_multiple_non_conforming_neighbors =
722 : total_neighbors > 1 and not neighbors.are_conforming();
723 : size_t neighbor_count = 1;
724 : for (const auto& neighbor : neighbors) {
725 : const auto& orientation = neighbors.orientation(neighbor);
726 : const auto direction_from_neighbor = orientation(direction.opposite());
727 : const DirectionalId<Dim> mortar_id{
728 : direction,
729 : has_multiple_non_conforming_neighbors ? element.id() : neighbor};
730 : const Mesh<Dim - 1>& mortar_mesh = mortar_meshes.at(mortar_id);
731 : auto volume_mesh_for_neighbor = volume_mesh;
732 : auto mortar_mesh_for_neighbor = mortar_mesh;
733 : DataVector neighbor_boundary_data_on_mortar{};
734 : std::optional<InterpolatedBoundaryData<Dim>> interpolated_boundary_data{
735 : std::nullopt};
736 :
737 : switch (mortar_info.at(mortar_id).interface_data_policy()) {
738 : case InterfaceDataPolicy::CopyProject:
739 : [[fallthrough]];
740 : case InterfaceDataPolicy::NonconformingNeighborInterpolates:
741 : neighbor_boundary_data_on_mortar =
742 : *all_mortar_data.at(mortar_id).local().mortar_data.value();
743 : break;
744 : case InterfaceDataPolicy::OrientCopyProject: {
745 : volume_mesh_for_neighbor = orientation(volume_mesh);
746 : mortar_mesh_for_neighbor = orient_mesh_on_slice(
747 : mortar_mesh, direction.dimension(), orientation);
748 : const auto& slice_extents = mortar_mesh.extents();
749 : neighbor_boundary_data_on_mortar = orient_variables_on_slice(
750 : all_mortar_data.at(mortar_id).local().mortar_data.value(),
751 : slice_extents, direction.dimension(), orientation);
752 : break;
753 : }
754 : case InterfaceDataPolicy::NonconformingSelfInterpolates: {
755 : if constexpr (Dim > 1) {
756 : neighbor_boundary_data_on_mortar =
757 : *all_mortar_data.at(mortar_id).local().mortar_data.value();
758 : const auto& interpolator =
759 : mortar_info.at(mortar_id).interpolator().value();
760 : interpolated_boundary_data = InterpolatedBoundaryData<Dim>{
761 : {.data = interpolator.interpolate_to_neighbor(
762 : neighbor_boundary_data_on_mortar),
763 : .target_mesh = interpolator.neighbor_mortar_mesh(),
764 : .offsets = interpolator.interpolated_neighbor_data_offsets()}};
765 : } else {
766 : ERROR("Cannot have non-conforming neighbors in 1D");
767 : }
768 : break;
769 : }
770 : default:
771 : ERROR("InterfaceDataPolicy "
772 : << mortar_info.at(mortar_id).interface_data_policy()
773 : << " is not handled yet, id = " << mortar_id);
774 : }
775 :
776 : const TimeStepId& next_time_step_id =
777 : db::get<::Tags::Next<::Tags::TimeStepId>>(*box);
778 :
779 : using SendData = evolution::dg::BoundaryData<Dim>;
780 : SendData data{};
781 :
782 : if (neighbor_count == total_neighbors) {
783 : data = SendData{volume_mesh_for_neighbor,
784 : ghost_data_mesh,
785 : mortar_mesh_for_neighbor,
786 : std::move(ghost_and_subcell_data),
787 : {std::move(neighbor_boundary_data_on_mortar)},
788 : next_time_step_id,
789 : tci_decision,
790 : integration_order,
791 : interpolated_boundary_data};
792 : } else {
793 : data = SendData{volume_mesh_for_neighbor,
794 : ghost_data_mesh,
795 : mortar_mesh_for_neighbor,
796 : ghost_and_subcell_data,
797 : {std::move(neighbor_boundary_data_on_mortar)},
798 : next_time_step_id,
799 : tci_decision,
800 : integration_order,
801 : interpolated_boundary_data};
802 : }
803 :
804 : // Send mortar data (the `std::tuple` named `data`) to neighbor
805 : if constexpr (Parallel::is_dg_element_collection_v<ParallelComponent>) {
806 : Parallel::local_synchronous_action<
807 : Parallel::Actions::SendDataToElement>(
808 : receiver_proxy, cache,
809 : evolution::dg::Tags::BoundaryCorrectionAndGhostCellsInbox<
810 : Dim, UseNodegroupDgElements>{},
811 : neighbor, time_step_id,
812 : std::make_pair(DirectionalId{direction_from_neighbor, element.id()},
813 : std::move(data)));
814 : } else {
815 : Parallel::receive_data<
816 : evolution::dg::Tags::BoundaryCorrectionAndGhostCellsInbox<
817 : Dim, UseNodegroupDgElements>>(
818 : receiver_proxy[neighbor], time_step_id,
819 : std::make_pair(DirectionalId{direction_from_neighbor, element.id()},
820 : std::move(data)));
821 : }
822 : ++neighbor_count;
823 : }
824 : }
825 :
826 : // We treat this as a set, but use a map because we don't have a
827 : // non-allocating set type.
828 : DirectionMap<Dim, bool> mortar_history_directions{};
829 : for (const auto& [mortar, info] : mortar_info) {
830 : if (info.time_stepping_policy() == TimeSteppingPolicy::Conservative) {
831 : mortar_history_directions.emplace(mortar.direction(), true);
832 : }
833 : }
834 :
835 : if (not mortar_history_directions.empty()) {
836 : // Need volume Jacobian for any face whose normal direction uses Gauss
837 : // points. This means mixed-quadrature non-hypercube elements (e.g.
838 : // full_cylinder) where some directions have collocated face points and
839 : // others do not.
840 : const bool any_direction_uses_gauss =
841 : alg::any_of(volume_mesh.quadrature(), [](const Spectral::Quadrature q) {
842 : return q == Spectral::Quadrature::Gauss;
843 : });
844 :
845 : const Scalar<DataVector> volume_det_inv_jacobian{};
846 : if (any_direction_uses_gauss) {
847 : // NOLINTNEXTLINE
848 : const_cast<DataVector&>(get(volume_det_inv_jacobian))
849 : .set_data_ref(make_not_null(&const_cast<DataVector&>( // NOLINT
850 : get(db::get<domain::Tags::DetInvJacobian<
851 : Frame::ElementLogical, Frame::Inertial>>(*box)))));
852 : }
853 :
854 : // Add face normal and Jacobian determinants to the local mortar data. We
855 : // only need the Jacobians for directions using Gauss points. Then copy
856 : // over into the boundary history, since that's what the LTS steppers use.
857 : //
858 : // The boundary history coupling computation (which computes the _lifted_
859 : // boundary correction) returns a Variables<dt<EvolvedVars>> instead of
860 : // using the `NormalDotNumericalFlux` prefix tag. This is because the
861 : // returned quantity is more a `dt` quantity than a
862 : // `NormalDotNormalDotFlux` since it's been lifted to the volume.
863 : db::mutate<evolution::dg::Tags::MortarData<Dim>,
864 : evolution::dg::Tags::MortarDataHistory<Dim>>(
865 : [&element, integration_order, &mortar_history_directions, &mortar_info,
866 : &time_step_id, any_direction_uses_gauss, &volume_det_inv_jacobian,
867 : &volume_mesh](
868 : const gsl::not_null<
869 : DirectionalIdMap<Dim, evolution::dg::MortarDataHolder<Dim>>*>
870 : mortar_data,
871 : const gsl::not_null<DirectionalIdMap<
872 : Dim, TimeSteppers::BoundaryHistory<
873 : evolution::dg::MortarData<Dim>,
874 : evolution::dg::MortarData<Dim>, DataVector>>*>
875 : boundary_data_history,
876 : const DirectionMap<Dim,
877 : std::optional<Variables<tmpl::list<
878 : evolution::dg::Tags::MagnitudeOfNormal,
879 : evolution::dg::Tags::NormalCovector<Dim>>>>>&
880 : normal_covector_and_magnitude) {
881 : Scalar<DataVector> volume_det_jacobian{};
882 : Scalar<DataVector> face_det_jacobian{};
883 : if (any_direction_uses_gauss) {
884 : get(volume_det_jacobian) = 1.0 / get(volume_det_inv_jacobian);
885 : }
886 : for (const auto& [direction, neighbors_in_direction] :
887 : element.neighbors()) {
888 : if (not mortar_history_directions.contains(direction)) {
889 : continue;
890 : }
891 : const size_t total_neighbors = neighbors_in_direction.size();
892 : // If there are multiple non-conforming neighbors, we only create a
893 : // single mortar labeled by the host ElementId. This is done
894 : // because the data from all neighbors will be combined onto a
895 : // single mortar as it makes no sense to have multiple mortars
896 : // between non-conforming Elements.
897 : const bool has_multiple_non_conforming_neighbors =
898 : total_neighbors > 1 and
899 : not neighbors_in_direction.are_conforming();
900 : // We can perform projections once for all neighbors in the
901 : // direction because we care about the _face_ mesh, not the mortar
902 : // mesh.
903 : ASSERT(normal_covector_and_magnitude.at(direction).has_value(),
904 : "The normal covector and magnitude have not been computed.");
905 : const Scalar<DataVector>& face_normal_magnitude =
906 : get<evolution::dg::Tags::MagnitudeOfNormal>(
907 : *normal_covector_and_magnitude.at(direction));
908 : if (volume_mesh.quadrature(direction.dimension()) ==
909 : Spectral::Quadrature::Gauss) {
910 : const Matrix identity{};
911 : auto interpolation_matrices =
912 : make_array<Dim>(std::cref(identity));
913 : const std::pair<Matrix, Matrix>& matrices =
914 : Spectral::boundary_interpolation_matrices(
915 : volume_mesh.slice_through(direction.dimension()));
916 : gsl::at(interpolation_matrices, direction.dimension()) =
917 : direction.side() == Side::Upper ? matrices.second
918 : : matrices.first;
919 : if (get(face_det_jacobian).size() !=
920 : get(face_normal_magnitude).size()) {
921 : get(face_det_jacobian) =
922 : DataVector{get(face_normal_magnitude).size()};
923 : }
924 : apply_matrices(make_not_null(&get(face_det_jacobian)),
925 : interpolation_matrices, get(volume_det_jacobian),
926 : volume_mesh.extents());
927 : }
928 :
929 : for (const auto& neighbor : neighbors_in_direction) {
930 : const DirectionalId<Dim> mortar_id{
931 : direction, has_multiple_non_conforming_neighbors
932 : ? element.id()
933 : : neighbor};
934 : if (mortar_info.at(mortar_id).time_stepping_policy() !=
935 : TimeSteppingPolicy::Conservative) {
936 : continue;
937 : }
938 : auto& local_mortar_data = mortar_data->at(mortar_id).local();
939 : local_mortar_data.face_normal_magnitude = face_normal_magnitude;
940 : if (volume_mesh.quadrature(direction.dimension()) ==
941 : Spectral::Quadrature::Gauss) {
942 : local_mortar_data.volume_mesh = volume_mesh;
943 : local_mortar_data.volume_det_inv_jacobian =
944 : volume_det_inv_jacobian;
945 : local_mortar_data.face_det_jacobian = face_det_jacobian;
946 : }
947 : ASSERT(boundary_data_history->count(mortar_id) != 0,
948 : "Could not insert the mortar data for "
949 : << mortar_id
950 : << " because the unordered map has not been "
951 : "initialized "
952 : "to have the mortar id.");
953 : boundary_data_history->at(mortar_id).local().insert(
954 : time_step_id, integration_order,
955 : std::move(mortar_data->at(mortar_id).local()));
956 : mortar_data->at(mortar_id) = MortarDataHolder<Dim>{};
957 : }
958 : }
959 : },
960 : box,
961 : db::get<evolution::dg::Tags::NormalCovectorAndMagnitude<Dim>>(*box));
962 : }
963 : }
964 : } // namespace evolution::dg::Actions
|
