8#include <compile_time_options.h>
50 template <
int dim,
typename Number>
51 DEAL_II_ALWAYS_INLINE
inline Number
70 const Number first_radicand = (
ScalarNumber(3.) * gamma + Number(11.)) /
73 const Number second_radicand =
74 Number(5. / 6.) + slope * (gamma - Number(3.));
76 Number radicand = std::min(first_radicand, second_radicand);
77 radicand = std::min(Number(1.), radicand);
78 radicand = std::max(Number(1. / 2.), radicand);
80 return std::sqrt(radicand);
84 template <
int dim,
typename Number>
86 const Number &rho,
const Number &gamma,
const Number &a)
const
88 const auto view = hyperbolic_system.view<dim, Number>();
89 const auto interpolation_b = view.eos_interpolation_b();
91 const Number numerator =
92 ScalarNumber(2.) * a * (Number(1.) - interpolation_b * rho);
94 const Number denominator = gamma - Number(1.);
96 return numerator / denominator;
100 template <
int dim,
typename Number>
101 DEAL_II_ALWAYS_INLINE
inline Number
106 const auto &[rho_i, u_i, p_i, gamma_i, a_i] = riemann_data_i;
107 const auto &[rho_j, u_j, p_j, gamma_j, a_j] = riemann_data_j;
108 const auto alpha_i = alpha(rho_i, gamma_i, a_i);
109 const auto alpha_j = alpha(rho_j, gamma_j, a_j);
119 const Number p_min = std::min(p_i, p_j);
120 const Number p_max = std::max(p_i, p_j);
122 const Number gamma_min =
123 dealii::compare_and_apply_mask<dealii::SIMDComparison::less_than>(
124 p_i, p_j, gamma_i, gamma_j);
126 const Number alpha_min =
127 dealii::compare_and_apply_mask<dealii::SIMDComparison::less_than>(
128 p_i, p_j, alpha_i, alpha_j);
130 const Number alpha_hat_min = c(gamma_min) * alpha_min;
132 const Number alpha_max = dealii::compare_and_apply_mask<
133 dealii::SIMDComparison::greater_than_or_equal>(
134 p_i, p_j, alpha_i, alpha_j);
136 const Number gamma_m = std::min(gamma_i, gamma_j);
137 const Number gamma_M = std::max(gamma_i, gamma_j);
139 const Number numerator =
142 const Number p_ratio = p_min / p_max;
150 const Number r_exponent =
151 (gamma_M - gamma_min) / (
ScalarNumber(2.) * gamma_min * gamma_M);
158 const Number first_exponent =
160 const Number first_exponent_inverse = Number(1.) / first_exponent;
162 const Number first_denom =
163 alpha_hat_min *
ryujin::pow(p_ratio, r_exponent - first_exponent) +
166 const Number p_1_tilde =
167 p_max *
ryujin::pow(numerator / first_denom, first_exponent_inverse);
168#ifdef DEBUG_RIEMANN_SOLVER
169 std::cout <<
"RS p_1_tilde = " << p_1_tilde <<
"\n";
177 const Number second_exponent =
179 const Number second_exponent_inverse = Number(1.) / second_exponent;
181 Number second_denom =
182 alpha_hat_min *
ryujin::pow(p_ratio, -second_exponent) +
185 const Number p_2_tilde = p_max *
ryujin::pow(numerator / second_denom,
186 second_exponent_inverse);
188#ifdef DEBUG_RIEMANN_SOLVER
189 std::cout <<
"RS p_2_tilde = " << p_2_tilde <<
"\n";
192 return std::min(p_1_tilde, p_2_tilde);
196 template <
int dim,
typename Number>
197 DEAL_II_ALWAYS_INLINE
inline Number
202 const auto &[rho_i, u_i, p_i, gamma_i, a_i] = riemann_data_i;
203 const auto &[rho_j, u_j, p_j, gamma_j, a_j] = riemann_data_j;
205 const Number gamma_m = std::min(gamma_i, gamma_j);
207 const Number alpha_hat_i = c(gamma_i) * alpha(rho_i, gamma_i, a_i);
208 const Number alpha_hat_j = c(gamma_j) * alpha(rho_j, gamma_j, a_j);
216 const Number exponent =
218 const Number exponent_inverse = Number(1.) / exponent;
220 const Number numerator =
223 const Number denominator =
224 alpha_hat_i *
ryujin::pow(p_i / p_j, -exponent) + alpha_hat_j;
226 const Number p_1_tilde =
227 p_j *
ryujin::pow(numerator / denominator, exponent_inverse);
229#ifdef DEBUG_RIEMANN_SOLVER
230 std::cout <<
"SS p_1_tilde = " << p_1_tilde <<
"\n";
233 const auto p_2_tilde = p_star_failsafe(riemann_data_i, riemann_data_j);
235 return std::min(p_1_tilde, p_2_tilde);
239 template <
int dim,
typename Number>
240 DEAL_II_ALWAYS_INLINE
inline Number
245 const auto view = hyperbolic_system.view<dim, Number>();
246 const auto interpolation_b = view.eos_interpolation_b();
248 const auto &[rho_i, u_i, p_i, gamma_i, a_i] = riemann_data_i;
249 const auto &[rho_j, u_j, p_j, gamma_j, a_j] = riemann_data_j;
257 const Number p_max = std::max(p_i, p_j);
260 ScalarNumber(2.) * (Number(1.) - interpolation_b * rho_i) * p_max;
261 radicand_i /= rho_i * ((gamma_i + Number(1.)) * p_max +
262 (gamma_i - Number(1.)) * p_i);
264 const Number x_i = std::sqrt(radicand_i);
267 ScalarNumber(2.) * (Number(1.) - interpolation_b * rho_j) * p_max;
268 radicand_j /= rho_j * ((gamma_j + Number(1.)) * p_max +
269 (gamma_j - Number(1.)) * p_j);
271 const Number x_j = std::sqrt(radicand_j);
273 const Number a = x_i + x_j;
274 const Number b = u_j - u_i;
275 const Number c = -p_i * x_i - p_j * x_j;
277 const Number base = (-b + std::sqrt(b * b -
ScalarNumber(4.) * a * c)) /
279 const Number p_2_tilde = base * base;
281#ifdef DEBUG_RIEMANN_SOLVER
282 std::cout <<
"SS p_2_tilde = " << p_2_tilde <<
"\n";
288 template <
int dim,
typename Number>
289 DEAL_II_ALWAYS_INLINE
inline Number
294 const auto &[rho_i, u_i, p_i, gamma_i, a_i] = riemann_data_i;
295 const auto &[rho_j, u_j, p_j, gamma_j, a_j] = riemann_data_j;
296 const auto alpha_i = alpha(rho_i, gamma_i, a_i);
297 const auto alpha_j = alpha(rho_j, gamma_j, a_j);
307 const Number p_min = std::min(p_i, p_j);
308 const Number p_max = std::max(p_i, p_j);
310 const Number gamma_min =
311 dealii::compare_and_apply_mask<dealii::SIMDComparison::less_than>(
312 p_i, p_j, gamma_i, gamma_j);
314 const Number alpha_min =
315 dealii::compare_and_apply_mask<dealii::SIMDComparison::less_than>(
316 p_i, p_j, alpha_i, alpha_j);
318 const Number alpha_hat_min = c(gamma_min) * alpha_min;
320 const Number gamma_max = dealii::compare_and_apply_mask<
321 dealii::SIMDComparison::greater_than_or_equal>(
322 p_i, p_j, gamma_i, gamma_j);
324 const Number alpha_max = dealii::compare_and_apply_mask<
325 dealii::SIMDComparison::greater_than_or_equal>(
326 p_i, p_j, alpha_i, alpha_j);
328 const Number alpha_hat_max = c(gamma_max) * alpha_max;
330 const Number gamma_m = std::min(gamma_i, gamma_j);
331 const Number gamma_M = std::max(gamma_i, gamma_j);
333 const Number p_ratio = p_min / p_max;
341 const Number r_exponent =
342 (gamma_M - gamma_min) / (
ScalarNumber(2.) * gamma_min * gamma_M);
351 const Number exponent =
353 const Number exponent_inverse = Number(1.) / exponent;
355 const Number numerator =
358 Number denominator = alpha_hat_min *
ryujin::pow(p_ratio, -exponent) +
361 const Number p_tilde =
362 p_max *
ryujin::pow(numerator / denominator, exponent_inverse);
364#ifdef DEBUG_RIEMANN_SOLVER
365 std::cout <<
"IN p_*_tilde = " << p_tilde <<
"\n";
372 template <
int dim,
typename Number>
373 DEAL_II_ALWAYS_INLINE
inline Number
375 const Number p_star)
const
377 const auto view = hyperbolic_system.view<dim, Number>();
378 const auto interpolation_b = view.eos_interpolation_b();
380 const auto &[rho, u, p, gamma, a] = riemann_data;
382 const Number one_minus_b_rho = Number(1.) - interpolation_b * rho;
385 ScalarNumber(2.) * one_minus_b_rho / (rho * (gamma + Number(1.)));
387 const Number Bz = (gamma - Number(1.)) / (gamma + Number(1.)) * p;
389 const Number radicand = Az / (p_star + Bz);
391 const Number true_value = (p_star - p) * std::sqrt(radicand);
393 const auto exponent = ScalarNumber(0.5) * (gamma - Number(1.)) / gamma;
394 const Number factor =
ryujin::pow(p_star / p, exponent) - Number(1.);
395 const auto false_value = ScalarNumber(2.) * a * one_minus_b_rho * factor /
396 (gamma - Number(1.));
398 return dealii::compare_and_apply_mask<
399 dealii::SIMDComparison::greater_than_or_equal>(
400 p_star, p, true_value, false_value);
404 template <
int dim,
typename Number>
405 DEAL_II_ALWAYS_INLINE
inline Number
407 const primitive_type &riemann_data_j,
408 const Number p_in)
const
410 const Number &u_i = riemann_data_i[1];
411 const Number &u_j = riemann_data_j[1];
413 return f(riemann_data_i, p_in) + f(riemann_data_j, p_in) + u_j - u_i;
417 template <
int dim,
typename Number>
418 DEAL_II_ALWAYS_INLINE
inline Number
423 const auto view = hyperbolic_system.view<dim, Number>();
424 const auto interpolation_b = view.eos_interpolation_b();
426 const auto &[rho_i, u_i, p_i, gamma_i, a_i] = riemann_data_i;
427 const auto &[rho_j, u_j, p_j, gamma_j, a_j] = riemann_data_j;
429 const Number p_max = std::max(p_i, p_j);
431 const Number radicand_inverse_i =
432 ScalarNumber(0.5) * rho_i / (Number(1.) - interpolation_b * rho_i) *
433 ((gamma_i + Number(1.)) * p_max + (gamma_i - Number(1.)) * p_i);
435 const Number value_i = (p_max - p_i) / std::sqrt(radicand_inverse_i);
437 const Number radicand_jnverse_j =
438 ScalarNumber(0.5) * rho_j / (Number(1.) - interpolation_b * rho_j) *
439 ((gamma_j + Number(1.)) * p_max + (gamma_j - Number(1.)) * p_j);
441 const Number value_j = (p_max - p_j) / std::sqrt(radicand_jnverse_j);
443 return value_i + value_j + u_j - u_i;
447 template <
int dim,
typename Number>
448 DEAL_II_ALWAYS_INLINE
inline Number
452 const auto &[rho, u, p, gamma, a] = riemann_data;
459 return u - a * std::sqrt(Number(1.) + factor * tmp);
463 template <
int dim,
typename Number>
464 DEAL_II_ALWAYS_INLINE
inline Number
466 const Number p_star)
const
468 const auto &[rho, u, p, gamma, a] = riemann_data;
475 return u + a * std::sqrt(Number(1.) + factor * tmp);
479 template <
int dim,
typename Number>
480 DEAL_II_ALWAYS_INLINE
inline Number
484 const Number p_star)
const
486 const Number nu_11 = lambda1_minus(riemann_data_i, p_star);
487 const Number nu_32 = lambda3_plus(riemann_data_j, p_star);
493 template <
int dim,
typename Number>
494 DEAL_II_ALWAYS_INLINE
inline auto
498 const dealii::Tensor<1, dim, Number> &n_ij)
const ->
primitive_type
500 const auto view = hyperbolic_system.view<dim, Number>();
502 const auto rho = view.density(U);
505 const auto m = view.momentum(U);
506 const auto proj_m = n_ij * m;
508 const auto gamma = view.surrogate_gamma(U, p);
510 const auto interpolation_b = view.eos_interpolation_b();
511 const auto x = Number(1.) - interpolation_b * rho;
512 const auto a = std::sqrt(gamma * p / (rho * x));
518 dealii::ExcMessage(
"Internal error: p <= 0."));
523 dealii::ExcMessage(
"Internal error: 1. - b * rho <= 0."));
528 dealii::ExcMessage(
"Internal error: gamma <= 1."));
531 return {{rho, proj_m * rho_inverse, p, gamma, a}};
535 template <
int dim,
typename Number>
540 const auto view = hyperbolic_system.view<dim, Number>();
542 const auto &[rho_i, u_i, p_i, gamma_i, a_i] = riemann_data_i;
543 const auto &[rho_j, u_j, p_j, gamma_j, a_j] = riemann_data_j;
545#ifdef DEBUG_RIEMANN_SOLVER
546 std::cout <<
"rho_left: " << rho_i << std::endl;
547 std::cout <<
"u_left: " << u_i << std::endl;
548 std::cout <<
"p_left: " << p_i << std::endl;
549 std::cout <<
"gamma_left: " << gamma_i << std::endl;
550 std::cout <<
"a_left: " << a_i << std::endl;
551 std::cout <<
"rho_right: " << rho_j << std::endl;
552 std::cout <<
"u_right: " << u_j << std::endl;
553 std::cout <<
"p_right: " << p_j << std::endl;
554 std::cout <<
"gamma_right: " << gamma_j << std::endl;
555 std::cout <<
"a_right: " << a_j << std::endl;
558 const Number p_max = std::max(p_i, p_j);
559 const Number phi_p_max = phi_of_p_max(riemann_data_i, riemann_data_j);
561 if (!view.compute_strict_bounds()) {
562#ifdef DEBUG_RIEMANN_SOLVER
563 const Number p_star_RS = p_star_RS_full(riemann_data_i, riemann_data_j);
564 const Number p_star_SS = p_star_SS_full(riemann_data_i, riemann_data_j);
565 const Number p_debug =
566 dealii::compare_and_apply_mask<dealii::SIMDComparison::less_than>(
567 phi_p_max, Number(0.), p_star_SS, std::min(p_max, p_star_RS));
568 std::cout <<
" p^*_debug = " << p_debug <<
"\n";
569 std::cout <<
" phi(p_*_d) = "
570 << phi(riemann_data_i, riemann_data_j, p_debug) <<
"\n";
571 std::cout <<
"-> lambda_deb = "
572 << compute_lambda(riemann_data_i, riemann_data_j, p_debug)
576 const Number p_star_tilde =
577 p_star_interpolated(riemann_data_i, riemann_data_j);
578 const Number p_star_backup =
579 p_star_failsafe(riemann_data_i, riemann_data_j);
582 dealii::compare_and_apply_mask<dealii::SIMDComparison::less_than>(
585 std::min(p_star_tilde, p_star_backup),
586 std::min(p_max, p_star_tilde));
588#ifdef DEBUG_RIEMANN_SOLVER
589 std::cout <<
" p^*_tilde = " << p_2 <<
"\n";
590 std::cout <<
" phi(p_*_t) = "
591 << phi(riemann_data_i, riemann_data_j, p_2) <<
"\n";
592 std::cout <<
"-> lambda_max = "
593 << compute_lambda(riemann_data_i, riemann_data_j, p_2) <<
"\n"
597 return compute_lambda(riemann_data_i, riemann_data_j, p_2);
600 const Number p_star_RS = p_star_RS_full(riemann_data_i, riemann_data_j);
601 const Number p_star_SS = p_star_SS_full(riemann_data_i, riemann_data_j);
604 dealii::compare_and_apply_mask<dealii::SIMDComparison::less_than>(
605 phi_p_max, Number(0.), p_star_SS, std::min(p_max, p_star_RS));
607#ifdef DEBUG_RIEMANN_SOLVER
608 std::cout <<
" p^*_tilde = " << p_2 <<
"\n";
609 std::cout <<
" phi(p_*_t) = "
610 << phi(riemann_data_i, riemann_data_j, p_2) <<
"\n";
611 std::cout <<
"-> lambda_max = "
612 << compute_lambda(riemann_data_i, riemann_data_j, p_2)
616 return compute_lambda(riemann_data_i, riemann_data_j, p_2);
620 template <
int dim,
typename Number>
624 const unsigned int i,
625 const unsigned int *js,
626 const dealii::Tensor<1, dim, Number> &n_ij)
const
628 const auto &[p_i, unused_i, s_i, eta_i] =
630 .template get_tensor<Number, precomputed_state_type>(i);
632 const auto &[p_j, unused_j, s_j, eta_j] =
634 .template get_tensor<Number, precomputed_state_type>(js);
636 const auto riemann_data_i = riemann_data_from_state(U_i, p_i, n_ij);
637 const auto riemann_data_j = riemann_data_from_state(U_j, p_j, n_ij);
639 return compute(riemann_data_i, riemann_data_j);
Number p_star_RS_full(const primitive_type &riemann_data_i, const primitive_type &riemann_data_j) const
Number lambda1_minus(const primitive_type &riemann_data, const Number p_star) const
Number alpha(const Number &rho, const Number &gamma, const Number &a) const
primitive_type riemann_data_from_state(const state_type &U, const Number &p, const dealii::Tensor< 1, dim, Number > &n_ij) const
Number p_star_failsafe(const primitive_type &riemann_data_i, const primitive_type &riemann_data_j) const
Number p_star_SS_full(const primitive_type &riemann_data_i, const primitive_type &riemann_data_j) const
Number lambda3_plus(const primitive_type &primitive_state, const Number p_star) const
typename View::state_type state_type
Number p_star_interpolated(const primitive_type &riemann_data_i, const primitive_type &riemann_data_j) const
Number phi_of_p_max(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 p_star) const
std::array< Number, riemann_data_size > primitive_type
typename View::ScalarNumber ScalarNumber
Number c(const Number gamma_Z) const
Number compute(const primitive_type &riemann_data_i, const primitive_type &riemann_data_j) const
#define AssertThrowSIMD(variable, condition, exception)
T pow(const T x, const T b)
DEAL_II_ALWAYS_INLINE Number negative_part(const Number number)
DEAL_II_ALWAYS_INLINE Number positive_part(const Number number)