AoWeno.hpp
1 // Distributed under the MIT License.
2 // See LICENSE.txt for details.
3 
4 #pragma once
5 
6 #include <array>
7 #include <cmath>
8 #include <cstddef>
9 #include <utility>
10 
11 #include "NumericalAlgorithms/FiniteDifference/Reconstruct.hpp"
14 #include "Utilities/Gsl.hpp"
15 #include "Utilities/Math.hpp"
16 #include "Utilities/TMPL.hpp"
17 
18 /// \cond
19 template <size_t Dim>
20 class Direction;
21 template <size_t Dim>
22 class Index;
23 /// \endcond
24 
25 namespace fd::reconstruction {
26 namespace detail {
27 template <size_t NonlinearWeightExponent>
28 struct AoWeno53Reconstructor {
30  const double* const u, const int stride, const double gamma_hi,
31  const double gamma_lo, const double epsilon) noexcept {
32  ASSERT(gamma_hi <= 1.0 and gamma_hi >= 0.0,
33  "gamma_hi must be in [0.0, 1.0] but is " << gamma_hi);
34  ASSERT(gamma_lo <= 1.0 and gamma_lo >= 0.0,
35  "gamma_lo must be in [0.0, 1.0] but is " << gamma_lo);
36  ASSERT(epsilon > 0.0,
37  "epsilon must be greater than zero but is " << epsilon);
38 
39  const std::array moments_sr3_1{
40  1.041666666666666 * u[0] - 0.08333333333333333 * u[-stride] +
41  0.04166666666666666 * u[-2 * stride],
42  0.5 * u[-2 * stride] - 2.0 * u[-stride] + 1.5 * u[0],
43  0.5 * u[-2 * stride] - u[-stride] + 0.5 * u[0]};
44  const std::array moments_sr3_2{0.04166666666666666 * u[stride] +
45  0.9166666666666666 * u[0] +
46  0.04166666666666666 * u[-stride],
47  0.5 * (u[stride] - u[-stride]),
48  0.5 * u[-stride] - u[0] + 0.5 * u[stride]};
49  const std::array moments_sr3_3{
50  0.04166666666666666 * u[2 * stride] - 0.08333333333333333 * u[stride] +
51  1.04166666666666666 * u[0],
52  -1.5 * u[0] + 2.0 * u[stride] - 0.5 * u[2 * stride],
53  0.5 * u[0] - u[stride] + 0.5 * u[2 * stride]};
54  const std::array moments_sr5{-2.95138888888888881e-03 * u[-2 * stride] +
55  5.34722222222222196e-02 * u[-stride] +
56  8.98958333333333304e-01 * u[0] +
57  5.34722222222222196e-02 * u[stride] +
58  -2.95138888888888881e-03 * u[2 * stride],
59  7.08333333333333315e-02 * u[-2 * stride] +
60  -6.41666666666666718e-01 * u[-stride] +
61  6.41666666666666718e-01 * u[stride] +
62  -7.08333333333333315e-02 * u[2 * stride],
63  -3.27380952380952397e-02 * u[-2 * stride] +
64  6.30952380952380931e-01 * u[-stride] +
65  -1.19642857142857140e+00 * u[0] +
66  6.30952380952380931e-01 * u[stride] +
67  -3.27380952380952397e-02 * u[2 * stride],
68  -8.33333333333333287e-02 * u[-2 * stride] +
69  1.66666666666666657e-01 * u[-stride] +
70  -1.66666666666666657e-01 * u[stride] +
71  8.33333333333333287e-02 * u[2 * stride],
72  4.16666666666666644e-02 * u[-2 * stride] +
73  -1.66666666666666657e-01 * u[-stride] +
74  2.50000000000000000e-01 * u[0] +
75  -1.66666666666666657e-01 * u[stride] +
76  4.16666666666666644e-02 * u[2 * stride]};
77 
78  // These are the "first alternative" oscillation indicators
79  constexpr double beta_r3_factor = 37.0 / 3.0;
80  const std::array beta_r3{
81  square(moments_sr3_1[1]) + beta_r3_factor * square(moments_sr3_1[2]),
82  square(moments_sr3_2[1]) + beta_r3_factor * square(moments_sr3_2[2]),
83  square(moments_sr3_3[1]) + beta_r3_factor * square(moments_sr3_3[2])};
84  const double beta_sr5 = square(moments_sr5[1]) +
85  61.0 / 5.0 * moments_sr5[1] * moments_sr5[3] +
86  37.0 / 3.0 * square(moments_sr5[2]) +
87  1538.0 / 7.0 * moments_sr5[2] * moments_sr5[4] +
88  8973.0 / 50.0 * square(moments_sr5[3]) +
89  167158.0 / 49.0 * square(moments_sr5[4]);
90 
91  // Compute linear and normalized nonlinear weights
92  const std::array linear_weights{
93  gamma_hi, 0.5 * (1.0 - gamma_hi) * (1.0 - gamma_lo),
94  (1.0 - gamma_hi) * gamma_lo, 0.5 * (1.0 - gamma_hi) * (1.0 - gamma_lo)};
95  std::array nonlinear_weights{
96  linear_weights[0] / pow<NonlinearWeightExponent>(beta_sr5 + epsilon),
97  linear_weights[1] / pow<NonlinearWeightExponent>(beta_r3[0] + epsilon),
98  linear_weights[2] / pow<NonlinearWeightExponent>(beta_r3[1] + epsilon),
99  linear_weights[3] / pow<NonlinearWeightExponent>(beta_r3[2] + epsilon)};
100  const double normalization = nonlinear_weights[0] + nonlinear_weights[1] +
101  nonlinear_weights[2] + nonlinear_weights[3];
102  for (double& nw : nonlinear_weights) {
103  nw /= normalization;
104  }
105 
106  const std::array<double, 5> moments{
107  {nonlinear_weights[0] / linear_weights[0] *
108  (moments_sr5[0] - linear_weights[1] * moments_sr3_1[0] -
109  linear_weights[2] * moments_sr3_2[0] -
110  linear_weights[3] * moments_sr3_3[0]) +
111  nonlinear_weights[1] * moments_sr3_1[0] +
112  nonlinear_weights[2] * moments_sr3_2[0] +
113  nonlinear_weights[3] * moments_sr3_3[0],
114  nonlinear_weights[0] / linear_weights[0] *
115  (moments_sr5[1] - linear_weights[1] * moments_sr3_1[1] -
116  linear_weights[2] * moments_sr3_2[1] -
117  linear_weights[3] * moments_sr3_3[1]) +
118  nonlinear_weights[1] * moments_sr3_1[1] +
119  nonlinear_weights[2] * moments_sr3_2[1] +
120  nonlinear_weights[3] * moments_sr3_3[1],
121  nonlinear_weights[0] / linear_weights[0] *
122  (moments_sr5[2] - linear_weights[1] * moments_sr3_1[2] -
123  linear_weights[2] * moments_sr3_2[2] -
124  linear_weights[3] * moments_sr3_3[2]) +
125  nonlinear_weights[1] * moments_sr3_1[2] +
126  nonlinear_weights[2] * moments_sr3_2[2] +
127  nonlinear_weights[3] * moments_sr3_3[2],
128  nonlinear_weights[0] / linear_weights[0] * moments_sr5[3],
129  nonlinear_weights[0] / linear_weights[0] * moments_sr5[4]}};
130 
131  // The polynomial values (excluding the L_0(x)) at the lower/upper faces.
132  // The polynomials are (x in [-0.5,0.5]):
133  // L0(x)=1.0 (not in array)
134  // L1(x)=x
135  // L2(x)=x^2-1/12
136  // L3(x)=x^3 - (3/20) x
137  // L4(x)=x^4 - (3/14) x^2 + 3/560
138  const std::array polys_at_plus_half{0.5, 0.16666666666666666, 0.05,
139  0.014285714285714289};
140 
141  return {{moments[0] - polys_at_plus_half[0] * moments[1] +
142  polys_at_plus_half[1] * moments[2] -
143  polys_at_plus_half[2] * moments[3] +
144  polys_at_plus_half[3] * moments[4],
145  moments[0] + polys_at_plus_half[0] * moments[1] +
146  polys_at_plus_half[1] * moments[2] +
147  polys_at_plus_half[2] * moments[3] +
148  polys_at_plus_half[3] * moments[4]}};
149  }
150 
151  SPECTRE_ALWAYS_INLINE static constexpr size_t stencil_width() noexcept {
152  return 5;
153  }
154 };
155 } // namespace detail
156 
157 /*!
158  * \ingroup FiniteDifferenceGroup
159  * \brief Performs an adaptive order (AO) WENO reconstruction using a single 5th
160  * order stencil and a 3rd-order CWENO scheme.
161  *
162  * The AO-WENO(5,3) scheme is based on the scheme presented in
163  * \cite Balsara2016780 but adjusted to do reconstruction on variables instead
164  * of fluxes. The Legendre basis functions on the domain \f$\xi\in[-1/2,1/2]\f$
165  * are given by:
166  *
167  * \f{align*}{
168  * L_0(\xi) &= 1 \\
169  * L_1(\xi) &= \xi \\
170  * L_2(\xi) &= \xi^2-\frac{1}{12} \\
171  * L_3(\xi) &= \xi^3 - \frac{3}{20} \xi \\
172  * L_4(\xi) &= \xi^4 - \frac{3}{14} \xi^2 + \frac{3}{560}
173  * \f}
174  *
175  * The oscillation indicators are given by
176  *
177  * \f{align*}{
178  * \beta_l = \sum_{k=1}^{p}\sum_{m=1}^{p}\sigma_{km} u_{k,l} u_{m,l}
179  * \f}
180  *
181  * where \f$p\f$ is the maximum degree of the basis function used for the
182  * stencil \f$l\f$, and
183  *
184  * \f{align*}{
185  * \sigma_{km}=\sum_{i=1}^{p}\int_{-1/2}^{1/2}
186  * \frac{d^i L_k(\xi)}{d\xi^i}\frac{d^i L_m(\xi)}{d\xi^i}d\xi
187  * \f}
188  *
189  * We write the 3rd-order reconstructed polynomial \f$P^{r3}_l(\xi)\f$
190  * associated with the stencil \f$S^{r3}_l\f$ as
191  *
192  * \f{align*}{
193  * P^{r3}_l(\xi) = u_0 +u_{\xi} L_1(\xi) + u_{\xi 2} L_2(\xi)
194  * \f}
195  *
196  * For the stencil \f$S^{r3}_1\f$ we get
197  *
198  * \f{align*}{
199  * u_0 &= \frac{25}{24} u_{j} -\frac{1}{12} u_{j-1} + \frac{1}{24} u_{j-2} \\
200  * u_{\xi} &= \frac{1}{2}u_{j-2} - 2 u_{j-1} + \frac{3}{2} u_j \\
201  * u_{\xi2} &= \frac{1}{2} u_{j-2} - u_{j-1} + \frac{1}{2} u_j
202  * \f}
203  *
204  * For the stencil \f$S^{r3}_2\f$ we get
205  *
206  * \f{align*}{
207  * u_0 &= \frac{1}{24} u_{j-1} + \frac{11}{12} u_{j} + \frac{1}{24} u_{j+1} \\
208  * u_{\xi} &= \frac{1}{2}(u_{j+1} - u_{j-1}) \\
209  * u_{\xi2} &= \frac{1}{2} u_{j-1} - u_{j} + \frac{1}{2} u_{j+1}
210  * \f}
211  *
212  * For the stencil \f$S^{r3}_3\f$ we get
213  *
214  * \f{align*}{
215  * u_0 &= \frac{25}{24} u_{j} -\frac{1}{12} u_{j+1} + \frac{1}{24} u_{j+2} \\
216  * u_{\xi} &= -\frac{1}{2}u_{j+2} + 2 u_{j+1} - \frac{3}{2} u_j \\
217  * u_{\xi2} &= \frac{1}{2} u_{j+2} - u_{j+1} + \frac{1}{2} u_j
218  * \f}
219  *
220  * The oscillation indicator for the 3rd-order stencils is given by
221  *
222  * \f{align*}{
223  * \beta^{r3}_l = \left(u_{\xi}\right)^2+\frac{13}{3}\left(u_{\xi2}\right)^2
224  * \f}
225  *
226  * We write the 5th-order reconstructed polynomial \f$P^{r5}(\xi)\f$ as
227  *
228  * \f{align*}{
229  * P^{r5}_l(\xi) = u_0 +u_{\xi} L_1(\xi) + u_{\xi 2} L_2(\xi)
230  * + u_{\xi3} L_3(\xi) + u_{\xi4} L_4(\xi)
231  * \f}
232  *
233  * with
234  *
235  * \f{align*}{
236  * u_0 &= -\frac{17}{5760} u_{j-2} + \frac{77}{1440} u_{j-1} +
237  * \frac{863}{960} u_{j} + \frac{77}{1440} u_{j+1} -
238  * \frac{17}{5760} u_{j+2} \\
239  * u_{\xi} &= \frac{17}{240} u_{j-2} - \frac{77}{120} u_{j-1} +
240  * \frac{77}{120} u_{j+1} - \frac{17}{240} u_{j+2} \\
241  * u_{\xi2} &= -\frac{11}{336} u_{j-2} + \frac{53}{84} u_{j-1} -
242  * \frac{67}{56} u_{j} + \frac{53}{84} u_{j+1} -
243  * \frac{11}{336} u_{j+2} \\
244  * u_{\xi3} &= -\frac{1}{12} u_{j-2} + \frac{1}{6} u_{j-1}
245  * -\frac{1}{6} u_{j+1} + \frac{1}{12} u_{j+2} \\
246  * u_{\xi4} &= \frac{1}{24} u_{j-2} -\frac{1}{6} u_{j-1}
247  * +\frac{1}{4} u_{j} - \frac{1}{6}u_{j+1}+\frac{1}{24}u_{j+2}
248  * \f}
249  *
250  * The oscillation indicator is given by
251  *
252  * \f{align*}{
253  * \beta^{r5}&=\left(u_{\xi}+\frac{1}{10}u_{\xi3}\right)^2 \\
254  * &+\frac{13}{3}\left(u_{\xi2}+\frac{123}{455}u_{\xi4}\right)^2\\
255  * &+\frac{781}{20}(u_{\xi3})^2+\frac{1421461}{2275}(u_{\xi4})^2
256  * \f}
257  *
258  * There are two linear weights \f$\gamma_{\mathrm{hi}}\f$ and
259  * \f$\gamma_{\mathrm{lo}}\f$. \f$\gamma_{\mathrm{hi}}\f$ controls how much of
260  * the 5th-order stencil is used in smooth regions, while
261  * \f$\gamma_{\mathrm{lo}}\f$ controls the linear weight of the central
262  * 3rd-order stencil. For larger \f$\gamma_{\mathrm{lo}}\f$, the 3rd-order
263  * method is more centrally weighted. The linear weights for the stencils
264  * are given by
265  *
266  * \f{align*}{
267  * \gamma^{r5}&=\gamma_{\mathrm{hi}} \\
268  * \gamma^{r3}_1& = (1-\gamma_{\mathrm{hi}})(1-\gamma_{\mathrm{lo}})/2 \\
269  * \gamma^{r3}_2& = (1-\gamma_{\mathrm{hi}})\gamma_{\mathrm{lo}} \\
270  * \gamma^{r3}_3& = (1-\gamma_{\mathrm{hi}})(1-\gamma_{\mathrm{lo}})/2
271  * \f}
272  *
273  * We use the standard nonlinear weights instead of the "Z" weights of
274  * \cite Borges20083191
275  *
276  * \f{align*}{
277  * w^{r5}&=\frac{\gamma^{r5}}{(\beta^{r5}+\epsilon)^q} \\
278  * w^{r3}_l&=\frac{\gamma^{r3}_l}{(\beta^{r3}_l+\epsilon)^q}
279  * \f}
280  *
281  * where \f$\epsilon\f$ is a small number used to avoid division by zero. The
282  * normalized nonlinear weights are denoted by \f$\bar{w}^{r5}\f$ and
283  * \f$\bar{w}_l^{r3}\f$. The final reconstructed polynomial \f$P(\xi)\f$ is
284  * given by
285  *
286  * \f{align*}{
287  * P(\xi)&=\frac{\bar{w}^{r5}}{\gamma^{r5}}
288  * \left(P^{r5}(\xi)-\gamma^{r3}_1P^{r3}_1(\xi)
289  * -\gamma^{r3}_2P^{r3}_2(\xi)-\gamma^{r3}_3P^{r3}_3(\xi)\right) \\
290  * &+\bar{w}^{r3}_1P^{r3}_1(\xi)+\bar{w}^{r3}_2P^{r3}_2(\xi)
291  * +\bar{w}^{r3}_3P^{r3}_3(\xi)
292  * \f}
293  *
294  * The weights \f$\gamma_{\mathrm{hi}}\f$ and \f$\gamma_{\mathrm{lo}}\f$
295  * are typically chosen to be in the range \f$[0.85,0.95]\f$.
296  *
297  * ### First alternative oscillation indicators
298  *
299  * Instead of integrating over just the cell, we can instead integrate the basis
300  * functions over the entire fit interval, \f$[-5/2,5/2]\f$. Using this interval
301  * for the \f$S^{r3}_l\f$ and the \f$S^{r5}\f$ stencils we get
302  *
303  * \f{align*}{
304  * \beta^{r3}_l&=(u_{\xi})^2 + \frac{37}{3} (u_{\xi2})^2 \\
305  * \beta^{r5} &=\left(u_{\xi}+\frac{61}{10}u_{\xi3}\right)^2
306  * + \frac{569}{4}u_{\xi3}^2 \\
307  * &+ \frac{1}{8190742}\left(5383 u_{\xi2} + 167158 u_{\xi4}\right)^2
308  * +\frac{4410763}{501474}(u_{\xi2})^2 \\
309  * &=u_{\xi}^2 + \frac{61}{5}u_{\xi}u_{\xi3} + \frac{37}{3}u_{\xi2}^2
310  * + \frac{1538}{7}u_{\xi2}u_{\xi4} \\
311  * &+ \frac{8973}{50}u_{\xi3}^2 + \frac{167158}{49}u_{\xi4}^2
312  * \f}
313  *
314  * Note that the indicator is manifestly non-negative, a required property of
315  * oscillation indicators. These indicators weight high modes more, which means
316  * the scheme is more sensitive to high-frequency features in the solution.
317  *
318  * \note currently it is the alternative indicators that are used. However, an
319  * option to control which are used can readily be added, probably best done
320  * as a template parameter with `if constexpr` to avoid conditionals inside
321  * tight loops.
322  */
323 template <size_t NonlinearWeightExponent, size_t Dim>
326  reconstructed_upper_side_of_face_vars,
328  reconstructed_lower_side_of_face_vars,
329  const gsl::span<const double>& volume_vars,
330  const DirectionMap<Dim, gsl::span<const double>>& ghost_cell_vars,
331  const Index<Dim>& volume_extents, const size_t number_of_variables,
332  const double gamma_hi, const double gamma_lo,
333  const double epsilon) noexcept {
334  detail::reconstruct<detail::AoWeno53Reconstructor<NonlinearWeightExponent>>(
335  reconstructed_upper_side_of_face_vars,
336  reconstructed_lower_side_of_face_vars, volume_vars, ghost_cell_vars,
337  volume_extents, number_of_variables, gamma_hi, gamma_lo, epsilon);
338 }
339 
340 /*!
341  * \brief Returns function pointers to the `aoweno_53` function, lower neighbor
342  * reconstruction, and upper neighbor reconstruction.
343  *
344  * This is useful for controlling template parameters like the
345  * `NonlinearWeightExponent` from an input file by setting a function pointer.
346  * Note that the reason the reconstruction functions instead of say the
347  * `pointwise` member function is returned is to avoid function pointers inside
348  * tight loops.
349  */
350 template <size_t Dim>
351 std::tuple<
356  const Index<Dim>&, size_t, double, double, double) noexcept,
357  void (*)(gsl::not_null<DataVector*>, const DataVector&, const DataVector&,
358  const Index<Dim>&, const Index<Dim>&, const Direction<Dim>&,
359  const double&, const double&, const double&) noexcept,
360  void (*)(gsl::not_null<DataVector*>, const DataVector&, const DataVector&,
361  const Index<Dim>&, const Index<Dim>&, const Direction<Dim>&,
362  const double&, const double&, const double&) noexcept>
363 aoweno_53_function_pointers(size_t nonlinear_weight_exponent) noexcept;
364 } // namespace fd::reconstruction
fd::reconstruction::aoweno_53
void aoweno_53(const gsl::not_null< std::array< gsl::span< double >, Dim > * > reconstructed_upper_side_of_face_vars, const gsl::not_null< std::array< gsl::span< double >, Dim > * > reconstructed_lower_side_of_face_vars, const gsl::span< const double > &volume_vars, const DirectionMap< Dim, gsl::span< const double >> &ghost_cell_vars, const Index< Dim > &volume_extents, const size_t number_of_variables, const double gamma_hi, const double gamma_lo, const double epsilon) noexcept
Performs an adaptive order (AO) WENO reconstruction using a single 5th order stencil and a 3rd-order ...
Definition: AoWeno.hpp:324
utility
fd::reconstruction
Variable and flux vector splitting reconstruction schemes for finite difference methods.
Definition: AoWeno.hpp:25
fd::reconstruction::aoweno_53_function_pointers
std::tuple< void(*)(gsl::not_null< std::array< gsl::span< double >, Dim > * >, gsl::not_null< std::array< gsl::span< double >, Dim > * >, const gsl::span< const double > &, const DirectionMap< Dim, gsl::span< const double >> &, const Index< Dim > &, size_t, double, double, double) noexcept, void(*)(gsl::not_null< DataVector * >, const DataVector &, const DataVector &, const Index< Dim > &, const Index< Dim > &, const Direction< Dim > &, const double &, const double &, const double &) noexcept, void(*)(gsl::not_null< DataVector * >, const DataVector &, const DataVector &, const Index< Dim > &, const Index< Dim > &, const Direction< Dim > &, const double &, const double &, const double &) noexcept > aoweno_53_function_pointers(size_t nonlinear_weight_exponent) noexcept
Returns function pointers to the aoweno_53 function, lower neighbor reconstruction,...
Definition: AoWeno.hpp:363
square
constexpr decltype(auto) square(const T &x)
Compute the square of x
Definition: ConstantExpressions.hpp:55
std::tuple
Index
Definition: Index.hpp:31
cmath
Direction
Definition: Direction.hpp:23
SPECTRE_ALWAYS_INLINE
#define SPECTRE_ALWAYS_INLINE
Definition: ForceInline.hpp:16
cstddef
array
ConstantExpressions.hpp
DataVector
Stores a collection of function values.
Definition: DataVector.hpp:46
ASSERT
#define ASSERT(a, m)
Assert that an expression should be true.
Definition: Assert.hpp:49
gsl::span
Create a span/view on a range, which is cheap to copy (one pointer).
Definition: Gsl.hpp:292
DirectionMap
Definition: DirectionMap.hpp:15
Gsl.hpp
ForceInline.hpp
TMPL.hpp
gsl::not_null
Require a pointer to not be a nullptr
Definition: ReadSpecPiecewisePolynomial.hpp:13