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 : public Riemann
30 : {
31 : public:
32 :
33 :
34 : static constexpr const char* name = "roe";
35 5 : Roe (IO::ParmParse &pp, std::string name)
36 5 : {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 5 : static void Parse(Roe & value, IO::ParmParse & pp)
51 : {
52 : // enable to dump diagnostic data if the roe solver fails
53 10 : pp.query_default("verbose", value.verbose, 1);
54 : // apply entropy fix if tru
55 10 : 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 10 : pp.query_default("lowmach",value.lowmach,false);
60 :
61 :
62 5 : if (value.entropy_fix == 1)
63 4 : Util::Warning(INFO,"The entropy fix is experimental and should be used with caution");
64 4 : else if (value.entropy_fix == 2)
65 4 : Util::Warning(INFO,"The entropy fix is experimental and should be used with caution. Has previously caused errors with FlowDrivenCavity regression test");
66 5 : }
67 :
68 36044800 : 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
69 : {
70 36044800 : Set::Scalar rho_L = lo.rho , rho_R = hi.rho;
71 36044800 : Set::Scalar Mn_L = lo.M_normal , Mn_R = hi.M_normal ;
72 36044800 : Set::Scalar Mt_L = lo.M_tangent , Mt_R = hi.M_tangent ;
73 36044800 : Set::Scalar E_L = lo.E , E_R = hi.E ;
74 :
75 : // Ensure no negative densities
76 36044800 : rho_L = std::max(small,rho_L);
77 36044800 : rho_R = std::max(small,rho_R);
78 :
79 : // STEP 1: Compute fluid primitives
80 :
81 : // This is a modification to Laney allowing for an arbitrary equation of state
82 : //////////////////////////////////////////////////////////////////////////////////////////////////
83 36044800 : Set::Scalar T_L=0, T_R=0, p_L=0, p_R=0, gamma_L=0, gamma_R=0;
84 36044800 : switch (side) {
85 9011200 : case 0: // xlo
86 9011200 : T_L = gas.ComputeT(rho_L, lo.M_normal, lo.M_tangent, E_L, 500.0, X, i-1, j ,k);
87 9011200 : p_L = gas.ComputeP(rho_L, T_L, X, i-1, j, k);
88 9011200 : gamma_L = gas.gamma(T_L, X, i-1, j, k);
89 9011200 : T_R = gas.ComputeT(rho_R, hi.M_normal, hi.M_tangent, E_R, 500.0, X, i, j ,k);
90 9011200 : p_R = gas.ComputeP(rho_R, T_R, X, i, j, k);
91 9011200 : gamma_R = gas.gamma(T_R, X, i, j, k);
92 9011200 : break;
93 9011200 : case 1: // xhi
94 9011200 : T_L = gas.ComputeT(rho_L, lo.M_normal, lo.M_tangent, E_L, 500.0, X, i, j ,k);
95 9011200 : p_L = gas.ComputeP(rho_L, T_L, X, i, j, k);
96 9011200 : gamma_L = gas.gamma(T_L, X, i, j, k);
97 9011200 : T_R = gas.ComputeT(rho_R, hi.M_normal, hi.M_tangent, E_R, 500.0, X, i+1, j ,k);
98 9011200 : p_R = gas.ComputeP(rho_R, T_R, X, i+1, j, k);
99 9011200 : gamma_R = gas.gamma(T_R, X, i+1, j, k);
100 9011200 : break;
101 9011200 : case 2: // ylo
102 9011200 : T_L = gas.ComputeT(rho_L, lo.M_normal, lo.M_tangent, E_L, 500.0, X, i, j-1 ,k);
103 9011200 : p_L = gas.ComputeP(rho_L, T_L, X, i, j-1, k);
104 9011200 : gamma_L = gas.gamma(T_L, X, i, j-1, k);
105 9011200 : T_R = gas.ComputeT(rho_R, hi.M_normal, hi.M_tangent, E_R, 500.0, X, i, j ,k);
106 9011200 : p_R = gas.ComputeP(rho_R, T_R, X, i, j, k);
107 9011200 : gamma_R = gas.gamma(T_R, X, i, j, k);
108 9011200 : break;
109 9011200 : case 3: // yhi
110 9011200 : T_L = gas.ComputeT(rho_L, lo.M_normal, lo.M_tangent, E_L, 500.0, X, i, j ,k);
111 9011200 : p_L = gas.ComputeP(rho_L, T_L, X, i, j, k);
112 9011200 : gamma_L = gas.gamma(T_L, X, i, j, k);
113 9011200 : T_R = gas.ComputeT(rho_R, hi.M_normal, hi.M_tangent, E_R, 500.0, X, i, j+1 ,k);
114 9011200 : p_R = gas.ComputeP(rho_R, T_R, X, i, j+1, k);
115 9011200 : gamma_R = gas.gamma(T_R, X, i, j+1, k);
116 9011200 : break;
117 0 : case 4: // zlo
118 0 : T_L = gas.ComputeT(rho_L, lo.M_normal, lo.M_tangent, E_L, 500.0, X, i, j ,k-1);
119 0 : p_L = gas.ComputeP(rho_L, T_L, X, i, j, k-1);
120 0 : gamma_L = gas.gamma(T_L, X, i, j, k-1);
121 0 : T_R = gas.ComputeT(rho_R, hi.M_normal, hi.M_tangent, E_R, 500.0, X, i, j ,k);
122 0 : p_R = gas.ComputeP(rho_R, T_R, X, i, j, k);
123 0 : gamma_R = gas.gamma(T_R, X, i, j, k);
124 0 : break;
125 0 : case 5: // zhi
126 0 : T_L = gas.ComputeT(rho_L, lo.M_normal, lo.M_tangent, E_L, 500.0, X, i, j ,k);
127 0 : p_L = gas.ComputeP(rho_L, T_L, X, i, j, k);
128 0 : gamma_L = gas.gamma(T_L, X, i, j, k);
129 0 : T_R = gas.ComputeT(rho_R, hi.M_normal, hi.M_tangent, E_R, 500.0, X, i, j ,k+1);
130 0 : p_R = gas.ComputeP(rho_R, T_R, X, i, j, k+1);
131 0 : gamma_R = gas.gamma(T_R, X, i, j, k+1);
132 0 : break;
133 0 : default:
134 0 : Util::Abort(INFO, "Unknown side in Riemann solver, ",side);
135 : }
136 : //////////////////////////////////////////////////////////////////////////////////////////////////
137 :
138 36044800 : Set::Scalar ke_L = 0.5 * (Mn_L * Mn_L /*+ Mt_L * Mt_L*/) / (rho_L+small); // KE per unit volume
139 36044800 : Set::Scalar ue_L = E_L - ke_L; // IE per unit volume
140 36044800 : Set::Scalar h_TL = (ke_L + ue_L + p_L) / (rho_L+small); // specific stagnation enthalpy (per unit mass)
141 :
142 36044800 : Set::Scalar ke_R = 0.5 * (Mn_R * Mn_R /*+ Mt_R * Mt_R*/) / (rho_R+small);
143 36044800 : Set::Scalar ue_R = E_R - ke_R;
144 36044800 : Set::Scalar h_TR = (ke_R + ue_R + p_R) / (rho_R+small);
145 :
146 36044800 : Set::Scalar u_L = Mn_L/(rho_L+small), u_R = Mn_R/(rho_R+small);
147 36044800 : Set::Scalar v_L = Mt_L/(rho_L+small), v_R = Mt_R/(rho_R+small);
148 :
149 : //
150 : // STEP 2: Compute Roe-averaged quantities
151 : //
152 36044800 : Set::Scalar rho_RL = std::sqrt(rho_L * rho_R);
153 36044800 : 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);
154 36044800 : 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);
155 : // This is not thermodynamically consistent to compute gamma_RL - need to compute a Roe-Average composition, temperature, and then gamma for consistency
156 : // But is should be a semi-close appproximation for now.
157 36044800 : Set::Scalar gamma_RL = (std::sqrt(rho_L) * gamma_L + std::sqrt(rho_R) * gamma_R) / (std::sqrt(rho_L) + std::sqrt(rho_R) + small);
158 36044800 : Set::Scalar a_RL_sq = std::max(0.0,(gamma_RL - 1.0) * (h_RL - 0.5 * u_RL * u_RL));
159 :
160 : #ifdef AMREX_DEBUG
161 : if (verbose && ((a_RL_sq<0) || (a_RL_sq!=a_RL_sq)))
162 : {
163 : Util::Message(INFO, "sound speed ", a_RL_sq);
164 :
165 : Util::Message(INFO, "mixed rho ", lo.rho, " ", hi.rho);
166 : Util::Message(INFO, "mixed Mn ", lo.M_normal, " ", hi.M_normal);
167 : Util::Message(INFO, "mixed Mt ", lo.M_tangent, " ", hi.M_tangent);
168 : Util::Message(INFO, "mixed E ", lo.E, " ", hi.E);
169 :
170 : Util::Message(INFO, "fluid rho ", rho_L, " ", rho_R);
171 : Util::Message(INFO, "fluid Mn ", Mn_L, " ", Mn_R);
172 : Util::Message(INFO, "fluid Mt ", Mt_L, " ", Mt_R);
173 : Util::Message(INFO, "fluid E ", E_L, " ", E_R);
174 :
175 : Util::Message(INFO, "fluid rho ", rho_L, " ", rho_R);
176 : Util::Message(INFO, "fluid u ", u_L, " ", u_R);
177 : Util::Message(INFO, "fluid v ", v_L, " ", v_R);
178 : Util::Message(INFO, "fluid p ", p_L, " ", p_R);
179 : }
180 : Util::AssertException(INFO,TEST(a_RL_sq==a_RL_sq)," a_RL_sq is nan/inf; (a_RL_sq=", a_RL_sq,")");
181 : Util::AssertException(INFO,TEST(a_RL_sq>=0), " a_RL_sq is negative; (a_RL_sq=(",a_RL_sq,")");
182 : #endif
183 :
184 36044800 : Set::Scalar a_RL = std::sqrt(a_RL_sq) + small;
185 :
186 : //
187 : // STEP 3: Compute Roe-averaged wave speeds
188 : //
189 36044800 : Set::Scalar lambda1 = u_RL; // 5.53a
190 36044800 : Set::Scalar lambda2 = u_RL + a_RL; // 5.53b
191 36044800 : Set::Scalar lambda3 = u_RL - a_RL; // 5.53c
192 :
193 : //
194 : // STEP 4: Compute wave strengths
195 : //
196 36044800 : Set::Scalar deltarho= rho_R - rho_L;
197 36044800 : Set::Scalar deltap = p_R - p_L;
198 36044800 : Set::Scalar deltau = u_R - u_L;
199 :
200 36044800 : if (lowmach)
201 : {
202 0 : Set::Scalar v_RL = (std::sqrt(rho_L) * v_L + std::sqrt(rho_R) * v_R ) /
203 0 : (std::sqrt(rho_L) + std::sqrt(rho_R) + small);
204 0 : Set::Scalar Ma = (std::abs(u_RL) + std::abs(v_RL)) / (a_RL + small);
205 0 : Ma = std::min(Ma, 1.0);
206 0 : deltau *= Ma;
207 : }
208 :
209 36044800 : Set::Scalar deltav1 = deltarho - deltap / (a_RL_sq + small); // 5.54a
210 36044800 : Set::Scalar deltav2 = deltau + deltap / (rho_RL * a_RL + small); // 5.54b
211 36044800 : Set::Scalar deltav3 = deltau - deltap / (rho_RL * a_RL + small); // 5.54c
212 :
213 : //
214 : // STEP 5: Compute the right eigenvectors
215 : //
216 36044800 : Set::Scalar r11 = 1.0;
217 36044800 : Set::Scalar r12 = u_RL;
218 36044800 : Set::Scalar r13 = 0.5*u_RL*u_RL;
219 36044800 : Set::Scalar r21 = 0.5*rho_RL/a_RL;
220 36044800 : Set::Scalar r22 = 0.5*rho_RL/a_RL * ( u_RL + a_RL );
221 36044800 : Set::Scalar r23 = 0.5*rho_RL/a_RL * ( h_RL + a_RL*u_RL );
222 36044800 : Set::Scalar r31 = -0.5*rho_RL/a_RL;
223 36044800 : Set::Scalar r32 = -0.5*rho_RL/a_RL * ( u_RL - a_RL );
224 36044800 : Set::Scalar r33 = -0.5*rho_RL/a_RL * ( h_RL - a_RL*u_RL );
225 :
226 : //
227 : // STEP 6: Compute solution - not needed since fluxes will be computed in STEP 7
228 : //
229 :
230 : //
231 : // ROE ENTROPY FIX (Source cited in header comments)
232 : //
233 36044800 : if (entropy_fix == 1) { // chimeracfd
234 8192000 : lambda1 = fabs(lambda1);
235 8192000 : lambda2 = fabs(lambda2);
236 8192000 : lambda3 = fabs(lambda3);
237 8192000 : if ( lambda1 < deltau ) lambda1 = 0.5*(lambda1*lambda1 + deltau*deltau)/deltau;
238 8192000 : if ( lambda2 < deltau ) lambda2 = 0.5*(lambda2*lambda2 + deltau*deltau)/deltau;
239 8192000 : if ( lambda3 < deltau ) lambda3 = 0.5*(lambda3*lambda3 + deltau*deltau)/deltau;
240 : }
241 27852800 : else if (entropy_fix == 2) { // Jayanti
242 8192000 : Set::Scalar a_L = std::sqrt(gamma_L * p_L / (rho_L + small)); // sound speed
243 8192000 : Set::Scalar a_R = std::sqrt(gamma_R * p_R / (rho_R + small));
244 8192000 : Set::Scalar lambda1_L = u_L; Set::Scalar lambda1_R = u_R; // eigenvalues
245 8192000 : Set::Scalar lambda2_L = u_L + a_L; Set::Scalar lambda2_R = u_R + a_R;
246 8192000 : Set::Scalar lambda3_L = u_L - a_L; Set::Scalar lambda3_R = u_R - a_R;
247 8192000 : Set::Scalar fix1 = std::max(0.0, std::max(lambda1 - lambda1_L, lambda1_R - lambda1));
248 8192000 : Set::Scalar fix2 = std::max(0.0, std::max(lambda2 - lambda2_L, lambda2_R - lambda2));
249 8192000 : Set::Scalar fix3 = std::max(0.0, std::max(lambda3 - lambda3_L, lambda3_R - lambda3));
250 8192000 : if ( lambda1 < fix1 ) lambda1 = fix1;
251 8192000 : if ( lambda2 < fix2 ) lambda2 = fix2;
252 8192000 : if ( lambda3 < fix3 ) lambda3 = fix3;
253 : }
254 :
255 : //
256 : // STEP 7: Compute fluxes
257 : //
258 36044800 : Flux fl;
259 :
260 36044800 : fl.mass = (0.5*(rho_L*u_L + rho_R*u_R) - 0.5*(
261 36044800 : r11*fabs(lambda1)*deltav1 +
262 36044800 : r21*fabs(lambda2)*deltav2 +
263 36044800 : r31*fabs(lambda3)*deltav3)
264 : );
265 :
266 36044800 : if (fl.mass != fl.mass)
267 : {
268 0 : if (verbose)
269 : {
270 0 : Util::ParallelMessage(INFO,"hi ", hi);
271 0 : Util::ParallelMessage(INFO,"lo ", lo);
272 0 : Util::ParallelMessage(INFO,"rho_R ", rho_R);
273 0 : Util::ParallelMessage(INFO,"rho_L ", rho_L);
274 0 : Util::ParallelMessage(INFO,"rho_RL ", rho_RL);
275 0 : Util::ParallelMessage(INFO,"u_R ", u_R);
276 0 : Util::ParallelMessage(INFO,"u_L ", u_L);
277 0 : Util::ParallelMessage(INFO,"u_RL ", u_RL);
278 0 : Util::ParallelMessage(INFO,"a_RL ", a_RL);
279 0 : Util::ParallelMessage(INFO,"lambda1 ", lambda1);
280 0 : Util::ParallelMessage(INFO,"lambda2 ", lambda2);
281 0 : Util::ParallelMessage(INFO,"lambda3 ", lambda3);
282 0 : Util::ParallelMessage(INFO,"deltav1 ", deltav1);
283 0 : Util::ParallelMessage(INFO,"deltav2 ", deltav2);
284 0 : Util::ParallelMessage(INFO,"deltav3 ", deltav3);
285 : }
286 0 : Util::Exception(INFO);
287 : }
288 :
289 :
290 36044800 : fl.momentum_normal = ( 0.5*(rho_L*u_L*u_L + p_L + rho_R*u_R*u_R + p_R) - 0.5*(
291 36044800 : r12*fabs(lambda1)*deltav1 +
292 36044800 : r22*fabs(lambda2)*deltav2 +
293 36044800 : r32*fabs(lambda3)*deltav3)
294 : );
295 :
296 :
297 36044800 : fl.energy = ( 0.5*(u_L*(ke_L + p_L + ue_L) + u_R*(ke_R + p_R + ue_R)) - 0.5*
298 : (
299 36044800 : r13*fabs(lambda1)*deltav1 +
300 36044800 : r23*fabs(lambda2)*deltav2 +
301 36044800 : r33*fabs(lambda3)*deltav3)
302 : );
303 :
304 : //
305 : // (Update the tangential momentum flux)
306 : //
307 36044800 : fl.momentum_tangent = 0.5 * (rho_L * u_L * v_L + rho_R * u_R * v_R);
308 :
309 72089600 : return fl;
310 : }
311 :
312 :
313 : //static int Test()
314 : //{
315 : // Roe solver;
316 : //
317 :
318 : // int failed = 0;
319 :
320 : // Set::Scalar gamma = 1.4;
321 : // Set::Scalar pref = 10.0;
322 : // Set::Scalar small = 1E-10;
323 :
324 : // // Test 1: Tangential Velocity Difference - No Normal Flux
325 : // try {
326 : // State left (1.0, 1.0, 0.0, 1.0);
327 : // State center(1.0, 1.0, 1.0, 1.0);
328 : // State right (1.0, 1.0, 2.0, 1.0);
329 : // Flux fluxlo = solver.Solve(center, right, gamma, pref, small);
330 : // Flux fluxhi = solver.Solve(left, center, gamma, pref, small);
331 :
332 : // if (fabs(fluxhi.mass - fluxlo.mass) > 1E-10
333 : // || fabs(fluxhi.momentum_normal - fluxlo.momentum_normal) > 1E-10
334 : // || fabs(fluxhi.energy - fluxlo.energy) > 1E-10) {
335 : // Util::Warning(INFO, "left: ",left);
336 : // Util::Warning(INFO, "center: ",center);
337 : // Util::Warning(INFO, "right: ",right);
338 : // Util::Warning(INFO, "Fluxlo: ",fluxlo);
339 : // Util::Warning(INFO, "Fluxhi: ",fluxhi);
340 : // Util::Exception(INFO, "Tangential velocity difference incorrectly affecting normal flux.");
341 : // }
342 : // Util::Test::SubMessage("Test 1: Tangential velocity should induce no normal flux",0);
343 : // } catch (const std::runtime_error& e)
344 : // {
345 : // failed++;
346 : // Util::Test::SubMessage("Test 1: Tangential velocity should induce no normal flux",1);
347 : // }
348 :
349 : // // Test 2: Pure Transverse Velocity Difference
350 : // try {
351 : // State left (1.0, 0.0, 0.0, 1.0);
352 : // State center(1.0, 0.0, 1.0, 1.0);
353 : // State right (1.0, 0.0, 2.0, 1.0);
354 : // Flux fluxlo = solver.Solve(left, center, gamma, pref, small);
355 : // Flux fluxhi = solver.Solve(center, right, gamma, pref, small);
356 : // if (fabs(fluxhi.mass - fluxlo.mass) > 1E-10
357 : // || fabs(fluxhi.momentum_normal - fluxlo.momentum_normal) > 1E-10
358 : // || fabs(fluxhi.energy - fluxlo.energy) > 1E-10) {
359 : // Util::Warning(INFO, "left: ",left);
360 : // Util::Warning(INFO, "center: ",center);
361 : // Util::Warning(INFO, "right: ",right);
362 : // Util::Warning(INFO, "Fluxhi: ",fluxhi);
363 : // Util::Warning(INFO, "Fluxlo: ",fluxlo);
364 : // Util::Exception(INFO, "Pure transverse velocity difference affecting normal flux.");
365 : // }
366 : // Util::Test::SubMessage("Test 2: Pure transverse velocity difference",0);
367 : // } catch (const std::runtime_error& e)
368 : // {
369 : // failed++;
370 : // Util::Test::SubMessage("Test 2: Pure transverse velocity difference",1);
371 : // }
372 :
373 : // // Test 3: Symmetry Test (no flux across identical states)
374 : // try {
375 : // State left(1.0, 0.0, 0.0, 1.0);
376 : // State center(1.0, 0.0, 0.0, 1.0);
377 : // State right(1.0, 0.0, 0.0, 1.0);
378 : // Flux fluxhi = solver.Solve(center, right, gamma, pref, small);
379 : // Flux fluxlo = solver.Solve(left, center, gamma, pref, small);
380 : // if (fabs(fluxhi.mass - fluxlo.mass) > 1E-10 // no change in mass flux
381 : // || fabs(fluxhi.momentum_normal - fluxlo.momentum_normal) > 1E-10 // no change in momentum flux
382 : // || fabs(fluxhi.momentum_tangent) > 1E-10 // zero tangent flux
383 : // || fabs(fluxlo.momentum_tangent) > 1E-10 // zero tangent flux
384 : // || fabs(fluxhi.energy-fluxlo.energy) > 1E-10 // no change in energy flux
385 : // ) {
386 : // Util::Warning(INFO, "left: ",left);
387 : // Util::Warning(INFO, "right: ",right);
388 : // Util::Warning(INFO, "Fluxhi: ",fluxhi);
389 : // Util::Warning(INFO, "Fluxlo: ",fluxlo);
390 : // Util::Exception(INFO, "Symmetric states should result in zero flux.");
391 : // }
392 : // Util::Test::SubMessage("Test 3: Constant states induces no flux difference",0);
393 : // } catch (const std::runtime_error& e)
394 : // {
395 : // failed++;
396 : // Util::Test::SubMessage("Test 3: Constant states induces no flux difference",1);
397 : // }
398 :
399 : // // Test 4: Uniform Flow Test (no flux across uniform flow)
400 : // try {
401 : // State left (1.0, 1.0, 0.5, 1.0);
402 : // State center(1.0, 1.0, 0.5, 1.0);
403 : // State right (1.0, 1.0, 0.5, 1.0);
404 : // Flux fluxhi = solver.Solve(center, right, gamma, pref, small);
405 : // Flux fluxlo = solver.Solve(left, center, gamma, pref, small);
406 : // if (fabs(fluxhi.mass - fluxlo.mass) > 1E-10 ||
407 : // fabs(fluxhi.momentum_normal - fluxlo.momentum_normal) > 1E-10 ||
408 : // fabs(fluxhi.energy - fluxlo.energy) > 1E-10) {
409 : // Util::Warning(INFO, "left: ",left);
410 : // Util::Warning(INFO, "center: ",center);
411 : // Util::Warning(INFO, "right: ",right);
412 : // Util::Warning(INFO, "Fluxlo: ",fluxlo);
413 : // Util::Warning(INFO, "Fluxhi: ",fluxhi);
414 : // Util::Exception(INFO, "Uniform flow should result in no flux.");
415 : // }
416 : // Util::Test::SubMessage("Test 4: Uniform flow should maintain constant flux",0);
417 : // } catch (const std::runtime_error& e)
418 : // {
419 : // failed++;
420 : // Util::Test::SubMessage("Test 4: Uniform flow should maintain constant flux",1);
421 : // }
422 :
423 : // return failed;
424 : //}
425 : };
426 : }
427 : }
428 : }
429 :
430 : #endif
|
