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