ryujin 2.1.1 revision 9dcb748690310d6a540ebb8b066d1a0834fc7604
limiter.template.h
Go to the documentation of this file.
1//
2// SPDX-License-Identifier: MIT
3// Copyright (C) 2020 - 2023 by the ryujin authors
4//
5
6#pragma once
7
8#include "limiter.h"
9
10namespace ryujin
11{
12 namespace Euler
13 {
14 template <int dim, typename Number>
15 std::tuple<Number, bool>
17 const state_type &U,
18 const state_type &P,
19 const Number t_min /* = Number(0.) */,
20 const Number t_max /* = Number(1.) */)
21 {
22 bool success = true;
23 Number t_r = t_max;
24
25 constexpr ScalarNumber eps = std::numeric_limits<ScalarNumber>::epsilon();
26 const ScalarNumber relax_small = ScalarNumber(1. + 10. * eps);
27 const ScalarNumber relax =
28 ScalarNumber(1. + hyperbolic_system.vacuum_state_relaxation() * eps);
29
30 /*
31 * First limit the density rho.
32 *
33 * See [Guermond, Nazarov, Popov, Thomas] (4.8):
34 */
35
36 {
37 const auto &rho_U = hyperbolic_system.density(U);
38 const auto &rho_P = hyperbolic_system.density(P);
39
40 const auto &rho_min = std::get<0>(bounds);
41 const auto &rho_max = std::get<1>(bounds);
42
43 /*
44 * Verify that rho_U is within bounds. This property might be
45 * violated for relative CFL numbers larger than 1.
46 */
47 const auto test_min = hyperbolic_system.filter_vacuum_density(
48 std::max(Number(0.), rho_U - relax * rho_max));
49 const auto test_max = hyperbolic_system.filter_vacuum_density(
50 std::max(Number(0.), rho_min - relax * rho_U));
51 if (!(test_min == Number(0.) && test_max == Number(0.))) {
52#ifdef DEBUG_OUTPUT
53 std::cout << std::fixed << std::setprecision(16);
54 std::cout << "Bounds violation: low-order density (critical)!"
55 << "\n\t\trho min: " << rho_min
56 << "\n\t\trho min (delta): "
57 << negative_part(rho_U - rho_min)
58 << "\n\t\trho: " << rho_U
59 << "\n\t\trho max (delta): "
60 << positive_part(rho_U - rho_max)
61 << "\n\t\trho max: " << rho_max << "\n"
62 << std::endl;
63#endif
64 success = false;
65 }
66
67 const Number denominator =
68 ScalarNumber(1.) / (std::abs(rho_P) + eps * rho_max);
69
70 t_r = dealii::compare_and_apply_mask<dealii::SIMDComparison::less_than>(
71 rho_max,
72 rho_U + t_r * rho_P,
73 /*
74 * rho_P is positive.
75 *
76 * Note: Do not take an absolute value here. If we are out of
77 * bounds we have to ensure that t_r is set to t_min.
78 */
79 (rho_max - rho_U) * denominator,
80 t_r);
81
82 t_r = dealii::compare_and_apply_mask<dealii::SIMDComparison::less_than>(
83 rho_U + t_r * rho_P,
84 rho_min,
85 /*
86 * rho_P is negative.
87 *
88 * Note: Do not take an absolute value here. If we are out of
89 * bounds we have to ensure that t_r is set to t_min.
90 */
91 (rho_U - rho_min) * denominator,
92 t_r);
93
94 /*
95 * Ensure that t_min <= t <= t_max. This might not be the case if
96 * rho_U is outside the interval [rho_min, rho_max]. Furthermore,
97 * the quotient we take above is prone to numerical cancellation in
98 * particular in the second pass of the limiter when rho_P might be
99 * small.
100 */
101 t_r = std::min(t_r, t_max);
102 t_r = std::max(t_r, t_min);
103
104#ifdef CHECK_BOUNDS
105 /*
106 * Verify that the new state is within bounds:
107 */
108 const auto rho_new = hyperbolic_system.density(U + t_r * P);
109 const auto test_new_min = hyperbolic_system.filter_vacuum_density(
110 std::max(Number(0.), rho_new - relax * rho_max));
111 const auto test_new_max = hyperbolic_system.filter_vacuum_density(
112 std::max(Number(0.), rho_min - relax * rho_new));
113 if (!(test_new_min == Number(0.) && test_new_max == Number(0.))) {
114#ifdef DEBUG_OUTPUT
115 std::cout << std::fixed << std::setprecision(16);
116 std::cout << "Bounds violation: high-order density!"
117 << "\n\t\trho min: " << rho_min
118 << "\n\t\trho min (delta): "
119 << negative_part(rho_new - rho_min)
120 << "\n\t\trho: " << rho_new
121 << "\n\t\trho max (delta): "
122 << positive_part(rho_new - rho_max)
123 << "\n\t\trho max: " << rho_max << "\n"
124 << std::endl;
125#endif
126 success = false;
127 }
128#endif
129 }
130
131 /*
132 * Then limit the specific entropy:
133 *
134 * See [Guermond, Nazarov, Popov, Thomas], Section 4.6 + Section 5.1:
135 */
136
137 Number t_l = t_min; // good state
138
139 const ScalarNumber gamma = hyperbolic_system.gamma();
140 const ScalarNumber gp1 = gamma + ScalarNumber(1.);
141
142 {
143 /*
144 * Prepare a quadratic Newton method:
145 *
146 * Given initial limiter values t_l and t_r with psi(t_l) > 0 and
147 * psi(t_r) < 0 we try to find t^\ast with psi(t^\ast) \approx 0.
148 *
149 * Here, psi is a 3-convex function obtained by scaling the specific
150 * entropy s:
151 *
152 * psi = \rho ^ {\gamma + 1} s
153 *
154 * (s in turn was defined as s =\varepsilon \rho ^{-\gamma}, where
155 * \varepsilon = (\rho e) is the internal energy.)
156 */
157
158 const auto &s_min = std::get<2>(bounds);
159
160 for (unsigned int n = 0; n < newton_max_iter; ++n) {
161
162 const auto U_r = U + t_r * P;
163 const auto rho_r = hyperbolic_system.density(U_r);
164 const auto rho_r_gamma = ryujin::pow(rho_r, gamma);
165 const auto rho_e_r = hyperbolic_system.internal_energy(U_r);
166
167 auto psi_r =
168 relax_small * rho_r * rho_e_r - s_min * rho_r * rho_r_gamma;
169
170#ifndef CHECK_BOUNDS
171 /*
172 * If psi_r > 0 the right state is fine, force returning t_r by
173 * setting t_l = t_r:
174 */
175 t_l = dealii::compare_and_apply_mask<
176 dealii::SIMDComparison::greater_than>(
177 psi_r, Number(0.), t_r, t_l);
178
179 /*
180 * If we have set t_l = t_r everywhere then all states state U_r
181 * with t_r obey the specific entropy inequality and we can
182 * break.
183 *
184 * This is a very important optimization: Only for 1 in (25 to
185 * 50) cases do we actually need to limit on the specific entropy
186 * because one of the right states failed. So we can skip
187 * constructing the left state U_l, which is expensive.
188 *
189 * This implies unfortunately that we might not accurately report
190 * whether the low_order update U itself obeyed bounds because
191 * U_r = U + t_r * P pushed us back into bounds. We thus skip
192 * this shortcut if `CHECK_BOUNDS` is set.
193 */
194 if (t_l == t_r)
195 break;
196#endif
197
198#ifdef DEBUG_OUTPUT_LIMITER
199 if (n == 0) {
200 std::cout << std::endl;
201 std::cout << std::fixed << std::setprecision(16);
202 std::cout << "t_l: (start) " << t_l << std::endl;
203 std::cout << "t_r: (start) " << t_r << std::endl;
204 }
205#endif
206
207 const auto U_l = U + t_l * P;
208 const auto rho_l = hyperbolic_system.density(U_l);
209 const auto rho_l_gamma = ryujin::pow(rho_l, gamma);
210 const auto rho_e_l = hyperbolic_system.internal_energy(U_l);
211
212 auto psi_l =
213 relax_small * rho_l * rho_e_l - s_min * rho_l * rho_l_gamma;
214
215 /*
216 * Verify that the left state is within bounds. This property might
217 * be violated for relative CFL numbers larger than 1.
218 */
219 const auto lower_bound =
220 (ScalarNumber(1.) - relax) * s_min * rho_l * rho_l_gamma;
221 if (n == 0 &&
222 !(std::min(Number(0.), psi_l - lower_bound) == Number(0.))) {
223#ifdef DEBUG_OUTPUT
224 std::cout << std::fixed << std::setprecision(16);
225 std::cout
226 << "Bounds violation: low-order specific entropy (critical)!\n";
227 std::cout << "\t\tPsi left: 0 <= " << psi_l << "\n" << std::endl;
228#endif
229 success = false;
230 }
231
232#ifdef CHECK_BOUNDS
233 /*
234 * If psi_r > 0 the right state is fine, force returning t_r by
235 * setting t_l = t_r:
236 */
237 t_l = dealii::compare_and_apply_mask<
238 dealii::SIMDComparison::greater_than>(
239 psi_r, Number(0.), t_r, t_l);
240#endif
241
242 /*
243 * Break if the window between t_l and t_r is within the prescribed
244 * tolerance:
245 */
246 if (std::max(Number(0.), t_r - t_l - newton_tolerance) == Number(0.))
247 break;
248
249 /* We got unlucky and have to perform a Newton step: */
250
251 const auto drho = hyperbolic_system.density(P);
252 const auto drho_e_l =
253 hyperbolic_system.internal_energy_derivative(U_l) * P;
254 const auto drho_e_r =
255 hyperbolic_system.internal_energy_derivative(U_r) * P;
256 const auto dpsi_l =
257 rho_l * drho_e_l + (rho_e_l - gp1 * s_min * rho_l_gamma) * drho;
258 const auto dpsi_r =
259 rho_r * drho_e_r + (rho_e_r - gp1 * s_min * rho_r_gamma) * drho;
260
262 t_l, t_r, psi_l, psi_r, dpsi_l, dpsi_r, Number(-1.));
263
264#ifdef DEBUG_OUTPUT_LIMITER
265 std::cout << "psi_l: " << psi_l << std::endl;
266 std::cout << "psi_r: " << psi_r << std::endl;
267 std::cout << "dpsi_l: " << dpsi_l << std::endl;
268 std::cout << "dpsi_r: " << dpsi_r << std::endl;
269 std::cout << "t_l: ( " << n << " ) " << t_l << std::endl;
270 std::cout << "t_r: ( " << n << " ) " << t_r << std::endl;
271#endif
272 }
273
274#ifdef CHECK_BOUNDS
275 /*
276 * Verify that the new state is within bounds:
277 */
278 {
279 const auto U_new = U + t_l * P;
280 const auto rho_new = hyperbolic_system.density(U_new);
281 const auto rho_new_gamma = ryujin::pow(rho_new, gamma);
282 const auto rho_e_new = hyperbolic_system.internal_energy(U_new);
283
284 auto psi_new = relax_small * rho_new * rho_e_new -
285 s_min * rho_new * rho_new_gamma;
286
287 const auto lower_bound =
288 (ScalarNumber(1.) - relax) * s_min * rho_new * rho_new_gamma;
289
290 const bool e_valid = std::min(Number(0.), rho_e_new) == Number(0.);
291 const bool psi_valid =
292 std::min(Number(0.), psi_new - lower_bound) == Number(0.);
293
294 if (!e_valid || !psi_valid) {
295#ifdef DEBUG_OUTPUT
296 std::cout << std::fixed << std::setprecision(16);
297 std::cout << "Bounds violation: high-order specific entropy!\n";
298 std::cout << "\t\trho e: 0 <= " << rho_e_new << "\n";
299 std::cout << "\t\tPsi: 0 <= " << psi_new << "\n" << std::endl;
300#endif
301 success = false;
302 }
303 }
304#endif
305 }
306
307 return {t_l, success};
308 }
309
310 } // namespace Euler
311} // namespace ryujin
typename HyperbolicSystemView::state_type state_type
Definition: limiter.h:67
std::array< Number, n_bounds > Bounds
Definition: limiter.h:119
typename HyperbolicSystemView::ScalarNumber ScalarNumber
Definition: limiter.h:90
std::tuple< Number, bool > limit(const Bounds &bounds, const state_type &U, const state_type &P, const Number t_min=Number(0.), const Number t_max=Number(1.))
DEAL_II_ALWAYS_INLINE void quadratic_newton_step(Number &p_1, Number &p_2, const Number phi_p_1, const Number phi_p_2, const Number dphi_p_1, const Number dphi_p_2, const Number sign=Number(1.0))
Definition: newton.h:39
T pow(const T x, const T b)
DEAL_II_ALWAYS_INLINE Number negative_part(const Number number)
Definition: simd.h:123
DEAL_II_ALWAYS_INLINE Number positive_part(const Number number)
Definition: simd.h:111