Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : #pragma once
5 :
6 : #include <array>
7 : #include <boost/functional/hash.hpp> // IWYU pragma: keep
8 : #include <cmath>
9 : #include <cstddef>
10 : #include <limits>
11 : #include <ostream>
12 : #include <pup.h> // IWYU pragma: keep
13 : #include <type_traits>
14 : #include <unordered_map>
15 : #include <utility> // for pair
16 :
17 : #include "DataStructures/DataVector.hpp" // IWYU pragma: keep
18 : #include "DataStructures/Index.hpp"
19 : #include "DataStructures/ModalVector.hpp" // IWYU pragma: keep
20 : #include "DataStructures/Tags.hpp"
21 : #include "DataStructures/Tensor/Metafunctions.hpp" // IWYU pragma: keep
22 : #include "DataStructures/Tensor/Tensor.hpp" // IWYU pragma: keep
23 : #include "DataStructures/Variables.hpp" // IWYU pragma: keep
24 : #include "Domain/Structure/Direction.hpp"
25 : #include "Domain/Structure/DirectionalId.hpp"
26 : #include "Domain/Structure/Element.hpp" // IWYU pragma: keep
27 : #include "Domain/Structure/ElementId.hpp" // IWYU pragma: keep
28 : #include "Domain/Structure/OrientationMap.hpp" // IWYU pragma: keep
29 : #include "Domain/Structure/OrientationMapHelpers.hpp"
30 : #include "Domain/Tags.hpp" // IWYU pragma: keep
31 : #include "NumericalAlgorithms/LinearOperators/CoefficientTransforms.hpp"
32 : #include "NumericalAlgorithms/Spectral/Mesh.hpp"
33 : #include "NumericalAlgorithms/Spectral/Spectral.hpp"
34 : #include "Options/Context.hpp"
35 : #include "Options/ParseError.hpp"
36 : #include "Options/String.hpp"
37 : #include "Utilities/Algorithm.hpp"
38 : #include "Utilities/ErrorHandling/Error.hpp"
39 : #include "Utilities/Gsl.hpp"
40 : #include "Utilities/MakeArray.hpp"
41 : #include "Utilities/Math.hpp"
42 : #include "Utilities/TMPL.hpp"
43 :
44 : // IWYU pragma: no_include <algorithm>
45 : // IWYU pragma: no_forward_declare Variables
46 :
47 : /// \cond
48 : namespace Limiters {
49 : template <size_t VolumeDim, typename TagsToLimit>
50 : class Krivodonova;
51 : } // namespace Limiters
52 : /// \endcond
53 :
54 : namespace Limiters {
55 : /*!
56 : * \ingroup LimitersGroup
57 : * \brief An implementation of the Krivodonova limiter.
58 : *
59 : * The limiter is described in \cite Krivodonova2007. The Krivodonova limiter
60 : * works by limiting the highest derivatives/modal coefficients using an
61 : * aggressive minmod approach, decreasing in derivative/modal coefficient order
62 : * until no more limiting is necessary. In 3d, the function being limited is
63 : * expanded as:
64 : *
65 : * \f{align}{
66 : * u^{l,m,n}=\sum_{i,j,k=0,0,0}^{N_i,N_j,N_k}c^{l,m,n}_{i,j,k}
67 : * P_{i}(\xi)P_{j}(\eta)P_{k}(\zeta)
68 : * \f}
69 : *
70 : * where \f$\left\{\xi, \eta, \zeta\right\}\f$ are the logical coordinates,
71 : * \f$P_{i}\f$ are the Legendre polynomials, the superscript \f$\{l,m,n\}\f$
72 : * represents the element indexed by \f$l,m,n\f$, and \f$N_i,N_j\f$ and
73 : * \f$N_k\f$ are the number of collocation points minus one in the
74 : * \f$\xi,\eta,\f$ and \f$\zeta\f$ direction, respectively. The coefficients are
75 : * limited according to:
76 : *
77 : * \f{align}{
78 : * \tilde{c}^{l,m,n}_{i,j,k}=\mathrm{minmod}
79 : * &\left(c_{i,j,k}^{l,m,n},
80 : * \alpha_i\left(c^{l+1,m,n}_{i-1,j,k}-c^{l,m,n}_{i-1,j,k}\right),
81 : * \alpha_i\left(c^{l,m,n}_{i-1,j,k}-c^{l-1,m,n}_{i-1,j,k}\right),
82 : * \right.\notag \\
83 : * &\;\;\;\;
84 : * \alpha_j\left(c^{l,m+1,n}_{i,j-1,k}-c^{l,m,n}_{i,j-1,k}\right),
85 : * \alpha_j\left(c^{l,m,n}_{i,j-1,k}-c^{l,m-1,n}_{i,j-1,k}\right),
86 : * \notag \\
87 : * &\;\;\;\;\left.
88 : * \alpha_k\left(c^{l,m,n+1}_{i,j,k-1}-c^{l,m,n}_{i,j,k-1}\right),
89 : * \alpha_k\left(c^{l,m,n}_{i,j,k-1}-c^{l,m,n-1}_{i,j,k-1}\right)
90 : * \right),
91 : * \label{eq:krivodonova 3d minmod}
92 : * \f}
93 : *
94 : * where \f$\mathrm{minmod}\f$ is the minmod function defined as
95 : *
96 : * \f{align}{
97 : * \mathrm{minmod}(a,b,c,\ldots)=
98 : * \left\{
99 : * \begin{array}{ll}
100 : * \mathrm{sgn}(a)\min(\lvert a\rvert, \lvert b\rvert,
101 : * \lvert c\rvert, \ldots) & \mathrm{if} \;
102 : * \mathrm{sgn}(a)=\mathrm{sgn}(b)=\mathrm{sgn}(c)=\mathrm{sgn}(\ldots) \\
103 : * 0 & \mathrm{otherwise}
104 : * \end{array}\right.
105 : * \f}
106 : *
107 : * Krivodonova \cite Krivodonova2007 requires \f$\alpha_i\f$ to be in the range
108 : *
109 : * \f{align*}{
110 : * \frac{1}{2(2i-1)}\le \alpha_i \le 1
111 : * \f}
112 : *
113 : * where the lower bound comes from finite differencing the coefficients between
114 : * neighbor elements when using Legendre polynomials (see \cite Krivodonova2007
115 : * for details). Note that we normalize our Legendre polynomials by \f$P_i(1) =
116 : * 1\f$; this is the normalization \cite Krivodonova2007 uses in 1D, (but not in
117 : * 2D), which is why our bounds on \f$\alpha_i\f$ match Eq. 14 of
118 : * \cite Krivodonova2007 (but not Eq. 23). We relax the lower bound:
119 : *
120 : * \f{align*}{
121 : * 0 \le \alpha_i \le 1
122 : * \f}
123 : *
124 : * to allow different basis functions (e.g. Chebyshev polynomials) and to allow
125 : * the limiter to be more dissipative if necessary. The same \f$\alpha_i\f$s are
126 : * used in all dimensions.
127 : *
128 : * \note The only place where the specific choice of 1d basis
129 : * comes in is the lower bound for the \f$\alpha_i\f$s, and so in general the
130 : * limiter can be applied to any 1d or tensor product of 1d basis functions.
131 : *
132 : * The limiting procedure must be applied from the highest derivatives to the
133 : * lowest, i.e. the highest coefficients to the lowest. Let us consider a 3d
134 : * element with \f$N+1\f$ coefficients in each dimension and denote the
135 : * coefficients as \f$c_{i,j,k}\f$. Then the limiting procedure starts at
136 : * \f$c_{N,N,N}\f$, followed by \f$c_{N,N,N-1}\f$, \f$c_{N,N-1,N}\f$, and
137 : * \f$c_{N-1,N,N}\f$. A detailed example is given below. Limiting is stopped if
138 : * all symmetric pairs of coefficients are left unchanged, i.e.
139 : * \f$c_{i,j,k}=\tilde{c}_{i,j,k}\f$. By all symmetric coefficients we mean
140 : * that, for example, \f$c_{N-i,N-j,N-k}\f$, \f$c_{N-j,N-i,N-k}\f$,
141 : * \f$c_{N-k,N-j,N-i}\f$, \f$c_{N-j,N-k,N-i}\f$, \f$c_{N-i,N-k,N-j}\f$, and
142 : * \f$c_{N-k,N-i,N-j}\f$ are not limited. As a concrete example, consider a 3d
143 : * element with 3 collocation points per dimension. Each limited coefficient is
144 : * defined as (though only computed if needed):
145 : *
146 : * \f{align*}{
147 : * \tilde{c}^{l,m,n}_{2,2,2}=\mathrm{minmod}
148 : * &\left(c_{2,2,2}^{l,m,n},
149 : * \alpha_2\left(c^{l+1,m,n}_{1,2,2}-c^{l,m,n}_{1,2,2}\right),
150 : * \alpha_2\left(c^{l,m,n}_{1,2,2}-c^{l-1,m,n}_{1,2,2}\right),
151 : * \right.\\
152 : * &\;\;\;\;
153 : * \alpha_2\left(c^{l,m+1,n}_{2,1,2}-c^{l,m,n}_{2,1,2}\right),
154 : * \alpha_2\left(c^{l,m,n}_{2,1,2}-c^{l,m-1,n}_{2,1,2}\right),\\
155 : * &\;\;\;\;\left.
156 : * \alpha_2\left(c^{l,m,n+1}_{2,2,1}-c^{l,m,n}_{2,2,1}\right),
157 : * \alpha_2\left(c^{l,m,n}_{2,2,1}-c^{l,m,n-1}_{2,2,1}\right)
158 : * \right),\\
159 : * \tilde{c}^{l,m,n}_{2,2,1}=\mathrm{minmod}
160 : * &\left(c_{2,2,1}^{l,m,n},
161 : * \alpha_2\left(c^{l+1,m,n}_{1,2,1}-c^{l,m,n}_{1,2,1}\right),
162 : * \alpha_2\left(c^{l,m,n}_{1,2,1}-c^{l-1,m,n}_{1,2,1}\right),
163 : * \right.\\
164 : * &\;\;\;\;
165 : * \alpha_2\left(c^{l,m+1,n}_{2,1,1}-c^{l,m,n}_{2,1,1}\right),
166 : * \alpha_2\left(c^{l,m,n}_{2,1,1}-c^{l,m-1,n}_{2,1,1}\right),\\
167 : * &\;\;\;\;\left.
168 : * \alpha_1\left(c^{l,m,n+1}_{2,2,0}-c^{l,m,n}_{2,2,0}\right),
169 : * \alpha_1\left(c^{l,m,n}_{2,2,0}-c^{l,m,n-1}_{2,2,0}\right)
170 : * \right),\\
171 : * \tilde{c}^{l,m,n}_{2,1,2}=\mathrm{minmod}
172 : * &\left(c_{2,1,2}^{l,m,n},
173 : * \alpha_2\left(c^{l+1,m,n}_{1,1,2}-c^{l,m,n}_{1,1,2}\right),
174 : * \alpha_2\left(c^{l,m,n}_{1,1,2}-c^{l-1,m,n}_{1,1,2}\right),
175 : * \right.\\
176 : * &\;\;\;\;
177 : * \alpha_1\left(c^{l,m+1,n}_{2,0,2}-c^{l,m,n}_{2,0,2}\right),
178 : * \alpha_1\left(c^{l,m,n}_{2,0,2}-c^{l,m-1,n}_{2,0,2}\right),\\
179 : * &\;\;\;\;\left.
180 : * \alpha_2\left(c^{l,m,n+1}_{2,1,1}-c^{l,m,n}_{2,1,1}\right),
181 : * \alpha_2\left(c^{l,m,n}_{2,1,1}-c^{l,m,n-1}_{2,1,1}\right)
182 : * \right),\\
183 : * \tilde{c}^{l,m,n}_{1,2,2}=\mathrm{minmod}
184 : * &\left(c_{1,2,2}^{l,m,n},
185 : * \alpha_1\left(c^{l+1,m,n}_{0,2,2}-c^{l,m,n}_{0,2,2}\right),
186 : * \alpha_1\left(c^{l,m,n}_{0,2,2}-c^{l-1,m,n}_{0,2,2}\right),
187 : * \right.\\
188 : * &\;\;\;\;
189 : * \alpha_2\left(c^{l,m+1,n}_{1,1,2}-c^{l,m,n}_{1,1,2}\right),
190 : * \alpha_2\left(c^{l,m,n}_{1,1,2}-c^{l,m-1,n}_{1,1,2}\right),\\
191 : * &\;\;\;\;\left.
192 : * \alpha_2\left(c^{l,m,n+1}_{1,2,1}-c^{l,m,n}_{1,2,1}\right),
193 : * \alpha_2\left(c^{l,m,n}_{1,2,1}-c^{l,m,n-1}_{1,2,1}\right)
194 : * \right),\\
195 : * \tilde{c}^{l,m,n}_{2,2,0}=\mathrm{minmod}
196 : * &\left(c_{2,2,0}^{l,m,n},
197 : * \alpha_2\left(c^{l+1,m,n}_{1,2,0}-c^{l,m,n}_{1,2,0}\right),
198 : * \alpha_2\left(c^{l,m,n}_{1,2,0}-c^{l-1,m,n}_{1,2,0}\right),
199 : * \right.\\
200 : * &\;\;\;\;\left.
201 : * \alpha_2\left(c^{l,m+1,n}_{2,1,0}-c^{l,m,n}_{2,1,0}\right),
202 : * \alpha_2\left(c^{l,m,n}_{2,1,0}-c^{l,m-1,n}_{2,1,0}\right)
203 : * \right),\\
204 : * \tilde{c}^{l,m,n}_{2,0,2}=\mathrm{minmod}
205 : * &\left(c_{2,0,2}^{l,m,n},
206 : * \alpha_2\left(c^{l+1,m,n}_{1,0,2}-c^{l,m,n}_{1,0,2}\right),
207 : * \alpha_2\left(c^{l,m,n}_{1,0,2}-c^{l-1,m,n}_{1,0,2}\right),
208 : * \right.\\
209 : * &\;\;\;\;\left.
210 : * \alpha_2\left(c^{l,m,n+1}_{2,0,1}-c^{l,m,n}_{2,0,1}\right),
211 : * \alpha_2\left(c^{l,m,n}_{2,0,1}-c^{l,m,n-1}_{2,0,1}\right)
212 : * \right),\\
213 : * \tilde{c}^{l,m,n}_{0,2,2}=\mathrm{minmod}
214 : * &\left(c_{0,2,2}^{l,m,n},
215 : * \alpha_2\left(c^{l,m+1,n}_{0,1,2}-c^{l,m,n}_{0,1,2}\right),
216 : * \alpha_2\left(c^{l,m,n}_{0,1,2}-c^{l,m-1,n}_{0,1,2}\right),
217 : * \right.\\
218 : * &\;\;\;\;\left.
219 : * \alpha_2\left(c^{l,m,n+1}_{0,2,1}-c^{l,m,n}_{0,2,1}\right),
220 : * \alpha_2\left(c^{l,m,n}_{0,2,1}-c^{l,m,n-1}_{0,2,1}\right)
221 : * \right),\\
222 : * \tilde{c}^{l,m,n}_{2,1,1}=\mathrm{minmod}
223 : * &\left(c_{2,1,1}^{l,m,n},
224 : * \alpha_2\left(c^{l+1,m,n}_{1,1,1}-c^{l,m,n}_{1,1,1}\right),
225 : * \alpha_2\left(c^{l,m,n}_{1,1,1}-c^{l-1,m,n}_{1,1,1}\right),
226 : * \right.\\
227 : * &\;\;\;\;
228 : * \alpha_1\left(c^{l,m+1,n}_{2,0,1}-c^{l,m,n}_{2,0,1}\right),
229 : * \alpha_1\left(c^{l,m,n}_{2,0,1}-c^{l,m-1,n}_{2,0,1}\right),\\
230 : * &\;\;\;\;\left.
231 : * \alpha_1\left(c^{l,m,n+1}_{2,1,0}-c^{l,m,n}_{2,1,0}\right),
232 : * \alpha_1\left(c^{l,m,n}_{2,1,0}-c^{l,m,n-1}_{2,1,0}\right)
233 : * \right),\\
234 : * \tilde{c}^{l,m,n}_{1,2,1}=\mathrm{minmod}
235 : * &\left(c_{1,2,1}^{l,m,n},
236 : * \alpha_1\left(c^{l+1,m,n}_{0,2,1}-c^{l,m,n}_{0,2,1}\right),
237 : * \alpha_1\left(c^{l,m,n}_{0,2,1}-c^{l-1,m,n}_{0,2,1}\right),
238 : * \right.\\
239 : * &\;\;\;\;
240 : * \alpha_2\left(c^{l,m+1,n}_{1,1,1}-c^{l,m,n}_{1,1,1}\right),
241 : * \alpha_2\left(c^{l,m,n}_{1,1,1}-c^{l,m-1,n}_{1,1,1}\right),\\
242 : * &\;\;\;\;\left.
243 : * \alpha_1\left(c^{l,m,n+1}_{1,2,0}-c^{l,m,n}_{1,2,0}\right),
244 : * \alpha_1\left(c^{l,m,n}_{1,2,0}-c^{l,m,n-1}_{1,2,0}\right)
245 : * \right),\\
246 : * \tilde{c}^{l,m,n}_{1,1,2}=\mathrm{minmod}
247 : * &\left(c_{1,1,2}^{l,m,n},
248 : * \alpha_1\left(c^{l+1,m,n}_{0,1,2}-c^{l,m,n}_{0,1,2}\right),
249 : * \alpha_1\left(c^{l,m,n}_{0,1,2}-c^{l-1,m,n}_{0,1,2}\right),
250 : * \right.\\
251 : * &\;\;\;\;
252 : * \alpha_1\left(c^{l,m+1,n}_{1,0,2}-c^{l,m,n}_{1,0,2}\right),
253 : * \alpha_1\left(c^{l,m,n}_{1,0,2}-c^{l,m-1,n}_{1,0,2}\right),\\
254 : * &\;\;\;\;\left.
255 : * \alpha_2\left(c^{l,m,n+1}_{1,1,1}-c^{l,m,n}_{1,1,1}\right),
256 : * \alpha_2\left(c^{l,m,n}_{1,1,1}-c^{l,m,n-1}_{1,1,1}\right)
257 : * \right),
258 : * \f}
259 : * \f{align*}{
260 : * \tilde{c}^{l,m,n}_{2,1,0}=\mathrm{minmod}
261 : * &\left(c_{2,1,0}^{l,m,n},
262 : * \alpha_2\left(c^{l+1,m,n}_{1,1,0}-c^{l,m,n}_{1,1,0}\right),
263 : * \alpha_2\left(c^{l,m,n}_{1,1,0}-c^{l-1,m,n}_{1,1,0}\right),
264 : * \right.\\
265 : * &\;\;\;\;\left.
266 : * \alpha_1\left(c^{l,m+1,n}_{2,0,0}-c^{l,m,n}_{2,0,0}\right),
267 : * \alpha_1\left(c^{l,m,n}_{2,0,0}-c^{l,m-1,n}_{2,0,0}\right)
268 : * \right),\\
269 : * \tilde{c}^{l,m,n}_{2,0,1}=\mathrm{minmod}
270 : * &\left(c_{2,0,1}^{l,m,n},
271 : * \alpha_2\left(c^{l+1,m,n}_{1,0,1}-c^{l,m,n}_{1,0,1}\right),
272 : * \alpha_2\left(c^{l,m,n}_{1,0,1}-c^{l-1,m,n}_{1,0,1}\right),
273 : * \right.\\
274 : * &\;\;\;\;\left.
275 : * \alpha_1\left(c^{l,m,n+1}_{2,0,0}-c^{l,m,n}_{2,0,0}\right),
276 : * \alpha_1\left(c^{l,m,n}_{2,0,0}-c^{l,m,n-1}_{2,0,0}\right)
277 : * \right),\\
278 : * \tilde{c}^{l,m,n}_{1,2,0}=\mathrm{minmod}
279 : * &\left(c_{1,2,0}^{l,m,n},
280 : * \alpha_1\left(c^{l+1,m,n}_{0,2,0}-c^{l,m,n}_{0,2,0}\right),
281 : * \alpha_1\left(c^{l,m,n}_{0,2,0}-c^{l-1,m,n}_{0,2,0}\right),
282 : * \right.\\
283 : * &\;\;\;\;\left.
284 : * \alpha_2\left(c^{l,m+1,n}_{1,1,0}-c^{l,m,n}_{1,1,0}\right),
285 : * \alpha_2\left(c^{l,m,n}_{1,1,0}-c^{l,m-1,n}_{1,1,0}\right)
286 : * \right),\\
287 : * \tilde{c}^{l,m,n}_{1,0,2}=\mathrm{minmod}
288 : * &\left(c_{1,0,2}^{l,m,n},
289 : * \alpha_1\left(c^{l+1,m,n}_{0,0,2}-c^{l,m,n}_{0,0,2}\right),
290 : * \alpha_1\left(c^{l,m,n}_{0,0,2}-c^{l-1,m,n}_{0,0,2}\right),
291 : * \right.\\
292 : * &\;\;\;\;\left.
293 : * \alpha_2\left(c^{l,m,n+1}_{1,0,1}-c^{l,m,n}_{1,0,1}\right),
294 : * \alpha_2\left(c^{l,m,n}_{1,0,1}-c^{l,m,n-1}_{1,0,1}\right)
295 : * \right),\\
296 : * \tilde{c}^{l,m,n}_{0,1,2}=\mathrm{minmod}
297 : * &\left(c_{0,1,2}^{l,m,n},
298 : * \alpha_1\left(c^{l,m+1,n}_{0,0,2}-c^{l,m,n}_{0,0,2}\right),
299 : * \alpha_1\left(c^{l,m,n}_{0,0,2}-c^{l,m-1,n}_{0,0,2}\right),
300 : * \right. \\
301 : * &\;\;\;\;\left.
302 : * \alpha_2\left(c^{l,m,n+1}_{0,1,1}-c^{l,m,n}_{0,1,1}\right),
303 : * \alpha_2\left(c^{l,m,n}_{0,1,1}-c^{l,m,n-1}_{0,1,1}\right)
304 : * \right),\\
305 : * \tilde{c}^{l,m,n}_{0,2,1}=\mathrm{minmod}
306 : * &\left(c_{0,2,1}^{l,m,n},
307 : * \alpha_2\left(c^{l,m+1,n}_{0,1,1}-c^{l,m,n}_{0,1,1}\right),
308 : * \alpha_2\left(c^{l,m,n}_{0,1,1}-c^{l,m-1,n}_{0,1,1}\right),
309 : * \right.\\
310 : * &\;\;\;\;\left.
311 : * \alpha_1\left(c^{l,m,n+1}_{0,2,0}-c^{l,m,n}_{0,2,0}\right),
312 : * \alpha_1\left(c^{l,m,n}_{0,2,0}-c^{l,m,n-1}_{0,2,0}\right)
313 : * \right),
314 : * \f}
315 : * \f{align*}{
316 : * \tilde{c}^{l,m,n}_{2,0,0}=\mathrm{minmod}
317 : * &\left(c_{2,0,0}^{l,m,n},
318 : * \alpha_2\left(c^{l+1,m,n}_{1,0,0}-c^{l,m,n}_{1,0,0}\right),
319 : * \alpha_2\left(c^{l,m,n}_{1,0,0}-c^{l-1,m,n}_{1,0,0}\right)
320 : * \right),\\
321 : * \tilde{c}^{l,m,n}_{0,2,0}=\mathrm{minmod}
322 : * &\left(c_{0,2,0}^{l,m,n},
323 : * \alpha_2\left(c^{l,m+1,n}_{0,1,0}-c^{l,m,n}_{0,1,0}\right),
324 : * \alpha_2\left(c^{l,m,n}_{0,1,0}-c^{l,m-1,n}_{0,1,0}\right)
325 : * \right),\\
326 : * \tilde{c}^{l,m,n}_{0,0,2}=\mathrm{minmod}
327 : * &\left(c_{0,0,2}^{l,m,n},
328 : * \alpha_2\left(c^{l,m,n+1}_{0,0,1}-c^{l,m,n}_{0,0,1}\right),
329 : * \alpha_2\left(c^{l,m,n}_{0,0,1}-c^{l,m,n-1}_{0,0,1}\right)
330 : * \right),\\
331 : * \tilde{c}^{l,m,n}_{1,1,1}=\mathrm{minmod}
332 : * &\left(c_{1,1,1}^{l,m,n},
333 : * \alpha_1\left(c^{l+1,m,n}_{0,1,1}-c^{l,m,n}_{0,1,1}\right),
334 : * \alpha_1\left(c^{l,m,n}_{0,1,1}-c^{l-1,m,n}_{0,1,1}\right),
335 : * \right.\\
336 : * &\;\;\;\;
337 : * \alpha_1\left(c^{l,m+1,n}_{1,0,1}-c^{l,m,n}_{1,0,1}\right),
338 : * \alpha_1\left(c^{l,m,n}_{1,0,1}-c^{l,m-1,n}_{1,0,1}\right),\\
339 : * &\;\;\;\;\left.
340 : * \alpha_1\left(c^{l,m,n+1}_{1,1,0}-c^{l,m,n}_{1,1,0}\right),
341 : * \alpha_1\left(c^{l,m,n}_{1,1,0}-c^{l,m,n-1}_{1,1,0}\right)
342 : * \right),\\
343 : * \tilde{c}^{l,m,n}_{1,1,0}=\mathrm{minmod}
344 : * &\left(c_{1,1,0}^{l,m,n},
345 : * \alpha_1\left(c^{l+1,m,n}_{0,1,0}-c^{l,m,n}_{0,1,0}\right),
346 : * \alpha_1\left(c^{l,m,n}_{0,1,0}-c^{l-1,m,n}_{0,1,0}\right),
347 : * \right.\\
348 : * &\;\;\;\;\left.
349 : * \alpha_1\left(c^{l,m+1,n}_{1,0,0}-c^{l,m,n}_{1,0,0}\right),
350 : * \alpha_1\left(c^{l,m,n}_{1,0,0}-c^{l,m-1,n}_{1,0,0}\right),
351 : * \right),\\
352 : * \tilde{c}^{l,m,n}_{1,0,1}=\mathrm{minmod}
353 : * &\left(c_{1,0,1}^{l,m,n},
354 : * \alpha_1\left(c^{l+1,m,n}_{0,0,1}-c^{l,m,n}_{0,0,1}\right),
355 : * \alpha_1\left(c^{l,m,n}_{0,0,1}-c^{l-1,m,n}_{0,0,1}\right),
356 : * \right.\\
357 : * &\;\;\;\;\left.
358 : * \alpha_1\left(c^{l,m,n+1}_{1,0,0}-c^{l,m,n}_{1,0,0}\right),
359 : * \alpha_1\left(c^{l,m,n}_{1,0,0}-c^{l,m,n-1}_{1,0,0}\right)
360 : * \right),\\
361 : * \tilde{c}^{l,m,n}_{0,1,1}=\mathrm{minmod}
362 : * &\left(c_{0,1,1}^{l,m,n},
363 : * \alpha_1\left(c^{l,m+1,n}_{0,0,1}-c^{l,m,n}_{0,0,1}\right),
364 : * \alpha_1\left(c^{l,m,n}_{0,0,1}-c^{l,m-1,n}_{0,0,1}\right),
365 : * \right.\\
366 : * &\;\;\;\;\left.
367 : * \alpha_1\left(c^{l,m,n+1}_{0,1,0}-c^{l,m,n}_{0,1,0}\right),
368 : * \alpha_1\left(c^{l,m,n}_{0,1,0}-c^{l,m,n-1}_{0,1,0}\right)
369 : * \right),\\
370 : * \tilde{c}^{l,m,n}_{1,0,0}=\mathrm{minmod}
371 : * &\left(c_{1,0,0}^{l,m,n},
372 : * \alpha_1\left(c^{l+1,m,n}_{0,0,0}-c^{l,m,n}_{0,0,0}\right),
373 : * \alpha_1\left(c^{l,m,n}_{0,0,0}-c^{l-1,m,n}_{0,0,0}\right)
374 : * \right),\\
375 : * \tilde{c}^{l,m,n}_{0,1,0}=\mathrm{minmod}
376 : * &\left(c_{0,1,0}^{l,m,n},
377 : * \alpha_1\left(c^{l,m+1,n}_{0,0,0}-c^{l,m,n}_{0,0,0}\right),
378 : * \alpha_1\left(c^{l,m,n}_{0,0,0}-c^{l,m-1,n}_{0,0,0}\right)
379 : * \right),\\
380 : * \tilde{c}^{l,m,n}_{0,0,1}=\mathrm{minmod}
381 : * &\left(c_{0,0,1}^{l,m,n},
382 : * \alpha_1\left(c^{l,m,n+1}_{0,0,0}-c^{l,m,n}_{0,0,0}\right),
383 : * \alpha_1\left(c^{l,m,n}_{0,0,0}-c^{l,m,n-1}_{0,0,0}\right)
384 : * \right),
385 : * \f}
386 : *
387 : *
388 : * The algorithm to perform the limiting is as follows:
389 : *
390 : * - limit \f$c_{2,2,2}\f$ (i.e. \f$c_{2,2,2}\leftarrow\tilde{c}_{2,2,2}\f$), if
391 : * not changed, stop
392 : * - limit \f$c_{2,2,1}\f$, \f$c_{2,1,2}\f$, and \f$c_{1,2,2}\f$, if all not
393 : * changed, stop
394 : * - limit \f$c_{2,2,0}\f$, \f$c_{2,0,2}\f$, and \f$c_{0,2,2}\f$, if all not
395 : * changed, stop
396 : * - limit \f$c_{2,1,1}\f$, \f$c_{1,2,1}\f$, and \f$c_{1,1,2}\f$, if all not
397 : * changed, stop
398 : * - limit \f$c_{2,1,0}\f$, \f$c_{2,0,1}\f$, \f$c_{1,2,0}\f$, \f$c_{1,0,2}\f$,
399 : * \f$c_{0,1,2}\f$, and \f$c_{0,2,1}\f$, if all not changed, stop
400 : * - limit \f$c_{2,0,0}\f$, \f$c_{0,2,0}\f$, and \f$c_{0,0,2}\f$, if all not
401 : * changed, stop
402 : * - limit \f$c_{1,1,1}\f$, if not changed, stop
403 : * - limit \f$c_{1,1,0}\f$, \f$c_{1,0,1}\f$, and \f$c_{0,1,1}\f$, if all not
404 : * changed, stop
405 : * - limit \f$c_{1,0,0}\f$, \f$c_{0,1,0}\f$, and \f$c_{0,0,1}\f$, if all not
406 : * changed, stop
407 : *
408 : * The 1d and 2d implementations are straightforward restrictions of the
409 : * described algorithm.
410 : *
411 : * #### Limitations:
412 : *
413 : * - We currently recompute the spectral coefficients for local use, after
414 : * having computed these same coefficients for sending to the neighbors. We
415 : * should be able to avoid this either by storing the coefficients in the
416 : * DataBox or by allowing the limiters' `packaged_data` function to return an
417 : * object to be passed as an additional argument to the `operator()` (still
418 : * would need to be stored in the DataBox).
419 : * - We cannot handle the case where neighbors have more/fewer
420 : * coefficients. In the case that the neighbor has more coefficients we could
421 : * just ignore the higher coefficients. In the case that the neighbor has
422 : * fewer coefficients we have a few choices.
423 : * - Having a different number of collocation points in different directions is
424 : * not supported. However, it is straightforward to handle this case. The
425 : * outermost loop should be in the direction with the most collocation points,
426 : * while the inner most loop should be over the direction with the fewest
427 : * collocation points. The highest to lowest coefficients can then be limited
428 : * appropriately again.
429 : * - h-refinement is not supported, but there is one reasonably straightforward
430 : * implementation that may work. In this case we would ignore refinement that
431 : * is not in the direction of the neighbor, treating the element as
432 : * simply having multiple neighbors in that direction. The only change would
433 : * be accounting for different refinement in the direction of the neighor,
434 : * which should be easy to add since the differences in coefficients in
435 : * Eq.\f$\ref{eq:krivodonova 3d minmod}\f$ will just be multiplied by
436 : * non-unity factors.
437 : */
438 : template <size_t VolumeDim, typename... Tags>
439 1 : class Krivodonova<VolumeDim, tmpl::list<Tags...>> {
440 : public:
441 : /*!
442 : * \brief The \f$\alpha_i\f$ values in the Krivodonova algorithm.
443 : */
444 1 : struct Alphas {
445 0 : using type = std::array<
446 : double, Spectral::maximum_number_of_points<Spectral::Basis::Legendre>>;
447 0 : static constexpr Options::String help = {
448 : "The alpha parameters of the Krivodonova limiter"};
449 : };
450 : /*!
451 : * \brief Turn the limiter off
452 : *
453 : * This option exists to temporarily disable the limiter for debugging
454 : * purposes. For problems where limiting is not needed, the preferred
455 : * approach is to not compile the limiter into the executable.
456 : */
457 1 : struct DisableForDebugging {
458 0 : using type = bool;
459 0 : static type suggested_value() { return false; }
460 0 : static constexpr Options::String help = {"Disable the limiter"};
461 : };
462 :
463 0 : using options = tmpl::list<Alphas, DisableForDebugging>;
464 0 : static constexpr Options::String help = {
465 : "The hierarchical limiter of Krivodonova.\n\n"
466 : "This limiter works by limiting the highest modal "
467 : "coefficients/derivatives using an aggressive minmod approach, "
468 : "decreasing in modal coefficient order until no more limiting is "
469 : "necessary."};
470 :
471 0 : explicit Krivodonova(
472 : std::array<double,
473 : Spectral::maximum_number_of_points<Spectral::Basis::Legendre>>
474 : alphas,
475 : bool disable_for_debugging = false, const Options::Context& context = {});
476 :
477 0 : Krivodonova() = default;
478 0 : Krivodonova(const Krivodonova&) = default;
479 0 : Krivodonova& operator=(const Krivodonova&) = default;
480 0 : Krivodonova(Krivodonova&&) = default;
481 0 : Krivodonova& operator=(Krivodonova&&) = default;
482 0 : ~Krivodonova() = default;
483 :
484 : // NOLINTNEXTLINE(google-runtime-references)
485 0 : void pup(PUP::er& p);
486 :
487 0 : bool operator==(const Krivodonova& rhs) const;
488 :
489 0 : struct PackagedData {
490 0 : Variables<tmpl::list<::Tags::Modal<Tags>...>> modal_volume_data;
491 0 : Mesh<VolumeDim> mesh;
492 :
493 : // NOLINTNEXTLINE(google-runtime-references)
494 0 : void pup(PUP::er& p) {
495 : p | modal_volume_data;
496 : p | mesh;
497 : }
498 : };
499 :
500 0 : using package_argument_tags =
501 : tmpl::list<Tags..., domain::Tags::Mesh<VolumeDim>>;
502 :
503 : /// \brief Package data for sending to neighbor elements.
504 1 : void package_data(gsl::not_null<PackagedData*> packaged_data,
505 : const typename Tags::type&... tensors,
506 : const Mesh<VolumeDim>& mesh,
507 : const OrientationMap<VolumeDim>& orientation_map) const;
508 :
509 0 : using limit_tags = tmpl::list<Tags...>;
510 0 : using limit_argument_tags = tmpl::list<domain::Tags::Element<VolumeDim>,
511 : domain::Tags::Mesh<VolumeDim>>;
512 :
513 0 : bool operator()(
514 : const gsl::not_null<std::add_pointer_t<typename Tags::type>>... tensors,
515 : const Element<VolumeDim>& element, const Mesh<VolumeDim>& mesh,
516 : const std::unordered_map<DirectionalId<VolumeDim>, PackagedData,
517 : boost::hash<DirectionalId<VolumeDim>>>&
518 : neighbor_data) const;
519 :
520 : private:
521 : template <typename Tag>
522 0 : char limit_one_tensor(
523 : gsl::not_null<Variables<tmpl::list<::Tags::Modal<Tags>...>>*> coeffs_self,
524 : gsl::not_null<bool*> limited_any_component, const Mesh<1>& mesh,
525 : const std::unordered_map<DirectionalId<1>, PackagedData,
526 : boost::hash<DirectionalId<1>>>& neighbor_data)
527 : const;
528 : template <typename Tag>
529 0 : char limit_one_tensor(
530 : gsl::not_null<Variables<tmpl::list<::Tags::Modal<Tags>...>>*> coeffs_self,
531 : gsl::not_null<bool*> limited_any_component, const Mesh<2>& mesh,
532 : const std::unordered_map<DirectionalId<2>, PackagedData,
533 : boost::hash<DirectionalId<2>>>& neighbor_data)
534 : const;
535 : template <typename Tag>
536 0 : char limit_one_tensor(
537 : gsl::not_null<Variables<tmpl::list<::Tags::Modal<Tags>...>>*> coeffs_self,
538 : gsl::not_null<bool*> limited_any_component, const Mesh<3>& mesh,
539 : const std::unordered_map<DirectionalId<3>, PackagedData,
540 : boost::hash<DirectionalId<3>>>& neighbor_data)
541 : const;
542 :
543 : template <typename Tag, size_t Dim>
544 0 : char fill_variables_tag_with_spectral_coeffs(
545 : gsl::not_null<Variables<tmpl::list<::Tags::Modal<Tags>...>>*>
546 : modal_coeffs,
547 : const typename Tag::type& nodal_tensor, const Mesh<Dim>& mesh) const;
548 :
549 : std::array<double,
550 : Spectral::maximum_number_of_points<Spectral::Basis::Legendre>>
551 0 : alphas_ = make_array<
552 : Spectral::maximum_number_of_points<Spectral::Basis::Legendre>>(
553 : std::numeric_limits<double>::signaling_NaN());
554 0 : bool disable_for_debugging_{false};
555 : };
556 :
557 : template <size_t VolumeDim, typename... Tags>
558 : Krivodonova<VolumeDim, tmpl::list<Tags...>>::Krivodonova(
559 : std::array<double,
560 : Spectral::maximum_number_of_points<Spectral::Basis::Legendre>>
561 : alphas,
562 : bool disable_for_debugging, const Options::Context& context)
563 : : alphas_(alphas), disable_for_debugging_(disable_for_debugging) {
564 : // See the main documentation for an explanation of why these bounds are
565 : // different from those of Krivodonova 2007
566 : if (alg::any_of(alphas_,
567 : [](const double t) { return t > 1.0 or t <= 0.0; })) {
568 : PARSE_ERROR(context,
569 : "The alphas in the Krivodonova limiter must be in the range "
570 : "(0,1].");
571 : }
572 : }
573 :
574 : template <size_t VolumeDim, typename... Tags>
575 : void Krivodonova<VolumeDim, tmpl::list<Tags...>>::package_data(
576 : const gsl::not_null<PackagedData*> packaged_data,
577 : const typename Tags::type&... tensors, const Mesh<VolumeDim>& mesh,
578 : const OrientationMap<VolumeDim>& orientation_map) const {
579 : if (UNLIKELY(disable_for_debugging_)) {
580 : // Do not initialize packaged_data
581 : return;
582 : }
583 :
584 : packaged_data->modal_volume_data.initialize(mesh.number_of_grid_points());
585 : // perform nodal coefficients to modal coefficients transformation on each
586 : // tensor component
587 : expand_pack(fill_variables_tag_with_spectral_coeffs<Tags>(
588 : &(packaged_data->modal_volume_data), tensors, mesh)...);
589 :
590 : packaged_data->modal_volume_data = orient_variables(
591 : packaged_data->modal_volume_data, mesh.extents(), orientation_map);
592 :
593 : // This doesn't currently do anything because the mesh is assumed to be the
594 : // same in all dimensions
595 : packaged_data->mesh = orientation_map(mesh);
596 : }
597 :
598 : template <size_t VolumeDim, typename... Tags>
599 : bool Krivodonova<VolumeDim, tmpl::list<Tags...>>::operator()(
600 : const gsl::not_null<std::add_pointer_t<typename Tags::type>>... tensors,
601 : const Element<VolumeDim>& element, const Mesh<VolumeDim>& mesh,
602 : const std::unordered_map<DirectionalId<VolumeDim>, PackagedData,
603 : boost::hash<DirectionalId<VolumeDim>>>&
604 : neighbor_data) const {
605 : if (UNLIKELY(disable_for_debugging_)) {
606 : // Do not modify input tensors
607 : return false;
608 : }
609 : if (UNLIKELY(mesh != Mesh<VolumeDim>(mesh.extents()[0], mesh.basis()[0],
610 : mesh.quadrature()[0]))) {
611 : ERROR(
612 : "The Krivodonova limiter does not yet support non-uniform number of "
613 : "collocation points, bases, and quadrature in each direction. The "
614 : "mesh is: "
615 : << mesh);
616 : }
617 : if (UNLIKELY(
618 : alg::any_of(element.neighbors(), [](const auto& direction_neighbors) {
619 : return direction_neighbors.second.size() != 1;
620 : }))) {
621 : ERROR("The Krivodonova limiter does not yet support h-refinement");
622 : }
623 : alg::for_each(neighbor_data, [&mesh](const auto& id_packaged_data) {
624 : if (UNLIKELY(id_packaged_data.second.mesh != mesh)) {
625 : ERROR(
626 : "The Krivodonova limiter does not yet support differing meshes "
627 : "between neighbors. Self mesh is: "
628 : << mesh << " neighbor mesh is: " << id_packaged_data.second.mesh);
629 : }
630 : });
631 :
632 : // Compute local modal coefficients
633 : Variables<tmpl::list<::Tags::Modal<Tags>...>> coeffs_self(
634 : mesh.number_of_grid_points(), 0.0);
635 : expand_pack(fill_variables_tag_with_spectral_coeffs<Tags>(&coeffs_self,
636 : *tensors, mesh)...);
637 :
638 : // Perform the limiting on the modal coefficients
639 : bool limited_any_component = false;
640 : expand_pack(limit_one_tensor<::Tags::Modal<Tags>>(
641 : make_not_null(&coeffs_self), make_not_null(&limited_any_component), mesh,
642 : neighbor_data)...);
643 :
644 : // transform back to nodal coefficients
645 : const auto wrap_copy_nodal_coeffs =
646 : [&mesh, &coeffs_self](auto tag, const auto tensor) {
647 : auto& coeffs_tensor = get<decltype(tag)>(coeffs_self);
648 : auto tensor_it = tensor->begin();
649 : for (auto coeffs_it = coeffs_tensor.begin();
650 : coeffs_it != coeffs_tensor.end();
651 : (void)++coeffs_it, (void)++tensor_it) {
652 : to_nodal_coefficients(make_not_null(&*tensor_it), *coeffs_it, mesh);
653 : }
654 : return '0';
655 : };
656 : expand_pack(wrap_copy_nodal_coeffs(::Tags::Modal<Tags>{}, tensors)...);
657 :
658 : return limited_any_component;
659 : }
660 :
661 : template <size_t VolumeDim, typename... Tags>
662 : template <typename Tag>
663 0 : char Krivodonova<VolumeDim, tmpl::list<Tags...>>::limit_one_tensor(
664 : const gsl::not_null<Variables<tmpl::list<::Tags::Modal<Tags>...>>*>
665 : coeffs_self,
666 : const gsl::not_null<bool*> limited_any_component, const Mesh<1>& mesh,
667 : const std::unordered_map<DirectionalId<1>, PackagedData,
668 : boost::hash<DirectionalId<1>>>& neighbor_data)
669 : const {
670 : using tensor_type = typename Tag::type;
671 : for (size_t storage_index = 0; storage_index < tensor_type::size();
672 : ++storage_index) {
673 : // for each coefficient...
674 : for (size_t i = mesh.extents()[0] - 1; i > 0; --i) {
675 : auto& self_coeffs = get<Tag>(*coeffs_self)[storage_index];
676 : double min_abs_coeff = std::abs(self_coeffs[i]);
677 : const double sgn_of_coeff = sgn(self_coeffs[i]);
678 : bool sgns_all_equal = true;
679 : for (const auto& kv : neighbor_data) {
680 : const auto& neighbor_coeffs =
681 : get<Tag>(kv.second.modal_volume_data)[storage_index];
682 : const double tmp = kv.first.direction.sign() * gsl::at(alphas_, i) *
683 : (neighbor_coeffs[i - 1] - self_coeffs[i - 1]);
684 :
685 : min_abs_coeff = std::min(min_abs_coeff, std::abs(tmp));
686 : sgns_all_equal &= sgn(tmp) == sgn_of_coeff;
687 : if (not sgns_all_equal) {
688 : self_coeffs[i] = 0.0;
689 : break;
690 : }
691 : }
692 : if (sgns_all_equal) {
693 : const double tmp = sgn_of_coeff * min_abs_coeff;
694 : if (tmp == self_coeffs[i]) {
695 : break;
696 : }
697 : *limited_any_component |= true;
698 : self_coeffs[i] = tmp;
699 : }
700 : }
701 : }
702 : return '0';
703 : }
704 :
705 : template <size_t VolumeDim, typename... Tags>
706 : template <typename Tag>
707 0 : char Krivodonova<VolumeDim, tmpl::list<Tags...>>::limit_one_tensor(
708 : const gsl::not_null<Variables<tmpl::list<::Tags::Modal<Tags>...>>*>
709 : coeffs_self,
710 : const gsl::not_null<bool*> limited_any_component, const Mesh<2>& mesh,
711 : const std::unordered_map<DirectionalId<2>, PackagedData,
712 : boost::hash<DirectionalId<2>>>& neighbor_data)
713 : const {
714 : using tensor_type = typename Tag::type;
715 : const auto minmod = [&coeffs_self, &mesh, &neighbor_data, this](
716 : const size_t local_i, const size_t local_j,
717 : const size_t local_tensor_storage_index) {
718 : const auto& self_coeffs =
719 : get<Tag>(*coeffs_self)[local_tensor_storage_index];
720 : double min_abs_coeff = std::abs(
721 : self_coeffs[mesh.storage_index(Index<VolumeDim>{local_i, local_j})]);
722 : const double sgn_of_coeff = sgn(
723 : self_coeffs[mesh.storage_index(Index<VolumeDim>{local_i, local_j})]);
724 : for (const auto& kv : neighbor_data) {
725 : const Direction<2>& dir = kv.first.direction;
726 : const auto& neighbor_coeffs =
727 : get<Tag>(kv.second.modal_volume_data)[local_tensor_storage_index];
728 : const size_t index_i =
729 : dir.axis() == Direction<VolumeDim>::Axis::Xi ? local_i - 1 : local_i;
730 : const size_t index_j =
731 : dir.axis() == Direction<VolumeDim>::Axis::Eta ? local_j - 1 : local_j;
732 : // skip neighbors where we cannot compute a finite difference in that
733 : // direction because we are already at the lowest coefficient.
734 : if (index_i == std::numeric_limits<size_t>::max() or
735 : index_j == std::numeric_limits<size_t>::max()) {
736 : continue;
737 : }
738 : const size_t alpha_index =
739 : dir.axis() == Direction<VolumeDim>::Axis::Xi ? local_i : local_j;
740 : const double tmp =
741 : dir.sign() * gsl::at(alphas_, alpha_index) *
742 : (neighbor_coeffs[mesh.storage_index(
743 : Index<VolumeDim>{index_i, index_j})] -
744 : self_coeffs[mesh.storage_index(Index<VolumeDim>{index_i, index_j})]);
745 :
746 : min_abs_coeff = std::min(min_abs_coeff, std::abs(tmp));
747 : if (sgn(tmp) != sgn_of_coeff) {
748 : return 0.0;
749 : }
750 : }
751 : return sgn_of_coeff * min_abs_coeff;
752 : };
753 :
754 : for (size_t tensor_storage_index = 0;
755 : tensor_storage_index < tensor_type::size(); ++tensor_storage_index) {
756 : // for each coefficient...
757 : for (size_t i = mesh.extents()[0] - 1; i > 0; --i) {
758 : for (size_t j = i; j < mesh.extents()[1]; --j) {
759 : // Check if we are done limiting, and if so we move on to the next
760 : // tensor component.
761 : auto& self_coeffs = get<Tag>(*coeffs_self)[tensor_storage_index];
762 : // We treat the different cases separately to reduce the number of
763 : // times we call minmod, not because it is required for correctness.
764 : if (UNLIKELY(i == j)) {
765 : const double tmp = minmod(i, j, tensor_storage_index);
766 : if (tmp == self_coeffs[mesh.storage_index(Index<VolumeDim>{i, j})]) {
767 : goto next_tensor_index;
768 : }
769 : *limited_any_component |= true;
770 : self_coeffs[mesh.storage_index(Index<VolumeDim>{i, j})] = tmp;
771 : } else {
772 : const double tmp_ij = minmod(i, j, tensor_storage_index);
773 : const double tmp_ji = minmod(j, i, tensor_storage_index);
774 : if (tmp_ij ==
775 : self_coeffs[mesh.storage_index(Index<VolumeDim>{i, j})] and
776 : tmp_ji ==
777 : self_coeffs[mesh.storage_index(Index<VolumeDim>{j, i})]) {
778 : goto next_tensor_index;
779 : }
780 : *limited_any_component |= true;
781 : self_coeffs[mesh.storage_index(Index<VolumeDim>{i, j})] = tmp_ij;
782 : self_coeffs[mesh.storage_index(Index<VolumeDim>{j, i})] = tmp_ji;
783 : }
784 : }
785 : }
786 : next_tensor_index:
787 : continue;
788 : }
789 : return '0';
790 : }
791 :
792 : template <size_t VolumeDim, typename... Tags>
793 : template <typename Tag>
794 0 : char Krivodonova<VolumeDim, tmpl::list<Tags...>>::limit_one_tensor(
795 : const gsl::not_null<Variables<tmpl::list<::Tags::Modal<Tags>...>>*>
796 : coeffs_self,
797 : const gsl::not_null<bool*> limited_any_component, const Mesh<3>& mesh,
798 : const std::unordered_map<DirectionalId<3>, PackagedData,
799 : boost::hash<DirectionalId<3>>>& neighbor_data)
800 : const {
801 : using tensor_type = typename Tag::type;
802 : const auto minmod = [&coeffs_self, &mesh, &neighbor_data, this](
803 : const size_t local_i, const size_t local_j,
804 : const size_t local_k,
805 : const size_t local_tensor_storage_index) {
806 : const auto& self_coeffs =
807 : get<Tag>(*coeffs_self)[local_tensor_storage_index];
808 : double min_abs_coeff = std::abs(self_coeffs[mesh.storage_index(
809 : Index<VolumeDim>{local_i, local_j, local_k})]);
810 : const double sgn_of_coeff = sgn(self_coeffs[mesh.storage_index(
811 : Index<VolumeDim>{local_i, local_j, local_k})]);
812 : for (const auto& kv : neighbor_data) {
813 : const Direction<3>& dir = kv.first.direction;
814 : const auto& neighbor_coeffs =
815 : get<Tag>(kv.second.modal_volume_data)[local_tensor_storage_index];
816 : const size_t index_i =
817 : dir.axis() == Direction<VolumeDim>::Axis::Xi ? local_i - 1 : local_i;
818 : const size_t index_j =
819 : dir.axis() == Direction<VolumeDim>::Axis::Eta ? local_j - 1 : local_j;
820 : const size_t index_k = dir.axis() == Direction<VolumeDim>::Axis::Zeta
821 : ? local_k - 1
822 : : local_k;
823 : // skip neighbors where we cannot compute a finite difference in that
824 : // direction because we are already at the lowest coefficient.
825 : if (index_i == std::numeric_limits<size_t>::max() or
826 : index_j == std::numeric_limits<size_t>::max() or
827 : index_k == std::numeric_limits<size_t>::max()) {
828 : continue;
829 : }
830 : const size_t alpha_index =
831 : dir.axis() == Direction<VolumeDim>::Axis::Xi
832 : ? local_i
833 : : dir.axis() == Direction<VolumeDim>::Axis::Eta ? local_j
834 : : local_k;
835 : const double tmp = dir.sign() * gsl::at(alphas_, alpha_index) *
836 : (neighbor_coeffs[mesh.storage_index(
837 : Index<VolumeDim>{index_i, index_j, index_k})] -
838 : self_coeffs[mesh.storage_index(
839 : Index<VolumeDim>{index_i, index_j, index_k})]);
840 :
841 : min_abs_coeff = std::min(min_abs_coeff, std::abs(tmp));
842 : if (sgn(tmp) != sgn_of_coeff) {
843 : return 0.0;
844 : }
845 : }
846 : return sgn_of_coeff * min_abs_coeff;
847 : };
848 :
849 : for (size_t tensor_storage_index = 0;
850 : tensor_storage_index < tensor_type::size(); ++tensor_storage_index) {
851 : // for each coefficient...
852 : for (size_t i = mesh.extents()[0] - 1; i > 0; --i) {
853 : for (size_t j = i; j < mesh.extents()[1]; --j) {
854 : for (size_t k = j; k < mesh.extents()[2]; --k) {
855 : // Check if we are done limiting, and if so we move on to the next
856 : // tensor component.
857 : auto& self_coeffs = get<Tag>(*coeffs_self)[tensor_storage_index];
858 : // We treat the different cases separately to reduce the number of
859 : // times we call minmod, not because it is required for correctness.
860 : //
861 : // Note that the case `i == k and i != j` cannot be encountered since
862 : // the loop bounds are `i >= j >= k`.
863 : if (UNLIKELY(i == j and i == k)) {
864 : const double tmp = minmod(i, j, k, tensor_storage_index);
865 : if (tmp ==
866 : self_coeffs[mesh.storage_index(Index<VolumeDim>{i, j, k})]) {
867 : goto next_tensor_index;
868 : }
869 : *limited_any_component |= true;
870 : self_coeffs[mesh.storage_index(Index<VolumeDim>{i, j, k})] = tmp;
871 : } else if (i == j and i != k) {
872 : const double tmp_ijk = minmod(i, j, k, tensor_storage_index);
873 : const double tmp_ikj = minmod(i, k, j, tensor_storage_index);
874 : const double tmp_kij = minmod(k, i, j, tensor_storage_index);
875 : if (tmp_ijk == self_coeffs[mesh.storage_index(
876 : Index<VolumeDim>{i, j, k})] and
877 : tmp_ikj == self_coeffs[mesh.storage_index(
878 : Index<VolumeDim>{i, k, j})] and
879 : tmp_kij == self_coeffs[mesh.storage_index(
880 : Index<VolumeDim>{k, i, j})]) {
881 : goto next_tensor_index;
882 : }
883 : *limited_any_component |= true;
884 : self_coeffs[mesh.storage_index(Index<VolumeDim>{i, j, k})] =
885 : tmp_ijk;
886 : self_coeffs[mesh.storage_index(Index<VolumeDim>{i, k, j})] =
887 : tmp_ikj;
888 : self_coeffs[mesh.storage_index(Index<VolumeDim>{k, i, j})] =
889 : tmp_kij;
890 : } else if (i != j and j == k) {
891 : const double tmp_ijk = minmod(i, j, k, tensor_storage_index);
892 : const double tmp_kij = minmod(k, i, j, tensor_storage_index);
893 : const double tmp_kji = minmod(k, j, i, tensor_storage_index);
894 : if (tmp_ijk == self_coeffs[mesh.storage_index(
895 : Index<VolumeDim>{i, j, k})] and
896 : tmp_kij == self_coeffs[mesh.storage_index(
897 : Index<VolumeDim>{k, i, j})] and
898 : tmp_kji == self_coeffs[mesh.storage_index(
899 : Index<VolumeDim>{k, j, i})]) {
900 : goto next_tensor_index;
901 : }
902 : *limited_any_component |= true;
903 : self_coeffs[mesh.storage_index(Index<VolumeDim>{i, j, k})] =
904 : tmp_ijk;
905 : self_coeffs[mesh.storage_index(Index<VolumeDim>{k, i, j})] =
906 : tmp_kij;
907 : self_coeffs[mesh.storage_index(Index<VolumeDim>{k, j, i})] =
908 : tmp_kji;
909 : } else {
910 : const double tmp_ijk = minmod(i, j, k, tensor_storage_index);
911 : const double tmp_jik = minmod(j, i, k, tensor_storage_index);
912 : const double tmp_ikj = minmod(i, k, j, tensor_storage_index);
913 : const double tmp_jki = minmod(j, k, i, tensor_storage_index);
914 : const double tmp_kij = minmod(k, i, j, tensor_storage_index);
915 : const double tmp_kji = minmod(k, j, i, tensor_storage_index);
916 : if (tmp_ijk == self_coeffs[mesh.storage_index(
917 : Index<VolumeDim>{i, j, k})] and
918 : tmp_jik == self_coeffs[mesh.storage_index(
919 : Index<VolumeDim>{j, i, k})] and
920 : tmp_ikj == self_coeffs[mesh.storage_index(
921 : Index<VolumeDim>{i, k, j})] and
922 : tmp_jki == self_coeffs[mesh.storage_index(
923 : Index<VolumeDim>{j, k, i})] and
924 : tmp_kij == self_coeffs[mesh.storage_index(
925 : Index<VolumeDim>{k, i, j})] and
926 : tmp_kji == self_coeffs[mesh.storage_index(
927 : Index<VolumeDim>{k, j, i})]) {
928 : goto next_tensor_index;
929 : }
930 : *limited_any_component |= true;
931 : self_coeffs[mesh.storage_index(Index<VolumeDim>{i, j, k})] =
932 : tmp_ijk;
933 : self_coeffs[mesh.storage_index(Index<VolumeDim>{j, i, k})] =
934 : tmp_jik;
935 : self_coeffs[mesh.storage_index(Index<VolumeDim>{i, k, j})] =
936 : tmp_ikj;
937 : self_coeffs[mesh.storage_index(Index<VolumeDim>{j, k, i})] =
938 : tmp_jki;
939 : self_coeffs[mesh.storage_index(Index<VolumeDim>{k, i, j})] =
940 : tmp_kij;
941 : self_coeffs[mesh.storage_index(Index<VolumeDim>{k, j, i})] =
942 : tmp_kji;
943 : }
944 : }
945 : }
946 : }
947 : next_tensor_index:
948 : continue;
949 : }
950 : return '0';
951 : }
952 :
953 : template <size_t VolumeDim, typename... Tags>
954 : template <typename Tag, size_t Dim>
955 0 : char Krivodonova<VolumeDim, tmpl::list<Tags...>>::
956 : fill_variables_tag_with_spectral_coeffs(
957 : const gsl::not_null<Variables<tmpl::list<::Tags::Modal<Tags>...>>*>
958 : modal_coeffs,
959 : const typename Tag::type& nodal_tensor, const Mesh<Dim>& mesh) const {
960 : auto& coeffs_tensor = get<::Tags::Modal<Tag>>(*modal_coeffs);
961 : auto tensor_it = nodal_tensor.begin();
962 : for (auto coeffs_it = coeffs_tensor.begin(); coeffs_it != coeffs_tensor.end();
963 : (void)++coeffs_it, (void)++tensor_it) {
964 : to_modal_coefficients(make_not_null(&*coeffs_it), *tensor_it, mesh);
965 : }
966 : return '0';
967 : }
968 :
969 : template <size_t VolumeDim, typename... Tags>
970 : void Krivodonova<VolumeDim, tmpl::list<Tags...>>::pup(PUP::er& p) {
971 : p | alphas_;
972 : p | disable_for_debugging_;
973 : }
974 :
975 : template <size_t VolumeDim, typename... Tags>
976 : bool Krivodonova<VolumeDim, tmpl::list<Tags...>>::operator==(
977 : const Krivodonova<VolumeDim, tmpl::list<Tags...>>& rhs) const {
978 : return alphas_ == rhs.alphas_ and
979 : disable_for_debugging_ == rhs.disable_for_debugging_;
980 : }
981 :
982 : template <size_t VolumeDim, typename... Tags>
983 0 : bool operator!=(const Krivodonova<VolumeDim, tmpl::list<Tags...>>& lhs,
984 : const Krivodonova<VolumeDim, tmpl::list<Tags...>>& rhs) {
985 : return not(lhs == rhs);
986 : }
987 : } // namespace Limiters
|
