RdmpTci.hpp
1 // Distributed under the MIT License.
2 // See LICENSE.txt for details.
3 
4 #pragma once
5 
6 #include <algorithm>
7 #include <cstddef>
8 #include <vector>
9 
10 #include "DataStructures/DataVector.hpp"
13 #include "Evolution/DgSubcell/Tags/Inactive.hpp"
15 #include "Utilities/TMPL.hpp"
16 
17 namespace evolution::dg::subcell {
18 /*!
19  * \brief Troubled cell indicator using a relaxed discrete maximum principle,
20  * comparing the candidate solution with the past solution in the element and
21  * its neighbors.
22  *
23  * Let the candidate solution be denoted by \f$u^\star_{\alpha}(t^{n+1})\f$.
24  * Then the RDMP requires that
25  *
26  * \f{align*}{
27  * \min_{\forall\mathcal{N}}\left(u_{\alpha}(t^n)\right)
28  * - \delta_\alpha
29  * \le
30  * u^\star_{\alpha}(t^{n+1})
31  * \le
32  * \max_{\forall\mathcal{N}}
33  * \left(u_{\alpha}(t^n)\right) + \delta_\alpha
34  * \f}
35  *
36  * where \f$\mathcal{N}\f$ are either the Neumann or Voronoi neighbors and the
37  * element itself, and \f$\delta_\alpha\f$ is a parameter defined below that
38  * relaxes the discrete maximum principle (DMP). When computing
39  * \f$\max(u_\alpha)\f$ and \f$\min(u_\alpha)\f$ over a DG element that is not
40  * using subcells we first project the DG solution to the subcells and then
41  * compute the maximum and minimum over *both* the DG grid and the subcell grid.
42  * However, when a DG element is using subcells we compute the maximum and
43  * minimum of \f$u_\alpha(t^n)\f$ over the subcells only. Note that the maximum
44  * and minimum values of \f$u^\star_\alpha\f$ are always computed over both the
45  * DG and the subcell grids, even when using the RDMP to check if the
46  * reconstructed DG solution would be admissible.
47  *
48  * The parameter \f$\delta_\alpha\f$ is given by:
49  *
50  * \f{align*}{
51  * \delta_\alpha =
52  * \max\left(\delta_{0},\epsilon
53  * \left(\max_{\forall\mathcal{N}}\left(u_{\alpha}(t^n)\right)
54  * - \min_{\forall\mathcal{N}}\left(u_{\alpha}(t^n)\right)\right)
55  * \right),
56  * \f}
57  *
58  * where we typically take \f$\delta_{0}=10^{-4}\f$ and \f$\epsilon=10^{-3}\f$.
59  */
60 template <typename... EvolvedVarsTags>
61 bool rdmp_tci(const Variables<tmpl::list<EvolvedVarsTags...>>&
62  active_grid_candidate_evolved_vars,
63  const Variables<tmpl::list<Tags::Inactive<EvolvedVarsTags>...>>&
64  inactive_grid_candidate_evolved_vars,
65  const std::vector<double>& max_of_past_variables,
66  const std::vector<double>& min_of_past_variables,
67  const double rdmp_delta0, const double rdmp_epsilon) noexcept {
68  bool cell_is_troubled = false;
69  size_t component_index = 0;
70  tmpl::for_each<tmpl::list<EvolvedVarsTags...>>(
71  [&active_grid_candidate_evolved_vars, &cell_is_troubled, &component_index,
72  &inactive_grid_candidate_evolved_vars, &max_of_past_variables,
73  &min_of_past_variables, rdmp_delta0, rdmp_epsilon](auto tag_v) noexcept {
74  if (cell_is_troubled) {
75  return;
76  }
77  using std::max;
78  using std::min;
79 
80  using tag = tmpl::type_from<decltype(tag_v)>;
81  using inactive_tag = Tags::Inactive<tag>;
82  const auto& active_var = get<tag>(active_grid_candidate_evolved_vars);
83  const auto& inactive_var =
84  get<inactive_tag>(inactive_grid_candidate_evolved_vars);
85  const size_t number_of_components = active_var.size();
86  ASSERT(number_of_components == inactive_var.size(),
87  "The active and inactive vars must have the same type of tensor "
88  "and therefore the same number of components.");
89 
90  for (size_t tensor_storage_index = 0;
91  tensor_storage_index < number_of_components;
92  ++tensor_storage_index) {
93  ASSERT(not cell_is_troubled,
94  "If a cell has already been marked as troubled during the "
95  "RDMP TCI, we should not be continuing to check other "
96  "variables.");
97  const double max_active = max(active_var[tensor_storage_index]);
98  const double min_active = min(active_var[tensor_storage_index]);
99  const double max_inactive = max(inactive_var[tensor_storage_index]);
100  const double min_inactive = min(inactive_var[tensor_storage_index]);
101  const double delta =
102  max(rdmp_delta0,
103  rdmp_epsilon * (max_of_past_variables[component_index] -
104  min_of_past_variables[component_index]));
105  cell_is_troubled =
106  max(max_active, max_inactive) >
107  max_of_past_variables[component_index] + delta or
108  min(min_active, min_inactive) <
109  min_of_past_variables[component_index] - delta;
110  if (cell_is_troubled) {
111  return;
112  }
113  ++component_index;
114  }
115  });
116  return cell_is_troubled;
117 }
118 
119 /*!
120  * \brief get the max and min of each component of the active and inactive
121  * variables. If `include_inactive_grid` is `false` then only the max over the
122  * `active_grid_evolved_vars` for each component is returned.
123  */
124 template <typename... EvolvedVarsTags>
126  const Variables<tmpl::list<EvolvedVarsTags...>>& active_grid_evolved_vars,
127  const Variables<tmpl::list<Tags::Inactive<EvolvedVarsTags>...>>&
128  inactive_grid_evolved_vars,
129  const bool include_inactive_grid) noexcept {
130  std::vector<double> max_of_vars(
131  active_grid_evolved_vars.number_of_independent_components,
133  std::vector<double> min_of_vars(
134  active_grid_evolved_vars.number_of_independent_components,
136  size_t component_index = 0;
137  tmpl::for_each<tmpl::list<EvolvedVarsTags...>>(
138  [&active_grid_evolved_vars, &component_index, &inactive_grid_evolved_vars,
139  &include_inactive_grid, &max_of_vars,
140  &min_of_vars](auto tag_v) noexcept {
141  using std::max;
142  using std::min;
143 
144  using tag = tmpl::type_from<decltype(tag_v)>;
145  const auto& active_var = get<tag>(active_grid_evolved_vars);
146  const size_t number_of_components_in_tensor = active_var.size();
147  for (size_t tensor_storage_index = 0;
148  tensor_storage_index < number_of_components_in_tensor;
149  ++tensor_storage_index) {
150  ASSERT(component_index < max_of_vars.size() and
151  component_index < min_of_vars.size(),
152  "The component index into the variables is out of bounds.");
153  max_of_vars[component_index] = max(active_var[tensor_storage_index]);
154  min_of_vars[component_index] = min(active_var[tensor_storage_index]);
155  if (include_inactive_grid) {
156  using inactive_tag = Tags::Inactive<tag>;
157  const auto& inactive_var =
158  get<inactive_tag>(inactive_grid_evolved_vars);
159  max_of_vars[component_index] =
160  max(max_of_vars[component_index],
161  max(inactive_var[tensor_storage_index]));
162  min_of_vars[component_index] =
163  min(min_of_vars[component_index],
164  min(inactive_var[tensor_storage_index]));
165  }
166  ++component_index;
167  }
168  });
169  return {std::move(max_of_vars), std::move(min_of_vars)};
170 }
171 } // namespace evolution::dg::subcell
evolution::dg::subcell::rdmp_max_min
std::pair< std::vector< double >, std::vector< double > > rdmp_max_min(const Variables< tmpl::list< EvolvedVarsTags... >> &active_grid_evolved_vars, const Variables< tmpl::list< Tags::Inactive< EvolvedVarsTags >... >> &inactive_grid_evolved_vars, const bool include_inactive_grid) noexcept
get the max and min of each component of the active and inactive variables. If include_inactive_grid ...
Definition: RdmpTci.hpp:125
std::pair
vector
algorithm
cstddef
Assert.hpp
ASSERT
#define ASSERT(a, m)
Assert that an expression should be true.
Definition: Assert.hpp:49
Variables.hpp
evolution::dg::subcell
Implementation of a generic finite volume/conservative finite difference subcell limiter.
Definition: Actions.hpp:6
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
Tensor.hpp
evolution::dg::subcell::Tags::Inactive
Mark a tag as holding data for the inactive grid.
Definition: Inactive.hpp:23
std::numeric_limits
TMPL.hpp