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