Line data Source code
1 : #ifndef SOLVER_LOCAL_RIEMANN_HLLE_H
2 : #define SOLVER_LOCAL_RIEMANN_HLLE_H
3 :
4 : #include "IO/ParmParse.H"
5 : #include "Solver/Local/Riemann/Riemann.H"
6 :
7 : /// A bunch of solvers
8 : namespace Solver
9 : {
10 : /// Local solvers
11 : namespace Local
12 : {
13 :
14 : namespace Riemann
15 : {
16 :
17 : class HLLE : public Riemann
18 : {
19 : public:
20 :
21 :
22 : static constexpr const char* name = "hlle";
23 1 : HLLE (IO::ParmParse &pp, std::string name)
24 1 : {pp_queryclass(name,*this);}
25 : HLLE (IO::ParmParse &pp)
26 : {pp_queryclass(*this);}
27 : HLLE ()
28 : {
29 : IO::ParmParse pp;
30 : pp_queryclass(*this);
31 : }
32 :
33 : int lowmach = 0;
34 : Set::Scalar cutoff = NAN;
35 :
36 1 : static void Parse(HLLE & value, IO::ParmParse & pp)
37 : {
38 : // toggle to apply a low mach correction approximation
39 2 : pp.query_default("lowmach",value.lowmach,false);
40 : // cutoff value if using the low mach approximation
41 2 : pp.query_default("cutoff",value.cutoff,0.1);
42 1 : }
43 :
44 8192000 : virtual Flux Solve(State lo, State hi, Model::Gas::Gas& gas, Set::Patch<const Set::Scalar>& X, int i, int j, int k, int side, Set::Scalar small) override
45 : {
46 : // Grab and label the input states
47 8192000 : Set::Scalar rho_L = lo.rho , rho_R = hi.rho;
48 8192000 : Set::Scalar Mn_L = lo.M_normal , Mn_R = hi.M_normal ;
49 : //Set::Scalar Mt_L = lo.M_tangent , Mt_R = hi.M_tangent ;
50 8192000 : Set::Scalar E_L = lo.E , E_R = hi.E ;
51 :
52 : // Ensure no negative densities
53 8192000 : rho_L = std::max(small,rho_L);
54 8192000 : rho_R = std::max(small,rho_R);
55 :
56 :
57 : // Calculate primitives
58 8192000 : Set::Scalar u_L = Mn_L/(rho_L+small), u_R = Mn_R/(rho_R+small);
59 :
60 8192000 : Set::Scalar T_L=0, T_R=0, p_L=0, p_R=0, gamma_L=0, gamma_R=0;
61 8192000 : switch (side) {
62 2048000 : case 0: // xlo
63 2048000 : T_L = gas.ComputeT(rho_L, lo.M_normal, lo.M_tangent, E_L, 500.0, X, i-1, j ,k);
64 2048000 : p_L = gas.ComputeP(rho_L, T_L, X, i-1, j, k);
65 2048000 : gamma_L = gas.gamma(T_L, X, i-1, j, k);
66 2048000 : T_R = gas.ComputeT(rho_R, hi.M_normal, hi.M_tangent, E_R, 500.0, X, i, j ,k);
67 2048000 : p_R = gas.ComputeP(rho_R, T_R, X, i, j, k);
68 2048000 : gamma_R = gas.gamma(T_R, X, i, j, k);
69 2048000 : break;
70 2048000 : case 1: // xhi
71 2048000 : T_L = gas.ComputeT(rho_L, lo.M_normal, lo.M_tangent, E_L, 500.0, X, i, j ,k);
72 2048000 : p_L = gas.ComputeP(rho_L, T_L, X, i, j, k);
73 2048000 : gamma_L = gas.gamma(T_L, X, i, j, k);
74 2048000 : T_R = gas.ComputeT(rho_R, hi.M_normal, hi.M_tangent, E_R, 500.0, X, i+1, j ,k);
75 2048000 : p_R = gas.ComputeP(rho_R, T_R, X, i+1, j, k);
76 2048000 : gamma_R = gas.gamma(T_R, X, i+1, j, k);
77 2048000 : break;
78 2048000 : case 2: // ylo
79 2048000 : T_L = gas.ComputeT(rho_L, lo.M_normal, lo.M_tangent, E_L, 500.0, X, i, j-1 ,k);
80 2048000 : p_L = gas.ComputeP(rho_L, T_L, X, i, j-1, k);
81 2048000 : gamma_L = gas.gamma(T_L, X, i, j-1, k);
82 2048000 : T_R = gas.ComputeT(rho_R, hi.M_normal, hi.M_tangent, E_R, 500.0, X, i, j ,k);
83 2048000 : p_R = gas.ComputeP(rho_R, T_R, X, i, j, k);
84 2048000 : gamma_R = gas.gamma(T_R, X, i, j, k);
85 2048000 : break;
86 2048000 : case 3: // yhi
87 2048000 : T_L = gas.ComputeT(rho_L, lo.M_normal, lo.M_tangent, E_L, 500.0, X, i, j ,k);
88 2048000 : p_L = gas.ComputeP(rho_L, T_L, X, i, j, k);
89 2048000 : gamma_L = gas.gamma(T_L, X, i, j, k);
90 2048000 : T_R = gas.ComputeT(rho_R, hi.M_normal, hi.M_tangent, E_R, 500.0, X, i, j+1 ,k);
91 2048000 : p_R = gas.ComputeP(rho_R, T_R, X, i, j+1, k);
92 2048000 : gamma_R = gas.gamma(T_R, X, i, j+1, k);
93 2048000 : break;
94 0 : case 4: // zlo
95 0 : T_L = gas.ComputeT(rho_L, lo.M_normal, lo.M_tangent, E_L, 500.0, X, i, j ,k-1);
96 0 : p_L = gas.ComputeP(rho_L, T_L, X, i, j, k-1);
97 0 : gamma_L = gas.gamma(T_L, X, i, j, k-1);
98 0 : T_R = gas.ComputeT(rho_R, hi.M_normal, hi.M_tangent, E_R, 500.0, X, i, j ,k);
99 0 : p_R = gas.ComputeP(rho_R, T_R, X, i, j, k);
100 0 : gamma_R = gas.gamma(T_R, X, i, j, k);
101 0 : break;
102 0 : case 5: // zhi
103 0 : T_L = gas.ComputeT(rho_L, lo.M_normal, lo.M_tangent, E_L, 500.0, X, i, j ,k);
104 0 : p_L = gas.ComputeP(rho_L, T_L, X, i, j, k);
105 0 : gamma_L = gas.gamma(T_L, X, i, j, k);
106 0 : T_R = gas.ComputeT(rho_R, hi.M_normal, hi.M_tangent, E_R, 500.0, X, i, j ,k+1);
107 0 : p_R = gas.ComputeP(rho_R, T_R, X, i, j, k+1);
108 0 : gamma_R = gas.gamma(T_R, X, i, j, k+1);
109 0 : break;
110 0 : default:
111 0 : Util::Abort(INFO, "Unknown side in Riemann solver, ",side);
112 : }
113 :
114 8192000 : Set::Scalar c_L = std::sqrt(gamma_L * p_L / (rho_L + small));
115 8192000 : Set::Scalar c_R = std::sqrt(gamma_R * p_R / (rho_R + small));
116 :
117 : // Max wave speeds
118 8192000 : Set::Scalar S_L = std::min(u_L-c_L, u_R - c_R);
119 8192000 : Set::Scalar S_R = std::max(u_L+c_L, u_R + c_R);
120 :
121 :
122 8192000 : if (lowmach)
123 : {
124 0 : Set::Scalar u_RL = 0.5 * (u_L + u_R);
125 0 : Set::Scalar a_RL = 0.5 * (c_L + c_R);
126 0 : Set::Scalar Mach = std::abs(u_RL) / (a_RL + small);
127 0 : Set::Scalar psi = std::max(Mach, cutoff); // Mach cutoff
128 0 : S_L *= psi;
129 0 : S_R *= psi;
130 : }
131 :
132 : Flux F_L (
133 : lo.M_normal,
134 8192000 : lo.M_normal * u_L + p_L,
135 8192000 : lo.M_tangent * u_L,
136 8192000 : u_L * (lo.E + p_L)
137 8192000 : );
138 : Flux F_R (
139 : hi.M_normal,
140 8192000 : hi.M_normal * u_R + p_R,
141 8192000 : hi.M_tangent * u_R,
142 8192000 : u_R * (hi.E + p_R)
143 8192000 : );
144 :
145 8192000 : if (S_L >= 0.0)
146 : {
147 287544 : return F_L;
148 : }
149 7904456 : else if (S_R <= 0.0)
150 : {
151 0 : return F_R;
152 : }
153 : else
154 : {
155 7904456 : State dState = hi - lo;
156 7904456 : Set::Scalar dS = S_R - S_L;
157 7904456 : return (S_R*F_L - S_L*F_R + S_L*S_R*Flux(dState)) / dS;
158 : }
159 : }
160 : };
161 : }
162 : }
163 : }
164 :
165 : #endif
|
