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 <optional>
10 : #include <tuple>
11 : #include <utility>
12 :
13 : #include "NumericalAlgorithms/FiniteDifference/FallbackReconstructorType.hpp"
14 : #include "NumericalAlgorithms/FiniteDifference/Reconstruct.hpp"
15 : #include "NumericalAlgorithms/FiniteDifference/Unlimited.hpp"
16 : #include "Utilities/ConstantExpressions.hpp"
17 : #include "Utilities/ForceInline.hpp"
18 : #include "Utilities/Gsl.hpp"
19 : #include "Utilities/Math.hpp"
20 : #include "Utilities/TMPL.hpp"
21 :
22 : /// \cond
23 : class DataVector;
24 : template <size_t Dim>
25 : class Direction;
26 : template <size_t Dim, typename T>
27 : class DirectionMap;
28 : template <size_t Dim>
29 : class Index;
30 : /// \endcond
31 :
32 : namespace fd::reconstruction {
33 : namespace detail {
34 : template <typename LowOrderReconstructor, bool PositivityPreserving,
35 : bool Use9thOrder, bool Use7thOrder>
36 : struct PositivityPreservingAdaptiveOrderReconstructor {
37 : using ReturnType = std::tuple<double, double, std::uint8_t>;
38 : SPECTRE_ALWAYS_INLINE static ReturnType pointwise(
39 : const double* const u, const int stride, const double four_to_the_alpha_5,
40 : // GCC9 complains that six_to_the_alpha_7 and eight_to_the_alpha_9
41 : // are unused because if-constexpr
42 : [[maybe_unused]] const double six_to_the_alpha_7,
43 : [[maybe_unused]] const double eight_to_the_alpha_9) {
44 : using std::get;
45 : if constexpr (Use9thOrder) {
46 : const auto unlimited_9 = UnlimitedReconstructor<8>::pointwise(u, stride);
47 : const ReturnType order_9_result{get<0>(unlimited_9), get<1>(unlimited_9),
48 : 9};
49 :
50 : if (not PositivityPreserving or LIKELY(get<0>(order_9_result) > 0.0 and
51 : get<1>(order_9_result) > 0.0)) {
52 : const double order_9_norm_of_top_modal_coefficient = square(
53 : -1.593380762005595 * u[stride] +
54 : 0.7966903810027975 * u[2 * stride] -
55 : 0.22762582314365648 * u[3 * stride] +
56 : 0.02845322789295706 * u[4 * stride] -
57 : 1.593380762005595 * u[-stride] +
58 : 0.7966903810027975 * u[-2 * stride] -
59 : 0.22762582314365648 * u[-3 * stride] +
60 : 0.02845322789295706 * u[-4 * stride] + 1.991725952506994 * u[0]);
61 :
62 : const double order_9_norm_of_polynomial =
63 : u[stride] * (25.393963433621668 * u[stride] -
64 : 31.738453392103736 * u[2 * stride] +
65 : 14.315575523531798 * u[3 * stride] -
66 : 5.422933317103013 * u[4 * stride] +
67 : 45.309550145164756 * u[-stride] -
68 : 25.682667845756164 * u[-2 * stride] +
69 : 10.394184200706238 * u[-3 * stride] -
70 : 3.5773996341558414 * u[-4 * stride] -
71 : 56.63693768145594 * u[0]) +
72 : u[2 * stride] * (10.664627625179254 * u[2 * stride] -
73 : 9.781510753231265 * u[3 * stride] +
74 : 3.783820939683476 * u[4 * stride] -
75 : 25.682667845756164 * u[-stride] +
76 : 13.59830711617153 * u[-2 * stride] -
77 : 5.064486634342602 * u[-3 * stride] +
78 : 1.5850428636128617 * u[-4 * stride] +
79 : 33.99576779042882 * u[0]) +
80 : u[3 * stride] * (2.5801312593878514 * u[3 * stride] -
81 : 1.812843724346584 * u[4 * stride] +
82 : 10.394184200706238 * u[-stride] -
83 : 5.064486634342602 * u[-2 * stride] +
84 : 1.6716163773782988 * u[-3 * stride] -
85 : 0.4380794296257583 * u[-4 * stride] -
86 : 14.626643302060115 * u[0]) +
87 : u[4 * stride] * (0.5249097623867759 * u[4 * stride] -
88 : 3.5773996341558414 * u[-stride] +
89 : 1.5850428636128617 * u[-2 * stride] -
90 : 0.4380794296257583 * u[-3 * stride] +
91 : 0.07624062080823268 * u[-4 * stride] +
92 : 5.336843456576288 * u[0]) +
93 : u[-stride] * (25.393963433621668 * u[-stride] -
94 : 31.738453392103736 * u[-2 * stride] +
95 : 14.315575523531798 * u[-3 * stride] -
96 : 5.422933317103013 * u[-4 * stride] -
97 : 56.63693768145594 * u[0]) +
98 : u[-2 * stride] * (10.664627625179254 * u[-2 * stride] -
99 : 9.781510753231265 * u[-3 * stride] +
100 : 3.783820939683476 * u[-4 * stride] +
101 : 33.99576779042882 * u[0]) +
102 : u[-3 * stride] * (2.5801312593878514 * u[-3 * stride] -
103 : 1.812843724346584 * u[-4 * stride] -
104 : 14.626643302060115 * u[0]) +
105 : u[-4 * stride] * (0.5249097623867759 * u[-4 * stride] +
106 : 5.336843456576288 * u[0]) +
107 : 33.758463458609164 * square(u[0]);
108 : if (square(eight_to_the_alpha_9) *
109 : order_9_norm_of_top_modal_coefficient <=
110 : order_9_norm_of_polynomial) {
111 : return order_9_result;
112 : }
113 : }
114 : }
115 :
116 : if constexpr (Use7thOrder) {
117 : const auto unlimited_7 = UnlimitedReconstructor<6>::pointwise(u, stride);
118 : const ReturnType order_7_result{get<0>(unlimited_7), get<1>(unlimited_7),
119 : 7};
120 :
121 : if (not PositivityPreserving or LIKELY(get<0>(order_7_result) > 0.0 and
122 : get<1>(order_7_result) > 0.0)) {
123 : const double order_7_norm_of_top_modal_coefficient =
124 : square(0.06936287633138594 * u[-3 * stride] -
125 : 0.4161772579883155 * u[-2 * stride] +
126 : 1.040443144970789 * u[-stride] - //
127 : 1.3872575266277185 * u[0] + //
128 : 1.040443144970789 * u[stride] -
129 : 0.4161772579883155 * u[2 * stride] + //
130 : 0.06936287633138594 * u[3 * stride]);
131 :
132 : const double order_7_norm_of_polynomial =
133 : u[stride] * (3.93094886671763 * u[stride] -
134 : 4.4887583031366605 * u[2 * stride] +
135 : 2.126671427664419 * u[3 * stride] +
136 : 6.081742742499426 * u[-stride] -
137 : 3.1180508323787337 * u[-2 * stride] +
138 : 1.2660604719155235 * u[-3 * stride] -
139 : 8.108990323332568 * u[0]) +
140 : u[2 * stride] * (1.7504056695205172 * u[2 * stride] -
141 : 1.402086588589091 * u[3 * stride] -
142 : 3.1180508323787337 * u[-stride] +
143 : 1.384291080027286 * u[-2 * stride] -
144 : 0.46498946172145633 * u[-3 * stride] +
145 : 4.614303600090953 * u[0]) +
146 : u[3 * stride] * (0.5786954880513824 * u[3 * stride] +
147 : 1.2660604719155235 * u[-stride] -
148 : 0.46498946172145633 * u[-2 * stride] +
149 : 0.10352871936656591 * u[-3 * stride] -
150 : 2.0705743873313183 * u[0]) +
151 : u[-stride] * (3.93094886671763 * u[-stride] -
152 : 4.4887583031366605 * u[-2 * stride] +
153 : 2.126671427664419 * u[-3 * stride] -
154 : 8.108990323332568 * u[0]) +
155 : u[-2 * stride] * (1.7504056695205172 * u[-2 * stride] -
156 : 1.402086588589091 * u[-3 * stride] +
157 : 4.614303600090953 * u[0]) +
158 : u[-3 * stride] * (0.5786954880513824 * u[-3 * stride] -
159 : 2.0705743873313183 * u[0]) +
160 : 5.203166203165525 * square(u[0]);
161 : if (square(six_to_the_alpha_7) *
162 : order_7_norm_of_top_modal_coefficient <=
163 : order_7_norm_of_polynomial) {
164 : return order_7_result;
165 : }
166 : }
167 : }
168 : const auto unlimited_5 = UnlimitedReconstructor<4>::pointwise(u, stride);
169 : const ReturnType order_5_result{get<0>(unlimited_5), get<1>(unlimited_5),
170 : 5};
171 : if (not PositivityPreserving or
172 : LIKELY(get<0>(order_5_result) > 0.0 and get<1>(order_5_result) > 0.0)) {
173 : // The Persson sensor is 4^alpha L2(\hat{u}) <= L2(u)
174 : const double order_5_norm_of_top_modal_coefficient =
175 : 0.2222222222222222 * square(-1.4880952380952381 * u[stride] +
176 : 0.37202380952380953 * u[2 * stride] -
177 : 1.4880952380952381 * u[-stride] +
178 : 0.37202380952380953 * u[-2 * stride] +
179 : 2.232142857142857 * u[0]);
180 :
181 : // Potential optimization: Simple approximation to the integral since we
182 : // only really need about 1 digit of accuracy. This does eliminate aliases
183 : // and so, while reducing the number of FLOPs, might not be accurate
184 : // enough in the cases we're interested in.
185 : //
186 : // const double order_5_norm_of_polynomial =
187 : // 1.1935763888888888 * square(u[-2 * stride]) +
188 : // 0.4340277777777778 * square(u[-stride]) +
189 : // 1.7447916666666667 * square(u[0]) +
190 : // 0.4340277777777778 * square(u[stride]) +
191 : // 1.1935763888888888 * square(u[2 * stride]);
192 : const double order_5_norm_of_polynomial =
193 : (u[stride] * (1.179711612654321 * u[stride] -
194 : 0.963946414792769 * u[2 * stride] +
195 : 1.0904086750440918 * u[-stride] -
196 : 0.5030502507716049 * u[-2 * stride] -
197 : 1.6356130125661377 * u[0]) +
198 : u[2 * stride] *
199 : (0.6699388830329586 * u[2 * stride] -
200 : 0.5030502507716049 * u[-stride] +
201 : 0.154568572944224 * u[-2 * stride] + 0.927411437665344 * u[0]) +
202 : u[-stride] * (1.179711612654321 * u[-stride] -
203 : 0.963946414792769 * u[-2 * stride] -
204 : 1.6356130125661377 * u[0]) +
205 : u[-2 * stride] * (0.6699388830329586 * u[-2 * stride] +
206 : 0.927411437665344 * u[0]) +
207 : 1.4061182415674602 * square(u[0]));
208 : if (square(four_to_the_alpha_5) * order_5_norm_of_top_modal_coefficient <=
209 : order_5_norm_of_polynomial) {
210 : return order_5_result;
211 : }
212 : }
213 : // Drop to low-order reconstructor
214 : const auto low_order = LowOrderReconstructor::pointwise(u, stride);
215 : const ReturnType low_order_result{get<0>(low_order), get<1>(low_order), 2};
216 : if (not PositivityPreserving or LIKELY(get<0>(low_order_result) > 0.0 and
217 : get<1>(low_order_result) > 0.0)) {
218 : return low_order_result;
219 : }
220 : // 1st-order reconstruction to guarantee positivity
221 : return {u[0], u[0], 1};
222 : }
223 :
224 : SPECTRE_ALWAYS_INLINE static constexpr size_t stencil_width() {
225 : return Use9thOrder ? 9 : (Use7thOrder ? 7 : 5);
226 : }
227 : };
228 : } // namespace detail
229 :
230 : /// @{
231 : /*!
232 : * \ingroup FiniteDifferenceGroup
233 : * \brief Performs positivity-preserving adaptive-order FD reconstruction.
234 : *
235 : * Performs a fifth-order unlimited reconstruction. If the reconstructed
236 : * values at the interfaces aren't positive (when `PositivityPreserving` is
237 : * `true`) or when the Persson TCI condition:
238 : *
239 : * \f{align}{
240 : * 4^\alpha \int_{x_{i-5/2}}^{x_{i+5/2}} \hat{u}^2(x) dx
241 : * > \int_{x_{i-5/2}}^{x_{i+5/2}} u^2(x) dx
242 : * \f}
243 : *
244 : * is satisfied, where \f$\hat{u}\f$ is the polynomial with only the largest
245 : * modal coefficient non-zero, then the `LowOrderReconstructor` is used.
246 : *
247 : * If `PositivityPreserving` is `true` then if the low-order reconstructed
248 : * solution isn't positive we use first-order (constant in space)
249 : * reconstruction.
250 : *
251 : * If `Use9thOrder` is `true` then first a ninth-order reconstruction is used,
252 : * followed by fifth-order. If `Use7thOrder` is `true` then seventh-order
253 : * reconstruction is used before fifth-order (but after ninth-order if
254 : * `Use9thOrder` is also `true`). This allows using the highest possible order
255 : * locally for reconstruction.
256 : *
257 : * Fifth order unlimited reconstruction is:
258 : *
259 : * \f{align}{
260 : * \hat{u}_{i+1/2}=\frac{3}{128}u_{i-2} - \frac{5}{32}u_{i-1}
261 : * + \frac{45}{64}u_{i} + \frac{15}{32}u_{i+1} - \frac{5}{128}u_{i+2}
262 : * \f}
263 : *
264 : * Seventh order unlimited reconstruction is:
265 : *
266 : * \f{align}{
267 : * \hat{u}_{i+1/2}&=-\frac{5}{1024}u_{i-3} + \frac{21}{512}u_{i-2}
268 : * - \frac{175}{1024}u_{i-1} + \frac{175}{256}u_{i} \\
269 : * &+ \frac{525}{1024}u_{i+1} - \frac{35}{512}u_{i+2}
270 : * + \frac{7}{1024}u_{i+3}
271 : * \f}
272 : *
273 : * Ninth order unlimited reconstruction is:
274 : *
275 : * \f{align}{
276 : * \hat{u}_{i+1/2}&=\frac{35}{32768}u_{i-4} - \frac{45}{4096}u_{i-3}
277 : * + \frac{441}{8291}u_{i-2} - \frac{735}{4096}u_{i-1} \\
278 : * &+ \frac{11025}{16384}u_{i}
279 : * + \frac{2205}{4096}u_{i+1} - \frac{735}{8192}u_{i+2}
280 : * + \frac{63}{4096}u_{i+3} \\
281 : * &- \frac{45}{32768}u_{i+4}
282 : * \f}
283 : *
284 : */
285 : template <typename LowOrderReconstructor, bool PositivityPreserving,
286 : bool Use9thOrder, bool Use7thOrder, size_t Dim>
287 1 : void positivity_preserving_adaptive_order(
288 : const gsl::not_null<std::array<gsl::span<double>, Dim>*>
289 : reconstructed_upper_side_of_face_vars,
290 : const gsl::not_null<std::array<gsl::span<double>, Dim>*>
291 : reconstructed_lower_side_of_face_vars,
292 : const gsl::span<const double>& volume_vars,
293 : const DirectionMap<Dim, gsl::span<const double>>& ghost_cell_vars,
294 : const Index<Dim>& volume_extents, const size_t number_of_variables,
295 : const double four_to_the_alpha_5, const double six_to_the_alpha_7,
296 : const double eight_to_the_alpha_9) {
297 : detail::reconstruct<detail::PositivityPreservingAdaptiveOrderReconstructor<
298 : LowOrderReconstructor, PositivityPreserving, Use9thOrder, Use7thOrder>>(
299 : reconstructed_upper_side_of_face_vars,
300 : reconstructed_lower_side_of_face_vars, volume_vars, ghost_cell_vars,
301 : volume_extents, number_of_variables, four_to_the_alpha_5,
302 : six_to_the_alpha_7, eight_to_the_alpha_9);
303 : }
304 :
305 : template <typename LowOrderReconstructor, bool PositivityPreserving,
306 : bool Use9thOrder, bool Use7thOrder, size_t Dim>
307 1 : void positivity_preserving_adaptive_order(
308 : const gsl::not_null<std::array<gsl::span<double>, Dim>*>
309 : reconstructed_upper_side_of_face_vars,
310 : const gsl::not_null<std::array<gsl::span<double>, Dim>*>
311 : reconstructed_lower_side_of_face_vars,
312 : const gsl::not_null<
313 : std::optional<std::array<gsl::span<std::uint8_t>, Dim>>*>
314 : reconstruction_order,
315 : const gsl::span<const double>& volume_vars,
316 : const DirectionMap<Dim, gsl::span<const double>>& ghost_cell_vars,
317 : const Index<Dim>& volume_extents, const size_t number_of_variables,
318 : const double four_to_the_alpha_5, const double six_to_the_alpha_7,
319 : const double eight_to_the_alpha_9) {
320 : detail::reconstruct<detail::PositivityPreservingAdaptiveOrderReconstructor<
321 : LowOrderReconstructor, PositivityPreserving, Use9thOrder, Use7thOrder>>(
322 : reconstructed_upper_side_of_face_vars,
323 : reconstructed_lower_side_of_face_vars, reconstruction_order, volume_vars,
324 : ghost_cell_vars, volume_extents, number_of_variables, four_to_the_alpha_5,
325 : six_to_the_alpha_7, eight_to_the_alpha_9);
326 : }
327 : /// @}
328 :
329 : namespace detail {
330 : template <size_t Dim, bool ReturnReconstructionOrder>
331 : using ppao_recons_type = tmpl::conditional_t<
332 : ReturnReconstructionOrder,
333 : void (*)(
334 : gsl::not_null<std::array<gsl::span<double>, Dim>*>,
335 : gsl::not_null<std::array<gsl::span<double>, Dim>*>,
336 : gsl::not_null<std::optional<std::array<gsl::span<std::uint8_t>, Dim>>*>,
337 : const gsl::span<const double>&,
338 : const DirectionMap<Dim, gsl::span<const double>>&, const Index<Dim>&,
339 : size_t, double, double, double),
340 : void (*)(gsl::not_null<std::array<gsl::span<double>, Dim>*>,
341 : gsl::not_null<std::array<gsl::span<double>, Dim>*>,
342 : const gsl::span<const double>&,
343 : const DirectionMap<Dim, gsl::span<const double>>&,
344 : const Index<Dim>&, size_t, double, double, double)>;
345 : }
346 :
347 : /*!
348 : * \brief Returns function pointers to the
349 : * `positivity_preserving_adaptive_order` function, lower neighbor
350 : * reconstruction, and upper neighbor reconstruction.
351 : */
352 : template <size_t Dim, bool ReturnReconstructionOrder>
353 1 : auto positivity_preserving_adaptive_order_function_pointers(
354 : bool positivity_preserving, bool use_9th_order, bool use_7th_order,
355 : FallbackReconstructorType fallback_recons)
356 : -> std::tuple<detail::ppao_recons_type<Dim, ReturnReconstructionOrder>,
357 : void (*)(gsl::not_null<DataVector*>, const DataVector&,
358 : const DataVector&, const Index<Dim>&,
359 : const Index<Dim>&, const Direction<Dim>&,
360 : const double&, const double&, const double&),
361 : void (*)(gsl::not_null<DataVector*>, const DataVector&,
362 : const DataVector&, const Index<Dim>&,
363 : const Index<Dim>&, const Direction<Dim>&,
364 : const double&, const double&, const double&)>;
365 : } // namespace fd::reconstruction
|
