ryujin 2.1.1 revision 9dcb748690310d6a540ebb8b066d1a0834fc7604
limiter.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 "limiter.h"
9
10namespace ryujin
11{
12 namespace ScalarConservation
13 {
14 template <int dim, typename Number>
15 std::tuple<Number, bool>
17 const state_type &U,
18 const state_type &P,
19 const Number t_min /* = Number(0.) */,
20 const Number t_max /* = Number(1.) */)
21 {
22 bool success = true;
23 Number t_r = t_max;
24
25 constexpr ScalarNumber eps = std::numeric_limits<ScalarNumber>::epsilon();
26 const ScalarNumber relax = ScalarNumber(1. + 10000. * eps);
27
28 const auto &u_U = hyperbolic_system.state(U);
29 const auto &u_P = hyperbolic_system.state(P);
30
31 const auto &u_min = std::get<0>(bounds);
32 const auto &u_max = std::get<1>(bounds);
33
34 /*
35 * Verify that u_U is within bounds. This property might be
36 * violated for relative CFL numbers larger than 1.
37 */
38 const auto test_min = std::max(Number(0.), u_U - relax * u_max);
39 const auto test_max = std::max(Number(0.), u_min - relax * u_U);
40 if (!(test_min == Number(0.) && test_max == Number(0.))) {
41#ifdef DEBUG_OUTPUT
42 std::cout << std::fixed << std::setprecision(16);
43 std::cout << "Bounds violation: low-order state (critical)!"
44 << "\n\t\tu min: " << u_min
45 << "\n\t\tu min (delta): " << negative_part(u_U - u_min)
46 << "\n\t\tu: " << u_U
47 << "\n\t\tu max (delta): " << positive_part(u_U - u_max)
48 << "\n\t\tu max: " << u_max << "\n"
49 << std::endl;
50#endif
51 success = false;
52 }
53
54 const Number denominator =
55 ScalarNumber(1.) / (std::abs(u_P) + eps * u_max);
56
57 t_r = dealii::compare_and_apply_mask<dealii::SIMDComparison::less_than>(
58 u_max,
59 u_U + t_r * u_P,
60 /*
61 * u_P is positive.
62 *
63 * Note: Do not take an absolute value here. If we are out of
64 * bounds we have to ensure that t_r is set to t_min.
65 */
66 (u_max - u_U) * denominator,
67 t_r);
68
69 t_r = dealii::compare_and_apply_mask<dealii::SIMDComparison::less_than>(
70 u_U + t_r * u_P,
71 u_min,
72 /*
73 * u_P is negative.
74 *
75 * Note: Do not take an absolute value here. If we are out of
76 * bounds we have to ensure that t_r is set to t_min.
77 */
78 (u_U - u_min) * denominator,
79 t_r);
80
81 /*
82 * Ensure that t_min <= t <= t_max. This might not be the case if
83 * u_U is outside the interval [u_min, u_max]. Furthermore,
84 * the quotient we take above is prone to numerical cancellation in
85 * particular in the second pass of the limiter when u_P might be
86 * small.
87 */
88 t_r = std::min(t_r, t_max);
89 t_r = std::max(t_r, t_min);
90
91#ifdef CHECK_BOUNDS
92 /*
93 * Verify that the new state is within bounds:
94 */
95 const auto u_new = hyperbolic_system.state(U + t_r * P);
96 const auto test_new_min = std::max(Number(0.), u_new - relax * u_max);
97 const auto test_new_max = std::max(Number(0.), u_min - relax * u_new);
98 if (!(test_new_min == Number(0.) && test_new_max == Number(0.))) {
99#ifdef DEBUG_OUTPUT
100 std::cout << std::fixed << std::setprecision(16);
101 std::cout << "Bounds violation: high-order state!"
102 << "\n\t\tu min: " << u_min
103 << "\n\t\tu min (delta): " << negative_part(u_new - u_min)
104 << "\n\t\tu: " << u_new
105 << "\n\t\tu max (delta): " << positive_part(u_new - u_max)
106 << "\n\t\tu max: " << u_max << "\n"
107 << std::endl;
108#endif
109 success = false;
110 }
111#endif
112
113 return {t_r, success};
114 }
115
116 } // namespace ScalarConservation
117} // namespace ryujin
typename HyperbolicSystemView::state_type state_type
Definition: limiter.h:36
std::array< Number, n_bounds > Bounds
Definition: limiter.h:82
std::tuple< Number, bool > limit(const Bounds &bounds, const state_type &U, const state_type &P, const Number t_min=Number(0.), const Number t_max=Number(1.))
typename get_value_type< Number >::type ScalarNumber
Definition: limiter.h:53
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