Alamo
HLLC.H
Go to the documentation of this file.
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
16namespace Solver {
17namespace Local {
18namespace Riemann {
19
20class HLLC : public Riemann
21{
22public:
23 static constexpr const char* name = "hllc";
24 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 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 }
33
34 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 Set::Scalar rho_L = std::max(lo.rho, small);
38 Set::Scalar rho_R = std::max(hi.rho, small);
39 // density average
40 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 Set::Scalar u_L = lo.M_normal / rho_L;
46 Set::Scalar u_R = hi.M_normal / rho_R;
47 Set::Scalar v_L = lo.M_tangent / rho_L;
48 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 Set::Scalar E_L = lo.E, E_R = hi.E;
54
55 // primitive pressures
56 Set::Scalar T_L=0, T_R=0, p_L=0, p_R=0, gamma_L=0, gamma_R=0;
57 switch (side) {
58 case 0: // xlo
59 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 p_L = std::max(small, gas.ComputeP(rho_L, T_L, X, i-1, j, k));
61 gamma_L = gas.gamma(T_L, X, i-1, j, k);
62 T_R = std::max(small, gas.ComputeT(rho_R, hi.M_normal, hi.M_tangent, E_R, 500.0, X, i, j ,k));
63 p_R = std::max(small, gas.ComputeP(rho_R, T_R, X, i, j, k));
64 gamma_R = gas.gamma(T_R, X, i, j, k);
65 break;
66 case 1: // xhi
67 T_L = std::max(small, gas.ComputeT(rho_L, lo.M_normal, lo.M_tangent, E_L, 500.0, X, i, j ,k));
68 p_L = std::max(small, gas.ComputeP(rho_L, T_L, X, i, j, k));
69 gamma_L = gas.gamma(T_L, X, i, j, k);
70 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 p_R = std::max(small, gas.ComputeP(rho_R, T_R, X, i+1, j, k));
72 gamma_R = gas.gamma(T_R, X, i+1, j, k);
73 break;
74 case 2: // ylo
75 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 p_L = std::max(small, gas.ComputeP(rho_L, T_L, X, i, j-1, k));
77 gamma_L = gas.gamma(T_L, X, i, j-1, k);
78 T_R = std::max(small, gas.ComputeT(rho_R, hi.M_normal, hi.M_tangent, E_R, 500.0, X, i, j ,k));
79 p_R = std::max(small, gas.ComputeP(rho_R, T_R, X, i, j, k));
80 gamma_R = gas.gamma(T_R, X, i, j, k);
81 break;
82 case 3: // yhi
83 T_L = std::max(small, gas.ComputeT(rho_L, lo.M_normal, lo.M_tangent, E_L, 500.0, X, i, j ,k));
84 p_L = std::max(small, gas.ComputeP(rho_L, T_L, X, i, j, k));
85 gamma_L = gas.gamma(T_L, X, i, j, k);
86 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 p_R = std::max(small, gas.ComputeP(rho_R, T_R, X, i, j+1, k));
88 gamma_R = gas.gamma(T_R, X, i, j+1, k);
89 break;
90 case 4: // zlo
91 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 p_L = std::max(small, gas.ComputeP(rho_L, T_L, X, i, j, k-1));
93 gamma_L = gas.gamma(T_L, X, i, j, k-1);
94 T_R = std::max(small, gas.ComputeT(rho_R, hi.M_normal, hi.M_tangent, E_R, 500.0, X, i, j ,k));
95 p_R = std::max(small, gas.ComputeP(rho_R, T_R, X, i, j, k));
96 gamma_R = gas.gamma(T_R, X, i, j, k);
97 break;
98 case 5: // zhi
99 T_L = std::max(small, gas.ComputeT(rho_L, lo.M_normal, lo.M_tangent, E_L, 500.0, X, i, j ,k));
100 p_L = std::max(small, gas.ComputeP(rho_L, T_L, X, i, j, k));
101 gamma_L = gas.gamma(T_L, X, i, j, k);
102 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 p_R = std::max(small, gas.ComputeP(rho_R, T_R, X, i, j, k+1));
104 gamma_R = gas.gamma(T_R, X, i, j, k+1);
105 break;
106 default:
107 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 Set::Scalar a_L = std::sqrt(gamma_L*p_L / rho_L);
118 Set::Scalar a_R = std::sqrt(gamma_R*p_R / rho_R);
119 // sound speed average
120 Set::Scalar abar = 0.5*(a_L + a_R);
121
122 Set::Scalar pstar = std::max(0.0, 0.5*(p_L + p_R) - 0.5*(u_R - u_L)*rhobar*abar);
123
124 Set::Scalar q_L = pstar <= p_L ?
125 1.0 :
126 std::sqrt(1.0 + ((gamma_L+1.0)/(2.0*gamma_L)) * (pstar/p_L - 1.0));
127 Set::Scalar q_R = pstar <= p_R ?
128 1.0 :
129 std::sqrt(1.0 + ((gamma_R+1.0)/(2.0*gamma_R)) * (pstar/p_R - 1.0));
130
131 // wave speed estimate
132 Set::Scalar S_L = u_L - a_L*q_L;
133 Set::Scalar S_R = u_R + a_R*q_R;
134
135 // intermediate speed
136 Set::Scalar Sstar =
137 (p_R - p_L + rho_L*u_L*(S_L-u_L) - rho_R*u_R*(S_R-u_R)) /
138 ( rho_L*(S_L-u_L) - rho_R*(S_R - u_R) );
139
140
141 if (0.0 < S_L) // F_L region
142 {
143 Flux F_L ( rho_L*u_L,
144 rho_L*u_L*u_L + p_L,
145 rho_L*u_L*v_L,
146 u_L * (E_L + p_L) );
147 return F_L;
148 }
149 else if (S_L <= 0.0 && 0.0 < Sstar) // Fstar_L region
150 {
151 Flux F_L ( rho_L*u_L,
152 rho_L*u_L*u_L + p_L,
153 rho_L*u_L*v_L,
154 u_L * (E_L + p_L) );
155 State U_L ( rho_L,
156 rho_L*u_L,
157 rho_L*v_L,
158 E_L);
159
160 State Ustar_L ( 1.0,
161 Sstar,
162 v_L,
163 (E_L/rho_L) + (Sstar-u_L) * (Sstar + p_L/(rho_L*(S_L-u_L))));
164 Ustar_L *= rho_L*(S_L - u_L)/(S_L-Sstar);
165
166 return F_L + S_L*(Ustar_L - U_L);
167 }
168 else if (Sstar <= 0.0 && 0.0 <= S_R )
169 {
170 Flux F_R ( rho_R*u_R,
171 rho_R*u_R*u_R + p_R,
172 rho_R*u_R*v_R,
173 u_R * (E_R + p_R) );
174 State U_R ( rho_R,
175 rho_R*u_R,
176 rho_R*v_R,
177 E_R);
178
179
180 State Ustar_R ( 1.0,
181 Sstar,
182 v_R,
183 (E_R/rho_R) + (Sstar - u_R)*(Sstar + p_R/(rho_R*(S_R-u_R))));
184 Ustar_R *= rho_R*(S_R-u_R)/(S_R-Sstar);
185
186 return F_R + S_R*(Ustar_R - U_R);
187 }
188 else if (S_R < 0.0)
189 {
190 Flux F_R ( rho_R*u_R,
191 rho_R*u_R*u_R + p_R,
192 rho_R*u_R*v_R,
193 u_R * (E_R + p_R) );
194 return F_R;
195 }
196 Util::Message(INFO,lo);
197 Util::Message(INFO,hi);
198 Util::Message(INFO,S_R);
199 Util::Message(INFO,Sstar);
200 Util::Message(INFO,S_L);
201 Util::Exception(INFO);
202 return Flux(NAN,NAN,NAN,NAN);
203 }
204};
205
206} // namespace Riemann
207} // namespace Local
208} // namespace Solver
209
210#endif
X
#define X(name)
pp_queryclass
#define pp_queryclass(...)
Definition ParmParse.H:109
INFO
#define INFO
Definition Util.H:24
Model::Gas::Gas::ComputeP
double ComputeP(double density, double T, Set::Patch< const Set::Scalar > &X, int i, int j, int k) const
Definition Gas.cpp:66
Model::Gas::Gas::ComputeT
double ComputeT(double density, double momentumx, double momentumy, double E, double Tguess, Set::Patch< const Set::Scalar > &X, int i, int j, int k, double rtol=1e-12) const
Definition Gas.cpp:52
Model::Gas::Gas::gamma
double gamma(double T, Set::Patch< const Set::Scalar > &X, int i, int j, int k) const
Definition Gas.H:129
Solver::Local::Riemann::HLLC::Solve
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
Definition HLLC.H:34
Solver::Local::Riemann::HLLC::HLLC
HLLC(IO::ParmParse &pp, std::string name)
Definition HLLC.H:24
Solver::Local::Riemann::HLLC::name
static constexpr const char * name
Definition HLLC.H:23
Solver::Local::Riemann::HLLC::Parse
static void Parse(HLLC &, IO::ParmParse &)
Definition HLLC.H:28
Solver::Local::Riemann::HLLC::HLLC
HLLC(IO::ParmParse &pp)
Definition HLLC.H:25
Set::Scalar
amrex::Real Scalar
Definition Base.H:18
Set::Vector
Eigen::Matrix< amrex::Real, AMREX_SPACEDIM, 1 > Vector
Definition Base.H:19
Solver
A bunch of solvers.
Definition CG.H:6
Util::Abort
void Abort(const char *msg)
Definition Util.cpp:268
Util::Message
void Message(std::string file, std::string func, int line, Args const &... args)
Definition Util.H:129
Util::Exception
void Exception(std::string file, std::string func, int line, Args const &... args)
Definition Util.H:226
Generated by doxygen 1.9.8