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 found = static_cast<size_t>(
447 : alg::count_if(received->second,
448 : [&expected_keys](const auto& message) {
449 : return expected_keys.contains(message.first);
450 : }));
451 : found < expected_keys.size()) {
452 : inbox.set_missing_messages(expected_keys.size() - found);
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 : const auto received_data_for_direction_it =
548 : alg::find_if(received_data, [&](const auto& entry) {
549 : return entry.first == DirectionalId<Dim>{direction_to_extend,
550 : relevant_neighbor_id};
551 : });
552 : // Received data must have entry for direction to extend.
553 : ASSERT(received_data_for_direction_it != received_data.end(),
554 : "Received data missing entry for direction to extend."
555 : << " direction_to_extend = " << direction_to_extend
556 : << ", relevant_neighbor_id = " << relevant_neighbor_id
557 : << ", problematic_direction = " << problematic_direction
558 : << ", element = " << element.id());
559 : const auto& received_data_for_direction =
560 : received_data_for_direction_it->second;
561 :
562 : // Received data must have received ghost data for this extension
563 : // direction.
564 : ASSERT(received_data_for_direction.ghost_cell_data.has_value(),
565 : "Ghost data missing for this extension direction."
566 : << " direction_to_extend = " << direction_to_extend
567 : << ", relevant_neighbor_id = " << relevant_neighbor_id
568 : << ", problematic_direction = " << problematic_direction
569 : << ", element = " << element.id());
570 :
571 : const ElementId<Dim> problematic_neighbor_id =
572 : *((element.neighbors()).at(problematic_direction).ids().begin());
573 :
574 : const auto& orientation = (element.neighbors())
575 : .at(problematic_direction)
576 : .orientation(problematic_neighbor_id);
577 : const auto direction_from_neighbor =
578 : orientation(problematic_direction.opposite());
579 :
580 : const size_t rdmp_size = rdmp_tci_data.max_variables_values.size() +
581 : rdmp_tci_data.min_variables_values.size();
582 :
583 : const DataVector& full_ghost_cell_data =
584 : received_data_for_direction.ghost_cell_data.value();
585 : const size_t relevant_ghost_data_size =
586 : full_ghost_cell_data.size() - rdmp_size;
587 :
588 : const DataVector relevant_ghost_data;
589 : make_const_view(make_not_null(&relevant_ghost_data), full_ghost_cell_data,
590 : 0, relevant_ghost_data_size);
591 :
592 : const auto& interpolant =
593 : (interpolants.at(DirectionalId<Dim>{problematic_direction,
594 : problematic_neighbor_id}))
595 : .value();
596 : const DataVector combined_data = combine_volume_ghost_data(
597 : volume_data_to_slice, relevant_ghost_data, subcell_extents,
598 : ghost_zone_size, direction_to_extend);
599 : const size_t result_size =
600 : ghost_zone_size * subcell_mesh.extents()
601 : .slice_away(problematic_direction.dimension())
602 : .product();
603 : const size_t span_size = result_size * number_of_components;
604 :
605 : DataVector subcell_data_to_send{span_size + rdmp_size};
606 :
607 : auto result_span = gsl::make_span(subcell_data_to_send.data(), span_size);
608 : interpolant.interpolate(
609 : make_not_null(&result_span),
610 : gsl::make_span(combined_data.data(), combined_data.size()));
611 :
612 : // Copy rdmp data to end of subcell_data_to_send
613 : std::copy(
614 : rdmp_tci_data.max_variables_values.cbegin(),
615 : rdmp_tci_data.max_variables_values.cend(),
616 : std::prev(subcell_data_to_send.end(), static_cast<int>(rdmp_size)));
617 : std::copy(rdmp_tci_data.min_variables_values.cbegin(),
618 : rdmp_tci_data.min_variables_values.cend(),
619 : std::prev(subcell_data_to_send.end(),
620 : static_cast<int>(
621 : rdmp_tci_data.min_variables_values.size())));
622 : evolution::dg::BoundaryData<Dim> data{
623 : dg_mesh, subcell_mesh,
624 : std::nullopt, std::move(subcell_data_to_send),
625 : std::nullopt, next_time_step_id,
626 : tci_decision, integration_order};
627 : Parallel::receive_data<
628 : evolution::dg::Tags::BoundaryCorrectionAndGhostCellsInbox<
629 : Dim, Parallel::is_dg_element_collection_v<ParallelComponent>>>(
630 : receiver_proxy[problematic_neighbor_id], time_step_id,
631 : std::pair{DirectionalId<Dim>{direction_from_neighbor, element.id()},
632 : std::move(data)});
633 : }
634 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
635 : }
636 : };
637 : /*!
638 : * \brief Receive the subcell data from our neighbor, and accumulate the data
639 : * from the relaxed discrete maximum principle troubled-cell indicator.
640 : *
641 : * Note:
642 : * - Since we only care about the min/max over all neighbors and ourself at the
643 : * past time, we accumulate all data immediately into the `RdmpTciData`.
644 : * - If the neighbor is using DG and therefore sends boundary correction data
645 : * then that is added into the `evolution::dg::Tags::MortarData` tag
646 : * - The next `TimeStepId` is recorded, but we do not yet support local time
647 : * stepping.
648 : * - This action will never care about what variables are sent for
649 : * reconstruction. It is only responsible for receiving the data and storing
650 : * it in the `NeighborData`.
651 : *
652 : * GlobalCache:
653 : * -Uses: nothing
654 : *
655 : * DataBox:
656 : * - Uses:
657 : * - `domain::Tags::Element<Dim>`
658 : * - `Tags::TimeStepId`
659 : * - `domain::Tags::Mesh<Dim>`
660 : * - `subcell::Tags::Mesh<Dim>`
661 : * - `domain::Tags::Element<Dim>`
662 : * - `Tags::Next<Tags::TimeStepId>`
663 : * - `subcell::Tags::ActiveGrid`
664 : * - `System::variables_tag`
665 : * - Adds: nothing
666 : * - Removes: nothing
667 : * - Modifies:
668 : * - `subcell::Tags::GhostDataForReconstruction<Dim>`
669 : * - `subcell::Tags::DataForRdmpTci`
670 : * - `evolution::dg::Tags::MortarData`
671 : * - `evolution::dg::Tags::MortarNextTemporalId`
672 : */
673 : template <size_t Dim>
674 1 : struct ReceiveDataForReconstruction {
675 : template <typename DbTags, typename... InboxTags, typename ArrayIndex,
676 : typename ActionList, typename ParallelComponent,
677 : typename Metavariables>
678 0 : static Parallel::iterable_action_return_t apply(
679 : db::DataBox<DbTags>& box, tuples::TaggedTuple<InboxTags...>& inboxes,
680 : const Parallel::GlobalCache<Metavariables>& /*cache*/,
681 : const ArrayIndex& /*array_index*/, const ActionList /*meta*/,
682 : const ParallelComponent* const /*meta*/) {
683 : const Element<Dim>& element = db::get<::domain::Tags::Element<Dim>>(box);
684 : const auto number_of_expected_messages = element.neighbors().size();
685 : if (UNLIKELY(number_of_expected_messages == 0)) {
686 : // We have no neighbors, so just continue without doing any work
687 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
688 : }
689 :
690 : using ::operator<<;
691 : const auto& current_time_step_id = db::get<::Tags::TimeStepId>(box);
692 : auto& inbox =
693 : tuples::get<evolution::dg::Tags::BoundaryCorrectionAndGhostCellsInbox<
694 : Metavariables::volume_dim,
695 : Parallel::is_dg_element_collection_v<ParallelComponent>>>(inboxes);
696 : inbox.collect_messages();
697 : const auto received = inbox.messages.find(current_time_step_id);
698 : // Check we have at least some data from correct time, and then check that
699 : // we have received all data
700 : if (received == inbox.messages.end()) {
701 : inbox.set_missing_messages(number_of_expected_messages);
702 : return {Parallel::AlgorithmExecution::Retry, std::nullopt};
703 : }
704 : if (received->second.size() != number_of_expected_messages) {
705 : inbox.set_missing_messages(number_of_expected_messages -
706 : received->second.size());
707 : return {Parallel::AlgorithmExecution::Retry, std::nullopt};
708 : }
709 :
710 : // Now that we have received all the data, copy it over as needed.
711 : auto received_data = std::move(received->second);
712 : inbox.messages.erase(received);
713 :
714 : const Mesh<Dim>& subcell_mesh = db::get<Tags::Mesh<Dim>>(box);
715 : const auto& mortar_meshes = get<evolution::dg::Tags::MortarMesh<Dim>>(box);
716 :
717 : db::mutate<Tags::GhostDataForReconstruction<Dim>, Tags::DataForRdmpTci,
718 : evolution::dg::Tags::MortarData<Dim>,
719 : evolution::dg::Tags::MortarNextTemporalId<Dim>,
720 : domain::Tags::NeighborMesh<Dim>,
721 : evolution::dg::subcell::Tags::MeshForGhostData<Dim>,
722 : evolution::dg::subcell::Tags::NeighborTciDecisions<Dim>>(
723 : [&element,
724 : ghost_zone_size = Metavariables::SubcellOptions::ghost_zone_size(box),
725 : &received_data, &subcell_mesh, &mortar_meshes](
726 : const gsl::not_null<DirectionalIdMap<Dim, GhostData>*>
727 : ghost_data_ptr,
728 : const gsl::not_null<RdmpTciData*> rdmp_tci_data_ptr,
729 : const gsl::not_null<
730 : DirectionalIdMap<Dim, evolution::dg::MortarDataHolder<Dim>>*>
731 : mortar_data,
732 : const gsl::not_null<DirectionalIdMap<Dim, TimeStepId>*>
733 : mortar_next_time_step_id,
734 : const gsl::not_null<DirectionalIdMap<Dim, Mesh<Dim>>*>
735 : neighbor_mesh,
736 : const gsl::not_null<DirectionalIdMap<Dim, Mesh<Dim>>*>
737 : mesh_for_ghost_data,
738 : const auto neighbor_tci_decisions,
739 : const DirectionalIdMap<Dim, std::optional<intrp::Irregular<Dim>>>&
740 : neighbor_dg_to_fd_interpolants) {
741 : // Remove neighbor meshes for neighbors that don't exist anymore
742 : domain::remove_nonexistent_neighbors(neighbor_mesh, element);
743 : domain::remove_nonexistent_neighbors(mesh_for_ghost_data, element);
744 :
745 : // Get the next time step id, and also the fluxes data if the neighbor
746 : // is doing DG.
747 : for (auto& received_mortar_data : received_data) {
748 : const auto& mortar_id = received_mortar_data.first;
749 : try {
750 : mortar_next_time_step_id->at(mortar_id) =
751 : received_mortar_data.second.validity_range;
752 : } catch (std::exception& e) {
753 : ERROR("Failed retrieving the MortarId: ("
754 : << mortar_id.direction() << ',' << mortar_id.id()
755 : << ") from the mortar_next_time_step_id. Got exception: "
756 : << e.what());
757 : }
758 : if (received_mortar_data.second.boundary_correction_data
759 : .has_value()) {
760 : mortar_data->at(mortar_id).neighbor().face_mesh =
761 : received_mortar_data.second.volume_mesh.slice_away(
762 : mortar_id.direction().dimension());
763 : mortar_data->at(mortar_id).neighbor().mortar_mesh =
764 : mortar_meshes.at(mortar_id);
765 : mortar_data->at(mortar_id).neighbor().mortar_data = std::move(
766 : *received_mortar_data.second.boundary_correction_data);
767 : }
768 : // Set new neighbor mesh
769 : neighbor_mesh->insert_or_assign(
770 : mortar_id, received_mortar_data.second.volume_mesh);
771 : mesh_for_ghost_data->insert_or_assign(
772 : mortar_id, received_mortar_data.second
773 : .volume_mesh_ghost_cell_data.value());
774 : }
775 :
776 : ASSERT(ghost_data_ptr->empty(),
777 : "Should have no elements in the neighbor data when "
778 : "receiving neighbor data");
779 : const size_t number_of_rdmp_vars =
780 : rdmp_tci_data_ptr->max_variables_values.size();
781 : ASSERT(rdmp_tci_data_ptr->min_variables_values.size() ==
782 : number_of_rdmp_vars,
783 : "The number of RDMP variables for which we have a maximum "
784 : "and minimum should be the same, but we have "
785 : << number_of_rdmp_vars << " for the max and "
786 : << rdmp_tci_data_ptr->min_variables_values.size()
787 : << " for the min.");
788 :
789 : for (const auto& [directional_element_id, boundary_data] :
790 : received_data) {
791 : ASSERT(ghost_data_ptr->count(directional_element_id) == 0,
792 : "Found neighbor already inserted in direction "
793 : << directional_element_id.direction()
794 : << " with ElementId " << directional_element_id.id());
795 : ASSERT(boundary_data.ghost_cell_data.has_value(),
796 : "Received subcell data message that does not contain any "
797 : "actual subcell data for reconstruction.");
798 : // Collect the max/min of u(t^n) for the RDMP as we receive data.
799 : // This reduces the memory footprint.
800 :
801 : evolution::dg::subcell::insert_neighbor_rdmp_and_volume_data(
802 : rdmp_tci_data_ptr, ghost_data_ptr,
803 : *boundary_data.ghost_cell_data, number_of_rdmp_vars,
804 : directional_element_id,
805 : mesh_for_ghost_data->at(directional_element_id), element,
806 : subcell_mesh, ghost_zone_size, neighbor_dg_to_fd_interpolants);
807 : ASSERT(neighbor_tci_decisions->contains(directional_element_id),
808 : "The NeighorTciDecisions should contain the neighbor ("
809 : << directional_element_id.direction() << ", "
810 : << directional_element_id.id() << ") but doesn't");
811 : neighbor_tci_decisions->at(directional_element_id) =
812 : boundary_data.tci_status;
813 : }
814 : },
815 : make_not_null(&box),
816 : db::get<
817 : evolution::dg::subcell::Tags::InterpolatorsFromNeighborDgToFd<Dim>>(
818 : box));
819 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
820 : }
821 : };
822 : } // namespace evolution::dg::subcell::Actions
|
