Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : #pragma once
5 :
6 : #include <algorithm>
7 : #include <array>
8 : #include <cstddef>
9 : #include <iterator>
10 : #include <limits>
11 : #include <map>
12 : #include <optional>
13 : #include <tuple>
14 : #include <unordered_set>
15 : #include <utility>
16 :
17 : #include "DataStructures/DataBox/DataBox.hpp"
18 : #include "DataStructures/DataBox/Prefixes.hpp"
19 : #include "DataStructures/DataVector.hpp"
20 : #include "DataStructures/Index.hpp"
21 : #include "DataStructures/Tensor/Tensor.hpp"
22 : #include "DataStructures/Variables.hpp"
23 : #include "DataStructures/VariablesTag.hpp"
24 : #include "Domain/Structure/Direction.hpp"
25 : #include "Domain/Structure/DirectionalId.hpp"
26 : #include "Domain/Structure/DirectionalIdMap.hpp"
27 : #include "Domain/Structure/Element.hpp"
28 : #include "Domain/Structure/ElementId.hpp"
29 : #include "Domain/Structure/OrientationMapHelpers.hpp"
30 : #include "Domain/Structure/TrimMap.hpp"
31 : #include "Domain/Tags.hpp"
32 : #include "Domain/Tags/NeighborMesh.hpp"
33 : #include "Evolution/DgSubcell/ActiveGrid.hpp"
34 : #include "Evolution/DgSubcell/CombineVolumeGhostData.hpp"
35 : #include "Evolution/DgSubcell/GhostData.hpp"
36 : #include "Evolution/DgSubcell/NeighborRdmpAndVolumeData.hpp"
37 : #include "Evolution/DgSubcell/Projection.hpp"
38 : #include "Evolution/DgSubcell/RdmpTci.hpp"
39 : #include "Evolution/DgSubcell/RdmpTciData.hpp"
40 : #include "Evolution/DgSubcell/SliceData.hpp"
41 : #include "Evolution/DgSubcell/SubcellOptions.hpp"
42 : #include "Evolution/DgSubcell/Tags/ActiveGrid.hpp"
43 : #include "Evolution/DgSubcell/Tags/CellCenteredFlux.hpp"
44 : #include "Evolution/DgSubcell/Tags/DataForRdmpTci.hpp"
45 : #include "Evolution/DgSubcell/Tags/GhostDataForReconstruction.hpp"
46 : #include "Evolution/DgSubcell/Tags/Interpolators.hpp"
47 : #include "Evolution/DgSubcell/Tags/Mesh.hpp"
48 : #include "Evolution/DgSubcell/Tags/MeshForGhostData.hpp"
49 : #include "Evolution/DgSubcell/Tags/SubcellOptions.hpp"
50 : #include "Evolution/DgSubcell/Tags/TciStatus.hpp"
51 : #include "Evolution/DiscontinuousGalerkin/BoundaryData.hpp"
52 : #include "Evolution/DiscontinuousGalerkin/InboxTags.hpp"
53 : #include "Evolution/DiscontinuousGalerkin/MortarData.hpp"
54 : #include "Evolution/DiscontinuousGalerkin/MortarDataHolder.hpp"
55 : #include "Evolution/DiscontinuousGalerkin/MortarTags.hpp"
56 : #include "Evolution/DiscontinuousGalerkin/TimeSteppingPolicy.hpp"
57 : #include "NumericalAlgorithms/Interpolation/IrregularInterpolant.hpp"
58 : #include "NumericalAlgorithms/Spectral/Mesh.hpp"
59 : #include "Parallel/AlgorithmExecution.hpp"
60 : #include "Parallel/ArrayCollection/IsDgElementCollection.hpp"
61 : #include "Parallel/GlobalCache.hpp"
62 : #include "Time/Tags/HistoryEvolvedVariables.hpp"
63 : #include "Time/TimeStepId.hpp"
64 : #include "Utilities/Algorithm.hpp"
65 : #include "Utilities/ErrorHandling/Assert.hpp"
66 : #include "Utilities/Gsl.hpp"
67 : #include "Utilities/Literals.hpp"
68 : #include "Utilities/MakeArray.hpp"
69 : #include "Utilities/TMPL.hpp"
70 :
71 : /// \cond
72 : namespace Tags {
73 : struct TimeStepId;
74 : } // namespace Tags
75 : namespace evolution::dg::Tags {
76 : template <size_t Dim>
77 : struct MortarInfo;
78 : } // namespace evolution::dg::Tags
79 : /// \endcond
80 :
81 : namespace evolution::dg::subcell::Actions {
82 : /*!
83 : * \brief Sets the local data from the relaxed discrete maximum principle
84 : * troubled-cell indicator and sends ghost zone data to neighboring elements.
85 : *
86 : * The action proceeds as follows:
87 : *
88 : * 1. Determine in which directions we have neighbors
89 : * 2. Slice the variables provided by GhostDataMutator to send to our neighbors
90 : * for ghost zones
91 : * 3. Send the ghost zone data, appending the max/min for the TCI at the end of
92 : * the `DataVector` we are sending.
93 : *
94 : * \warning This assumes the RDMP TCI data in the DataBox has been set, it does
95 : * not calculate it automatically. The reason is this way we can only calculate
96 : * the RDMP data when it's needed since computing it can be pretty expensive.
97 : *
98 : * Some notes:
99 : * - In the future we will need to send the cell-centered fluxes to do
100 : * high-order FD without additional reconstruction being necessary.
101 : *
102 : * GlobalCache:
103 : * - Uses:
104 : * - `ParallelComponent` proxy
105 : *
106 : * DataBox:
107 : * - Uses:
108 : * - `domain::Tags::Mesh<Dim>`
109 : * - `subcell::Tags::Mesh<Dim>`
110 : * - `domain::Tags::Element<Dim>`
111 : * - `Tags::TimeStepId`
112 : * - `Tags::Next<Tags::TimeStepId>`
113 : * - `subcell::Tags::ActiveGrid`
114 : * - `System::variables_tag`
115 : * - `subcell::Tags::DataForRdmpTci`
116 : * - Adds: nothing
117 : * - Removes: nothing
118 : * - Modifies:
119 : * - `subcell::Tags::GhostDataForReconstruction<Dim>`
120 : */
121 : template <size_t Dim, typename GhostDataMutator, bool UseNodegroupDgElements>
122 1 : struct SendDataForReconstruction {
123 0 : using inbox_tags =
124 : tmpl::list<evolution::dg::Tags::BoundaryCorrectionAndGhostCellsInbox<
125 : Dim, UseNodegroupDgElements>>;
126 :
127 : template <typename DbTags, typename... InboxTags, typename ArrayIndex,
128 : typename ActionList, typename ParallelComponent,
129 : typename Metavariables>
130 0 : static Parallel::iterable_action_return_t apply(
131 : db::DataBox<DbTags>& box, tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
132 : Parallel::GlobalCache<Metavariables>& cache,
133 : const ArrayIndex& /*array_index*/, const ActionList /*meta*/,
134 : const ParallelComponent* const /*meta*/) {
135 : static_assert(UseNodegroupDgElements ==
136 : Parallel::is_dg_element_collection_v<ParallelComponent>,
137 : "The action SendDataForReconstruction is told by the "
138 : "template parameter UseNodegroupDgElements that it is being "
139 : "used with a DgElementCollection, but the ParallelComponent "
140 : "is not a DgElementCollection. You need to change the "
141 : "template parameter on the SendDataForReconstruction action "
142 : "in your action list.");
143 :
144 : ASSERT(db::get<Tags::ActiveGrid>(box) == ActiveGrid::Subcell,
145 : "The SendDataForReconstruction action can only be called when "
146 : "Subcell is the active scheme.");
147 :
148 : const Element<Dim>& element = db::get<::domain::Tags::Element<Dim>>(box);
149 : ASSERT(alg::all_of(db::get<evolution::dg::Tags::MortarInfo<Dim>>(box),
150 : [](const auto& mortar) {
151 : return mortar.second.time_stepping_policy() ==
152 : TimeSteppingPolicy::EqualRate;
153 : }),
154 : "Cannot send subcell data from "
155 : << element.id() << " across an LTS mortar: "
156 : << db::get<evolution::dg::Tags::MortarInfo<Dim>>(box));
157 :
158 : using flux_variables = typename Metavariables::system::flux_variables;
159 :
160 : db::mutate<Tags::GhostDataForReconstruction<Dim>>(
161 : [](const auto ghost_data_ptr) {
162 : // Clear the previous neighbor data and add current local data
163 : ghost_data_ptr->clear();
164 : },
165 : make_not_null(&box));
166 :
167 : const Mesh<Dim>& dg_mesh = db::get<::domain::Tags::Mesh<Dim>>(box);
168 : const Mesh<Dim>& subcell_mesh = db::get<Tags::Mesh<Dim>>(box);
169 : const size_t ghost_zone_size =
170 : Metavariables::SubcellOptions::ghost_zone_size(box);
171 :
172 : // Optimization note: could save a copy+allocation if we moved
173 : // all_sliced_data when possible before sending.
174 : //
175 : // Note: RDMP size doesn't help here since we need to slice data after
176 : // anyway, so no way to save an allocation through that.
177 : const auto& cell_centered_flux =
178 : db::get<Tags::CellCenteredFlux<flux_variables, Dim>>(box);
179 : DataVector volume_data_to_slice = db::mutate_apply(
180 : GhostDataMutator{}, make_not_null(&box),
181 : cell_centered_flux.has_value() ? cell_centered_flux.value().size()
182 : : 0_st);
183 : if (cell_centered_flux.has_value()) {
184 : std::copy(
185 : cell_centered_flux.value().data(),
186 : std::next(
187 : cell_centered_flux.value().data(),
188 : static_cast<std::ptrdiff_t>(cell_centered_flux.value().size())),
189 : std::next(
190 : volume_data_to_slice.data(),
191 : static_cast<std::ptrdiff_t>(volume_data_to_slice.size() -
192 : cell_centered_flux.value().size())));
193 : }
194 :
195 : // When using enable_extension_directions, we send the ghost data
196 : // for problematic directions separately with the new action
197 : // ReceiveAndSendDataForReconstruction, so we only slice the data
198 : // for non-problematic directions here. (see the documentation of
199 : // ReceiveAndSendDataForReconstruction for what "problematic" means)
200 : const auto& extension_directions =
201 : db::get<evolution::dg::subcell::Tags::ExtensionDirections<Dim>>(box);
202 : std::unordered_set<Direction<Dim>> directions_to_work;
203 : for (const auto& internal_direction : element.internal_boundaries()) {
204 : if (extension_directions.contains(internal_direction)) {
205 : continue;
206 : } else {
207 : directions_to_work.insert(internal_direction);
208 : }
209 : }
210 :
211 : const DirectionMap<Dim, DataVector> all_sliced_data = slice_data(
212 : volume_data_to_slice, subcell_mesh.extents(), ghost_zone_size,
213 : directions_to_work, 0,
214 : db::get<
215 : evolution::dg::subcell::Tags::InterpolatorsFromFdToNeighborFd<Dim>>(
216 : box));
217 :
218 : auto& receiver_proxy =
219 : Parallel::get_parallel_component<ParallelComponent>(cache);
220 : const RdmpTciData& rdmp_tci_data = db::get<Tags::DataForRdmpTci>(box);
221 : const TimeStepId& time_step_id = db::get<::Tags::TimeStepId>(box);
222 : const TimeStepId& next_time_step_id =
223 : db::get<::Tags::Next<::Tags::TimeStepId>>(box);
224 :
225 : const int tci_decision =
226 : db::get<evolution::dg::subcell::Tags::TciDecision>(box);
227 : using history_tags = ::Tags::get_all_history_tags<DbTags>;
228 : static_assert(tmpl::size<history_tags>::value == 1);
229 : const auto& integration_order =
230 : db::get<tmpl::front<history_tags>>(box).integration_order();
231 : // Compute and send actual variables
232 : for (const auto& [direction, neighbors_in_direction] :
233 : element.neighbors()) {
234 : // Only need to send data for directions that are not
235 : // problematic directions (keys of extension_directions).
236 : if (not extension_directions.contains(direction)) {
237 : ASSERT(neighbors_in_direction.size() == 1,
238 : "AMR is not yet supported when using DG-subcell. Note that this "
239 : "condition could be relaxed to support AMR only where the "
240 : "evolution is using DG without any changes to subcell.");
241 :
242 : for (const ElementId<Dim>& neighbor : neighbors_in_direction) {
243 : const auto& orientation =
244 : neighbors_in_direction.orientation(neighbor);
245 : const auto direction_from_neighbor =
246 : orientation(direction.opposite());
247 : const size_t rdmp_size = rdmp_tci_data.max_variables_values.size() +
248 : rdmp_tci_data.min_variables_values.size();
249 : const auto& sliced_data_in_direction = all_sliced_data.at(direction);
250 :
251 : // Allocate with subcell data and rdmp data
252 : DataVector subcell_data_to_send{sliced_data_in_direction.size() +
253 : rdmp_size};
254 : // Note: Currently we interpolate our solution to our neighbor FD grid
255 : // even when grid points align but are oriented differently. There's a
256 : // possible optimization for the rare (almost never?) edge case where
257 : // two blocks have the same ghost zone coordinates but have different
258 : // orientations (e.g. RotatedBricks). Since this shouldn't ever happen
259 : // outside of tests, we currently don't bother with it. If we wanted
260 : // to, here's the code:
261 : //
262 : // if (not orientation.is_aligned()) {
263 : // std::array<size_t, Dim> slice_extents{};
264 : // for (size_t d = 0; d < Dim; ++d) {
265 : // gsl::at(slice_extents, d) = subcell_mesh.extents(d);
266 : // }
267 : // gsl::at(slice_extents, direction.dimension()) = ghost_zone_size;
268 : // // Need a view so we only get the subcell data and not the rdmp
269 : // // data
270 : // DataVector subcell_data_to_send_view{
271 : // subcell_data_to_send.data(),
272 : // subcell_data_to_send.size() - rdmp_size};
273 : // orient_variables(make_not_null(&subcell_data_to_send_view),
274 : // sliced_data_in_direction,
275 : // Index<Dim>{slice_extents}, orientation);
276 : // } else { std::copy(...); }
277 : //
278 :
279 : // Copy over data since it's already oriented from interpolation
280 : std::copy(sliced_data_in_direction.begin(),
281 : sliced_data_in_direction.end(),
282 : subcell_data_to_send.begin());
283 : // Copy rdmp data to end of subcell_data_to_send
284 : std::copy(rdmp_tci_data.max_variables_values.cbegin(),
285 : rdmp_tci_data.max_variables_values.cend(),
286 : std::prev(subcell_data_to_send.end(),
287 : static_cast<int>(rdmp_size)));
288 : std::copy(rdmp_tci_data.min_variables_values.cbegin(),
289 : rdmp_tci_data.min_variables_values.cend(),
290 : std::prev(subcell_data_to_send.end(),
291 : static_cast<int>(
292 : rdmp_tci_data.min_variables_values.size())));
293 :
294 : evolution::dg::BoundaryData<Dim> data{
295 : dg_mesh, subcell_mesh,
296 : std::nullopt, std::move(subcell_data_to_send),
297 : std::nullopt, next_time_step_id,
298 : tci_decision, integration_order};
299 :
300 : Parallel::receive_data<
301 : evolution::dg::Tags::BoundaryCorrectionAndGhostCellsInbox<
302 : Dim,
303 : Parallel::is_dg_element_collection_v<ParallelComponent>>>(
304 : receiver_proxy[neighbor], time_step_id,
305 : std::pair{
306 : DirectionalId<Dim>{direction_from_neighbor, element.id()},
307 : std::move(data)});
308 : }
309 : }
310 : }
311 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
312 : }
313 : };
314 :
315 : /*!
316 : * \brief Handles ghost data at block boundaries to ensure interpolation
317 : * is (almost) always used instead of extrapolation.
318 : *
319 : * Since the coordinate maps are only continuous and not smooth at block
320 : * boundaries, the logical coordinate axes, and therefore grid point axes,
321 : * do not necessarily align. An example is given in the image below, which is
322 : * a snapshot of 2d circular domain built from one central square surrounded
323 : * by four deformed wedges. We zoom in on the upper right corner of the
324 : * cube for illustration purposes.
325 : *
326 : * \image html curved_mesh_illustration.png width=600px
327 : *
328 : * Blue circles denote the cell-centered FD points in the two elements whose
329 : * ghost points are being exchanged, bright red diamonds denote the ghost
330 : * points needed for reconstruction in the element on the right, and dark red
331 : * squares denote the cell-centered FD points in a neighboring element not
332 : * directly participating in the exchange. The dashed blue and dash-dotted
333 : * red lines show lines of constant logical coordinates in the left and
334 : * right elements, respectively. Notice that they intersect on the boundary,
335 : * but do not align.
336 : *
337 : * In this example, three lowest ghost points cannot be directly
338 : * interpolated from the element on the left. In such a case,
339 : * we flag the direction (from the perspective of the element on the left,
340 : * in this example, the direction the green arrow points to) as
341 : * problematic and the ghost points may be filled either
342 : * by extrapolation or interpolation. Extrapolation can lead to an unphysical
343 : * state like negative densities, so interpolation is generally preferred.
344 : * To enable interpolation, set `EnableExtensionDirections` to true in
345 : * `SubcellOptions` part of the input file. This enables the use of
346 : * ghost data from the neighboring element in the extension direction
347 : * (in this example, the direction the grey arrow points to) to fill the
348 : * ghost points.
349 : *
350 : * \warning This option is only available in the case where we are only
351 : * using FD scheme for the evolution, i.e. we are using true for
352 : * `AlwaysUseSubcell` in the `SubcellOptions` part of the input file.
353 : * Currently, the following cases are not supported:
354 : * 1. There are multiple neighbors in a direction
355 : * 2. There are multiple extension directions required for a
356 : * single problematic direction.
357 : * 3. The extension direction is itself a problematic direction, which can
358 : * result in a deadlock.
359 : *
360 : * When disabled, this action does nothing.
361 : *
362 : * When enabled, the action:
363 : * 1. Receives subcell ghost data from neighbors for all non-problematic
364 : * directions.
365 : * 2. For “problematic directions” it extends the element’s volume data
366 : * using ghost data from another neighbor (in the "extension direction").
367 : * This extension ensures that the ghost data can be filled using
368 : * interpolation.
369 : * 3. Interpolates to the requested ghost points, and then sends the completed
370 : * ghost data to the original neighbor.
371 : */
372 : template <size_t Dim, typename GhostDataMutator, bool UseNodegroupDgElements>
373 1 : struct ReceiveAndSendDataForReconstruction {
374 0 : using inbox_tags =
375 : tmpl::list<evolution::dg::Tags::BoundaryCorrectionAndGhostCellsInbox<
376 : Dim, UseNodegroupDgElements>>;
377 : template <typename DbTags, typename... InboxTags, typename ArrayIndex,
378 : typename ActionList, typename ParallelComponent,
379 : typename Metavariables>
380 0 : static Parallel::iterable_action_return_t apply(
381 : db::DataBox<DbTags>& box, tuples::TaggedTuple<InboxTags...>& inboxes,
382 : Parallel::GlobalCache<Metavariables>& cache,
383 : const ArrayIndex& /*array_index*/, const ActionList /*meta*/,
384 : const ParallelComponent* const /*meta*/) {
385 : if (not db::get<evolution::dg::subcell::Tags::SubcellOptions<Dim>>(box)
386 : .enable_extension_directions()) {
387 : // We are not using extension directions, so just continue without doing
388 : // any work.
389 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
390 : }
391 :
392 : const auto& extension_directions =
393 : db::get<evolution::dg::subcell::Tags::ExtensionDirections<Dim>>(box);
394 : if (extension_directions.empty()) {
395 : // For this element, we have no extension directions, so just
396 : // continue without doing any work.
397 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
398 : }
399 :
400 : const Element<Dim>& element = db::get<::domain::Tags::Element<Dim>>(box);
401 : // Need to subtract number of problematic directions.
402 : const auto number_of_expected_messages =
403 : element.neighbors().size() - extension_directions.size();
404 :
405 : if (UNLIKELY(number_of_expected_messages == 0)) {
406 : // We have no neighbors, so just continue without doing any work.
407 : // Technically, this could also happen if all of the neighbors are in
408 : // problematic directions, but this case is unlikely to ever happen.
409 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
410 : }
411 :
412 : std::unordered_set<DirectionalId<Dim>> expected_keys;
413 : for (const auto& internal_direction : element.internal_boundaries()) {
414 : if (not extension_directions.contains(internal_direction)) {
415 : // Note here, we assume that we have only one neighbor per
416 : // direction, so we can just take the first one.
417 : ASSERT(element.neighbors().at(internal_direction).ids().size() == 1,
418 : "Assumption one neighbor per direction failed. "
419 : "direction: "
420 : << internal_direction << ", neighbors.size() = "
421 : << element.neighbors().at(internal_direction).ids().size()
422 : << ", element: " << element.id());
423 : expected_keys.emplace(
424 : internal_direction,
425 : *element.neighbors().at(internal_direction).ids().begin());
426 : }
427 : }
428 :
429 : const auto& interpolants = db::get<
430 : evolution::dg::subcell::Tags::InterpolatorsFromFdToNeighborFd<Dim>>(
431 : box);
432 :
433 : const auto& current_time_step_id = db::get<::Tags::TimeStepId>(box);
434 : auto& inbox =
435 : tuples::get<evolution::dg::Tags::BoundaryCorrectionAndGhostCellsInbox<
436 : Metavariables::volume_dim,
437 : Parallel::is_dg_element_collection_v<ParallelComponent>>>(inboxes);
438 : inbox.collect_messages();
439 : const auto received = inbox.messages.find(current_time_step_id);
440 : // Check we have at least some data from correct time, and then check
441 : // we have received all data
442 : if (received == inbox.messages.end()) {
443 : inbox.set_missing_messages(expected_keys.size());
444 : return {Parallel::AlgorithmExecution::Retry, std::nullopt};
445 : }
446 : if (const auto missing = std::count_if(
447 : expected_keys.begin(), expected_keys.end(),
448 : [&received](const auto& key) {
449 : return received->second.find(key) == received->second.end();
450 : });
451 : missing != 0) {
452 : inbox.set_missing_messages(static_cast<size_t>(missing));
453 : return {Parallel::AlgorithmExecution::Retry, std::nullopt};
454 : }
455 :
456 : const size_t ghost_zone_size =
457 : Metavariables::SubcellOptions::ghost_zone_size(box);
458 : const Mesh<Dim>& dg_mesh = db::get<::domain::Tags::Mesh<Dim>>(box);
459 : const Mesh<Dim>& subcell_mesh = db::get<Tags::Mesh<Dim>>(box);
460 : const Index<Dim>& subcell_extents = subcell_mesh.extents();
461 :
462 : const auto& received_data = received->second;
463 : ASSERT(received_data.size() >= number_of_expected_messages,
464 : "received_data size: " << received_data.size()
465 : << " less than expected number of messages: "
466 : << number_of_expected_messages << " !");
467 :
468 : using flux_variables = typename Metavariables::system::flux_variables;
469 : const auto& cell_centered_flux =
470 : db::get<Tags::CellCenteredFlux<flux_variables, Dim>>(box);
471 : DataVector volume_data_to_slice = db::mutate_apply(
472 : GhostDataMutator{}, make_not_null(&box),
473 : cell_centered_flux.has_value() ? cell_centered_flux.value().size()
474 : : 0_st);
475 : if (cell_centered_flux.has_value()) {
476 : std::copy(
477 : cell_centered_flux.value().data(),
478 : std::next(
479 : cell_centered_flux.value().data(),
480 : static_cast<std::ptrdiff_t>(cell_centered_flux.value().size())),
481 : std::next(
482 : volume_data_to_slice.data(),
483 : static_cast<std::ptrdiff_t>(volume_data_to_slice.size() -
484 : cell_centered_flux.value().size())));
485 : }
486 :
487 : auto& receiver_proxy =
488 : Parallel::get_parallel_component<ParallelComponent>(cache);
489 : const RdmpTciData& rdmp_tci_data = db::get<Tags::DataForRdmpTci>(box);
490 : const TimeStepId& time_step_id = db::get<::Tags::TimeStepId>(box);
491 : const TimeStepId& next_time_step_id =
492 : db::get<::Tags::Next<::Tags::TimeStepId>>(box);
493 :
494 : const int tci_decision =
495 : db::get<evolution::dg::subcell::Tags::TciDecision>(box);
496 : using history_tags = ::Tags::get_all_history_tags<DbTags>;
497 : static_assert(tmpl::size<history_tags>::value == 1);
498 : const auto& integration_order =
499 : db::get<tmpl::front<history_tags>>(box).integration_order();
500 :
501 : const size_t number_of_points = subcell_mesh.extents().product();
502 : // Number of independent components per grid point.
503 : const size_t number_of_components =
504 : volume_data_to_slice.size() / number_of_points;
505 :
506 : // For each problematic direction (a direction for which the element
507 : // cannot directly provide ghost data):
508 : // 1. Receive ghost data from our neighbor in the direction specified by
509 : // extension_direction.direction_to_extend.
510 : // 2. Use this data to extend our own volume data, then interpolate to the
511 : // required ghost points.
512 : // 3. Send the final ghost data to the neighbor at problematic_direction.
513 : for (const auto& [problematic_direction, extension_direction] :
514 : extension_directions) {
515 : // Direction to extend the volume data.
516 : const Direction<Dim> direction_to_extend =
517 : extension_direction.direction_to_extend;
518 :
519 : // Check that direction_to_extend is not a problematic direction.
520 : ASSERT(problematic_direction != direction_to_extend,
521 : "The direction to extend must not be a problematic direction. "
522 : "problematic_direction: "
523 : << problematic_direction
524 : << ", direction_to_extend: " << direction_to_extend);
525 :
526 : // Only one neighbor per direction.
527 : ASSERT(element.neighbors().at(direction_to_extend).ids().size() == 1,
528 : "Assumption one neighbor per direction failed. "
529 : "direction_to_extend: "
530 : << direction_to_extend << ", neighbors.size() = "
531 : << element.neighbors().at(direction_to_extend).ids().size()
532 : << ", element: " << element.id());
533 :
534 : // Only internal boundary for the extension.
535 : ASSERT(element.internal_boundaries().contains(direction_to_extend),
536 : "Direction to extend not an internal boundary! "
537 : "dir_to_extend = "
538 : << direction_to_extend << ", element.internal_boundaries() = "
539 : << element.internal_boundaries()
540 : << ", element = " << element.id());
541 :
542 : // Again, we are assuming that we have only one neighbor per direction,
543 : // and we are just taking the first one.
544 : const ElementId<Dim> relevant_neighbor_id =
545 : *((element.neighbors()).at(direction_to_extend).ids().begin());
546 :
547 : // Received data must have entry for direction to extend.
548 : ASSERT((received_data.contains(DirectionalId<Dim>{direction_to_extend,
549 : relevant_neighbor_id})),
550 : "Received data missing entry for direction to extend."
551 : << " direction_to_extend = " << direction_to_extend
552 : << ", relevant_neighbor_id = " << relevant_neighbor_id
553 : << ", problematic_direction = " << problematic_direction
554 : << ", element = " << element.id());
555 :
556 : // Received data must have received ghost data for this extension
557 : // direction.
558 : ASSERT(((received_data.at(DirectionalId<Dim>{direction_to_extend,
559 : relevant_neighbor_id}))
560 : .ghost_cell_data.has_value()),
561 : "Ghost data missing for this extension direction."
562 : << " direction_to_extend = " << direction_to_extend
563 : << ", relevant_neighbor_id = " << relevant_neighbor_id
564 : << ", problematic_direction = " << problematic_direction
565 : << ", element = " << element.id());
566 :
567 : const ElementId<Dim> problematic_neighbor_id =
568 : *((element.neighbors()).at(problematic_direction).ids().begin());
569 :
570 : const auto& orientation = (element.neighbors())
571 : .at(problematic_direction)
572 : .orientation(problematic_neighbor_id);
573 : const auto direction_from_neighbor =
574 : orientation(problematic_direction.opposite());
575 :
576 : const size_t rdmp_size = rdmp_tci_data.max_variables_values.size() +
577 : rdmp_tci_data.min_variables_values.size();
578 :
579 : const DataVector& full_ghost_cell_data =
580 : received_data
581 : .at(DirectionalId<Dim>{direction_to_extend, relevant_neighbor_id})
582 : .ghost_cell_data.value();
583 : const size_t relevant_ghost_data_size =
584 : full_ghost_cell_data.size() - rdmp_size;
585 :
586 : const DataVector relevant_ghost_data;
587 : make_const_view(make_not_null(&relevant_ghost_data), full_ghost_cell_data,
588 : 0, relevant_ghost_data_size);
589 :
590 : const auto& interpolant =
591 : (interpolants.at(DirectionalId<Dim>{problematic_direction,
592 : problematic_neighbor_id}))
593 : .value();
594 : const DataVector combined_data = combine_volume_ghost_data(
595 : volume_data_to_slice, relevant_ghost_data, subcell_extents,
596 : ghost_zone_size, direction_to_extend);
597 : const size_t result_size =
598 : ghost_zone_size * subcell_mesh.extents()
599 : .slice_away(problematic_direction.dimension())
600 : .product();
601 : const size_t span_size = result_size * number_of_components;
602 :
603 : DataVector subcell_data_to_send{span_size + rdmp_size};
604 :
605 : auto result_span = gsl::make_span(subcell_data_to_send.data(), span_size);
606 : interpolant.interpolate(
607 : make_not_null(&result_span),
608 : gsl::make_span(combined_data.data(), combined_data.size()));
609 :
610 : // Copy rdmp data to end of subcell_data_to_send
611 : std::copy(
612 : rdmp_tci_data.max_variables_values.cbegin(),
613 : rdmp_tci_data.max_variables_values.cend(),
614 : std::prev(subcell_data_to_send.end(), static_cast<int>(rdmp_size)));
615 : std::copy(rdmp_tci_data.min_variables_values.cbegin(),
616 : rdmp_tci_data.min_variables_values.cend(),
617 : std::prev(subcell_data_to_send.end(),
618 : static_cast<int>(
619 : rdmp_tci_data.min_variables_values.size())));
620 : evolution::dg::BoundaryData<Dim> data{
621 : dg_mesh, subcell_mesh,
622 : std::nullopt, std::move(subcell_data_to_send),
623 : std::nullopt, next_time_step_id,
624 : tci_decision, integration_order};
625 : Parallel::receive_data<
626 : evolution::dg::Tags::BoundaryCorrectionAndGhostCellsInbox<
627 : Dim, Parallel::is_dg_element_collection_v<ParallelComponent>>>(
628 : receiver_proxy[problematic_neighbor_id], time_step_id,
629 : std::pair{DirectionalId<Dim>{direction_from_neighbor, element.id()},
630 : std::move(data)});
631 : }
632 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
633 : }
634 : };
635 : /*!
636 : * \brief Receive the subcell data from our neighbor, and accumulate the data
637 : * from the relaxed discrete maximum principle troubled-cell indicator.
638 : *
639 : * Note:
640 : * - Since we only care about the min/max over all neighbors and ourself at the
641 : * past time, we accumulate all data immediately into the `RdmpTciData`.
642 : * - If the neighbor is using DG and therefore sends boundary correction data
643 : * then that is added into the `evolution::dg::Tags::MortarData` tag
644 : * - The next `TimeStepId` is recorded, but we do not yet support local time
645 : * stepping.
646 : * - This action will never care about what variables are sent for
647 : * reconstruction. It is only responsible for receiving the data and storing
648 : * it in the `NeighborData`.
649 : *
650 : * GlobalCache:
651 : * -Uses: nothing
652 : *
653 : * DataBox:
654 : * - Uses:
655 : * - `domain::Tags::Element<Dim>`
656 : * - `Tags::TimeStepId`
657 : * - `domain::Tags::Mesh<Dim>`
658 : * - `subcell::Tags::Mesh<Dim>`
659 : * - `domain::Tags::Element<Dim>`
660 : * - `Tags::Next<Tags::TimeStepId>`
661 : * - `subcell::Tags::ActiveGrid`
662 : * - `System::variables_tag`
663 : * - Adds: nothing
664 : * - Removes: nothing
665 : * - Modifies:
666 : * - `subcell::Tags::GhostDataForReconstruction<Dim>`
667 : * - `subcell::Tags::DataForRdmpTci`
668 : * - `evolution::dg::Tags::MortarData`
669 : * - `evolution::dg::Tags::MortarNextTemporalId`
670 : */
671 : template <size_t Dim>
672 1 : struct ReceiveDataForReconstruction {
673 : template <typename DbTags, typename... InboxTags, typename ArrayIndex,
674 : typename ActionList, typename ParallelComponent,
675 : typename Metavariables>
676 0 : static Parallel::iterable_action_return_t apply(
677 : db::DataBox<DbTags>& box, tuples::TaggedTuple<InboxTags...>& inboxes,
678 : const Parallel::GlobalCache<Metavariables>& /*cache*/,
679 : const ArrayIndex& /*array_index*/, const ActionList /*meta*/,
680 : const ParallelComponent* const /*meta*/) {
681 : const Element<Dim>& element = db::get<::domain::Tags::Element<Dim>>(box);
682 : const auto number_of_expected_messages = element.neighbors().size();
683 : if (UNLIKELY(number_of_expected_messages == 0)) {
684 : // We have no neighbors, so just continue without doing any work
685 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
686 : }
687 :
688 : using ::operator<<;
689 : const auto& current_time_step_id = db::get<::Tags::TimeStepId>(box);
690 : auto& inbox =
691 : tuples::get<evolution::dg::Tags::BoundaryCorrectionAndGhostCellsInbox<
692 : Metavariables::volume_dim,
693 : Parallel::is_dg_element_collection_v<ParallelComponent>>>(inboxes);
694 : inbox.collect_messages();
695 : const auto received = inbox.messages.find(current_time_step_id);
696 : // Check we have at least some data from correct time, and then check that
697 : // we have received all data
698 : if (received == inbox.messages.end()) {
699 : inbox.set_missing_messages(number_of_expected_messages);
700 : return {Parallel::AlgorithmExecution::Retry, std::nullopt};
701 : }
702 : if (received->second.size() != number_of_expected_messages) {
703 : inbox.set_missing_messages(number_of_expected_messages -
704 : received->second.size());
705 : return {Parallel::AlgorithmExecution::Retry, std::nullopt};
706 : }
707 :
708 : // Now that we have received all the data, copy it over as needed.
709 : DirectionalIdMap<Dim, evolution::dg::BoundaryData<Dim>> received_data =
710 : std::move(received->second);
711 : inbox.messages.erase(received);
712 :
713 : const Mesh<Dim>& subcell_mesh = db::get<Tags::Mesh<Dim>>(box);
714 : const auto& mortar_meshes = get<evolution::dg::Tags::MortarMesh<Dim>>(box);
715 :
716 : db::mutate<Tags::GhostDataForReconstruction<Dim>, Tags::DataForRdmpTci,
717 : evolution::dg::Tags::MortarData<Dim>,
718 : evolution::dg::Tags::MortarNextTemporalId<Dim>,
719 : domain::Tags::NeighborMesh<Dim>,
720 : evolution::dg::subcell::Tags::MeshForGhostData<Dim>,
721 : evolution::dg::subcell::Tags::NeighborTciDecisions<Dim>>(
722 : [&element,
723 : ghost_zone_size = Metavariables::SubcellOptions::ghost_zone_size(box),
724 : &received_data, &subcell_mesh, &mortar_meshes](
725 : const gsl::not_null<DirectionalIdMap<Dim, GhostData>*>
726 : ghost_data_ptr,
727 : const gsl::not_null<RdmpTciData*> rdmp_tci_data_ptr,
728 : const gsl::not_null<
729 : DirectionalIdMap<Dim, evolution::dg::MortarDataHolder<Dim>>*>
730 : mortar_data,
731 : const gsl::not_null<DirectionalIdMap<Dim, TimeStepId>*>
732 : mortar_next_time_step_id,
733 : const gsl::not_null<DirectionalIdMap<Dim, Mesh<Dim>>*>
734 : neighbor_mesh,
735 : const gsl::not_null<DirectionalIdMap<Dim, Mesh<Dim>>*>
736 : mesh_for_ghost_data,
737 : const auto neighbor_tci_decisions,
738 : const DirectionalIdMap<Dim, std::optional<intrp::Irregular<Dim>>>&
739 : neighbor_dg_to_fd_interpolants) {
740 : // Remove neighbor meshes for neighbors that don't exist anymore
741 : domain::remove_nonexistent_neighbors(neighbor_mesh, element);
742 : domain::remove_nonexistent_neighbors(mesh_for_ghost_data, element);
743 :
744 : // Get the next time step id, and also the fluxes data if the neighbor
745 : // is doing DG.
746 : for (auto& received_mortar_data : received_data) {
747 : const auto& mortar_id = received_mortar_data.first;
748 : try {
749 : mortar_next_time_step_id->at(mortar_id) =
750 : received_mortar_data.second.validity_range;
751 : } catch (std::exception& e) {
752 : ERROR("Failed retrieving the MortarId: ("
753 : << mortar_id.direction() << ',' << mortar_id.id()
754 : << ") from the mortar_next_time_step_id. Got exception: "
755 : << e.what());
756 : }
757 : if (received_mortar_data.second.boundary_correction_data
758 : .has_value()) {
759 : mortar_data->at(mortar_id).neighbor().face_mesh =
760 : received_mortar_data.second.volume_mesh.slice_away(
761 : mortar_id.direction().dimension());
762 : mortar_data->at(mortar_id).neighbor().mortar_mesh =
763 : mortar_meshes.at(mortar_id);
764 : mortar_data->at(mortar_id).neighbor().mortar_data = std::move(
765 : *received_mortar_data.second.boundary_correction_data);
766 : }
767 : // Set new neighbor mesh
768 : neighbor_mesh->insert_or_assign(
769 : mortar_id, received_mortar_data.second.volume_mesh);
770 : mesh_for_ghost_data->insert_or_assign(
771 : mortar_id, received_mortar_data.second
772 : .volume_mesh_ghost_cell_data.value());
773 : }
774 :
775 : ASSERT(ghost_data_ptr->empty(),
776 : "Should have no elements in the neighbor data when "
777 : "receiving neighbor data");
778 : const size_t number_of_rdmp_vars =
779 : rdmp_tci_data_ptr->max_variables_values.size();
780 : ASSERT(rdmp_tci_data_ptr->min_variables_values.size() ==
781 : number_of_rdmp_vars,
782 : "The number of RDMP variables for which we have a maximum "
783 : "and minimum should be the same, but we have "
784 : << number_of_rdmp_vars << " for the max and "
785 : << rdmp_tci_data_ptr->min_variables_values.size()
786 : << " for the min.");
787 :
788 : for (const auto& [direction, neighbors_in_direction] :
789 : element.neighbors()) {
790 : for (const auto& neighbor : neighbors_in_direction) {
791 : DirectionalId<Dim> directional_element_id{direction, neighbor};
792 : ASSERT(ghost_data_ptr->count(directional_element_id) == 0,
793 : "Found neighbor already inserted in direction "
794 : << direction << " with ElementId " << neighbor);
795 : ASSERT(received_data[directional_element_id]
796 : .ghost_cell_data.has_value(),
797 : "Received subcell data message that does not contain any "
798 : "actual subcell data for reconstruction.");
799 : // Collect the max/min of u(t^n) for the RDMP as we receive data.
800 : // This reduces the memory footprint.
801 :
802 : evolution::dg::subcell::insert_neighbor_rdmp_and_volume_data(
803 : rdmp_tci_data_ptr, ghost_data_ptr,
804 : *received_data[directional_element_id].ghost_cell_data,
805 : number_of_rdmp_vars, directional_element_id,
806 : mesh_for_ghost_data->at(directional_element_id), element,
807 : subcell_mesh, ghost_zone_size,
808 : neighbor_dg_to_fd_interpolants);
809 : ASSERT(neighbor_tci_decisions->contains(directional_element_id),
810 : "The NeighorTciDecisions should contain the neighbor ("
811 : << directional_element_id.direction() << ", "
812 : << directional_element_id.id() << ") but doesn't");
813 : neighbor_tci_decisions->at(directional_element_id) =
814 : received_data[directional_element_id].tci_status;
815 : }
816 : }
817 : },
818 : make_not_null(&box),
819 : db::get<
820 : evolution::dg::subcell::Tags::InterpolatorsFromNeighborDgToFd<Dim>>(
821 : box));
822 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
823 : }
824 : };
825 : } // namespace evolution::dg::subcell::Actions
|
