TciAndRollback.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 <iterator>
9 #include <tuple>
10 #include <type_traits>
11 #include <utility>
12 
14 #include "DataStructures/DataBox/PrefixHelpers.hpp"
17 #include "Domain/Tags.hpp"
18 #include "Evolution/DgSubcell/Actions/Labels.hpp"
19 #include "Evolution/DgSubcell/ActiveGrid.hpp"
20 #include "Evolution/DgSubcell/NeighborData.hpp"
21 #include "Evolution/DgSubcell/Projection.hpp"
22 #include "Evolution/DgSubcell/RdmpTci.hpp"
23 #include "Evolution/DgSubcell/SubcellOptions.hpp"
24 #include "Evolution/DgSubcell/Tags/ActiveGrid.hpp"
25 #include "Evolution/DgSubcell/Tags/Coordinates.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"
32 #include "Evolution/DiscontinuousGalerkin/InboxTags.hpp"
34 #include "Parallel/Actions/Goto.hpp"
35 #include "Parallel/GlobalCache.hpp"
36 #include "Time/Actions/SelfStartActions.hpp"
37 #include "Time/History.hpp"
38 #include "Time/Tags.hpp"
40 #include "Utilities/TMPL.hpp"
41 #include "Utilities/TaggedTuple.hpp"
42 
44 /*!
45  * \brief Run the troubled-cell indicator on the candidate solution and perform
46  * the time step rollback if needed.
47  *
48  * The troubled-cell indicator (TCI) is given by the mutator `TciMutator` and
49  * can mutate tags in the DataBox, but should do so cautiously. The main reason
50  * that this is a mutator is because primitive variables, such as the pressure,
51  * are used to check if the solution is physical. In the relativistic case, even
52  * just whether or not the primitive variables can be recovered is used as a
53  * condition. Note that the evolved variables are projected to the subcells
54  * _before_ the TCI is called, and so `Tags::Inactive<variables_tag>` can be
55  * used to retrieve the candidate solution on the subcell (e.g. for an RDMP
56  * TCI).
57  *
58  * After rollback, the subcell scheme must project the DG boundary corrections
59  * \f$G\f$ to the subcells for the scheme to be conservative. The subcell
60  * actions know if a rollback was done because the local mortar data would
61  * already be computed.
62  */
63 template <typename TciMutator>
65  template <typename DbTags, typename... InboxTags, typename Metavariables,
66  typename ArrayIndex, typename ActionList,
67  typename ParallelComponent, size_t Dim = Metavariables::volume_dim>
68  static std::tuple<db::DataBox<DbTags>&&, bool, size_t> apply(
69  db::DataBox<DbTags>& box,
70  const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
71  const Parallel::GlobalCache<Metavariables>& /*cache*/,
72  const ArrayIndex& /*array_index*/, const ActionList /*meta*/,
73  const ParallelComponent* const /*meta*/) noexcept {
74  static_assert(
75  tmpl::count_if<
76  ActionList,
77  std::is_same<tmpl::_1, tmpl::pin<TciAndRollback>>>::value == 1,
78  "Must have the TciAndRollback action exactly once in the action list "
79  "of a phase.");
80  static_assert(
81  tmpl::count_if<
82  ActionList,
83  std::is_same<tmpl::_1,
84  tmpl::pin<::Actions::Label<
86  BeginSubcellAfterDgRollback>>>>::value == 1,
87  "Must have the BeginSubcellAfterDgRollback label exactly once in the "
88  "action list of a phase.");
89 
90  using variables_tag = typename Metavariables::system::variables_tag;
91  using dt_variables_tag = db::add_tag_prefix<::Tags::dt, variables_tag>;
92 
93  const ActiveGrid active_grid = db::get<Tags::ActiveGrid>(box);
94  ASSERT(active_grid == ActiveGrid::Dg,
95  "Must be using DG when calling TciAndRollback action.");
96 
97  const Mesh<Dim>& dg_mesh = db::get<::domain::Tags::Mesh<Dim>>(box);
98  const Mesh<Dim>& subcell_mesh = db::get<Tags::Mesh<Dim>>(box);
99 
100  // Project candidate solution to subcells
101  db::mutate<Tags::Inactive<variables_tag>>(
102  make_not_null(&box),
103  [&dg_mesh, &subcell_mesh](const auto inactive_vars_ptr,
104  const auto& active_vars) noexcept {
105  // Note: strictly speaking, to be conservative this should project uJ
106  // instead of u.
107  fd::project(inactive_vars_ptr, active_vars, dg_mesh,
108  subcell_mesh.extents());
109  },
110  db::get<variables_tag>(box));
111 
112  // Run RDMP TCI since no user info beyond the input file options are needed
113  // for that.
114  std::pair self_id{
117  ASSERT(
119  self_id) != 0,
120  "The self ID is not in the NeighborData.");
121  const NeighborData& self_neighbor_data =
122  db::get<Tags::NeighborDataForReconstructionAndRdmpTci<Dim>>(box).at(
123  self_id);
124  const SubcellOptions& subcell_options = db::get<Tags::SubcellOptions>(box);
125  // Note: we assume the max/min over all neighbors and ourselves at the past
126  // time step has been collected into
127  // `self_neighbor_data.max/min_variables_values`
128  bool cell_is_troubled =
129  rdmp_tci(db::get<variables_tag>(box),
131  self_neighbor_data.max_variables_values,
132  self_neighbor_data.min_variables_values,
133  subcell_options.rdmp_delta0(), subcell_options.rdmp_epsilon());
134 
135  // If the RDMP TCI marked the candidate as acceptable, check with the
136  // user-specified TCI, since that could be stricter.
137  if (not cell_is_troubled) {
138  // The reason we pass in the persson_exponent explicitly instead of
139  // leaving it to the user is because the value of the exponent that should
140  // be used to decide if it is safe to switch back to DG should be
141  // `persson_exponent+1` to prevent the code from rapidly switching back
142  // and forth between DG and subcell. Rather than trying to enforce this by
143  // documentation, the switching back to DG TCI gets passed in the exponent
144  // it should use, and to keep the interface between the TCIs consistent,
145  // we also pass the exponent in separately here.
146  cell_is_troubled = db::mutate_apply<TciMutator>(
147  make_not_null(&box), subcell_options.persson_exponent());
148  }
149  if (cell_is_troubled or
150  (subcell_options.always_use_subcells() and
152  .external_boundaries()
153  .empty())) {
154  db::mutate<variables_tag, Tags::Inactive<variables_tag>,
157  make_not_null(&box),
158  [&dg_mesh, &subcell_mesh](
159  const auto active_vars_ptr, const auto inactive_vars_ptr,
160  const auto active_history_ptr,
161  const gsl::not_null<ActiveGrid*> active_grid_ptr,
162  const gsl::not_null<bool*> did_rollback_ptr) noexcept {
163  ASSERT(
164  active_history_ptr->size() > 0,
165  "We cannot have an empty history when unwinding, that's just "
166  "nutty. Did you call the action too early in the action list?");
167  // Rollback u^{n+1}* to u^n (undoing the candidate solution) by
168  // using the time stepper history.
169  *active_vars_ptr = std::prev(active_history_ptr->end()).value();
170  fd::project(inactive_vars_ptr, *active_vars_ptr, dg_mesh,
171  subcell_mesh.extents());
172  using std::swap;
173  swap(*active_vars_ptr, *inactive_vars_ptr);
174 
175  // Project the time stepper history to the subcells, excluding the
176  // most recent inadmissible history.
177  TimeSteppers::History<typename variables_tag::type,
178  typename dt_variables_tag::type>
179  subcell_history{active_history_ptr->integration_order()};
180  const auto end_it = std::prev(active_history_ptr->end());
181  for (auto it = active_history_ptr->begin(); it != end_it; ++it) {
182  subcell_history.insert(
183  it.time_step_id(),
184  fd::project(it.value(), dg_mesh, subcell_mesh.extents()),
185  fd::project(it.derivative(), dg_mesh,
186  subcell_mesh.extents()));
187  }
188  *active_history_ptr = std::move(subcell_history);
189  *active_grid_ptr = ActiveGrid::Subcell;
190  *did_rollback_ptr = true;
191  // Note: We do _not_ project the boundary history here because that
192  // needs to be done at the lifting stage of the subcell method,
193  // since we need to lift G+D instead of the ingredients that go into
194  // G+D, which is what we would be projecting here.
195  });
196 
197  if (UNLIKELY(db::get<::Tags::TimeStepId>(box).slab_number() < 0)) {
198  // If we are doing self start, then we need to project the initial guess
199  // to the subcells as well.
200  //
201  // Warning: this unfortunately needs to be kept in sync with the
202  // self-start procedure.
203  //
204  // Note: if we switch to the subcells then we might have an inconsistent
205  // state between the primitive and conservative variables on the
206  // subcells. The most correct thing is to re-compute the primitive
207  // variables on the subcells, since projecting the conservative
208  // variables is conservative.
209  if constexpr (Metavariables::system::
210  has_primitive_and_conservative_vars) {
211  db::mutate<
214  typename Metavariables::system::primitive_variables_tag>>(
215  make_not_null(&box),
216  [&dg_mesh, &subcell_mesh](
217  const auto initial_vars_ptr,
218  const auto initial_prim_vars_ptr) noexcept {
219  // Note: for strict conservation, we need to project uJ instead
220  // of just u.
221  std::get<0>(*initial_vars_ptr) =
222  fd::project(std::get<0>(*initial_vars_ptr), dg_mesh,
223  subcell_mesh.extents());
224  std::get<0>(*initial_prim_vars_ptr) =
225  fd::project(std::get<0>(*initial_prim_vars_ptr), dg_mesh,
226  subcell_mesh.extents());
227  });
228  } else {
229  db::mutate<SelfStart::Tags::InitialValue<variables_tag>>(
230  make_not_null(&box),
231  [&dg_mesh, &subcell_mesh](const auto initial_vars_ptr) noexcept {
232  // Note: for strict conservation, we need to project uJ instead
233  // of just u.
234  std::get<0>(*initial_vars_ptr) =
235  fd::project(std::get<0>(*initial_vars_ptr), dg_mesh,
236  subcell_mesh.extents());
237  });
238  }
239  }
240 
241  return {std::move(box), false,
242  tmpl::index_of<
243  ActionList,
244  ::Actions::Label<evolution::dg::subcell::Actions::Labels::
245  BeginSubcellAfterDgRollback>>::value +
246  1};
247  }
248  // The unlimited DG solver has passed, so we can remove the current neighbor
249  // data.
250  db::mutate<subcell::Tags::NeighborDataForReconstructionAndRdmpTci<Dim>>(
251  make_not_null(&box), [](const auto neighbor_data_ptr) noexcept {
252  neighbor_data_ptr->clear();
253  });
254 
255  return {std::move(box), false,
256  tmpl::index_of<ActionList, TciAndRollback>::value + 1};
257  }
258 };
259 } // 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
evolution::dg::subcell::Actions::Labels
Labels used to navigate the action list when using a DG-subcell scheme.
Definition: Labels.hpp:8
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
evolution::dg::subcell::fd::project
DataVector project(const DataVector &dg_u, const Mesh< Dim > &dg_mesh, const Index< Dim > &subcell_extents) noexcept
Project the variable dg_u onto the subcell grid with extents subcell_extents.
UNLIKELY
#define UNLIKELY(x)
Definition: Gsl.hpp:73
Parallel::GlobalCache
Definition: ElementReceiveInterpPoints.hpp:15
std::pair
GlobalCache.hpp
Tags.hpp
domain::Tags::Element
Definition: Tags.hpp:97
iterator
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 ...
SelfStart::Tags::InitialValue
Definition: SelfStartActions.hpp:105
evolution::dg::subcell::ActiveGrid
ActiveGrid
Definition: ActiveGrid.hpp:11
db::mutate
decltype(auto) mutate(const gsl::not_null< DataBox< TagList > * > box, Invokable &&invokable, Args &&... args) noexcept
Allows changing the state of one or more non-computed elements in the DataBox.
Definition: DataBox.hpp:641
evolution::dg::subcell::Actions::TciAndRollback
Run the troubled-cell indicator on the candidate solution and perform the time step rollback if neede...
Definition: TciAndRollback.hpp:64
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
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
ASSERT
#define ASSERT(a, m)
Assert that an expression should be true.
Definition: Assert.hpp:49
Element.hpp
Mesh
Holds the number of grid points, basis, and quadrature in each direction of the computational grid.
Definition: Mesh.hpp:48
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:100
evolution::dg::subcell::NeighborData
Holds neighbor data needed for the TCI and reconstruction.
Definition: NeighborData.hpp:24
TimeSteppers::History
Definition: History.hpp:33
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
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
Actions::Label
Definition: Goto.hpp:38
evolution::dg::subcell::Tags::DidRollback
Tag indicating whether we are retrying a step after a rollback of a failed DG step.
Definition: DidRollback.hpp:16
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
type_traits
TMPL.hpp
Mesh.hpp
gsl::not_null
Require a pointer to not be a nullptr
Definition: ReadSpecThirdOrderPiecewisePolynomial.hpp:13