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