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 <deque>
8 : #include <optional>
9 : #include <tuple>
10 : #include <type_traits>
11 : #include <utility>
12 :
13 : #include "DataStructures/DataBox/DataBox.hpp"
14 : #include "DataStructures/DataBox/PrefixHelpers.hpp"
15 : #include "DataStructures/DataBox/Prefixes.hpp"
16 : #include "DataStructures/Variables.hpp"
17 : #include "DataStructures/VariablesTag.hpp"
18 : #include "Domain/Structure/Direction.hpp"
19 : #include "Domain/Structure/ElementId.hpp"
20 : #include "Domain/Tags.hpp"
21 : #include "Evolution/DgSubcell/ActiveGrid.hpp"
22 : #include "Evolution/DgSubcell/RdmpTci.hpp"
23 : #include "Evolution/DgSubcell/RdmpTciData.hpp"
24 : #include "Evolution/DgSubcell/Reconstruction.hpp"
25 : #include "Evolution/DgSubcell/ReconstructionMethod.hpp"
26 : #include "Evolution/DgSubcell/SubcellOptions.hpp"
27 : #include "Evolution/DgSubcell/Tags/ActiveGrid.hpp"
28 : #include "Evolution/DgSubcell/Tags/CellCenteredFlux.hpp"
29 : #include "Evolution/DgSubcell/Tags/DataForRdmpTci.hpp"
30 : #include "Evolution/DgSubcell/Tags/DidRollback.hpp"
31 : #include "Evolution/DgSubcell/Tags/GhostDataForReconstruction.hpp"
32 : #include "Evolution/DgSubcell/Tags/Mesh.hpp"
33 : #include "Evolution/DgSubcell/Tags/StepsSinceTciCall.hpp"
34 : #include "Evolution/DgSubcell/Tags/SubcellOptions.hpp"
35 : #include "Evolution/DgSubcell/Tags/TciCallsSinceRollback.hpp"
36 : #include "Evolution/DgSubcell/Tags/TciGridHistory.hpp"
37 : #include "Evolution/DgSubcell/Tags/TciStatus.hpp"
38 : #include "NumericalAlgorithms/Spectral/Mesh.hpp"
39 : #include "Parallel/AlgorithmExecution.hpp"
40 : #include "Time/History.hpp"
41 : #include "Time/TimeStepId.hpp"
42 : #include "Utilities/ContainerHelpers.hpp"
43 : #include "Utilities/ErrorHandling/Assert.hpp"
44 : #include "Utilities/Gsl.hpp"
45 : #include "Utilities/TMPL.hpp"
46 :
47 : /// \cond
48 : namespace Parallel {
49 : template <typename Metavariables>
50 : class GlobalCache;
51 : } // namespace Parallel
52 : namespace Tags {
53 : template <typename Tag>
54 : struct HistoryEvolvedVariables;
55 : struct TimeStepId;
56 : template <typename StepperInterface>
57 : struct TimeStepper;
58 : } // namespace Tags
59 : class TimeStepper;
60 : namespace tuples {
61 : template <typename...>
62 : class TaggedTuple;
63 : } // namespace tuples
64 : /// \endcond
65 :
66 : namespace evolution::dg::subcell::Actions {
67 : /*!
68 : * \brief Run the troubled-cell indicator on the subcell solution to see if it
69 : * is safe to switch back to DG.
70 : *
71 : * In terms of the DG-subcell/FD hybrid solver, this action is run after the FD
72 : * step has calculated the solution at \f$t^{n+1}\f$. At this point we check if
73 : * the FD solution at the new time \f$t^{n+1}\f$ is representable on the DG
74 : * grid.
75 : *
76 : * The algorithm proceeds as follows:
77 : * 1. If we are using a substep time integrator and are not at the end of a
78 : * step, or we are in the self-starting stage of a multistep method, or the
79 : * `subcell_options.always_use_subcells() == true`, then we do not run any
80 : * TCI or try to go back to DG. We need to avoid reconstructing (in the sense
81 : * of the inverse of projecting the DG solution to the subcells) the time
82 : * stepper history if there are shocks present in the history, and for
83 : * substep methods this is most easily handled by only switching back at the
84 : * end of a full time step. During the self-start phase of the multistep time
85 : * integrators we integrate over the same region of time at increasingly
86 : * higher order, which means if we were on subcell "previously" (since we use
87 : * a forward-in-time self-start method the time history is actually in the
88 : * future of the current step) then we will very likely need to again switch
89 : * to subcell.
90 : * 2. Reconstruct the subcell solution to the DG grid.
91 : * 3. Run the relaxed discrete maximum principle (RDMP) troubled-cell indicator
92 : * (TCI), checking both the subcell solution at \f$t^{n+1}\f$ and the
93 : * reconstructed DG solution at \f$t^{n+1}\f$ to make sure they are
94 : * admissible.
95 : * 4. If the RDMP TCI marked the DG solution as admissible, run the
96 : * user-specified mutator TCI `TciMutator`.
97 : * 5. If the cell is not troubled, and the time integrator type is substep or
98 : * the TCI history indicates the entire history for the multistep method is
99 : * free of discontinuities, then we can switch back to DG. Switching back to
100 : * DG requires reconstructing dg variables, reconstructing the time stepper
101 : * history, marking the active grid as `ActiveGrid::Dg`, and clearing the
102 : * subcell neighbor data.
103 : * 6. If we are not using a substep method, then record the TCI decision in the
104 : * `subcell::Tags::TciGridHistory`.
105 : *
106 : * \note Unlike `Actions::TciAndRollback`, this action does _not_ jump back to
107 : * `Labels::BeginDg`. This is because users may add actions after a time step
108 : * has been completed. In that sense, it may be more proper to actually check
109 : * the TCI and switch back to DG at the start of the step rather than the end.
110 : *
111 : * \note This action always sets `subcell::Tags::DidRollback` to `false` at the
112 : * very beginning since this action is called after an FD step has completed.
113 : *
114 : * \note If `subcell::Tags::DidRollback=True` i.e. the grid has been switched
115 : * from DG to FD in the current time step by preceding actions, this action is
116 : * skipped except setting `DidRollback` to `false`. Stated differently, if an
117 : * element switched from DG to FD it needs to remain at least one time step on
118 : * the FD grid.
119 : *
120 : * GlobalCache:
121 : * - Uses:
122 : * - `subcell::Tags::SubcellOptions`
123 : *
124 : * DataBox:
125 : * - Uses:
126 : * - `domain::Tags::Mesh<Dim>`
127 : * - `subcell::Tags::Mesh<Dim>`
128 : * - `Tags::TimeStepId`
129 : * - `subcell::Tags::ActiveGrid`
130 : * - `subcell::Tags::DataForRdmpTci`
131 : * - `subcell::Tags::TciGridHistory`
132 : * - Adds: nothing
133 : * - Removes: nothing
134 : * - Modifies:
135 : * - `System::variables_tag` if the cell is not troubled
136 : * - `Tags::HistoryEvolvedVariables` if the cell is not troubled
137 : * - `subcell::Tags::ActiveGrid` if the cell is not troubled
138 : * - `subcell::Tags::DidRollback` sets to `false`
139 : * - `subcell::Tags::TciDecision` is set to an integer value according to the
140 : * return of TciMutator.
141 : * - `subcell::Tags::GhostDataForReconstruction<Dim>`
142 : * if the cell is not troubled
143 : * - `subcell::Tags::TciGridHistory` if the time stepper is a multistep method
144 : */
145 : template <typename TciMutator>
146 1 : struct TciAndSwitchToDg {
147 : template <typename DbTags, typename... InboxTags, typename Metavariables,
148 : typename ArrayIndex, typename ActionList,
149 : typename ParallelComponent, size_t Dim = Metavariables::volume_dim>
150 0 : static Parallel::iterable_action_return_t apply(
151 : db::DataBox<DbTags>& box,
152 : const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
153 : const Parallel::GlobalCache<Metavariables>& /*cache*/,
154 : const ArrayIndex& /*array_index*/, const ActionList /*meta*/,
155 : const ParallelComponent* const /*meta*/) {
156 : static_assert(
157 : tmpl::count_if<
158 : ActionList,
159 : std::is_same<tmpl::_1, tmpl::pin<TciAndSwitchToDg>>>::value == 1,
160 : "Must have the TciAndSwitchToDg action exactly once in the action list "
161 : "of a phase.");
162 :
163 : using variables_tag = typename Metavariables::system::variables_tag;
164 : using flux_variables = typename Metavariables::system::flux_variables;
165 :
166 : ASSERT(db::get<Tags::ActiveGrid>(box) == ActiveGrid::Subcell,
167 : "Must be using subcells when calling TciAndSwitchToDg action.");
168 : const Mesh<Dim>& dg_mesh = db::get<::domain::Tags::Mesh<Dim>>(box);
169 : const Mesh<Dim>& subcell_mesh = db::get<subcell::Tags::Mesh<Dim>>(box);
170 : const SubcellOptions& subcell_options =
171 : db::get<Tags::SubcellOptions<Dim>>(box);
172 : const TimeStepId& time_step_id = db::get<::Tags::TimeStepId>(box);
173 :
174 : // This should never be run if we are prohibited from using subcell on this
175 : // element.
176 : ASSERT(not subcell_options.only_dg_block_ids().contains(
177 : db::get<domain::Tags::Element<Dim>>(box).id().block_id()),
178 : "Should never use subcell on element "
179 : << db::get<domain::Tags::Element<Dim>>(box).id());
180 :
181 : // The first condition is that for substep time integrators we only allow
182 : // switching back to DG on step boundaries. This is the easiest way to
183 : // avoid having a shock in the time stepper history, since there is no
184 : // history at step boundaries.
185 : //
186 : // The second condition is that if we are in the self-start procedure of
187 : // the time stepper, and we don't want to switch from subcell back to DG
188 : // during self-start since we integrate over the same temporal region at
189 : // increasingly higher order.
190 : //
191 : // The third condition is that the user has requested we always do
192 : // subcell, so effectively a finite difference/volume code.
193 : const bool only_need_rdmp_data =
194 : db::get<subcell::Tags::DidRollback>(box) or
195 : (time_step_id.substep() != 0 or time_step_id.slab_number() < 0)
196 :
197 : or db::get<evolution::dg::subcell::Tags::StepsSinceTciCall>(box) + 1 <
198 : subcell_options.number_of_steps_between_tci_calls();
199 : if (UNLIKELY(db::get<subcell::Tags::DidRollback>(box))) {
200 : db::mutate<subcell::Tags::DidRollback>(
201 : [](const gsl::not_null<bool*> did_rollback) {
202 : *did_rollback = false;
203 : },
204 : make_not_null(&box));
205 : }
206 :
207 : if (subcell_options.always_use_subcells()) {
208 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
209 : }
210 :
211 : std::tuple<int, evolution::dg::subcell::RdmpTciData> tci_result =
212 : db::mutate_apply<TciMutator>(make_not_null(&box),
213 : subcell_options.persson_exponent() + 1.0,
214 : only_need_rdmp_data);
215 :
216 : db::mutate<evolution::dg::subcell::Tags::DataForRdmpTci,
217 : evolution::dg::subcell::Tags::TciCallsSinceRollback,
218 : evolution::dg::subcell::Tags::StepsSinceTciCall>(
219 : [only_need_rdmp_data, &subcell_options, &tci_result](
220 : const auto rdmp_data_ptr,
221 : const gsl::not_null<size_t*> tci_calls_since_rollback_ptr,
222 : const gsl::not_null<size_t*> steps_since_tci_call_ptr) {
223 : *rdmp_data_ptr = std::move(std::get<1>(std::move(tci_result)));
224 : (*tci_calls_since_rollback_ptr) += (only_need_rdmp_data ? 0 : 1);
225 : if ((*steps_since_tci_call_ptr) + 1 <
226 : subcell_options.number_of_steps_between_tci_calls()) {
227 : (*steps_since_tci_call_ptr) += 1;
228 : } else {
229 : (*steps_since_tci_call_ptr) = 0;
230 : }
231 : },
232 : make_not_null(&box));
233 :
234 : // Use <= for TciCallsSinceRollback since we've already incremented the
235 : // count in a mutate call above. This means that if this is the first TCI
236 : // call after a rollback, TciCallsSinceRollback will be 1, not 0. We also
237 : // require that `subcell_options.min_tci_calls_after_rollback() >= 1`, so
238 : // 1 TCI call means only one step was taken after a rollback.
239 : if (only_need_rdmp_data or
240 : db::get<evolution::dg::subcell::Tags::TciCallsSinceRollback>(box) <=
241 : subcell_options.min_tci_calls_after_rollback()) {
242 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
243 : }
244 :
245 : const int tci_decision = std::get<0>(tci_result);
246 : const bool cell_is_troubled =
247 : tci_decision != 0 or (subcell_options.use_halo() and [&box]() -> bool {
248 : for (const auto& [_, neighbor_decision] :
249 : db::get<evolution::dg::subcell::Tags::NeighborTciDecisions<Dim>>(
250 : box)) {
251 : if (neighbor_decision != 0) {
252 : return true;
253 : }
254 : }
255 : return false;
256 : }());
257 :
258 : db::mutate<Tags::TciDecision>(
259 : [&tci_decision](const gsl::not_null<int*> tci_decision_ptr) {
260 : *tci_decision_ptr = tci_decision;
261 : },
262 : make_not_null(&box));
263 :
264 : // If the cell is not troubled, then we _might_ be able to switch back to
265 : // DG. This depends on the type of time stepper we are using:
266 : // - ADER: Not yet implemented, but here the TCI history is irrelevant
267 : // because it is a one-step scheme, so we can always switch.
268 : // - Multistep: the TCI must have marked the entire evolved variable history
269 : // as free of shocks. In practice for LMMs this means the TCI
270 : // history is as long as the evolved variables history and that
271 : // the entire TCI history is `ActiveGrid::Dg`.
272 : // - Substep: the easiest is to restrict switching back to DG to step
273 : // boundaries where there is no history.
274 : const auto& time_stepper = db::get<::Tags::TimeStepper<TimeStepper>>(box);
275 : const bool is_substep_method = time_stepper.number_of_substeps() != 1;
276 : const bool is_multistep_method = time_stepper.number_of_past_steps() != 0;
277 : ASSERT(time_stepper.number_of_substeps() != 0,
278 : "Don't know how to handle a time stepper with zero substeps. This "
279 : "might be totally fine, but the case should be thought about.");
280 : ASSERT(subcell_options.min_clear_tci_before_dg() > 0,
281 : "The value of 'subcell_options.min_clear_tci_before_dg()' must be "
282 : "at least 1");
283 : if (const auto& tci_history = db::get<subcell::Tags::TciGridHistory>(box);
284 : not cell_is_troubled and
285 :
286 : (((is_substep_method and
287 : tci_history.size() >= subcell_options.min_clear_tci_before_dg() - 1)
288 :
289 : or
290 :
291 : (is_multistep_method and
292 : tci_history.size() ==
293 : std::max(subcell_options.min_clear_tci_before_dg(),
294 : // Use `number_of_past_steps()` instead of
295 : // `number_of_past_steps()+number_of_substeps()`
296 : // since we only perform this check if we are on
297 : // substep 0 (i.e. starting a new step).
298 : time_stepper.number_of_past_steps() + 1)))
299 :
300 : and alg::all_of(tci_history,
301 : [](const ActiveGrid tci_grid_decision) {
302 : return tci_grid_decision == ActiveGrid::Dg;
303 : })
304 :
305 : )) {
306 : db::mutate<
307 : variables_tag, ::Tags::HistoryEvolvedVariables<variables_tag>,
308 : Tags::ActiveGrid, subcell::Tags::GhostDataForReconstruction<Dim>,
309 : evolution::dg::subcell::Tags::TciGridHistory,
310 : evolution::dg::subcell::Tags::TciCallsSinceRollback,
311 : evolution::dg::subcell::Tags::CellCenteredFlux<flux_variables, Dim>>(
312 : [&dg_mesh, &subcell_mesh, &subcell_options](
313 : const auto active_vars_ptr, const auto active_history_ptr,
314 : const gsl::not_null<ActiveGrid*> active_grid_ptr,
315 : const auto subcell_ghost_data_ptr,
316 : const gsl::not_null<
317 : std::deque<evolution::dg::subcell::ActiveGrid>*>
318 : tci_grid_history_ptr,
319 : const gsl::not_null<size_t*> tci_calls_since_rollback_ptr,
320 : const auto subcell_cell_centered_fluxes) {
321 : // Note: strictly speaking, to be conservative this should
322 : // reconstruct uJ instead of u.
323 : *active_vars_ptr = fd::reconstruct(
324 : *active_vars_ptr, dg_mesh, subcell_mesh.extents(),
325 : subcell_options.reconstruction_method());
326 :
327 : // Reconstruct the DG solution for each time in the time stepper
328 : // history
329 : active_history_ptr->map_entries(
330 : [&dg_mesh, &subcell_mesh, &subcell_options](const auto entry) {
331 : *entry =
332 : fd::reconstruct(*entry, dg_mesh, subcell_mesh.extents(),
333 : subcell_options.reconstruction_method());
334 : });
335 : *active_grid_ptr = ActiveGrid::Dg;
336 :
337 : // Clear the neighbor data needed for subcell reconstruction since
338 : // we have now completed the time step.
339 : subcell_ghost_data_ptr->clear();
340 :
341 : // Clear the TCI grid history since we don't need to use it when on
342 : // the DG grid.
343 : tci_grid_history_ptr->clear();
344 :
345 : // Reset tci_calls_since_rollback
346 : *tci_calls_since_rollback_ptr = 0;
347 :
348 : // Clear the allocation for the cell-centered fluxes.
349 : *subcell_cell_centered_fluxes = std::nullopt;
350 : },
351 : make_not_null(&box));
352 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
353 : }
354 :
355 : // We record the TCI history for a couple reasons.
356 : //
357 : // For multistep methods we track the TCI decision because we need the
358 : // discontinuity to clear the entire history before we can switch back to
359 : // DG.
360 : //
361 : // For substep methods tracking the TCI history allows us to stay on FD a
362 : // bit longer to allow the solution to smoothen a bit more before going
363 : // back to DG.
364 : db::mutate<evolution::dg::subcell::Tags::TciGridHistory>(
365 : [cell_is_troubled, &subcell_options, &time_stepper](
366 : const gsl::not_null<std::deque<evolution::dg::subcell::ActiveGrid>*>
367 : tci_grid_history) {
368 : tci_grid_history->push_front(cell_is_troubled ? ActiveGrid::Subcell
369 : : ActiveGrid::Dg);
370 : if (tci_grid_history->size() >
371 : std::max(subcell_options.min_clear_tci_before_dg() - 1,
372 : // Use `number_of_past_steps()` instead of
373 : // `number_of_past_steps()+number_of_substeps()`
374 : // since we only perform this check if we are on
375 : // substep 0 (i.e. starting a new step).
376 : time_stepper.number_of_past_steps() + 1)) {
377 : tci_grid_history->pop_back();
378 : }
379 : ASSERT(tci_grid_history->size() <=
380 : std::max(subcell_options.min_clear_tci_before_dg() - 1,
381 : time_stepper.number_of_past_steps() + 1),
382 : "The TCI history is not being correctly cleared. This is an "
383 : "internal bug. Please file an issue.\n"
384 : "tci_grid_history->size(): "
385 : << tci_grid_history->size() << "\nExpected: "
386 : << std::max(subcell_options.min_clear_tci_before_dg() - 1,
387 : time_stepper.number_of_past_steps() + 1));
388 : },
389 : make_not_null(&box));
390 :
391 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
392 : }
393 : };
394 : } // namespace evolution::dg::subcell::Actions
|
