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