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