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