Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : #pragma once
5 :
6 : #include <cstddef>
7 : #include <optional>
8 : #include <tuple>
9 : #include <type_traits>
10 : #include <utility>
11 :
12 : #include "DataStructures/DataBox/DataBox.hpp"
13 : #include "Domain/Structure/DirectionalId.hpp"
14 : #include "Domain/Structure/Element.hpp"
15 : #include "Domain/Tags.hpp"
16 : #include "Evolution/DgSubcell/ActiveGrid.hpp"
17 : #include "Evolution/DgSubcell/GhostData.hpp"
18 : #include "Evolution/DgSubcell/Mesh.hpp"
19 : #include "Evolution/DgSubcell/Projection.hpp"
20 : #include "Evolution/DgSubcell/RdmpTciData.hpp"
21 : #include "Evolution/DgSubcell/Reconstruction.hpp"
22 : #include "Evolution/DgSubcell/ReconstructionMethod.hpp"
23 : #include "Evolution/DgSubcell/Tags/ActiveGrid.hpp"
24 : #include "Evolution/DgSubcell/Tags/CellCenteredFlux.hpp"
25 : #include "Evolution/DgSubcell/Tags/Coordinates.hpp"
26 : #include "Evolution/DgSubcell/Tags/DataForRdmpTci.hpp"
27 : #include "Evolution/DgSubcell/Tags/DidRollback.hpp"
28 : #include "Evolution/DgSubcell/Tags/GhostDataForReconstruction.hpp"
29 : #include "Evolution/DgSubcell/Tags/InitialTciData.hpp"
30 : #include "Evolution/DgSubcell/Tags/Interpolators.hpp"
31 : #include "Evolution/DgSubcell/Tags/Jacobians.hpp"
32 : #include "Evolution/DgSubcell/Tags/Mesh.hpp"
33 : #include "Evolution/DgSubcell/Tags/MeshForGhostData.hpp"
34 : #include "Evolution/DgSubcell/Tags/ReconstructionOrder.hpp"
35 : #include "Evolution/DgSubcell/Tags/StepsSinceTciCall.hpp"
36 : #include "Evolution/DgSubcell/Tags/SubcellOptions.hpp"
37 : #include "Evolution/DgSubcell/Tags/TciCallsSinceRollback.hpp"
38 : #include "Evolution/DgSubcell/Tags/TciGridHistory.hpp"
39 : #include "Evolution/DgSubcell/Tags/TciStatus.hpp"
40 : #include "Evolution/Initialization/SetVariables.hpp"
41 : #include "NumericalAlgorithms/Interpolation/IrregularInterpolant.hpp"
42 : #include "NumericalAlgorithms/Spectral/Mesh.hpp"
43 : #include "Parallel/AlgorithmExecution.hpp"
44 : #include "Parallel/GlobalCache.hpp"
45 : #include "Utilities/CallWithDynamicType.hpp"
46 : #include "Utilities/ContainerHelpers.hpp"
47 : #include "Utilities/ErrorHandling/Error.hpp"
48 : #include "Utilities/TMPL.hpp"
49 : #include "Utilities/TaggedTuple.hpp"
50 :
51 : /// \cond
52 : namespace Tags {
53 : template <typename Tag>
54 : struct HistoryEvolvedVariables;
55 : } // namespace Tags
56 : /// \endcond
57 :
58 : namespace evolution::dg::subcell::Actions {
59 : /*!
60 : * \brief Initialize the subcell grid, including the size of the evolved
61 : * `Variables` and, if present, primitive `Variables`.
62 : *
63 : * By default sets the element to `subcell::ActiveGrid::Subcell` unless it
64 : * is not allowed to use subcell either because it is at an external boundary
65 : * or because it or one of its neighbors has been marked as DG-only.
66 : *
67 : * GlobalCache:
68 : * - Uses:
69 : * - `subcell::Tags::SubcellOptions`
70 : *
71 : * DataBox:
72 : * - Uses:
73 : * - `domain::Tags::Mesh<Dim>`
74 : * - `domain::Tags::Element<Dim>`
75 : * - `System::variables_tag`
76 : * - Adds:
77 : * - `subcell::Tags::Mesh<Dim>`
78 : * - `subcell::Tags::MeshForGhostData<Dim>`
79 : * - `subcell::Tags::ActiveGrid`
80 : * - `subcell::Tags::DidRollback`
81 : * - `subcell::Tags::TciGridHistory`
82 : * - `subcell::Tags::TciCallsSinceRollback`
83 : * - `subcell::Tags::GhostDataForReconstruction<Dim>`
84 : * - `subcell::Tags::TciDecision`
85 : * - `subcell::Tags::DataForRdmpTci`
86 : * - `subcell::fd::Tags::InverseJacobianLogicalToGrid<Dim>`
87 : * - `subcell::fd::Tags::DetInverseJacobianLogicalToGrid`
88 : * - `subcell::Tags::LogicalCoordinates<Dim>`
89 : * - `subcell::Tags::ReconstructionOrder<Dim>` (set as `std::nullopt`)
90 : * - `subcell::Tags::Coordinates<Dim, Frame::Grid>` (as compute tag)
91 : * - `subcell::Tags::Coordinates<Dim, Frame::Inertial>` (as compute tag)
92 : * - Removes: nothing
93 : * - Modifies:
94 : * - `System::variables_tag` and `System::primitive_variables_tag` if the cell
95 : * is troubled
96 : * - `Tags::dt<System::variables_tag>` if the cell is troubled
97 : */
98 : template <size_t Dim, typename System, bool UseNumericInitialData>
99 1 : struct SetSubcellGrid {
100 0 : using const_global_cache_tags = tmpl::list<Tags::SubcellOptions<Dim>>;
101 :
102 0 : using simple_tags = tmpl::list<
103 : Tags::ActiveGrid, Tags::DidRollback, Tags::TciGridHistory,
104 : Tags::TciCallsSinceRollback, Tags::StepsSinceTciCall,
105 : evolution::dg::subcell::Tags::MeshForGhostData<Dim>,
106 : Tags::GhostDataForReconstruction<Dim>, Tags::TciDecision,
107 : Tags::NeighborTciDecisions<Dim>, Tags::DataForRdmpTci,
108 : subcell::Tags::CellCenteredFlux<typename System::flux_variables, Dim>,
109 : subcell::Tags::ReconstructionOrder<Dim>,
110 : evolution::dg::subcell::Tags::InterpolatorsFromFdToNeighborFd<Dim>,
111 : evolution::dg::subcell::Tags::InterpolatorsFromDgToNeighborFd<Dim>,
112 : evolution::dg::subcell::Tags::InterpolatorsFromNeighborDgToFd<Dim>,
113 : typename System::variables_tag,
114 : evolution::dg::subcell::Tags::ExtensionDirections<Dim>>;
115 0 : using compute_tags =
116 : tmpl::list<Tags::MeshCompute<Dim>, Tags::LogicalCoordinatesCompute<Dim>,
117 : ::domain::Tags::MappedCoordinates<
118 : ::domain::Tags::ElementMap<Dim, Frame::Grid>,
119 : subcell::Tags::Coordinates<Dim, Frame::ElementLogical>,
120 : subcell::Tags::Coordinates>,
121 : Tags::InertialCoordinatesCompute<
122 : ::domain::CoordinateMaps::Tags::CoordinateMap<
123 : Dim, Frame::Grid, Frame::Inertial>>,
124 : fd::Tags::InverseJacobianLogicalToGridCompute<
125 : ::domain::Tags::ElementMap<Dim, Frame::Grid>, Dim>,
126 : fd::Tags::DetInverseJacobianLogicalToGridCompute<Dim>,
127 : fd::Tags::InverseJacobianLogicalToInertialCompute<
128 : ::domain::CoordinateMaps::Tags::CoordinateMap<
129 : Dim, Frame::Grid, Frame::Inertial>,
130 : Dim>,
131 : fd::Tags::DetInverseJacobianLogicalToInertialCompute<
132 : ::domain::CoordinateMaps::Tags::CoordinateMap<
133 : Dim, Frame::Grid, Frame::Inertial>,
134 : Dim>>;
135 :
136 : template <typename DbTagsList, typename... InboxTags, typename ArrayIndex,
137 : typename ActionList, typename ParallelComponent,
138 : typename Metavariables>
139 0 : static Parallel::iterable_action_return_t apply(
140 : db::DataBox<DbTagsList>& box,
141 : [[maybe_unused]] const tuples::TaggedTuple<InboxTags...>& inboxes,
142 : [[maybe_unused]] const Parallel::GlobalCache<Metavariables>& cache,
143 : [[maybe_unused]] const ArrayIndex& array_index, ActionList /*meta*/,
144 : const ParallelComponent* const /*meta*/) {
145 : const SubcellOptions& subcell_options =
146 : db::get<Tags::SubcellOptions<Dim>>(box);
147 : const Mesh<Dim>& dg_mesh = db::get<::domain::Tags::Mesh<Dim>>(box);
148 : const Mesh<Dim>& subcell_mesh = db::get<subcell::Tags::Mesh<Dim>>(box);
149 : const Element<Dim>& element = db::get<::domain::Tags::Element<Dim>>(box);
150 :
151 : for (size_t d = 0; d < Dim; ++d) {
152 : if (subcell_options.persson_num_highest_modes() >= dg_mesh.extents(d)) {
153 : ERROR("Number of the highest modes to be monitored by the Persson TCI ("
154 : << subcell_options.persson_num_highest_modes()
155 : << ") must be smaller than the extent of the DG mesh ("
156 : << dg_mesh.extents(d) << ").");
157 : }
158 : }
159 :
160 : // Loop over block neighbors and if neighbor id is inside of
161 : // subcell_options.only_dg_block_ids(), then bordering DG-only block
162 : const bool bordering_dg_block = alg::any_of(
163 : element.neighbors(),
164 : [&subcell_options](const auto& direction_and_neighbor) {
165 : const size_t first_block_id =
166 : direction_and_neighbor.second.ids().begin()->block_id();
167 : return std::binary_search(subcell_options.only_dg_block_ids().begin(),
168 : subcell_options.only_dg_block_ids().end(),
169 : first_block_id);
170 : });
171 :
172 : const bool subcell_allowed_in_element =
173 : not std::binary_search(subcell_options.only_dg_block_ids().begin(),
174 : subcell_options.only_dg_block_ids().end(),
175 : element.id().block_id()) and
176 : not bordering_dg_block;
177 : const bool cell_is_not_on_external_boundary =
178 : db::get<::domain::Tags::Element<Dim>>(box)
179 : .external_boundaries()
180 : .empty();
181 :
182 : constexpr bool subcell_enabled_at_external_boundary =
183 : Metavariables::SubcellOptions::subcell_enabled_at_external_boundary;
184 :
185 : db::mutate<Tags::NeighborTciDecisions<Dim>>(
186 : [&element](const auto neighbor_decisions_ptr) {
187 : neighbor_decisions_ptr->clear();
188 : for (const auto& [direction, neighbors_in_direction] :
189 : element.neighbors()) {
190 : for (const auto& neighbor : neighbors_in_direction.ids()) {
191 : neighbor_decisions_ptr->insert(
192 : std::pair{DirectionalId<Dim>{direction, neighbor}, 0});
193 : }
194 : }
195 : },
196 : make_not_null(&box));
197 :
198 : db::mutate_apply<
199 : tmpl::list<Tags::ActiveGrid, Tags::DidRollback,
200 : typename System::variables_tag, subcell::Tags::TciDecision,
201 : subcell::Tags::TciCallsSinceRollback,
202 : subcell::Tags::StepsSinceTciCall>,
203 : tmpl::list<>>(
204 : [&cell_is_not_on_external_boundary, &dg_mesh,
205 : subcell_allowed_in_element, &subcell_mesh](
206 : const gsl::not_null<ActiveGrid*> active_grid_ptr,
207 : const gsl::not_null<bool*> did_rollback_ptr,
208 : const auto active_vars_ptr,
209 : const gsl::not_null<int*> tci_decision_ptr,
210 : const gsl::not_null<size_t*> tci_calls_since_rollback_ptr,
211 : const gsl::not_null<size_t*> steps_since_tci_call_ptr) {
212 : // We don't consider setting the initial grid to subcell as rolling
213 : // back. Since no time step is undone, we just continue on the
214 : // subcells as a normal solve.
215 : *did_rollback_ptr = false;
216 :
217 : if ((cell_is_not_on_external_boundary or
218 : subcell_enabled_at_external_boundary) and
219 : subcell_allowed_in_element) {
220 : *active_grid_ptr = ActiveGrid::Subcell;
221 : active_vars_ptr->initialize(subcell_mesh.number_of_grid_points());
222 : } else {
223 : *active_grid_ptr = ActiveGrid::Dg;
224 : active_vars_ptr->initialize(dg_mesh.number_of_grid_points());
225 : }
226 :
227 : *tci_decision_ptr = 0;
228 : *tci_calls_since_rollback_ptr = 0;
229 : *steps_since_tci_call_ptr = 0;
230 : },
231 : make_not_null(&box));
232 : if constexpr (System::has_primitive_and_conservative_vars) {
233 : db::mutate<typename System::primitive_variables_tag>(
234 : [&dg_mesh, &subcell_mesh](const auto prim_vars_ptr,
235 : const auto active_grid) {
236 : if (active_grid == ActiveGrid::Dg) {
237 : prim_vars_ptr->initialize(dg_mesh.number_of_grid_points());
238 : } else {
239 : prim_vars_ptr->initialize(subcell_mesh.number_of_grid_points());
240 : }
241 : },
242 : make_not_null(&box), db::get<Tags::ActiveGrid>(box));
243 : }
244 : if constexpr (not UseNumericInitialData) {
245 : if (db::get<Tags::ActiveGrid>(box) ==
246 : evolution::dg::subcell::ActiveGrid::Dg) {
247 : evolution::Initialization::Actions::SetVariables<
248 : ::domain::Tags::Coordinates<Dim, Frame::ElementLogical>>::
249 : apply(box, inboxes, cache, array_index, ActionList{},
250 : std::add_pointer_t<ParallelComponent>{nullptr});
251 : } else {
252 : evolution::Initialization::Actions::
253 : SetVariables<Tags::Coordinates<Dim, Frame::ElementLogical>>::apply(
254 : box, inboxes, cache, array_index, ActionList{},
255 : std::add_pointer_t<ParallelComponent>{nullptr});
256 : }
257 : }
258 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
259 : }
260 : };
261 :
262 : /*!
263 : * \brief Sets the RDMP data from the initial data and sends it to neighboring
264 : * elements.
265 : *
266 : * GlobalCache:
267 : * - Uses:
268 : * - `ParallelComponent` proxy
269 : *
270 : * DataBox:
271 : * - Uses:
272 : * - `domain::Tags::Element<Dim>`
273 : * - `subcell::Tags::DataForRdmpTci`
274 : * - `subcell::Tags::InitialTciData`
275 : * - whatever `SetInitialRdmpData` uses
276 : * - Adds: nothing
277 : * - Removes: nothing
278 : * - Modifies:
279 : * - whatever `SetInitialRdmpData` mutates
280 : */
281 : template <size_t Dim, typename SetInitialRdmpData>
282 1 : struct SetAndCommunicateInitialRdmpData {
283 0 : using inbox_tags =
284 : tmpl::list<evolution::dg::subcell::Tags::InitialTciData<Dim>>;
285 :
286 : template <typename DbTagsList, typename... InboxTags, typename ArrayIndex,
287 : typename ActionList, typename ParallelComponent,
288 : typename Metavariables>
289 0 : static Parallel::iterable_action_return_t apply(
290 : db::DataBox<DbTagsList>& box,
291 : const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
292 : Parallel::GlobalCache<Metavariables>& cache,
293 : const ArrayIndex& /*array_index*/, ActionList /*meta*/,
294 : const ParallelComponent* const /*meta*/) {
295 : // Get the RDMP data on this element and then initialize it.
296 : db::mutate_apply<SetInitialRdmpData>(make_not_null(&box));
297 :
298 : // Send RDMP data to neighbors
299 : const auto& element = db::get<domain::Tags::Element<Dim>>(box);
300 : const auto& rdmp_data =
301 : db::get<evolution::dg::subcell::Tags::DataForRdmpTci>(box);
302 : auto& receiver_proxy =
303 : Parallel::get_parallel_component<ParallelComponent>(cache);
304 : for (const auto& [direction, neighbors] : element.neighbors()) {
305 : for (const auto& neighbor : neighbors) {
306 : const auto& orientation = neighbors.orientation(neighbor);
307 : const auto direction_from_neighbor = orientation(direction.opposite());
308 : evolution::dg::subcell::InitialTciData data{{}, rdmp_data};
309 : // We use temporal ID 0 for sending RDMP data
310 : const int temporal_id = 0;
311 : Parallel::receive_data<
312 : evolution::dg::subcell::Tags::InitialTciData<Dim>>(
313 : receiver_proxy[neighbor], temporal_id,
314 : std::make_pair(
315 : DirectionalId<Dim>{direction_from_neighbor, element.id()},
316 : std::move(data)));
317 : }
318 : }
319 :
320 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
321 : }
322 : };
323 :
324 : /*!
325 : * \brief Apply the TCI on the FD grid to the initial data and send the TCI
326 : * decision to neighboring elements.
327 : *
328 : * GlobalCache:
329 : * - Uses:
330 : * - `ParallelComponent` proxy
331 : *
332 : * DataBox:
333 : * - Uses:
334 : * - `domain::Tags::Element<Dim>`
335 : * - `subcell::Tags::DataForRdmpTci`
336 : * - `subcell::Tags::InitialTciData`
337 : * - `subcell::Tags::SubcellOptions`
338 : * - `subcell::Tags::ActiveGrid`
339 : * - whatever `TciOnFdGridMutator` uses
340 : * - Adds: nothing
341 : * - Removes: nothing
342 : * - Modifies:
343 : * - `subcell::Tags::DataForRdmpTci`
344 : * - `subcell::Tags::TciDecision`
345 : */
346 : template <size_t Dim, typename System, typename TciOnFdGridMutator>
347 1 : struct ComputeAndSendTciOnInitialGrid {
348 0 : using inbox_tags =
349 : tmpl::list<evolution::dg::subcell::Tags::InitialTciData<Dim>>;
350 :
351 : template <typename DbTagsList, typename... InboxTags, typename ArrayIndex,
352 : typename ActionList, typename ParallelComponent,
353 : typename Metavariables>
354 0 : static Parallel::iterable_action_return_t apply(
355 : db::DataBox<DbTagsList>& box, tuples::TaggedTuple<InboxTags...>& inboxes,
356 : Parallel::GlobalCache<Metavariables>& cache,
357 : const ArrayIndex& /*array_index*/, ActionList /*meta*/,
358 : const ParallelComponent* const /*meta*/) {
359 : const Element<Dim>& element = db::get<::domain::Tags::Element<Dim>>(box);
360 :
361 : // Check if we have received all RDMP data.
362 : if (LIKELY(element.number_of_neighbors() != 0)) {
363 : auto& inbox =
364 : tuples::get<evolution::dg::subcell::Tags::InitialTciData<Dim>>(
365 : inboxes);
366 : const auto& received = inbox.find(0);
367 : if (received == inbox.end() or
368 : received->second.size() != element.number_of_neighbors()) {
369 : return {Parallel::AlgorithmExecution::Retry, std::nullopt};
370 : }
371 :
372 : db::mutate<evolution::dg::subcell::Tags::DataForRdmpTci>(
373 : [&element, &received](const auto rdmp_tci_data_ptr) {
374 : (void)element;
375 : const size_t number_of_rdmp_vars =
376 : rdmp_tci_data_ptr->max_variables_values.size();
377 : ASSERT(rdmp_tci_data_ptr->max_variables_values.size() ==
378 : number_of_rdmp_vars,
379 : "The number of local max vars is "
380 : << number_of_rdmp_vars
381 : << " while the number of local min vars is "
382 : << rdmp_tci_data_ptr->max_variables_values.size()
383 : << " the local element ID is " << element.id());
384 : for (const auto& [direction_and_neighbor_element_id,
385 : neighbor_initial_tci_data] : received->second) {
386 : ASSERT(neighbor_initial_tci_data.initial_rdmp_data.has_value(),
387 : "Neighbor in direction "
388 : << direction_and_neighbor_element_id.direction()
389 : << " with element ID "
390 : << direction_and_neighbor_element_id.id() << " of "
391 : << element.id()
392 : << " didn't send initial TCI data correctly");
393 : ASSERT(
394 : neighbor_initial_tci_data.initial_rdmp_data.value()
395 : .max_variables_values.size() == number_of_rdmp_vars,
396 : "The number of local RDMP vars is "
397 : << number_of_rdmp_vars
398 : << " while the number of remote max vars is "
399 : << neighbor_initial_tci_data.initial_rdmp_data.value()
400 : .max_variables_values.size()
401 : << " the local element ID is " << element.id()
402 : << " and the remote id is "
403 : << direction_and_neighbor_element_id.id());
404 : ASSERT(
405 : neighbor_initial_tci_data.initial_rdmp_data.value()
406 : .min_variables_values.size() == number_of_rdmp_vars,
407 : "The number of local RDMP vars is "
408 : << number_of_rdmp_vars
409 : << " while the number of remote min vars is "
410 : << neighbor_initial_tci_data.initial_rdmp_data.value()
411 : .min_variables_values.size()
412 : << " the local element ID is " << element.id()
413 : << " and the remote id is "
414 : << direction_and_neighbor_element_id.id());
415 : for (size_t var_index = 0; var_index < number_of_rdmp_vars;
416 : ++var_index) {
417 : rdmp_tci_data_ptr->max_variables_values[var_index] =
418 : std::max(rdmp_tci_data_ptr->max_variables_values[var_index],
419 : neighbor_initial_tci_data.initial_rdmp_data.value()
420 : .max_variables_values[var_index]);
421 : rdmp_tci_data_ptr->min_variables_values[var_index] =
422 : std::min(rdmp_tci_data_ptr->min_variables_values[var_index],
423 : neighbor_initial_tci_data.initial_rdmp_data.value()
424 : .min_variables_values[var_index]);
425 : }
426 : }
427 : },
428 : make_not_null(&box));
429 : inbox.erase(received);
430 : }
431 :
432 : const auto send_tci_decision = [&cache, &element](const int tci_decision) {
433 : if (UNLIKELY(element.number_of_neighbors() == 0)) {
434 : return;
435 : }
436 : auto& receiver_proxy =
437 : Parallel::get_parallel_component<ParallelComponent>(cache);
438 : for (const auto& [direction, neighbors] : element.neighbors()) {
439 : for (const auto& neighbor : neighbors) {
440 : const auto& orientation = neighbors.orientation(neighbor);
441 : const auto direction_from_neighbor =
442 : orientation(direction.opposite());
443 : evolution::dg::subcell::InitialTciData data{tci_decision, {}};
444 : // We use temporal ID 1 for ending the TCI decision.
445 : const int temporal_id = 1;
446 : Parallel::receive_data<
447 : evolution::dg::subcell::Tags::InitialTciData<Dim>>(
448 : receiver_proxy[neighbor], temporal_id,
449 : std::make_pair(
450 : DirectionalId<Dim>{direction_from_neighbor, element.id()},
451 : std::move(data)));
452 : }
453 : }
454 : };
455 :
456 : const SubcellOptions& subcell_options =
457 : db::get<Tags::SubcellOptions<Dim>>(box);
458 :
459 : if (subcell_options.always_use_subcells() or
460 : get<Tags::ActiveGrid>(box) == ActiveGrid::Dg) {
461 : db::mutate<Tags::TciDecision>(
462 : [](const gsl::not_null<int*> tci_decision_ptr) {
463 : *tci_decision_ptr = 0;
464 : },
465 : make_not_null(&box));
466 : send_tci_decision(0);
467 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
468 : }
469 :
470 : // Now run the TCI to see if we could switch back to DG.
471 : const std::tuple<int, evolution::dg::subcell::RdmpTciData> tci_result =
472 : db::mutate_apply<TciOnFdGridMutator>(
473 : make_not_null(&box), subcell_options.persson_exponent() + 1.0,
474 : false);
475 :
476 : db::mutate<Tags::TciDecision>(
477 : [&tci_result](const gsl::not_null<int*> tci_decision_ptr) {
478 : *tci_decision_ptr = std::get<0>(tci_result);
479 : },
480 : make_not_null(&box));
481 : send_tci_decision(std::get<0>(tci_result));
482 :
483 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
484 : }
485 : };
486 :
487 : /*!
488 : * \brief Using the local and neighboring TCI decisions, switches the element to
489 : * DG if the DG solution was determined to be admissible.
490 : *
491 : * GlobalCache:
492 : * - Uses:
493 : * - `ParallelComponent` proxy
494 : *
495 : * DataBox:
496 : * - Uses:
497 : * - `domain::Tags::Element<Dim>`
498 : * - `subcell::Tags::DataForRdmpTci`
499 : * - `subcell::Tags::InitialTciData`
500 : * - `subcell::Tags::SubcellOptions`
501 : * - `subcell::Tags::ActiveGrid`
502 : * - whatever `TciOnFdGridMutator` uses
503 : * - Adds: nothing
504 : * - Removes: nothing
505 : * - Modifies:
506 : * - `subcell::Tags::NeighborTciDecisions`
507 : * - `System::variables_tag`
508 : * - `Tags::HistoryEvolvedVariables<System::variables_tag>`
509 : * - `subcell::Tags::GhostDataForReconstruction`
510 : * - `subcell::Tags::TciGridHistory`
511 : * - `subcell::Tags::CellCenteredFlux`
512 : */
513 : template <size_t Dim, typename System>
514 1 : struct SetInitialGridFromTciData {
515 : template <typename DbTagsList, typename... InboxTags, typename ArrayIndex,
516 : typename ActionList, typename ParallelComponent,
517 : typename Metavariables>
518 0 : static Parallel::iterable_action_return_t apply(
519 : db::DataBox<DbTagsList>& box, tuples::TaggedTuple<InboxTags...>& inboxes,
520 : const Parallel::GlobalCache<Metavariables>& /*cache*/,
521 : const ArrayIndex& /*array_index*/, ActionList /*meta*/,
522 : const ParallelComponent* const /*meta*/) {
523 : const Element<Dim>& element = db::get<::domain::Tags::Element<Dim>>(box);
524 : if (LIKELY(element.number_of_neighbors() != 0)) {
525 : auto& inbox =
526 : tuples::get<evolution::dg::subcell::Tags::InitialTciData<Dim>>(
527 : inboxes);
528 : const auto& received = inbox.find(1);
529 : // Check if we have received all TCI decisions.
530 : if (received == inbox.end() or
531 : received->second.size() != element.number_of_neighbors()) {
532 : return {Parallel::AlgorithmExecution::Retry, std::nullopt};
533 : }
534 :
535 : db::mutate<evolution::dg::subcell::Tags::NeighborTciDecisions<Dim>>(
536 : [&element, &received](const auto neighbor_tci_decisions_ptr) {
537 : for (const auto& [directional_element_id,
538 : neighbor_initial_tci_data] : received->second) {
539 : (void)element;
540 : ASSERT(neighbor_initial_tci_data.tci_status.has_value(),
541 : "Neighbor in direction "
542 : << directional_element_id.direction()
543 : << " with element ID " << directional_element_id.id()
544 : << " of " << element.id()
545 : << " didn't send initial TCI decision correctly");
546 : neighbor_tci_decisions_ptr->at(directional_element_id) =
547 : neighbor_initial_tci_data.tci_status.value();
548 : }
549 : },
550 : make_not_null(&box));
551 : inbox.erase(received);
552 : }
553 :
554 : if (get<Tags::ActiveGrid>(box) == ActiveGrid::Dg) {
555 : // In this case we are allowed to only do DG in this element. No need to
556 : // even do any checks.
557 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
558 : }
559 :
560 : const SubcellOptions& subcell_options =
561 : db::get<Tags::SubcellOptions<Dim>>(box);
562 :
563 : bool cell_is_troubled =
564 : subcell_options.always_use_subcells() or
565 : (subcell_options.use_halo() and [&box]() -> bool {
566 : for (const auto& [_, neighbor_decision] :
567 : db::get<evolution::dg::subcell::Tags::NeighborTciDecisions<Dim>>(
568 : box)) {
569 : if (neighbor_decision != 0) {
570 : return true;
571 : }
572 : }
573 : return false;
574 : }()) or
575 : (db::get<Tags::TciDecision>(box) != 0);
576 :
577 : if (not cell_is_troubled) {
578 : using variables_tag = typename System::variables_tag;
579 : using flux_variables = typename System::flux_variables;
580 :
581 : const Mesh<Dim>& dg_mesh = db::get<::domain::Tags::Mesh<Dim>>(box);
582 : const Mesh<Dim>& subcell_mesh = db::get<subcell::Tags::Mesh<Dim>>(box);
583 : db::mutate<
584 : variables_tag, ::Tags::HistoryEvolvedVariables<variables_tag>,
585 : Tags::ActiveGrid, subcell::Tags::GhostDataForReconstruction<Dim>,
586 : evolution::dg::subcell::Tags::TciGridHistory,
587 : evolution::dg::subcell::Tags::CellCenteredFlux<flux_variables, Dim>>(
588 : [&dg_mesh, &subcell_mesh, &subcell_options](
589 : const auto active_vars_ptr, const auto active_history_ptr,
590 : const gsl::not_null<ActiveGrid*> active_grid_ptr,
591 : const auto subcell_ghost_data_ptr,
592 : const gsl::not_null<
593 : std::deque<evolution::dg::subcell::ActiveGrid>*>
594 : tci_grid_history_ptr,
595 : const auto subcell_cell_centered_fluxes) {
596 : // Note: strictly speaking, to be conservative this should
597 : // reconstruct uJ instead of u.
598 : *active_vars_ptr = fd::reconstruct(
599 : *active_vars_ptr, dg_mesh, subcell_mesh.extents(),
600 : subcell_options.reconstruction_method());
601 :
602 : // Reconstruct the DG solution for each time in the time stepper
603 : // history
604 : active_history_ptr->map_entries(
605 : [&dg_mesh, &subcell_mesh, &subcell_options](const auto entry) {
606 : *entry =
607 : fd::reconstruct(*entry, dg_mesh, subcell_mesh.extents(),
608 : subcell_options.reconstruction_method());
609 : });
610 : *active_grid_ptr = ActiveGrid::Dg;
611 :
612 : // Clear the neighbor data needed for subcell reconstruction since
613 : // we have now completed the time step.
614 : subcell_ghost_data_ptr->clear();
615 :
616 : // Clear the TCI grid history since we don't need to use it when on
617 : // the DG grid.
618 : tci_grid_history_ptr->clear();
619 :
620 : // Clear the allocation for the cell-centered fluxes.
621 : *subcell_cell_centered_fluxes = std::nullopt;
622 : },
623 : make_not_null(&box));
624 : }
625 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
626 : }
627 : };
628 : } // namespace evolution::dg::subcell::Actions
|
