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