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