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