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 std::binary_search(
177 : subcell_options.only_dg_block_ids().begin(),
178 : subcell_options.only_dg_block_ids().end(),
179 : db::get<domain::Tags::Element<Dim>>(box).id().block_id()),
180 : "Should never use subcell on element "
181 : << db::get<domain::Tags::Element<Dim>>(box).id());
182 :
183 : // The first condition is that for substep time integrators we only allow
184 : // switching back to DG on step boundaries. This is the easiest way to
185 : // avoid having a shock in the time stepper history, since there is no
186 : // history at step boundaries.
187 : //
188 : // The second condition is that if we are in the self-start procedure of
189 : // the time stepper, and we don't want to switch from subcell back to DG
190 : // during self-start since we integrate over the same temporal region at
191 : // increasingly higher order.
192 : //
193 : // The third condition is that the user has requested we always do
194 : // subcell, so effectively a finite difference/volume code.
195 : const bool only_need_rdmp_data =
196 : db::get<subcell::Tags::DidRollback>(box) or
197 : (time_step_id.substep() != 0 or time_step_id.slab_number() < 0)
198 :
199 : or db::get<evolution::dg::subcell::Tags::StepsSinceTciCall>(box) + 1 <
200 : subcell_options.number_of_steps_between_tci_calls();
201 : if (UNLIKELY(db::get<subcell::Tags::DidRollback>(box))) {
202 : db::mutate<subcell::Tags::DidRollback>(
203 : [](const gsl::not_null<bool*> did_rollback) {
204 : *did_rollback = false;
205 : },
206 : make_not_null(&box));
207 : }
208 :
209 : if (subcell_options.always_use_subcells()) {
210 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
211 : }
212 :
213 : std::tuple<int, evolution::dg::subcell::RdmpTciData> tci_result =
214 : db::mutate_apply<TciMutator>(make_not_null(&box),
215 : subcell_options.persson_exponent() + 1.0,
216 : only_need_rdmp_data);
217 :
218 : db::mutate<evolution::dg::subcell::Tags::DataForRdmpTci,
219 : evolution::dg::subcell::Tags::TciCallsSinceRollback,
220 : evolution::dg::subcell::Tags::StepsSinceTciCall>(
221 : [only_need_rdmp_data, &subcell_options, &tci_result](
222 : const auto rdmp_data_ptr,
223 : const gsl::not_null<size_t*> tci_calls_since_rollback_ptr,
224 : const gsl::not_null<size_t*> steps_since_tci_call_ptr) {
225 : *rdmp_data_ptr = std::move(std::get<1>(std::move(tci_result)));
226 : (*tci_calls_since_rollback_ptr) += (only_need_rdmp_data ? 0 : 1);
227 : if ((*steps_since_tci_call_ptr) + 1 <
228 : subcell_options.number_of_steps_between_tci_calls()) {
229 : (*steps_since_tci_call_ptr) += 1;
230 : } else {
231 : (*steps_since_tci_call_ptr) = 0;
232 : }
233 : },
234 : make_not_null(&box));
235 :
236 : // Use <= for TciCallsSinceRollback since we've already incremented the
237 : // count in a mutate call above. This means that if this is the first TCI
238 : // call after a rollback, TciCallsSinceRollback will be 1, not 0. We also
239 : // require that `subcell_options.min_tci_calls_after_rollback() >= 1`, so
240 : // 1 TCI call means only one step was taken after a rollback.
241 : if (only_need_rdmp_data or
242 : db::get<evolution::dg::subcell::Tags::TciCallsSinceRollback>(box) <=
243 : subcell_options.min_tci_calls_after_rollback()) {
244 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
245 : }
246 :
247 : const int tci_decision = std::get<0>(tci_result);
248 : const bool cell_is_troubled =
249 : tci_decision != 0 or (subcell_options.use_halo() and [&box]() -> bool {
250 : for (const auto& [_, neighbor_decision] :
251 : db::get<evolution::dg::subcell::Tags::NeighborTciDecisions<Dim>>(
252 : box)) {
253 : if (neighbor_decision != 0) {
254 : return true;
255 : }
256 : }
257 : return false;
258 : }());
259 :
260 : db::mutate<Tags::TciDecision>(
261 : [&tci_decision](const gsl::not_null<int*> tci_decision_ptr) {
262 : *tci_decision_ptr = tci_decision;
263 : },
264 : make_not_null(&box));
265 :
266 : // If the cell is not troubled, then we _might_ be able to switch back to
267 : // DG. This depends on the type of time stepper we are using:
268 : // - ADER: Not yet implemented, but here the TCI history is irrelevant
269 : // because it is a one-step scheme, so we can always switch.
270 : // - Multistep: the TCI must have marked the entire evolved variable history
271 : // as free of shocks. In practice for LMMs this means the TCI
272 : // history is as long as the evolved variables history and that
273 : // the entire TCI history is `ActiveGrid::Dg`.
274 : // - Substep: the easiest is to restrict switching back to DG to step
275 : // boundaries where there is no history.
276 : const auto& time_stepper = db::get<::Tags::TimeStepper<TimeStepper>>(box);
277 : const bool is_substep_method = time_stepper.number_of_substeps() != 1;
278 : const bool is_multistep_method = time_stepper.number_of_past_steps() != 0;
279 : ASSERT(time_stepper.number_of_substeps() != 0,
280 : "Don't know how to handle a time stepper with zero substeps. This "
281 : "might be totally fine, but the case should be thought about.");
282 : ASSERT(subcell_options.min_clear_tci_before_dg() > 0,
283 : "The value of 'subcell_options.min_clear_tci_before_dg()' must be "
284 : "at least 1");
285 : if (const auto& tci_history = db::get<subcell::Tags::TciGridHistory>(box);
286 : not cell_is_troubled and
287 :
288 : (((is_substep_method and
289 : tci_history.size() >= subcell_options.min_clear_tci_before_dg() - 1)
290 :
291 : or
292 :
293 : (is_multistep_method and
294 : tci_history.size() ==
295 : std::max(subcell_options.min_clear_tci_before_dg(),
296 : // Use `number_of_past_steps()` instead of
297 : // `number_of_past_steps()+number_of_substeps()`
298 : // since we only perform this check if we are on
299 : // substep 0 (i.e. starting a new step).
300 : time_stepper.number_of_past_steps() + 1)))
301 :
302 : and alg::all_of(tci_history,
303 : [](const ActiveGrid tci_grid_decision) {
304 : return tci_grid_decision == ActiveGrid::Dg;
305 : })
306 :
307 : )) {
308 : db::mutate<
309 : variables_tag, ::Tags::HistoryEvolvedVariables<variables_tag>,
310 : Tags::ActiveGrid, subcell::Tags::GhostDataForReconstruction<Dim>,
311 : evolution::dg::subcell::Tags::TciGridHistory,
312 : evolution::dg::subcell::Tags::TciCallsSinceRollback,
313 : evolution::dg::subcell::Tags::CellCenteredFlux<flux_variables, Dim>>(
314 : [&dg_mesh, &subcell_mesh, &subcell_options](
315 : const auto active_vars_ptr, const auto active_history_ptr,
316 : const gsl::not_null<ActiveGrid*> active_grid_ptr,
317 : const auto subcell_ghost_data_ptr,
318 : const gsl::not_null<
319 : std::deque<evolution::dg::subcell::ActiveGrid>*>
320 : tci_grid_history_ptr,
321 : const gsl::not_null<size_t*> tci_calls_since_rollback_ptr,
322 : const auto subcell_cell_centered_fluxes) {
323 : // Note: strictly speaking, to be conservative this should
324 : // reconstruct uJ instead of u.
325 : *active_vars_ptr = fd::reconstruct(
326 : *active_vars_ptr, dg_mesh, subcell_mesh.extents(),
327 : subcell_options.reconstruction_method());
328 :
329 : // Reconstruct the DG solution for each time in the time stepper
330 : // history
331 : active_history_ptr->map_entries(
332 : [&dg_mesh, &subcell_mesh, &subcell_options](const auto entry) {
333 : *entry =
334 : fd::reconstruct(*entry, dg_mesh, subcell_mesh.extents(),
335 : subcell_options.reconstruction_method());
336 : });
337 : *active_grid_ptr = ActiveGrid::Dg;
338 :
339 : // Clear the neighbor data needed for subcell reconstruction since
340 : // we have now completed the time step.
341 : subcell_ghost_data_ptr->clear();
342 :
343 : // Clear the TCI grid history since we don't need to use it when on
344 : // the DG grid.
345 : tci_grid_history_ptr->clear();
346 :
347 : // Reset tci_calls_since_rollback
348 : *tci_calls_since_rollback_ptr = 0;
349 :
350 : // Clear the allocation for the cell-centered fluxes.
351 : *subcell_cell_centered_fluxes = std::nullopt;
352 : },
353 : make_not_null(&box));
354 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
355 : }
356 :
357 : // We record the TCI history for a couple reasons.
358 : //
359 : // For multistep methods we track the TCI decision because we need the
360 : // discontinuity to clear the entire history before we can switch back to
361 : // DG.
362 : //
363 : // For substep methods tracking the TCI history allows us to stay on FD a
364 : // bit longer to allow the solution to smoothen a bit more before going
365 : // back to DG.
366 : db::mutate<evolution::dg::subcell::Tags::TciGridHistory>(
367 : [cell_is_troubled, &subcell_options, &time_stepper](
368 : const gsl::not_null<std::deque<evolution::dg::subcell::ActiveGrid>*>
369 : tci_grid_history) {
370 : tci_grid_history->push_front(cell_is_troubled ? ActiveGrid::Subcell
371 : : ActiveGrid::Dg);
372 : if (tci_grid_history->size() >
373 : std::max(subcell_options.min_clear_tci_before_dg() - 1,
374 : // Use `number_of_past_steps()` instead of
375 : // `number_of_past_steps()+number_of_substeps()`
376 : // since we only perform this check if we are on
377 : // substep 0 (i.e. starting a new step).
378 : time_stepper.number_of_past_steps() + 1)) {
379 : tci_grid_history->pop_back();
380 : }
381 : ASSERT(tci_grid_history->size() <=
382 : std::max(subcell_options.min_clear_tci_before_dg() - 1,
383 : time_stepper.number_of_past_steps() + 1),
384 : "The TCI history is not being correctly cleared. This is an "
385 : "internal bug. Please file an issue.\n"
386 : "tci_grid_history->size(): "
387 : << tci_grid_history->size() << "\nExpected: "
388 : << std::max(subcell_options.min_clear_tci_before_dg() - 1,
389 : time_stepper.number_of_past_steps() + 1));
390 : },
391 : make_not_null(&box));
392 :
393 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
394 : }
395 : };
396 : } // namespace evolution::dg::subcell::Actions
|
