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 : // Systems may declare an `auxiliary_variables` type alias whose first
420 : // derivatives are needed by the volume terms. The detect-or-default
421 : // metafunction yields an empty list for systems without it.
422 : using auxiliary_variables =
423 : detail::get_auxiliary_variables_or_default_t<EvolutionSystem,
424 : tmpl::list<>>;
425 :
426 : const Mesh<Dim>& mesh = db::get<::domain::Tags::Mesh<Dim>>(box);
427 : const Element<Dim>& element = db::get<domain::Tags::Element<Dim>>(box);
428 : const ::dg::Formulation dg_formulation =
429 : db::get<::dg::Tags::Formulation>(box);
430 : ASSERT(alg::all_of(mesh.basis(),
431 : [&mesh](const Spectral::Basis current_basis) {
432 : return current_basis == mesh.basis(0);
433 : }) or
434 : element.topologies() != domain::topologies::hypercube<Dim>,
435 : "An isotropic basis must be used in the evolution code. While "
436 : "theoretically this restriction could be lifted, the simplification "
437 : "it offers are quite substantial. Relaxing this assumption is likely "
438 : "to require quite a bit of careful code refactoring and debugging.");
439 : ASSERT(alg::all_of(mesh.quadrature(),
440 : [&mesh](const Spectral::Quadrature current_quadrature) {
441 : return current_quadrature == mesh.quadrature(0);
442 : }) or
443 : element.topologies() != domain::topologies::hypercube<Dim>,
444 : "An isotropic quadrature must be used in the evolution code. While "
445 : "theoretically this restriction could be lifted, the simplification "
446 : "it offers are quite substantial. Relaxing this assumption is likely "
447 : "to require quite a bit of careful code refactoring and debugging.");
448 :
449 : const auto& boundary_correction =
450 : db::get<evolution::Tags::BoundaryCorrection>(box);
451 : using derived_boundary_corrections =
452 : tmpl::at<typename Metavariables::factory_creation::factory_classes,
453 : evolution::BoundaryCorrection>;
454 :
455 : // To avoid a second allocation in internal_mortar_data, we allocate the
456 : // variables needed to construct the fields on the faces here along with
457 : // everything else. This requires us to know all the tags necessary to apply
458 : // boundary corrections. However, since we pick boundary corrections at
459 : // runtime, we just gather all possible tags from all possible boundary
460 : // corrections and lump them into the allocation. This may result in a
461 : // larger-than-necessary allocation, but it won't be that much larger.
462 : using all_dg_package_temporary_tags =
463 : tmpl::transform<derived_boundary_corrections,
464 : detail::get_dg_package_temporary_tags<tmpl::_1>>;
465 : using all_primitive_tags_for_face =
466 : tmpl::transform<derived_boundary_corrections,
467 : detail::get_primitive_tags_for_face<
468 : tmpl::pin<EvolutionSystem>, tmpl::_1>>;
469 : using fluxes_tags = db::wrap_tags_in<::Tags::Flux, flux_variables,
470 : tmpl::size_t<Dim>, Frame::Inertial>;
471 : using dg_package_data_projected_tags =
472 : tmpl::list<typename variables_tag::tags_list, fluxes_tags,
473 : all_dg_package_temporary_tags, all_primitive_tags_for_face>;
474 : using all_face_temporary_tags =
475 : tmpl::remove_duplicates<tmpl::flatten<tmpl::push_back<
476 : tmpl::list<dg_package_data_projected_tags,
477 : detail::inverse_spatial_metric_tag<EvolutionSystem>>,
478 : detail::OneOverNormalVectorMagnitude, detail::NormalVector<Dim>>>>;
479 : // To avoid additional allocations in internal_mortar_data, we provide a
480 : // buffer used to compute the packaged data before it has to be projected to
481 : // the mortar. We get all mortar tags for similar reasons as described above
482 : using all_mortar_tags = tmpl::remove_duplicates<tmpl::flatten<
483 : tmpl::transform<derived_boundary_corrections,
484 : detail::get_dg_package_field_tags<tmpl::_1>>>>;
485 :
486 : // We also don't use the number of volume mesh grid points. We instead use the
487 : // max number of grid points from each face. That way, our allocation will be
488 : // large enough to hold any face and we can reuse the allocation for each face
489 : // without having to resize it.
490 : size_t num_face_temporary_grid_points = 0;
491 : {
492 : for (const auto& [direction, neighbors_in_direction] :
493 : element.neighbors()) {
494 : (void)neighbors_in_direction;
495 : const auto face_mesh = mesh.slice_away(direction.dimension());
496 : num_face_temporary_grid_points = std::max(
497 : num_face_temporary_grid_points, face_mesh.number_of_grid_points());
498 : }
499 : }
500 :
501 : // Allocate the Variables classes needed for the time derivative
502 : // computation.
503 : //
504 : // This is factored out so that we will be able to do ADER-DG/CG where a
505 : // spacetime polynomial is constructed by solving implicit equations in time
506 : // using a Picard iteration. A high-order initial guess is needed to
507 : // efficiently construct the ADER spacetime solution. This initial guess is
508 : // obtained using continuous RK methods, and so we will want to reuse
509 : // buffers. Thus, the volume_terms function returns by reference rather than
510 : // by value.
511 : using VarsTemporaries =
512 : Variables<typename compute_volume_time_derivative_terms::temporary_tags>;
513 : using VarsFluxes =
514 : Variables<db::wrap_tags_in<::Tags::Flux, flux_variables,
515 : tmpl::size_t<Dim>, Frame::Inertial>>;
516 : using VarsPartialDerivatives =
517 : Variables<db::wrap_tags_in<::Tags::deriv, partial_derivative_tags,
518 : tmpl::size_t<Dim>, Frame::Inertial>>;
519 : using VarsDivFluxes = Variables<db::wrap_tags_in<
520 : ::Tags::div, db::wrap_tags_in<::Tags::Flux, flux_variables,
521 : tmpl::size_t<Dim>, Frame::Inertial>>>;
522 : using VarsFaceTemporaries = Variables<all_face_temporary_tags>;
523 : using DgPackagedDataVarsOnFace = Variables<all_mortar_tags>;
524 : const size_t number_of_grid_points = mesh.number_of_grid_points();
525 : const size_t buffer_size =
526 : (VarsTemporaries::number_of_independent_components +
527 : VarsFluxes::number_of_independent_components +
528 : VarsPartialDerivatives::number_of_independent_components +
529 : VarsDivFluxes::number_of_independent_components) *
530 : number_of_grid_points +
531 : // Different number of grid points. See explanation above where
532 : // num_face_temporary_grid_points is defined
533 : (VarsFaceTemporaries::number_of_independent_components +
534 : DgPackagedDataVarsOnFace::number_of_independent_components) *
535 : num_face_temporary_grid_points;
536 : auto buffer = cpp20::make_unique_for_overwrite<double[]>(buffer_size);
537 : #ifdef SPECTRE_NAN_INIT
538 : std::fill(&buffer[0], &buffer[buffer_size],
539 : std::numeric_limits<double>::signaling_NaN());
540 : #endif
541 : VarsTemporaries temporaries{
542 : &buffer[0], VarsTemporaries::number_of_independent_components *
543 : number_of_grid_points};
544 : VarsFluxes volume_fluxes{
545 : &buffer[VarsTemporaries::number_of_independent_components *
546 : number_of_grid_points],
547 : VarsFluxes::number_of_independent_components * number_of_grid_points};
548 : VarsPartialDerivatives partial_derivs{
549 : &buffer[(VarsTemporaries::number_of_independent_components +
550 : VarsFluxes::number_of_independent_components) *
551 : number_of_grid_points],
552 : VarsPartialDerivatives::number_of_independent_components *
553 : number_of_grid_points};
554 : VarsDivFluxes div_fluxes{
555 : &buffer[(VarsTemporaries::number_of_independent_components +
556 : VarsFluxes::number_of_independent_components +
557 : VarsPartialDerivatives::number_of_independent_components) *
558 : number_of_grid_points],
559 : VarsDivFluxes::number_of_independent_components * number_of_grid_points};
560 : // Lighter weight data structure than a Variables to avoid passing even more
561 : // templates to internal_mortar_data.
562 : gsl::span<double> face_temporaries = gsl::make_span<double>(
563 : &buffer[(VarsTemporaries::number_of_independent_components +
564 : VarsFluxes::number_of_independent_components +
565 : VarsPartialDerivatives::number_of_independent_components +
566 : VarsDivFluxes::number_of_independent_components) *
567 : number_of_grid_points],
568 : // Different number of grid points. See explanation above where
569 : // num_face_temporary_grid_points is defined
570 : VarsFaceTemporaries::number_of_independent_components *
571 : num_face_temporary_grid_points);
572 : gsl::span<double> packaged_data_buffer = gsl::make_span<double>(
573 : &buffer[(VarsTemporaries::number_of_independent_components +
574 : VarsFluxes::number_of_independent_components +
575 : VarsPartialDerivatives::number_of_independent_components +
576 : VarsDivFluxes::number_of_independent_components) *
577 : number_of_grid_points +
578 : VarsFaceTemporaries::number_of_independent_components *
579 : num_face_temporary_grid_points],
580 : // Different number of grid points. See explanation above where
581 : // num_face_temporary_grid_points is defined
582 : DgPackagedDataVarsOnFace::number_of_independent_components *
583 : num_face_temporary_grid_points);
584 :
585 : const Scalar<DataVector>* det_inverse_jacobian = nullptr;
586 : if constexpr (tmpl::size<flux_variables>::value != 0) {
587 : if (dg_formulation == ::dg::Formulation::WeakInertial) {
588 : det_inverse_jacobian = &db::get<
589 : domain::Tags::DetInvJacobian<Frame::ElementLogical, Frame::Inertial>>(
590 : box);
591 : }
592 : }
593 : if constexpr (tmpl::size<auxiliary_variables>::value != 0) {
594 : static_assert(
595 : tmpl::size<tmpl::list_difference<
596 : partial_derivative_tags,
597 : tmpl::append<typename variables_tag::tags_list,
598 : auxiliary_variables>>>::value == 0,
599 : "Every gradient variable must be an evolved variable (in "
600 : "variables_tag) or an auxiliary variable (in auxiliary_variables); "
601 : "otherwise it is not populated in the combined differentiation "
602 : "source.");
603 : Variables<detail::evolved_and_auxiliary_vars_tags<EvolutionSystem>>
604 : evolved_and_auxiliary_vars{mesh.number_of_grid_points()};
605 : evolved_and_auxiliary_vars.assign_subset(db::get<variables_tag>(box));
606 : evolved_and_auxiliary_vars.assign_subset(
607 : db::get<::Tags::Variables<auxiliary_variables>>(box));
608 : db::mutate_apply<
609 : tmpl::list<dt_variables_tag>,
610 : typename compute_volume_time_derivative_terms::argument_tags>(
611 : [&dg_formulation, &div_fluxes, &det_inverse_jacobian,
612 : &div_mesh_velocity = db::get<::domain::Tags::DivMeshVelocity>(box),
613 : &evolved_and_auxiliary_vars,
614 : &inertial_coordinates =
615 : db::get<domain::Tags::Coordinates<Dim, Frame::Inertial>>(box),
616 : &logical_to_inertial_inv_jacobian =
617 : db::get<::domain::Tags::InverseJacobian<Dim, Frame::ElementLogical,
618 : Frame::Inertial>>(box),
619 : &mesh,
620 : &mesh_velocity = db::get<::domain::Tags::MeshVelocity<Dim>>(box),
621 : &partial_derivs, &temporaries,
622 : &volume_fluxes](const gsl::not_null<Variables<db::wrap_tags_in<
623 : ::Tags::dt, typename variables_tag::tags_list>>*>
624 : dt_vars_ptr,
625 : const auto&... time_derivative_args) {
626 : detail::volume_terms<compute_volume_time_derivative_terms>(
627 : dt_vars_ptr, make_not_null(&volume_fluxes),
628 : make_not_null(&partial_derivs), make_not_null(&temporaries),
629 : make_not_null(&div_fluxes), evolved_and_auxiliary_vars,
630 : dg_formulation, mesh, inertial_coordinates,
631 : logical_to_inertial_inv_jacobian, det_inverse_jacobian,
632 : mesh_velocity, div_mesh_velocity, time_derivative_args...);
633 : },
634 : make_not_null(&box));
635 : } else {
636 : db::mutate_apply<
637 : tmpl::list<dt_variables_tag>,
638 : typename compute_volume_time_derivative_terms::argument_tags>(
639 : [&dg_formulation, &div_fluxes, &det_inverse_jacobian,
640 : &div_mesh_velocity = db::get<::domain::Tags::DivMeshVelocity>(box),
641 : &evolved_variables = db::get<variables_tag>(box),
642 : &inertial_coordinates =
643 : db::get<domain::Tags::Coordinates<Dim, Frame::Inertial>>(box),
644 : &logical_to_inertial_inv_jacobian =
645 : db::get<::domain::Tags::InverseJacobian<Dim, Frame::ElementLogical,
646 : Frame::Inertial>>(box),
647 : &mesh,
648 : &mesh_velocity = db::get<::domain::Tags::MeshVelocity<Dim>>(box),
649 : &partial_derivs, &temporaries,
650 : &volume_fluxes](const gsl::not_null<Variables<db::wrap_tags_in<
651 : ::Tags::dt, typename variables_tag::tags_list>>*>
652 : dt_vars_ptr,
653 : const auto&... time_derivative_args) {
654 : detail::volume_terms<compute_volume_time_derivative_terms>(
655 : dt_vars_ptr, make_not_null(&volume_fluxes),
656 : make_not_null(&partial_derivs), make_not_null(&temporaries),
657 : make_not_null(&div_fluxes), evolved_variables, dg_formulation,
658 : mesh, inertial_coordinates, logical_to_inertial_inv_jacobian,
659 : det_inverse_jacobian, mesh_velocity, div_mesh_velocity,
660 : time_derivative_args...);
661 : },
662 : make_not_null(&box));
663 : }
664 :
665 : const Variables<detail::get_primitive_vars_tags_from_system<EvolutionSystem>>*
666 : primitive_vars{nullptr};
667 : if constexpr (EvolutionSystem::has_primitive_and_conservative_vars) {
668 : primitive_vars =
669 : &db::get<typename EvolutionSystem::primitive_variables_tag>(box);
670 : }
671 :
672 : static_assert(
673 : tmpl::all<derived_boundary_corrections, std::is_final<tmpl::_1>>::value,
674 : "All createable classes for boundary corrections must be marked "
675 : "final.");
676 : tmpl::for_each<derived_boundary_corrections>(
677 : [&boundary_correction, &box, &partial_derivs, &primitive_vars,
678 : &temporaries, &volume_fluxes, &packaged_data_buffer,
679 : &face_temporaries](auto derived_correction_v) {
680 : using DerivedCorrection =
681 : tmpl::type_from<decltype(derived_correction_v)>;
682 : if (typeid(boundary_correction) == typeid(DerivedCorrection)) {
683 : // Compute internal boundary quantities on the mortar for sides
684 : // of the element that have neighbors, i.e. they are not an
685 : // external side.
686 : // Note: this call mutates:
687 : // - evolution::dg::Tags::NormalCovectorAndMagnitude<Dim>,
688 : // - evolution::dg::Tags::MortarData<Dim>
689 : detail::internal_mortar_data<EvolutionSystem, Dim>(
690 : make_not_null(&box), make_not_null(&face_temporaries),
691 : make_not_null(&packaged_data_buffer),
692 : dynamic_cast<const DerivedCorrection&>(boundary_correction),
693 : db::get<variables_tag>(box), volume_fluxes, temporaries,
694 : primitive_vars,
695 : typename DerivedCorrection::dg_package_data_volume_tags{});
696 :
697 : detail::apply_boundary_conditions_on_all_external_faces<
698 : EvolutionSystem, Dim>(
699 : make_not_null(&box),
700 : dynamic_cast<const DerivedCorrection&>(boundary_correction),
701 : temporaries, volume_fluxes, partial_derivs, primitive_vars);
702 : }
703 : });
704 :
705 : if constexpr (LocalTimeStepping) {
706 : db::mutate_apply<ChangeStepSize<DgStepChoosers>>(make_not_null(&box));
707 : }
708 :
709 : send_data_for_fluxes<ParallelComponent>(make_not_null(&cache),
710 : make_not_null(&box), volume_fluxes);
711 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
712 : }
713 :
714 : template <size_t Dim, typename EvolutionSystem, typename DgStepChoosers,
715 : bool LocalTimeStepping, bool UseNodegroupDgElements>
716 : template <typename ParallelComponent, typename DbTagsList,
717 : typename Metavariables>
718 : void ComputeTimeDerivative<Dim, EvolutionSystem, DgStepChoosers,
719 0 : LocalTimeStepping, UseNodegroupDgElements>::
720 : send_data_for_fluxes(
721 : const gsl::not_null<Parallel::GlobalCache<Metavariables>*> cache,
722 : const gsl::not_null<db::DataBox<DbTagsList>*> box,
723 : [[maybe_unused]] const Variables<db::wrap_tags_in<
724 : ::Tags::Flux, typename EvolutionSystem::flux_variables,
725 : tmpl::size_t<Dim>, Frame::Inertial>>& volume_fluxes) {
726 : using variables_tag = typename EvolutionSystem::variables_tag;
727 :
728 : auto& receiver_proxy =
729 : Parallel::get_parallel_component<ParallelComponent>(*cache);
730 : const auto& element = db::get<domain::Tags::Element<Dim>>(*box);
731 :
732 : const auto& time_step_id = db::get<::Tags::TimeStepId>(*box);
733 : const auto integration_order =
734 : db::get<::Tags::HistoryEvolvedVariables<variables_tag>>(*box)
735 : .integration_order();
736 : const auto& all_mortar_data =
737 : db::get<evolution::dg::Tags::MortarData<Dim>>(*box);
738 : const auto& mortar_meshes = get<evolution::dg::Tags::MortarMesh<Dim>>(*box);
739 : const auto& mortar_info = get<evolution::dg::Tags::MortarInfo<Dim>>(*box);
740 :
741 : std::optional<DirectionMap<Dim, DataVector>>
742 : all_neighbor_data_for_reconstruction = std::nullopt;
743 : int tci_decision = 0;
744 : const Mesh<Dim>& volume_mesh = db::get<domain::Tags::Mesh<Dim>>(*box);
745 : std::optional<Mesh<Dim>> ghost_data_mesh = std::nullopt;
746 : if constexpr (using_subcell_v<Metavariables>) {
747 : if (not all_neighbor_data_for_reconstruction.has_value()) {
748 : all_neighbor_data_for_reconstruction = DirectionMap<Dim, DataVector>{};
749 : }
750 :
751 : evolution::dg::subcell::prepare_neighbor_data<Metavariables>(
752 : make_not_null(&all_neighbor_data_for_reconstruction.value()),
753 : make_not_null(&ghost_data_mesh), box, volume_fluxes);
754 : tci_decision = evolution::dg::subcell::get_tci_decision(*box);
755 : }
756 :
757 : for (const auto& [direction, neighbors] : element.neighbors()) {
758 : std::optional<DataVector> ghost_and_subcell_data = std::nullopt;
759 : if constexpr (using_subcell_v<Metavariables>) {
760 : ASSERT(all_neighbor_data_for_reconstruction.has_value(),
761 : "Trying to do DG-subcell but the ghost and subcell data for the "
762 : "neighbor has not been set.");
763 : ghost_and_subcell_data =
764 : std::move(all_neighbor_data_for_reconstruction.value()[direction]);
765 : }
766 :
767 : const size_t total_neighbors = neighbors.size();
768 : // If there are multiple non-conforming neighbors, we only create a single
769 : // mortar labeled by the host ElementId. This is done because the data
770 : // from all neighbors will be combined onto a single mortar as it makes no
771 : // sense to have multiple mortars between non-conforming Elements.
772 : const bool has_multiple_non_conforming_neighbors =
773 : total_neighbors > 1 and not neighbors.are_conforming();
774 : size_t neighbor_count = 1;
775 : for (const auto& neighbor : neighbors) {
776 : const auto& orientation = neighbors.orientation(neighbor);
777 : const auto direction_from_neighbor = orientation(direction.opposite());
778 : const DirectionalId<Dim> mortar_id{
779 : direction,
780 : has_multiple_non_conforming_neighbors ? element.id() : neighbor};
781 : const Mesh<Dim - 1>& mortar_mesh = mortar_meshes.at(mortar_id);
782 : auto volume_mesh_for_neighbor = volume_mesh;
783 : auto mortar_mesh_for_neighbor = mortar_mesh;
784 : DataVector neighbor_boundary_data_on_mortar{};
785 : std::optional<InterpolatedBoundaryData<Dim>> interpolated_boundary_data{
786 : std::nullopt};
787 :
788 : switch (mortar_info.at(mortar_id).interface_data_policy()) {
789 : case InterfaceDataPolicy::CopyProject:
790 : [[fallthrough]];
791 : case InterfaceDataPolicy::NonconformingNeighborInterpolates:
792 : neighbor_boundary_data_on_mortar =
793 : *all_mortar_data.at(mortar_id).local().mortar_data.value();
794 : break;
795 : case InterfaceDataPolicy::OrientCopyProject: {
796 : volume_mesh_for_neighbor = orientation(volume_mesh);
797 : mortar_mesh_for_neighbor = orient_mesh_on_slice(
798 : mortar_mesh, direction.dimension(), orientation);
799 : const auto& slice_extents = mortar_mesh.extents();
800 : neighbor_boundary_data_on_mortar = orient_variables_on_slice(
801 : all_mortar_data.at(mortar_id).local().mortar_data.value(),
802 : slice_extents, direction.dimension(), orientation);
803 : break;
804 : }
805 : case InterfaceDataPolicy::NonconformingSelfInterpolates: {
806 : if constexpr (Dim > 1) {
807 : neighbor_boundary_data_on_mortar =
808 : *all_mortar_data.at(mortar_id).local().mortar_data.value();
809 : const auto& interpolator =
810 : mortar_info.at(mortar_id).interpolator().value();
811 : interpolated_boundary_data = InterpolatedBoundaryData<Dim>{
812 : {.data = interpolator.interpolate_to_neighbor(
813 : neighbor_boundary_data_on_mortar),
814 : .target_mesh = interpolator.neighbor_mortar_mesh(),
815 : .offsets = interpolator.interpolated_neighbor_data_offsets()}};
816 : } else {
817 : ERROR("Cannot have non-conforming neighbors in 1D");
818 : }
819 : break;
820 : }
821 : default:
822 : ERROR("InterfaceDataPolicy "
823 : << mortar_info.at(mortar_id).interface_data_policy()
824 : << " is not handled yet, id = " << mortar_id);
825 : }
826 :
827 : const TimeStepId& next_time_step_id =
828 : db::get<::Tags::Next<::Tags::TimeStepId>>(*box);
829 :
830 : using SendData = evolution::dg::BoundaryData<Dim>;
831 : SendData data{};
832 :
833 : if (neighbor_count == total_neighbors) {
834 : data = SendData{volume_mesh_for_neighbor,
835 : ghost_data_mesh,
836 : mortar_mesh_for_neighbor,
837 : std::move(ghost_and_subcell_data),
838 : {std::move(neighbor_boundary_data_on_mortar)},
839 : next_time_step_id,
840 : tci_decision,
841 : integration_order,
842 : interpolated_boundary_data};
843 : } else {
844 : data = SendData{volume_mesh_for_neighbor,
845 : ghost_data_mesh,
846 : mortar_mesh_for_neighbor,
847 : ghost_and_subcell_data,
848 : {std::move(neighbor_boundary_data_on_mortar)},
849 : next_time_step_id,
850 : tci_decision,
851 : integration_order,
852 : interpolated_boundary_data};
853 : }
854 :
855 : // Send mortar data (the `std::tuple` named `data`) to neighbor
856 : if constexpr (Parallel::is_dg_element_collection_v<ParallelComponent>) {
857 : Parallel::local_synchronous_action<
858 : Parallel::Actions::SendDataToElement>(
859 : receiver_proxy, cache,
860 : evolution::dg::Tags::BoundaryCorrectionAndGhostCellsInbox<
861 : Dim, UseNodegroupDgElements>{},
862 : neighbor, time_step_id,
863 : std::make_pair(DirectionalId{direction_from_neighbor, element.id()},
864 : std::move(data)));
865 : } else {
866 : Parallel::receive_data<
867 : evolution::dg::Tags::BoundaryCorrectionAndGhostCellsInbox<
868 : Dim, UseNodegroupDgElements>>(
869 : receiver_proxy[neighbor], time_step_id,
870 : std::make_pair(DirectionalId{direction_from_neighbor, element.id()},
871 : std::move(data)));
872 : }
873 : ++neighbor_count;
874 : }
875 : }
876 :
877 : // We treat this as a set, but use a map because we don't have a
878 : // non-allocating set type.
879 : DirectionMap<Dim, bool> mortar_history_directions{};
880 : for (const auto& [mortar, info] : mortar_info) {
881 : if (info.time_stepping_policy() == TimeSteppingPolicy::Conservative) {
882 : mortar_history_directions.emplace(mortar.direction(), true);
883 : }
884 : }
885 :
886 : if (not mortar_history_directions.empty()) {
887 : // Need volume Jacobian for any face whose normal direction uses Gauss
888 : // points. This means mixed-quadrature non-hypercube elements (e.g.
889 : // full_cylinder) where some directions have collocated face points and
890 : // others do not.
891 : const bool any_direction_uses_gauss =
892 : alg::any_of(volume_mesh.quadrature(), [](const Spectral::Quadrature q) {
893 : return q == Spectral::Quadrature::Gauss;
894 : });
895 :
896 : const Scalar<DataVector> volume_det_inv_jacobian{};
897 : if (any_direction_uses_gauss) {
898 : // NOLINTNEXTLINE
899 : const_cast<DataVector&>(get(volume_det_inv_jacobian))
900 : .set_data_ref(make_not_null(&const_cast<DataVector&>( // NOLINT
901 : get(db::get<domain::Tags::DetInvJacobian<
902 : Frame::ElementLogical, Frame::Inertial>>(*box)))));
903 : }
904 :
905 : // Add face normal and Jacobian determinants to the local mortar data. We
906 : // only need the Jacobians for directions using Gauss points. Then copy
907 : // over into the boundary history, since that's what the LTS steppers use.
908 : //
909 : // The boundary history coupling computation (which computes the _lifted_
910 : // boundary correction) returns a Variables<dt<EvolvedVars>> instead of
911 : // using the `NormalDotNumericalFlux` prefix tag. This is because the
912 : // returned quantity is more a `dt` quantity than a
913 : // `NormalDotNormalDotFlux` since it's been lifted to the volume.
914 : db::mutate<evolution::dg::Tags::MortarData<Dim>,
915 : evolution::dg::Tags::MortarDataHistory<Dim>>(
916 : [&element, integration_order, &mortar_history_directions, &mortar_info,
917 : &time_step_id, any_direction_uses_gauss, &volume_det_inv_jacobian,
918 : &volume_mesh](
919 : const gsl::not_null<
920 : DirectionalIdMap<Dim, evolution::dg::MortarDataHolder<Dim>>*>
921 : mortar_data,
922 : const gsl::not_null<DirectionalIdMap<
923 : Dim, TimeSteppers::BoundaryHistory<
924 : evolution::dg::MortarData<Dim>,
925 : evolution::dg::MortarData<Dim>, DataVector>>*>
926 : boundary_data_history,
927 : const DirectionMap<Dim,
928 : std::optional<Variables<tmpl::list<
929 : evolution::dg::Tags::MagnitudeOfNormal,
930 : evolution::dg::Tags::NormalCovector<Dim>>>>>&
931 : normal_covector_and_magnitude) {
932 : Scalar<DataVector> volume_det_jacobian{};
933 : Scalar<DataVector> face_det_jacobian{};
934 : if (any_direction_uses_gauss) {
935 : get(volume_det_jacobian) = 1.0 / get(volume_det_inv_jacobian);
936 : }
937 : for (const auto& [direction, neighbors_in_direction] :
938 : element.neighbors()) {
939 : if (not mortar_history_directions.contains(direction)) {
940 : continue;
941 : }
942 : const size_t total_neighbors = neighbors_in_direction.size();
943 : // If there are multiple non-conforming neighbors, we only create a
944 : // single mortar labeled by the host ElementId. This is done
945 : // because the data from all neighbors will be combined onto a
946 : // single mortar as it makes no sense to have multiple mortars
947 : // between non-conforming Elements.
948 : const bool has_multiple_non_conforming_neighbors =
949 : total_neighbors > 1 and
950 : not neighbors_in_direction.are_conforming();
951 : // We can perform projections once for all neighbors in the
952 : // direction because we care about the _face_ mesh, not the mortar
953 : // mesh.
954 : ASSERT(normal_covector_and_magnitude.at(direction).has_value(),
955 : "The normal covector and magnitude have not been computed.");
956 : const Scalar<DataVector>& face_normal_magnitude =
957 : get<evolution::dg::Tags::MagnitudeOfNormal>(
958 : *normal_covector_and_magnitude.at(direction));
959 : if (volume_mesh.quadrature(direction.dimension()) ==
960 : Spectral::Quadrature::Gauss) {
961 : const Matrix identity{};
962 : auto interpolation_matrices =
963 : make_array<Dim>(std::cref(identity));
964 : const std::pair<Matrix, Matrix>& matrices =
965 : Spectral::boundary_interpolation_matrices(
966 : volume_mesh.slice_through(direction.dimension()));
967 : gsl::at(interpolation_matrices, direction.dimension()) =
968 : direction.side() == Side::Upper ? matrices.second
969 : : matrices.first;
970 : if (get(face_det_jacobian).size() !=
971 : get(face_normal_magnitude).size()) {
972 : get(face_det_jacobian) =
973 : DataVector{get(face_normal_magnitude).size()};
974 : }
975 : apply_matrices(make_not_null(&get(face_det_jacobian)),
976 : interpolation_matrices, get(volume_det_jacobian),
977 : volume_mesh.extents());
978 : }
979 :
980 : for (const auto& neighbor : neighbors_in_direction) {
981 : const DirectionalId<Dim> mortar_id{
982 : direction, has_multiple_non_conforming_neighbors
983 : ? element.id()
984 : : neighbor};
985 : if (mortar_info.at(mortar_id).time_stepping_policy() !=
986 : TimeSteppingPolicy::Conservative) {
987 : continue;
988 : }
989 : auto& local_mortar_data = mortar_data->at(mortar_id).local();
990 : local_mortar_data.face_normal_magnitude = face_normal_magnitude;
991 : if (volume_mesh.quadrature(direction.dimension()) ==
992 : Spectral::Quadrature::Gauss) {
993 : local_mortar_data.volume_mesh = volume_mesh;
994 : local_mortar_data.volume_det_inv_jacobian =
995 : volume_det_inv_jacobian;
996 : local_mortar_data.face_det_jacobian = face_det_jacobian;
997 : }
998 : ASSERT(boundary_data_history->count(mortar_id) != 0,
999 : "Could not insert the mortar data for "
1000 : << mortar_id
1001 : << " because the unordered map has not been "
1002 : "initialized "
1003 : "to have the mortar id.");
1004 : boundary_data_history->at(mortar_id).local().insert(
1005 : time_step_id, integration_order,
1006 : std::move(mortar_data->at(mortar_id).local()));
1007 : mortar_data->at(mortar_id) = MortarDataHolder<Dim>{};
1008 : }
1009 : }
1010 : },
1011 : box,
1012 : db::get<evolution::dg::Tags::NormalCovectorAndMagnitude<Dim>>(*box));
1013 : }
1014 : }
1015 : } // namespace evolution::dg::Actions
|
