12 namespace ScalarConservation
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();
30 const auto &u_U = view.state(U);
31 const auto &u_P = view.state(P);
33 const auto &u_min = std::get<0>(bounds);
34 const auto &u_max = std::get<1>(bounds);
40 const auto test_min = std::max(Number(0.), u_U - relax * u_max);
41 const auto test_max = std::max(Number(0.), u_min - relax * u_U);
42 if (!(test_min == Number(0.) && test_max == Number(0.))) {
44 std::cout << std::fixed << std::setprecision(16);
45 std::cout <<
"Bounds violation: low-order state (critical)!"
46 <<
"\n\t\tu min: " << u_min
50 <<
"\n\t\tu max: " << u_max <<
"\n"
56 const auto regularization =
57 Number(100. * std::numeric_limits<ScalarNumber>::min());
59 const Number denominator =
61 std::max(regularization, std::abs(u_P) + eps * u_max);
63 t_r = dealii::compare_and_apply_mask<dealii::SIMDComparison::less_than>(
72 (u_max - u_U) * denominator,
75 t_r = dealii::compare_and_apply_mask<dealii::SIMDComparison::less_than>(
84 (u_U - u_min) * denominator,
94 t_r = std::min(t_r, t_max);
95 t_r = std::max(t_r, t_min);
101 const auto u_new = view.state(U + t_r * P);
102 const auto test_new_min = std::max(Number(0.), u_new - relax * u_max);
103 const auto test_new_max = std::max(Number(0.), u_min - relax * u_new);
104 if (!(test_new_min == Number(0.) && test_new_max == Number(0.))) {
106 std::cout << std::fixed << std::setprecision(16);
107 std::cout <<
"Bounds violation: high-order state!"
108 <<
"\n\t\tu min: " << u_min
110 <<
"\n\t\tu: " << u_new
112 <<
"\n\t\tu max: " << u_max <<
"\n"
119 return {t_r, success};
std::array< Number, n_bounds > Bounds
typename View::state_type state_type
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.))
typename View::ScalarNumber ScalarNumber
DEAL_II_ALWAYS_INLINE Number negative_part(const Number number)
DEAL_II_ALWAYS_INLINE Number positive_part(const Number number)