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