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