14 template <
int dim,
typename Number>
15 std::tuple<Number, bool>
22 const auto view = hyperbolic_system.view<dim, Number>();
27 constexpr ScalarNumber eps = std::numeric_limits<ScalarNumber>::epsilon();
28 const auto small = view.vacuum_state_relaxation_small();
29 const auto large = view.vacuum_state_relaxation_large();
40 const auto &rho_U = view.density(U);
41 const auto &rho_P = view.density(P);
43 const auto &rho_min = std::get<0>(bounds);
44 const auto &rho_max = std::get<1>(bounds);
50 const auto test_min = view.filter_vacuum_density(
51 std::max(Number(0.), rho_U - relax * rho_max));
52 const auto test_max = view.filter_vacuum_density(
53 std::max(Number(0.), rho_min - relax * rho_U));
54 if (!(test_min == Number(0.) && test_max == Number(0.))) {
56 std::cout << std::fixed << std::setprecision(16);
57 std::cout <<
"Bounds violation: low-order density (critical)!"
58 <<
"\n\t\trho min: " << rho_min
59 <<
"\n\t\trho min (delta): "
61 <<
"\n\t\trho: " << rho_U
62 <<
"\n\t\trho max (delta): "
64 <<
"\n\t\trho max: " << rho_max <<
"\n"
70 const Number denominator =
73 t_r = dealii::compare_and_apply_mask<dealii::SIMDComparison::less_than>(
82 (rho_max - rho_U) * denominator,
85 t_r = dealii::compare_and_apply_mask<dealii::SIMDComparison::less_than>(
94 (rho_U - rho_min) * denominator,
104 t_r = std::min(t_r, t_max);
105 t_r = std::max(t_r, t_min);
111 const auto rho_new = view.density(U + t_r * P);
112 const auto test_new_min = view.filter_vacuum_density(
113 std::max(Number(0.), rho_new - relax * rho_max));
114 const auto test_new_max = view.filter_vacuum_density(
115 std::max(Number(0.), rho_min - relax * rho_new));
116 if (!(test_new_min == Number(0.) && test_new_max == Number(0.))) {
118 std::cout << std::fixed << std::setprecision(16);
119 std::cout <<
"Bounds violation: high-order density!"
120 <<
"\n\t\trho min: " << rho_min
121 <<
"\n\t\trho min (delta): "
123 <<
"\n\t\trho: " << rho_new
124 <<
"\n\t\trho max (delta): "
126 <<
"\n\t\trho max: " << rho_max <<
"\n"
142 const auto &gamma = std::get<3>(bounds) ;
143 const Number gm1 = gamma - Number(1.);
145 const auto interpolation_b = view.eos_interpolation_b();
165 const auto &s_min = std::get<2>(bounds);
167 for (
unsigned int n = 0; n < parameters.newton_max_iterations(); ++n) {
169 const auto U_r = U + t_r * P;
170 const auto rho_r = view.density(U_r);
171 const auto rho_r_gamma =
ryujin::pow(rho_r, gamma);
172 const auto rho_e_r = view.internal_energy(U_r);
173 const auto covolume_r = Number(1.) - interpolation_b * rho_r;
176 relax_small * rho_r * rho_e_r -
177 s_min * rho_r * rho_r_gamma *
ryujin::pow(covolume_r, -gm1);
184 t_l = dealii::compare_and_apply_mask<
185 dealii::SIMDComparison::greater_than>(
186 psi_r, Number(0.), t_r, t_l);
207#ifdef DEBUG_OUTPUT_LIMITER
209 std::cout << std::endl;
210 std::cout << std::fixed << std::setprecision(16);
211 std::cout <<
"t_l: (start) " << t_l << std::endl;
212 std::cout <<
"t_r: (start) " << t_r << std::endl;
216 const auto U_l = U + t_l * P;
217 const auto rho_l = view.density(U_l);
218 const auto rho_l_gamma =
ryujin::pow(rho_l, gamma);
219 const auto rho_e_l = view.internal_energy(U_l);
220 const auto covolume_l = Number(1.) - interpolation_b * rho_l;
223 relax_small * rho_l * rho_e_l -
224 s_min * rho_l * rho_l_gamma *
ryujin::pow(covolume_l, -gm1);
230 const auto lower_bound = (
ScalarNumber(1.) - relax) * s_min * rho_l *
233 !(std::min(Number(0.), psi_l - lower_bound) == Number(0.))) {
235 std::cout << std::fixed << std::setprecision(16);
237 <<
"Bounds violation: low-order specific entropy (critical)!\n";
238 std::cout <<
"\t\tPsi left: 0 <= " << psi_l <<
"\n" << std::endl;
248 t_l = dealii::compare_and_apply_mask<
249 dealii::SIMDComparison::greater_than>(
250 psi_r, Number(0.), t_r, t_l);
257 const Number tolerance(parameters.newton_tolerance());
258 if (std::max(Number(0.), t_r - t_l - tolerance) == Number(0.))
263 const auto drho = view.density(P);
264 const auto drho_e_l = view.internal_energy_derivative(U_l) * P;
265 const auto drho_e_r = view.internal_energy_derivative(U_r) * P;
267 const auto extra_term_l =
269 (covolume_l + gamma - interpolation_b * rho_l);
270 const auto extra_term_r =
272 (covolume_r + gamma - interpolation_b * rho_r);
275 rho_l * drho_e_l + (rho_e_l - extra_term_l) * drho;
277 rho_r * drho_e_r + (rho_e_r - extra_term_r) * drho;
280 t_l, t_r, psi_l, psi_r, dpsi_l, dpsi_r, Number(-1.));
282#ifdef DEBUG_OUTPUT_LIMITER
283 std::cout <<
"psi_l: " << psi_l << std::endl;
284 std::cout <<
"psi_r: " << psi_r << std::endl;
285 std::cout <<
"dpsi_l: " << dpsi_l << std::endl;
286 std::cout <<
"dpsi_r: " << dpsi_r << std::endl;
287 std::cout <<
"t_l: ( " << n <<
" ) " << t_l << std::endl;
288 std::cout <<
"t_r: ( " << n <<
" ) " << t_r << std::endl;
297 const auto U_new = U + t_l * P;
299 const auto rho_new = view.density(U_new);
300 const auto rho_new_gamma =
ryujin::pow(rho_new, gamma);
301 const auto rho_e_new = view.internal_energy(U_new);
302 const auto covolume_new = Number(1.) - interpolation_b * rho_new;
305 relax_small * rho_new * rho_e_new -
306 s_min * rho_new * rho_new_gamma *
ryujin::pow(covolume_new, -gm1);
308 const auto lower_bound = (
ScalarNumber(1.) - relax) * s_min *
309 rho_new * rho_new_gamma *
312 const bool e_valid = std::min(Number(0.), rho_e_new) == Number(0.);
313 const bool psi_valid =
314 std::min(Number(0.), psi_new - lower_bound) == Number(0.);
316 if (!e_valid || !psi_valid) {
318 std::cout << std::fixed << std::setprecision(16);
319 std::cout <<
"Bounds violation: high-order specific entropy!\n";
320 std::cout <<
"\t\trho e: 0 <= " << rho_e_new <<
"\n";
321 std::cout <<
"\t\tPsi: 0 <= " << psi_new <<
"\n" << std::endl;
329 return {t_l, success};
typename View::state_type state_type
typename View::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.))
std::array< Number, n_bounds > Bounds
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)