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