Skip to content

Commit ee2b2c7

Browse files
tomidakntomidakn
authored
Merge pull request #620 from PrincetonUniversity/HLLD_Fix
Removed a degeneracy check from the HLLD and LHLLD Riemann Solvers.
2 parents 1591aab + 0e3ca3d commit ee2b2c7

File tree

2 files changed

+91
-101
lines changed

2 files changed

+91
-101
lines changed

‎src/hydro/rsolvers/mhd/hlld.cpp

+45-50
Original file line numberDiff line numberDiff line change
@@ -237,46 +237,41 @@ void Hydro::RiemannSolver(const int k, const int j, const int il, const int iu,
237237
urst.e = (sdr*ur.e - ptr*wri[IVX] + ptst*spd[2] +
238238
bxi*(wri[IVX]*bxi + (wri[IVY]*ur.by + wri[IVZ]*ur.bz) - vbstr))*sdmr_inv;
239239
// ul** and ur** - if Bx is near zero, same as *-states
240-
if (0.5*bxsq < (SMALL_NUMBER)*ptst) {
241-
uldst = ulst;
242-
urdst = urst;
243-
} else {
244-
Real invsumd = 1.0/(sqrtdl + sqrtdr);
245-
Real bxsig = (bxi > 0.0 ? 1.0 : -1.0);
246-
247-
uldst.d = ulst.d;
248-
urdst.d = urst.d;
249-
250-
uldst.mx = ulst.mx;
251-
urdst.mx = urst.mx;
252-
253-
// eqn (59) of M&K
254-
Real tmp = invsumd*(sqrtdl*(ulst.my*ulst_d_inv) + sqrtdr*(urst.my*urst_d_inv) +
255-
bxsig*(urst.by - ulst.by));
256-
uldst.my = uldst.d * tmp;
257-
urdst.my = urdst.d * tmp;
258-
259-
// eqn (60) of M&K
260-
tmp = invsumd*(sqrtdl*(ulst.mz*ulst_d_inv) + sqrtdr*(urst.mz*urst_d_inv) +
261-
bxsig*(urst.bz - ulst.bz));
262-
uldst.mz = uldst.d * tmp;
263-
urdst.mz = urdst.d * tmp;
264-
265-
// eqn (61) of M&K
266-
tmp = invsumd*(sqrtdl*urst.by + sqrtdr*ulst.by +
267-
bxsig*sqrtdl*sqrtdr*((urst.my*urst_d_inv) - (ulst.my*ulst_d_inv)));
268-
uldst.by = urdst.by = tmp;
269-
270-
// eqn (62) of M&K
271-
tmp = invsumd*(sqrtdl*urst.bz + sqrtdr*ulst.bz +
272-
bxsig*sqrtdl*sqrtdr*((urst.mz*urst_d_inv) - (ulst.mz*ulst_d_inv)));
273-
uldst.bz = urdst.bz = tmp;
274-
275-
// eqn (63) of M&K
276-
tmp = spd[2]*bxi + (uldst.my*uldst.by + uldst.mz*uldst.bz)/uldst.d;
277-
uldst.e = ulst.e - sqrtdl*bxsig*(vbstl - tmp);
278-
urdst.e = urst.e + sqrtdr*bxsig*(vbstr - tmp);
279-
}
240+
Real invsumd = 1.0/(sqrtdl + sqrtdr);
241+
Real bxsig = (bxi > 0.0 ? 1.0 : -1.0);
242+
243+
uldst.d = ulst.d;
244+
urdst.d = urst.d;
245+
246+
uldst.mx = ulst.mx;
247+
urdst.mx = urst.mx;
248+
249+
// eqn (59) of M&K
250+
Real tmp = invsumd*(sqrtdl*(ulst.my*ulst_d_inv) + sqrtdr*(urst.my*urst_d_inv) +
251+
bxsig*(urst.by - ulst.by));
252+
uldst.my = uldst.d * tmp;
253+
urdst.my = urdst.d * tmp;
254+
255+
// eqn (60) of M&K
256+
tmp = invsumd*(sqrtdl*(ulst.mz*ulst_d_inv) + sqrtdr*(urst.mz*urst_d_inv) +
257+
bxsig*(urst.bz - ulst.bz));
258+
uldst.mz = uldst.d * tmp;
259+
urdst.mz = urdst.d * tmp;
260+
261+
// eqn (61) of M&K
262+
tmp = invsumd*(sqrtdl*urst.by + sqrtdr*ulst.by +
263+
bxsig*sqrtdl*sqrtdr*((urst.my*urst_d_inv) - (ulst.my*ulst_d_inv)));
264+
uldst.by = urdst.by = tmp;
265+
266+
// eqn (62) of M&K
267+
tmp = invsumd*(sqrtdl*urst.bz + sqrtdr*ulst.bz +
268+
bxsig*sqrtdl*sqrtdr*((urst.mz*urst_d_inv) - (ulst.mz*ulst_d_inv)));
269+
uldst.bz = urdst.bz = tmp;
270+
271+
// eqn (63) of M&K
272+
tmp = spd[2]*bxi + (uldst.my*uldst.by + uldst.mz*uldst.bz)/uldst.d;
273+
uldst.e = ulst.e - sqrtdl*bxsig*(vbstl - tmp);
274+
urdst.e = urst.e + sqrtdr*bxsig*(vbstr - tmp);
280275

281276
//--- Step 6. Compute flux
282277
uldst.d = spd[1] * (uldst.d - ulst.d);
@@ -338,6 +333,15 @@ void Hydro::RiemannSolver(const int k, const int j, const int il, const int iu,
338333
flxi[IEN] = fl.e + ulst.e;
339334
flxi[IBY] = fl.by + ulst.by;
340335
flxi[IBZ] = fl.bz + ulst.bz;
336+
} else if (spd[3] <= 0.0) {
337+
// return Fr*
338+
flxi[IDN] = fr.d + urst.d;
339+
flxi[IVX] = fr.mx + urst.mx;
340+
flxi[IVY] = fr.my + urst.my;
341+
flxi[IVZ] = fr.mz + urst.mz;
342+
flxi[IEN] = fr.e + urst.e;
343+
flxi[IBY] = fr.by + urst.by;
344+
flxi[IBZ] = fr.bz + urst.bz;
341345
} else if (spd[2] >= 0.0) {
342346
// return Fl**
343347
flxi[IDN] = fl.d + ulst.d + uldst.d;
@@ -347,7 +351,7 @@ void Hydro::RiemannSolver(const int k, const int j, const int il, const int iu,
347351
flxi[IEN] = fl.e + ulst.e + uldst.e;
348352
flxi[IBY] = fl.by + ulst.by + uldst.by;
349353
flxi[IBZ] = fl.bz + ulst.bz + uldst.bz;
350-
} else if (spd[3] > 0.0) {
354+
} else {
351355
// return Fr**
352356
flxi[IDN] = fr.d + urst.d + urdst.d;
353357
flxi[IVX] = fr.mx + urst.mx + urdst.mx;
@@ -356,15 +360,6 @@ void Hydro::RiemannSolver(const int k, const int j, const int il, const int iu,
356360
flxi[IEN] = fr.e + urst.e + urdst.e;
357361
flxi[IBY] = fr.by + urst.by + urdst.by;
358362
flxi[IBZ] = fr.bz + urst.bz + urdst.bz;
359-
} else {
360-
// return Fr*
361-
flxi[IDN] = fr.d + urst.d;
362-
flxi[IVX] = fr.mx + urst.mx;
363-
flxi[IVY] = fr.my + urst.my;
364-
flxi[IVZ] = fr.mz + urst.mz;
365-
flxi[IEN] = fr.e + urst.e;
366-
flxi[IBY] = fr.by + urst.by;
367-
flxi[IBZ] = fr.bz + urst.bz;
368363
}
369364

370365
flx(IDN,k,j,i) = flxi[IDN];

‎src/hydro/rsolvers/mhd/lhlld.cpp

+46-51
Original file line numberDiff line numberDiff line change
@@ -242,47 +242,42 @@ void Hydro::RiemannSolver(const int k, const int j, const int il, const int iu,
242242
// (KGF): group transverse by, bz terms for floating-point associativity symmetry
243243
urst.e = (sdr*ur.e - ptr*wri[IVX] + ptst*spd[2] +
244244
bxi*(wri[IVX]*bxi + (wri[IVY]*ur.by + wri[IVZ]*ur.bz) - vbstr))*sdmr_inv;
245-
// ul** and ur** - if Bx is near zero, same as *-states
246-
if (0.5*bxsq < (SMALL_NUMBER)*ptst) {
247-
uldst = ulst;
248-
urdst = urst;
249-
} else {
250-
Real invsumd = 1.0/(sqrtdl + sqrtdr);
251-
Real bxsig = (bxi > 0.0 ? 1.0 : -1.0);
252-
253-
uldst.d = ulst.d;
254-
urdst.d = urst.d;
255-
256-
uldst.mx = ulst.mx;
257-
urdst.mx = urst.mx;
258-
259-
// eqn (59) of M&K
260-
Real tmp = invsumd*(sqrtdl*(ulst.my*ulst_d_inv) + sqrtdr*(urst.my*urst_d_inv) +
261-
bxsig*(urst.by - ulst.by));
262-
uldst.my = uldst.d * tmp;
263-
urdst.my = urdst.d * tmp;
264-
265-
// eqn (60) of M&K
266-
tmp = invsumd*(sqrtdl*(ulst.mz*ulst_d_inv) + sqrtdr*(urst.mz*urst_d_inv) +
267-
bxsig*(urst.bz - ulst.bz));
268-
uldst.mz = uldst.d * tmp;
269-
urdst.mz = urdst.d * tmp;
270-
271-
// eqn (61) of M&K
272-
tmp = invsumd*(sqrtdl*urst.by + sqrtdr*ulst.by +
273-
bxsig*sqrtdl*sqrtdr*((urst.my*urst_d_inv) - (ulst.my*ulst_d_inv)));
274-
uldst.by = urdst.by = tmp;
275-
276-
// eqn (62) of M&K
277-
tmp = invsumd*(sqrtdl*urst.bz + sqrtdr*ulst.bz +
278-
bxsig*sqrtdl*sqrtdr*((urst.mz*urst_d_inv) - (ulst.mz*ulst_d_inv)));
279-
uldst.bz = urdst.bz = tmp;
280-
281-
// eqn (63) of M&K
282-
tmp = spd[2]*bxi + (uldst.my*uldst.by + uldst.mz*uldst.bz)/uldst.d;
283-
uldst.e = ulst.e - sqrtdl*bxsig*(vbstl - tmp);
284-
urdst.e = urst.e + sqrtdr*bxsig*(vbstr - tmp);
285-
}
245+
246+
Real invsumd = 1.0/(sqrtdl + sqrtdr);
247+
Real bxsig = (bxi > 0.0 ? 1.0 : -1.0);
248+
249+
uldst.d = ulst.d;
250+
urdst.d = urst.d;
251+
252+
uldst.mx = ulst.mx;
253+
urdst.mx = urst.mx;
254+
255+
// eqn (59) of M&K
256+
Real tmp = invsumd*(sqrtdl*(ulst.my*ulst_d_inv) + sqrtdr*(urst.my*urst_d_inv) +
257+
bxsig*(urst.by - ulst.by));
258+
uldst.my = uldst.d * tmp;
259+
urdst.my = urdst.d * tmp;
260+
261+
// eqn (60) of M&K
262+
tmp = invsumd*(sqrtdl*(ulst.mz*ulst_d_inv) + sqrtdr*(urst.mz*urst_d_inv) +
263+
bxsig*(urst.bz - ulst.bz));
264+
uldst.mz = uldst.d * tmp;
265+
urdst.mz = urdst.d * tmp;
266+
267+
// eqn (61) of M&K
268+
tmp = invsumd*(sqrtdl*urst.by + sqrtdr*ulst.by +
269+
bxsig*sqrtdl*sqrtdr*((urst.my*urst_d_inv) - (ulst.my*ulst_d_inv)));
270+
uldst.by = urdst.by = tmp;
271+
272+
// eqn (62) of M&K
273+
tmp = invsumd*(sqrtdl*urst.bz + sqrtdr*ulst.bz +
274+
bxsig*sqrtdl*sqrtdr*((urst.mz*urst_d_inv) - (ulst.mz*ulst_d_inv)));
275+
uldst.bz = urdst.bz = tmp;
276+
277+
// eqn (63) of M&K
278+
tmp = spd[2]*bxi + (uldst.my*uldst.by + uldst.mz*uldst.bz)/uldst.d;
279+
uldst.e = ulst.e - sqrtdl*bxsig*(vbstl - tmp);
280+
urdst.e = urst.e + sqrtdr*bxsig*(vbstr - tmp);
286281

287282
//--- Step 6. Compute flux
288283
uldst.d = spd[1] * (uldst.d - ulst.d);
@@ -344,6 +339,15 @@ void Hydro::RiemannSolver(const int k, const int j, const int il, const int iu,
344339
flxi[IEN] = fl.e + ulst.e;
345340
flxi[IBY] = fl.by + ulst.by;
346341
flxi[IBZ] = fl.bz + ulst.bz;
342+
} else if (spd[3] <= 0.0) {
343+
// return Fr*
344+
flxi[IDN] = fr.d + urst.d;
345+
flxi[IVX] = fr.mx + urst.mx;
346+
flxi[IVY] = fr.my + urst.my;
347+
flxi[IVZ] = fr.mz + urst.mz;
348+
flxi[IEN] = fr.e + urst.e;
349+
flxi[IBY] = fr.by + urst.by;
350+
flxi[IBZ] = fr.bz + urst.bz;
347351
} else if (spd[2] >= 0.0) {
348352
// return Fl**
349353
flxi[IDN] = fl.d + ulst.d + uldst.d;
@@ -353,7 +357,7 @@ void Hydro::RiemannSolver(const int k, const int j, const int il, const int iu,
353357
flxi[IEN] = fl.e + ulst.e + uldst.e;
354358
flxi[IBY] = fl.by + ulst.by + uldst.by;
355359
flxi[IBZ] = fl.bz + ulst.bz + uldst.bz;
356-
} else if (spd[3] > 0.0) {
360+
} else {
357361
// return Fr**
358362
flxi[IDN] = fr.d + urst.d + urdst.d;
359363
flxi[IVX] = fr.mx + urst.mx + urdst.mx;
@@ -362,15 +366,6 @@ void Hydro::RiemannSolver(const int k, const int j, const int il, const int iu,
362366
flxi[IEN] = fr.e + urst.e + urdst.e;
363367
flxi[IBY] = fr.by + urst.by + urdst.by;
364368
flxi[IBZ] = fr.bz + urst.bz + urdst.bz;
365-
} else {
366-
// return Fr*
367-
flxi[IDN] = fr.d + urst.d;
368-
flxi[IVX] = fr.mx + urst.mx;
369-
flxi[IVY] = fr.my + urst.my;
370-
flxi[IVZ] = fr.mz + urst.mz;
371-
flxi[IEN] = fr.e + urst.e;
372-
flxi[IBY] = fr.by + urst.by;
373-
flxi[IBZ] = fr.bz + urst.bz;
374369
}
375370

376371
flx(IDN,k,j,i) = flxi[IDN];

0 commit comments

Comments
 (0)