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 constexpr auto lt = dealii::SIMDComparison::less_than;
75 t_r = dealii::compare_and_apply_mask<lt>(
84 (rho_max - rho_U) * denominator,
87 t_r = dealii::compare_and_apply_mask<lt>(
96 (rho_U - rho_min) * denominator,
106 t_r = std::min(t_r, t_max);
107 t_r = std::max(t_r, t_min);
113 const auto rho_new = view.density(U + t_r * P);
114 const auto test_new_min = view.filter_vacuum_density(
115 std::max(Number(0.), rho_new - relax * rho_max));
116 const auto test_new_max = view.filter_vacuum_density(
117 std::max(Number(0.), rho_min - relax * rho_new));
118 if (!(test_new_min == Number(0.) && test_new_max == Number(0.))) {
120 std::cout << std::fixed << std::setprecision(16);
121 std::cout <<
"Bounds violation: high-order density!"
122 <<
"\n\t\trho min: " << rho_min
123 <<
"\n\t\trho min (delta): "
125 <<
"\n\t\trho: " << rho_new
126 <<
"\n\t\trho max (delta): "
128 <<
"\n\t\trho max: " << rho_max <<
"\n"
163 const auto &s_min = std::get<2>(bounds);
165 for (
unsigned int n = 0; n < parameters.newton_max_iterations(); ++n) {
167 const auto U_r = U + t_r * P;
168 const auto rho_r = view.density(U_r);
169 const auto rho_r_gamma =
ryujin::pow(rho_r, gamma);
170 const auto rho_e_r = view.internal_energy(U_r);
173 relax_small * rho_r * rho_e_r - s_min * rho_r * rho_r_gamma;
180 t_l = dealii::compare_and_apply_mask<
181 dealii::SIMDComparison::greater_than>(
182 psi_r, Number(0.), t_r, t_l);
203#ifdef DEBUG_OUTPUT_LIMITER
205 std::cout << std::endl;
206 std::cout << std::fixed << std::setprecision(16);
207 std::cout <<
"t_l: (start) " << t_l << std::endl;
208 std::cout <<
"t_r: (start) " << t_r << std::endl;
212 const auto U_l = U + t_l * P;
213 const auto rho_l = view.density(U_l);
214 const auto rho_l_gamma =
ryujin::pow(rho_l, gamma);
215 const auto rho_e_l = view.internal_energy(U_l);
218 relax_small * rho_l * rho_e_l - s_min * rho_l * rho_l_gamma;
224 const auto lower_bound =
225 (
ScalarNumber(1.) - relax) * s_min * rho_l * rho_l_gamma;
227 !(std::min(Number(0.), psi_l - lower_bound) == Number(0.))) {
229 std::cout << std::fixed << std::setprecision(16);
231 <<
"Bounds violation: low-order specific entropy (critical)!\n";
232 std::cout <<
"\t\tPsi left: 0 <= " << psi_l <<
"\n" << std::endl;
242 t_l = dealii::compare_and_apply_mask<
243 dealii::SIMDComparison::greater_than>(
244 psi_r, Number(0.), t_r, t_l);
251 const Number tolerance(parameters.newton_tolerance());
252 if (std::max(Number(0.), t_r - t_l - tolerance) == Number(0.))
257 const auto drho = view.density(P);
258 const auto drho_e_l = view.internal_energy_derivative(U_l) * P;
259 const auto drho_e_r = view.internal_energy_derivative(U_r) * P;
261 rho_l * drho_e_l + (rho_e_l - gp1 * s_min * rho_l_gamma) * drho;
263 rho_r * drho_e_r + (rho_e_r - gp1 * s_min * rho_r_gamma) * drho;
266 t_l, t_r, psi_l, psi_r, dpsi_l, dpsi_r, Number(-1.));
268#ifdef DEBUG_OUTPUT_LIMITER
269 std::cout <<
"psi_l: " << psi_l << std::endl;
270 std::cout <<
"psi_r: " << psi_r << std::endl;
271 std::cout <<
"dpsi_l: " << dpsi_l << std::endl;
272 std::cout <<
"dpsi_r: " << dpsi_r << std::endl;
273 std::cout <<
"t_l: ( " << n <<
" ) " << t_l << std::endl;
274 std::cout <<
"t_r: ( " << n <<
" ) " << t_r << std::endl;
283 const auto U_new = U + t_l * P;
284 const auto rho_new = view.density(U_new);
285 const auto rho_new_gamma =
ryujin::pow(rho_new, gamma);
286 const auto rho_e_new = view.internal_energy(U_new);
288 auto psi_new = relax_small * rho_new * rho_e_new -
289 s_min * rho_new * rho_new_gamma;
291 const auto lower_bound =
292 (
ScalarNumber(1.) - relax) * s_min * rho_new * rho_new_gamma;
294 const bool e_valid = std::min(Number(0.), rho_e_new) == Number(0.);
295 const bool psi_valid =
296 std::min(Number(0.), psi_new - lower_bound) == Number(0.);
298 if (!e_valid || !psi_valid) {
300 std::cout << std::fixed << std::setprecision(16);
301 std::cout <<
"Bounds violation: high-order specific entropy!\n";
302 std::cout <<
"\t\trho e: 0 <= " << rho_e_new <<
"\n";
303 std::cout <<
"\t\tPsi: 0 <= " << psi_new <<
"\n" << std::endl;
311 return {t_l, success};
typename View::ScalarNumber ScalarNumber
typename View::state_type state_type
std::array< Number, n_bounds > Bounds
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)