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