TciAndSwitchToDg.hpp
1 // Distributed under the MIT License.
2 // See LICENSE.txt for details.
3 
4 #pragma once
5 
6 #include <cstddef>
7 #include <deque>
8 #include <tuple>
9 #include <type_traits>
10 #include <utility>
11 
13 #include "DataStructures/DataBox/PrefixHelpers.hpp"
16 #include "DataStructures/VariablesTag.hpp"
19 #include "Domain/Tags.hpp"
20 #include "Evolution/DgSubcell/ActiveGrid.hpp"
21 #include "Evolution/DgSubcell/NeighborData.hpp"
22 #include "Evolution/DgSubcell/RdmpTci.hpp"
23 #include "Evolution/DgSubcell/Reconstruction.hpp"
24 #include "Evolution/DgSubcell/SubcellOptions.hpp"
25 #include "Evolution/DgSubcell/Tags/ActiveGrid.hpp"
26 #include "Evolution/DgSubcell/Tags/DidRollback.hpp"
27 #include "Evolution/DgSubcell/Tags/Inactive.hpp"
28 #include "Evolution/DgSubcell/Tags/Mesh.hpp"
29 #include "Evolution/DgSubcell/Tags/NeighborData.hpp"
30 #include "Evolution/DgSubcell/Tags/SubcellOptions.hpp"
31 #include "Evolution/DgSubcell/Tags/TciGridHistory.hpp"
33 #include "Time/History.hpp"
34 #include "Time/Tags.hpp"
35 #include "Time/TimeStepId.hpp"
36 #include "Time/TimeSteppers/TimeStepper.hpp"
38 #include "Utilities/Gsl.hpp"
39 #include "Utilities/TMPL.hpp"
40 
41 /// \cond
42 namespace Parallel {
43 template <typename Metavariables>
44 class GlobalCache;
45 } // namespace Parallel
46 namespace tuples {
47 template <typename...>
48 class TaggedTuple;
49 } // namespace tuples
50 /// \endcond
51 
53 /*!
54  * \brief Run the troubled-cell indicator on the subcell solution to see if it
55  * is safe to switch back to DG.
56  *
57  * In terms of the DG-subcell/FD hybrid solver, this action is run after the FD
58  * step has calculated the solution at \f$t^{n+1}\f$. At this point we check if
59  * the FD solution at the new time \f$t^{n+1}\f$ is representable on the DG
60  * grid.
61  *
62  * The algorithm proceeds as follows:
63  * 1. If we are using a substep time integrator and are not at the end of a
64  * step, or we are in the self-starting stage of a multistep method, or the
65  * `subcell_options.always_use_subcells() == true`, then we do not run any
66  * TCI or try to go back to DG. We need to avoid reconstructing (in the sense
67  * of the inverse of projecting the DG solution to the subcells) the time
68  * stepper history if there are shocks present in the history, and for
69  * substep methods this is most easily handled by only switching back at the
70  * end of a full time step. During the self-start phase of the multistep time
71  * integrators we integrate over the same region of time at increasingly
72  * higher order, which means if we were on subcell "previously" (since we use
73  * a forward-in-time self-start method the time history is actually in the
74  * future of the current step) then we will very likely need to again switch
75  * to subcell.
76  * 2. Reconstruct the subcell solution to the DG grid.
77  * 3. Run the relaxed discrete maximum principle (RDMP) troubled-cell indicator
78  * (TCI), checking both the subcell solution at \f$t^{n+1}\f$ and the
79  * reconstructed DG solution at \f$t^{n+1}\f$ to make sure they are
80  * admissible.
81  * 4. If the RDMP TCI marked the DG solution as admissible, run the
82  * user-specified mutator TCI `TciMutator`.
83  * 5. If the cell is not troubled, and the time integrator type is substep or
84  * the TCI history indicates the entire history for the multistep method is
85  * free of discontinuities, then we can switch back to DG. Switching back to
86  * DG requires swapping the active and inactive evolved variables,
87  * reconstructing the time stepper history, marking the active grid as
88  * `ActiveGrid::Dg`, and clearing the subcell neighbor data.
89  * 6. If we are not using a substep method, then record the TCI decision in the
90  * `subcell::Tags::TciGridHistory`.
91  *
92  * \note Unlike `Actions::TciAndRollback`, this action does _not_ jump back to
93  * `Labels::BeginDg`. This is because users may add actions after a time step
94  * has been completed. In that sense, it may be more proper to actually check
95  * the TCI and switch back to DG at the start of the step rather than the end.
96  *
97  * \note This action always sets `subcell::Tags::DidRollback` to `false` at the
98  * very beginning since this action is called after an FD step has completed.
99  *
100  * GlobalCache:
101  * - Uses:
102  * - `subcell::Tags::SubcellOptions`
103  *
104  * DataBox:
105  * - Uses:
106  * - `domain::Tags::Mesh<Dim>`
107  * - `subcell::Tags::Mesh<Dim>`
108  * - `Tags::TimeStepId`
109  * - `subcell::Tags::ActiveGrid`
110  * - `subcell::Tags::NeighborDataForReconstructionAndRdmpTci<Dim>`
111  * - `subcell::Tags::TciGridHistory`
112  * - Adds: nothing
113  * - Removes: nothing
114  * - Modifies:
115  * - `subcell::Tags::Inactive<System::variables_tag>`
116  * - `System::variables_tag` if the cell is not troubled
117  * - `Tags::HistoryEvolvedVariables` if the cell is not troubled
118  * - `subcell::Tags::ActiveGrid` if the cell is not troubled
119  * - `subcell::Tags::DidRollback` sets to `false`
120  * - `subcell::Tags::NeighborDataForReconstructionAndRdmpTci<Dim>`
121  * if the cell is not troubled
122  * - `subcell::Tags::TciGridHistory` if the time stepper is a multistep method
123  */
124 template <typename TciMutator>
126  template <typename DbTags, typename... InboxTags, typename Metavariables,
127  typename ArrayIndex, typename ActionList,
128  typename ParallelComponent, size_t Dim = Metavariables::volume_dim>
129  static std::tuple<db::DataBox<DbTags>&&> apply(
130  db::DataBox<DbTags>& box,
131  const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
132  const Parallel::GlobalCache<Metavariables>& /*cache*/,
133  const ArrayIndex& /*array_index*/, const ActionList /*meta*/,
134  const ParallelComponent* const /*meta*/) noexcept {
135  static_assert(
136  tmpl::count_if<
137  ActionList,
138  std::is_same<tmpl::_1, tmpl::pin<TciAndSwitchToDg>>>::value == 1,
139  "Must have the TciAndSwitchToDg action exactly once in the action list "
140  "of a phase.");
141 
142  db::mutate<subcell::Tags::DidRollback>(
143  make_not_null(&box),
144  [](const gsl::not_null<bool*> did_rollback) noexcept {
145  *did_rollback = false;
146  });
147 
148  const TimeStepId& time_step_id = db::get<::Tags::TimeStepId>(box);
149  const SubcellOptions& subcell_options = db::get<Tags::SubcellOptions>(box);
150  if (time_step_id.substep() != 0 or
151  UNLIKELY(time_step_id.slab_number() < 0) or
152  UNLIKELY(subcell_options.always_use_subcells())) {
153  // The first condition is that for substep time integrators we only allow
154  // switching back to DG on step boundaries. This is the easiest way to
155  // avoid having a shock in the time stepper history, since there is no
156  // history at step boundaries.
157  //
158  // The second condition is that if we are in the self-start procedure of
159  // the time stepper, and we don't want to switch from subcell back to DG
160  // during self-start since we integrate over the same temporal region at
161  // increasingly higher order.
162  //
163  // The third condition is that the user has requested we always do
164  // subcell, so effectively a finite difference/volume code.
165  return {std::move(box)};
166  }
167 
168  using variables_tag = typename Metavariables::system::variables_tag;
169 
170  ASSERT(db::get<Tags::ActiveGrid>(box) == ActiveGrid::Subcell,
171  "Must be using subcells when calling TciAndSwitchToDg action.");
172  const Mesh<Dim>& dg_mesh = db::get<::domain::Tags::Mesh<Dim>>(box);
173  const Mesh<Dim>& subcell_mesh = db::get<subcell::Tags::Mesh<Dim>>(box);
174 
175  db::mutate<Tags::Inactive<variables_tag>>(
176  make_not_null(&box),
177  [&dg_mesh, &subcell_mesh](const auto inactive_vars_ptr,
178  const auto& active_vars) noexcept {
179  // Note: strictly speaking, to be conservative this should reconstruct
180  // uJ instead of u.
181  fd::reconstruct(inactive_vars_ptr, active_vars, dg_mesh,
182  subcell_mesh.extents());
183  },
184  db::get<variables_tag>(box));
185 
186  // Run RDMP TCI since no user info beyond the input file options are needed
187  // for that.
188  const std::pair self_id{Direction<Dim>::lower_xi(),
190  ASSERT(
192  self_id) != 0,
193  "The self ID is not in the NeighborData but should have been added "
194  "before TciAndSwitchToDg was called.");
195  const NeighborData& self_neighbor_data =
196  db::get<Tags::NeighborDataForReconstructionAndRdmpTci<Dim>>(box).at(
197  self_id);
198 
199  // Note: we assume the max/min over all neighbors and ourselves at the past
200  // time step has been collected into
201  // `self_neighbor_data.max/min_variables_values`
202  bool cell_is_troubled =
203  rdmp_tci(db::get<variables_tag>(box),
205  self_neighbor_data.max_variables_values,
206  self_neighbor_data.min_variables_values,
207  subcell_options.rdmp_delta0(), subcell_options.rdmp_epsilon());
208 
209  // If the RDMP TCI marked the candidate as acceptable, check with the
210  // user-specified TCI, since that could be stricter. We pass in the Persson
211  // exponent with one added in order to avoid flip-flopping between DG and
212  // subcell. That is, `persson_exponent+1.0` is stricter than
213  // `persson_exponent`.
214  if (not cell_is_troubled) {
215  cell_is_troubled = db::mutate_apply<TciMutator>(
216  make_not_null(&box), subcell_options.persson_exponent() + 1.0);
217  }
218 
219  // If the cell is not troubled, then we _might_ be able to switch back to
220  // DG. This depends on the type of time stepper we are using:
221  // - ADER: Not yet implemented, but here the TCI history is irrelevant
222  // because it is a one-step scheme, so we can always switch.
223  // - Multistep: the TCI must have marked the entire evolved variable history
224  // as free of shocks. In practice for LMMs this means the TCI
225  // history is as long as the evolved variables history and that
226  // the entire TCI history is `ActiveGrid::Dg`.
227  // - Substep: the easiest is to restrict switching back to DG to step
228  // boundaries where there is no history.
229  const auto& time_stepper = db::get<::Tags::TimeStepper<>>(box);
230  const bool is_substep_method = time_stepper.number_of_substeps() != 1;
231  ASSERT(time_stepper.number_of_substeps() != 0,
232  "Don't know how to handle a time stepper with zero substeps. This "
233  "might be totally fine, but the case should be thought about.");
234  if (const auto& tci_history = db::get<subcell::Tags::TciGridHistory>(box);
235  not cell_is_troubled and
236  (is_substep_method or
237  (tci_history.size() == time_stepper.order() and
238  alg::all_of(tci_history, [](const ActiveGrid tci_decision) noexcept {
239  return tci_decision == ActiveGrid::Dg;
240  })))) {
241  db::mutate<variables_tag, Tags::Inactive<variables_tag>,
246  make_not_null(&box),
247  [&dg_mesh, &subcell_mesh](
248  const auto active_vars_ptr, const auto inactive_vars_ptr,
249  const auto active_history_ptr,
250  const gsl::not_null<ActiveGrid*> active_grid_ptr,
251  const auto subcell_neighbor_data_ptr,
252  const gsl::not_null<
254  tci_grid_history_ptr) noexcept {
255  using std::swap;
256  swap(*active_vars_ptr, *inactive_vars_ptr);
257 
258  // Reconstruct the DG solution for each time in the time stepper
259  // history
260  using dt_variables_tag =
262  TimeSteppers::History<typename variables_tag::type,
263  typename dt_variables_tag::type>
264  dg_history{active_history_ptr->integration_order()};
265  const auto end_it = active_history_ptr->end();
266  for (auto it = active_history_ptr->begin(); it != end_it; ++it) {
267  dg_history.insert(it.time_step_id(),
268  fd::reconstruct(it.derivative(), dg_mesh,
269  subcell_mesh.extents()));
270  }
271  dg_history.most_recent_value() =
272  fd::reconstruct(active_history_ptr->most_recent_value(),
273  dg_mesh, subcell_mesh.extents()),
274  *active_history_ptr = std::move(dg_history);
275  *active_grid_ptr = ActiveGrid::Dg;
276 
277  // Clear the neighbor data needed for subcell reconstruction since
278  // we have now completed the time step.
279  subcell_neighbor_data_ptr->clear();
280 
281  // Clear the TCI grid history since we don't need to use it when on
282  // the DG grid.
283  tci_grid_history_ptr->clear();
284  });
285  return {std::move(box)};
286  }
287 
288  if (not is_substep_method) {
289  // For multistep methods we need to record the TCI decision history.
290  // We track the TCI decision, not which grid we are on because for
291  // multistep methods we need the discontinuity to clear the entire
292  // history before we can switch back to DG.
293  db::mutate<evolution::dg::subcell::Tags::TciGridHistory>(
294  make_not_null(&box),
295  [cell_is_troubled,
296  &time_stepper](const gsl::not_null<
298  tci_grid_history) noexcept {
299  tci_grid_history->push_front(cell_is_troubled ? ActiveGrid::Subcell
300  : ActiveGrid::Dg);
301  if (tci_grid_history->size() > time_stepper.order()) {
302  tci_grid_history->pop_back();
303  }
304  });
305  }
306  return {std::move(box)};
307  }
308 };
309 } // namespace evolution::dg::subcell::Actions
gsl::at
constexpr T & at(std::array< T, N > &arr, Size index)
Retrieve a entry from a container, with checks in Debug mode that the index being retrieved is valid.
Definition: Gsl.hpp:125
std::is_same
ElementId::external_boundary_id
static ElementId< VolumeDim > external_boundary_id() noexcept
Returns an ElementId meant for identifying data on external boundaries, which should never correspond...
evolution::dg::subcell::Actions
Actions for the DG-subcell hybrid solver.
Definition: Actions.hpp:9
utility
UNLIKELY
#define UNLIKELY(x)
Definition: Gsl.hpp:73
Parallel::GlobalCache
Definition: ElementReceiveInterpPoints.hpp:15
std::pair
Tags.hpp
db::add_tag_prefix
typename detail::add_tag_prefix_impl< Prefix, Tag, Args... >::type add_tag_prefix
Definition: PrefixHelpers.hpp:51
tuple
db::get
const auto & get(const DataBox< TagList > &box) noexcept
Retrieve the item with tag Tag from the DataBox.
Definition: DataBox.hpp:791
Direction::lower_xi
static Direction< VolumeDim > lower_xi() noexcept
Helper functions for creating specific Directions. These are labeled by the logical-coordinate names ...
evolution::dg::subcell::ActiveGrid
ActiveGrid
Definition: ActiveGrid.hpp:11
ElementId.hpp
evolution::dg::subcell::Tags::NeighborDataForReconstructionAndRdmpTci
The neighbor data for reconstruction and the RDMP troubled-cell indicator.
Definition: NeighborData.hpp:25
DataBox.hpp
cstddef
Assert.hpp
evolution::dg::subcell::Tags::TciGridHistory
A record of which grid the TCI requested we use.
Definition: TciGridHistory.hpp:20
deque
tuples::TaggedTuple
An associative container that is indexed by structs.
Definition: TaggedTuple.hpp:271
evolution::dg::subcell::SubcellOptions
Holds the system-agnostic subcell parameters, such as numbers controlling when to switch between DG a...
Definition: SubcellOptions.hpp:23
TimeStepId
Definition: TimeStepId.hpp:25
ASSERT
#define ASSERT(a, m)
Assert that an expression should be true.
Definition: Assert.hpp:49
Variables.hpp
Mesh
Holds the number of grid points, basis, and quadrature in each direction of the computational grid.
Definition: Mesh.hpp:49
TimeSteppers::History::integration_order
size_t integration_order() const noexcept
Get or set the current order of integration. TimeSteppers may impose restrictions on the valid values...
Definition: History.hpp:105
evolution::dg::subcell::NeighborData
Holds neighbor data needed for the TCI and reconstruction.
Definition: NeighborData.hpp:24
TimeSteppers::History
Definition: History.hpp:33
TimeStepId.hpp
Gsl.hpp
evolution::dg::subcell::rdmp_tci
bool rdmp_tci(const Variables< tmpl::list< EvolvedVarsTags... >> &active_grid_candidate_evolved_vars, const Variables< tmpl::list< Tags::Inactive< EvolvedVarsTags >... >> &inactive_grid_candidate_evolved_vars, const std::vector< double > &max_of_past_variables, const std::vector< double > &min_of_past_variables, const double rdmp_delta0, const double rdmp_epsilon) noexcept
Troubled cell indicator using a relaxed discrete maximum principle, comparing the candidate solution ...
Definition: RdmpTci.hpp:61
evolution::dg::subcell::Actions::TciAndSwitchToDg
Run the troubled-cell indicator on the subcell solution to see if it is safe to switch back to DG.
Definition: TciAndSwitchToDg.hpp:125
make_not_null
gsl::not_null< T * > make_not_null(T *ptr) noexcept
Construct a not_null from a pointer. Often this will be done as an implicit conversion,...
Definition: Gsl.hpp:880
Direction.hpp
Tags::HistoryEvolvedVariables
Definition: Tags.hpp:85
evolution::dg::subcell::Tags::Inactive
Mark a tag as holding data for the inactive grid.
Definition: Inactive.hpp:23
Prefixes.hpp
alg::all_of
decltype(auto) all_of(const Container &c, UnaryPredicate &&unary_predicate)
Convenience wrapper around std::all_of.
Definition: Algorithm.hpp:181
Parallel
Functionality for parallelization.
Definition: ElementReceiveInterpPoints.hpp:13
type_traits
TMPL.hpp
evolution::dg::subcell::fd::reconstruct
DataVector reconstruct(const DataVector &subcell_u_times_projected_det_jac, const Mesh< Dim > &dg_mesh, const Index< Dim > &subcell_extents) noexcept
reconstruct the variable subcell_u_times_projected_det_jac onto the DG grid dg_mesh.
Mesh.hpp
gsl::not_null
Require a pointer to not be a nullptr
Definition: ReadSpecPiecewisePolynomial.hpp:13