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/BoundaryCorrection.hpp"
28 : #include "Evolution/BoundaryCorrectionTags.hpp"
29 : #include "Evolution/DiscontinuousGalerkin/Actions/BoundaryConditionsImpl.hpp"
30 : #include "Evolution/DiscontinuousGalerkin/Actions/ComputeTimeDerivativeHelpers.hpp"
31 : #include "Evolution/DiscontinuousGalerkin/Actions/InternalMortarDataImpl.hpp"
32 : #include "Evolution/DiscontinuousGalerkin/Actions/NormalCovectorAndMagnitude.hpp"
33 : #include "Evolution/DiscontinuousGalerkin/Actions/PackageDataImpl.hpp"
34 : #include "Evolution/DiscontinuousGalerkin/Actions/VolumeTermsImpl.hpp"
35 : #include "Evolution/DiscontinuousGalerkin/BoundaryData.hpp"
36 : #include "Evolution/DiscontinuousGalerkin/InboxTags.hpp"
37 : #include "Evolution/DiscontinuousGalerkin/MortarData.hpp"
38 : #include "Evolution/DiscontinuousGalerkin/MortarDataHolder.hpp"
39 : #include "Evolution/DiscontinuousGalerkin/MortarTags.hpp"
40 : #include "Evolution/DiscontinuousGalerkin/NormalVectorTags.hpp"
41 : #include "Evolution/DiscontinuousGalerkin/UsingSubcell.hpp"
42 : #include "NumericalAlgorithms/DiscontinuousGalerkin/Formulation.hpp"
43 : #include "NumericalAlgorithms/DiscontinuousGalerkin/MortarHelpers.hpp"
44 : #include "NumericalAlgorithms/DiscontinuousGalerkin/Tags.hpp"
45 : #include "NumericalAlgorithms/DiscontinuousGalerkin/Tags/Formulation.hpp"
46 : #include "NumericalAlgorithms/LinearOperators/Divergence.hpp"
47 : #include "NumericalAlgorithms/LinearOperators/PartialDerivatives.hpp"
48 : #include "NumericalAlgorithms/Spectral/Basis.hpp"
49 : #include "NumericalAlgorithms/Spectral/BoundaryInterpolationMatrices.hpp"
50 : #include "NumericalAlgorithms/Spectral/Mesh.hpp"
51 : #include "NumericalAlgorithms/Spectral/Quadrature.hpp"
52 : #include "Parallel/AlgorithmExecution.hpp"
53 : #include "Parallel/ArrayCollection/IsDgElementCollection.hpp"
54 : #include "Parallel/ArrayCollection/SendDataToElement.hpp"
55 : #include "Parallel/GlobalCache.hpp"
56 : #include "Parallel/Invoke.hpp"
57 : #include "Time/BoundaryHistory.hpp"
58 : #include "Time/ChangeStepSize.hpp"
59 : #include "Utilities/Algorithm.hpp"
60 : #include "Utilities/Gsl.hpp"
61 : #include "Utilities/TMPL.hpp"
62 :
63 : /// \cond
64 : namespace Tags {
65 : template <typename Tag>
66 : struct HistoryEvolvedVariables;
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 must be derived from
251 : * `evolution::BoundaryCorrection`. Each concrete boundary correction
252 : * must specify:
253 : *
254 : * - type alias `dg_package_field_tags`. These are what will be returned by
255 : * `gsl::not_null` from the `dg_package_data` member function.
256 : *
257 : * - type alias `dg_package_data_temporary_tags`. These are temporary tags
258 : * that are projected to the face and then passed to the `dg_package_data`
259 : * function.
260 : *
261 : * - type alias `dg_package_data_primitive_tags`. These are the primitive
262 : * variables (if any) that are projected to the face and then passed to
263 : * `dg_package_data`.
264 : *
265 : * - type alias `dg_package_data_volume_tags`. These are tags that are not
266 : * projected to the interface and are retrieved directly from the `DataBox`.
267 : * The equation of state for hydrodynamics systems is an example of what
268 : * would be a "volume tag".
269 : *
270 : * A `static constexpr bool need_normal_vector` must be specified. If `true`
271 : * then the normal vector is computed from the normal covector. This is
272 : * currently not implemented.
273 : *
274 : * The `dg_package_data` function takes as arguments `gsl::not_null` of the
275 : * `dg_package_field_tags`, then the projected evolved variables, the
276 : * projected fluxes, the projected temporaries, the projected primitives, the
277 : * unit normal covector, mesh velocity, normal dotted into the mesh velocity,
278 : * the `volume_tags`, and finally the `dg::Formulation`. The `dg_package_data`
279 : * function must compute all ingredients for the boundary correction, including
280 : * mesh-velocity-corrected characteristic speeds. However, the projected fluxes
281 : * passed in are \f$F^i - u v^i_g\f$ (the mesh velocity term is already
282 : * included). The `dg_package_data` function must also return a `double` that is
283 : * the maximum absolute characteristic speed over the entire face. This will be
284 : * used for checking that the time step doesn't violate the CFL condition.
285 : *
286 : * Here is an example of the type aliases and `bool`:
287 : *
288 : * \snippet ComputeTimeDerivativeImpl.tpp bt_ta
289 : *
290 : * The normal vector requirement is:
291 : *
292 : * \snippet ComputeTimeDerivativeImpl.tpp bt_nnv
293 : *
294 : * For a conservative system with primitive variables and using the `TimeStepId`
295 : * as a volume tag the `dg_package_data` function looks like:
296 : *
297 : * \snippet ComputeTimeDerivativeImpl.tpp bt_cp
298 : *
299 : * For a mixed conservative-nonconservative system with primitive variables and
300 : * using the `TimeStepId` as a volume tag the `dg_package_data` function looks
301 : * like:
302 : *
303 : * \snippet ComputeTimeDerivativeImpl.tpp bt_mp
304 : *
305 : * Uses:
306 : * - System:
307 : * - `boundary_correction`
308 : * - `variables_tag`
309 : * - `flux_variables`
310 : * - `gradients_tags`
311 : * - `compute_volume_time_derivative`
312 : * - `has_primitive_and_conservative_vars`
313 : * - `primitive_variables_tag` if system has primitive variables
314 : *
315 : * - DataBox:
316 : * - `domain::Tags::Element<Dim>`
317 : * - `domain::Tags::Mesh<Dim>`
318 : * - `evolution::dg::Tags::MortarMesh<Dim>`
319 : * - `evolution::dg::Tags::MortarData<Dim>`
320 : * - `Tags::TimeStepId`
321 : * - \code{.cpp}
322 : * domain::Tags::Interface<domain::Tags::InternalDirections<Dim>,
323 : * domain::Tags::Mesh<Dim - 1>>
324 : * \endcode
325 : * - \code{.cpp}
326 : * domain::Tags::Interface<
327 : * domain::Tags::InternalDirections<Dim>,
328 : * ::Tags::Normalized<
329 : * domain::Tags::UnnormalizedFaceNormal<Dim, Frame::Inertial>>>
330 : * \endcode
331 : * - \code{.cpp}
332 : * domain::Tags::Interface<
333 : * domain::Tags::InternalDirections<Dim>,
334 : * domain::Tags::MeshVelocity<Dim, Frame::Inertial>>
335 : * \endcode
336 : * - `Metavariables::system::variables_tag`
337 : * - `Metavariables::system::flux_variables`
338 : * - `Metavariables::system::primitive_tags` if exists
339 : * - boundary correction `dg_package_data_volume_tags`
340 : *
341 : * DataBox changes:
342 : * - Adds: nothing
343 : * - Removes: nothing
344 : * - Modifies:
345 : * - `evolution::dg::Tags::MortarData<Dim>`
346 : */
347 : template <size_t Dim, typename EvolutionSystem, typename DgStepChoosers,
348 : bool LocalTimeStepping, bool UseNodegroupDgElements>
349 1 : struct ComputeTimeDerivative {
350 0 : using inbox_tags =
351 : tmpl::list<evolution::dg::Tags::BoundaryCorrectionAndGhostCellsInbox<
352 : Dim, UseNodegroupDgElements>>;
353 0 : using const_global_cache_tags = tmpl::append<
354 : tmpl::list<::dg::Tags::Formulation, evolution::Tags::BoundaryCorrection,
355 : domain::Tags::ExternalBoundaryConditions<Dim>>,
356 : tmpl::conditional_t<
357 : LocalTimeStepping,
358 : typename ChangeStepSize<DgStepChoosers>::const_global_cache_tags,
359 : 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, bool UseNodegroupDgElements>
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 : UseNodegroupDgElements>::
390 : apply(db::DataBox<DbTagsList>& box,
391 : tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
392 : Parallel::GlobalCache<Metavariables>& cache,
393 : const ArrayIndex& /*array_index*/, ActionList /*meta*/,
394 : const ParallelComponent* const /*meta*/) { // NOLINT const
395 : static_assert(UseNodegroupDgElements ==
396 : Parallel::is_dg_element_collection_v<ParallelComponent>,
397 : "The action ComputeTimeDerivative is told by the "
398 : "template parameter UseNodegroupDgElements that it is being "
399 : "used with a DgElementCollection, but the ParallelComponent "
400 : "is not a DgElementCollection. You need to change the "
401 : "template parameter on the ComputeTimeDerivative action "
402 : "in your action list.");
403 :
404 : using variables_tag = typename EvolutionSystem::variables_tag;
405 : using dt_variables_tag = db::add_tag_prefix<::Tags::dt, variables_tag>;
406 : using partial_derivative_tags = typename EvolutionSystem::gradient_variables;
407 : using flux_variables = typename EvolutionSystem::flux_variables;
408 : using compute_volume_time_derivative_terms =
409 : typename EvolutionSystem::compute_volume_time_derivative_terms;
410 :
411 : const Mesh<Dim>& mesh = db::get<::domain::Tags::Mesh<Dim>>(box);
412 : const Element<Dim>& element = db::get<domain::Tags::Element<Dim>>(box);
413 : const ::dg::Formulation dg_formulation =
414 : db::get<::dg::Tags::Formulation>(box);
415 : ASSERT(alg::all_of(mesh.basis(),
416 : [&mesh](const Spectral::Basis current_basis) {
417 : return current_basis == mesh.basis(0);
418 : }) or
419 : element.topologies() != domain::topologies::hypercube<Dim>,
420 : "An isotropic basis must be used in the evolution code. While "
421 : "theoretically this restriction could be lifted, the simplification "
422 : "it offers are quite substantial. Relaxing this assumption is likely "
423 : "to require quite a bit of careful code refactoring and debugging.");
424 : ASSERT(alg::all_of(mesh.quadrature(),
425 : [&mesh](const Spectral::Quadrature current_quadrature) {
426 : return current_quadrature == mesh.quadrature(0);
427 : }) or
428 : element.topologies() != domain::topologies::hypercube<Dim>,
429 : "An isotropic quadrature must be used in the evolution code. While "
430 : "theoretically this restriction could be lifted, the simplification "
431 : "it offers are quite substantial. Relaxing this assumption is likely "
432 : "to require quite a bit of careful code refactoring and debugging.");
433 :
434 : const auto& boundary_correction =
435 : db::get<evolution::Tags::BoundaryCorrection>(box);
436 : using derived_boundary_corrections =
437 : tmpl::at<typename Metavariables::factory_creation::factory_classes,
438 : evolution::BoundaryCorrection>;
439 :
440 : // To avoid a second allocation in internal_mortar_data, we allocate the
441 : // variables needed to construct the fields on the faces here along with
442 : // everything else. This requires us to know all the tags necessary to apply
443 : // boundary corrections. However, since we pick boundary corrections at
444 : // runtime, we just gather all possible tags from all possible boundary
445 : // corrections and lump them into the allocation. This may result in a
446 : // larger-than-necessary allocation, but it won't be that much larger.
447 : using all_dg_package_temporary_tags =
448 : tmpl::transform<derived_boundary_corrections,
449 : detail::get_dg_package_temporary_tags<tmpl::_1>>;
450 : using all_primitive_tags_for_face =
451 : tmpl::transform<derived_boundary_corrections,
452 : detail::get_primitive_tags_for_face<
453 : tmpl::pin<EvolutionSystem>, tmpl::_1>>;
454 : using fluxes_tags = db::wrap_tags_in<::Tags::Flux, flux_variables,
455 : tmpl::size_t<Dim>, Frame::Inertial>;
456 : using dg_package_data_projected_tags =
457 : tmpl::list<typename variables_tag::tags_list, fluxes_tags,
458 : all_dg_package_temporary_tags, all_primitive_tags_for_face>;
459 : using all_face_temporary_tags =
460 : tmpl::remove_duplicates<tmpl::flatten<tmpl::push_back<
461 : tmpl::list<dg_package_data_projected_tags,
462 : detail::inverse_spatial_metric_tag<EvolutionSystem>>,
463 : detail::OneOverNormalVectorMagnitude, detail::NormalVector<Dim>>>>;
464 : // To avoid additional allocations in internal_mortar_data, we provide a
465 : // buffer used to compute the packaged data before it has to be projected to
466 : // the mortar. We get all mortar tags for similar reasons as described above
467 : using all_mortar_tags = tmpl::remove_duplicates<tmpl::flatten<
468 : tmpl::transform<derived_boundary_corrections,
469 : detail::get_dg_package_field_tags<tmpl::_1>>>>;
470 :
471 : // We also don't use the number of volume mesh grid points. We instead use the
472 : // max number of grid points from each face. That way, our allocation will be
473 : // large enough to hold any face and we can reuse the allocation for each face
474 : // without having to resize it.
475 : size_t num_face_temporary_grid_points = 0;
476 : {
477 : for (const auto& [direction, neighbors_in_direction] :
478 : element.neighbors()) {
479 : (void)neighbors_in_direction;
480 : const auto face_mesh = mesh.slice_away(direction.dimension());
481 : num_face_temporary_grid_points = std::max(
482 : num_face_temporary_grid_points, face_mesh.number_of_grid_points());
483 : }
484 : }
485 :
486 : // Allocate the Variables classes needed for the time derivative
487 : // computation.
488 : //
489 : // This is factored out so that we will be able to do ADER-DG/CG where a
490 : // spacetime polynomial is constructed by solving implicit equations in time
491 : // using a Picard iteration. A high-order initial guess is needed to
492 : // efficiently construct the ADER spacetime solution. This initial guess is
493 : // obtained using continuous RK methods, and so we will want to reuse
494 : // buffers. Thus, the volume_terms function returns by reference rather than
495 : // by value.
496 : using VarsTemporaries =
497 : Variables<typename compute_volume_time_derivative_terms::temporary_tags>;
498 : using VarsFluxes =
499 : Variables<db::wrap_tags_in<::Tags::Flux, flux_variables,
500 : tmpl::size_t<Dim>, Frame::Inertial>>;
501 : using VarsPartialDerivatives =
502 : Variables<db::wrap_tags_in<::Tags::deriv, partial_derivative_tags,
503 : tmpl::size_t<Dim>, Frame::Inertial>>;
504 : using VarsDivFluxes = Variables<db::wrap_tags_in<
505 : ::Tags::div, db::wrap_tags_in<::Tags::Flux, flux_variables,
506 : tmpl::size_t<Dim>, Frame::Inertial>>>;
507 : using VarsFaceTemporaries = Variables<all_face_temporary_tags>;
508 : using DgPackagedDataVarsOnFace = Variables<all_mortar_tags>;
509 : const size_t number_of_grid_points = mesh.number_of_grid_points();
510 : const size_t buffer_size =
511 : (VarsTemporaries::number_of_independent_components +
512 : VarsFluxes::number_of_independent_components +
513 : VarsPartialDerivatives::number_of_independent_components +
514 : VarsDivFluxes::number_of_independent_components) *
515 : number_of_grid_points +
516 : // Different number of grid points. See explanation above where
517 : // num_face_temporary_grid_points is defined
518 : (VarsFaceTemporaries::number_of_independent_components +
519 : DgPackagedDataVarsOnFace::number_of_independent_components) *
520 : num_face_temporary_grid_points;
521 : auto buffer = cpp20::make_unique_for_overwrite<double[]>(buffer_size);
522 : #ifdef SPECTRE_DEBUG
523 : std::fill(&buffer[0], &buffer[buffer_size],
524 : std::numeric_limits<double>::signaling_NaN());
525 : #endif
526 : VarsTemporaries temporaries{
527 : &buffer[0], VarsTemporaries::number_of_independent_components *
528 : number_of_grid_points};
529 : VarsFluxes volume_fluxes{
530 : &buffer[VarsTemporaries::number_of_independent_components *
531 : number_of_grid_points],
532 : VarsFluxes::number_of_independent_components * number_of_grid_points};
533 : VarsPartialDerivatives partial_derivs{
534 : &buffer[(VarsTemporaries::number_of_independent_components +
535 : VarsFluxes::number_of_independent_components) *
536 : number_of_grid_points],
537 : VarsPartialDerivatives::number_of_independent_components *
538 : number_of_grid_points};
539 : VarsDivFluxes div_fluxes{
540 : &buffer[(VarsTemporaries::number_of_independent_components +
541 : VarsFluxes::number_of_independent_components +
542 : VarsPartialDerivatives::number_of_independent_components) *
543 : number_of_grid_points],
544 : VarsDivFluxes::number_of_independent_components * number_of_grid_points};
545 : // Lighter weight data structure than a Variables to avoid passing even more
546 : // templates to internal_mortar_data.
547 : gsl::span<double> face_temporaries = gsl::make_span<double>(
548 : &buffer[(VarsTemporaries::number_of_independent_components +
549 : VarsFluxes::number_of_independent_components +
550 : VarsPartialDerivatives::number_of_independent_components +
551 : VarsDivFluxes::number_of_independent_components) *
552 : number_of_grid_points],
553 : // Different number of grid points. See explanation above where
554 : // num_face_temporary_grid_points is defined
555 : VarsFaceTemporaries::number_of_independent_components *
556 : num_face_temporary_grid_points);
557 : gsl::span<double> packaged_data_buffer = gsl::make_span<double>(
558 : &buffer[(VarsTemporaries::number_of_independent_components +
559 : VarsFluxes::number_of_independent_components +
560 : VarsPartialDerivatives::number_of_independent_components +
561 : VarsDivFluxes::number_of_independent_components) *
562 : number_of_grid_points +
563 : VarsFaceTemporaries::number_of_independent_components *
564 : num_face_temporary_grid_points],
565 : // Different number of grid points. See explanation above where
566 : // num_face_temporary_grid_points is defined
567 : DgPackagedDataVarsOnFace::number_of_independent_components *
568 : num_face_temporary_grid_points);
569 :
570 : const Scalar<DataVector>* det_inverse_jacobian = nullptr;
571 : if constexpr (tmpl::size<flux_variables>::value != 0) {
572 : if (dg_formulation == ::dg::Formulation::WeakInertial) {
573 : det_inverse_jacobian = &db::get<
574 : domain::Tags::DetInvJacobian<Frame::ElementLogical, Frame::Inertial>>(
575 : box);
576 : }
577 : }
578 : db::mutate_apply<
579 : tmpl::list<dt_variables_tag>,
580 : typename compute_volume_time_derivative_terms::argument_tags>(
581 : [&dg_formulation, &div_fluxes, &det_inverse_jacobian,
582 : &div_mesh_velocity = db::get<::domain::Tags::DivMeshVelocity>(box),
583 : &evolved_variables = db::get<variables_tag>(box),
584 : &inertial_coordinates =
585 : db::get<domain::Tags::Coordinates<Dim, Frame::Inertial>>(box),
586 : &logical_to_inertial_inv_jacobian =
587 : db::get<::domain::Tags::InverseJacobian<Dim, Frame::ElementLogical,
588 : Frame::Inertial>>(box),
589 : &mesh, &mesh_velocity = db::get<::domain::Tags::MeshVelocity<Dim>>(box),
590 : &partial_derivs, &temporaries, &volume_fluxes](
591 : const gsl::not_null<Variables<
592 : db::wrap_tags_in<::Tags::dt, typename variables_tag::tags_list>>*>
593 : dt_vars_ptr,
594 : const auto&... time_derivative_args) {
595 : detail::volume_terms<compute_volume_time_derivative_terms>(
596 : dt_vars_ptr, make_not_null(&volume_fluxes),
597 : make_not_null(&partial_derivs), make_not_null(&temporaries),
598 : make_not_null(&div_fluxes), evolved_variables, dg_formulation, mesh,
599 : inertial_coordinates, logical_to_inertial_inv_jacobian,
600 : det_inverse_jacobian, mesh_velocity, div_mesh_velocity,
601 : time_derivative_args...);
602 : },
603 : make_not_null(&box));
604 :
605 : const Variables<detail::get_primitive_vars_tags_from_system<EvolutionSystem>>*
606 : primitive_vars{nullptr};
607 : if constexpr (EvolutionSystem::has_primitive_and_conservative_vars) {
608 : primitive_vars =
609 : &db::get<typename EvolutionSystem::primitive_variables_tag>(box);
610 : }
611 :
612 : static_assert(
613 : tmpl::all<derived_boundary_corrections, std::is_final<tmpl::_1>>::value,
614 : "All createable classes for boundary corrections must be marked "
615 : "final.");
616 : tmpl::for_each<derived_boundary_corrections>(
617 : [&boundary_correction, &box, &partial_derivs, &primitive_vars,
618 : &temporaries, &volume_fluxes, &packaged_data_buffer,
619 : &face_temporaries](auto derived_correction_v) {
620 : using DerivedCorrection =
621 : tmpl::type_from<decltype(derived_correction_v)>;
622 : if (typeid(boundary_correction) == typeid(DerivedCorrection)) {
623 : // Compute internal boundary quantities on the mortar for sides
624 : // of the element that have neighbors, i.e. they are not an
625 : // external side.
626 : // Note: this call mutates:
627 : // - evolution::dg::Tags::NormalCovectorAndMagnitude<Dim>,
628 : // - evolution::dg::Tags::MortarData<Dim>
629 : detail::internal_mortar_data<EvolutionSystem, Dim>(
630 : make_not_null(&box), make_not_null(&face_temporaries),
631 : make_not_null(&packaged_data_buffer),
632 : dynamic_cast<const DerivedCorrection&>(boundary_correction),
633 : db::get<variables_tag>(box), volume_fluxes, temporaries,
634 : primitive_vars,
635 : typename DerivedCorrection::dg_package_data_volume_tags{});
636 :
637 : detail::apply_boundary_conditions_on_all_external_faces<
638 : EvolutionSystem, Dim>(
639 : make_not_null(&box),
640 : dynamic_cast<const DerivedCorrection&>(boundary_correction),
641 : temporaries, volume_fluxes, partial_derivs, primitive_vars);
642 : }
643 : });
644 :
645 : if constexpr (LocalTimeStepping) {
646 : db::mutate_apply<ChangeStepSize<DgStepChoosers>>(make_not_null(&box));
647 : }
648 :
649 : send_data_for_fluxes<ParallelComponent>(make_not_null(&cache),
650 : make_not_null(&box), volume_fluxes);
651 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
652 : }
653 :
654 : template <size_t Dim, typename EvolutionSystem, typename DgStepChoosers,
655 : bool LocalTimeStepping, bool UseNodegroupDgElements>
656 : template <typename ParallelComponent, typename DbTagsList,
657 : typename Metavariables>
658 : void ComputeTimeDerivative<Dim, EvolutionSystem, DgStepChoosers,
659 0 : LocalTimeStepping, UseNodegroupDgElements>::
660 : send_data_for_fluxes(
661 : const gsl::not_null<Parallel::GlobalCache<Metavariables>*> cache,
662 : const gsl::not_null<db::DataBox<DbTagsList>*> box,
663 : [[maybe_unused]] const Variables<db::wrap_tags_in<
664 : ::Tags::Flux, typename EvolutionSystem::flux_variables,
665 : tmpl::size_t<Dim>, Frame::Inertial>>& volume_fluxes) {
666 : using variables_tag = typename EvolutionSystem::variables_tag;
667 :
668 : auto& receiver_proxy =
669 : Parallel::get_parallel_component<ParallelComponent>(*cache);
670 : const auto& element = db::get<domain::Tags::Element<Dim>>(*box);
671 :
672 : const auto& time_step_id = db::get<::Tags::TimeStepId>(*box);
673 : const auto integration_order =
674 : db::get<::Tags::HistoryEvolvedVariables<variables_tag>>(*box)
675 : .integration_order();
676 : const auto& all_mortar_data =
677 : db::get<evolution::dg::Tags::MortarData<Dim>>(*box);
678 : const auto& mortar_meshes = get<evolution::dg::Tags::MortarMesh<Dim>>(*box);
679 :
680 : std::optional<DirectionMap<Dim, DataVector>>
681 : all_neighbor_data_for_reconstruction = std::nullopt;
682 : int tci_decision = 0;
683 : const Mesh<Dim>& volume_mesh = db::get<domain::Tags::Mesh<Dim>>(*box);
684 : std::optional<Mesh<Dim>> ghost_data_mesh = std::nullopt;
685 : if constexpr (using_subcell_v<Metavariables>) {
686 : if (not all_neighbor_data_for_reconstruction.has_value()) {
687 : all_neighbor_data_for_reconstruction = DirectionMap<Dim, DataVector>{};
688 : }
689 :
690 : evolution::dg::subcell::prepare_neighbor_data<Metavariables>(
691 : make_not_null(&all_neighbor_data_for_reconstruction.value()),
692 : make_not_null(&ghost_data_mesh), box, volume_fluxes);
693 : tci_decision = evolution::dg::subcell::get_tci_decision(*box);
694 : }
695 :
696 : for (const auto& [direction, neighbors] : element.neighbors()) {
697 : DataVector ghost_and_subcell_data{};
698 : if constexpr (using_subcell_v<Metavariables>) {
699 : ASSERT(all_neighbor_data_for_reconstruction.has_value(),
700 : "Trying to do DG-subcell but the ghost and subcell data for the "
701 : "neighbor has not been set.");
702 : ghost_and_subcell_data =
703 : std::move(all_neighbor_data_for_reconstruction.value()[direction]);
704 : }
705 :
706 : const size_t total_neighbors = neighbors.size();
707 : size_t neighbor_count = 1;
708 : for (const auto& neighbor : neighbors) {
709 : const auto& orientation = neighbors.orientation(neighbor);
710 : const auto direction_from_neighbor = orientation(direction.opposite());
711 : const DirectionalId<Dim> mortar_id{direction, neighbor};
712 :
713 : const Mesh<Dim - 1>& mortar_mesh = mortar_meshes.at(mortar_id);
714 : const auto volume_mesh_for_neighbor = orientation(volume_mesh);
715 : const auto mortar_mesh_for_neighbor =
716 : orient_mesh_on_slice(mortar_mesh, direction.dimension(), orientation);
717 : DataVector neighbor_boundary_data_on_mortar{};
718 :
719 : if (LIKELY(orientation.is_aligned())) {
720 : neighbor_boundary_data_on_mortar =
721 : *all_mortar_data.at(mortar_id).local().mortar_data.value();
722 : } else {
723 : const auto& slice_extents = mortar_mesh.extents();
724 : neighbor_boundary_data_on_mortar = orient_variables_on_slice(
725 : all_mortar_data.at(mortar_id).local().mortar_data.value(),
726 : slice_extents, direction.dimension(), orientation);
727 : }
728 :
729 : const TimeStepId& next_time_step_id =
730 : db::get<::Tags::Next<::Tags::TimeStepId>>(*box);
731 :
732 : using SendData = evolution::dg::BoundaryData<Dim>;
733 : SendData data{};
734 :
735 : if (neighbor_count == total_neighbors) {
736 : data = SendData{volume_mesh_for_neighbor,
737 : ghost_data_mesh,
738 : mortar_mesh_for_neighbor,
739 : std::move(ghost_and_subcell_data),
740 : {std::move(neighbor_boundary_data_on_mortar)},
741 : next_time_step_id,
742 : tci_decision,
743 : integration_order,
744 : std::nullopt};
745 : } else {
746 : data = SendData{volume_mesh_for_neighbor,
747 : ghost_data_mesh,
748 : mortar_mesh_for_neighbor,
749 : ghost_and_subcell_data,
750 : {std::move(neighbor_boundary_data_on_mortar)},
751 : next_time_step_id,
752 : tci_decision,
753 : integration_order,
754 : std::nullopt};
755 : }
756 :
757 : // Send mortar data (the `std::tuple` named `data`) to neighbor
758 : if constexpr (Parallel::is_dg_element_collection_v<ParallelComponent>) {
759 : Parallel::local_synchronous_action<
760 : Parallel::Actions::SendDataToElement>(
761 : receiver_proxy, cache,
762 : evolution::dg::Tags::BoundaryCorrectionAndGhostCellsInbox<
763 : Dim, UseNodegroupDgElements>{},
764 : neighbor, time_step_id,
765 : std::make_pair(DirectionalId{direction_from_neighbor, element.id()},
766 : std::move(data)));
767 : } else {
768 : Parallel::receive_data<
769 : evolution::dg::Tags::BoundaryCorrectionAndGhostCellsInbox<
770 : Dim, UseNodegroupDgElements>>(
771 : receiver_proxy[neighbor], time_step_id,
772 : std::make_pair(DirectionalId{direction_from_neighbor, element.id()},
773 : std::move(data)));
774 : }
775 : ++neighbor_count;
776 : }
777 : }
778 :
779 : if constexpr (LocalTimeStepping) {
780 : using dt_variables_tag = db::add_tag_prefix<::Tags::dt, variables_tag>;
781 : // We assume isotropic quadrature, i.e. the quadrature is the same in
782 : // all directions.
783 : const bool using_gauss_points =
784 : db::get<domain::Tags::Mesh<Dim>>(*box).quadrature() ==
785 : make_array<Dim>(Spectral::Quadrature::Gauss);
786 :
787 : const Scalar<DataVector> volume_det_inv_jacobian{};
788 : if (using_gauss_points) {
789 : // NOLINTNEXTLINE
790 : const_cast<DataVector&>(get(volume_det_inv_jacobian))
791 : .set_data_ref(make_not_null(&const_cast<DataVector&>( // NOLINT
792 : get(db::get<domain::Tags::DetInvJacobian<
793 : Frame::ElementLogical, Frame::Inertial>>(*box)))));
794 : }
795 :
796 : // Add face normal and Jacobian determinants to the local mortar data. We
797 : // only need the Jacobians if we are using Gauss points. Then copy over
798 : // into the boundary history, since that's what the LTS steppers use.
799 : //
800 : // The boundary history coupling computation (which computes the _lifted_
801 : // boundary correction) returns a Variables<dt<EvolvedVars>> instead of
802 : // using the `NormalDotNumericalFlux` prefix tag. This is because the
803 : // returned quantity is more a `dt` quantity than a
804 : // `NormalDotNormalDotFlux` since it's been lifted to the volume.
805 : db::mutate<evolution::dg::Tags::MortarData<Dim>,
806 : evolution::dg::Tags::MortarDataHistory<
807 : Dim, typename dt_variables_tag::type>>(
808 : [&element, integration_order, &time_step_id, using_gauss_points,
809 : &volume_det_inv_jacobian, &volume_mesh](
810 : const gsl::not_null<
811 : DirectionalIdMap<Dim, evolution::dg::MortarDataHolder<Dim>>*>
812 : mortar_data,
813 : const gsl::not_null<
814 : DirectionalIdMap<Dim, TimeSteppers::BoundaryHistory<
815 : evolution::dg::MortarData<Dim>,
816 : evolution::dg::MortarData<Dim>,
817 : typename dt_variables_tag::type>>*>
818 : boundary_data_history,
819 : const DirectionMap<Dim,
820 : std::optional<Variables<tmpl::list<
821 : evolution::dg::Tags::MagnitudeOfNormal,
822 : evolution::dg::Tags::NormalCovector<Dim>>>>>&
823 : normal_covector_and_magnitude) {
824 : Scalar<DataVector> volume_det_jacobian{};
825 : Scalar<DataVector> face_det_jacobian{};
826 : if (using_gauss_points) {
827 : get(volume_det_jacobian) = 1.0 / get(volume_det_inv_jacobian);
828 : }
829 : for (const auto& [direction, neighbors_in_direction] :
830 : element.neighbors()) {
831 : // We can perform projections once for all neighbors in the
832 : // direction because we care about the _face_ mesh, not the mortar
833 : // mesh.
834 : ASSERT(normal_covector_and_magnitude.at(direction).has_value(),
835 : "The normal covector and magnitude have not been computed.");
836 : const Scalar<DataVector>& face_normal_magnitude =
837 : get<evolution::dg::Tags::MagnitudeOfNormal>(
838 : *normal_covector_and_magnitude.at(direction));
839 : if (using_gauss_points) {
840 : const Matrix identity{};
841 : auto interpolation_matrices =
842 : make_array<Dim>(std::cref(identity));
843 : const std::pair<Matrix, Matrix>& matrices =
844 : Spectral::boundary_interpolation_matrices(
845 : volume_mesh.slice_through(direction.dimension()));
846 : gsl::at(interpolation_matrices, direction.dimension()) =
847 : direction.side() == Side::Upper ? matrices.second
848 : : matrices.first;
849 : if (get(face_det_jacobian).size() !=
850 : get(face_normal_magnitude).size()) {
851 : get(face_det_jacobian) =
852 : DataVector{get(face_normal_magnitude).size()};
853 : }
854 : apply_matrices(make_not_null(&get(face_det_jacobian)),
855 : interpolation_matrices, get(volume_det_jacobian),
856 : volume_mesh.extents());
857 : }
858 :
859 : for (const auto& neighbor : neighbors_in_direction) {
860 : const DirectionalId<Dim> mortar_id{direction, neighbor};
861 : auto& local_mortar_data = mortar_data->at(mortar_id).local();
862 : local_mortar_data.face_normal_magnitude = face_normal_magnitude;
863 : if (using_gauss_points) {
864 : local_mortar_data.volume_mesh = volume_mesh;
865 : local_mortar_data.volume_det_inv_jacobian =
866 : volume_det_inv_jacobian;
867 : local_mortar_data.face_det_jacobian = face_det_jacobian;
868 : }
869 : ASSERT(boundary_data_history->count(mortar_id) != 0,
870 : "Could not insert the mortar data for "
871 : << mortar_id
872 : << " because the unordered map has not been "
873 : "initialized "
874 : "to have the mortar id.");
875 : boundary_data_history->at(mortar_id).local().insert(
876 : time_step_id, integration_order,
877 : std::move(mortar_data->at(mortar_id).local()));
878 : mortar_data->at(mortar_id) = MortarDataHolder<Dim>{};
879 : }
880 : }
881 : },
882 : box,
883 : db::get<evolution::dg::Tags::NormalCovectorAndMagnitude<Dim>>(*box));
884 : }
885 : }
886 : } // namespace evolution::dg::Actions
|
