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 <cmath>
8 : #include <cstddef>
9 : #include <utility>
10 :
11 : #include "NumericalAlgorithms/FiniteDifference/Reconstruct.hpp"
12 : #include "Utilities/ConstantExpressions.hpp"
13 : #include "Utilities/ForceInline.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 1 : namespace fd::reconstruction {
26 : namespace detail {
27 : template <size_t NonlinearWeightExponent>
28 : struct AoWeno53Reconstructor {
29 : SPECTRE_ALWAYS_INLINE static std::array<double, 2> pointwise(
30 : const double* const u, const int stride, const double gamma_hi,
31 : const double gamma_lo, const double epsilon) {
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() { return 5; }
152 : };
153 : } // namespace detail
154 :
155 : /*!
156 : * \ingroup FiniteDifferenceGroup
157 : * \brief Performs an adaptive order (AO) WENO reconstruction using a single 5th
158 : * order stencil and a 3rd-order CWENO scheme.
159 : *
160 : * The AO-WENO(5,3) scheme is based on the scheme presented in
161 : * \cite Balsara2016780 but adjusted to do reconstruction on variables instead
162 : * of fluxes. The Legendre basis functions on the domain \f$\xi\in[-1/2,1/2]\f$
163 : * are given by:
164 : *
165 : * \f{align*}{
166 : * L_0(\xi) &= 1 \\
167 : * L_1(\xi) &= \xi \\
168 : * L_2(\xi) &= \xi^2-\frac{1}{12} \\
169 : * L_3(\xi) &= \xi^3 - \frac{3}{20} \xi \\
170 : * L_4(\xi) &= \xi^4 - \frac{3}{14} \xi^2 + \frac{3}{560}
171 : * \f}
172 : *
173 : * The oscillation indicators are given by
174 : *
175 : * \f{align*}{
176 : * \beta_l = \sum_{k=1}^{p}\sum_{m=1}^{p}\sigma_{km} u_{k,l} u_{m,l}
177 : * \f}
178 : *
179 : * where \f$p\f$ is the maximum degree of the basis function used for the
180 : * stencil \f$l\f$, and
181 : *
182 : * \f{align*}{
183 : * \sigma_{km}=\sum_{i=1}^{p}\int_{-1/2}^{1/2}
184 : * \frac{d^i L_k(\xi)}{d\xi^i}\frac{d^i L_m(\xi)}{d\xi^i}d\xi
185 : * \f}
186 : *
187 : * We write the 3rd-order reconstructed polynomial \f$P^{r3}_l(\xi)\f$
188 : * associated with the stencil \f$S^{r3}_l\f$ as
189 : *
190 : * \f{align*}{
191 : * P^{r3}_l(\xi) = u_0 +u_{\xi} L_1(\xi) + u_{\xi 2} L_2(\xi)
192 : * \f}
193 : *
194 : * For the stencil \f$S^{r3}_1\f$ we get
195 : *
196 : * \f{align*}{
197 : * u_0 &= \frac{25}{24} u_{j} -\frac{1}{12} u_{j-1} + \frac{1}{24} u_{j-2} \\
198 : * u_{\xi} &= \frac{1}{2}u_{j-2} - 2 u_{j-1} + \frac{3}{2} u_j \\
199 : * u_{\xi2} &= \frac{1}{2} u_{j-2} - u_{j-1} + \frac{1}{2} u_j
200 : * \f}
201 : *
202 : * For the stencil \f$S^{r3}_2\f$ we get
203 : *
204 : * \f{align*}{
205 : * u_0 &= \frac{1}{24} u_{j-1} + \frac{11}{12} u_{j} + \frac{1}{24} u_{j+1} \\
206 : * u_{\xi} &= \frac{1}{2}(u_{j+1} - u_{j-1}) \\
207 : * u_{\xi2} &= \frac{1}{2} u_{j-1} - u_{j} + \frac{1}{2} u_{j+1}
208 : * \f}
209 : *
210 : * For the stencil \f$S^{r3}_3\f$ we get
211 : *
212 : * \f{align*}{
213 : * u_0 &= \frac{25}{24} u_{j} -\frac{1}{12} u_{j+1} + \frac{1}{24} u_{j+2} \\
214 : * u_{\xi} &= -\frac{1}{2}u_{j+2} + 2 u_{j+1} - \frac{3}{2} u_j \\
215 : * u_{\xi2} &= \frac{1}{2} u_{j+2} - u_{j+1} + \frac{1}{2} u_j
216 : * \f}
217 : *
218 : * The oscillation indicator for the 3rd-order stencils is given by
219 : *
220 : * \f{align*}{
221 : * \beta^{r3}_l = \left(u_{\xi}\right)^2+\frac{13}{3}\left(u_{\xi2}\right)^2
222 : * \f}
223 : *
224 : * We write the 5th-order reconstructed polynomial \f$P^{r5}(\xi)\f$ as
225 : *
226 : * \f{align*}{
227 : * P^{r5}_l(\xi) = u_0 +u_{\xi} L_1(\xi) + u_{\xi 2} L_2(\xi)
228 : * + u_{\xi3} L_3(\xi) + u_{\xi4} L_4(\xi)
229 : * \f}
230 : *
231 : * with
232 : *
233 : * \f{align*}{
234 : * u_0 &= -\frac{17}{5760} u_{j-2} + \frac{77}{1440} u_{j-1} +
235 : * \frac{863}{960} u_{j} + \frac{77}{1440} u_{j+1} -
236 : * \frac{17}{5760} u_{j+2} \\
237 : * u_{\xi} &= \frac{17}{240} u_{j-2} - \frac{77}{120} u_{j-1} +
238 : * \frac{77}{120} u_{j+1} - \frac{17}{240} u_{j+2} \\
239 : * u_{\xi2} &= -\frac{11}{336} u_{j-2} + \frac{53}{84} u_{j-1} -
240 : * \frac{67}{56} u_{j} + \frac{53}{84} u_{j+1} -
241 : * \frac{11}{336} u_{j+2} \\
242 : * u_{\xi3} &= -\frac{1}{12} u_{j-2} + \frac{1}{6} u_{j-1}
243 : * -\frac{1}{6} u_{j+1} + \frac{1}{12} u_{j+2} \\
244 : * u_{\xi4} &= \frac{1}{24} u_{j-2} -\frac{1}{6} u_{j-1}
245 : * +\frac{1}{4} u_{j} - \frac{1}{6}u_{j+1}+\frac{1}{24}u_{j+2}
246 : * \f}
247 : *
248 : * The oscillation indicator is given by
249 : *
250 : * \f{align*}{
251 : * \beta^{r5}&=\left(u_{\xi}+\frac{1}{10}u_{\xi3}\right)^2 \\
252 : * &+\frac{13}{3}\left(u_{\xi2}+\frac{123}{455}u_{\xi4}\right)^2\\
253 : * &+\frac{781}{20}(u_{\xi3})^2+\frac{1421461}{2275}(u_{\xi4})^2
254 : * \f}
255 : *
256 : * There are two linear weights \f$\gamma_{\mathrm{hi}}\f$ and
257 : * \f$\gamma_{\mathrm{lo}}\f$. \f$\gamma_{\mathrm{hi}}\f$ controls how much of
258 : * the 5th-order stencil is used in smooth regions, while
259 : * \f$\gamma_{\mathrm{lo}}\f$ controls the linear weight of the central
260 : * 3rd-order stencil. For larger \f$\gamma_{\mathrm{lo}}\f$, the 3rd-order
261 : * method is more centrally weighted. The linear weights for the stencils
262 : * are given by
263 : *
264 : * \f{align*}{
265 : * \gamma^{r5}&=\gamma_{\mathrm{hi}} \\
266 : * \gamma^{r3}_1& = (1-\gamma_{\mathrm{hi}})(1-\gamma_{\mathrm{lo}})/2 \\
267 : * \gamma^{r3}_2& = (1-\gamma_{\mathrm{hi}})\gamma_{\mathrm{lo}} \\
268 : * \gamma^{r3}_3& = (1-\gamma_{\mathrm{hi}})(1-\gamma_{\mathrm{lo}})/2
269 : * \f}
270 : *
271 : * We use the standard nonlinear weights instead of the "Z" weights of
272 : * \cite Borges20083191
273 : *
274 : * \f{align*}{
275 : * w^{r5}&=\frac{\gamma^{r5}}{(\beta^{r5}+\epsilon)^q} \\
276 : * w^{r3}_l&=\frac{\gamma^{r3}_l}{(\beta^{r3}_l+\epsilon)^q}
277 : * \f}
278 : *
279 : * where \f$\epsilon\f$ is a small number used to avoid division by zero. The
280 : * normalized nonlinear weights are denoted by \f$\bar{w}^{r5}\f$ and
281 : * \f$\bar{w}_l^{r3}\f$. The final reconstructed polynomial \f$P(\xi)\f$ is
282 : * given by
283 : *
284 : * \f{align*}{
285 : * P(\xi)&=\frac{\bar{w}^{r5}}{\gamma^{r5}}
286 : * \left(P^{r5}(\xi)-\gamma^{r3}_1P^{r3}_1(\xi)
287 : * -\gamma^{r3}_2P^{r3}_2(\xi)-\gamma^{r3}_3P^{r3}_3(\xi)\right) \\
288 : * &+\bar{w}^{r3}_1P^{r3}_1(\xi)+\bar{w}^{r3}_2P^{r3}_2(\xi)
289 : * +\bar{w}^{r3}_3P^{r3}_3(\xi)
290 : * \f}
291 : *
292 : * The weights \f$\gamma_{\mathrm{hi}}\f$ and \f$\gamma_{\mathrm{lo}}\f$
293 : * are typically chosen to be in the range \f$[0.85,0.95]\f$.
294 : *
295 : * ### First alternative oscillation indicators
296 : *
297 : * Instead of integrating over just the cell, we can instead integrate the basis
298 : * functions over the entire fit interval, \f$[-5/2,5/2]\f$. Using this interval
299 : * for the \f$S^{r3}_l\f$ and the \f$S^{r5}\f$ stencils we get
300 : *
301 : * \f{align*}{
302 : * \beta^{r3}_l&=(u_{\xi})^2 + \frac{37}{3} (u_{\xi2})^2 \\
303 : * \beta^{r5} &=\left(u_{\xi}+\frac{61}{10}u_{\xi3}\right)^2
304 : * + \frac{569}{4}u_{\xi3}^2 \\
305 : * &+ \frac{1}{8190742}\left(5383 u_{\xi2} + 167158 u_{\xi4}\right)^2
306 : * +\frac{4410763}{501474}(u_{\xi2})^2 \\
307 : * &=u_{\xi}^2 + \frac{61}{5}u_{\xi}u_{\xi3} + \frac{37}{3}u_{\xi2}^2
308 : * + \frac{1538}{7}u_{\xi2}u_{\xi4} \\
309 : * &+ \frac{8973}{50}u_{\xi3}^2 + \frac{167158}{49}u_{\xi4}^2
310 : * \f}
311 : *
312 : * Note that the indicator is manifestly non-negative, a required property of
313 : * oscillation indicators. These indicators weight high modes more, which means
314 : * the scheme is more sensitive to high-frequency features in the solution.
315 : *
316 : * \note currently it is the alternative indicators that are used. However, an
317 : * option to control which are used can readily be added, probably best done
318 : * as a template parameter with `if constexpr` to avoid conditionals inside
319 : * tight loops.
320 : */
321 : template <size_t NonlinearWeightExponent, size_t Dim>
322 1 : void aoweno_53(
323 : const gsl::not_null<std::array<gsl::span<double>, Dim>*>
324 : reconstructed_upper_side_of_face_vars,
325 : const gsl::not_null<std::array<gsl::span<double>, Dim>*>
326 : reconstructed_lower_side_of_face_vars,
327 : const gsl::span<const double>& volume_vars,
328 : const DirectionMap<Dim, gsl::span<const double>>& ghost_cell_vars,
329 : const Index<Dim>& volume_extents, const size_t number_of_variables,
330 : const double gamma_hi, const double gamma_lo, const double epsilon) {
331 : detail::reconstruct<detail::AoWeno53Reconstructor<NonlinearWeightExponent>>(
332 : reconstructed_upper_side_of_face_vars,
333 : reconstructed_lower_side_of_face_vars, volume_vars, ghost_cell_vars,
334 : volume_extents, number_of_variables, gamma_hi, gamma_lo, epsilon);
335 : }
336 :
337 : /*!
338 : * \brief Returns function pointers to the `aoweno_53` function, lower neighbor
339 : * reconstruction, and upper neighbor reconstruction.
340 : *
341 : * This is useful for controlling template parameters like the
342 : * `NonlinearWeightExponent` from an input file by setting a function pointer.
343 : * Note that the reason the reconstruction functions instead of say the
344 : * `pointwise` member function is returned is to avoid function pointers inside
345 : * tight loops.
346 : */
347 : template <size_t Dim>
348 1 : auto aoweno_53_function_pointers(size_t nonlinear_weight_exponent)
349 : -> std::tuple<void (*)(gsl::not_null<std::array<gsl::span<double>, Dim>*>,
350 : gsl::not_null<std::array<gsl::span<double>, Dim>*>,
351 : const gsl::span<const double>&,
352 : const DirectionMap<Dim, gsl::span<const double>>&,
353 : const Index<Dim>&, size_t, double, double, double),
354 : void (*)(gsl::not_null<DataVector*>, const DataVector&,
355 : const DataVector&, const Index<Dim>&,
356 : const Index<Dim>&, const Direction<Dim>&,
357 : const double&, const double&, const double&),
358 : void (*)(gsl::not_null<DataVector*>, const DataVector&,
359 : const DataVector&, const Index<Dim>&,
360 : const Index<Dim>&, const Direction<Dim>&,
361 : const double&, const double&, const double&)>;
362 : } // namespace fd::reconstruction
|
