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