Alamo
Roe.H
Go to the documentation of this file.
1//
2// This implements the Riemann Roe solver.
3//
4// Notation and algorithm follow the presentation in Section 5.3.3
5// of *Computational Gasdynamics* by Culbert B. Laney (page 88)
6//
7// This solver uses an optional entropy fix
8// Option 1: chimeracfd method https://chimeracfd.com/programming/gryphon/fluxroe.html
9// Option 2: Eq. 4.3.67 in *Computational Fluid Dynamics for Engineers and Scientists* by Sreenivas Jayanti
10//
11
12#ifndef SOLVER_LOCAL_RIEMANN_ROE_H
13#define SOLVER_LOCAL_RIEMANN_ROE_H
14
15#include "IO/ParmParse.H"
16#include "Solver/Local/Riemann/Riemann.H"
17
18/// A bunch of solvers
19namespace Solver
20{
21/// Local solvers
22namespace Local
23{
24
25namespace Riemann
26{
27
28/// Roe Riemann Solver based on Gas Dynamics - Culbert B. Laney
29class Roe : public Riemann
30{
31public:
32
33
34 static constexpr const char* name = "roe";
35 Roe (IO::ParmParse &pp, std::string name)
36 {pp_queryclass(name,*this);}
37 Roe (IO::ParmParse &pp)
38 {pp_queryclass(*this);}
39 Roe ()
40 {
41 IO::ParmParse pp;
42 pp_queryclass(*this);
43 }
44
45 int verbose = 0;
46 int entropy_fix = 0;
47 int lowmach = 0;
48 Set::Scalar phi = NAN;
49
50 static void Parse(Roe & value, IO::ParmParse & pp)
51 {
52 // enable to dump diagnostic data if the roe solver fails
53 pp.query_default("verbose", value.verbose, 1);
54 // apply entropy fix if tru
55 pp.query_default("entropy_fix", value.entropy_fix, false);
56
57 // Apply the lowmach fix descripte in Rieper 2010
58 // "A low-Mach number fix for Roe’s approximate Riemann solver"
59 pp.query_default("lowmach",value.lowmach,false);
60
61
62 if (value.entropy_fix == 1)
63 Util::Warning(INFO,"The entropy fix is experimental and should be used with caution");
64 else if (value.entropy_fix == 2)
65 Util::Warning(INFO,"The entropy fix is experimental and should be used with caution. Has previously caused errors with FlowDrivenCavity regression test");
66 }
67
68 virtual Flux Solve(State lo, State hi, Set::Scalar gamma, Set::Scalar p_ref, Set::Scalar small) override
69 {
70 Set::Scalar rho_L = lo.rho , rho_R = hi.rho;
71 Set::Scalar Mn_L = lo.M_normal , Mn_R = hi.M_normal ;
72 Set::Scalar Mt_L = lo.M_tangent , Mt_R = hi.M_tangent ;
73 Set::Scalar E_L = lo.E , E_R = hi.E ;
74
75 // Ensure no negative densities
76 rho_L = std::max(0.0,rho_L);
77 rho_R = std::max(0.0,rho_R);
78
79 // STEP 1: Compute fluid primitives
80 Set::Scalar ke_L = 0.5 * (Mn_L * Mn_L /*+ Mt_L * Mt_L*/) / (rho_L+small); // KE per unit volume
81 Set::Scalar ue_L = E_L - ke_L; // IE per unit volume
82 Set::Scalar p_L = (gamma - 1.0) * ue_L + p_ref; // pressure
83 Set::Scalar h_TL = (ke_L + ue_L + p_L) / (rho_L+small); // specific stagnation enthalpy (per unit mass)
84
85 Set::Scalar ke_R = 0.5 * (Mn_R * Mn_R /*+ Mt_R * Mt_R*/) / (rho_R+small);
86 Set::Scalar ue_R = E_R - ke_R;
87 Set::Scalar p_R = (gamma - 1.0) * ue_R + p_ref;
88 Set::Scalar h_TR = (ke_R + ue_R + p_R) / (rho_R+small);
89
90 Set::Scalar u_L = Mn_L/(rho_L+small), u_R = Mn_R/(rho_R+small);
91 Set::Scalar v_L = Mt_L/(rho_L+small), v_R = Mt_R/(rho_R+small);
92
93 //
94 // STEP 2: Compute Roe-averaged quantities
95 //
96 Set::Scalar rho_RL = std::sqrt(rho_L * rho_R);
97 Set::Scalar u_RL = (std::sqrt(rho_L) * u_L + std::sqrt(rho_R) * u_R ) / (std::sqrt(rho_L) + std::sqrt(rho_R) + small);
98 Set::Scalar h_RL = (std::sqrt(rho_L) * h_TL + std::sqrt(rho_R) * h_TR) / (std::sqrt(rho_L) + std::sqrt(rho_R) + small);
99 Set::Scalar a_RL_sq = std::max(0.0,(gamma - 1.0) * (h_RL - 0.5 * u_RL * u_RL));
100
101
102#ifdef AMREX_DEBUG
103 if (verbose && ((a_RL_sq<0) || (a_RL_sq!=a_RL_sq)))
104 {
105 Util::Message(INFO, "sound speed ", a_RL_sq);
106
107 Util::Message(INFO, "mixed rho ", lo.rho, " ", hi.rho);
108 Util::Message(INFO, "mixed Mn ", lo.M_normal, " ", hi.M_normal);
109 Util::Message(INFO, "mixed Mt ", lo.M_tangent, " ", hi.M_tangent);
110 Util::Message(INFO, "mixed E ", lo.E, " ", hi.E);
111
112 Util::Message(INFO, "fluid rho ", rho_L, " ", rho_R);
113 Util::Message(INFO, "fluid Mn ", Mn_L, " ", Mn_R);
114 Util::Message(INFO, "fluid Mt ", Mt_L, " ", Mt_R);
115 Util::Message(INFO, "fluid E ", E_L, " ", E_R);
116
117 Util::Message(INFO, "fluid rho ", rho_L, " ", rho_R);
118 Util::Message(INFO, "fluid u ", u_L, " ", u_R);
119 Util::Message(INFO, "fluid v ", v_L, " ", v_R);
120 Util::Message(INFO, "fluid p ", p_L, " ", p_R);
121 }
122 Util::AssertException(INFO,TEST(a_RL_sq==a_RL_sq)," a_RL_sq is nan/inf; (a_RL_sq=", a_RL_sq,")");
123 Util::AssertException(INFO,TEST(a_RL_sq>=0), " a_RL_sq is negative; (a_RL_sq=(",a_RL_sq,")");
124#endif
125
126 Set::Scalar a_RL = std::sqrt(a_RL_sq) + small;
127
128 //
129 // STEP 3: Compute Roe-averaged wave speeds
130 //
131 Set::Scalar lambda1 = u_RL; // 5.53a
132 Set::Scalar lambda2 = u_RL + a_RL; // 5.53b
133 Set::Scalar lambda3 = u_RL - a_RL; // 5.53c
134
135 //
136 // STEP 4: Compute wave strengths
137 //
138 Set::Scalar deltarho= rho_R - rho_L;
139 Set::Scalar deltap = p_R - p_L;
140 Set::Scalar deltau = u_R - u_L;
141
142 if (lowmach)
143 {
144 Set::Scalar v_RL = (std::sqrt(rho_L) * v_L + std::sqrt(rho_R) * v_R ) /
145 (std::sqrt(rho_L) + std::sqrt(rho_R) + small);
146 Set::Scalar Ma = (std::abs(u_RL) + std::abs(v_RL)) / (a_RL + small);
147 Ma = std::min(Ma, 1.0);
148 deltau *= Ma;
149 }
150
151 Set::Scalar deltav1 = deltarho - deltap / (a_RL_sq + small); // 5.54a
152 Set::Scalar deltav2 = deltau + deltap / (rho_RL * a_RL + small); // 5.54b
153 Set::Scalar deltav3 = deltau - deltap / (rho_RL * a_RL + small); // 5.54c
154
155 //
156 // STEP 5: Compute the right eigenvectors
157 //
158 Set::Scalar r11 = 1.0;
159 Set::Scalar r12 = u_RL;
160 Set::Scalar r13 = 0.5*u_RL*u_RL;
161 Set::Scalar r21 = 0.5*rho_RL/a_RL;
162 Set::Scalar r22 = 0.5*rho_RL/a_RL * ( u_RL + a_RL );
163 Set::Scalar r23 = 0.5*rho_RL/a_RL * ( h_RL + a_RL*u_RL );
164 Set::Scalar r31 = -0.5*rho_RL/a_RL;
165 Set::Scalar r32 = -0.5*rho_RL/a_RL * ( u_RL - a_RL );
166 Set::Scalar r33 = -0.5*rho_RL/a_RL * ( h_RL - a_RL*u_RL );
167
168 //
169 // STEP 6: Compute solution - not needed since fluxes will be computed in STEP 7
170 //
171
172 //
173 // ROE ENTROPY FIX (Source cited in header comments)
174 //
175 if (entropy_fix == 1) { // chimeracfd
176 lambda1 = fabs(lambda1);
177 lambda2 = fabs(lambda2);
178 lambda3 = fabs(lambda3);
179 if ( lambda1 < deltau ) lambda1 = 0.5*(lambda1*lambda1 + deltau*deltau)/deltau;
180 if ( lambda2 < deltau ) lambda2 = 0.5*(lambda2*lambda2 + deltau*deltau)/deltau;
181 if ( lambda3 < deltau ) lambda3 = 0.5*(lambda3*lambda3 + deltau*deltau)/deltau;
182 }
183 else if (entropy_fix == 2) { // Jayanti
184 Set::Scalar a_L = std::sqrt(gamma * p_L / (rho_L + small)); // sound speed
185 Set::Scalar a_R = std::sqrt(gamma * p_R / (rho_R + small));
186 Set::Scalar lambda1_L = u_L; Set::Scalar lambda1_R = u_R; // eigenvalues
187 Set::Scalar lambda2_L = u_L + a_L; Set::Scalar lambda2_R = u_R + a_R;
188 Set::Scalar lambda3_L = u_L - a_L; Set::Scalar lambda3_R = u_R - a_R;
189 Set::Scalar fix1 = std::max(0.0, std::max(lambda1 - lambda1_L, lambda1_R - lambda1));
190 Set::Scalar fix2 = std::max(0.0, std::max(lambda2 - lambda2_L, lambda2_R - lambda2));
191 Set::Scalar fix3 = std::max(0.0, std::max(lambda3 - lambda3_L, lambda3_R - lambda3));
192 if ( lambda1 < fix1 ) lambda1 = fix1;
193 if ( lambda2 < fix2 ) lambda2 = fix2;
194 if ( lambda3 < fix3 ) lambda3 = fix3;
195 }
196
197 //
198 // STEP 7: Compute fluxes
199 //
200 Flux fl;
201
202 fl.mass = (0.5*(rho_L*u_L + rho_R*u_R) - 0.5*(
203 r11*fabs(lambda1)*deltav1 +
204 r21*fabs(lambda2)*deltav2 +
205 r31*fabs(lambda3)*deltav3)
206 );
207
208 if (fl.mass != fl.mass)
209 {
210 if (verbose)
211 {
212 Util::ParallelMessage(INFO,"hi ", hi);
213 Util::ParallelMessage(INFO,"lo ", lo);
214 Util::ParallelMessage(INFO,"rho_R ", rho_R);
215 Util::ParallelMessage(INFO,"rho_L ", rho_L);
216 Util::ParallelMessage(INFO,"rho_RL ", rho_RL);
217 Util::ParallelMessage(INFO,"u_R ", u_R);
218 Util::ParallelMessage(INFO,"u_L ", u_L);
219 Util::ParallelMessage(INFO,"u_RL ", u_RL);
220 Util::ParallelMessage(INFO,"a_RL ", a_RL);
221 Util::ParallelMessage(INFO,"lambda1 ", lambda1);
222 Util::ParallelMessage(INFO,"lambda2 ", lambda2);
223 Util::ParallelMessage(INFO,"lambda3 ", lambda3);
224 Util::ParallelMessage(INFO,"deltav1 ", deltav1);
225 Util::ParallelMessage(INFO,"deltav2 ", deltav2);
226 Util::ParallelMessage(INFO,"deltav3 ", deltav3);
227 }
228 Util::Exception(INFO);
229 }
230
231
232 fl.momentum_normal = ( 0.5*(rho_L*u_L*u_L + p_L + rho_R*u_R*u_R + p_R) - 0.5*(
233 r12*fabs(lambda1)*deltav1 +
234 r22*fabs(lambda2)*deltav2 +
235 r32*fabs(lambda3)*deltav3)
236 );
237
238
239 fl.energy = ( 0.5*(u_L*(ke_L + p_L + ue_L) + u_R*(ke_R + p_R + ue_R)) - 0.5*
240 (
241 r13*fabs(lambda1)*deltav1 +
242 r23*fabs(lambda2)*deltav2 +
243 r33*fabs(lambda3)*deltav3)
244 );
245
246 //
247 // (Update the tangential momentum flux)
248 //
249 fl.momentum_tangent = 0.5 * (rho_L * u_L * v_L + rho_R * u_R * v_R);
250
251 return fl;
252 }
253
254
255 static int Test()
256 {
257 Roe solver;
258
259
260 int failed = 0;
261
262 Set::Scalar gamma = 1.4;
263 Set::Scalar pref = 10.0;
264 Set::Scalar small = 1E-10;
265
266 // Test 1: Tangential Velocity Difference - No Normal Flux
267 try {
268 State left (1.0, 1.0, 0.0, 1.0);
269 State center(1.0, 1.0, 1.0, 1.0);
270 State right (1.0, 1.0, 2.0, 1.0);
271 Flux fluxlo = solver.Solve(center, right, gamma, pref, small);
272 Flux fluxhi = solver.Solve(left, center, gamma, pref, small);
273
274 if (fabs(fluxhi.mass - fluxlo.mass) > 1E-10
275 || fabs(fluxhi.momentum_normal - fluxlo.momentum_normal) > 1E-10
276 || fabs(fluxhi.energy - fluxlo.energy) > 1E-10) {
277 Util::Warning(INFO, "left: ",left);
278 Util::Warning(INFO, "center: ",center);
279 Util::Warning(INFO, "right: ",right);
280 Util::Warning(INFO, "Fluxlo: ",fluxlo);
281 Util::Warning(INFO, "Fluxhi: ",fluxhi);
282 Util::Exception(INFO, "Tangential velocity difference incorrectly affecting normal flux.");
283 }
284 Util::Test::SubMessage("Test 1: Tangential velocity should induce no normal flux",0);
285 } catch (const std::runtime_error& e)
286 {
287 failed++;
288 Util::Test::SubMessage("Test 1: Tangential velocity should induce no normal flux",1);
289 }
290
291 // Test 2: Pure Transverse Velocity Difference
292 try {
293 State left (1.0, 0.0, 0.0, 1.0);
294 State center(1.0, 0.0, 1.0, 1.0);
295 State right (1.0, 0.0, 2.0, 1.0);
296 Flux fluxlo = solver.Solve(left, center, gamma, pref, small);
297 Flux fluxhi = solver.Solve(center, right, gamma, pref, small);
298 if (fabs(fluxhi.mass - fluxlo.mass) > 1E-10
299 || fabs(fluxhi.momentum_normal - fluxlo.momentum_normal) > 1E-10
300 || fabs(fluxhi.energy - fluxlo.energy) > 1E-10) {
301 Util::Warning(INFO, "left: ",left);
302 Util::Warning(INFO, "center: ",center);
303 Util::Warning(INFO, "right: ",right);
304 Util::Warning(INFO, "Fluxhi: ",fluxhi);
305 Util::Warning(INFO, "Fluxlo: ",fluxlo);
306 Util::Exception(INFO, "Pure transverse velocity difference affecting normal flux.");
307 }
308 Util::Test::SubMessage("Test 2: Pure transverse velocity difference",0);
309 } catch (const std::runtime_error& e)
310 {
311 failed++;
312 Util::Test::SubMessage("Test 2: Pure transverse velocity difference",1);
313 }
314
315 // Test 3: Symmetry Test (no flux across identical states)
316 try {
317 State left(1.0, 0.0, 0.0, 1.0);
318 State center(1.0, 0.0, 0.0, 1.0);
319 State right(1.0, 0.0, 0.0, 1.0);
320 Flux fluxhi = solver.Solve(center, right, gamma, pref, small);
321 Flux fluxlo = solver.Solve(left, center, gamma, pref, small);
322 if (fabs(fluxhi.mass - fluxlo.mass) > 1E-10 // no change in mass flux
323 || fabs(fluxhi.momentum_normal - fluxlo.momentum_normal) > 1E-10 // no change in momentum flux
324 || fabs(fluxhi.momentum_tangent) > 1E-10 // zero tangent flux
325 || fabs(fluxlo.momentum_tangent) > 1E-10 // zero tangent flux
326 || fabs(fluxhi.energy-fluxlo.energy) > 1E-10 // no change in energy flux
327 ) {
328 Util::Warning(INFO, "left: ",left);
329 Util::Warning(INFO, "right: ",right);
330 Util::Warning(INFO, "Fluxhi: ",fluxhi);
331 Util::Warning(INFO, "Fluxlo: ",fluxlo);
332 Util::Exception(INFO, "Symmetric states should result in zero flux.");
333 }
334 Util::Test::SubMessage("Test 3: Constant states induces no flux difference",0);
335 } catch (const std::runtime_error& e)
336 {
337 failed++;
338 Util::Test::SubMessage("Test 3: Constant states induces no flux difference",1);
339 }
340
341 // Test 4: Uniform Flow Test (no flux across uniform flow)
342 try {
343 State left (1.0, 1.0, 0.5, 1.0);
344 State center(1.0, 1.0, 0.5, 1.0);
345 State right (1.0, 1.0, 0.5, 1.0);
346 Flux fluxhi = solver.Solve(center, right, gamma, pref, small);
347 Flux fluxlo = solver.Solve(left, center, gamma, pref, small);
348 if (fabs(fluxhi.mass - fluxlo.mass) > 1E-10 ||
349 fabs(fluxhi.momentum_normal - fluxlo.momentum_normal) > 1E-10 ||
350 fabs(fluxhi.energy - fluxlo.energy) > 1E-10) {
351 Util::Warning(INFO, "left: ",left);
352 Util::Warning(INFO, "center: ",center);
353 Util::Warning(INFO, "right: ",right);
354 Util::Warning(INFO, "Fluxlo: ",fluxlo);
355 Util::Warning(INFO, "Fluxhi: ",fluxhi);
356 Util::Exception(INFO, "Uniform flow should result in no flux.");
357 }
358 Util::Test::SubMessage("Test 4: Uniform flow should maintain constant flux",0);
359 } catch (const std::runtime_error& e)
360 {
361 failed++;
362 Util::Test::SubMessage("Test 4: Uniform flow should maintain constant flux",1);
363 }
364
365 return failed;
366 }
367};
368}
369}
370}
371
372#endif
pp_queryclass
#define pp_queryclass(...)
Definition ParmParse.H:107
TEST
#define TEST(x)
Definition Util.H:21
INFO
#define INFO
Definition Util.H:20
IO::ParmParse::query_default
int query_default(std::string name, T &value, T defaultvalue, std::string="", std::string="", int=-1)
Definition ParmParse.H:203
Solver::Local::Riemann::Roe
Roe Riemann Solver based on Gas Dynamics - Culbert B. Laney.
Definition Roe.H:30
Solver::Local::Riemann::Roe::Roe
Roe(IO::ParmParse &pp)
Definition Roe.H:37
Solver::Local::Riemann::Roe::Parse
static void Parse(Roe &value, IO::ParmParse &pp)
Definition Roe.H:50
Solver::Local::Riemann::Roe::Solve
virtual Flux Solve(State lo, State hi, Set::Scalar gamma, Set::Scalar p_ref, Set::Scalar small) override
Definition Roe.H:68
Solver::Local::Riemann::Roe::Roe
Roe(IO::ParmParse &pp, std::string name)
Definition Roe.H:35
Solver::Local::Riemann::Roe::name
static constexpr const char * name
Definition Roe.H:34
Solver::Local::Riemann::Roe::Test
static int Test()
Definition Roe.H:255
Set::Scalar
amrex::Real Scalar
Definition Base.H:19
Solver
A bunch of solvers.
Definition CG.H:6
Util::Test::SubMessage
int SubMessage(std::string testname, int failed)
Definition Util.cpp:369
Util::ParallelMessage
void ParallelMessage(std::string file, std::string func, int line, Args const &... args)
Definition Util.H:164
Util::AssertException
AMREX_FORCE_INLINE void AssertException(std::string file, std::string func, int line, std::string smt, bool pass, Args const &... args)
Definition Util.H:233
Util::Warning
void Warning(std::string file, std::string func, int line, Args const &... args)
Definition Util.H:181
Util::Message
void Message(std::string file, std::string func, int line, Args const &... args)
Definition Util.H:141
Util::Exception
void Exception(std::string file, std::string func, int line, Args const &... args)
Definition Util.H:205
Solver::Local::Riemann::Flux::momentum_normal
Set::Scalar momentum_normal
Definition Riemann.H:92
Solver::Local::Riemann::Flux::momentum_tangent
Set::Scalar momentum_tangent
Definition Riemann.H:93
Generated by doxygen 1.9.8