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<
443 : double, Spectral::maximum_number_of_points<Spectral::Basis::Legendre>>;
444 0 : static constexpr Options::String help = {
445 : "The alpha parameters of the Krivodonova limiter"};
446 : };
447 : /*!
448 : * \brief Turn the limiter off
449 : *
450 : * This option exists to temporarily disable the limiter for debugging
451 : * purposes. For problems where limiting is not needed, the preferred
452 : * approach is to not compile the limiter into the executable.
453 : */
454 1 : struct DisableForDebugging {
455 0 : using type = bool;
456 0 : static type suggested_value() { return false; }
457 0 : static constexpr Options::String help = {"Disable the limiter"};
458 : };
459 :
460 0 : using options = tmpl::list<Alphas, DisableForDebugging>;
461 0 : static constexpr Options::String help = {
462 : "The hierarchical limiter of Krivodonova.\n\n"
463 : "This limiter works by limiting the highest modal "
464 : "coefficients/derivatives using an aggressive minmod approach, "
465 : "decreasing in modal coefficient order until no more limiting is "
466 : "necessary."};
467 :
468 0 : explicit Krivodonova(
469 : std::array<double,
470 : Spectral::maximum_number_of_points<Spectral::Basis::Legendre>>
471 : alphas,
472 : bool disable_for_debugging = false, const Options::Context& context = {});
473 :
474 0 : Krivodonova() = default;
475 0 : Krivodonova(const Krivodonova&) = default;
476 0 : Krivodonova& operator=(const Krivodonova&) = default;
477 0 : Krivodonova(Krivodonova&&) = default;
478 0 : Krivodonova& operator=(Krivodonova&&) = default;
479 0 : ~Krivodonova() = default;
480 :
481 : // NOLINTNEXTLINE(google-runtime-references)
482 0 : void pup(PUP::er& p);
483 :
484 0 : bool operator==(const Krivodonova& rhs) const;
485 :
486 0 : struct PackagedData {
487 0 : Variables<tmpl::list<::Tags::Modal<Tags>...>> modal_volume_data;
488 0 : Mesh<VolumeDim> mesh;
489 :
490 : // NOLINTNEXTLINE(google-runtime-references)
491 0 : void pup(PUP::er& p) {
492 : p | modal_volume_data;
493 : p | mesh;
494 : }
495 : };
496 :
497 0 : using package_argument_tags =
498 : tmpl::list<Tags..., domain::Tags::Mesh<VolumeDim>>;
499 :
500 : /// \brief Package data for sending to neighbor elements.
501 1 : void package_data(gsl::not_null<PackagedData*> packaged_data,
502 : const typename Tags::type&... tensors,
503 : const Mesh<VolumeDim>& mesh,
504 : const OrientationMap<VolumeDim>& orientation_map) const;
505 :
506 0 : using limit_tags = tmpl::list<Tags...>;
507 0 : using limit_argument_tags = tmpl::list<domain::Tags::Element<VolumeDim>,
508 : domain::Tags::Mesh<VolumeDim>>;
509 :
510 0 : bool operator()(
511 : const gsl::not_null<std::add_pointer_t<typename Tags::type>>... tensors,
512 : const Element<VolumeDim>& element, const Mesh<VolumeDim>& mesh,
513 : const std::unordered_map<DirectionalId<VolumeDim>, PackagedData,
514 : boost::hash<DirectionalId<VolumeDim>>>&
515 : neighbor_data) const;
516 :
517 : private:
518 : template <typename Tag>
519 0 : char limit_one_tensor(
520 : gsl::not_null<Variables<tmpl::list<::Tags::Modal<Tags>...>>*> coeffs_self,
521 : gsl::not_null<bool*> limited_any_component, const Mesh<1>& mesh,
522 : const std::unordered_map<DirectionalId<1>, PackagedData,
523 : boost::hash<DirectionalId<1>>>& neighbor_data)
524 : const;
525 : template <typename Tag>
526 0 : char limit_one_tensor(
527 : gsl::not_null<Variables<tmpl::list<::Tags::Modal<Tags>...>>*> coeffs_self,
528 : gsl::not_null<bool*> limited_any_component, const Mesh<2>& mesh,
529 : const std::unordered_map<DirectionalId<2>, PackagedData,
530 : boost::hash<DirectionalId<2>>>& neighbor_data)
531 : const;
532 : template <typename Tag>
533 0 : char limit_one_tensor(
534 : gsl::not_null<Variables<tmpl::list<::Tags::Modal<Tags>...>>*> coeffs_self,
535 : gsl::not_null<bool*> limited_any_component, const Mesh<3>& mesh,
536 : const std::unordered_map<DirectionalId<3>, PackagedData,
537 : boost::hash<DirectionalId<3>>>& neighbor_data)
538 : const;
539 :
540 : template <typename Tag, size_t Dim>
541 0 : char fill_variables_tag_with_spectral_coeffs(
542 : gsl::not_null<Variables<tmpl::list<::Tags::Modal<Tags>...>>*>
543 : modal_coeffs,
544 : const typename Tag::type& nodal_tensor, const Mesh<Dim>& mesh) const;
545 :
546 : std::array<double,
547 : Spectral::maximum_number_of_points<Spectral::Basis::Legendre>>
548 0 : alphas_ = make_array<
549 : Spectral::maximum_number_of_points<Spectral::Basis::Legendre>>(
550 : std::numeric_limits<double>::signaling_NaN());
551 0 : bool disable_for_debugging_{false};
552 : };
553 :
554 : template <size_t VolumeDim, typename... Tags>
555 : Krivodonova<VolumeDim, tmpl::list<Tags...>>::Krivodonova(
556 : std::array<double,
557 : Spectral::maximum_number_of_points<Spectral::Basis::Legendre>>
558 : alphas,
559 : bool disable_for_debugging, const Options::Context& context)
560 : : alphas_(alphas), disable_for_debugging_(disable_for_debugging) {
561 : // See the main documentation for an explanation of why these bounds are
562 : // different from those of Krivodonova 2007
563 : if (alg::any_of(alphas_,
564 : [](const double t) { return t > 1.0 or t <= 0.0; })) {
565 : PARSE_ERROR(context,
566 : "The alphas in the Krivodonova limiter must be in the range "
567 : "(0,1].");
568 : }
569 : }
570 :
571 : template <size_t VolumeDim, typename... Tags>
572 : void Krivodonova<VolumeDim, tmpl::list<Tags...>>::package_data(
573 : const gsl::not_null<PackagedData*> packaged_data,
574 : const typename Tags::type&... tensors, const Mesh<VolumeDim>& mesh,
575 : const OrientationMap<VolumeDim>& orientation_map) const {
576 : if (UNLIKELY(disable_for_debugging_)) {
577 : // Do not initialize packaged_data
578 : return;
579 : }
580 :
581 : packaged_data->modal_volume_data.initialize(mesh.number_of_grid_points());
582 : // perform nodal coefficients to modal coefficients transformation on each
583 : // tensor component
584 : expand_pack(fill_variables_tag_with_spectral_coeffs<Tags>(
585 : &(packaged_data->modal_volume_data), tensors, mesh)...);
586 :
587 : packaged_data->modal_volume_data = orient_variables(
588 : packaged_data->modal_volume_data, mesh.extents(), orientation_map);
589 :
590 : // This doesn't currently do anything because the mesh is assumed to be the
591 : // same in all dimensions
592 : packaged_data->mesh = orientation_map(mesh);
593 : }
594 :
595 : template <size_t VolumeDim, typename... Tags>
596 : bool Krivodonova<VolumeDim, tmpl::list<Tags...>>::operator()(
597 : const gsl::not_null<std::add_pointer_t<typename Tags::type>>... tensors,
598 : const Element<VolumeDim>& element, const Mesh<VolumeDim>& mesh,
599 : const std::unordered_map<DirectionalId<VolumeDim>, PackagedData,
600 : boost::hash<DirectionalId<VolumeDim>>>&
601 : neighbor_data) const {
602 : if (UNLIKELY(disable_for_debugging_)) {
603 : // Do not modify input tensors
604 : return false;
605 : }
606 : if (UNLIKELY(mesh != Mesh<VolumeDim>(mesh.extents()[0], mesh.basis()[0],
607 : mesh.quadrature()[0]))) {
608 : ERROR(
609 : "The Krivodonova limiter does not yet support non-uniform number of "
610 : "collocation points, bases, and quadrature in each direction. The "
611 : "mesh is: "
612 : << mesh);
613 : }
614 : if (UNLIKELY(
615 : alg::any_of(element.neighbors(), [](const auto& direction_neighbors) {
616 : return direction_neighbors.second.size() != 1;
617 : }))) {
618 : ERROR("The Krivodonova limiter does not yet support h-refinement");
619 : }
620 : alg::for_each(neighbor_data, [&mesh](const auto& id_packaged_data) {
621 : if (UNLIKELY(id_packaged_data.second.mesh != mesh)) {
622 : ERROR(
623 : "The Krivodonova limiter does not yet support differing meshes "
624 : "between neighbors. Self mesh is: "
625 : << mesh << " neighbor mesh is: " << id_packaged_data.second.mesh);
626 : }
627 : });
628 :
629 : // Compute local modal coefficients
630 : Variables<tmpl::list<::Tags::Modal<Tags>...>> coeffs_self(
631 : mesh.number_of_grid_points(), 0.0);
632 : expand_pack(fill_variables_tag_with_spectral_coeffs<Tags>(&coeffs_self,
633 : *tensors, mesh)...);
634 :
635 : // Perform the limiting on the modal coefficients
636 : bool limited_any_component = false;
637 : expand_pack(limit_one_tensor<::Tags::Modal<Tags>>(
638 : make_not_null(&coeffs_self), make_not_null(&limited_any_component), mesh,
639 : neighbor_data)...);
640 :
641 : // transform back to nodal coefficients
642 : const auto wrap_copy_nodal_coeffs =
643 : [&mesh, &coeffs_self](auto tag, const auto tensor) {
644 : auto& coeffs_tensor = get<decltype(tag)>(coeffs_self);
645 : auto tensor_it = tensor->begin();
646 : for (auto coeffs_it = coeffs_tensor.begin();
647 : coeffs_it != coeffs_tensor.end();
648 : (void)++coeffs_it, (void)++tensor_it) {
649 : to_nodal_coefficients(make_not_null(&*tensor_it), *coeffs_it, mesh);
650 : }
651 : return '0';
652 : };
653 : expand_pack(wrap_copy_nodal_coeffs(::Tags::Modal<Tags>{}, tensors)...);
654 :
655 : return limited_any_component;
656 : }
657 :
658 : template <size_t VolumeDim, typename... Tags>
659 : template <typename Tag>
660 0 : char Krivodonova<VolumeDim, tmpl::list<Tags...>>::limit_one_tensor(
661 : const gsl::not_null<Variables<tmpl::list<::Tags::Modal<Tags>...>>*>
662 : coeffs_self,
663 : const gsl::not_null<bool*> limited_any_component, const Mesh<1>& mesh,
664 : const std::unordered_map<DirectionalId<1>, PackagedData,
665 : boost::hash<DirectionalId<1>>>& neighbor_data)
666 : const {
667 : using tensor_type = typename Tag::type;
668 : for (size_t storage_index = 0; storage_index < tensor_type::size();
669 : ++storage_index) {
670 : // for each coefficient...
671 : for (size_t i = mesh.extents()[0] - 1; i > 0; --i) {
672 : auto& self_coeffs = get<Tag>(*coeffs_self)[storage_index];
673 : double min_abs_coeff = std::abs(self_coeffs[i]);
674 : const double sgn_of_coeff = sgn(self_coeffs[i]);
675 : bool sgns_all_equal = true;
676 : for (const auto& kv : neighbor_data) {
677 : const auto& neighbor_coeffs =
678 : get<Tag>(kv.second.modal_volume_data)[storage_index];
679 : const double tmp = kv.first.direction().sign() * gsl::at(alphas_, i) *
680 : (neighbor_coeffs[i - 1] - self_coeffs[i - 1]);
681 :
682 : min_abs_coeff = std::min(min_abs_coeff, std::abs(tmp));
683 : sgns_all_equal &= sgn(tmp) == sgn_of_coeff;
684 : if (not sgns_all_equal) {
685 : self_coeffs[i] = 0.0;
686 : break;
687 : }
688 : }
689 : if (sgns_all_equal) {
690 : const double tmp = sgn_of_coeff * min_abs_coeff;
691 : if (tmp == self_coeffs[i]) {
692 : break;
693 : }
694 : *limited_any_component |= true;
695 : self_coeffs[i] = tmp;
696 : }
697 : }
698 : }
699 : return '0';
700 : }
701 :
702 : template <size_t VolumeDim, typename... Tags>
703 : template <typename Tag>
704 0 : char Krivodonova<VolumeDim, tmpl::list<Tags...>>::limit_one_tensor(
705 : const gsl::not_null<Variables<tmpl::list<::Tags::Modal<Tags>...>>*>
706 : coeffs_self,
707 : const gsl::not_null<bool*> limited_any_component, const Mesh<2>& mesh,
708 : const std::unordered_map<DirectionalId<2>, PackagedData,
709 : boost::hash<DirectionalId<2>>>& neighbor_data)
710 : const {
711 : using tensor_type = typename Tag::type;
712 : const auto minmod = [&coeffs_self, &mesh, &neighbor_data, this](
713 : const size_t local_i, const size_t local_j,
714 : const size_t local_tensor_storage_index) {
715 : const auto& self_coeffs =
716 : get<Tag>(*coeffs_self)[local_tensor_storage_index];
717 : double min_abs_coeff = std::abs(
718 : self_coeffs[mesh.storage_index(Index<VolumeDim>{local_i, local_j})]);
719 : const double sgn_of_coeff = sgn(
720 : self_coeffs[mesh.storage_index(Index<VolumeDim>{local_i, local_j})]);
721 : for (const auto& kv : neighbor_data) {
722 : const Direction<2>& dir = kv.first.direction();
723 : const auto& neighbor_coeffs =
724 : get<Tag>(kv.second.modal_volume_data)[local_tensor_storage_index];
725 : const size_t index_i =
726 : dir.axis() == Direction<VolumeDim>::Axis::Xi ? local_i - 1 : local_i;
727 : const size_t index_j =
728 : dir.axis() == Direction<VolumeDim>::Axis::Eta ? local_j - 1 : local_j;
729 : // skip neighbors where we cannot compute a finite difference in that
730 : // direction because we are already at the lowest coefficient.
731 : if (index_i == std::numeric_limits<size_t>::max() or
732 : index_j == std::numeric_limits<size_t>::max()) {
733 : continue;
734 : }
735 : const size_t alpha_index =
736 : dir.axis() == Direction<VolumeDim>::Axis::Xi ? local_i : local_j;
737 : const double tmp =
738 : dir.sign() * gsl::at(alphas_, alpha_index) *
739 : (neighbor_coeffs[mesh.storage_index(
740 : Index<VolumeDim>{index_i, index_j})] -
741 : self_coeffs[mesh.storage_index(Index<VolumeDim>{index_i, index_j})]);
742 :
743 : min_abs_coeff = std::min(min_abs_coeff, std::abs(tmp));
744 : if (sgn(tmp) != sgn_of_coeff) {
745 : return 0.0;
746 : }
747 : }
748 : return sgn_of_coeff * min_abs_coeff;
749 : };
750 :
751 : for (size_t tensor_storage_index = 0;
752 : tensor_storage_index < tensor_type::size(); ++tensor_storage_index) {
753 : // for each coefficient...
754 : for (size_t i = mesh.extents()[0] - 1; i > 0; --i) {
755 : for (size_t j = i; j < mesh.extents()[1]; --j) {
756 : // Check if we are done limiting, and if so we move on to the next
757 : // tensor component.
758 : auto& self_coeffs = get<Tag>(*coeffs_self)[tensor_storage_index];
759 : // We treat the different cases separately to reduce the number of
760 : // times we call minmod, not because it is required for correctness.
761 : if (UNLIKELY(i == j)) {
762 : const double tmp = minmod(i, j, tensor_storage_index);
763 : if (tmp == self_coeffs[mesh.storage_index(Index<VolumeDim>{i, j})]) {
764 : goto next_tensor_index;
765 : }
766 : *limited_any_component |= true;
767 : self_coeffs[mesh.storage_index(Index<VolumeDim>{i, j})] = tmp;
768 : } else {
769 : const double tmp_ij = minmod(i, j, tensor_storage_index);
770 : const double tmp_ji = minmod(j, i, tensor_storage_index);
771 : if (tmp_ij ==
772 : self_coeffs[mesh.storage_index(Index<VolumeDim>{i, j})] and
773 : tmp_ji ==
774 : self_coeffs[mesh.storage_index(Index<VolumeDim>{j, i})]) {
775 : goto next_tensor_index;
776 : }
777 : *limited_any_component |= true;
778 : self_coeffs[mesh.storage_index(Index<VolumeDim>{i, j})] = tmp_ij;
779 : self_coeffs[mesh.storage_index(Index<VolumeDim>{j, i})] = tmp_ji;
780 : }
781 : }
782 : }
783 : next_tensor_index:
784 : continue;
785 : }
786 : return '0';
787 : }
788 :
789 : template <size_t VolumeDim, typename... Tags>
790 : template <typename Tag>
791 0 : char Krivodonova<VolumeDim, tmpl::list<Tags...>>::limit_one_tensor(
792 : const gsl::not_null<Variables<tmpl::list<::Tags::Modal<Tags>...>>*>
793 : coeffs_self,
794 : const gsl::not_null<bool*> limited_any_component, const Mesh<3>& mesh,
795 : const std::unordered_map<DirectionalId<3>, PackagedData,
796 : boost::hash<DirectionalId<3>>>& neighbor_data)
797 : const {
798 : using tensor_type = typename Tag::type;
799 : const auto minmod = [&coeffs_self, &mesh, &neighbor_data, this](
800 : const size_t local_i, const size_t local_j,
801 : const size_t local_k,
802 : const size_t local_tensor_storage_index) {
803 : const auto& self_coeffs =
804 : get<Tag>(*coeffs_self)[local_tensor_storage_index];
805 : double min_abs_coeff = std::abs(self_coeffs[mesh.storage_index(
806 : Index<VolumeDim>{local_i, local_j, local_k})]);
807 : const double sgn_of_coeff = sgn(self_coeffs[mesh.storage_index(
808 : Index<VolumeDim>{local_i, local_j, local_k})]);
809 : for (const auto& kv : neighbor_data) {
810 : const Direction<3>& dir = kv.first.direction();
811 : const auto& neighbor_coeffs =
812 : get<Tag>(kv.second.modal_volume_data)[local_tensor_storage_index];
813 : const size_t index_i =
814 : dir.axis() == Direction<VolumeDim>::Axis::Xi ? local_i - 1 : local_i;
815 : const size_t index_j =
816 : dir.axis() == Direction<VolumeDim>::Axis::Eta ? local_j - 1 : local_j;
817 : const size_t index_k = dir.axis() == Direction<VolumeDim>::Axis::Zeta
818 : ? local_k - 1
819 : : local_k;
820 : // skip neighbors where we cannot compute a finite difference in that
821 : // direction because we are already at the lowest coefficient.
822 : if (index_i == std::numeric_limits<size_t>::max() or
823 : index_j == std::numeric_limits<size_t>::max() or
824 : index_k == std::numeric_limits<size_t>::max()) {
825 : continue;
826 : }
827 : const size_t alpha_index =
828 : dir.axis() == Direction<VolumeDim>::Axis::Xi
829 : ? local_i
830 : : dir.axis() == Direction<VolumeDim>::Axis::Eta ? local_j
831 : : local_k;
832 : const double tmp = dir.sign() * gsl::at(alphas_, alpha_index) *
833 : (neighbor_coeffs[mesh.storage_index(
834 : Index<VolumeDim>{index_i, index_j, index_k})] -
835 : self_coeffs[mesh.storage_index(
836 : Index<VolumeDim>{index_i, index_j, index_k})]);
837 :
838 : min_abs_coeff = std::min(min_abs_coeff, std::abs(tmp));
839 : if (sgn(tmp) != sgn_of_coeff) {
840 : return 0.0;
841 : }
842 : }
843 : return sgn_of_coeff * min_abs_coeff;
844 : };
845 :
846 : for (size_t tensor_storage_index = 0;
847 : tensor_storage_index < tensor_type::size(); ++tensor_storage_index) {
848 : // for each coefficient...
849 : for (size_t i = mesh.extents()[0] - 1; i > 0; --i) {
850 : for (size_t j = i; j < mesh.extents()[1]; --j) {
851 : for (size_t k = j; k < mesh.extents()[2]; --k) {
852 : // Check if we are done limiting, and if so we move on to the next
853 : // tensor component.
854 : auto& self_coeffs = get<Tag>(*coeffs_self)[tensor_storage_index];
855 : // We treat the different cases separately to reduce the number of
856 : // times we call minmod, not because it is required for correctness.
857 : //
858 : // Note that the case `i == k and i != j` cannot be encountered since
859 : // the loop bounds are `i >= j >= k`.
860 : if (UNLIKELY(i == j and i == k)) {
861 : const double tmp = minmod(i, j, k, tensor_storage_index);
862 : if (tmp ==
863 : self_coeffs[mesh.storage_index(Index<VolumeDim>{i, j, k})]) {
864 : goto next_tensor_index;
865 : }
866 : *limited_any_component |= true;
867 : self_coeffs[mesh.storage_index(Index<VolumeDim>{i, j, k})] = tmp;
868 : } else if (i == j and i != k) {
869 : const double tmp_ijk = minmod(i, j, k, tensor_storage_index);
870 : const double tmp_ikj = minmod(i, k, j, tensor_storage_index);
871 : const double tmp_kij = minmod(k, i, j, tensor_storage_index);
872 : if (tmp_ijk == self_coeffs[mesh.storage_index(
873 : Index<VolumeDim>{i, j, k})] and
874 : tmp_ikj == self_coeffs[mesh.storage_index(
875 : Index<VolumeDim>{i, k, j})] and
876 : tmp_kij == self_coeffs[mesh.storage_index(
877 : Index<VolumeDim>{k, i, j})]) {
878 : goto next_tensor_index;
879 : }
880 : *limited_any_component |= true;
881 : self_coeffs[mesh.storage_index(Index<VolumeDim>{i, j, k})] =
882 : tmp_ijk;
883 : self_coeffs[mesh.storage_index(Index<VolumeDim>{i, k, j})] =
884 : tmp_ikj;
885 : self_coeffs[mesh.storage_index(Index<VolumeDim>{k, i, j})] =
886 : tmp_kij;
887 : } else if (i != j and j == k) {
888 : const double tmp_ijk = minmod(i, j, k, tensor_storage_index);
889 : const double tmp_kij = minmod(k, i, j, tensor_storage_index);
890 : const double tmp_kji = minmod(k, j, i, tensor_storage_index);
891 : if (tmp_ijk == self_coeffs[mesh.storage_index(
892 : Index<VolumeDim>{i, j, k})] and
893 : tmp_kij == self_coeffs[mesh.storage_index(
894 : Index<VolumeDim>{k, i, j})] and
895 : tmp_kji == self_coeffs[mesh.storage_index(
896 : Index<VolumeDim>{k, j, i})]) {
897 : goto next_tensor_index;
898 : }
899 : *limited_any_component |= true;
900 : self_coeffs[mesh.storage_index(Index<VolumeDim>{i, j, k})] =
901 : tmp_ijk;
902 : self_coeffs[mesh.storage_index(Index<VolumeDim>{k, i, j})] =
903 : tmp_kij;
904 : self_coeffs[mesh.storage_index(Index<VolumeDim>{k, j, i})] =
905 : tmp_kji;
906 : } else {
907 : const double tmp_ijk = minmod(i, j, k, tensor_storage_index);
908 : const double tmp_jik = minmod(j, i, k, tensor_storage_index);
909 : const double tmp_ikj = minmod(i, k, j, tensor_storage_index);
910 : const double tmp_jki = minmod(j, k, i, tensor_storage_index);
911 : const double tmp_kij = minmod(k, i, j, tensor_storage_index);
912 : const double tmp_kji = minmod(k, j, i, tensor_storage_index);
913 : if (tmp_ijk == self_coeffs[mesh.storage_index(
914 : Index<VolumeDim>{i, j, k})] and
915 : tmp_jik == self_coeffs[mesh.storage_index(
916 : Index<VolumeDim>{j, i, k})] and
917 : tmp_ikj == self_coeffs[mesh.storage_index(
918 : Index<VolumeDim>{i, k, j})] and
919 : tmp_jki == self_coeffs[mesh.storage_index(
920 : Index<VolumeDim>{j, k, i})] and
921 : tmp_kij == self_coeffs[mesh.storage_index(
922 : Index<VolumeDim>{k, i, j})] and
923 : tmp_kji == self_coeffs[mesh.storage_index(
924 : Index<VolumeDim>{k, j, i})]) {
925 : goto next_tensor_index;
926 : }
927 : *limited_any_component |= true;
928 : self_coeffs[mesh.storage_index(Index<VolumeDim>{i, j, k})] =
929 : tmp_ijk;
930 : self_coeffs[mesh.storage_index(Index<VolumeDim>{j, i, k})] =
931 : tmp_jik;
932 : self_coeffs[mesh.storage_index(Index<VolumeDim>{i, k, j})] =
933 : tmp_ikj;
934 : self_coeffs[mesh.storage_index(Index<VolumeDim>{j, k, i})] =
935 : tmp_jki;
936 : self_coeffs[mesh.storage_index(Index<VolumeDim>{k, i, j})] =
937 : tmp_kij;
938 : self_coeffs[mesh.storage_index(Index<VolumeDim>{k, j, i})] =
939 : tmp_kji;
940 : }
941 : }
942 : }
943 : }
944 : next_tensor_index:
945 : continue;
946 : }
947 : return '0';
948 : }
949 :
950 : template <size_t VolumeDim, typename... Tags>
951 : template <typename Tag, size_t Dim>
952 0 : char Krivodonova<VolumeDim, tmpl::list<Tags...>>::
953 : fill_variables_tag_with_spectral_coeffs(
954 : const gsl::not_null<Variables<tmpl::list<::Tags::Modal<Tags>...>>*>
955 : modal_coeffs,
956 : const typename Tag::type& nodal_tensor, const Mesh<Dim>& mesh) const {
957 : auto& coeffs_tensor = get<::Tags::Modal<Tag>>(*modal_coeffs);
958 : auto tensor_it = nodal_tensor.begin();
959 : for (auto coeffs_it = coeffs_tensor.begin(); coeffs_it != coeffs_tensor.end();
960 : (void)++coeffs_it, (void)++tensor_it) {
961 : to_modal_coefficients(make_not_null(&*coeffs_it), *tensor_it, mesh);
962 : }
963 : return '0';
964 : }
965 :
966 : template <size_t VolumeDim, typename... Tags>
967 : void Krivodonova<VolumeDim, tmpl::list<Tags...>>::pup(PUP::er& p) {
968 : p | alphas_;
969 : p | disable_for_debugging_;
970 : }
971 :
972 : template <size_t VolumeDim, typename... Tags>
973 : bool Krivodonova<VolumeDim, tmpl::list<Tags...>>::operator==(
974 : const Krivodonova<VolumeDim, tmpl::list<Tags...>>& rhs) const {
975 : return alphas_ == rhs.alphas_ and
976 : disable_for_debugging_ == rhs.disable_for_debugging_;
977 : }
978 :
979 : template <size_t VolumeDim, typename... Tags>
980 0 : bool operator!=(const Krivodonova<VolumeDim, tmpl::list<Tags...>>& lhs,
981 : const Krivodonova<VolumeDim, tmpl::list<Tags...>>& rhs) {
982 : return not(lhs == rhs);
983 : }
984 : } // namespace Limiters
|
