14 template <
int dim,
typename Number>
15 std::tuple<Number, bool>
25 constexpr ScalarNumber eps = std::numeric_limits<ScalarNumber>::epsilon();
28 ScalarNumber(1. + hyperbolic_system.vacuum_state_relaxation() * eps);
37 const auto &rho_U = hyperbolic_system.density(U);
38 const auto &rho_P = hyperbolic_system.density(P);
40 const auto &rho_min = std::get<0>(bounds);
41 const auto &rho_max = std::get<1>(bounds);
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.))) {
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): "
58 <<
"\n\t\trho: " << rho_U
59 <<
"\n\t\trho max (delta): "
61 <<
"\n\t\trho max: " << rho_max <<
"\n"
67 const Number denominator =
70 t_r = dealii::compare_and_apply_mask<dealii::SIMDComparison::less_than>(
79 (rho_max - rho_U) * denominator,
82 t_r = dealii::compare_and_apply_mask<dealii::SIMDComparison::less_than>(
91 (rho_U - rho_min) * denominator,
101 t_r = std::min(t_r, t_max);
102 t_r = std::max(t_r, t_min);
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.))) {
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): "
120 <<
"\n\t\trho: " << rho_new
121 <<
"\n\t\trho max (delta): "
123 <<
"\n\t\trho max: " << rho_max <<
"\n"
158 const auto &s_min = std::get<2>(bounds);
160 for (
unsigned int n = 0; n < newton_max_iter; ++n) {
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);
168 relax_small * rho_r * rho_e_r - s_min * rho_r * rho_r_gamma;
175 t_l = dealii::compare_and_apply_mask<
176 dealii::SIMDComparison::greater_than>(
177 psi_r, Number(0.), t_r, t_l);
198#ifdef DEBUG_OUTPUT_LIMITER
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;
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);
213 relax_small * rho_l * rho_e_l - s_min * rho_l * rho_l_gamma;
219 const auto lower_bound =
220 (
ScalarNumber(1.) - relax) * s_min * rho_l * rho_l_gamma;
222 !(std::min(Number(0.), psi_l - lower_bound) == Number(0.))) {
224 std::cout << std::fixed << std::setprecision(16);
226 <<
"Bounds violation: low-order specific entropy (critical)!\n";
227 std::cout <<
"\t\tPsi left: 0 <= " << psi_l <<
"\n" << std::endl;
237 t_l = dealii::compare_and_apply_mask<
238 dealii::SIMDComparison::greater_than>(
239 psi_r, Number(0.), t_r, t_l);
246 if (std::max(Number(0.), t_r - t_l - newton_tolerance) == Number(0.))
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;
257 rho_l * drho_e_l + (rho_e_l - gp1 * s_min * rho_l_gamma) * drho;
259 rho_r * drho_e_r + (rho_e_r - gp1 * s_min * rho_r_gamma) * drho;
262 t_l, t_r, psi_l, psi_r, dpsi_l, dpsi_r, Number(-1.));
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;
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);
284 auto psi_new = relax_small * rho_new * rho_e_new -
285 s_min * rho_new * rho_new_gamma;
287 const auto lower_bound =
288 (
ScalarNumber(1.) - relax) * s_min * rho_new * rho_new_gamma;
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.);
294 if (!e_valid || !psi_valid) {
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;
307 return {t_l, success};
typename HyperbolicSystemView::state_type state_type
std::array< Number, n_bounds > Bounds
typename HyperbolicSystemView::ScalarNumber ScalarNumber
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))
T pow(const T x, const T b)
DEAL_II_ALWAYS_INLINE Number negative_part(const Number number)
DEAL_II_ALWAYS_INLINE Number positive_part(const Number number)