ryujin 2.1.1 revision fd1a8e8ebb4535655802e400781bce67403e1e49
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 <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 EulerAEOS
20 {
21 /*
22 * The RiemannSolver is a guaranteed maximal wavespeed (GMS) estimate
23 * for the extended Riemann problem outlined in
24 * @cite ClaytonGuermondPopov-2022.
25 *
26 * In contrast to the algorithm outlined in above reference the
27 * algorithm takes a couple of shortcuts to significantly decrease the
28 * computational footprint. These simplifications still guarantee that
29 * we have an upper bound on the maximal wavespeed - but the number
30 * bound might be larger. In particular:
31 *
32 * - We do not check and treat the case phi(p_min) > 0. This
33 * corresponds to two expansion waves, see ยง5.2 in the reference. In
34 * this case we have
35 *
36 * 0 < p_star < p_min <= p_max.
37 *
38 * And due to the fact that p_star < p_min the wavespeeds reduce to
39 * a left wavespeed v_L - a_L and right wavespeed v_R + a_R. This
40 * implies that it is sufficient to set p_2 to ANY value provided
41 * that p_2 <= p_min hold true in order to compute the correct
42 * wavespeed.
43 *
44 * If p_2 > p_min then a more pessimistic bound is computed.
45 *
46 * - FIXME: Simplification in p_star_RS
47 */
48
49
50 template <int dim, typename Number>
51 DEAL_II_ALWAYS_INLINE inline Number
52 RiemannSolver<dim, Number>::c(const Number gamma) const
53 {
54 /*
55 * We implement the continuous and monotonic function c(gamma) as
56 * defined in (A.3) on page A469 of @cite ClaytonGuermondPopov-2022.
57 * But with a simplified quick cut-off for the case gamma > 3:
58 *
59 * c(gamma)^2 = 1 for gamma <= 5 / 3
60 * c(gamma)^2 = (3 * gamma + 11) / (6 * gamma + 6) in between
61 * c(gamma)^2 = max(1/2, 5 / 6 - slope (gamma - 3)) for gamma > 3
62 *
63 * Due to the fact that the function is monotonic we can simply clip
64 * the values without checking the conditions:
65 */
66
67 constexpr ScalarNumber slope =
68 ScalarNumber(-0.34976871477801828189920753948709);
69
70 const Number first_radicand = (ScalarNumber(3.) * gamma + Number(11.)) /
71 (ScalarNumber(6.) * gamma + Number(6.));
72
73 const Number second_radicand =
74 Number(5. / 6.) + slope * (gamma - Number(3.));
75
76 Number radicand = std::min(first_radicand, second_radicand);
77 radicand = std::min(Number(1.), radicand);
78 radicand = std::max(Number(1. / 2.), radicand);
79
80 return std::sqrt(radicand);
81 }
82
83
84 template <int dim, typename Number>
85 DEAL_II_ALWAYS_INLINE inline Number RiemannSolver<dim, Number>::alpha(
86 const Number &rho, const Number &gamma, const Number &a) const
87 {
88 const auto view = hyperbolic_system.view<dim, Number>();
89 const auto interpolation_b = view.eos_interpolation_b();
90
91 const Number numerator =
92 ScalarNumber(2.) * a * (Number(1.) - interpolation_b * rho);
93
94 const Number denominator = gamma - Number(1.);
95
96 return numerator / denominator;
97 }
98
99
100 template <int dim, typename Number>
101 DEAL_II_ALWAYS_INLINE inline Number
103 const primitive_type &riemann_data_i,
104 const primitive_type &riemann_data_j) const
105 {
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);
110
111 /*
112 * First get p_min, p_max.
113 *
114 * Then, we get gamma_min/max, and alpha_min/max. Note that the
115 * *_min/max values are associated with p_min/max and are not
116 * necessarily the minimum/maximum of *_i vs *_j.
117 */
118
119 const Number p_min = std::min(p_i, p_j);
120 const Number p_max = std::max(p_i, p_j);
121
122 const Number gamma_min =
123 dealii::compare_and_apply_mask<dealii::SIMDComparison::less_than>(
124 p_i, p_j, gamma_i, gamma_j);
125
126 const Number alpha_min =
127 dealii::compare_and_apply_mask<dealii::SIMDComparison::less_than>(
128 p_i, p_j, alpha_i, alpha_j);
129
130 const Number alpha_hat_min = c(gamma_min) * alpha_min;
131
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);
135
136 const Number gamma_m = std::min(gamma_i, gamma_j);
137 const Number gamma_M = std::max(gamma_i, gamma_j);
138
139 const Number numerator =
140 positive_part(alpha_hat_min + alpha_max - (u_j - u_i));
141
142 const Number p_ratio = p_min / p_max;
143
144 /*
145 * Here, we use a trick: The r-factor only shows up in the formula
146 * for the case \gamma_min = \gamma_m, otherwise the r-factor
147 * vanishes. We can accomplish this by using the following modified
148 * exponent (where we substitute gamma_m by gamma_min):
149 */
150 const Number r_exponent =
151 (gamma_M - gamma_min) / (ScalarNumber(2.) * gamma_min * gamma_M);
152
153 /*
154 * Compute (5.7) first formula for \tilde p_1^\ast and (5.8)
155 * second formula for \tilde p_2^\ast at the same time:
156 */
157
158 const Number first_exponent =
159 (gamma_M - Number(1.)) / (ScalarNumber(2.) * gamma_M);
160 const Number first_exponent_inverse = Number(1.) / first_exponent;
161
162 const Number first_denom =
163 alpha_hat_min * ryujin::pow(p_ratio, r_exponent - first_exponent) +
164 alpha_max;
165
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";
170#endif
171
172 /*
173 * Compute (5.7) second formula for \tilde p_2^\ast and (5.8) first
174 * formula for \tilde p_1^\ast at the same time:
175 */
176
177 const Number second_exponent =
178 (gamma_m - Number(1.)) / (ScalarNumber(2.) * gamma_m);
179 const Number second_exponent_inverse = Number(1.) / second_exponent;
180
181 Number second_denom =
182 alpha_hat_min * ryujin::pow(p_ratio, -second_exponent) +
183 alpha_max * ryujin::pow(p_ratio, r_exponent);
184
185 const Number p_2_tilde = p_max * ryujin::pow(numerator / second_denom,
186 second_exponent_inverse);
187
188#ifdef DEBUG_RIEMANN_SOLVER
189 std::cout << "RS p_2_tilde = " << p_2_tilde << "\n";
190#endif
191
192 return std::min(p_1_tilde, p_2_tilde);
193 }
194
195
196 template <int dim, typename Number>
197 DEAL_II_ALWAYS_INLINE inline Number
199 const primitive_type &riemann_data_i,
200 const primitive_type &riemann_data_j) const
201 {
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;
204
205 const Number gamma_m = std::min(gamma_i, gamma_j);
206
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);
209
210 /*
211 * Compute (5.10) formula for \tilde p_1^\ast:
212 *
213 * Cost: 2x pow, 4x division, 0x sqrt
214 */
215
216 const Number exponent =
217 (gamma_m - Number(1.)) / (ScalarNumber(2.) * gamma_m);
218 const Number exponent_inverse = Number(1.) / exponent;
219
220 const Number numerator =
221 positive_part(alpha_hat_i + alpha_hat_j - (u_j - u_i));
222
223 const Number denominator =
224 alpha_hat_i * ryujin::pow(p_i / p_j, -exponent) + alpha_hat_j;
225
226 const Number p_1_tilde =
227 p_j * ryujin::pow(numerator / denominator, exponent_inverse);
228
229#ifdef DEBUG_RIEMANN_SOLVER
230 std::cout << "SS p_1_tilde = " << p_1_tilde << "\n";
231#endif
232
233 const auto p_2_tilde = p_star_failsafe(riemann_data_i, riemann_data_j);
234
235 return std::min(p_1_tilde, p_2_tilde);
236 }
237
238
239 template <int dim, typename Number>
240 DEAL_II_ALWAYS_INLINE inline Number
242 const primitive_type &riemann_data_i,
243 const primitive_type &riemann_data_j) const
244 {
245 const auto view = hyperbolic_system.view<dim, Number>();
246 const auto interpolation_b = view.eos_interpolation_b();
247
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;
250
251 /*
252 * Compute (5.11) formula for \tilde p_2^\ast:
253 *
254 * Cost: 0x pow, 3x division, 3x sqrt
255 */
256
257 const Number p_max = std::max(p_i, p_j);
258
259 Number radicand_i =
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);
263
264 const Number x_i = std::sqrt(radicand_i);
265
266 Number radicand_j =
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);
270
271 const Number x_j = std::sqrt(radicand_j);
272
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;
276
277 const Number base = (-b + std::sqrt(b * b - ScalarNumber(4.) * a * c)) /
278 (ScalarNumber(2.) * a);
279 const Number p_2_tilde = base * base;
280
281#ifdef DEBUG_RIEMANN_SOLVER
282 std::cout << "SS p_2_tilde = " << p_2_tilde << "\n";
283#endif
284 return p_2_tilde;
285 }
286
287
288 template <int dim, typename Number>
289 DEAL_II_ALWAYS_INLINE inline Number
291 const primitive_type &riemann_data_i,
292 const primitive_type &riemann_data_j) const
293 {
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);
298
299 /*
300 * First get p_min, p_max.
301 *
302 * Then, we get gamma_min/max, and alpha_min/max. Note that the
303 * *_min/max values are associated with p_min/max and are not
304 * necessarily the minimum/maximum of *_i vs *_j.
305 */
306
307 const Number p_min = std::min(p_i, p_j);
308 const Number p_max = std::max(p_i, p_j);
309
310 const Number gamma_min =
311 dealii::compare_and_apply_mask<dealii::SIMDComparison::less_than>(
312 p_i, p_j, gamma_i, gamma_j);
313
314 const Number alpha_min =
315 dealii::compare_and_apply_mask<dealii::SIMDComparison::less_than>(
316 p_i, p_j, alpha_i, alpha_j);
317
318 const Number alpha_hat_min = c(gamma_min) * alpha_min;
319
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);
323
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);
327
328 const Number alpha_hat_max = c(gamma_max) * alpha_max;
329
330 const Number gamma_m = std::min(gamma_i, gamma_j);
331 const Number gamma_M = std::max(gamma_i, gamma_j);
332
333 const Number p_ratio = p_min / p_max;
334
335 /*
336 * Here, we use a trick: The r-factor only shows up in the formula
337 * for the case \gamma_min = \gamma_m, otherwise the r-factor
338 * vanishes. We can accomplish this by using the following modified
339 * exponent (where we substitute gamma_m by gamma_min):
340 */
341 const Number r_exponent =
342 (gamma_M - gamma_min) / (ScalarNumber(2.) * gamma_min * gamma_M);
343
344 /*
345 * Compute a simultaneous upper bound on
346 * (5.7) second formula for \tilde p_2^\ast
347 * (5.8) first formula for \tilde p_1^\ast
348 * (5.11) formula for \tilde p_2^\ast
349 */
350
351 const Number exponent =
352 (gamma_m - Number(1.)) / (ScalarNumber(2.) * gamma_m);
353 const Number exponent_inverse = Number(1.) / exponent;
354
355 const Number numerator =
356 positive_part(alpha_hat_min + /*SIC!*/ alpha_max - (u_j - u_i));
357
358 Number denominator = alpha_hat_min * ryujin::pow(p_ratio, -exponent) +
359 alpha_hat_max * ryujin::pow(p_ratio, r_exponent);
360
361 const Number p_tilde =
362 p_max * ryujin::pow(numerator / denominator, exponent_inverse);
363
364#ifdef DEBUG_RIEMANN_SOLVER
365 std::cout << "IN p_*_tilde = " << p_tilde << "\n";
366#endif
367
368 return p_tilde;
369 }
370
371
372 template <int dim, typename Number>
373 DEAL_II_ALWAYS_INLINE inline Number
374 RiemannSolver<dim, Number>::f(const primitive_type &riemann_data,
375 const Number p_star) const
376 {
377 const auto view = hyperbolic_system.view<dim, Number>();
378 const auto interpolation_b = view.eos_interpolation_b();
379
380 const auto &[rho, u, p, gamma, a] = riemann_data;
381
382 const Number one_minus_b_rho = Number(1.) - interpolation_b * rho;
383
384 const Number Az =
385 ScalarNumber(2.) * one_minus_b_rho / (rho * (gamma + Number(1.)));
386
387 const Number Bz = (gamma - Number(1.)) / (gamma + Number(1.)) * p;
388
389 const Number radicand = Az / (p_star + Bz);
390
391 const Number true_value = (p_star - p) * std::sqrt(radicand);
392
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.));
397
398 return dealii::compare_and_apply_mask<
399 dealii::SIMDComparison::greater_than_or_equal>(
400 p_star, p, true_value, false_value);
401 }
402
403
404 template <int dim, typename Number>
405 DEAL_II_ALWAYS_INLINE inline Number
406 RiemannSolver<dim, Number>::phi(const primitive_type &riemann_data_i,
407 const primitive_type &riemann_data_j,
408 const Number p_in) const
409 {
410 const Number &u_i = riemann_data_i[1];
411 const Number &u_j = riemann_data_j[1];
412
413 return f(riemann_data_i, p_in) + f(riemann_data_j, p_in) + u_j - u_i;
414 }
415
416
417 template <int dim, typename Number>
418 DEAL_II_ALWAYS_INLINE inline Number
420 const primitive_type &riemann_data_i,
421 const primitive_type &riemann_data_j) const
422 {
423 const auto view = hyperbolic_system.view<dim, Number>();
424 const auto interpolation_b = view.eos_interpolation_b();
425
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;
428
429 const Number p_max = std::max(p_i, p_j);
430
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);
434
435 const Number value_i = (p_max - p_i) / std::sqrt(radicand_inverse_i);
436
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);
440
441 const Number value_j = (p_max - p_j) / std::sqrt(radicand_jnverse_j);
442
443 return value_i + value_j + u_j - u_i;
444 }
445
446
447 template <int dim, typename Number>
448 DEAL_II_ALWAYS_INLINE inline Number
450 const primitive_type &riemann_data, const Number p_star) const
451 {
452 const auto &[rho, u, p, gamma, a] = riemann_data;
453
454 const auto factor =
455 ScalarNumber(0.5) * (gamma + ScalarNumber(1.)) / gamma;
456
457 const Number tmp = positive_part((p_star - p) / p);
458
459 return u - a * std::sqrt(Number(1.) + factor * tmp);
460 }
461
462
463 template <int dim, typename Number>
464 DEAL_II_ALWAYS_INLINE inline Number
466 const Number p_star) const
467 {
468 const auto &[rho, u, p, gamma, a] = riemann_data;
469
470 const auto factor =
471 ScalarNumber(0.5) * (gamma + ScalarNumber(1.)) / gamma;
472
473 const Number tmp = positive_part((p_star - p) / p);
474
475 return u + a * std::sqrt(Number(1.) + factor * tmp);
476 }
477
478
479 template <int dim, typename Number>
480 DEAL_II_ALWAYS_INLINE inline Number
482 const primitive_type &riemann_data_i,
483 const primitive_type &riemann_data_j,
484 const Number p_star) const
485 {
486 const Number nu_11 = lambda1_minus(riemann_data_i, p_star);
487 const Number nu_32 = lambda3_plus(riemann_data_j, p_star);
488
489 return std::max(positive_part(nu_32), negative_part(nu_11));
490 }
491
492
493 template <int dim, typename Number>
494 DEAL_II_ALWAYS_INLINE inline auto
496 const state_type &U,
497 const Number &p,
498 const dealii::Tensor<1, dim, Number> &n_ij) const -> primitive_type
499 {
500 const auto view = hyperbolic_system.view<dim, Number>();
501
502 const auto rho = view.density(U);
503 const auto rho_inverse = ScalarNumber(1.0) / rho;
504
505 const auto m = view.momentum(U);
506 const auto proj_m = n_ij * m;
507
508 const auto gamma = view.surrogate_gamma(U, p);
509
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));
513
514#ifdef CHECK_BOUNDS
516 Number(p),
517 [](auto val) { return val > ScalarNumber(0.); },
518 dealii::ExcMessage("Internal error: p <= 0."));
519
521 x,
522 [](auto val) { return val > ScalarNumber(0.); },
523 dealii::ExcMessage("Internal error: 1. - b * rho <= 0."));
524
526 gamma,
527 [](auto val) { return val > ScalarNumber(1.); },
528 dealii::ExcMessage("Internal error: gamma <= 1."));
529#endif
530
531 return {{rho, proj_m * rho_inverse, p, gamma, a}};
532 }
533
534
535 template <int dim, typename Number>
537 const primitive_type &riemann_data_i,
538 const primitive_type &riemann_data_j) const
539 {
540 const auto view = hyperbolic_system.view<dim, Number>();
541
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;
544
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;
556#endif
557
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);
560
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)
573 << std::endl;
574#endif
575
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);
580
581 const Number p_2 =
582 dealii::compare_and_apply_mask<dealii::SIMDComparison::less_than>(
583 phi_p_max,
584 Number(0.),
585 std::min(p_star_tilde, p_star_backup),
586 std::min(p_max, p_star_tilde));
587
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"
594 << std::endl;
595#endif
596
597 return compute_lambda(riemann_data_i, riemann_data_j, p_2);
598 }
599
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);
602
603 const Number p_2 =
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));
606
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)
613 << std::endl;
614#endif
615
616 return compute_lambda(riemann_data_i, riemann_data_j, p_2);
617 }
618
619
620 template <int dim, typename Number>
621 DEAL_II_ALWAYS_INLINE inline Number RiemannSolver<dim, Number>::compute(
622 const state_type &U_i,
623 const state_type &U_j,
624 const unsigned int i,
625 const unsigned int *js,
626 const dealii::Tensor<1, dim, Number> &n_ij) const
627 {
628 const auto &[p_i, unused_i, s_i, eta_i] =
629 precomputed_values.template get_tensor<Number, precomputed_type>(i);
630
631 const auto &[p_j, unused_j, s_j, eta_j] =
632 precomputed_values.template get_tensor<Number, precomputed_type>(js);
633
634 const auto riemann_data_i = riemann_data_from_state(U_i, p_i, n_ij);
635 const auto riemann_data_j = riemann_data_from_state(U_j, p_j, n_ij);
636
637 return compute(riemann_data_i, riemann_data_j);
638 }
639
640
641 } // namespace EulerAEOS
642} // 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