ryujin 2.1.1 revision ed06d5bf0c7d511bfbe1eea34219071b2f01d008
riemann_solver.template.h
Go to the documentation of this file.
1//
2// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
3// Copyright (C) 2020 - 2024 by the ryujin authors
4//
5
6#pragma once
7
8#include "riemann_solver.h"
9
10#include <newton.h>
11#include <simd.h>
12
13// #define DEBUG_RIEMANN_SOLVER
14
15namespace ryujin
16{
17 namespace EulerAEOS
18 {
19 /*
20 * The RiemannSolver is a guaranteed maximal wavespeed (GMS) estimate
21 * for the extended Riemann problem outlined in
22 * @cite ClaytonGuermondPopov-2022. For extenstions on handling negative
23 * pressures, we follow @cite clayton2023robust (see §4.6).
24 *
25 * In contrast to the algorithm outlined in above reference the
26 * algorithm takes a couple of shortcuts to significantly decrease the
27 * computational footprint. These simplifications still guarantee that
28 * we have an upper bound on the maximal wavespeed - but the number
29 * bound might be larger. In particular:
30 *
31 * - We do not check and treat the case phi(p_min) > 0. This
32 * corresponds to two expansion waves, see §5.2 in the reference. In
33 * this case we have
34 *
35 * 0 < p_star < p_min <= p_max.
36 *
37 * And due to the fact that p_star < p_min the wavespeeds reduce to
38 * a left wavespeed v_L - a_L and right wavespeed v_R + a_R. This
39 * implies that it is sufficient to set p_2 to ANY value provided
40 * that p_2 <= p_min hold true in order to compute the correct
41 * wavespeed.
42 *
43 * If p_2 > p_min then a more pessimistic bound is computed.
44 *
45 * - FIXME: Simplification in p_star_RS
46 */
47
48
49 template <int dim, typename Number>
50 DEAL_II_ALWAYS_INLINE inline Number
51 RiemannSolver<dim, Number>::c(const Number &gamma) const
52 {
53 /*
54 * We implement the continuous and monotonic function c(gamma) as
55 * defined in (A.3) on page A469 of @cite ClaytonGuermondPopov-2022.
56 * But with a simplified quick cut-off for the case gamma > 3:
57 *
58 * c(gamma)^2 = 1 for gamma <= 5 / 3
59 * c(gamma)^2 = (3 * gamma + 11) / (6 * gamma + 6) in between
60 * c(gamma)^2 = max(1/2, 5 / 6 - slope (gamma - 3)) for gamma > 3
61 *
62 * Due to the fact that the function is monotonic we can simply clip
63 * the values without checking the conditions:
64 */
65
66 constexpr ScalarNumber slope =
67 ScalarNumber(-0.34976871477801828189920753948709);
68
69 const Number first_radicand = (ScalarNumber(3.) * gamma + Number(11.)) /
70 (ScalarNumber(6.) * gamma + Number(6.));
71
72 const Number second_radicand =
73 Number(5. / 6.) + slope * (gamma - Number(3.));
74
75 Number radicand = std::min(first_radicand, second_radicand);
76 radicand = std::min(Number(1.), radicand);
77 radicand = std::max(Number(1. / 2.), radicand);
78
79 return std::sqrt(radicand);
80 }
81
82
83 template <int dim, typename Number>
84 DEAL_II_ALWAYS_INLINE inline Number RiemannSolver<dim, Number>::alpha(
85 const Number &rho, const Number &gamma, const Number &a) const
86 {
87 const auto view = hyperbolic_system.view<dim, Number>();
88 const auto interpolation_b = view.eos_interpolation_b();
89
90 const Number numerator =
91 ScalarNumber(2.) * a * (Number(1.) - interpolation_b * rho);
92
93 const Number denominator = gamma - Number(1.);
94
95 return safe_division(numerator, denominator);
96 }
97
98
99 template <int dim, typename Number>
100 DEAL_II_ALWAYS_INLINE inline Number
102 const primitive_type &riemann_data_i,
103 const primitive_type &riemann_data_j) const
104 {
105 const auto view = hyperbolic_system.view<dim, Number>();
106 const auto pinf = view.eos_interpolation_pinfty();
107
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);
112
113 /*
114 * First get p_min, p_max.
115 *
116 * Then, we get gamma_min/max, and alpha_min/max. Note that the
117 * *_min/max values are associated with p_min/max and are not
118 * necessarily the minimum/maximum of *_i vs *_j.
119 */
120
121 const Number p_min = std::min(p_i, p_j);
122 const Number p_max = std::max(p_i, p_j);
123
124 const Number gamma_min =
125 dealii::compare_and_apply_mask<dealii::SIMDComparison::less_than>(
126 p_i, p_j, gamma_i, gamma_j);
127
128 const Number alpha_min =
129 dealii::compare_and_apply_mask<dealii::SIMDComparison::less_than>(
130 p_i, p_j, alpha_i, alpha_j);
131
132 const Number alpha_hat_min = c(gamma_min) * alpha_min;
133
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);
137
138 const Number gamma_m = std::min(gamma_i, gamma_j);
139 const Number gamma_M = std::max(gamma_i, gamma_j);
140
141 const Number numerator =
142 dealii::compare_and_apply_mask<dealii::SIMDComparison::equal>(
143 p_max + pinf,
144 Number(0.),
145 Number(0.),
146 positive_part(alpha_hat_min + alpha_max - (u_j - u_i)));
147
148 /*
149 * The admissible set is p_min >= pinf. But numerically let's avoid
150 * division by zero and ensure positivity:
151 */
152 const Number p_ratio = safe_division(p_min + pinf, p_max + pinf);
153
154 /*
155 * Here, we use a trick: The r-factor only shows up in the formula
156 * for the case \gamma_min = \gamma_m, otherwise the r-factor
157 * vanishes. We can accomplish this by using the following modified
158 * exponent (where we substitute gamma_m by gamma_min):
159 */
160 const Number r_exponent =
161 (gamma_M - gamma_min) / (ScalarNumber(2.) * gamma_min * gamma_M);
162
163 /*
164 * Compute (5.7) first formula for \tilde p_1^\ast and (5.8)
165 * second formula for \tilde p_2^\ast at the same time:
166 */
167
168 const Number first_exponent =
169 (gamma_M - Number(1.)) / (ScalarNumber(2.) * gamma_M);
170
171 const Number first_exponent_inverse =
172 safe_division(Number(1.), first_exponent);
173
174 const Number first_denom =
175 alpha_hat_min * ryujin::pow(p_ratio, r_exponent - first_exponent) +
176 alpha_max;
177
178 const Number p_1_tilde =
179 (p_max + pinf) * ryujin::pow(safe_division(numerator, first_denom),
180 first_exponent_inverse) -
181 pinf;
182
183#ifdef DEBUG_RIEMANN_SOLVER
184 std::cout << "RS p_1_tilde = " << p_1_tilde << "\n";
185#endif
186
187 /*
188 * Compute (5.7) second formula for \tilde p_2^\ast and (5.8) first
189 * formula for \tilde p_1^\ast at the same time:
190 */
191
192 const Number second_exponent =
193 (gamma_m - Number(1.)) / (ScalarNumber(2.) * gamma_m);
194
195 const Number second_exponent_inverse =
196 safe_division(Number(1.), second_exponent);
197
198 Number second_denom =
199 alpha_hat_min * ryujin::pow(p_ratio, -second_exponent) +
200 alpha_max * ryujin::pow(p_ratio, r_exponent);
201
202 const Number p_2_tilde =
203 (p_max + pinf) * ryujin::pow(safe_division(numerator, second_denom),
204 second_exponent_inverse) -
205 pinf;
206
207#ifdef DEBUG_RIEMANN_SOLVER
208 std::cout << "RS p_2_tilde = " << p_2_tilde << "\n";
209#endif
210
211 return std::min(p_1_tilde, p_2_tilde);
212 }
213
214
215 template <int dim, typename Number>
216 DEAL_II_ALWAYS_INLINE inline Number
218 const primitive_type &riemann_data_i,
219 const primitive_type &riemann_data_j) const
220 {
221 const auto view = hyperbolic_system.view<dim, Number>();
222 const auto pinf = view.eos_interpolation_pinfty();
223
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;
226
227 const Number gamma_m = std::min(gamma_i, gamma_j);
228
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);
231
232 /*
233 * Compute (5.10) formula for \tilde p_1^\ast:
234 *
235 * Cost: 2x pow, 4x division, 0x sqrt
236 */
237
238 const Number exponent =
239 (gamma_m - Number(1.)) / (ScalarNumber(2.) * gamma_m);
240 const Number exponent_inverse = Number(1.) / exponent;
241
242 const Number numerator =
243 dealii::compare_and_apply_mask<dealii::SIMDComparison::equal>(
244 p_j + pinf,
245 Number(0.),
246 Number(0.),
247 positive_part(alpha_hat_i + alpha_hat_j - (u_j - u_i)));
248
249 const Number denominator =
250 alpha_hat_i *
251 ryujin::pow(safe_division(p_i + pinf, p_j + pinf), -exponent) +
252 alpha_hat_j;
253
254 const Number p_1_tilde =
255 (p_j + pinf) * ryujin::pow(safe_division(numerator, denominator),
256 exponent_inverse) -
257 pinf;
258
259#ifdef DEBUG_RIEMANN_SOLVER
260 std::cout << "SS p_1_tilde = " << p_1_tilde << "\n";
261#endif
262
263 const auto p_2_tilde = p_star_failsafe(riemann_data_i, riemann_data_j);
264
265 return std::min(p_1_tilde, p_2_tilde);
266 }
267
268
269 template <int dim, typename Number>
270 DEAL_II_ALWAYS_INLINE inline Number
272 const primitive_type &riemann_data_i,
273 const primitive_type &riemann_data_j) const
274 {
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();
278
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;
281
282 /*
283 * Compute (5.11) formula for \tilde p_2^\ast:
284 *
285 * Cost: 0x pow, 3x division, 3x sqrt
286 */
287
288 const Number p_max = std::max(p_i, p_j) + pinf;
289
290 const Number radicand_i = safe_division(
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)));
294
295 const Number x_i = std::sqrt(radicand_i);
296
297 const Number radicand_j = safe_division(
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)));
301
302 const Number x_j = std::sqrt(radicand_j);
303
304 const Number a = x_i + x_j;
305 const Number b =
306 dealii::compare_and_apply_mask<dealii::SIMDComparison::equal>(
307 a, Number(0.), Number(0.), u_j - u_i);
308
309 const Number c = -(p_i + pinf) * x_i - (p_j + pinf) * x_j;
310
311 const Number base = safe_division(
312 std::abs(-b +
313 std::sqrt(positive_part(b * b - ScalarNumber(4.) * a * c))),
314 std::abs(ScalarNumber(2.) * a));
315
316 const Number p_2_tilde = base * base - pinf;
317
318#ifdef DEBUG_RIEMANN_SOLVER
319 std::cout << "SS p_2_tilde = " << p_2_tilde << "\n";
320#endif
321 return p_2_tilde;
322 }
323
324
325 template <int dim, typename Number>
326 DEAL_II_ALWAYS_INLINE inline Number
328 const primitive_type &riemann_data_i,
329 const primitive_type &riemann_data_j) const
330 {
331 const auto view = hyperbolic_system.view<dim, Number>();
332 const auto pinf = view.eos_interpolation_pinfty();
333
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);
338
339 /*
340 * First get p_min, p_max.
341 *
342 * Then, we get gamma_min/max, and alpha_min/max. Note that the
343 * *_min/max values are associated with p_min/max and are not
344 * necessarily the minimum/maximum of *_i vs *_j.
345 */
346
347 const Number p_min = std::min(p_i, p_j) + pinf;
348 const Number p_max = std::max(p_i, p_j) + pinf;
349
350 const Number gamma_min =
351 dealii::compare_and_apply_mask<dealii::SIMDComparison::less_than>(
352 p_i, p_j, gamma_i, gamma_j);
353
354 const Number alpha_min =
355 dealii::compare_and_apply_mask<dealii::SIMDComparison::less_than>(
356 p_i, p_j, alpha_i, alpha_j);
357
358 const Number alpha_hat_min = c(gamma_min) * alpha_min;
359
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);
363
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);
367
368 const Number alpha_hat_max = c(gamma_max) * alpha_max;
369
370 const Number gamma_m = std::min(gamma_i, gamma_j);
371 const Number gamma_M = std::max(gamma_i, gamma_j);
372
373 const Number p_ratio = safe_division(p_min, p_max);
374
375 /*
376 * Here, we use a trick: The r-factor only shows up in the formula
377 * for the case \gamma_min = \gamma_m, otherwise the r-factor
378 * vanishes. We can accomplish this by using the following modified
379 * exponent (where we substitute gamma_m by gamma_min):
380 */
381 const Number r_exponent =
382 (gamma_M - gamma_min) / (ScalarNumber(2.) * gamma_min * gamma_M);
383
384 /*
385 * Compute a simultaneous upper bound on
386 * (5.7) second formula for \tilde p_2^\ast
387 * (5.8) first formula for \tilde p_1^\ast
388 * (5.11) formula for \tilde p_2^\ast
389 */
390
391 const Number exponent =
392 (gamma_m - Number(1.)) / (ScalarNumber(2.) * gamma_m);
393 const Number exponent_inverse = Number(1.) / exponent;
394
395 const Number numerator =
396 positive_part(alpha_hat_min + /*SIC!*/ alpha_max - (u_j - u_i));
397
398 Number denominator = alpha_hat_min * ryujin::pow(p_ratio, -exponent) +
399 alpha_hat_max * ryujin::pow(p_ratio, r_exponent);
400
401 const auto temp = safe_division(numerator, denominator);
402
403 const Number p_tilde = p_max * ryujin::pow(temp, exponent_inverse) - pinf;
404
405#ifdef DEBUG_RIEMANN_SOLVER
406 std::cout << "IN p_*_tilde = " << p_tilde << "\n";
407#endif
408
409 return p_tilde;
410 }
411
412
413 template <int dim, typename Number>
414 DEAL_II_ALWAYS_INLINE inline Number
415 RiemannSolver<dim, Number>::f(const primitive_type &riemann_data,
416 const Number p_star) const
417 {
418 constexpr ScalarNumber min = std::numeric_limits<ScalarNumber>::min();
419
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();
423
424 const auto &[rho, u, p, gamma, a] = riemann_data;
425
426 const Number one_minus_b_rho = Number(1.) - interpolation_b * rho;
427 const Number gamma_minus_one = gamma - Number(1.);
428
429 const Number Az =
430 ScalarNumber(2.) * one_minus_b_rho / (rho * (gamma + Number(1.)));
431
432 const Number Bz = gamma_minus_one / (gamma + Number(1.)) * (p + pinf);
433
434 const Number radicand = safe_division(Az, p_star + pinf + Bz);
435
436 /* true_value is shock case */
437 const Number true_value = (p_star - p) * std::sqrt(radicand);
438
439 const auto exponent = ScalarNumber(0.5) * gamma_minus_one / gamma;
440
441 const Number ratio = safe_division(p_star + pinf, p + pinf);
442 const Number factor = ryujin::pow(ratio, exponent) - Number(1.);
443
444 /* false_value is rarefaction case */
445 const auto false_value = ScalarNumber(2.) * a * one_minus_b_rho * factor /
446 std::max(gamma_minus_one, Number(min));
447
448 return dealii::compare_and_apply_mask<
449 dealii::SIMDComparison::greater_than_or_equal>(
450 p_star, p, true_value, false_value);
451 }
452
453
454 template <int dim, typename Number>
455 DEAL_II_ALWAYS_INLINE inline Number
456 RiemannSolver<dim, Number>::phi(const primitive_type &riemann_data_i,
457 const primitive_type &riemann_data_j,
458 const Number p_in) const
459 {
460 const Number &u_i = riemann_data_i[1];
461 const Number &u_j = riemann_data_j[1];
462
463 return f(riemann_data_i, p_in) + f(riemann_data_j, p_in) + u_j - u_i;
464 }
465
466
467 template <int dim, typename Number>
468 DEAL_II_ALWAYS_INLINE inline Number
470 const primitive_type &riemann_data_i,
471 const primitive_type &riemann_data_j) const
472 {
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();
476
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;
479
480 const Number p_max = std::max(p_i, p_j) + pinf;
481
482 const Number radicand_inverse_i =
483 safe_division(ScalarNumber(0.5) * rho_i,
484 Number(1.) - interpolation_b * rho_i) *
485 ((gamma_i + Number(1.)) * p_max +
486 (gamma_i - Number(1.)) * (p_i + pinf));
487
488 const Number value_i =
489 safe_division(p_max - p_i, std::sqrt(radicand_inverse_i));
490
491 const Number radicand_inverse_j =
492 safe_division(ScalarNumber(0.5) * rho_j,
493 Number(1.) - interpolation_b * rho_j) *
494 ((gamma_j + Number(1.)) * p_max +
495 (gamma_j - Number(1.)) * (p_j + pinf));
496
497 const Number value_j =
498 safe_division(p_max - p_j, std::sqrt(radicand_inverse_j));
499
500 return value_i + value_j + u_j - u_i;
501 }
502
503
504 template <int dim, typename Number>
505 DEAL_II_ALWAYS_INLINE inline Number
507 const primitive_type &riemann_data, const Number p_star) const
508 {
509 const auto view = hyperbolic_system.view<dim, Number>();
510 const auto pinf = view.eos_interpolation_pinfty();
511
512 const auto &[rho, u, p, gamma, a] = riemann_data;
513
514 const auto factor =
515 ScalarNumber(0.5) * (gamma + ScalarNumber(1.)) / gamma;
516
517 const Number tmp = safe_division(positive_part(p_star - p), p + pinf);
518
519 return u - a * std::sqrt(Number(1.) + factor * tmp);
520 }
521
522
523 template <int dim, typename Number>
524 DEAL_II_ALWAYS_INLINE inline Number
526 const Number p_star) const
527 {
528 const auto view = hyperbolic_system.view<dim, Number>();
529 const auto pinf = view.eos_interpolation_pinfty();
530
531 const auto &[rho, u, p, gamma, a] = riemann_data;
532
533 const auto factor =
534 ScalarNumber(0.5) * (gamma + ScalarNumber(1.)) / gamma;
535
536 const Number tmp = safe_division(positive_part(p_star - p), p + pinf);
537
538 return u + a * std::sqrt(Number(1.) + factor * tmp);
539 }
540
541
542 template <int dim, typename Number>
543 DEAL_II_ALWAYS_INLINE inline Number
545 const primitive_type &riemann_data_i,
546 const primitive_type &riemann_data_j,
547 const Number p_star) const
548 {
549 const Number nu_11 = lambda1_minus(riemann_data_i, p_star);
550 const Number nu_32 = lambda3_plus(riemann_data_j, p_star);
551
552 return std::max(positive_part(nu_32), negative_part(nu_11));
553 }
554
555
556 template <int dim, typename Number>
557 DEAL_II_ALWAYS_INLINE inline auto
559 const state_type &U,
560 const Number &p,
561 const dealii::Tensor<1, dim, Number> &n_ij) const -> primitive_type
562 {
563 const auto view = hyperbolic_system.view<dim, Number>();
564
565 const auto rho = view.density(U);
566 const auto rho_inverse = ScalarNumber(1.0) / rho;
567
568 const auto m = view.momentum(U);
569 const auto proj_m = n_ij * m;
570
571 const auto gamma = view.surrogate_gamma(U, p);
572
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));
577
578#ifdef DEBUG_EXPENSIVE_BOUNDS_CHECK
580 Number(p + pinf),
581 [](auto val) { return val >= ScalarNumber(0.); },
582 dealii::ExcMessage("Internal error: p + pinf < 0."));
583
585 x,
586 [](auto val) { return val > ScalarNumber(0.); },
587 dealii::ExcMessage("Internal error: 1. - b * rho <= 0."));
588
590 gamma,
591 [](auto val) { return val >= ScalarNumber(1.); },
592 dealii::ExcMessage("Internal error: gamma < 1."));
593#endif
594
595 return {{rho, proj_m * rho_inverse, p, gamma, a}};
596 }
597
598
599 template <int dim, typename Number>
601 const primitive_type &riemann_data_i,
602 const primitive_type &riemann_data_j) const
603 {
604 const auto view = hyperbolic_system.view<dim, Number>();
605 const auto pinf = view.eos_interpolation_pinfty();
606
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;
609
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;
621#endif
622
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);
625
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)
638 << std::endl;
639#endif
640
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);
645
646 const Number p_2 =
647 dealii::compare_and_apply_mask<dealii::SIMDComparison::less_than>(
648 phi_p_max,
649 Number(0.),
650 std::min(p_star_tilde, p_star_backup),
651 std::min(p_max, p_star_tilde));
652
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"
659 << std::endl;
660#endif
661
662 return compute_lambda(riemann_data_i, riemann_data_j, p_2);
663 }
664
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);
667
668 const Number p_2 =
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));
671
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)
678 << std::endl;
679#endif
680
681 return compute_lambda(riemann_data_i, riemann_data_j, p_2);
682 }
683
684
685 template <int dim, typename Number>
686 DEAL_II_ALWAYS_INLINE inline Number RiemannSolver<dim, Number>::compute(
687 const state_type &U_i,
688 const state_type &U_j,
689 const unsigned int i,
690 const unsigned int *js,
691 const dealii::Tensor<1, dim, Number> &n_ij) const
692 {
693 const auto &[p_i, unused_i, s_i, eta_i] =
694 precomputed_values.template get_tensor<Number, precomputed_type>(i);
695
696 const auto &[p_j, unused_j, s_j, eta_j] =
697 precomputed_values.template get_tensor<Number, precomputed_type>(js);
698
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);
701
702 return compute(riemann_data_i, riemann_data_j);
703 }
704
705
706 } // namespace EulerAEOS
707} // namespace ryujin
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)
Definition: simd.h:124
DEAL_II_ALWAYS_INLINE Number positive_part(const Number number)
Definition: simd.h:112
DEAL_II_ALWAYS_INLINE Number safe_division(const Number &numerator, const Number &denominator)