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