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