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/GhostData.hpp"
35 : #include "Evolution/DgSubcell/NeighborRdmpAndVolumeData.hpp"
36 : #include "Evolution/DgSubcell/Projection.hpp"
37 : #include "Evolution/DgSubcell/RdmpTci.hpp"
38 : #include "Evolution/DgSubcell/RdmpTciData.hpp"
39 : #include "Evolution/DgSubcell/SliceData.hpp"
40 : #include "Evolution/DgSubcell/Tags/ActiveGrid.hpp"
41 : #include "Evolution/DgSubcell/Tags/CellCenteredFlux.hpp"
42 : #include "Evolution/DgSubcell/Tags/DataForRdmpTci.hpp"
43 : #include "Evolution/DgSubcell/Tags/GhostDataForReconstruction.hpp"
44 : #include "Evolution/DgSubcell/Tags/Interpolators.hpp"
45 : #include "Evolution/DgSubcell/Tags/Mesh.hpp"
46 : #include "Evolution/DgSubcell/Tags/Reconstructor.hpp"
47 : #include "Evolution/DgSubcell/Tags/TciStatus.hpp"
48 : #include "Evolution/DiscontinuousGalerkin/InboxTags.hpp"
49 : #include "Evolution/DiscontinuousGalerkin/MortarData.hpp"
50 : #include "Evolution/DiscontinuousGalerkin/MortarTags.hpp"
51 : #include "NumericalAlgorithms/Interpolation/IrregularInterpolant.hpp"
52 : #include "NumericalAlgorithms/Spectral/Mesh.hpp"
53 : #include "NumericalAlgorithms/Spectral/Spectral.hpp"
54 : #include "Parallel/AlgorithmExecution.hpp"
55 : #include "Parallel/GlobalCache.hpp"
56 : #include "Time/TimeStepId.hpp"
57 : #include "Utilities/ErrorHandling/Assert.hpp"
58 : #include "Utilities/Gsl.hpp"
59 : #include "Utilities/Literals.hpp"
60 : #include "Utilities/MakeArray.hpp"
61 : #include "Utilities/TMPL.hpp"
62 :
63 : /// \cond
64 : namespace Tags {
65 : struct TimeStepId;
66 : } // namespace Tags
67 : /// \endcond
68 :
69 : namespace evolution::dg::subcell::Actions {
70 : /*!
71 : * \brief Sets the local data from the relaxed discrete maximum principle
72 : * troubled-cell indicator and sends ghost zone data to neighboring elements.
73 : *
74 : * The action proceeds as follows:
75 : *
76 : * 1. Determine in which directions we have neighbors
77 : * 2. Slice the variables provided by GhostDataMutator to send to our neighbors
78 : * for ghost zones
79 : * 3. Send the ghost zone data, appending the max/min for the TCI at the end of
80 : * the `DataVector` we are sending.
81 : *
82 : * \warning This assumes the RDMP TCI data in the DataBox has been set, it does
83 : * not calculate it automatically. The reason is this way we can only calculate
84 : * the RDMP data when it's needed since computing it can be pretty expensive.
85 : *
86 : * Some notes:
87 : * - In the future we will need to send the cell-centered fluxes to do
88 : * high-order FD without additional reconstruction being necessary.
89 : *
90 : * GlobalCache:
91 : * - Uses:
92 : * - `ParallelComponent` proxy
93 : *
94 : * DataBox:
95 : * - Uses:
96 : * - `domain::Tags::Mesh<Dim>`
97 : * - `subcell::Tags::Mesh<Dim>`
98 : * - `domain::Tags::Element<Dim>`
99 : * - `Tags::TimeStepId`
100 : * - `Tags::Next<Tags::TimeStepId>`
101 : * - `subcell::Tags::ActiveGrid`
102 : * - `System::variables_tag`
103 : * - `subcell::Tags::DataForRdmpTci`
104 : * - Adds: nothing
105 : * - Removes: nothing
106 : * - Modifies:
107 : * - `subcell::Tags::GhostDataForReconstruction<Dim>`
108 : */
109 : template <size_t Dim, typename GhostDataMutator, bool LocalTimeStepping>
110 1 : struct SendDataForReconstruction {
111 0 : using inbox_tags = tmpl::list<
112 : evolution::dg::Tags::BoundaryCorrectionAndGhostCellsInbox<Dim>>;
113 :
114 : template <typename DbTags, typename... InboxTags, typename ArrayIndex,
115 : typename ActionList, typename ParallelComponent,
116 : typename Metavariables>
117 0 : static Parallel::iterable_action_return_t apply(
118 : db::DataBox<DbTags>& box, tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
119 : Parallel::GlobalCache<Metavariables>& cache,
120 : const ArrayIndex& /*array_index*/, const ActionList /*meta*/,
121 : const ParallelComponent* const /*meta*/) {
122 : static_assert(
123 : not LocalTimeStepping,
124 : "DG-subcell does not yet support local time stepping. The "
125 : "reconstruction data must be sent using dense output sometimes, and "
126 : "not at all other times. However, the data for the RDMP TCI should be "
127 : "sent along with the data for reconstruction each time.");
128 :
129 : ASSERT(db::get<Tags::ActiveGrid>(box) == ActiveGrid::Subcell,
130 : "The SendDataForReconstruction action can only be called when "
131 : "Subcell is the active scheme.");
132 : using flux_variables = typename Metavariables::system::flux_variables;
133 :
134 : db::mutate<Tags::GhostDataForReconstruction<Dim>>(
135 : [](const auto ghost_data_ptr) {
136 : // Clear the previous neighbor data and add current local data
137 : ghost_data_ptr->clear();
138 : },
139 : make_not_null(&box));
140 :
141 : const Mesh<Dim>& dg_mesh = db::get<::domain::Tags::Mesh<Dim>>(box);
142 : const Mesh<Dim>& subcell_mesh = db::get<Tags::Mesh<Dim>>(box);
143 : const Element<Dim>& element = db::get<::domain::Tags::Element<Dim>>(box);
144 : const size_t ghost_zone_size =
145 : db::get<evolution::dg::subcell::Tags::Reconstructor>(box)
146 : .ghost_zone_size();
147 :
148 : // Optimization note: could save a copy+allocation if we moved
149 : // all_sliced_data when possible before sending.
150 : //
151 : // Note: RDMP size doesn't help here since we need to slice data after
152 : // anyway, so no way to save an allocation through that.
153 : const auto& cell_centered_flux =
154 : db::get<Tags::CellCenteredFlux<flux_variables, Dim>>(box);
155 : DataVector volume_data_to_slice = db::mutate_apply(
156 : GhostDataMutator{}, make_not_null(&box),
157 : cell_centered_flux.has_value() ? cell_centered_flux.value().size()
158 : : 0_st);
159 : if (cell_centered_flux.has_value()) {
160 : std::copy(
161 : cell_centered_flux.value().data(),
162 : std::next(
163 : cell_centered_flux.value().data(),
164 : static_cast<std::ptrdiff_t>(cell_centered_flux.value().size())),
165 : std::next(
166 : volume_data_to_slice.data(),
167 : static_cast<std::ptrdiff_t>(volume_data_to_slice.size() -
168 : cell_centered_flux.value().size())));
169 : }
170 : const DirectionMap<Dim, DataVector> all_sliced_data = slice_data(
171 : volume_data_to_slice, subcell_mesh.extents(), ghost_zone_size,
172 : element.internal_boundaries(), 0,
173 : db::get<
174 : evolution::dg::subcell::Tags::InterpolatorsFromFdToNeighborFd<Dim>>(
175 : box));
176 :
177 : auto& receiver_proxy =
178 : Parallel::get_parallel_component<ParallelComponent>(cache);
179 : const RdmpTciData& rdmp_tci_data = db::get<Tags::DataForRdmpTci>(box);
180 : const TimeStepId& time_step_id = db::get<::Tags::TimeStepId>(box);
181 : const TimeStepId& next_time_step_id = [&box]() {
182 : if (LocalTimeStepping) {
183 : return db::get<::Tags::Next<::Tags::TimeStepId>>(box);
184 : } else {
185 : return db::get<::Tags::TimeStepId>(box);
186 : }
187 : }();
188 :
189 : const int tci_decision =
190 : db::get<evolution::dg::subcell::Tags::TciDecision>(box);
191 : // Compute and send actual variables
192 : for (const auto& [direction, neighbors_in_direction] :
193 : element.neighbors()) {
194 : const auto& orientation = neighbors_in_direction.orientation();
195 : const auto direction_from_neighbor = orientation(direction.opposite());
196 : ASSERT(neighbors_in_direction.size() == 1,
197 : "AMR is not yet supported when using DG-subcell. Note that this "
198 : "condition could be relaxed to support AMR only where the "
199 : "evolution is using DG without any changes to subcell.");
200 :
201 : for (const ElementId<Dim>& neighbor : neighbors_in_direction) {
202 : const size_t rdmp_size = rdmp_tci_data.max_variables_values.size() +
203 : rdmp_tci_data.min_variables_values.size();
204 : const auto& sliced_data_in_direction = all_sliced_data.at(direction);
205 : // Allocate with subcell data and rdmp data
206 : DataVector subcell_data_to_send{sliced_data_in_direction.size() +
207 : rdmp_size};
208 : // Note: Currently we interpolate our solution to our neighbor FD grid
209 : // even when grid points align but are oriented differently. There's a
210 : // possible optimization for the rare (almost never?) edge case where
211 : // two blocks have the same ghost zone coordinates but have different
212 : // orientations (e.g. RotatedBricks). Since this shouldn't ever happen
213 : // outside of tests, we currently don't bother with it. If we wanted to,
214 : // here's the code:
215 : //
216 : // if (not orientation.is_aligned()) {
217 : // std::array<size_t, Dim> slice_extents{};
218 : // for (size_t d = 0; d < Dim; ++d) {
219 : // gsl::at(slice_extents, d) = subcell_mesh.extents(d);
220 : // }
221 : // gsl::at(slice_extents, direction.dimension()) = ghost_zone_size;
222 : // // Need a view so we only get the subcell data and not the rdmp
223 : // // data
224 : // DataVector subcell_data_to_send_view{
225 : // subcell_data_to_send.data(),
226 : // subcell_data_to_send.size() - rdmp_size};
227 : // orient_variables(make_not_null(&subcell_data_to_send_view),
228 : // sliced_data_in_direction, Index<Dim>{slice_extents},
229 : // orientation);
230 : // } else { std::copy(...); }
231 : //
232 : // Copy over data since it's already oriented from interpolation
233 : std::copy(sliced_data_in_direction.begin(),
234 : sliced_data_in_direction.end(), subcell_data_to_send.begin());
235 : // Copy rdmp data to end of subcell_data_to_send
236 : std::copy(
237 : rdmp_tci_data.max_variables_values.cbegin(),
238 : rdmp_tci_data.max_variables_values.cend(),
239 : std::prev(subcell_data_to_send.end(), static_cast<int>(rdmp_size)));
240 : std::copy(rdmp_tci_data.min_variables_values.cbegin(),
241 : rdmp_tci_data.min_variables_values.cend(),
242 : std::prev(subcell_data_to_send.end(),
243 : static_cast<int>(
244 : rdmp_tci_data.min_variables_values.size())));
245 :
246 : std::tuple<Mesh<Dim>, Mesh<Dim - 1>, std::optional<DataVector>,
247 : std::optional<DataVector>, ::TimeStepId, int>
248 : data{subcell_mesh,
249 : dg_mesh.slice_away(direction.dimension()),
250 : std::move(subcell_data_to_send),
251 : std::nullopt,
252 : next_time_step_id,
253 : tci_decision};
254 :
255 : Parallel::receive_data<
256 : evolution::dg::Tags::BoundaryCorrectionAndGhostCellsInbox<Dim>>(
257 : receiver_proxy[neighbor], time_step_id,
258 : std::pair{DirectionalId<Dim>{direction_from_neighbor, element.id()},
259 : std::move(data)});
260 : }
261 : }
262 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
263 : }
264 : };
265 :
266 : /*!
267 : * \brief Receive the subcell data from our neighbor, and accumulate the data
268 : * from the relaxed discrete maximum principle troubled-cell indicator.
269 : *
270 : * Note:
271 : * - Since we only care about the min/max over all neighbors and ourself at the
272 : * past time, we accumulate all data immediately into the `RdmpTciData`.
273 : * - If the neighbor is using DG and therefore sends boundary correction data
274 : * then that is added into the `evolution::dg::Tags::MortarData` tag
275 : * - The next `TimeStepId` is recorded, but we do not yet support local time
276 : * stepping.
277 : * - This action will never care about what variables are sent for
278 : * reconstruction. It is only responsible for receiving the data and storing
279 : * it in the `NeighborData`.
280 : *
281 : * GlobalCache:
282 : * -Uses: nothing
283 : *
284 : * DataBox:
285 : * - Uses:
286 : * - `domain::Tags::Element<Dim>`
287 : * - `Tags::TimeStepId`
288 : * - `domain::Tags::Mesh<Dim>`
289 : * - `subcell::Tags::Mesh<Dim>`
290 : * - `domain::Tags::Element<Dim>`
291 : * - `Tags::Next<Tags::TimeStepId>`
292 : * - `subcell::Tags::ActiveGrid`
293 : * - `System::variables_tag`
294 : * - Adds: nothing
295 : * - Removes: nothing
296 : * - Modifies:
297 : * - `subcell::Tags::GhostDataForReconstruction<Dim>`
298 : * - `subcell::Tags::DataForRdmpTci`
299 : * - `evolution::dg::Tags::MortarData`
300 : * - `evolution::dg::Tags::MortarNextTemporalId`
301 : */
302 : template <size_t Dim>
303 1 : struct ReceiveDataForReconstruction {
304 : template <typename DbTags, typename... InboxTags, typename ArrayIndex,
305 : typename ActionList, typename ParallelComponent,
306 : typename Metavariables>
307 0 : static Parallel::iterable_action_return_t apply(
308 : db::DataBox<DbTags>& box, tuples::TaggedTuple<InboxTags...>& inboxes,
309 : const Parallel::GlobalCache<Metavariables>& /*cache*/,
310 : const ArrayIndex& /*array_index*/, const ActionList /*meta*/,
311 : const ParallelComponent* const /*meta*/) {
312 : const Element<Dim>& element = db::get<::domain::Tags::Element<Dim>>(box);
313 : const auto number_of_expected_messages = element.neighbors().size();
314 : if (UNLIKELY(number_of_expected_messages == 0)) {
315 : // We have no neighbors, so just continue without doing any work
316 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
317 : }
318 :
319 : using ::operator<<;
320 : using Key = DirectionalId<Dim>;
321 : const auto& current_time_step_id = db::get<::Tags::TimeStepId>(box);
322 : std::map<TimeStepId,
323 : DirectionalIdMap<Dim, std::tuple<Mesh<Dim>, Mesh<Dim - 1>,
324 : std::optional<DataVector>,
325 : std::optional<DataVector>,
326 : ::TimeStepId, int>>>& inbox =
327 : tuples::get<evolution::dg::Tags::BoundaryCorrectionAndGhostCellsInbox<
328 : Metavariables::volume_dim>>(inboxes);
329 : const auto& received = inbox.find(current_time_step_id);
330 : // Check we have at least some data from correct time, and then check that
331 : // we have received all data
332 : if (received == inbox.end() or
333 : received->second.size() != number_of_expected_messages) {
334 : return {Parallel::AlgorithmExecution::Retry, std::nullopt};
335 : }
336 :
337 : // Now that we have received all the data, copy it over as needed.
338 : DirectionalIdMap<
339 : Dim, std::tuple<Mesh<Dim>, Mesh<Dim - 1>, std::optional<DataVector>,
340 : std::optional<DataVector>, ::TimeStepId, int>>
341 : received_data = std::move(inbox[current_time_step_id]);
342 : inbox.erase(current_time_step_id);
343 :
344 : const Mesh<Dim>& subcell_mesh = db::get<Tags::Mesh<Dim>>(box);
345 :
346 : db::mutate<Tags::GhostDataForReconstruction<Dim>, Tags::DataForRdmpTci,
347 : evolution::dg::Tags::MortarData<Dim>,
348 : evolution::dg::Tags::MortarNextTemporalId<Dim>,
349 : domain::Tags::NeighborMesh<Dim>,
350 : evolution::dg::subcell::Tags::NeighborTciDecisions<Dim>>(
351 : [¤t_time_step_id, &element,
352 : ghost_zone_size =
353 : db::get<evolution::dg::subcell::Tags::Reconstructor>(box)
354 : .ghost_zone_size(),
355 : &received_data, &subcell_mesh](
356 : const gsl::not_null<DirectionalIdMap<Dim, GhostData>*>
357 : ghost_data_ptr,
358 : const gsl::not_null<RdmpTciData*> rdmp_tci_data_ptr,
359 : const gsl::not_null<std::unordered_map<
360 : Key, evolution::dg::MortarData<Dim>, boost::hash<Key>>*>
361 : mortar_data,
362 : const gsl::not_null<
363 : std::unordered_map<Key, TimeStepId, boost::hash<Key>>*>
364 : mortar_next_time_step_id,
365 : const gsl::not_null<DirectionalIdMap<Dim, Mesh<Dim>>*>
366 : neighbor_mesh,
367 : const auto neighbor_tci_decisions,
368 : const DirectionalIdMap<Dim, std::optional<intrp::Irregular<Dim>>>&
369 : neighbor_dg_to_fd_interpolants) {
370 : // Remove neighbor meshes for neighbors that don't exist anymore
371 : domain::remove_nonexistent_neighbors(neighbor_mesh, element);
372 :
373 : // Get the next time step id, and also the fluxes data if the neighbor
374 : // is doing DG.
375 : for (auto& received_mortar_data : received_data) {
376 : const auto& mortar_id = received_mortar_data.first;
377 : try {
378 : mortar_next_time_step_id->at(mortar_id) =
379 : std::get<4>(received_mortar_data.second);
380 : } catch (std::exception& e) {
381 : ERROR("Failed retrieving the MortarId: ("
382 : << mortar_id.direction << ',' << mortar_id.id
383 : << ") from the mortar_next_time_step_id. Got exception: "
384 : << e.what());
385 : }
386 : if (std::get<3>(received_mortar_data.second).has_value()) {
387 : mortar_data->at(mortar_id).insert_neighbor_mortar_data(
388 : current_time_step_id,
389 : std::get<1>(received_mortar_data.second),
390 : std::move(*std::get<3>(received_mortar_data.second)));
391 : }
392 : // Set new neighbor mesh
393 : neighbor_mesh->insert_or_assign(
394 : mortar_id, std::get<0>(received_mortar_data.second));
395 : }
396 :
397 : ASSERT(ghost_data_ptr->empty(),
398 : "Should have no elements in the neighbor data when "
399 : "receiving neighbor data");
400 : const size_t number_of_rdmp_vars =
401 : rdmp_tci_data_ptr->max_variables_values.size();
402 : ASSERT(rdmp_tci_data_ptr->min_variables_values.size() ==
403 : number_of_rdmp_vars,
404 : "The number of RDMP variables for which we have a maximum "
405 : "and minimum should be the same, but we have "
406 : << number_of_rdmp_vars << " for the max and "
407 : << rdmp_tci_data_ptr->min_variables_values.size()
408 : << " for the min.");
409 :
410 : for (const auto& [direction, neighbors_in_direction] :
411 : element.neighbors()) {
412 : for (const auto& neighbor : neighbors_in_direction) {
413 : DirectionalId<Dim> directional_element_id{direction, neighbor};
414 : ASSERT(ghost_data_ptr->count(directional_element_id) == 0,
415 : "Found neighbor already inserted in direction "
416 : << direction << " with ElementId " << neighbor);
417 : ASSERT(std::get<2>(received_data[directional_element_id])
418 : .has_value(),
419 : "Received subcell data message that does not contain any "
420 : "actual subcell data for reconstruction.");
421 : // Collect the max/min of u(t^n) for the RDMP as we receive data.
422 : // This reduces the memory footprint.
423 :
424 : evolution::dg::subcell::insert_neighbor_rdmp_and_volume_data(
425 : rdmp_tci_data_ptr, ghost_data_ptr,
426 : *std::get<2>(received_data[directional_element_id]),
427 : number_of_rdmp_vars, directional_element_id,
428 : neighbor_mesh->at(directional_element_id), element,
429 : subcell_mesh, ghost_zone_size,
430 : neighbor_dg_to_fd_interpolants);
431 : ASSERT(neighbor_tci_decisions->contains(directional_element_id),
432 : "The NeighorTciDecisions should contain the neighbor ("
433 : << directional_element_id.direction << ", "
434 : << directional_element_id.id << ") but doesn't");
435 : neighbor_tci_decisions->at(directional_element_id) =
436 : std::get<5>(received_data[directional_element_id]);
437 : }
438 : }
439 : },
440 : make_not_null(&box),
441 : db::get<
442 : evolution::dg::subcell::Tags::InterpolatorsFromNeighborDgToFd<Dim>>(
443 : box));
444 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
445 : }
446 : };
447 : } // namespace evolution::dg::subcell::Actions
|
