14 namespace ShallowWater
16 template <
int dim,
typename Number>
17 std::tuple<Number, bool>
24 const auto view = hyperbolic_system.view<dim, Number>();
30 const auto &[h_min, h_max, h_small, kin_max, v2_max] = bounds;
33 constexpr ScalarNumber eps = std::numeric_limits<ScalarNumber>::epsilon();
34 const auto small = view.dry_state_relaxation_small();
35 const auto large = view.dry_state_relaxation_large();
46 auto h_U = view.water_depth(U);
47 const auto &h_P = view.water_depth(P);
49 const auto test_min = view.filter_dry_water_depth(
50 std::max(Number(0.), h_U - relax * h_max));
51 const auto test_max = view.filter_dry_water_depth(
52 std::max(Number(0.), h_min - relax * h_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 water depth (critical)!\n"
58 <<
"\n\t\th min: " << h_min
62 <<
"\n\t\th max: " << h_max <<
"\n"
68 const Number denominator =
71 constexpr auto lt = dealii::SIMDComparison::less_than;
73 t_r = dealii::compare_and_apply_mask<lt>(
82 (h_max - h_U) * denominator,
89 const auto h_min_tilde = std::max(h_small, h_min);
91 t_r = dealii::compare_and_apply_mask<lt>(
100 (h_U - h_min_tilde) * denominator,
110 t_r = std::min(t_r, t_max);
111 t_r = std::max(t_r, t_min);
114#ifdef EXPENSIVE_BOUNDS_CHECK
118 const auto h_new = view.water_depth(U + t_r * P);
119 const auto test_new_min = view.filter_dry_water_depth(
120 std::max(Number(0.), h_new - relax * h_max));
121 const auto test_new_max = view.filter_dry_water_depth(
122 std::max(Number(0.), h_min - relax * h_new));
124 if (!(test_new_min == Number(0.) && test_new_max == Number(0.))) {
126 std::cout << std::fixed << std::setprecision(30);
127 std::cout <<
"Bounds violation: high-order water depth!\n"
128 <<
"\n\t\th min: " << h_min
130 <<
"\n\t\th: " << h_new
132 <<
"\n\t\th max: " << h_max <<
"\n"
145 !parameters.limit_on_kinetic_energy()))
146 return {t_l, success};
161 if (parameters.limit_on_kinetic_energy()) {
165 const auto U_r = U + t_r * P;
166 const auto h_r = view.water_depth(U_r);
167 const auto q_r = view.momentum(U_r);
170 relax_small * h_r * kin_max -
ScalarNumber(0.5) * q_r.norm_square();
176 t_l = dealii::compare_and_apply_mask<
177 dealii::SIMDComparison::greater_than>(psi_r, Number(0.), t_r, t_l);
180 if (!parameters.limit_on_square_velocity() && t_l == t_r)
181 return {t_l, success};
183#ifdef DEBUG_OUTPUT_LIMITER
185 std::cout << std::endl;
186 std::cout << std::fixed << std::setprecision(16);
187 std::cout <<
"t_l: (start) " << t_l << std::endl;
188 std::cout <<
"t_r: (start) " << t_r << std::endl;
192 const auto U_l = U + t_l * P;
193 const auto h_l = view.water_depth(U_l);
194 const auto q_l = view.momentum(U_l);
197 relax_small * h_l * kin_max -
ScalarNumber(0.5) * q_l.norm_square();
207 const auto filtered_h_l = view.filter_dry_water_depth(h_l);
208 const auto lower_bound =
209 (
ScalarNumber(1.) - relax) * filtered_h_l * kin_max - eps;
210 if (!(std::min(Number(0.), psi_l - lower_bound) == Number(0.))) {
212 std::cout << std::fixed << std::setprecision(16);
214 <<
"Bounds violation: low-order kinetic energy (critical)!\n";
215 std::cout <<
"\t\tPsi left: 0 <= " << psi_l <<
"\n" << std::endl;
224 const Number tolerance(parameters.newton_tolerance());
225 if (!(std::max(Number(0.), t_r - t_l - tolerance) == Number(0.))) {
243 const auto &h_P = view.water_depth(P);
244 const auto &q_U = view.momentum(U);
245 const auto &q_P = view.momentum(P);
247 const auto dpsi_l = h_P * kin_max - (q_U * q_P) - q_P * q_P * t_l;
248 const auto dpsi_r = h_P * kin_max - (q_U * q_P) - q_P * q_P * t_r;
251 t_l, t_r, psi_l, psi_r, dpsi_l, dpsi_r, Number(-1.));
253#ifdef DEBUG_OUTPUT_LIMITER
254 if (std::max(Number(0.), psi_r + Number(eps)) == Number(0.)) {
255 std::cout <<
"psi_l: " << psi_l << std::endl;
256 std::cout <<
"psi_r: " << psi_r << std::endl;
257 std::cout <<
"dpsi_l: " << dpsi_l << std::endl;
258 std::cout <<
"dpsi_r: " << dpsi_r << std::endl;
259 std::cout <<
"t_l: (end) " << t_l << std::endl;
260 std::cout <<
"t_r: (end) " << t_r << std::endl;
265#ifdef EXPENSIVE_BOUNDS_CHECK
270 const auto U_new = U + t_l * P;
271 const auto h_new = view.water_depth(U_new);
272 const auto q_new = view.momentum(U_new);
274 const auto psi_new = relax_small * h_new * kin_max -
277 const auto lower_bound =
280 const bool psi_valid =
281 std::min(Number(0.), psi_new - lower_bound) == Number(0.);
284 std::cout << std::fixed << std::setprecision(16);
285 std::cout <<
"Bounds violation: high-order kinetic energy!\n";
286 std::cout <<
"\t\tPsi: 0 <= " << psi_new <<
"\n" << std::endl;
294 if (parameters.limit_on_square_velocity()) {
311 if (parameters.limit_on_square_velocity()) {
314 const auto U_r = U + t_r * P;
315 const auto h_r = view.water_depth(U_r);
316 const auto q_r = view.momentum(U_r);
318 const auto psi_r = relax_small * h_r * h_r * v2_max - q_r.norm_square();
324 t_l = dealii::compare_and_apply_mask<
325 dealii::SIMDComparison::greater_than>(psi_r, Number(0.), t_r, t_l);
329 return {t_l, success};
331#ifdef DEBUG_OUTPUT_LIMITER
333 std::cout << std::endl;
334 std::cout << std::fixed << std::setprecision(16);
335 std::cout <<
"t_l: (start) " << t_l << std::endl;
336 std::cout <<
"t_r: (start) " << t_r << std::endl;
340 const auto U_l = U + t_l * P;
341 const auto h_l = view.water_depth(U_l);
342 const auto q_l = view.momentum(U_l);
344 const auto psi_l = relax_small * h_l * h_l * v2_max - q_l.norm_square();
354 const auto filtered_h_l = view.filter_dry_water_depth(h_l);
355 const auto lower_bound =
356 (
ScalarNumber(1.) - relax) * filtered_h_l * filtered_h_l * v2_max -
358 if (!(std::min(Number(0.), psi_l - lower_bound) == Number(0.))) {
360 std::cout << std::fixed << std::setprecision(16);
362 <<
"Bounds violation: low-order square velocity (critical)!\n";
363 std::cout <<
"\t\tPsi left: 0 <= " << psi_l <<
"\n" << std::endl;
372 const Number tolerance(parameters.newton_tolerance());
373 if (!(std::max(Number(0.), t_r - t_l - tolerance) == Number(0.))) {
391 const auto &h_U = view.water_depth(U);
392 const auto &h_P = view.water_depth(P);
393 const auto &q_U = view.momentum(U);
394 const auto &q_P = view.momentum(P);
397 (h_U + t_l * h_P) * h_P * v2_max -
400 (h_U + t_r * h_P) * h_P * v2_max -
404 t_l, t_r, psi_l, psi_r, dpsi_l, dpsi_r, Number(-1.));
406#ifdef DEBUG_OUTPUT_LIMITER
407 if (std::max(Number(0.), psi_r + Number(eps)) == Number(0.)) {
408 std::cout <<
"psi_l: " << psi_l << std::endl;
409 std::cout <<
"psi_r: " << psi_r << std::endl;
410 std::cout <<
"dpsi_l: " << dpsi_l << std::endl;
411 std::cout <<
"dpsi_r: " << dpsi_r << std::endl;
412 std::cout <<
"t_l: (end) " << t_l << std::endl;
413 std::cout <<
"t_r: (end) " << t_r << std::endl;
418#ifdef EXPENSIVE_BOUNDS_CHECK
423 const auto U_new = U + t_l * P;
424 const auto h_new = view.water_depth(U_new);
425 const auto q_new = view.momentum(U_new);
428 relax_small * h_new * h_new * v2_max - q_new.norm_square();
430 const auto lower_bound =
434 const bool psi_valid =
435 std::min(Number(0.), psi_new - lower_bound) == Number(0.);
438 std::cout << std::fixed << std::setprecision(16);
439 std::cout <<
"Bounds violation: high-order square velocity!\n";
440 std::cout <<
"\t\tPsi: 0 <= " << psi_new <<
"\n" << std::endl;
448 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.))
typename View::ScalarNumber ScalarNumber
std::array< Number, n_bounds > Bounds
typename View::state_type state_type
#define RYUJIN_UNLIKELY(x)
DEAL_II_ALWAYS_INLINE Number negative_part(const Number number)
DEAL_II_ALWAYS_INLINE Number positive_part(const Number number)
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))