49 template <
int dim,
typename Number>
50 DEAL_II_ALWAYS_INLINE
inline Number
69 const Number first_radicand = (
ScalarNumber(3.) * gamma + Number(11.)) /
72 const Number second_radicand =
73 Number(5. / 6.) + slope * (gamma - Number(3.));
75 Number radicand = std::min(first_radicand, second_radicand);
76 radicand = std::min(Number(1.), radicand);
77 radicand = std::max(Number(1. / 2.), radicand);
79 return std::sqrt(radicand);
83 template <
int dim,
typename Number>
85 const Number &rho,
const Number &gamma,
const Number &a)
const
87 const auto view = hyperbolic_system.view<dim, Number>();
88 const auto interpolation_b = view.eos_interpolation_b();
90 const Number numerator =
91 ScalarNumber(2.) * a * (Number(1.) - interpolation_b * rho);
93 const Number denominator = gamma - Number(1.);
99 template <
int dim,
typename Number>
100 DEAL_II_ALWAYS_INLINE
inline Number
105 const auto view = hyperbolic_system.view<dim, Number>();
106 const auto pinf = view.eos_interpolation_pinfty();
108 const auto &[rho_i, u_i, p_i, gamma_i, a_i] = riemann_data_i;
109 const auto &[rho_j, u_j, p_j, gamma_j, a_j] = riemann_data_j;
110 const auto alpha_i = alpha(rho_i, gamma_i, a_i);
111 const auto alpha_j = alpha(rho_j, gamma_j, a_j);
121 const Number p_min = std::min(p_i, p_j);
122 const Number p_max = std::max(p_i, p_j);
124 const Number gamma_min =
125 dealii::compare_and_apply_mask<dealii::SIMDComparison::less_than>(
126 p_i, p_j, gamma_i, gamma_j);
128 const Number alpha_min =
129 dealii::compare_and_apply_mask<dealii::SIMDComparison::less_than>(
130 p_i, p_j, alpha_i, alpha_j);
132 const Number alpha_hat_min = c(gamma_min) * alpha_min;
134 const Number alpha_max = dealii::compare_and_apply_mask<
135 dealii::SIMDComparison::greater_than_or_equal>(
136 p_i, p_j, alpha_i, alpha_j);
138 const Number gamma_m = std::min(gamma_i, gamma_j);
139 const Number gamma_M = std::max(gamma_i, gamma_j);
141 const Number numerator =
142 dealii::compare_and_apply_mask<dealii::SIMDComparison::equal>(
152 const Number p_ratio =
safe_division(p_min + pinf, p_max + pinf);
160 const Number r_exponent =
161 (gamma_M - gamma_min) / (
ScalarNumber(2.) * gamma_min * gamma_M);
168 const Number first_exponent =
171 const Number first_exponent_inverse =
174 const Number first_denom =
175 alpha_hat_min *
ryujin::pow(p_ratio, r_exponent - first_exponent) +
178 const Number p_1_tilde =
180 first_exponent_inverse) -
183#ifdef DEBUG_RIEMANN_SOLVER
184 std::cout <<
"RS p_1_tilde = " << p_1_tilde <<
"\n";
192 const Number second_exponent =
195 const Number second_exponent_inverse =
198 Number second_denom =
199 alpha_hat_min *
ryujin::pow(p_ratio, -second_exponent) +
202 const Number p_2_tilde =
204 second_exponent_inverse) -
207#ifdef DEBUG_RIEMANN_SOLVER
208 std::cout <<
"RS p_2_tilde = " << p_2_tilde <<
"\n";
211 return std::min(p_1_tilde, p_2_tilde);
215 template <
int dim,
typename Number>
216 DEAL_II_ALWAYS_INLINE
inline Number
221 const auto view = hyperbolic_system.view<dim, Number>();
222 const auto pinf = view.eos_interpolation_pinfty();
224 const auto &[rho_i, u_i, p_i, gamma_i, a_i] = riemann_data_i;
225 const auto &[rho_j, u_j, p_j, gamma_j, a_j] = riemann_data_j;
227 const Number gamma_m = std::min(gamma_i, gamma_j);
229 const Number alpha_hat_i = c(gamma_i) * alpha(rho_i, gamma_i, a_i);
230 const Number alpha_hat_j = c(gamma_j) * alpha(rho_j, gamma_j, a_j);
238 const Number exponent =
240 const Number exponent_inverse = Number(1.) / exponent;
242 const Number numerator =
243 dealii::compare_and_apply_mask<dealii::SIMDComparison::equal>(
249 const Number denominator =
254 const Number p_1_tilde =
259#ifdef DEBUG_RIEMANN_SOLVER
260 std::cout <<
"SS p_1_tilde = " << p_1_tilde <<
"\n";
263 const auto p_2_tilde = p_star_failsafe(riemann_data_i, riemann_data_j);
265 return std::min(p_1_tilde, p_2_tilde);
269 template <
int dim,
typename Number>
270 DEAL_II_ALWAYS_INLINE
inline Number
275 const auto view = hyperbolic_system.view<dim, Number>();
276 const auto interpolation_b = view.eos_interpolation_b();
277 const auto pinf = view.eos_interpolation_pinfty();
279 const auto &[rho_i, u_i, p_i, gamma_i, a_i] = riemann_data_i;
280 const auto &[rho_j, u_j, p_j, gamma_j, a_j] = riemann_data_j;
288 const Number p_max = std::max(p_i, p_j) + pinf;
291 ScalarNumber(2.) * (Number(1.) - interpolation_b * rho_i) * p_max,
292 rho_i * ((gamma_i + Number(1.)) * p_max +
293 (gamma_i - Number(1.)) * (p_i + pinf)));
295 const Number x_i = std::sqrt(radicand_i);
298 ScalarNumber(2.) * (Number(1.) - interpolation_b * rho_j) * p_max,
299 rho_j * ((gamma_j + Number(1.)) * p_max +
300 (gamma_j - Number(1.)) * (p_j + pinf)));
302 const Number x_j = std::sqrt(radicand_j);
304 const Number a = x_i + x_j;
306 dealii::compare_and_apply_mask<dealii::SIMDComparison::equal>(
307 a, Number(0.), Number(0.), u_j - u_i);
309 const Number c = -(p_i + pinf) * x_i - (p_j + pinf) * x_j;
316 const Number p_2_tilde = base * base - pinf;
318#ifdef DEBUG_RIEMANN_SOLVER
319 std::cout <<
"SS p_2_tilde = " << p_2_tilde <<
"\n";
325 template <
int dim,
typename Number>
326 DEAL_II_ALWAYS_INLINE
inline Number
331 const auto view = hyperbolic_system.view<dim, Number>();
332 const auto pinf = view.eos_interpolation_pinfty();
334 const auto &[rho_i, u_i, p_i, gamma_i, a_i] = riemann_data_i;
335 const auto &[rho_j, u_j, p_j, gamma_j, a_j] = riemann_data_j;
336 const auto alpha_i = alpha(rho_i, gamma_i, a_i);
337 const auto alpha_j = alpha(rho_j, gamma_j, a_j);
347 const Number p_min = std::min(p_i, p_j) + pinf;
348 const Number p_max = std::max(p_i, p_j) + pinf;
350 const Number gamma_min =
351 dealii::compare_and_apply_mask<dealii::SIMDComparison::less_than>(
352 p_i, p_j, gamma_i, gamma_j);
354 const Number alpha_min =
355 dealii::compare_and_apply_mask<dealii::SIMDComparison::less_than>(
356 p_i, p_j, alpha_i, alpha_j);
358 const Number alpha_hat_min = c(gamma_min) * alpha_min;
360 const Number gamma_max = dealii::compare_and_apply_mask<
361 dealii::SIMDComparison::greater_than_or_equal>(
362 p_i, p_j, gamma_i, gamma_j);
364 const Number alpha_max = dealii::compare_and_apply_mask<
365 dealii::SIMDComparison::greater_than_or_equal>(
366 p_i, p_j, alpha_i, alpha_j);
368 const Number alpha_hat_max = c(gamma_max) * alpha_max;
370 const Number gamma_m = std::min(gamma_i, gamma_j);
371 const Number gamma_M = std::max(gamma_i, gamma_j);
381 const Number r_exponent =
382 (gamma_M - gamma_min) / (
ScalarNumber(2.) * gamma_min * gamma_M);
391 const Number exponent =
393 const Number exponent_inverse = Number(1.) / exponent;
395 const Number numerator =
398 Number denominator = alpha_hat_min *
ryujin::pow(p_ratio, -exponent) +
403 const Number p_tilde = p_max *
ryujin::pow(temp, exponent_inverse) - pinf;
405#ifdef DEBUG_RIEMANN_SOLVER
406 std::cout <<
"IN p_*_tilde = " << p_tilde <<
"\n";
413 template <
int dim,
typename Number>
414 DEAL_II_ALWAYS_INLINE
inline Number
416 const Number p_star)
const
418 constexpr ScalarNumber
min = std::numeric_limits<ScalarNumber>::min();
420 const auto view = hyperbolic_system.view<dim, Number>();
421 const auto interpolation_b = view.eos_interpolation_b();
422 const auto pinf = view.eos_interpolation_pinfty();
424 const auto &[rho, u, p, gamma, a] = riemann_data;
426 const Number one_minus_b_rho = Number(1.) - interpolation_b * rho;
427 const Number gamma_minus_one = gamma - Number(1.);
430 ScalarNumber(2.) * one_minus_b_rho / (rho * (gamma + Number(1.)));
432 const Number Bz = gamma_minus_one / (gamma + Number(1.)) * (p + pinf);
434 const Number radicand =
safe_division(Az, p_star + pinf + Bz);
437 const Number true_value = (p_star - p) * std::sqrt(radicand);
439 const auto exponent = ScalarNumber(0.5) * gamma_minus_one / gamma;
442 const Number factor =
ryujin::pow(ratio, exponent) - Number(1.);
445 const auto false_value = ScalarNumber(2.) * a * one_minus_b_rho * factor /
446 std::max(gamma_minus_one, Number(min));
448 return dealii::compare_and_apply_mask<
449 dealii::SIMDComparison::greater_than_or_equal>(
450 p_star, p, true_value, false_value);
454 template <
int dim,
typename Number>
455 DEAL_II_ALWAYS_INLINE
inline Number
457 const primitive_type &riemann_data_j,
458 const Number p_in)
const
460 const Number &u_i = riemann_data_i[1];
461 const Number &u_j = riemann_data_j[1];
463 return f(riemann_data_i, p_in) + f(riemann_data_j, p_in) + u_j - u_i;
467 template <
int dim,
typename Number>
468 DEAL_II_ALWAYS_INLINE
inline Number
473 const auto view = hyperbolic_system.view<dim, Number>();
474 const auto interpolation_b = view.eos_interpolation_b();
475 const auto pinf = view.eos_interpolation_pinfty();
477 const auto &[rho_i, u_i, p_i, gamma_i, a_i] = riemann_data_i;
478 const auto &[rho_j, u_j, p_j, gamma_j, a_j] = riemann_data_j;
480 const Number p_max = std::max(p_i, p_j) + pinf;
482 const Number radicand_inverse_i =
484 Number(1.) - interpolation_b * rho_i) *
485 ((gamma_i + Number(1.)) * p_max +
486 (gamma_i - Number(1.)) * (p_i + pinf));
488 const Number value_i =
491 const Number radicand_inverse_j =
493 Number(1.) - interpolation_b * rho_j) *
494 ((gamma_j + Number(1.)) * p_max +
495 (gamma_j - Number(1.)) * (p_j + pinf));
497 const Number value_j =
500 return value_i + value_j + u_j - u_i;
504 template <
int dim,
typename Number>
505 DEAL_II_ALWAYS_INLINE
inline Number
509 const auto view = hyperbolic_system.view<dim, Number>();
510 const auto pinf = view.eos_interpolation_pinfty();
512 const auto &[rho, u, p, gamma, a] = riemann_data;
519 return u - a * std::sqrt(Number(1.) + factor * tmp);
523 template <
int dim,
typename Number>
524 DEAL_II_ALWAYS_INLINE
inline Number
526 const Number p_star)
const
528 const auto view = hyperbolic_system.view<dim, Number>();
529 const auto pinf = view.eos_interpolation_pinfty();
531 const auto &[rho, u, p, gamma, a] = riemann_data;
538 return u + a * std::sqrt(Number(1.) + factor * tmp);
542 template <
int dim,
typename Number>
543 DEAL_II_ALWAYS_INLINE
inline Number
547 const Number p_star)
const
549 const Number nu_11 = lambda1_minus(riemann_data_i, p_star);
550 const Number nu_32 = lambda3_plus(riemann_data_j, p_star);
556 template <
int dim,
typename Number>
557 DEAL_II_ALWAYS_INLINE
inline auto
561 const dealii::Tensor<1, dim, Number> &n_ij)
const ->
primitive_type
563 const auto view = hyperbolic_system.view<dim, Number>();
565 const auto rho = view.density(U);
568 const auto m = view.momentum(U);
569 const auto proj_m = n_ij * m;
571 const auto gamma = view.surrogate_gamma(U, p);
573 const auto interpolation_b = view.eos_interpolation_b();
574 const auto pinf = view.eos_interpolation_pinfty();
575 const auto x = Number(1.) - interpolation_b * rho;
576 const auto a = std::sqrt(gamma * (p + pinf) / (rho * x));
578#ifdef DEBUG_EXPENSIVE_BOUNDS_CHECK
582 dealii::ExcMessage(
"Internal error: p + pinf < 0."));
587 dealii::ExcMessage(
"Internal error: 1. - b * rho <= 0."));
592 dealii::ExcMessage(
"Internal error: gamma < 1."));
595 return {{rho, proj_m * rho_inverse, p, gamma, a}};
599 template <
int dim,
typename Number>
604 const auto view = hyperbolic_system.view<dim, Number>();
605 const auto pinf = view.eos_interpolation_pinfty();
607 const auto &[rho_i, u_i, p_i, gamma_i, a_i] = riemann_data_i;
608 const auto &[rho_j, u_j, p_j, gamma_j, a_j] = riemann_data_j;
610#ifdef DEBUG_RIEMANN_SOLVER
611 std::cout <<
"rho_left: " << rho_i << std::endl;
612 std::cout <<
"u_left: " << u_i << std::endl;
613 std::cout <<
"p_left: " << p_i << std::endl;
614 std::cout <<
"gamma_left: " << gamma_i << std::endl;
615 std::cout <<
"a_left: " << a_i << std::endl;
616 std::cout <<
"rho_right: " << rho_j << std::endl;
617 std::cout <<
"u_right: " << u_j << std::endl;
618 std::cout <<
"p_right: " << p_j << std::endl;
619 std::cout <<
"gamma_right: " << gamma_j << std::endl;
620 std::cout <<
"a_right: " << a_j << std::endl;
623 const Number p_max = std::max(p_i, p_j) + pinf;
624 const Number phi_p_max = phi_of_p_max(riemann_data_i, riemann_data_j);
626 if (!view.compute_strict_bounds()) {
627#ifdef DEBUG_RIEMANN_SOLVER
628 const Number p_star_RS = p_star_RS_full(riemann_data_i, riemann_data_j);
629 const Number p_star_SS = p_star_SS_full(riemann_data_i, riemann_data_j);
630 const Number p_debug =
631 dealii::compare_and_apply_mask<dealii::SIMDComparison::less_than>(
632 phi_p_max, Number(0.), p_star_SS, std::min(p_max, p_star_RS));
633 std::cout <<
" p^*_debug = " << p_debug <<
"\n";
634 std::cout <<
" phi(p_*_d) = "
635 << phi(riemann_data_i, riemann_data_j, p_debug) <<
"\n";
636 std::cout <<
"-> lambda_deb = "
637 << compute_lambda(riemann_data_i, riemann_data_j, p_debug)
641 const Number p_star_tilde =
642 p_star_interpolated(riemann_data_i, riemann_data_j);
643 const Number p_star_backup =
644 p_star_failsafe(riemann_data_i, riemann_data_j);
647 dealii::compare_and_apply_mask<dealii::SIMDComparison::less_than>(
650 std::min(p_star_tilde, p_star_backup),
651 std::min(p_max, p_star_tilde));
653#ifdef DEBUG_RIEMANN_SOLVER
654 std::cout <<
" p^*_tilde = " << p_2 <<
"\n";
655 std::cout <<
" phi(p_*_t) = "
656 << phi(riemann_data_i, riemann_data_j, p_2) <<
"\n";
657 std::cout <<
"-> lambda_max = "
658 << compute_lambda(riemann_data_i, riemann_data_j, p_2) <<
"\n"
662 return compute_lambda(riemann_data_i, riemann_data_j, p_2);
665 const Number p_star_RS = p_star_RS_full(riemann_data_i, riemann_data_j);
666 const Number p_star_SS = p_star_SS_full(riemann_data_i, riemann_data_j);
669 dealii::compare_and_apply_mask<dealii::SIMDComparison::less_than>(
670 phi_p_max, Number(0.), p_star_SS, std::min(p_max, p_star_RS));
672#ifdef DEBUG_RIEMANN_SOLVER
673 std::cout <<
" p^*_tilde = " << p_2 <<
"\n";
674 std::cout <<
" phi(p_*_t) = "
675 << phi(riemann_data_i, riemann_data_j, p_2) <<
"\n";
676 std::cout <<
"-> lambda_max = "
677 << compute_lambda(riemann_data_i, riemann_data_j, p_2)
681 return compute_lambda(riemann_data_i, riemann_data_j, p_2);
685 template <
int dim,
typename Number>
689 const unsigned int i,
690 const unsigned int *js,
691 const dealii::Tensor<1, dim, Number> &n_ij)
const
693 const auto &[p_i, unused_i, s_i, eta_i] =
694 precomputed_values.template get_tensor<Number, precomputed_type>(i);
696 const auto &[p_j, unused_j, s_j, eta_j] =
697 precomputed_values.template get_tensor<Number, precomputed_type>(js);
699 const auto riemann_data_i = riemann_data_from_state(U_i, p_i, n_ij);
700 const auto riemann_data_j = riemann_data_from_state(U_j, p_j, n_ij);
702 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
typename std::array< Number, riemann_data_size > primitive_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
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)
DEAL_II_ALWAYS_INLINE Number safe_division(const Number &numerator, const Number &denominator)