Line data Source code
1 : //
2 : // HLLC flux calculation based on equations 10.67 - 10.73 of
3 : // _Riemann solvers and numerical methods for fluid dynamics_ by E.F.Toro.
4 : //
5 : //
6 :
7 : #ifndef SOLVER_LOCAL_RIEMANN_HLLC_H
8 : #define SOLVER_LOCAL_RIEMANN_HLLC_H
9 :
10 : #include "IO/ParmParse.H"
11 : #include "Solver/Local/Riemann/Riemann.H"
12 : #include "Util/Util.H"
13 : #include "Model/Gas/Gas.H"
14 : #include <memory>
15 :
16 : namespace Solver {
17 : namespace Local {
18 : namespace Riemann {
19 :
20 : class HLLC : public Riemann
21 : {
22 : public:
23 : static constexpr const char* name = "hllc";
24 2 : HLLC(IO::ParmParse &pp, std::string name) { pp_queryclass(name, *this); }
25 : HLLC(IO::ParmParse &pp) { pp_queryclass(*this); }
26 : HLLC() { IO::ParmParse pp; pp_queryclass(*this); }
27 :
28 2 : static void Parse(HLLC & /*value*/, IO::ParmParse & /*pp*/)
29 : {
30 : // pp.query_default("lowmach", value.lowmach, false);
31 : // pp.query_default("cutoff", value.cutoff, 0.1);
32 2 : }
33 :
34 8704000 : 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
35 : {
36 : // get densities (and regularize)
37 8704000 : Set::Scalar rho_L = std::max(lo.rho, small);
38 8704000 : Set::Scalar rho_R = std::max(hi.rho, small);
39 : // density average
40 8704000 : Set::Scalar rhobar = 0.5*(rho_L + rho_R);
41 : // density ratio
42 : /*Set::Scalar Rrho = std::sqrt(rho_R / rho_L);*/
43 :
44 : // primitive velocities
45 8704000 : Set::Scalar u_L = lo.M_normal / rho_L;
46 8704000 : Set::Scalar u_R = hi.M_normal / rho_R;
47 8704000 : Set::Scalar v_L = lo.M_tangent / rho_L;
48 8704000 : Set::Scalar v_R = hi.M_tangent / rho_R;
49 : // roe average velocity
50 : /*Set::Scalar ubar = (u_L + u_R*Rrho) / (1.0 + Rrho);*/
51 :
52 : // label the energies
53 8704000 : Set::Scalar E_L = lo.E, E_R = hi.E;
54 :
55 : // primitive pressures
56 8704000 : Set::Scalar T_L=0, T_R=0, p_L=0, p_R=0, gamma_L=0, gamma_R=0;
57 8704000 : switch (side) {
58 2176000 : case 0: // xlo
59 2176000 : T_L = std::max(small, gas.ComputeT(rho_L, lo.M_normal, lo.M_tangent, E_L, 500.0, X, i-1, j ,k));
60 2176000 : p_L = std::max(small, gas.ComputeP(rho_L, T_L, X, i-1, j, k));
61 2176000 : gamma_L = gas.gamma(T_L, X, i-1, j, k);
62 2176000 : T_R = std::max(small, gas.ComputeT(rho_R, hi.M_normal, hi.M_tangent, E_R, 500.0, X, i, j ,k));
63 2176000 : p_R = std::max(small, gas.ComputeP(rho_R, T_R, X, i, j, k));
64 2176000 : gamma_R = gas.gamma(T_R, X, i, j, k);
65 2176000 : break;
66 2176000 : case 1: // xhi
67 2176000 : T_L = std::max(small, gas.ComputeT(rho_L, lo.M_normal, lo.M_tangent, E_L, 500.0, X, i, j ,k));
68 2176000 : p_L = std::max(small, gas.ComputeP(rho_L, T_L, X, i, j, k));
69 2176000 : gamma_L = gas.gamma(T_L, X, i, j, k);
70 2176000 : T_R = std::max(small, gas.ComputeT(rho_R, hi.M_normal, hi.M_tangent, E_R, 500.0, X, i+1, j ,k));
71 2176000 : p_R = std::max(small, gas.ComputeP(rho_R, T_R, X, i+1, j, k));
72 2176000 : gamma_R = gas.gamma(T_R, X, i+1, j, k);
73 2176000 : break;
74 2176000 : case 2: // ylo
75 2176000 : T_L = std::max(small, gas.ComputeT(rho_L, lo.M_normal, lo.M_tangent, E_L, 500.0, X, i, j-1 ,k));
76 2176000 : p_L = std::max(small, gas.ComputeP(rho_L, T_L, X, i, j-1, k));
77 2176000 : gamma_L = gas.gamma(T_L, X, i, j-1, k);
78 2176000 : T_R = std::max(small, gas.ComputeT(rho_R, hi.M_normal, hi.M_tangent, E_R, 500.0, X, i, j ,k));
79 2176000 : p_R = std::max(small, gas.ComputeP(rho_R, T_R, X, i, j, k));
80 2176000 : gamma_R = gas.gamma(T_R, X, i, j, k);
81 2176000 : break;
82 2176000 : case 3: // yhi
83 2176000 : T_L = std::max(small, gas.ComputeT(rho_L, lo.M_normal, lo.M_tangent, E_L, 500.0, X, i, j ,k));
84 2176000 : p_L = std::max(small, gas.ComputeP(rho_L, T_L, X, i, j, k));
85 2176000 : gamma_L = gas.gamma(T_L, X, i, j, k);
86 2176000 : T_R = std::max(small, gas.ComputeT(rho_R, hi.M_normal, hi.M_tangent, E_R, 500.0, X, i, j+1 ,k));
87 2176000 : p_R = std::max(small, gas.ComputeP(rho_R, T_R, X, i, j+1, k));
88 2176000 : gamma_R = gas.gamma(T_R, X, i, j+1, k);
89 2176000 : break;
90 0 : case 4: // zlo
91 0 : T_L = std::max(small, gas.ComputeT(rho_L, lo.M_normal, lo.M_tangent, E_L, 500.0, X, i, j ,k-1));
92 0 : p_L = std::max(small, gas.ComputeP(rho_L, T_L, X, i, j, k-1));
93 0 : gamma_L = gas.gamma(T_L, X, i, j, k-1);
94 0 : T_R = std::max(small, gas.ComputeT(rho_R, hi.M_normal, hi.M_tangent, E_R, 500.0, X, i, j ,k));
95 0 : p_R = std::max(small, gas.ComputeP(rho_R, T_R, X, i, j, k));
96 0 : gamma_R = gas.gamma(T_R, X, i, j, k);
97 0 : break;
98 0 : case 5: // zhi
99 0 : T_L = std::max(small, gas.ComputeT(rho_L, lo.M_normal, lo.M_tangent, E_L, 500.0, X, i, j ,k));
100 0 : p_L = std::max(small, gas.ComputeP(rho_L, T_L, X, i, j, k));
101 0 : gamma_L = gas.gamma(T_L, X, i, j, k);
102 0 : T_R = std::max(small, gas.ComputeT(rho_R, hi.M_normal, hi.M_tangent, E_R, 500.0, X, i, j ,k+1));
103 0 : p_R = std::max(small, gas.ComputeP(rho_R, T_R, X, i, j, k+1));
104 0 : gamma_R = gas.gamma(T_R, X, i, j, k+1);
105 0 : break;
106 0 : default:
107 0 : Util::Abort(INFO, "Unknown side in Riemann solver, ",side);
108 : }
109 :
110 : //T_L = std::max(T_L, 1.0e-6);
111 : //T_R = std::max(T_R, 1.0e-6);
112 : //Util::Message(INFO,"Side: ", side);
113 : //Util::Message(INFO, "Left: rho/u/v/E/p/T: ", rho_L, " ", u_L, " " , v_L, " ", E_L, " ", p_L, " ", T_L);
114 : //Util::Message(INFO, "Right: rho/u/v/E/p/T: ", rho_R, " ", u_R, " " , v_R, " ", E_R, " ", p_R, " ", T_R);
115 :
116 : // sound speeds
117 8704000 : Set::Scalar a_L = std::sqrt(gamma_L*p_L / rho_L);
118 8704000 : Set::Scalar a_R = std::sqrt(gamma_R*p_R / rho_R);
119 : // sound speed average
120 8704000 : Set::Scalar abar = 0.5*(a_L + a_R);
121 :
122 8704000 : Set::Scalar pstar = std::max(0.0, 0.5*(p_L + p_R) - 0.5*(u_R - u_L)*rhobar*abar);
123 :
124 8704000 : Set::Scalar q_L = pstar <= p_L ?
125 : 1.0 :
126 1055824 : std::sqrt(1.0 + ((gamma_L+1.0)/(2.0*gamma_L)) * (pstar/p_L - 1.0));
127 8704000 : Set::Scalar q_R = pstar <= p_R ?
128 : 1.0 :
129 1515992 : std::sqrt(1.0 + ((gamma_R+1.0)/(2.0*gamma_R)) * (pstar/p_R - 1.0));
130 :
131 : // wave speed estimate
132 8704000 : Set::Scalar S_L = u_L - a_L*q_L;
133 8704000 : Set::Scalar S_R = u_R + a_R*q_R;
134 :
135 : // intermediate speed
136 8704000 : Set::Scalar Sstar =
137 8704000 : (p_R - p_L + rho_L*u_L*(S_L-u_L) - rho_R*u_R*(S_R-u_R)) /
138 8704000 : ( rho_L*(S_L-u_L) - rho_R*(S_R - u_R) );
139 :
140 :
141 8704000 : if (0.0 < S_L) // F_L region
142 : {
143 : Flux F_L ( rho_L*u_L,
144 302008 : rho_L*u_L*u_L + p_L,
145 302008 : rho_L*u_L*v_L,
146 302008 : u_L * (E_L + p_L) );
147 302008 : return F_L;
148 : }
149 8401992 : else if (S_L <= 0.0 && 0.0 < Sstar) // Fstar_L region
150 : {
151 : Flux F_L ( rho_L*u_L,
152 1886680 : rho_L*u_L*u_L + p_L,
153 1886680 : rho_L*u_L*v_L,
154 1886680 : u_L * (E_L + p_L) );
155 : State U_L ( rho_L,
156 : rho_L*u_L,
157 : rho_L*v_L,
158 1886680 : E_L);
159 :
160 : State Ustar_L ( 1.0,
161 : Sstar,
162 : v_L,
163 1886680 : (E_L/rho_L) + (Sstar-u_L) * (Sstar + p_L/(rho_L*(S_L-u_L))));
164 1886680 : Ustar_L *= rho_L*(S_L - u_L)/(S_L-Sstar);
165 :
166 1886680 : return F_L + S_L*(Ustar_L - U_L);
167 : }
168 6515312 : else if (Sstar <= 0.0 && 0.0 <= S_R )
169 : {
170 : Flux F_R ( rho_R*u_R,
171 6515312 : rho_R*u_R*u_R + p_R,
172 6515312 : rho_R*u_R*v_R,
173 6515312 : u_R * (E_R + p_R) );
174 : State U_R ( rho_R,
175 : rho_R*u_R,
176 : rho_R*v_R,
177 6515312 : E_R);
178 :
179 :
180 : State Ustar_R ( 1.0,
181 : Sstar,
182 : v_R,
183 6515312 : (E_R/rho_R) + (Sstar - u_R)*(Sstar + p_R/(rho_R*(S_R-u_R))));
184 6515312 : Ustar_R *= rho_R*(S_R-u_R)/(S_R-Sstar);
185 :
186 6515312 : return F_R + S_R*(Ustar_R - U_R);
187 : }
188 0 : else if (S_R < 0.0)
189 : {
190 : Flux F_R ( rho_R*u_R,
191 0 : rho_R*u_R*u_R + p_R,
192 0 : rho_R*u_R*v_R,
193 0 : u_R * (E_R + p_R) );
194 0 : return F_R;
195 : }
196 0 : Util::Message(INFO,lo);
197 0 : Util::Message(INFO,hi);
198 0 : Util::Message(INFO,S_R);
199 0 : Util::Message(INFO,Sstar);
200 0 : Util::Message(INFO,S_L);
201 0 : Util::Exception(INFO);
202 0 : return Flux(NAN,NAN,NAN,NAN);
203 : }
204 : };
205 :
206 : } // namespace Riemann
207 : } // namespace Local
208 : } // namespace Solver
209 :
210 : #endif
|
