Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : #pragma once
5 :
6 : #include <cstddef>
7 : #include <limits>
8 : #include <map>
9 : #include <optional>
10 : #include <tuple>
11 : #include <type_traits>
12 : #include <utility>
13 : #include <vector>
14 :
15 : #include "DataStructures/DataBox/AsAccess.hpp"
16 : #include "DataStructures/DataBox/DataBox.hpp"
17 : #include "DataStructures/DataBox/PrefixHelpers.hpp"
18 : #include "DataStructures/DataBox/Prefixes.hpp"
19 : #include "DataStructures/Tensor/EagerMath/Magnitude.hpp"
20 : #include "Domain/FaceNormal.hpp"
21 : #include "Domain/Structure/DirectionalIdMap.hpp"
22 : #include "Domain/Structure/Element.hpp"
23 : #include "Domain/Structure/ElementId.hpp"
24 : #include "Domain/Structure/TrimMap.hpp"
25 : #include "Domain/Tags.hpp"
26 : #include "Domain/Tags/NeighborMesh.hpp"
27 : #include "Evolution/BoundaryCorrectionTags.hpp"
28 : #include "Evolution/DiscontinuousGalerkin/InboxTags.hpp"
29 : #include "Evolution/DiscontinuousGalerkin/MortarData.hpp"
30 : #include "Evolution/DiscontinuousGalerkin/MortarTags.hpp"
31 : #include "Evolution/DiscontinuousGalerkin/NormalVectorTags.hpp"
32 : #include "Evolution/DiscontinuousGalerkin/UsingSubcell.hpp"
33 : #include "NumericalAlgorithms/DiscontinuousGalerkin/Formulation.hpp"
34 : #include "NumericalAlgorithms/DiscontinuousGalerkin/LiftFlux.hpp"
35 : #include "NumericalAlgorithms/DiscontinuousGalerkin/LiftFromBoundary.hpp"
36 : #include "NumericalAlgorithms/DiscontinuousGalerkin/MortarHelpers.hpp"
37 : #include "NumericalAlgorithms/DiscontinuousGalerkin/Tags/Formulation.hpp"
38 : #include "NumericalAlgorithms/Spectral/Mesh.hpp"
39 : #include "NumericalAlgorithms/Spectral/Projection.hpp"
40 : #include "NumericalAlgorithms/Spectral/Quadrature.hpp"
41 : #include "NumericalAlgorithms/Spectral/Spectral.hpp"
42 : #include "Parallel/AlgorithmExecution.hpp"
43 : #include "Parallel/GlobalCache.hpp"
44 : #include "Time/BoundaryHistory.hpp"
45 : #include "Time/EvolutionOrdering.hpp"
46 : #include "Time/Time.hpp"
47 : #include "Time/TimeStepId.hpp"
48 : #include "Time/TimeSteppers/LtsTimeStepper.hpp"
49 : #include "Time/TimeSteppers/TimeStepper.hpp"
50 : #include "Utilities/Algorithm.hpp"
51 : #include "Utilities/CallWithDynamicType.hpp"
52 : #include "Utilities/ErrorHandling/Error.hpp"
53 : #include "Utilities/Gsl.hpp"
54 : #include "Utilities/MakeArray.hpp"
55 : #include "Utilities/TMPL.hpp"
56 : #include "Utilities/TaggedTuple.hpp"
57 :
58 : /// \cond
59 : namespace Tags {
60 : struct Time;
61 : struct TimeStep;
62 : struct TimeStepId;
63 : template <typename StepperInterface>
64 : struct TimeStepper;
65 : } // namespace Tags
66 : /// \endcond
67 :
68 : /// \cond
69 : namespace evolution::dg::subcell {
70 : // We use a forward declaration instead of including a header file to avoid
71 : // coupling to the DG-subcell libraries for executables that don't use subcell.
72 : template <size_t VolumeDim, typename DgComputeSubcellNeighborPackagedData>
73 : void neighbor_reconstructed_face_solution(
74 : gsl::not_null<db::Access*> box,
75 : gsl::not_null<std::pair<
76 : const TimeStepId,
77 : DirectionalIdMap<
78 : VolumeDim,
79 : std::tuple<Mesh<VolumeDim>, Mesh<VolumeDim - 1>,
80 : std::optional<DataVector>, std::optional<DataVector>,
81 : ::TimeStepId, int>>>*>
82 : received_temporal_id_and_data);
83 : template <size_t Dim>
84 : void neighbor_tci_decision(
85 : gsl::not_null<db::Access*> box,
86 : const std::pair<
87 : const TimeStepId,
88 : DirectionalIdMap<
89 : Dim, std::tuple<Mesh<Dim>, Mesh<Dim - 1>, std::optional<DataVector>,
90 : std::optional<DataVector>, ::TimeStepId, int>>>&
91 : received_temporal_id_and_data);
92 : } // namespace evolution::dg::subcell
93 : /// \endcond
94 :
95 : namespace evolution::dg {
96 : namespace detail {
97 : template <typename BoundaryCorrectionClass>
98 : struct get_dg_boundary_terms {
99 : using type = typename BoundaryCorrectionClass::dg_boundary_terms_volume_tags;
100 : };
101 :
102 : template <typename Tag, typename Type = db::const_item_type<Tag, tmpl::list<>>>
103 : struct TemporaryReference {
104 : using tag = Tag;
105 : using type = const Type&;
106 : };
107 : } // namespace detail
108 :
109 : /// Receive boundary data for global time-stepping. Returns true if
110 : /// all necessary data has been received.
111 : template <typename Metavariables, typename DbTagsList, typename... InboxTags>
112 1 : bool receive_boundary_data_global_time_stepping(
113 : const gsl::not_null<db::DataBox<DbTagsList>*> box,
114 : const gsl::not_null<tuples::TaggedTuple<InboxTags...>*> inboxes) {
115 : constexpr size_t volume_dim = Metavariables::system::volume_dim;
116 :
117 : const TimeStepId& temporal_id = get<::Tags::TimeStepId>(*box);
118 : using Key = DirectionalId<volume_dim>;
119 : std::map<TimeStepId,
120 : DirectionalIdMap<
121 : volume_dim,
122 : std::tuple<Mesh<volume_dim>, Mesh<volume_dim - 1>,
123 : std::optional<DataVector>, std::optional<DataVector>,
124 : ::TimeStepId, int>>>& inbox =
125 : tuples::get<evolution::dg::Tags::BoundaryCorrectionAndGhostCellsInbox<
126 : volume_dim>>(*inboxes);
127 : const auto received_temporal_id_and_data = inbox.find(temporal_id);
128 : if (received_temporal_id_and_data == inbox.end()) {
129 : return false;
130 : }
131 : const auto& received_neighbor_data = received_temporal_id_and_data->second;
132 : const Element<volume_dim>& element =
133 : db::get<domain::Tags::Element<volume_dim>>(*box);
134 : for (const auto& [direction, neighbors] : element.neighbors()) {
135 : for (const auto& neighbor : neighbors) {
136 : const auto neighbor_received =
137 : received_neighbor_data.find(Key{direction, neighbor});
138 : if (neighbor_received == received_neighbor_data.end()) {
139 : return false;
140 : }
141 : }
142 : }
143 :
144 : // Move inbox contents into the DataBox
145 : if constexpr (using_subcell_v<Metavariables>) {
146 : evolution::dg::subcell::neighbor_reconstructed_face_solution<
147 : volume_dim, typename Metavariables::SubcellOptions::
148 : DgComputeSubcellNeighborPackagedData>(
149 : &db::as_access(*box), make_not_null(&*received_temporal_id_and_data));
150 : evolution::dg::subcell::neighbor_tci_decision<volume_dim>(
151 : make_not_null(&db::as_access(*box)), *received_temporal_id_and_data);
152 : }
153 :
154 : db::mutate<evolution::dg::Tags::MortarData<volume_dim>,
155 : evolution::dg::Tags::MortarNextTemporalId<volume_dim>,
156 : domain::Tags::NeighborMesh<volume_dim>>(
157 : [&received_temporal_id_and_data](
158 : const gsl::not_null<std::unordered_map<
159 : Key, evolution::dg::MortarData<volume_dim>, boost::hash<Key>>*>
160 : mortar_data,
161 : const gsl::not_null<
162 : std::unordered_map<Key, TimeStepId, boost::hash<Key>>*>
163 : mortar_next_time_step_id,
164 : const gsl::not_null<DirectionalIdMap<volume_dim, Mesh<volume_dim>>*>
165 : neighbor_mesh) {
166 : neighbor_mesh->clear();
167 : for (auto& received_mortar_data :
168 : received_temporal_id_and_data->second) {
169 : const auto& mortar_id = received_mortar_data.first;
170 : ASSERT(received_temporal_id_and_data->first ==
171 : mortar_data->at(mortar_id).time_step_id(),
172 : "Expected to receive mortar data on mortar "
173 : << mortar_id << " at time "
174 : << mortar_next_time_step_id->at(mortar_id)
175 : << " but actually received at time "
176 : << received_temporal_id_and_data->first);
177 : neighbor_mesh->insert_or_assign(
178 : mortar_id, std::get<0>(received_mortar_data.second));
179 : mortar_next_time_step_id->at(mortar_id) =
180 : std::get<4>(received_mortar_data.second);
181 : ASSERT(using_subcell_v<Metavariables> or
182 : std::get<3>(received_mortar_data.second).has_value(),
183 : "Must receive number boundary correction data when not using "
184 : "DG-subcell.");
185 : if (std::get<3>(received_mortar_data.second).has_value()) {
186 : mortar_data->at(mortar_id).insert_neighbor_mortar_data(
187 : received_temporal_id_and_data->first,
188 : std::get<1>(received_mortar_data.second),
189 : std::move(*std::get<3>(received_mortar_data.second)));
190 : }
191 : }
192 : },
193 : box);
194 : inbox.erase(received_temporal_id_and_data);
195 : return true;
196 : }
197 :
198 : /// Receive boundary data for local time-stepping. Returns true if
199 : /// all necessary data has been received.
200 : ///
201 : /// Setting \p DenseOutput to true receives data required for output
202 : /// at `::Tags::Time` instead of `::Tags::Next<::Tags::TimeStepId>`.
203 : template <typename System, size_t Dim, bool DenseOutput, typename DbTagsList,
204 : typename... InboxTags>
205 1 : bool receive_boundary_data_local_time_stepping(
206 : const gsl::not_null<db::DataBox<DbTagsList>*> box,
207 : const gsl::not_null<tuples::TaggedTuple<InboxTags...>*> inboxes) {
208 : using variables_tag = typename System::variables_tag;
209 : using dt_variables_tag = db::add_tag_prefix<::Tags::dt, variables_tag>;
210 :
211 : // The boundary history coupling computation (which computes the _lifted_
212 : // boundary correction) returns a Variables<dt<EvolvedVars>> instead of
213 : // using the `NormalDotNumericalFlux` prefix tag. This is because the
214 : // returned quantity is more a `dt` quantity than a
215 : // `NormalDotNormalDotFlux` since it's been lifted to the volume.
216 : using Key = DirectionalId<Dim>;
217 : std::map<TimeStepId,
218 : DirectionalIdMap<
219 : Dim,
220 : std::tuple<Mesh<Dim>, Mesh<Dim - 1>,
221 : std::optional<DataVector>, std::optional<DataVector>,
222 : ::TimeStepId, int>>>& inbox =
223 : tuples::get<evolution::dg::Tags::BoundaryCorrectionAndGhostCellsInbox<
224 : Dim>>(*inboxes);
225 :
226 : const auto needed_time = [&box]() {
227 : const LtsTimeStepper& time_stepper =
228 : db::get<::Tags::TimeStepper<LtsTimeStepper>>(*box);
229 : if constexpr (DenseOutput) {
230 : const auto& dense_output_time = db::get<::Tags::Time>(*box);
231 : return [&dense_output_time, &time_stepper](const TimeStepId& id) {
232 : return time_stepper.neighbor_data_required(dense_output_time, id);
233 : };
234 : } else {
235 : const auto& next_temporal_id =
236 : db::get<::Tags::Next<::Tags::TimeStepId>>(*box);
237 : return [&next_temporal_id, &time_stepper](const TimeStepId& id) {
238 : return time_stepper.neighbor_data_required(next_temporal_id, id);
239 : };
240 : }
241 : }();
242 :
243 : const bool have_all_intermediate_messages = db::mutate<
244 : evolution::dg::Tags::MortarDataHistory<Dim,
245 : typename dt_variables_tag::type>,
246 : evolution::dg::Tags::MortarNextTemporalId<Dim>,
247 : domain::Tags::NeighborMesh<Dim>>(
248 : [&inbox, &needed_time](
249 : const gsl::not_null<
250 : std::unordered_map<Key,
251 : TimeSteppers::BoundaryHistory<
252 : evolution::dg::MortarData<Dim>,
253 : evolution::dg::MortarData<Dim>,
254 : typename dt_variables_tag::type>,
255 : boost::hash<Key>>*>
256 : boundary_data_history,
257 : const gsl::not_null<
258 : std::unordered_map<Key, TimeStepId, boost::hash<Key>>*>
259 : mortar_next_time_step_ids,
260 : const gsl::not_null<DirectionalIdMap<Dim, Mesh<Dim>>*>
261 : neighbor_mesh,
262 : const Element<Dim>& element) {
263 : // Remove neighbor meshes for neighbors that don't exist anymore
264 : domain::remove_nonexistent_neighbors(neighbor_mesh, element);
265 :
266 : // Move received boundary data into boundary history.
267 : for (auto& [mortar_id, mortar_next_time_step_id] :
268 : *mortar_next_time_step_ids) {
269 : if (mortar_id.id == ElementId<Dim>::external_boundary_id()) {
270 : continue;
271 : }
272 : while (needed_time(mortar_next_time_step_id)) {
273 : const auto time_entry = inbox.find(mortar_next_time_step_id);
274 : if (time_entry == inbox.end()) {
275 : return false;
276 : }
277 : const auto received_mortar_data =
278 : time_entry->second.find(mortar_id);
279 : if (received_mortar_data == time_entry->second.end()) {
280 : return false;
281 : }
282 :
283 : MortarData<Dim> neighbor_mortar_data{};
284 : // Insert:
285 : // - the current TimeStepId of the neighbor
286 : // - the current face mesh of the neighbor
287 : // - the current boundary correction data of the neighbor
288 : ASSERT(std::get<3>(received_mortar_data->second).has_value(),
289 : "Did not receive boundary correction data from the "
290 : "neighbor\nMortarId: "
291 : << mortar_id
292 : << "\nTimeStepId: " << mortar_next_time_step_id);
293 : neighbor_mesh->insert_or_assign(
294 : mortar_id, std::get<0>(received_mortar_data->second));
295 : neighbor_mortar_data.insert_neighbor_mortar_data(
296 : mortar_next_time_step_id,
297 : std::get<1>(received_mortar_data->second),
298 : std::move(*std::get<3>(received_mortar_data->second)));
299 : // We don't yet communicate the integration order, because
300 : // we don't have any variable-order methods. The
301 : // fixed-order methods ignore the field.
302 : boundary_data_history->at(mortar_id).remote().insert(
303 : mortar_next_time_step_id, std::numeric_limits<size_t>::max(),
304 : std::move(neighbor_mortar_data));
305 : mortar_next_time_step_id =
306 : std::get<4>(received_mortar_data->second);
307 : time_entry->second.erase(received_mortar_data);
308 : if (time_entry->second.empty()) {
309 : inbox.erase(time_entry);
310 : }
311 : }
312 : }
313 : return true;
314 : },
315 : box, db::get<::domain::Tags::Element<Dim>>(*box));
316 :
317 : return have_all_intermediate_messages;
318 : }
319 :
320 : /// Apply corrections from boundary communication.
321 : ///
322 : /// If `LocalTimeStepping` is false, updates the derivative of the variables,
323 : /// which should be done before taking a time step. If
324 : /// `LocalTimeStepping` is true, updates the variables themselves, which should
325 : /// be done after the volume update.
326 : ///
327 : /// Setting \p DenseOutput to true receives data required for output
328 : /// at ::Tags::Time instead of performing a full step. This is only
329 : /// used for local time-stepping.
330 : template <bool LocalTimeStepping, typename System, size_t VolumeDim,
331 : bool DenseOutput>
332 1 : struct ApplyBoundaryCorrections {
333 0 : static constexpr bool local_time_stepping = LocalTimeStepping;
334 : static_assert(local_time_stepping or not DenseOutput,
335 : "GTS does not use ApplyBoundaryCorrections for dense output.");
336 :
337 0 : using system = System;
338 0 : static constexpr size_t volume_dim = VolumeDim;
339 0 : using variables_tag = typename system::variables_tag;
340 0 : using dt_variables_tag = db::add_tag_prefix<::Tags::dt, variables_tag>;
341 0 : using DtVariables = typename dt_variables_tag::type;
342 0 : using volume_tags_for_dg_boundary_terms =
343 : tmpl::remove_duplicates<tmpl::flatten<tmpl::transform<
344 : typename system::boundary_correction_base::creatable_classes,
345 : detail::get_dg_boundary_terms<tmpl::_1>>>>;
346 :
347 0 : using TimeStepperType =
348 : tmpl::conditional_t<local_time_stepping, LtsTimeStepper, TimeStepper>;
349 :
350 0 : using tag_to_update =
351 : tmpl::conditional_t<local_time_stepping, variables_tag, dt_variables_tag>;
352 0 : using mortar_data_tag = tmpl::conditional_t<
353 : local_time_stepping,
354 : evolution::dg::Tags::MortarDataHistory<volume_dim, DtVariables>,
355 : evolution::dg::Tags::MortarData<volume_dim>>;
356 0 : using MortarDataType =
357 : tmpl::conditional_t<DenseOutput, const typename mortar_data_tag::type,
358 : typename mortar_data_tag::type>;
359 :
360 0 : using return_tags =
361 : tmpl::conditional_t<DenseOutput, tmpl::list<tag_to_update>,
362 : tmpl::list<tag_to_update, mortar_data_tag>>;
363 0 : using argument_tags = tmpl::append<
364 : tmpl::flatten<tmpl::list<
365 : tmpl::conditional_t<DenseOutput, mortar_data_tag, tmpl::list<>>,
366 : domain::Tags::Mesh<volume_dim>, Tags::MortarMesh<volume_dim>,
367 : Tags::MortarSize<volume_dim>, ::dg::Tags::Formulation,
368 : evolution::dg::Tags::NormalCovectorAndMagnitude<volume_dim>,
369 : ::Tags::TimeStepper<TimeStepperType>,
370 : evolution::Tags::BoundaryCorrection<system>,
371 : tmpl::conditional_t<DenseOutput, ::Tags::Time, ::Tags::TimeStep>,
372 : tmpl::conditional_t<local_time_stepping, tmpl::list<>,
373 : domain::Tags::DetInvJacobian<
374 : Frame::ElementLogical, Frame::Inertial>>>>,
375 : volume_tags_for_dg_boundary_terms>;
376 :
377 : // full step
378 : template <typename... VolumeArgs>
379 0 : static void apply(
380 : const gsl::not_null<typename tag_to_update::type*> vars_to_update,
381 : const gsl::not_null<MortarDataType*> mortar_data,
382 : const Mesh<volume_dim>& volume_mesh,
383 : const typename Tags::MortarMesh<volume_dim>::type& mortar_meshes,
384 : const typename Tags::MortarSize<volume_dim>::type& mortar_sizes,
385 : const ::dg::Formulation dg_formulation,
386 : const DirectionMap<
387 : volume_dim, std::optional<Variables<tmpl::list<
388 : evolution::dg::Tags::MagnitudeOfNormal,
389 : evolution::dg::Tags::NormalCovector<volume_dim>>>>>&
390 : face_normal_covector_and_magnitude,
391 : const TimeStepperType& time_stepper,
392 : const typename system::boundary_correction_base& boundary_correction,
393 : const TimeDelta& time_step,
394 : const Scalar<DataVector>& gts_det_inv_jacobian,
395 : const VolumeArgs&... volume_args) {
396 : apply_impl(vars_to_update, mortar_data, volume_mesh, mortar_meshes,
397 : mortar_sizes, dg_formulation, face_normal_covector_and_magnitude,
398 : time_stepper, boundary_correction, time_step,
399 : std::numeric_limits<double>::signaling_NaN(),
400 : gts_det_inv_jacobian, volume_args...);
401 : }
402 :
403 : template <typename... VolumeArgs>
404 0 : static void apply(
405 : const gsl::not_null<typename tag_to_update::type*> vars_to_update,
406 : const gsl::not_null<MortarDataType*> mortar_data,
407 : const Mesh<volume_dim>& volume_mesh,
408 : const typename Tags::MortarMesh<volume_dim>::type& mortar_meshes,
409 : const typename Tags::MortarSize<volume_dim>::type& mortar_sizes,
410 : const ::dg::Formulation dg_formulation,
411 : const DirectionMap<
412 : volume_dim, std::optional<Variables<tmpl::list<
413 : evolution::dg::Tags::MagnitudeOfNormal,
414 : evolution::dg::Tags::NormalCovector<volume_dim>>>>>&
415 : face_normal_covector_and_magnitude,
416 : const TimeStepperType& time_stepper,
417 : const typename system::boundary_correction_base& boundary_correction,
418 : const TimeDelta& time_step, const VolumeArgs&... volume_args) {
419 : apply_impl(vars_to_update, mortar_data, volume_mesh, mortar_meshes,
420 : mortar_sizes, dg_formulation, face_normal_covector_and_magnitude,
421 : time_stepper, boundary_correction, time_step,
422 : std::numeric_limits<double>::signaling_NaN(), {},
423 : volume_args...);
424 : }
425 :
426 : // dense output (LTS only)
427 : template <typename... VolumeArgs>
428 0 : static void apply(
429 : const gsl::not_null<typename variables_tag::type*> vars_to_update,
430 : const MortarDataType& mortar_data, const Mesh<volume_dim>& volume_mesh,
431 : const typename Tags::MortarMesh<volume_dim>::type& mortar_meshes,
432 : const typename Tags::MortarSize<volume_dim>::type& mortar_sizes,
433 : const ::dg::Formulation dg_formulation,
434 : const DirectionMap<
435 : volume_dim, std::optional<Variables<tmpl::list<
436 : evolution::dg::Tags::MagnitudeOfNormal,
437 : evolution::dg::Tags::NormalCovector<volume_dim>>>>>&
438 : face_normal_covector_and_magnitude,
439 : const LtsTimeStepper& time_stepper,
440 : const typename system::boundary_correction_base& boundary_correction,
441 : const double dense_output_time, const VolumeArgs&... volume_args) {
442 : apply_impl(vars_to_update, &mortar_data, volume_mesh, mortar_meshes,
443 : mortar_sizes, dg_formulation, face_normal_covector_and_magnitude,
444 : time_stepper, boundary_correction, TimeDelta{},
445 : dense_output_time, {}, volume_args...);
446 : }
447 :
448 : template <typename DbTagsList, typename... InboxTags, typename Metavariables,
449 : typename ArrayIndex, typename ParallelComponent>
450 0 : static bool is_ready(
451 : const gsl::not_null<db::DataBox<DbTagsList>*> box,
452 : const gsl::not_null<tuples::TaggedTuple<InboxTags...>*> inboxes,
453 : Parallel::GlobalCache<Metavariables>& /*cache*/,
454 : const ArrayIndex& /*array_index*/,
455 : const ParallelComponent* const /*component*/) {
456 : if constexpr (local_time_stepping) {
457 : return receive_boundary_data_local_time_stepping<System, VolumeDim,
458 : DenseOutput>(box,
459 : inboxes);
460 : } else {
461 : return receive_boundary_data_global_time_stepping<Metavariables>(box,
462 : inboxes);
463 : }
464 : }
465 :
466 : private:
467 : template <typename... VolumeArgs>
468 0 : static void apply_impl(
469 : const gsl::not_null<typename tag_to_update::type*> vars_to_update,
470 : const gsl::not_null<MortarDataType*> mortar_data,
471 : const Mesh<volume_dim>& volume_mesh,
472 : const typename Tags::MortarMesh<volume_dim>::type& mortar_meshes,
473 : const typename Tags::MortarSize<volume_dim>::type& mortar_sizes,
474 : const ::dg::Formulation dg_formulation,
475 : const DirectionMap<
476 : volume_dim, std::optional<Variables<tmpl::list<
477 : evolution::dg::Tags::MagnitudeOfNormal,
478 : evolution::dg::Tags::NormalCovector<volume_dim>>>>>&
479 : face_normal_covector_and_magnitude,
480 : const TimeStepperType& time_stepper,
481 : const typename system::boundary_correction_base& boundary_correction,
482 : const TimeDelta& time_step, const double dense_output_time,
483 : const Scalar<DataVector>& gts_det_inv_jacobian,
484 : const VolumeArgs&... volume_args) {
485 : tuples::tagged_tuple_from_typelist<db::wrap_tags_in<
486 : detail::TemporaryReference, volume_tags_for_dg_boundary_terms>>
487 : volume_args_tuple{volume_args...};
488 :
489 : // Set up helper lambda that will compute and lift the boundary corrections
490 : ASSERT(
491 : volume_mesh.quadrature() ==
492 : make_array<volume_dim>(volume_mesh.quadrature(0)),
493 : "Must have isotropic quadrature, but got volume mesh: " << volume_mesh);
494 : const bool using_gauss_lobatto_points =
495 : volume_mesh.quadrature(0) == Spectral::Quadrature::GaussLobatto;
496 :
497 : Scalar<DataVector> volume_det_inv_jacobian{};
498 : Scalar<DataVector> volume_det_jacobian{};
499 : if constexpr (not local_time_stepping) {
500 : if (not using_gauss_lobatto_points) {
501 : get(volume_det_inv_jacobian)
502 : .set_data_ref(make_not_null(
503 : // NOLINTNEXTLINE(cppcoreguidelines-pro-type-const-cast)
504 : &const_cast<DataVector&>(get(gts_det_inv_jacobian))));
505 : get(volume_det_jacobian) = 1.0 / get(volume_det_inv_jacobian);
506 : }
507 : }
508 :
509 : using derived_boundary_corrections =
510 : typename std::decay_t<decltype(boundary_correction)>::creatable_classes;
511 :
512 : static_assert(
513 : tmpl::all<derived_boundary_corrections, std::is_final<tmpl::_1>>::value,
514 : "All createable classes for boundary corrections must be marked "
515 : "final.");
516 : call_with_dynamic_type<void, derived_boundary_corrections>(
517 : &boundary_correction,
518 : [&dense_output_time, &dg_formulation,
519 : &face_normal_covector_and_magnitude, &mortar_data, &mortar_meshes,
520 : &mortar_sizes, &time_step, &time_stepper, using_gauss_lobatto_points,
521 : &vars_to_update, &volume_args_tuple, &volume_det_jacobian,
522 : &volume_det_inv_jacobian,
523 : &volume_mesh](auto* typed_boundary_correction) {
524 : using BcType = std::decay_t<decltype(*typed_boundary_correction)>;
525 : // Compute internal boundary quantities on the mortar for sides of
526 : // the element that have neighbors, i.e. they are not an external
527 : // side.
528 : using mortar_tags_list = typename BcType::dg_package_field_tags;
529 :
530 : // Variables for reusing allocations. The actual values are
531 : // not reused.
532 : DtVariables dt_boundary_correction_on_mortar{};
533 : DtVariables volume_dt_correction{};
534 : // These variables may change size for each mortar and require
535 : // a new memory allocation, but they may also happen to need
536 : // to be the same size twice in a row, in which case holding
537 : // on to the allocation is a win.
538 : Scalar<DataVector> face_det_jacobian{};
539 : Variables<mortar_tags_list> local_data_on_mortar{};
540 : Variables<mortar_tags_list> neighbor_data_on_mortar{};
541 :
542 : for (auto& mortar_id_and_data : *mortar_data) {
543 : const auto& mortar_id = mortar_id_and_data.first;
544 : const auto& direction = mortar_id.direction;
545 : if (UNLIKELY(mortar_id.id ==
546 : ElementId<volume_dim>::external_boundary_id())) {
547 : ERROR(
548 : "Cannot impose boundary conditions on external boundary in "
549 : "direction "
550 : << direction
551 : << " in the ApplyBoundaryCorrections action. Boundary "
552 : "conditions are applied in the ComputeTimeDerivative "
553 : "action "
554 : "instead. You may have unintentionally added external "
555 : "mortars in one of the initialization actions.");
556 : }
557 :
558 : const Mesh<volume_dim - 1> face_mesh =
559 : volume_mesh.slice_away(direction.dimension());
560 :
561 : const auto compute_correction_coupling =
562 : [&typed_boundary_correction, &direction, dg_formulation,
563 : &dt_boundary_correction_on_mortar, &face_det_jacobian,
564 : &face_mesh, &face_normal_covector_and_magnitude,
565 : &local_data_on_mortar, &mortar_id, &mortar_meshes,
566 : &mortar_sizes, &neighbor_data_on_mortar,
567 : using_gauss_lobatto_points, &volume_args_tuple,
568 : &volume_det_jacobian, &volume_det_inv_jacobian,
569 : &volume_dt_correction, &volume_mesh](
570 : const MortarData<volume_dim>& local_mortar_data,
571 : const MortarData<volume_dim>& neighbor_mortar_data)
572 : -> DtVariables {
573 : if (local_time_stepping and not using_gauss_lobatto_points) {
574 : // This needs to be updated every call because the Jacobian
575 : // may be time-dependent. In the case of time-independent maps
576 : // and local time stepping we could first perform the integral
577 : // on the boundaries, and then lift to the volume. This is
578 : // left as a future optimization.
579 : local_mortar_data.get_local_volume_det_inv_jacobian(
580 : make_not_null(&volume_det_inv_jacobian));
581 : get(volume_det_jacobian) = 1.0 / get(volume_det_inv_jacobian);
582 : }
583 : const auto& mortar_mesh = mortar_meshes.at(mortar_id);
584 :
585 : // Extract local and neighbor data, copy into Variables because
586 : // we store them in a std::vector for type erasure.
587 : const std::pair<Mesh<volume_dim - 1>, DataVector>&
588 : local_mesh_and_data = *local_mortar_data.local_mortar_data();
589 : const std::pair<Mesh<volume_dim - 1>, DataVector>&
590 : neighbor_mesh_and_data =
591 : *neighbor_mortar_data.neighbor_mortar_data();
592 : local_data_on_mortar.set_data_ref(
593 : const_cast<double*>(std::get<1>(local_mesh_and_data).data()),
594 : std::get<1>(local_mesh_and_data).size());
595 : neighbor_data_on_mortar.set_data_ref(
596 : const_cast<double*>(
597 : std::get<1>(neighbor_mesh_and_data).data()),
598 : std::get<1>(neighbor_mesh_and_data).size());
599 :
600 : // The boundary computations and lifting can be further
601 : // optimized by in the h-refinement case having only one
602 : // allocation for the face and having the projection from the
603 : // mortar to the face be done in place. E.g.
604 : // local_data_on_mortar and neighbor_data_on_mortar could be
605 : // allocated fewer times, as well as `needs_projection` section
606 : // below could do an in-place projection.
607 : dt_boundary_correction_on_mortar.initialize(
608 : mortar_mesh.number_of_grid_points());
609 :
610 : call_boundary_correction(
611 : make_not_null(&dt_boundary_correction_on_mortar),
612 : local_data_on_mortar, neighbor_data_on_mortar,
613 : *typed_boundary_correction, dg_formulation, volume_args_tuple,
614 : typename BcType::dg_boundary_terms_volume_tags{});
615 :
616 : const std::array<Spectral::MortarSize, volume_dim - 1>&
617 : mortar_size = mortar_sizes.at(mortar_id);
618 :
619 : // This cannot reuse an allocation because it is initialized
620 : // via move-assignment. (If it is used at all.)
621 : DtVariables dt_boundary_correction_projected_onto_face{};
622 : auto& dt_boundary_correction =
623 : [&dt_boundary_correction_on_mortar,
624 : &dt_boundary_correction_projected_onto_face, &face_mesh,
625 : &mortar_mesh, &mortar_size]() -> DtVariables& {
626 : if (Spectral::needs_projection(face_mesh, mortar_mesh,
627 : mortar_size)) {
628 : dt_boundary_correction_projected_onto_face =
629 : ::dg::project_from_mortar(
630 : dt_boundary_correction_on_mortar, face_mesh,
631 : mortar_mesh, mortar_size);
632 : return dt_boundary_correction_projected_onto_face;
633 : }
634 : return dt_boundary_correction_on_mortar;
635 : }();
636 :
637 : // Both paths initialize this to be non-owning.
638 : Scalar<DataVector> magnitude_of_face_normal{};
639 : if constexpr (local_time_stepping) {
640 : (void)face_normal_covector_and_magnitude;
641 : local_mortar_data.get_local_face_normal_magnitude(
642 : &magnitude_of_face_normal);
643 : } else {
644 : ASSERT(
645 : face_normal_covector_and_magnitude.count(direction) == 1 and
646 : face_normal_covector_and_magnitude.at(direction)
647 : .has_value(),
648 : "Face normal covector and magnitude not set in "
649 : "direction: "
650 : << direction);
651 : get(magnitude_of_face_normal)
652 : .set_data_ref(make_not_null(&const_cast<DataVector&>(
653 : get(get<evolution::dg::Tags::MagnitudeOfNormal>(
654 : *face_normal_covector_and_magnitude.at(
655 : direction))))));
656 : }
657 :
658 : if (using_gauss_lobatto_points) {
659 : // The lift_flux function lifts only on the slice, it does not
660 : // add the contribution to the volume.
661 : ::dg::lift_flux(make_not_null(&dt_boundary_correction),
662 : volume_mesh.extents(direction.dimension()),
663 : magnitude_of_face_normal);
664 : return std::move(dt_boundary_correction);
665 : } else {
666 : // We are using Gauss points.
667 : //
668 : // Notes:
669 : // - We should really lift both sides simultaneously since this
670 : // reduces memory accesses. Lifting all sides at the same
671 : // time is unlikely to improve performance since we lift by
672 : // jumping through slices. There may also be compatibility
673 : // issues with local time stepping.
674 : // - If we lift both sides at the same time we first need to
675 : // deal with projecting from mortars to the face, then lift
676 : // off the faces. With non-owning Variables memory
677 : // allocations could be significantly reduced in this code.
678 : if constexpr (local_time_stepping) {
679 : ASSERT(get(volume_det_inv_jacobian).size() > 0,
680 : "For local time stepping the volume determinant of "
681 : "the inverse Jacobian has not been set.");
682 :
683 : local_mortar_data.get_local_face_det_jacobian(
684 : make_not_null(&face_det_jacobian));
685 : } else {
686 : // Project the determinant of the Jacobian to the face. This
687 : // could be optimized by caching in the time-independent case.
688 : get(face_det_jacobian)
689 : .destructive_resize(face_mesh.number_of_grid_points());
690 : const Matrix identity{};
691 : auto interpolation_matrices =
692 : make_array<volume_dim>(std::cref(identity));
693 : const std::pair<Matrix, Matrix>& matrices =
694 : Spectral::boundary_interpolation_matrices(
695 : volume_mesh.slice_through(direction.dimension()));
696 : gsl::at(interpolation_matrices, direction.dimension()) =
697 : direction.side() == Side::Upper ? matrices.second
698 : : matrices.first;
699 : apply_matrices(make_not_null(&get(face_det_jacobian)),
700 : interpolation_matrices,
701 : get(volume_det_jacobian),
702 : volume_mesh.extents());
703 : }
704 :
705 : volume_dt_correction.initialize(
706 : volume_mesh.number_of_grid_points(), 0.0);
707 : ::dg::lift_boundary_terms_gauss_points(
708 : make_not_null(&volume_dt_correction),
709 : volume_det_inv_jacobian, volume_mesh, direction,
710 : dt_boundary_correction, magnitude_of_face_normal,
711 : face_det_jacobian);
712 : return std::move(volume_dt_correction);
713 : }
714 : };
715 :
716 : if constexpr (local_time_stepping) {
717 : typename variables_tag::type lgl_lifted_data{};
718 : auto& lifted_data = using_gauss_lobatto_points ? lgl_lifted_data
719 : : *vars_to_update;
720 : if (using_gauss_lobatto_points) {
721 : lifted_data.initialize(face_mesh.number_of_grid_points(), 0.0);
722 : }
723 :
724 : auto& mortar_data_history = mortar_id_and_data.second;
725 : if constexpr (DenseOutput) {
726 : (void)time_step;
727 : time_stepper.boundary_dense_output(
728 : &lifted_data, mortar_data_history, dense_output_time,
729 : compute_correction_coupling);
730 : } else {
731 : (void)dense_output_time;
732 : time_stepper.add_boundary_delta(
733 : &lifted_data, make_not_null(&mortar_data_history),
734 : time_step, compute_correction_coupling);
735 : }
736 :
737 : if (using_gauss_lobatto_points) {
738 : // Add the flux contribution to the volume data
739 : add_slice_to_data(
740 : vars_to_update, lifted_data, volume_mesh.extents(),
741 : direction.dimension(),
742 : index_to_slice_at(volume_mesh.extents(), direction));
743 : }
744 : } else {
745 : (void)time_step;
746 : (void)time_stepper;
747 : (void)dense_output_time;
748 :
749 : // Choose an allocation cache that may be empty, so we
750 : // might be able to reuse the allocation obtained for the
751 : // lifted data. This may result in a self assignment,
752 : // depending on the code paths taken, but handling the
753 : // results this way makes the GTS and LTS paths more
754 : // similar because the LTS code always stores the result
755 : // in the history and so sometimes benefits from moving
756 : // into the return value of compute_correction_coupling.
757 : auto& lifted_data = using_gauss_lobatto_points
758 : ? dt_boundary_correction_on_mortar
759 : : volume_dt_correction;
760 : lifted_data = compute_correction_coupling(
761 : mortar_id_and_data.second, mortar_id_and_data.second);
762 : // Remove data since it's tagged with the time. In the future we
763 : // _might_ be able to reuse allocations, but this optimization
764 : // should only be done after profiling.
765 : mortar_id_and_data.second.extract();
766 :
767 : if (using_gauss_lobatto_points) {
768 : // Add the flux contribution to the volume data
769 : add_slice_to_data(
770 : vars_to_update, lifted_data, volume_mesh.extents(),
771 : direction.dimension(),
772 : index_to_slice_at(volume_mesh.extents(), direction));
773 : } else {
774 : *vars_to_update += lifted_data;
775 : }
776 : }
777 : }
778 : });
779 : }
780 :
781 : template <typename... BoundaryCorrectionTags, typename... Tags,
782 : typename BoundaryCorrection, typename... AllVolumeArgs,
783 : typename... VolumeTagsForCorrection>
784 0 : static void call_boundary_correction(
785 : const gsl::not_null<Variables<tmpl::list<BoundaryCorrectionTags...>>*>
786 : boundary_corrections_on_mortar,
787 : const Variables<tmpl::list<Tags...>>& local_boundary_data,
788 : const Variables<tmpl::list<Tags...>>& neighbor_boundary_data,
789 : const BoundaryCorrection& boundary_correction,
790 : const ::dg::Formulation dg_formulation,
791 : const tuples::TaggedTuple<detail::TemporaryReference<AllVolumeArgs>...>&
792 : volume_args_tuple,
793 : tmpl::list<VolumeTagsForCorrection...> /*meta*/) {
794 : boundary_correction.dg_boundary_terms(
795 : make_not_null(
796 : &get<BoundaryCorrectionTags>(*boundary_corrections_on_mortar))...,
797 : get<Tags>(local_boundary_data)..., get<Tags>(neighbor_boundary_data)...,
798 : dg_formulation,
799 : tuples::get<detail::TemporaryReference<VolumeTagsForCorrection>>(
800 : volume_args_tuple)...);
801 : }
802 : };
803 :
804 1 : namespace Actions {
805 : /*!
806 : * \brief Computes the boundary corrections for global time-stepping
807 : * and adds them to the time derivative.
808 : */
809 : template <typename System, size_t VolumeDim, bool DenseOutput>
810 1 : struct ApplyBoundaryCorrectionsToTimeDerivative {
811 0 : using inbox_tags = tmpl::list<
812 : evolution::dg::Tags::BoundaryCorrectionAndGhostCellsInbox<VolumeDim>>;
813 0 : using const_global_cache_tags =
814 : tmpl::list<evolution::Tags::BoundaryCorrection<System>,
815 : ::dg::Tags::Formulation>;
816 :
817 : template <typename DbTagsList, typename... InboxTags, typename Metavariables,
818 : typename ArrayIndex, typename ActionList,
819 : typename ParallelComponent>
820 0 : static Parallel::iterable_action_return_t apply(
821 : db::DataBox<DbTagsList>& box, tuples::TaggedTuple<InboxTags...>& inboxes,
822 : const Parallel::GlobalCache<Metavariables>& /*cache*/,
823 : const ArrayIndex& /*array_index*/, ActionList /*meta*/,
824 : const ParallelComponent* const /*meta*/) {
825 : constexpr size_t volume_dim = Metavariables::system::volume_dim;
826 : const Element<volume_dim>& element =
827 : db::get<domain::Tags::Element<volume_dim>>(box);
828 :
829 : if (UNLIKELY(element.number_of_neighbors() == 0)) {
830 : // We have no neighbors, yay!
831 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
832 : }
833 :
834 : if (not receive_boundary_data_global_time_stepping<Metavariables>(
835 : make_not_null(&box), make_not_null(&inboxes))) {
836 : return {Parallel::AlgorithmExecution::Retry, std::nullopt};
837 : }
838 :
839 : db::mutate_apply<
840 : ApplyBoundaryCorrections<false, System, VolumeDim, DenseOutput>>(
841 : make_not_null(&box));
842 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
843 : }
844 : };
845 :
846 : /*!
847 : * \brief Computes the boundary corrections for local time-stepping
848 : * and adds them to the variables.
849 : *
850 : * When using local time stepping the neighbor sends data at the neighbor's
851 : * current temporal id. Along with the boundary data, the next temporal id at
852 : * which the neighbor will send data is also sent. This is equal to the
853 : * neighbor's `::Tags::Next<::Tags::TimeStepId>`. When inserting into the mortar
854 : * data history, we insert the received temporal id, that is, the current time
855 : * of the neighbor, along with the boundary correction data.
856 : */
857 : template <typename System, size_t VolumeDim, bool DenseOutput>
858 1 : struct ApplyLtsBoundaryCorrections {
859 0 : using inbox_tags = tmpl::list<
860 : evolution::dg::Tags::BoundaryCorrectionAndGhostCellsInbox<VolumeDim>>;
861 0 : using const_global_cache_tags =
862 : tmpl::list<evolution::Tags::BoundaryCorrection<System>,
863 : ::dg::Tags::Formulation>;
864 :
865 : template <typename DbTagsList, typename... InboxTags, typename Metavariables,
866 : typename ArrayIndex, typename ActionList,
867 : typename ParallelComponent>
868 0 : static Parallel::iterable_action_return_t apply(
869 : db::DataBox<DbTagsList>& box, tuples::TaggedTuple<InboxTags...>& inboxes,
870 : const Parallel::GlobalCache<Metavariables>& /*cache*/,
871 : const ArrayIndex& /*array_index*/, ActionList /*meta*/,
872 : const ParallelComponent* const /*meta*/) {
873 : constexpr size_t volume_dim = Metavariables::system::volume_dim;
874 : const Element<volume_dim>& element =
875 : db::get<domain::Tags::Element<volume_dim>>(box);
876 :
877 : if (UNLIKELY(element.number_of_neighbors() == 0)) {
878 : // We have no neighbors, yay!
879 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
880 : }
881 :
882 : if (not receive_boundary_data_local_time_stepping<System, VolumeDim, false>(
883 : make_not_null(&box), make_not_null(&inboxes))) {
884 : return {Parallel::AlgorithmExecution::Retry, std::nullopt};
885 : }
886 :
887 : db::mutate_apply<
888 : ApplyBoundaryCorrections<true, System, VolumeDim, DenseOutput>>(
889 : make_not_null(&box));
890 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
891 : }
892 : };
893 : } // namespace Actions
894 : } // namespace evolution::dg
|
