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