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