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