8#include <compile_time_options.h>
19 namespace ShallowWater
21 using namespace dealii;
23 template <
int dim,
typename Number>
24 DEAL_II_ALWAYS_INLINE
inline Number
26 const Number &h)
const
29 const ScalarNumber gravity = hyperbolic_system.gravity();
31 const auto &[h_Z, u_Z, a_Z] = riemann_data_Z;
33 const auto left_value = ScalarNumber(2.) * (std::sqrt(gravity * h) - a_Z);
35 const Number radicand =
36 ScalarNumber(0.5) * gravity * (h + h_Z) / (h * h_Z);
37 const Number right_value = (h - h_Z) * std::sqrt(radicand);
39 return dealii::compare_and_apply_mask<
40 dealii::SIMDComparison::less_than_or_equal>(
41 h, h_Z, left_value, right_value);
45 template <
int dim,
typename Number>
46 DEAL_II_ALWAYS_INLINE
inline Number
48 const primitive_type &riemann_data_j,
49 const Number &h)
const
51 const Number &u_i = riemann_data_i[1];
52 const Number &u_j = riemann_data_j[1];
54#ifdef DEBUG_RIEMANN_SOLVER
55 std::cout <<
"f_L --> " << f(riemann_data_i, h) << std::endl;
56 std::cout <<
"f_R --> " << f(riemann_data_j, h) << std::endl;
58 return f(riemann_data_i, h) + f(riemann_data_j, h) + u_j - u_i;
62 template <
int dim,
typename Number>
63 DEAL_II_ALWAYS_INLINE
inline Number
69 const auto &[h, u, a] = riemann_data;
74 return u - a * std::sqrt((
ScalarNumber(1.) + half_factor) *
79 template <
int dim,
typename Number>
80 DEAL_II_ALWAYS_INLINE
inline Number
82 const Number h_star)
const
86 const auto &[h, u, a] = riemann_data;
91 return u + a * std::sqrt((
ScalarNumber(1.) + half_factor) *
96 template <
int dim,
typename Number>
97 DEAL_II_ALWAYS_INLINE
inline Number
101 const Number h_star)
const
103 const Number lambda1 = lambda1_minus(riemann_data_i, h_star);
104 const Number lambda3 = lambda3_plus(riemann_data_j, h_star);
110 template <
int dim,
typename Number>
111 DEAL_II_ALWAYS_INLINE
inline Number
117 const ScalarNumber gravity = hyperbolic_system.gravity();
118 const auto gravity_inverse =
ScalarNumber(1.) / gravity;
120 const auto &[h_i, u_i, a_i] = riemann_data_i;
121 const auto &[h_j, u_j, a_j] = riemann_data_j;
123 const Number h_min = std::min(h_i, h_j);
124 const Number h_max = std::max(h_i, h_j);
126#ifdef DEBUG_RIEMANN_SOLVER
127 std::cout << h_min <<
"<- h_min/max ->" << h_max << std::endl;
130 const Number a_min = std::sqrt(gravity * h_min);
131 const Number a_max = std::sqrt(gravity * h_max);
133#ifdef DEBUG_RIEMANN_SOLVER
134 std::cout << a_min <<
"<- a_min/max ->" << a_max << std::endl;
140 const Number x0 = Number(9.) -
ScalarNumber(4.) * sqrt_two;
142 const Number phi_value_min =
143 phi(riemann_data_i, riemann_data_j, x0 * h_min);
144#ifdef DEBUG_RIEMANN_SOLVER
145 std::cout <<
"phi_value_min ->" << phi_value_min << std::endl;
148 const Number phi_value_max =
149 phi(riemann_data_i, riemann_data_j, x0 * h_max);
150#ifdef DEBUG_RIEMANN_SOLVER
151 std::cout <<
"phi_value_max ->" << phi_value_max << std::endl;
162 const Number h_star_left =
165#ifdef DEBUG_RIEMANN_SOLVER
166 std::cout <<
"left: " << h_star_left << std::endl;
171 tmp = Number(1.) + sqrt_two * (u_i - u_j) / (a_min + a_max);
172 const Number h_star_middle = std::sqrt(h_min * h_max) * tmp;
174#ifdef DEBUG_RIEMANN_SOLVER
175 std::cout <<
"middle: " << h_star_middle << std::endl;
180 const auto left_radicand =
184 const auto right_radicand =
185 sqrt_two * std::sqrt(gravity_inverse * h_min) * (u_i - u_j);
187 tmp = std::sqrt(
positive_part(left_radicand + right_radicand));
188 tmp -= sqrt_two * std::sqrt(h_min);
190 const Number h_star_right = tmp * tmp;
192#ifdef DEBUG_RIEMANN_SOLVER
193 std::cout <<
"right: " << h_star_right << std::endl;
198 Number h_star = dealii::compare_and_apply_mask<
199 dealii::SIMDComparison::less_than_or_equal>(
200 Number(0.), phi_value_min, h_star_left, h_star_right);
203 dealii::compare_and_apply_mask<dealii::SIMDComparison::less_than>(
204 phi_value_max, Number(0.), h_star_middle, h_star_right);
210 template <
int dim,
typename Number>
215 const Number h_star =
216 h_star_two_rarefaction(riemann_data_i, riemann_data_j);
218 const Number lambda_max =
219 compute_lambda(riemann_data_i, riemann_data_j, h_star);
Number phi(const primitive_type &riemann_data_i, const primitive_type &riemann_data_j, const Number &h) const
Number lambda3_plus(const primitive_type &riemann_data, const Number h_star) const
Number lambda1_minus(const primitive_type &riemann_data, const Number h_star) const
Number h_star_two_rarefaction(const primitive_type &riemann_data_i, const primitive_type &riemann_data_j) const
Number compute_lambda(const primitive_type &riemann_data_i, const primitive_type &riemann_data_j, const Number h_star) const
std::array< Number, riemann_data_size > primitive_type
typename get_value_type< Number >::type ScalarNumber
Number compute(const primitive_type &riemann_data_i, const primitive_type &riemann_data_j) const
Number f(const primitive_type &primitive_state, const Number &h_star) const
DEAL_II_ALWAYS_INLINE Number negative_part(const Number number)
DEAL_II_ALWAYS_INLINE Number positive_part(const Number number)