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