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 : const bool bordering_dg_block =
173 : evolution::dg::borders_only_dg_block(element, only_dg_block_ids);
174 :
175 : // Non-hypercube topologies (e.g. spherical shells) can never use subcell
176 : // and are automatically treated as DG-only.
177 : const bool subcell_allowed_in_element =
178 : fd::dg_mesh_supports_subcell(dg_mesh) and
179 : not alg::found(only_dg_block_ids, element.id().block_id()) and
180 : not bordering_dg_block;
181 : const bool cell_is_not_on_external_boundary =
182 : db::get<::domain::Tags::Element<Dim>>(box)
183 : .external_boundaries()
184 : .empty();
185 :
186 : constexpr bool subcell_enabled_at_external_boundary =
187 : Metavariables::SubcellOptions::subcell_enabled_at_external_boundary;
188 :
189 : // Non-hypercube elements (e.g. spherical shells) cannot use subcell and may
190 : // have > 24 neighbors, which would overflow the fixed-size
191 : // DirectionalIdMap, so their map is left empty. Hypercube elements that
192 : // are forced to DG (e.g. because they border a DG-only block) still have
193 : // a bounded neighbor count and continue to track TCI decisions normally.
194 : //
195 : // Within a hypercube element, also skip MultipleNonconforming directions.
196 : // These arise on the "one" side of a many-to-one interface where
197 : // bordering a non-hypercube element
198 : if (fd::dg_mesh_supports_subcell(dg_mesh)) {
199 : db::mutate<Tags::NeighborTciDecisions<Dim>>(
200 : [&element](const auto neighbor_decisions_ptr) {
201 : neighbor_decisions_ptr->clear();
202 : for (const auto& [direction, neighbors_in_direction] :
203 : element.neighbors()) {
204 : if (element.face_types().at(direction) ==
205 : ::domain::FaceType::MultipleNonconforming) {
206 : continue;
207 : }
208 : for (const auto& neighbor : neighbors_in_direction.ids()) {
209 : neighbor_decisions_ptr->insert(
210 : std::pair{DirectionalId<Dim>{direction, neighbor}, 0});
211 : }
212 : }
213 : },
214 : make_not_null(&box));
215 : }
216 :
217 : db::mutate_apply<
218 : tmpl::list<Tags::ActiveGrid, Tags::DidRollback,
219 : typename System::variables_tag, subcell::Tags::TciDecision,
220 : subcell::Tags::TciCallsSinceRollback,
221 : subcell::Tags::StepsSinceTciCall>,
222 : tmpl::list<>>(
223 : [&cell_is_not_on_external_boundary, &dg_mesh,
224 : subcell_allowed_in_element, &subcell_mesh](
225 : const gsl::not_null<ActiveGrid*> active_grid_ptr,
226 : const gsl::not_null<bool*> did_rollback_ptr,
227 : const auto active_vars_ptr,
228 : const gsl::not_null<int*> tci_decision_ptr,
229 : const gsl::not_null<size_t*> tci_calls_since_rollback_ptr,
230 : const gsl::not_null<size_t*> steps_since_tci_call_ptr) {
231 : // We don't consider setting the initial grid to subcell as rolling
232 : // back. Since no time step is undone, we just continue on the
233 : // subcells as a normal solve.
234 : *did_rollback_ptr = false;
235 :
236 : if ((cell_is_not_on_external_boundary or
237 : subcell_enabled_at_external_boundary) and
238 : subcell_allowed_in_element) {
239 : *active_grid_ptr = ActiveGrid::Subcell;
240 : active_vars_ptr->initialize(subcell_mesh.number_of_grid_points());
241 : } else {
242 : *active_grid_ptr = ActiveGrid::Dg;
243 : active_vars_ptr->initialize(dg_mesh.number_of_grid_points());
244 : }
245 :
246 : *tci_decision_ptr = 0;
247 : *tci_calls_since_rollback_ptr = 0;
248 : *steps_since_tci_call_ptr = 0;
249 : },
250 : make_not_null(&box));
251 : if constexpr (System::has_primitive_and_conservative_vars) {
252 : db::mutate<typename System::primitive_variables_tag>(
253 : [&dg_mesh, &subcell_mesh](const auto prim_vars_ptr,
254 : const auto active_grid) {
255 : if (active_grid == ActiveGrid::Dg) {
256 : prim_vars_ptr->initialize(dg_mesh.number_of_grid_points());
257 : } else {
258 : prim_vars_ptr->initialize(subcell_mesh.number_of_grid_points());
259 : }
260 : },
261 : make_not_null(&box), db::get<Tags::ActiveGrid>(box));
262 : }
263 : if constexpr (not UseNumericInitialData) {
264 : if (db::get<Tags::ActiveGrid>(box) ==
265 : evolution::dg::subcell::ActiveGrid::Dg) {
266 : evolution::Initialization::Actions::SetVariables<
267 : ::domain::Tags::Coordinates<Dim, Frame::ElementLogical>>::
268 : apply(box, inboxes, cache, array_index, ActionList{},
269 : std::add_pointer_t<ParallelComponent>{nullptr});
270 : } else {
271 : evolution::Initialization::Actions::
272 : SetVariables<Tags::Coordinates<Dim, Frame::ElementLogical>>::apply(
273 : box, inboxes, cache, array_index, ActionList{},
274 : std::add_pointer_t<ParallelComponent>{nullptr});
275 : }
276 : }
277 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
278 : }
279 : };
280 :
281 : /*!
282 : * \brief Sets the RDMP data from the initial data and sends it to neighboring
283 : * elements.
284 : *
285 : * GlobalCache:
286 : * - Uses:
287 : * - `ParallelComponent` proxy
288 : *
289 : * DataBox:
290 : * - Uses:
291 : * - `domain::Tags::Element<Dim>`
292 : * - `subcell::Tags::DataForRdmpTci`
293 : * - `subcell::Tags::InitialTciData`
294 : * - whatever `SetInitialRdmpData` uses
295 : * - Adds: nothing
296 : * - Removes: nothing
297 : * - Modifies:
298 : * - whatever `SetInitialRdmpData` mutates
299 : */
300 : template <size_t Dim, typename SetInitialRdmpData>
301 1 : struct SetAndCommunicateInitialRdmpData {
302 0 : using inbox_tags =
303 : tmpl::list<evolution::dg::subcell::Tags::InitialTciData<Dim>>;
304 :
305 : template <typename DbTagsList, typename... InboxTags, typename ArrayIndex,
306 : typename ActionList, typename ParallelComponent,
307 : typename Metavariables>
308 0 : static Parallel::iterable_action_return_t apply(
309 : db::DataBox<DbTagsList>& box,
310 : const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
311 : Parallel::GlobalCache<Metavariables>& cache,
312 : const ArrayIndex& /*array_index*/, ActionList /*meta*/,
313 : const ParallelComponent* const /*meta*/) {
314 : // Get the RDMP data on this element and then initialize it.
315 : db::mutate_apply<SetInitialRdmpData>(make_not_null(&box));
316 :
317 : // Send RDMP data to neighbors
318 : const auto& element = db::get<domain::Tags::Element<Dim>>(box);
319 : const auto& rdmp_data =
320 : db::get<evolution::dg::subcell::Tags::DataForRdmpTci>(box);
321 : auto& receiver_proxy =
322 : Parallel::get_parallel_component<ParallelComponent>(cache);
323 : for (const auto& [direction, neighbors] : element.neighbors()) {
324 : for (const auto& neighbor : neighbors) {
325 : const auto& orientation = neighbors.orientation(neighbor);
326 : const auto direction_from_neighbor = orientation(direction.opposite());
327 : evolution::dg::subcell::InitialTciData data{{}, rdmp_data};
328 : // We use temporal ID 0 for sending RDMP data
329 : const int temporal_id = 0;
330 : Parallel::receive_data<
331 : evolution::dg::subcell::Tags::InitialTciData<Dim>>(
332 : receiver_proxy[neighbor], temporal_id,
333 : std::make_pair(
334 : DirectionalId<Dim>{direction_from_neighbor, element.id()},
335 : std::move(data)));
336 : }
337 : }
338 :
339 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
340 : }
341 : };
342 :
343 : /*!
344 : * \brief Apply the TCI on the FD grid to the initial data and send the TCI
345 : * decision to neighboring elements.
346 : *
347 : * GlobalCache:
348 : * - Uses:
349 : * - `ParallelComponent` proxy
350 : *
351 : * DataBox:
352 : * - Uses:
353 : * - `domain::Tags::Element<Dim>`
354 : * - `subcell::Tags::DataForRdmpTci`
355 : * - `subcell::Tags::InitialTciData`
356 : * - `subcell::Tags::SubcellOptions`
357 : * - `subcell::Tags::ActiveGrid`
358 : * - whatever `TciOnFdGridMutator` uses
359 : * - Adds: nothing
360 : * - Removes: nothing
361 : * - Modifies:
362 : * - `subcell::Tags::DataForRdmpTci`
363 : * - `subcell::Tags::TciDecision`
364 : */
365 : template <size_t Dim, typename System, typename TciOnFdGridMutator>
366 1 : struct ComputeAndSendTciOnInitialGrid {
367 0 : using inbox_tags =
368 : tmpl::list<evolution::dg::subcell::Tags::InitialTciData<Dim>>;
369 :
370 : template <typename DbTagsList, typename... InboxTags, typename ArrayIndex,
371 : typename ActionList, typename ParallelComponent,
372 : typename Metavariables>
373 0 : static Parallel::iterable_action_return_t apply(
374 : db::DataBox<DbTagsList>& box, tuples::TaggedTuple<InboxTags...>& inboxes,
375 : Parallel::GlobalCache<Metavariables>& cache,
376 : const ArrayIndex& /*array_index*/, ActionList /*meta*/,
377 : const ParallelComponent* const /*meta*/) {
378 : const Element<Dim>& element = db::get<::domain::Tags::Element<Dim>>(box);
379 :
380 : // Check if we have received all RDMP data.
381 : if (LIKELY(element.number_of_neighbors() != 0)) {
382 : auto& inbox =
383 : tuples::get<evolution::dg::subcell::Tags::InitialTciData<Dim>>(
384 : inboxes);
385 : const auto& received = inbox.find(0);
386 : if (received == inbox.end() or
387 : received->second.size() != element.number_of_neighbors()) {
388 : return {Parallel::AlgorithmExecution::Retry, std::nullopt};
389 : }
390 :
391 : db::mutate<evolution::dg::subcell::Tags::DataForRdmpTci>(
392 : [&element, &received](const auto rdmp_tci_data_ptr) {
393 : (void)element;
394 : const size_t number_of_rdmp_vars =
395 : rdmp_tci_data_ptr->max_variables_values.size();
396 : ASSERT(rdmp_tci_data_ptr->max_variables_values.size() ==
397 : number_of_rdmp_vars,
398 : "The number of local max vars is "
399 : << number_of_rdmp_vars
400 : << " while the number of local min vars is "
401 : << rdmp_tci_data_ptr->max_variables_values.size()
402 : << " the local element ID is " << element.id());
403 : for (const auto& [direction_and_neighbor_element_id,
404 : neighbor_initial_tci_data] : received->second) {
405 : ASSERT(neighbor_initial_tci_data.initial_rdmp_data.has_value(),
406 : "Neighbor in direction "
407 : << direction_and_neighbor_element_id.direction()
408 : << " with element ID "
409 : << direction_and_neighbor_element_id.id() << " of "
410 : << element.id()
411 : << " didn't send initial TCI data correctly");
412 : ASSERT(
413 : neighbor_initial_tci_data.initial_rdmp_data.value()
414 : .max_variables_values.size() == number_of_rdmp_vars,
415 : "The number of local RDMP vars is "
416 : << number_of_rdmp_vars
417 : << " while the number of remote max vars is "
418 : << neighbor_initial_tci_data.initial_rdmp_data.value()
419 : .max_variables_values.size()
420 : << " the local element ID is " << element.id()
421 : << " and the remote id is "
422 : << direction_and_neighbor_element_id.id());
423 : ASSERT(
424 : neighbor_initial_tci_data.initial_rdmp_data.value()
425 : .min_variables_values.size() == number_of_rdmp_vars,
426 : "The number of local RDMP vars is "
427 : << number_of_rdmp_vars
428 : << " while the number of remote min vars is "
429 : << neighbor_initial_tci_data.initial_rdmp_data.value()
430 : .min_variables_values.size()
431 : << " the local element ID is " << element.id()
432 : << " and the remote id is "
433 : << direction_and_neighbor_element_id.id());
434 : for (size_t var_index = 0; var_index < number_of_rdmp_vars;
435 : ++var_index) {
436 : rdmp_tci_data_ptr->max_variables_values[var_index] =
437 : std::max(rdmp_tci_data_ptr->max_variables_values[var_index],
438 : neighbor_initial_tci_data.initial_rdmp_data.value()
439 : .max_variables_values[var_index]);
440 : rdmp_tci_data_ptr->min_variables_values[var_index] =
441 : std::min(rdmp_tci_data_ptr->min_variables_values[var_index],
442 : neighbor_initial_tci_data.initial_rdmp_data.value()
443 : .min_variables_values[var_index]);
444 : }
445 : }
446 : },
447 : make_not_null(&box));
448 : inbox.erase(received);
449 : }
450 :
451 : const auto send_tci_decision = [&cache, &element](const int tci_decision) {
452 : if (UNLIKELY(element.number_of_neighbors() == 0)) {
453 : return;
454 : }
455 : auto& receiver_proxy =
456 : Parallel::get_parallel_component<ParallelComponent>(cache);
457 : for (const auto& [direction, neighbors] : element.neighbors()) {
458 : for (const auto& neighbor : neighbors) {
459 : const auto& orientation = neighbors.orientation(neighbor);
460 : const auto direction_from_neighbor =
461 : orientation(direction.opposite());
462 : evolution::dg::subcell::InitialTciData data{tci_decision, {}};
463 : // We use temporal ID 1 for ending the TCI decision.
464 : const int temporal_id = 1;
465 : Parallel::receive_data<
466 : evolution::dg::subcell::Tags::InitialTciData<Dim>>(
467 : receiver_proxy[neighbor], temporal_id,
468 : std::make_pair(
469 : DirectionalId<Dim>{direction_from_neighbor, element.id()},
470 : std::move(data)));
471 : }
472 : }
473 : };
474 :
475 : const SubcellOptions& subcell_options =
476 : db::get<Tags::SubcellOptions<Dim>>(box);
477 :
478 : if (subcell_options.always_use_subcells() or
479 : get<Tags::ActiveGrid>(box) == ActiveGrid::Dg) {
480 : db::mutate<Tags::TciDecision>(
481 : [](const gsl::not_null<int*> tci_decision_ptr) {
482 : *tci_decision_ptr = 0;
483 : },
484 : make_not_null(&box));
485 : send_tci_decision(0);
486 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
487 : }
488 :
489 : // Now run the TCI to see if we could switch back to DG.
490 : const std::tuple<int, evolution::dg::subcell::RdmpTciData> tci_result =
491 : db::mutate_apply<TciOnFdGridMutator>(
492 : make_not_null(&box), subcell_options.persson_exponent() + 1.0,
493 : false);
494 :
495 : db::mutate<Tags::TciDecision>(
496 : [&tci_result](const gsl::not_null<int*> tci_decision_ptr) {
497 : *tci_decision_ptr = std::get<0>(tci_result);
498 : },
499 : make_not_null(&box));
500 : send_tci_decision(std::get<0>(tci_result));
501 :
502 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
503 : }
504 : };
505 :
506 : /*!
507 : * \brief Using the local and neighboring TCI decisions, switches the element to
508 : * DG if the DG solution was determined to be admissible.
509 : *
510 : * GlobalCache:
511 : * - Uses:
512 : * - `ParallelComponent` proxy
513 : *
514 : * DataBox:
515 : * - Uses:
516 : * - `domain::Tags::Element<Dim>`
517 : * - `subcell::Tags::DataForRdmpTci`
518 : * - `subcell::Tags::InitialTciData`
519 : * - `subcell::Tags::SubcellOptions`
520 : * - `subcell::Tags::ActiveGrid`
521 : * - whatever `TciOnFdGridMutator` uses
522 : * - Adds: nothing
523 : * - Removes: nothing
524 : * - Modifies:
525 : * - `subcell::Tags::NeighborTciDecisions`
526 : * - `System::variables_tag`
527 : * - `Tags::HistoryEvolvedVariables<System::variables_tag>`
528 : * - `subcell::Tags::GhostDataForReconstruction`
529 : * - `subcell::Tags::TciGridHistory`
530 : * - `subcell::Tags::CellCenteredFlux`
531 : */
532 : template <size_t Dim, typename System>
533 1 : struct SetInitialGridFromTciData {
534 : template <typename DbTagsList, typename... InboxTags, typename ArrayIndex,
535 : typename ActionList, typename ParallelComponent,
536 : typename Metavariables>
537 0 : static Parallel::iterable_action_return_t apply(
538 : db::DataBox<DbTagsList>& box, tuples::TaggedTuple<InboxTags...>& inboxes,
539 : const Parallel::GlobalCache<Metavariables>& /*cache*/,
540 : const ArrayIndex& /*array_index*/, ActionList /*meta*/,
541 : const ParallelComponent* const /*meta*/) {
542 : const Element<Dim>& element = db::get<::domain::Tags::Element<Dim>>(box);
543 : if (LIKELY(element.number_of_neighbors() != 0)) {
544 : auto& inbox =
545 : tuples::get<evolution::dg::subcell::Tags::InitialTciData<Dim>>(
546 : inboxes);
547 : const auto& received = inbox.find(1);
548 : // Check if we have received all TCI decisions.
549 : if (received == inbox.end() or
550 : received->second.size() != element.number_of_neighbors()) {
551 : return {Parallel::AlgorithmExecution::Retry, std::nullopt};
552 : }
553 :
554 : db::mutate<evolution::dg::subcell::Tags::NeighborTciDecisions<Dim>>(
555 : [&element, &received](const auto neighbor_tci_decisions_ptr) {
556 : // Non-hypercube elements (e.g. spherical shells) have an empty
557 : // NeighborTciDecisions map (see SetSubcellGrid) and will remain
558 : // on DG regardless of neighbor TCI decisions.
559 : if (neighbor_tci_decisions_ptr->empty()) {
560 : return;
561 : }
562 : for (const auto& [directional_element_id,
563 : neighbor_initial_tci_data] : received->second) {
564 : ASSERT(neighbor_initial_tci_data.tci_status.has_value(),
565 : "Neighbor in direction "
566 : << directional_element_id.direction()
567 : << " with element ID " << directional_element_id.id()
568 : << " of " << element.id()
569 : << " didn't send initial TCI decision correctly");
570 : if (not neighbor_tci_decisions_ptr->contains(
571 : directional_element_id)) {
572 : // TCI decisions for MultipleNonconforming neighbors are not
573 : // tracked because those elements are forced to remain on DG.
574 : ASSERT(element.face_types().at(
575 : directional_element_id.direction()) ==
576 : ::domain::FaceType::MultipleNonconforming,
577 : "NeighborTciDecisions does not contain the neighbor "
578 : << directional_element_id
579 : << " but the face is not MultipleNonconforming. "
580 : "This indicates a bug in the initialization of "
581 : "NeighborTciDecisions.");
582 : continue;
583 : }
584 : neighbor_tci_decisions_ptr->at(directional_element_id) =
585 : neighbor_initial_tci_data.tci_status.value();
586 : }
587 : },
588 : make_not_null(&box));
589 : inbox.erase(received);
590 : }
591 :
592 : if (get<Tags::ActiveGrid>(box) == ActiveGrid::Dg) {
593 : // In this case we are allowed to only do DG in this element. No need to
594 : // even do any checks.
595 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
596 : }
597 :
598 : const SubcellOptions& subcell_options =
599 : db::get<Tags::SubcellOptions<Dim>>(box);
600 :
601 : bool cell_is_troubled =
602 : subcell_options.always_use_subcells() or
603 : (subcell_options.use_halo() and [&box]() -> bool {
604 : for (const auto& [_, neighbor_decision] :
605 : db::get<evolution::dg::subcell::Tags::NeighborTciDecisions<Dim>>(
606 : box)) {
607 : if (neighbor_decision != 0) {
608 : return true;
609 : }
610 : }
611 : return false;
612 : }()) or
613 : (db::get<Tags::TciDecision>(box) != 0);
614 :
615 : if (not cell_is_troubled) {
616 : using variables_tag = typename System::variables_tag;
617 : using flux_variables = typename System::flux_variables;
618 :
619 : const Mesh<Dim>& dg_mesh = db::get<::domain::Tags::Mesh<Dim>>(box);
620 : const Mesh<Dim>& subcell_mesh = db::get<subcell::Tags::Mesh<Dim>>(box);
621 : db::mutate<
622 : variables_tag, ::Tags::HistoryEvolvedVariables<variables_tag>,
623 : Tags::ActiveGrid, subcell::Tags::GhostDataForReconstruction<Dim>,
624 : evolution::dg::subcell::Tags::TciGridHistory,
625 : evolution::dg::subcell::Tags::CellCenteredFlux<flux_variables, Dim>>(
626 : [&dg_mesh, &subcell_mesh, &subcell_options](
627 : const auto active_vars_ptr, const auto active_history_ptr,
628 : const gsl::not_null<ActiveGrid*> active_grid_ptr,
629 : const auto subcell_ghost_data_ptr,
630 : const gsl::not_null<
631 : std::deque<evolution::dg::subcell::ActiveGrid>*>
632 : tci_grid_history_ptr,
633 : const auto subcell_cell_centered_fluxes) {
634 : // Note: strictly speaking, to be conservative this should
635 : // reconstruct uJ instead of u.
636 : *active_vars_ptr = fd::reconstruct(
637 : *active_vars_ptr, dg_mesh, subcell_mesh.extents(),
638 : subcell_options.reconstruction_method());
639 :
640 : // Reconstruct the DG solution for each time in the time stepper
641 : // history
642 : active_history_ptr->map_entries(
643 : [&dg_mesh, &subcell_mesh, &subcell_options](const auto entry) {
644 : *entry =
645 : fd::reconstruct(*entry, dg_mesh, subcell_mesh.extents(),
646 : subcell_options.reconstruction_method());
647 : });
648 : *active_grid_ptr = ActiveGrid::Dg;
649 :
650 : // Clear the neighbor data needed for subcell reconstruction since
651 : // we have now completed the time step.
652 : subcell_ghost_data_ptr->clear();
653 :
654 : // Clear the TCI grid history since we don't need to use it when on
655 : // the DG grid.
656 : tci_grid_history_ptr->clear();
657 :
658 : // Clear the allocation for the cell-centered fluxes.
659 : *subcell_cell_centered_fluxes = std::nullopt;
660 : },
661 : make_not_null(&box));
662 : }
663 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
664 : }
665 : };
666 : } // namespace evolution::dg::subcell::Actions
|
