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