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