ryujin 2.1.1 revision 9dcb748690310d6a540ebb8b066d1a0834fc7604
riemann_solver.template.h
Go to the documentation of this file.
1//
2// SPDX-License-Identifier: MIT
3// Copyright (C) 2020 - 2023 by the ryujin authors
4//
5
6#pragma once
7
8#include <compile_time_options.h>
9
10#include "riemann_solver.h"
11
12#include <newton.h>
13#include <simd.h>
14
15// #define DEBUG_RIEMANN_SOLVER
16
17namespace ryujin
18{
19 namespace ShallowWater
20 {
21 using namespace dealii;
22
23 template <int dim, typename Number>
24 DEAL_II_ALWAYS_INLINE inline Number
25 RiemannSolver<dim, Number>::f(const primitive_type &riemann_data_Z,
26 const Number &h) const
27 {
28 using ScalarNumber = typename get_value_type<Number>::type;
29 const ScalarNumber gravity = hyperbolic_system.gravity();
30
31 const auto &[h_Z, u_Z, a_Z] = riemann_data_Z;
32
33 const auto left_value = ScalarNumber(2.) * (std::sqrt(gravity * h) - a_Z);
34
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);
38
39 return dealii::compare_and_apply_mask<
40 dealii::SIMDComparison::less_than_or_equal>(
41 h, h_Z, left_value, right_value);
42 }
43
44
45 template <int dim, typename Number>
46 DEAL_II_ALWAYS_INLINE inline Number
47 RiemannSolver<dim, Number>::phi(const primitive_type &riemann_data_i,
48 const primitive_type &riemann_data_j,
49 const Number &h) const
50 {
51 const Number &u_i = riemann_data_i[1];
52 const Number &u_j = riemann_data_j[1];
53
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;
57#endif
58 return f(riemann_data_i, h) + f(riemann_data_j, h) + u_j - u_i;
59 }
60
61
62 template <int dim, typename Number>
63 DEAL_II_ALWAYS_INLINE inline Number
65 const primitive_type &riemann_data, const Number h_star) const
66 {
68
69 const auto &[h, u, a] = riemann_data;
70
71 const Number factor = positive_part((h_star - h) / h);
72 const Number half_factor = ScalarNumber(0.5) * factor;
73
74 return u - a * std::sqrt((ScalarNumber(1.) + half_factor) *
75 (ScalarNumber(1.) + factor));
76 }
77
78
79 template <int dim, typename Number>
80 DEAL_II_ALWAYS_INLINE inline Number
82 const Number h_star) const
83 {
85
86 const auto &[h, u, a] = riemann_data;
87
88 const Number factor = positive_part((h_star - h) / h);
89 const Number half_factor = ScalarNumber(0.5) * factor;
90
91 return u + a * std::sqrt((ScalarNumber(1.) + half_factor) *
92 (ScalarNumber(1.) + factor));
93 }
94
95
96 template <int dim, typename Number>
97 DEAL_II_ALWAYS_INLINE inline Number
99 const primitive_type &riemann_data_i,
100 const primitive_type &riemann_data_j,
101 const Number h_star) const
102 {
103 const Number lambda1 = lambda1_minus(riemann_data_i, h_star);
104 const Number lambda3 = lambda3_plus(riemann_data_j, h_star);
105
106 return std::max(negative_part(lambda1), positive_part(lambda3));
107 }
108
109
110 template <int dim, typename Number>
111 DEAL_II_ALWAYS_INLINE inline Number
113 const primitive_type &riemann_data_i,
114 const primitive_type &riemann_data_j) const
115 {
117 const ScalarNumber gravity = hyperbolic_system.gravity();
118 const auto gravity_inverse = ScalarNumber(1.) / gravity;
119
120 const auto &[h_i, u_i, a_i] = riemann_data_i;
121 const auto &[h_j, u_j, a_j] = riemann_data_j;
122
123 const Number h_min = std::min(h_i, h_j);
124 const Number h_max = std::max(h_i, h_j);
125
126#ifdef DEBUG_RIEMANN_SOLVER
127 std::cout << h_min << "<- h_min/max ->" << h_max << std::endl;
128#endif
129
130 const Number a_min = std::sqrt(gravity * h_min);
131 const Number a_max = std::sqrt(gravity * h_max);
132
133#ifdef DEBUG_RIEMANN_SOLVER
134 std::cout << a_min << "<- a_min/max ->" << a_max << std::endl;
135#endif
136
137 const Number sqrt_two = std::sqrt(ScalarNumber(2.));
138
139 /* x0 = (2 sqrt(2) - 1)^2 */
140 const Number x0 = Number(9.) - ScalarNumber(4.) * sqrt_two;
141
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;
146#endif
147
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;
152#endif
153
154
155 /* We compute the three h_star quantities */
156
157 Number tmp;
158
159 // h_star_left
160
161 tmp = positive_part(u_i - u_j + ScalarNumber(2.) * (a_i + a_j));
162 const Number h_star_left =
163 ScalarNumber(0.0625) * gravity_inverse * tmp * tmp;
164
165#ifdef DEBUG_RIEMANN_SOLVER
166 std::cout << "left: " << h_star_left << std::endl;
167#endif
168
169 // h_star_middle
170
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;
173
174#ifdef DEBUG_RIEMANN_SOLVER
175 std::cout << "middle: " << h_star_middle << std::endl;
176#endif
177
178 // h_star_right
179
180 const auto left_radicand =
181 ScalarNumber(3.) * h_min +
182 ScalarNumber(2.) * sqrt_two * std::sqrt(h_min * h_max);
183
184 const auto right_radicand =
185 sqrt_two * std::sqrt(gravity_inverse * h_min) * (u_i - u_j);
186
187 tmp = std::sqrt(positive_part(left_radicand + right_radicand));
188 tmp -= sqrt_two * std::sqrt(h_min);
189
190 const Number h_star_right = tmp * tmp;
191
192#ifdef DEBUG_RIEMANN_SOLVER
193 std::cout << "right: " << h_star_right << std::endl;
194#endif
195
196 /* Finally define h_star */
197
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);
201
202 h_star =
203 dealii::compare_and_apply_mask<dealii::SIMDComparison::less_than>(
204 phi_value_max, Number(0.), h_star_middle, h_star_right);
205
206 return h_star;
207 }
208
209
210 template <int dim, typename Number>
212 const primitive_type &riemann_data_i,
213 const primitive_type &riemann_data_j) const
214 {
215 const Number h_star =
216 h_star_two_rarefaction(riemann_data_i, riemann_data_j);
217
218 const Number lambda_max =
219 compute_lambda(riemann_data_i, riemann_data_j, h_star);
220
221 return lambda_max;
222 }
223
224 } // namespace ShallowWater
225} // namespace ryujin
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)
Definition: simd.h:123
DEAL_II_ALWAYS_INLINE Number positive_part(const Number number)
Definition: simd.h:111