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