12 namespace ShallowWater
14 template <
int dim,
typename Number>
15 std::tuple<Number, bool>
28 constexpr ScalarNumber eps = std::numeric_limits<ScalarNumber>::epsilon();
30 1. + hyperbolic_system.dry_state_relaxation_sharp() * eps);
32 ScalarNumber(1. + hyperbolic_system.dry_state_relaxation_mollified() *
41 auto h_U = hyperbolic_system.water_depth(U);
42 const auto &h_P = hyperbolic_system.water_depth(P);
43 const auto &h_min = std::get<0>(bounds);
44 const auto &h_max = std::get<1>(bounds);
46 const auto test_min = hyperbolic_system.filter_dry_water_depth(
47 std::max(Number(0.), h_U - relax * h_max));
48 const auto test_max = hyperbolic_system.filter_dry_water_depth(
49 std::max(Number(0.), h_min - relax * h_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 water depth (critical)!\n";
55 std::cout <<
"\t\th min: " << h_min <<
"\n";
56 std::cout <<
"\t\th: " << h_U <<
"\n";
57 std::cout <<
"\t\th max: " << h_max <<
"\n" << std::endl;
62 const Number denominator =
63 h_P / std::max(std::abs(h_P * h_P), Number(100. *
min));
64 constexpr auto lte = dealii::SIMDComparison::less_than_or_equal;
66 t_r = dealii::compare_and_apply_mask<lte>(h_max,
72 t_r = dealii::compare_and_apply_mask<lte>(h_U + t_r * h_P,
79 t_r = std::min(t_r, t_max);
80 t_r = std::max(t_r, t_min);
85 t_r = dealii::compare_and_apply_mask<dealii::SIMDComparison::equal>(
86 hyperbolic_system.filter_dry_water_depth(h_U),
95 const auto h_new = hyperbolic_system.water_depth(U + t_r * P);
96 const auto test_new_min = hyperbolic_system.filter_dry_water_depth(
97 std::max(Number(0.), h_new - relax * h_max));
98 const auto test_new_max = hyperbolic_system.filter_dry_water_depth(
99 std::max(Number(0.), h_min - relax * h_new));
101 if (!(test_new_min == Number(0.) && test_new_max == Number(0.))) {
103 std::cout << std::fixed << std::setprecision(30);
104 std::cout <<
"Bounds violation: high-order water depth!\n";
105 std::cout <<
"\t\th min: " << h_min <<
"\n";
106 std::cout <<
"\t\th: " << h_new <<
"\n";
107 std::cout <<
"\t\th max: " << h_max <<
"\n";
108 std::cout <<
"\t\th_U: " << h_U <<
"\n";
109 std::cout <<
"\t\th_P: " << h_P <<
"\n";
110 std::cout <<
"\t\tt_r: " << t_r <<
"\n" << std::endl;
131 const auto &kin_max = std::get<2>(bounds);
135 const auto U_r = U + t_r * P;
136 const auto h_r = hyperbolic_system.water_depth(U_r);
137 const auto q_r = hyperbolic_system.momentum(U_r);
140 relax_small * h_r * kin_max -
ScalarNumber(0.5) * q_r.norm_square();
146 t_l = dealii::compare_and_apply_mask<
147 dealii::SIMDComparison::greater_than>(psi_r, Number(0.), t_r, t_l);
151 return {t_l, success};
153#ifdef DEBUG_OUTPUT_LIMITER
155 std::cout << std::endl;
156 std::cout << std::fixed << std::setprecision(16);
157 std::cout <<
"t_l: (start) " << t_l << std::endl;
158 std::cout <<
"t_r: (start) " << t_r << std::endl;
162 const auto U_l = U + t_l * P;
163 const auto h_l = hyperbolic_system.water_depth(U_l);
164 const auto q_l = hyperbolic_system.momentum(U_l);
167 relax_small * h_l * kin_max -
ScalarNumber(0.5) * q_l.norm_square();
177 const auto filtered_h_l = hyperbolic_system.filter_dry_water_depth(h_l);
178 const auto lower_bound =
179 (
ScalarNumber(1.) - relax) * filtered_h_l * kin_max - eps;
180 if (!(std::min(Number(0.), psi_l - lower_bound +
min) == Number(0.))) {
182 std::cout << std::fixed << std::setprecision(16);
184 <<
"Bounds violation: low-order kinetic energy (critical)!\n";
185 std::cout <<
"\t\tPsi left: 0 <= " << psi_l <<
"\n" << std::endl;
194 if (std::max(Number(0.), t_r - t_l - newton_tolerance) == Number(0.))
195 return {t_l, success};
214 const auto &h_P = hyperbolic_system.water_depth(P);
215 const auto &q_U = hyperbolic_system.momentum(U);
216 const auto &q_P = hyperbolic_system.momentum(P);
218 const auto dpsi_l = h_P * kin_max - (q_U * q_P) - q_P * q_P * t_l;
219 const auto dpsi_r = h_P * kin_max - (q_U * q_P) - q_P * q_P * t_r;
222 t_l, t_r, psi_l, psi_r, dpsi_l, dpsi_r, Number(-1.));
224#ifdef DEBUG_OUTPUT_LIMITER
225 if (std::max(Number(0.), psi_r + Number(eps)) == Number(0.)) {
226 std::cout <<
"psi_l: " << psi_l << std::endl;
227 std::cout <<
"psi_r: " << psi_r << std::endl;
228 std::cout <<
"dpsi_l: " << dpsi_l << std::endl;
229 std::cout <<
"dpsi_r: " << dpsi_r << std::endl;
230 std::cout <<
"t_l: (end) " << t_l << std::endl;
231 std::cout <<
"t_r: (end) " << t_r << std::endl;
240 const auto U_new = U + t_l * P;
241 const auto h_new = hyperbolic_system.water_depth(U_new);
242 const auto q_new = hyperbolic_system.momentum(U_new);
244 const auto psi_new = relax_small * h_new * kin_max -
247 const auto lower_bound =
250 const bool psi_valid =
251 std::min(Number(0.), psi_new - lower_bound +
min) == Number(0.);
254 std::cout << std::fixed << std::setprecision(16);
255 std::cout <<
"High-order bounds violation!\n";
256 std::cout <<
"\t\tDepth = " << h_new <<
"\n";
257 std::cout <<
"\t\tPsi: 0 <= " << psi_new <<
"\n" << std::endl;
265 return {t_l, success};
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
typename get_value_type< Number >::type ScalarNumber
HyperbolicSystem::state_type< dim, Number > state_type
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))
DEAL_II_ALWAYS_INLINE Number positive_part(const Number number)