ryujin 2.1.1 revision 075c649df0ea80b3fbb39a7a6e10bf95d4e28390
limiter.template.h
Go to the documentation of this file.
1//
2// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
3// Copyright (C) 2023 - 2024 by the ryujin authors
4//
5
6#pragma once
7
8#include "limiter.h"
9
10namespace ryujin
11{
12 namespace ScalarConservation
13 {
14 template <int dim, typename Number>
15 std::tuple<Number, bool>
17 const state_type &U,
18 const state_type &P,
19 const Number t_min /* = Number(0.) */,
20 const Number t_max /* = Number(1.) */)
21 {
22 const auto view = hyperbolic_system.view<dim, Number>();
23
24 bool success = true;
25 Number t_r = t_max;
26
27 constexpr ScalarNumber eps = std::numeric_limits<ScalarNumber>::epsilon();
28 const ScalarNumber relax = ScalarNumber(1. + 10000. * eps);
29
30 const auto &u_U = view.state(U);
31 const auto &u_P = view.state(P);
32
33 const auto &u_min = std::get<0>(bounds);
34 const auto &u_max = std::get<1>(bounds);
35
36 /*
37 * Verify that u_U is within bounds. This property might be
38 * violated for relative CFL numbers larger than 1.
39 *
40 * u_min, u_U, u_max might be negative, thus relax in both directions.
41 */
42 const auto test_max = std::max(
43 Number(0.), std::min(u_U - relax * u_max, relax * u_U - u_max));
44 const auto test_min = std::max(
45 Number(0.), std::min(u_min - relax * u_U, relax * u_min - u_U));
46 if (!(test_max == Number(0.) && test_min == Number(0.))) {
47#ifdef DEBUG_OUTPUT
48 std::cout << std::fixed << std::setprecision(16);
49 std::cout << "Bounds violation: low-order state (critical)!"
50 << "\n\t\tu min: " << u_min
51 << "\n\t\tu min (delta): " << negative_part(u_U - u_min)
52 << "\n\t\tu: " << u_U
53 << "\n\t\tu max (delta): " << positive_part(u_U - u_max)
54 << "\n\t\tu max: " << u_max << "\n"
55 << std::endl;
56#endif
57 success = false;
58 }
59
60 const auto regularization =
61 Number(100. * std::numeric_limits<ScalarNumber>::min());
62
63 const Number denominator =
64 ScalarNumber(1.) /
65 std::max(regularization, std::abs(u_P) + eps * u_max);
66
67 t_r = dealii::compare_and_apply_mask<dealii::SIMDComparison::less_than>(
68 u_max,
69 u_U + t_r * u_P,
70 /*
71 * u_P is positive.
72 *
73 * Note: Do not take an absolute value here. If we are out of
74 * bounds we have to ensure that t_r is set to t_min.
75 */
76 (u_max - u_U) * denominator,
77 t_r);
78
79 t_r = dealii::compare_and_apply_mask<dealii::SIMDComparison::less_than>(
80 u_U + t_r * u_P,
81 u_min,
82 /*
83 * u_P is negative.
84 *
85 * Note: Do not take an absolute value here. If we are out of
86 * bounds we have to ensure that t_r is set to t_min.
87 */
88 (u_U - u_min) * denominator,
89 t_r);
90
91 /*
92 * Ensure that t_min <= t <= t_max. This might not be the case if
93 * u_U is outside the interval [u_min, u_max]. Furthermore,
94 * the quotient we take above is prone to numerical cancellation in
95 * particular in the second pass of the limiter when u_P might be
96 * small.
97 */
98 t_r = std::min(t_r, t_max);
99 t_r = std::max(t_r, t_min);
100
101#ifdef CHECK_BOUNDS
102 /*
103 * Verify that the new state is within bounds:
104 *
105 * u_min, u_U, u_max might be negative, thus relax in both directions.
106 */
107 const auto u_new = view.state(U + t_r * P);
108 const auto test_new_max = std::max(
109 Number(0.), std::min(u_new - relax * u_max, relax * u_new - u_max));
110 const auto test_new_min = std::max(
111 Number(0.), std::min(u_min - relax * u_new, relax * u_min - u_new));
112 if (!(test_new_max == Number(0.) && test_new_min == Number(0.))) {
113#ifdef DEBUG_OUTPUT
114 std::cout << std::fixed << std::setprecision(16);
115 std::cout << "Bounds violation: high-order state!"
116 << "\n\t\tu min: " << u_min
117 << "\n\t\tu min (delta): " << negative_part(u_new - u_min)
118 << "\n\t\tu: " << u_new
119 << "\n\t\tu max (delta): " << positive_part(u_new - u_max)
120 << "\n\t\tu max: " << u_max << "\n"
121 << std::endl;
122#endif
123 success = false;
124 }
125#endif
126
127 return {t_r, success};
128 }
129
130 } // namespace ScalarConservation
131} // namespace ryujin
std::array< Number, n_bounds > Bounds
Definition: limiter.h:103
typename View::state_type state_type
Definition: limiter.h:65
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
Definition: limiter.h:61
DEAL_II_ALWAYS_INLINE Number negative_part(const Number number)
Definition: simd.h:124
DEAL_II_ALWAYS_INLINE Number positive_part(const Number number)
Definition: simd.h:112