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 <tuple>
10 :
11 : #include "NumericalAlgorithms/FiniteDifference/FallbackReconstructorType.hpp"
12 : #include "NumericalAlgorithms/FiniteDifference/Reconstruct.hpp"
13 : #include "Utilities/ConstantExpressions.hpp"
14 : #include "Utilities/ErrorHandling/Assert.hpp"
15 : #include "Utilities/ForceInline.hpp"
16 : #include "Utilities/Gsl.hpp"
17 :
18 : /// \cond
19 : class DataVector;
20 : template <size_t>
21 : class Direction;
22 : template <size_t Dim, typename T>
23 : class DirectionMap;
24 : template <size_t Dim>
25 : class Index;
26 : /// \endcond
27 :
28 : namespace fd::reconstruction {
29 : namespace detail {
30 : // pointwise reconstruction routine for the original Wcns5z scheme
31 : template <size_t NonlinearWeightExponent>
32 : struct Wcns5zWork {
33 : SPECTRE_ALWAYS_INLINE static std::array<double, 2> pointwise(
34 : const double* const q, const int stride, const double epsilon) {
35 : ASSERT(epsilon > 0.0,
36 : "epsilon must be greater than zero but is " << epsilon);
37 :
38 : using std::abs;
39 :
40 : const std::array beta{
41 : 1.0833333333333333 * square(q[-2 * stride] - 2.0 * q[-stride] + q[0]) +
42 : 0.25 * square(q[-2 * stride] - 4.0 * q[-stride] + 3.0 * q[0]),
43 : 1.0833333333333333 * square(q[-stride] - 2.0 * q[0] + q[stride]) +
44 : 0.25 * square(q[stride] - q[-stride]),
45 : 1.0833333333333333 * square(q[2 * stride] - 2.0 * q[stride] + q[0]) +
46 : 0.25 * square(q[2 * stride] - 4.0 * q[stride] + 3.0 * q[0])};
47 :
48 : const double tau5{abs(beta[2] - beta[0])};
49 :
50 : const std::array epsilon_k{
51 : epsilon * (1.0 + abs(q[0]) + abs(q[-stride]) + abs(q[-2 * stride])),
52 : epsilon * (1.0 + abs(q[0]) + abs(q[-stride]) + abs(q[stride])),
53 : epsilon * (1.0 + abs(q[0]) + abs(q[stride]) + abs(q[2 * stride]))};
54 :
55 : const std::array nw_buffer{
56 : 1.0 + pow<NonlinearWeightExponent>(tau5 / (beta[0] + epsilon_k[0])),
57 : 1.0 + pow<NonlinearWeightExponent>(tau5 / (beta[1] + epsilon_k[1])),
58 : 1.0 + pow<NonlinearWeightExponent>(tau5 / (beta[2] + epsilon_k[2]))};
59 :
60 : // nonlinear weights for upper and lower reconstructions. The factor 1/16
61 : // for `alpha`s is omitted here since it is eventually canceled out by
62 : // denominator when evaluating modified nonlinear weight `omega`s (see the
63 : // documentation of the `wcns5z()` function below).
64 : const std::array alpha_upper{nw_buffer[0], 10.0 * nw_buffer[1],
65 : 5.0 * nw_buffer[2]};
66 : const std::array alpha_lower{nw_buffer[2], 10.0 * nw_buffer[1],
67 : 5.0 * nw_buffer[0]};
68 : const double alpha_norm_upper =
69 : alpha_upper[0] + alpha_upper[1] + alpha_upper[2];
70 : const double alpha_norm_lower =
71 : alpha_lower[0] + alpha_lower[1] + alpha_lower[2];
72 :
73 : // reconstruction stencils
74 : const std::array recons_stencils_upper{
75 : 0.375 * q[-2 * stride] - 1.25 * q[-stride] + 1.875 * q[0],
76 : -0.125 * q[-stride] + 0.75 * q[0] + 0.375 * q[stride],
77 : 0.375 * q[0] + 0.75 * q[stride] - 0.125 * q[2 * stride]};
78 : const std::array recons_stencils_lower{
79 : 0.375 * q[2 * stride] - 1.25 * q[stride] + 1.875 * q[0],
80 : -0.125 * q[stride] + 0.75 * q[0] + 0.375 * q[-stride],
81 : 0.375 * q[0] + 0.75 * q[-stride] - 0.125 * q[-2 * stride]};
82 :
83 : // reconstructed solutions
84 : return {{(alpha_lower[0] * recons_stencils_lower[0] +
85 : alpha_lower[1] * recons_stencils_lower[1] +
86 : alpha_lower[2] * recons_stencils_lower[2]) /
87 : alpha_norm_lower,
88 : (alpha_upper[0] * recons_stencils_upper[0] +
89 : alpha_upper[1] * recons_stencils_upper[1] +
90 : alpha_upper[2] * recons_stencils_upper[2]) /
91 : alpha_norm_upper}};
92 : }
93 : };
94 :
95 : template <size_t NonlinearWeightExponent, class FallbackReconstructor>
96 : struct Wcns5zReconstructor {
97 : SPECTRE_ALWAYS_INLINE static std::array<double, 2> pointwise(
98 : const double* const q, const int stride, const double epsilon,
99 : const size_t max_number_of_extrema) {
100 : // count the number of extrema in the given FD stencil
101 : size_t n_extrema{0};
102 : for (int i = -1; i < 2; ++i) {
103 : // check if q[i * stride] is local maximum
104 : n_extrema += (q[i * stride] > q[(i - 1) * stride]) and
105 : (q[i * stride] > q[(i + 1) * stride]);
106 : // check if q[i * stride] is local minimum
107 : n_extrema += (q[i * stride] < q[(i - 1) * stride]) and
108 : (q[i * stride] < q[(i + 1) * stride]);
109 : }
110 :
111 : // if `n_extrema` is equal or smaller than a specified number, use the
112 : // original Wcns5z reconstruction
113 : if (n_extrema < max_number_of_extrema + 1) {
114 : return Wcns5zWork<NonlinearWeightExponent>::pointwise(q, stride, epsilon);
115 : } else {
116 : // otherwise use a fallback reconstruction method
117 : return FallbackReconstructor::pointwise(q, stride);
118 : }
119 : }
120 :
121 : SPECTRE_ALWAYS_INLINE static constexpr size_t stencil_width() { return 5; }
122 : };
123 :
124 : template <size_t NonlinearWeightExponent>
125 : struct Wcns5zReconstructor<NonlinearWeightExponent, void> {
126 : SPECTRE_ALWAYS_INLINE static std::array<double, 2> pointwise(
127 : const double* const q, const int stride, const double epsilon,
128 : const size_t /*max_number_of_extrema*/) {
129 : return Wcns5zWork<NonlinearWeightExponent>::pointwise(q, stride, epsilon);
130 : }
131 : SPECTRE_ALWAYS_INLINE static constexpr size_t stencil_width() { return 5; }
132 : };
133 :
134 : } // namespace detail
135 :
136 : /*!
137 : * \ingroup FiniteDifferenceGroup
138 : * \brief Performs fifth order weighted compact nonlinear scheme reconstruction
139 : * based on the WENO-Z method (WCNS-5Z). Adaptive fallback combined with
140 : * an auxiliary reconstruction method (e.g. monotonised central) is also
141 : * supported.
142 : *
143 : * Using the WENO oscillation indicators given by \cite Jiang1996
144 : *
145 : * \f{align}
146 : * \beta_0 & = \frac{13}{12} (q_{i-2} - 2q_{i-1} + q_{i})^2
147 : * + \frac{1}{4} (q_{i-2} - 4q_{i-1} + 3q_{i})^2 \\
148 : * \beta_1 & = \frac{13}{12} (q_{i-1} - 2q_{i} + q_{i+1})^2
149 : * + \frac{1}{4} (q_{i+1} - q_{i-1})^2 \\
150 : * \beta_2 & = \frac{13}{12} (q_{i} - 2q_{i+1} + q_{i+2})^2
151 : * + \frac{1}{4} (3q_{i} - 4q_{i+1} + q_{i+2})^2 ,
152 : * \f}
153 : *
154 : * compute the modified nonlinear weights (cf. WENO-Z scheme: see Eq. 42 of
155 : * \cite Borges20083191)
156 : *
157 : * \f{align}
158 : * \omega_k^z = \frac{\alpha_k^z}{\sum_{l=0}^z \alpha_l^z}, \quad
159 : * \alpha_k^z = c_k \left(1 + \left(
160 : * \frac{\tau_5}{\beta_k + \epsilon_k} \right)^q \right),
161 : * \quad k = 0, 1, 2
162 : * \f}
163 : *
164 : * where \f$\epsilon_k\f$ is a factor used to avoid division by zero and is set
165 : * to
166 : *
167 : * \f{align}
168 : * \epsilon_k = \varepsilon (1 + |q_{i+k}| + |q_{i+k-1}| + |q_{i+k-2}|),
169 : * \quad k = 0, 1, 2,
170 : * \f}
171 : *
172 : * linear weights \f$c_k\f$ are adopted from \cite Nonomura2013 (see Table 17
173 : * of it)
174 : *
175 : * \f{align}
176 : * c_0 = 1/16, \quad c_1 = 10/16, \quad c_2 = 5/16,
177 : * \f}
178 : *
179 : * and \f$\tau_5 \equiv |\beta_2-\beta_0|\f$.
180 : *
181 : * Reconstruction stencils are given by Lagrange interpolations (e.g. see Eq.
182 : * 29d - 29f of \cite Brehm2015)
183 : * \f{align}
184 : * q_{i+1/2}^0 &= \frac{3}{8}q_{i-2} -\frac{5}{4}q_{i-1} +\frac{15}{8}q_{i} \\
185 : * q_{i+1/2}^1 &= -\frac{1}{8}q_{i-1} +\frac{3}{4}q_{i} +\frac{3}{8}q_{i+1} \\
186 : * q_{i+1/2}^2 &= \frac{3}{8}q_{i} +\frac{3}{4}q_{i+1} -\frac{1}{8}q_{i+2}
187 : * \f}
188 : *
189 : * and the final reconstructed solution is given by
190 : * \f{align}
191 : * q_{i+1/2} = \sum_{k=0}^2 \omega_k q_{i+1/2}^k .
192 : * \f}
193 : *
194 : * The nonlinear weights exponent \f$q (=1 \text{ or } 2)\f$ and the small
195 : * factor \f$\varepsilon\f$ can be chosen via an input file.
196 : *
197 : * If the template parameter `FallbackReconstructor` is set to one of
198 : * the FD reconstructor structs of which names are listed in
199 : * `fd::reconstruction::FallbackReconstructorType`, adaptive reconstruction
200 : * is performed as follows. For each finite difference stencils, first check how
201 : * many extrema are in the stencil. If the number of local extrema is less than
202 : * or equal to a non-negative integer `max_number_of_extrema` which is given as
203 : * an input parameter, perform the WCNS-5Z reconstruction; otherwise switch to
204 : * the given `FallbackReconstructor` for performing reconstruction. If
205 : * `FallbackReconstructor` is set to `void`, the adaptive method is disabled and
206 : * WCNS-5Z is always used.
207 : *
208 : */
209 : template <size_t NonlinearWeightExponent, class FallbackReconstructor,
210 : size_t Dim>
211 1 : void wcns5z(const gsl::not_null<std::array<gsl::span<double>, Dim>*>
212 : reconstructed_upper_side_of_face_vars,
213 : const gsl::not_null<std::array<gsl::span<double>, Dim>*>
214 : reconstructed_lower_side_of_face_vars,
215 : const gsl::span<const double>& volume_vars,
216 : const DirectionMap<Dim, gsl::span<const double>>& ghost_cell_vars,
217 : const Index<Dim>& volume_extents, const size_t number_of_variables,
218 : const double epsilon, const size_t max_number_of_extrema) {
219 : detail::reconstruct<detail::Wcns5zReconstructor<NonlinearWeightExponent,
220 : FallbackReconstructor>>(
221 : reconstructed_upper_side_of_face_vars,
222 : reconstructed_lower_side_of_face_vars, volume_vars, ghost_cell_vars,
223 : volume_extents, number_of_variables, epsilon, max_number_of_extrema);
224 : }
225 :
226 : /*!
227 : * \brief Returns function pointers to the `wcns5z` function, lower neighbor
228 : * reconstruction, and upper neighbor reconstruction.
229 : *
230 : */
231 : template <size_t Dim>
232 1 : auto wcns5z_function_pointers(size_t nonlinear_weight_exponent,
233 : FallbackReconstructorType fallback_recons)
234 : -> std::tuple<
235 : void (*)(gsl::not_null<std::array<gsl::span<double>, Dim>*>,
236 : gsl::not_null<std::array<gsl::span<double>, Dim>*>,
237 : const gsl::span<const double>&,
238 : const DirectionMap<Dim, gsl::span<const double>>&,
239 : const Index<Dim>&, size_t, double, size_t),
240 : void (*)(gsl::not_null<DataVector*>, const DataVector&,
241 : const DataVector&, const Index<Dim>&, const Index<Dim>&,
242 : const Direction<Dim>&, const double&, const size_t&),
243 : void (*)(gsl::not_null<DataVector*>, const DataVector&,
244 : const DataVector&, const Index<Dim>&, const Index<Dim>&,
245 : const Direction<Dim>&, const double&, const size_t&)>;
246 :
247 : } // namespace fd::reconstruction
|
