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