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 :
14 : namespace Solver {
15 : namespace Local {
16 : namespace Riemann {
17 :
18 : class HLLC : public Riemann
19 : {
20 : public:
21 : static constexpr const char* name = "hllc";
22 1 : HLLC(IO::ParmParse &pp, std::string name) { pp_queryclass(name, *this); }
23 : HLLC(IO::ParmParse &pp) { pp_queryclass(*this); }
24 : HLLC() { IO::ParmParse pp; pp_queryclass(*this); }
25 :
26 1 : static void Parse(HLLC & /*value*/, IO::ParmParse & /*pp*/)
27 : {
28 : // pp.query_default("lowmach", value.lowmach, false);
29 : // pp.query_default("cutoff", value.cutoff, 0.1);
30 1 : }
31 :
32 8192000 : virtual Flux Solve(State lo, State hi, Set::Scalar gamma, Set::Scalar p_ref, Set::Scalar small) override
33 : {
34 : // get densities (and regularize)
35 8192000 : Set::Scalar rho_L = std::max(lo.rho, small);
36 8192000 : Set::Scalar rho_R = std::max(hi.rho, small);
37 : // density average
38 8192000 : Set::Scalar rhobar = 0.5*(rho_L + rho_R);
39 : // density ratio
40 : /*Set::Scalar Rrho = std::sqrt(rho_R / rho_L);*/
41 :
42 : // primitive velocities
43 8192000 : Set::Scalar u_L = lo.M_normal / rho_L;
44 8192000 : Set::Scalar u_R = hi.M_normal / rho_R;
45 8192000 : Set::Scalar v_L = lo.M_tangent / rho_L;
46 8192000 : Set::Scalar v_R = hi.M_tangent / rho_R;
47 : // roe average velocity
48 : /*Set::Scalar ubar = (u_L + u_R*Rrho) / (1.0 + Rrho);*/
49 :
50 : // label the energies
51 8192000 : Set::Scalar E_L = lo.E, E_R = hi.E;
52 :
53 : // primitive pressures
54 8192000 : Set::Scalar p_L = (gamma-1.0) * (E_L - 0.5*rho_L*(u_L*u_L + v_L*v_L)) + p_ref;
55 8192000 : Set::Scalar p_R = (gamma-1.0) * (E_R - 0.5*rho_R*(u_R*u_R + v_R*v_R)) + p_ref;
56 :
57 : // sound speeds
58 8192000 : Set::Scalar a_L = std::sqrt(gamma*p_L / rho_L);
59 8192000 : Set::Scalar a_R = std::sqrt(gamma*p_R / rho_R);
60 : // sound speed average
61 8192000 : Set::Scalar abar = 0.5*(a_L + a_R);
62 :
63 8192000 : Set::Scalar pstar = std::max(0.0, 0.5*(p_L + p_R) - 0.5*(u_R - u_L)*rhobar*abar);
64 :
65 8192000 : Set::Scalar q_L = pstar <= p_L ?
66 : 1.0 :
67 816048 : std::sqrt(1.0 + ((gamma+1)/(2*gamma)) * (pstar/p_L - 1.0));
68 8192000 : Set::Scalar q_R = pstar <= p_R ?
69 : 1.0 :
70 1449064 : std::sqrt(1.0 + ((gamma+1)/(2*gamma)) * (pstar/p_R - 1.0));
71 :
72 : // wave speed estimate
73 8192000 : Set::Scalar S_L = u_L - a_L*q_L;
74 8192000 : Set::Scalar S_R = u_R + a_R*q_R;
75 :
76 : // intermediate speed
77 8192000 : Set::Scalar Sstar =
78 8192000 : (p_R - p_L + rho_L*u_L*(S_L-u_L) - rho_R*u_R*(S_R-u_R)) /
79 8192000 : ( rho_L*(S_L-u_L) - rho_R*(S_R - u_R) );
80 :
81 :
82 8192000 : if (0.0 < S_L) // F_L region
83 : {
84 : Flux F_L ( rho_L*u_L,
85 302032 : rho_L*u_L*u_L + p_L,
86 302032 : rho_L*u_L*v_L,
87 302032 : u_L * (E_L + p_L) );
88 302032 : return F_L;
89 : }
90 7889968 : else if (S_L <= 0.0 && 0.0 < Sstar) // Fstar_L region
91 : {
92 : Flux F_L ( rho_L*u_L,
93 1886656 : rho_L*u_L*u_L + p_L,
94 1886656 : rho_L*u_L*v_L,
95 1886656 : u_L * (E_L + p_L) );
96 : State U_L ( rho_L,
97 : rho_L*u_L,
98 : rho_L*v_L,
99 1886656 : E_L);
100 :
101 : State Ustar_L ( 1.0,
102 : Sstar,
103 : v_L,
104 1886656 : (E_L/rho_L) + (Sstar-u_L) * (Sstar + p_L/(rho_L*(S_L-u_L))));
105 1886656 : Ustar_L *= rho_L*(S_L - u_L)/(S_L-Sstar);
106 :
107 1886656 : return F_L + S_L*(Ustar_L - U_L);
108 : }
109 6003312 : else if (Sstar <= 0.0 && 0.0 <= S_R )
110 : {
111 : Flux F_R ( rho_R*u_R,
112 6003312 : rho_R*u_R*u_R + p_R,
113 6003312 : rho_R*u_R*v_R,
114 6003312 : u_R * (E_R + p_R) );
115 : State U_R ( rho_R,
116 : rho_R*u_R,
117 : rho_R*v_R,
118 6003312 : E_R);
119 :
120 :
121 : State Ustar_R ( 1.0,
122 : Sstar,
123 : v_R,
124 6003312 : (E_R/rho_R) + (Sstar - u_R)*(Sstar + p_R/(rho_R*(S_R-u_R))));
125 6003312 : Ustar_R *= rho_R*(S_R-u_R)/(S_R-Sstar);
126 :
127 6003312 : return F_R + S_R*(Ustar_R - U_R);
128 : }
129 0 : else if (S_R < 0.0)
130 : {
131 : Flux F_R ( rho_R*u_R,
132 0 : rho_R*u_R*u_R + p_R,
133 0 : rho_R*u_R*v_R,
134 0 : u_R * (E_R + p_R) );
135 0 : return F_R;
136 : }
137 0 : Util::Message(INFO,lo);
138 0 : Util::Message(INFO,hi);
139 0 : Util::Message(INFO,S_R);
140 0 : Util::Message(INFO,Sstar);
141 0 : Util::Message(INFO,S_L);
142 0 : Util::Exception(INFO);
143 0 : return Flux(NAN,NAN,NAN,NAN);
144 : }
145 : };
146 :
147 : } // namespace Riemann
148 : } // namespace Local
149 : } // namespace Solver
150 :
151 : #endif
|
