Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : #pragma once
5 :
6 : #include <algorithm>
7 : #include <cstddef>
8 :
9 : #include "DataStructures/DataVector.hpp"
10 : #include "DataStructures/Tensor/Tensor.hpp"
11 : #include "DataStructures/Variables.hpp"
12 : #include "Evolution/DgSubcell/Tags/Inactive.hpp"
13 : #include "Utilities/ErrorHandling/Assert.hpp"
14 : #include "Utilities/TMPL.hpp"
15 :
16 : namespace evolution::dg::subcell {
17 : /*!
18 : * \brief Troubled cell indicator using a relaxed discrete maximum principle,
19 : * comparing the candidate solution with the past solution in the element and
20 : * its neighbors.
21 : *
22 : * Let the candidate solution be denoted by \f$u^\star_{\alpha}(t^{n+1})\f$.
23 : * Then the RDMP requires that
24 : *
25 : * \f{align*}{
26 : * \min_{\forall\mathcal{N}}\left(u_{\alpha}(t^n)\right)
27 : * - \delta_\alpha
28 : * \le
29 : * u^\star_{\alpha}(t^{n+1})
30 : * \le
31 : * \max_{\forall\mathcal{N}}
32 : * \left(u_{\alpha}(t^n)\right) + \delta_\alpha
33 : * \f}
34 : *
35 : * where \f$\mathcal{N}\f$ are either the Neumann or Voronoi neighbors and the
36 : * element itself, and \f$\delta_\alpha\f$ is a parameter defined below that
37 : * relaxes the discrete maximum principle (DMP). When computing
38 : * \f$\max(u_\alpha)\f$ and \f$\min(u_\alpha)\f$ over a DG element that is not
39 : * using subcells we first project the DG solution to the subcells and then
40 : * compute the maximum and minimum over *both* the DG grid and the subcell grid.
41 : * However, when a DG element is using subcells we compute the maximum and
42 : * minimum of \f$u_\alpha(t^n)\f$ over the subcells only. Note that the maximum
43 : * and minimum values of \f$u^\star_\alpha\f$ are always computed over both the
44 : * DG and the subcell grids, even when using the RDMP to check if the
45 : * reconstructed DG solution would be admissible.
46 : *
47 : * The parameter \f$\delta_\alpha\f$ is given by:
48 : *
49 : * \f{align*}{
50 : * \delta_\alpha =
51 : * \max\left(\delta_{0},\epsilon
52 : * \left(\max_{\forall\mathcal{N}}\left(u_{\alpha}(t^n)\right)
53 : * - \min_{\forall\mathcal{N}}\left(u_{\alpha}(t^n)\right)\right)
54 : * \right),
55 : * \f}
56 : *
57 : * where we typically take \f$\delta_{0}=10^{-4}\f$ and \f$\epsilon=10^{-3}\f$.
58 : *
59 : * If all checks are passed and cell is not troubled, returns an integer `0`.
60 : * Otherwise returns 1-based index of the tag in the input Variables that fails
61 : * the check. For instance, if we have following two Variables objects as
62 : * candidate solutions on active and inactive grids
63 : *
64 : * - `Variables<tmpl::list<DgVar1, DgVar2, DgVar3>>`
65 : * - `Variables<tmpl::list<SubVar1, SubVar2, SubVar3>>`
66 : *
67 : * and TCI flags the second pair `DgVar2` and `SubVar2` not satisfying two-mesh
68 : * RDMP criteria, returned value is `2` since the second pair of tags failed the
69 : * check.
70 : *
71 : * \note Once a single pair of tags fails to satisfy the check, checks for the
72 : * remaining part of the input variables are skipped. In the example above, for
73 : * instance if the second pair (`DgVar2`,`SubVar2`) is flagged as troubled, the
74 : * third pair (`DgVar3`,`SubVar3`) is ignored and not checked.
75 : *
76 : */
77 : template <typename... EvolvedVarsTags>
78 1 : int rdmp_tci(const Variables<tmpl::list<EvolvedVarsTags...>>&
79 : active_grid_candidate_evolved_vars,
80 : const Variables<tmpl::list<Tags::Inactive<EvolvedVarsTags>...>>&
81 : inactive_grid_candidate_evolved_vars,
82 : const DataVector& max_of_past_variables,
83 : const DataVector& min_of_past_variables, const double rdmp_delta0,
84 : const double rdmp_epsilon) {
85 : bool cell_is_troubled = false;
86 : int rdmp_tci_status = 0;
87 : size_t component_index = 0;
88 :
89 : tmpl::for_each<tmpl::list<EvolvedVarsTags...>>(
90 : [&active_grid_candidate_evolved_vars, &cell_is_troubled, &component_index,
91 : &inactive_grid_candidate_evolved_vars, &max_of_past_variables,
92 : &min_of_past_variables, &rdmp_tci_status, rdmp_delta0,
93 : rdmp_epsilon](auto tag_v) {
94 : if (cell_is_troubled) {
95 : return;
96 : }
97 : using std::max;
98 : using std::min;
99 :
100 : using tag = tmpl::type_from<decltype(tag_v)>;
101 : using inactive_tag = Tags::Inactive<tag>;
102 : const auto& active_var = get<tag>(active_grid_candidate_evolved_vars);
103 : const auto& inactive_var =
104 : get<inactive_tag>(inactive_grid_candidate_evolved_vars);
105 : const size_t number_of_components = active_var.size();
106 : ASSERT(number_of_components == inactive_var.size(),
107 : "The active and inactive vars must have the same type of tensor "
108 : "and therefore the same number of components.");
109 :
110 : for (size_t tensor_storage_index = 0;
111 : tensor_storage_index < number_of_components;
112 : ++tensor_storage_index) {
113 : ASSERT(not cell_is_troubled,
114 : "If a cell has already been marked as troubled during the "
115 : "RDMP TCI, we should not be continuing to check other "
116 : "variables.");
117 : const double max_active = max(active_var[tensor_storage_index]);
118 : const double min_active = min(active_var[tensor_storage_index]);
119 : const double max_inactive = max(inactive_var[tensor_storage_index]);
120 : const double min_inactive = min(inactive_var[tensor_storage_index]);
121 : const double delta =
122 : max(rdmp_delta0,
123 : rdmp_epsilon * (max_of_past_variables[component_index] -
124 : min_of_past_variables[component_index]));
125 : cell_is_troubled =
126 : max(max_active, max_inactive) >
127 : max_of_past_variables[component_index] + delta or
128 : min(min_active, min_inactive) <
129 : min_of_past_variables[component_index] - delta;
130 : if (cell_is_troubled) {
131 : rdmp_tci_status = static_cast<int>(component_index + 1);
132 : return;
133 : }
134 : ++component_index;
135 : }
136 : });
137 : return rdmp_tci_status;
138 : }
139 :
140 : /*!
141 : * \brief get the max and min of each component of the active and inactive
142 : * variables. If `include_inactive_grid` is `false` then only the max over the
143 : * `active_grid_evolved_vars` for each component is returned.
144 : */
145 : template <typename... EvolvedVarsTags>
146 1 : std::pair<DataVector, DataVector> rdmp_max_min(
147 : const Variables<tmpl::list<EvolvedVarsTags...>>& active_grid_evolved_vars,
148 : const Variables<tmpl::list<Tags::Inactive<EvolvedVarsTags>...>>&
149 : inactive_grid_evolved_vars,
150 : const bool include_inactive_grid) {
151 : DataVector max_of_vars{
152 : active_grid_evolved_vars.number_of_independent_components,
153 : std::numeric_limits<double>::min()};
154 : DataVector min_of_vars{
155 : active_grid_evolved_vars.number_of_independent_components,
156 : std::numeric_limits<double>::max()};
157 : size_t component_index = 0;
158 : tmpl::for_each<tmpl::list<EvolvedVarsTags...>>(
159 : [&active_grid_evolved_vars, &component_index, &inactive_grid_evolved_vars,
160 : &include_inactive_grid, &max_of_vars, &min_of_vars](auto tag_v) {
161 : using std::max;
162 : using std::min;
163 :
164 : using tag = tmpl::type_from<decltype(tag_v)>;
165 : const auto& active_var = get<tag>(active_grid_evolved_vars);
166 : const size_t number_of_components_in_tensor = active_var.size();
167 : for (size_t tensor_storage_index = 0;
168 : tensor_storage_index < number_of_components_in_tensor;
169 : ++tensor_storage_index) {
170 : ASSERT(component_index < max_of_vars.size() and
171 : component_index < min_of_vars.size(),
172 : "The component index into the variables is out of bounds.");
173 : max_of_vars[component_index] = max(active_var[tensor_storage_index]);
174 : min_of_vars[component_index] = min(active_var[tensor_storage_index]);
175 : if (include_inactive_grid) {
176 : using inactive_tag = Tags::Inactive<tag>;
177 : const auto& inactive_var =
178 : get<inactive_tag>(inactive_grid_evolved_vars);
179 : max_of_vars[component_index] =
180 : max(max_of_vars[component_index],
181 : max(inactive_var[tensor_storage_index]));
182 : min_of_vars[component_index] =
183 : min(min_of_vars[component_index],
184 : min(inactive_var[tensor_storage_index]));
185 : }
186 : ++component_index;
187 : }
188 : });
189 : return {std::move(max_of_vars), std::move(min_of_vars)};
190 : }
191 :
192 : /*!
193 : * \brief Check if the current variables satisfy the RDMP. Returns an integer
194 : * `0` if cell is not troubled and an integer `i+1` if the `[i]`-th element of
195 : * the input vector is responsible for failing the RDMP.
196 : *
197 : * Let the candidate solution be denoted by \f$u^\star_{\alpha}(t^{n+1})\f$.
198 : * Then the RDMP requires that
199 : *
200 : * \f{align*}{
201 : * \min_{\forall\mathcal{N}}\left(u_{\alpha}(t^n)\right)
202 : * - \delta_\alpha
203 : * \le
204 : * u^\star_{\alpha}(t^{n+1})
205 : * \le
206 : * \max_{\forall\mathcal{N}}
207 : * \left(u_{\alpha}(t^n)\right) + \delta_\alpha
208 : * \f}
209 : *
210 : * where \f$\mathcal{N}\f$ are either the Neumann or Voronoi neighbors and the
211 : * element itself, and \f$\delta_\alpha\f$ is a parameter defined below that
212 : * relaxes the discrete maximum principle (DMP). When computing
213 : * \f$\max(u_\alpha)\f$ and \f$\min(u_\alpha)\f$ over a DG element that is not
214 : * using subcells we first project the DG solution to the subcells and then
215 : * compute the maximum and minimum over *both* the DG grid and the subcell grid.
216 : * However, when a DG element is using subcells we compute the maximum and
217 : * minimum of \f$u_\alpha(t^n)\f$ over the subcells only. Note that the maximum
218 : * and minimum values of \f$u^\star_\alpha\f$ are always computed over both the
219 : * DG and the subcell grids, even when using the RDMP to check if the
220 : * reconstructed DG solution would be admissible.
221 : *
222 : * The parameter \f$\delta_\alpha\f$ is given by:
223 : *
224 : * \f{align*}{
225 : * \delta_\alpha =
226 : * \max\left(\delta_{0},\epsilon
227 : * \left(\max_{\forall\mathcal{N}}\left(u_{\alpha}(t^n)\right)
228 : * - \min_{\forall\mathcal{N}}\left(u_{\alpha}(t^n)\right)\right)
229 : * \right),
230 : * \f}
231 : *
232 : * where we typically take \f$\delta_{0}=10^{-4}\f$ and \f$\epsilon=10^{-3}\f$.
233 : *
234 : * If all checks are passed and cell is not troubled, returns an integer `0`.
235 : * Otherwise returns an 1-based index of the element in the input
236 : * `DataVector` that fails the check.
237 : *
238 : * e.g. Suppose we have three variables to check RDMP so that
239 : * `max_of_current_variables.size() == 3`. If RDMP TCI flags
240 : * `max_of_current_variables[1]`, `min_of_current_variables[1]`, .. (and so on)
241 : * as troubled, returned integer value is `2`.
242 : *
243 : * Once cell is marked as troubled, checks for the remaining part of the input
244 : * `std::vector`s are skipped. In the example above, for instance if `[1]`-th
245 : * component of inputs is flagged as troubled, checking the remaining index
246 : * `[2]` is skipped.
247 : *
248 : */
249 1 : int rdmp_tci(const DataVector& max_of_current_variables,
250 : const DataVector& min_of_current_variables,
251 : const DataVector& max_of_past_variables,
252 : const DataVector& min_of_past_variables, double rdmp_delta0,
253 : double rdmp_epsilon);
254 : } // namespace evolution::dg::subcell
|
