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