Skip to content

Commit

Permalink
Merge pull request #620 from PrincetonUniversity/HLLD_Fix
Browse files Browse the repository at this point in the history 
Removed a degeneracy check from the HLLD and LHLLD Riemann Solvers.
  • Loading branch information
@tomidakn
tomidakn authored Oct 12, 2024
2 parents 1591aab + 0e3ca3d commit ee2b2c7
Show file tree
Hide file tree
Showing 2 changed files with 91 additions and 101 deletions.
95 changes: 45 additions & 50 deletions src/hydro/rsolvers/mhd/hlld.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -237,46 +237,41 @@ void Hydro::RiemannSolver(const int k, const int j, const int il, const int iu,
urst.e = (sdr*ur.e - ptr*wri[IVX] + ptst*spd[2] +
bxi*(wri[IVX]*bxi + (wri[IVY]*ur.by + wri[IVZ]*ur.bz) - vbstr))*sdmr_inv;
// ul** and ur** - if Bx is near zero, same as *-states
if (0.5*bxsq < (SMALL_NUMBER)*ptst) {
uldst = ulst;
urdst = urst;
} else {
Real invsumd = 1.0/(sqrtdl + sqrtdr);
Real bxsig = (bxi > 0.0 ? 1.0 : -1.0);

uldst.d = ulst.d;
urdst.d = urst.d;

uldst.mx = ulst.mx;
urdst.mx = urst.mx;

// eqn (59) of M&K
Real tmp = invsumd*(sqrtdl*(ulst.my*ulst_d_inv) + sqrtdr*(urst.my*urst_d_inv) +
bxsig*(urst.by - ulst.by));
uldst.my = uldst.d * tmp;
urdst.my = urdst.d * tmp;

// eqn (60) of M&K
tmp = invsumd*(sqrtdl*(ulst.mz*ulst_d_inv) + sqrtdr*(urst.mz*urst_d_inv) +
bxsig*(urst.bz - ulst.bz));
uldst.mz = uldst.d * tmp;
urdst.mz = urdst.d * tmp;

// eqn (61) of M&K
tmp = invsumd*(sqrtdl*urst.by + sqrtdr*ulst.by +
bxsig*sqrtdl*sqrtdr*((urst.my*urst_d_inv) - (ulst.my*ulst_d_inv)));
uldst.by = urdst.by = tmp;

// eqn (62) of M&K
tmp = invsumd*(sqrtdl*urst.bz + sqrtdr*ulst.bz +
bxsig*sqrtdl*sqrtdr*((urst.mz*urst_d_inv) - (ulst.mz*ulst_d_inv)));
uldst.bz = urdst.bz = tmp;

// eqn (63) of M&K
tmp = spd[2]*bxi + (uldst.my*uldst.by + uldst.mz*uldst.bz)/uldst.d;
uldst.e = ulst.e - sqrtdl*bxsig*(vbstl - tmp);
urdst.e = urst.e + sqrtdr*bxsig*(vbstr - tmp);
}
Real invsumd = 1.0/(sqrtdl + sqrtdr);
Real bxsig = (bxi > 0.0 ? 1.0 : -1.0);

uldst.d = ulst.d;
urdst.d = urst.d;

uldst.mx = ulst.mx;
urdst.mx = urst.mx;

// eqn (59) of M&K
Real tmp = invsumd*(sqrtdl*(ulst.my*ulst_d_inv) + sqrtdr*(urst.my*urst_d_inv) +
bxsig*(urst.by - ulst.by));
uldst.my = uldst.d * tmp;
urdst.my = urdst.d * tmp;

// eqn (60) of M&K
tmp = invsumd*(sqrtdl*(ulst.mz*ulst_d_inv) + sqrtdr*(urst.mz*urst_d_inv) +
bxsig*(urst.bz - ulst.bz));
uldst.mz = uldst.d * tmp;
urdst.mz = urdst.d * tmp;

// eqn (61) of M&K
tmp = invsumd*(sqrtdl*urst.by + sqrtdr*ulst.by +
bxsig*sqrtdl*sqrtdr*((urst.my*urst_d_inv) - (ulst.my*ulst_d_inv)));
uldst.by = urdst.by = tmp;

// eqn (62) of M&K
tmp = invsumd*(sqrtdl*urst.bz + sqrtdr*ulst.bz +
bxsig*sqrtdl*sqrtdr*((urst.mz*urst_d_inv) - (ulst.mz*ulst_d_inv)));
uldst.bz = urdst.bz = tmp;

// eqn (63) of M&K
tmp = spd[2]*bxi + (uldst.my*uldst.by + uldst.mz*uldst.bz)/uldst.d;
uldst.e = ulst.e - sqrtdl*bxsig*(vbstl - tmp);
urdst.e = urst.e + sqrtdr*bxsig*(vbstr - tmp);

//--- Step 6. Compute flux
uldst.d = spd[1] * (uldst.d - ulst.d);
Expand Down Expand Up @@ -338,6 +333,15 @@ void Hydro::RiemannSolver(const int k, const int j, const int il, const int iu,
flxi[IEN] = fl.e + ulst.e;
flxi[IBY] = fl.by + ulst.by;
flxi[IBZ] = fl.bz + ulst.bz;
} else if (spd[3] <= 0.0) {
// return Fr*
flxi[IDN] = fr.d + urst.d;
flxi[IVX] = fr.mx + urst.mx;
flxi[IVY] = fr.my + urst.my;
flxi[IVZ] = fr.mz + urst.mz;
flxi[IEN] = fr.e + urst.e;
flxi[IBY] = fr.by + urst.by;
flxi[IBZ] = fr.bz + urst.bz;
} else if (spd[2] >= 0.0) {
// return Fl**
flxi[IDN] = fl.d + ulst.d + uldst.d;
Expand All @@ -347,7 +351,7 @@ void Hydro::RiemannSolver(const int k, const int j, const int il, const int iu,
flxi[IEN] = fl.e + ulst.e + uldst.e;
flxi[IBY] = fl.by + ulst.by + uldst.by;
flxi[IBZ] = fl.bz + ulst.bz + uldst.bz;
} else if (spd[3] > 0.0) {
} else {
// return Fr**
flxi[IDN] = fr.d + urst.d + urdst.d;
flxi[IVX] = fr.mx + urst.mx + urdst.mx;
Expand All @@ -356,15 +360,6 @@ void Hydro::RiemannSolver(const int k, const int j, const int il, const int iu,
flxi[IEN] = fr.e + urst.e + urdst.e;
flxi[IBY] = fr.by + urst.by + urdst.by;
flxi[IBZ] = fr.bz + urst.bz + urdst.bz;
} else {
// return Fr*
flxi[IDN] = fr.d + urst.d;
flxi[IVX] = fr.mx + urst.mx;
flxi[IVY] = fr.my + urst.my;
flxi[IVZ] = fr.mz + urst.mz;
flxi[IEN] = fr.e + urst.e;
flxi[IBY] = fr.by + urst.by;
flxi[IBZ] = fr.bz + urst.bz;
}

flx(IDN,k,j,i) = flxi[IDN];
Expand Down
97 changes: 46 additions & 51 deletions src/hydro/rsolvers/mhd/lhlld.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -242,47 +242,42 @@ void Hydro::RiemannSolver(const int k, const int j, const int il, const int iu,
// (KGF): group transverse by, bz terms for floating-point associativity symmetry
urst.e = (sdr*ur.e - ptr*wri[IVX] + ptst*spd[2] +
bxi*(wri[IVX]*bxi + (wri[IVY]*ur.by + wri[IVZ]*ur.bz) - vbstr))*sdmr_inv;
// ul** and ur** - if Bx is near zero, same as *-states
if (0.5*bxsq < (SMALL_NUMBER)*ptst) {
uldst = ulst;
urdst = urst;
} else {
Real invsumd = 1.0/(sqrtdl + sqrtdr);
Real bxsig = (bxi > 0.0 ? 1.0 : -1.0);

uldst.d = ulst.d;
urdst.d = urst.d;

uldst.mx = ulst.mx;
urdst.mx = urst.mx;

// eqn (59) of M&K
Real tmp = invsumd*(sqrtdl*(ulst.my*ulst_d_inv) + sqrtdr*(urst.my*urst_d_inv) +
bxsig*(urst.by - ulst.by));
uldst.my = uldst.d * tmp;
urdst.my = urdst.d * tmp;

// eqn (60) of M&K
tmp = invsumd*(sqrtdl*(ulst.mz*ulst_d_inv) + sqrtdr*(urst.mz*urst_d_inv) +
bxsig*(urst.bz - ulst.bz));
uldst.mz = uldst.d * tmp;
urdst.mz = urdst.d * tmp;

// eqn (61) of M&K
tmp = invsumd*(sqrtdl*urst.by + sqrtdr*ulst.by +
bxsig*sqrtdl*sqrtdr*((urst.my*urst_d_inv) - (ulst.my*ulst_d_inv)));
uldst.by = urdst.by = tmp;

// eqn (62) of M&K
tmp = invsumd*(sqrtdl*urst.bz + sqrtdr*ulst.bz +
bxsig*sqrtdl*sqrtdr*((urst.mz*urst_d_inv) - (ulst.mz*ulst_d_inv)));
uldst.bz = urdst.bz = tmp;

// eqn (63) of M&K
tmp = spd[2]*bxi + (uldst.my*uldst.by + uldst.mz*uldst.bz)/uldst.d;
uldst.e = ulst.e - sqrtdl*bxsig*(vbstl - tmp);
urdst.e = urst.e + sqrtdr*bxsig*(vbstr - tmp);
}

Real invsumd = 1.0/(sqrtdl + sqrtdr);
Real bxsig = (bxi > 0.0 ? 1.0 : -1.0);

uldst.d = ulst.d;
urdst.d = urst.d;

uldst.mx = ulst.mx;
urdst.mx = urst.mx;

// eqn (59) of M&K
Real tmp = invsumd*(sqrtdl*(ulst.my*ulst_d_inv) + sqrtdr*(urst.my*urst_d_inv) +
bxsig*(urst.by - ulst.by));
uldst.my = uldst.d * tmp;
urdst.my = urdst.d * tmp;

// eqn (60) of M&K
tmp = invsumd*(sqrtdl*(ulst.mz*ulst_d_inv) + sqrtdr*(urst.mz*urst_d_inv) +
bxsig*(urst.bz - ulst.bz));
uldst.mz = uldst.d * tmp;
urdst.mz = urdst.d * tmp;

// eqn (61) of M&K
tmp = invsumd*(sqrtdl*urst.by + sqrtdr*ulst.by +
bxsig*sqrtdl*sqrtdr*((urst.my*urst_d_inv) - (ulst.my*ulst_d_inv)));
uldst.by = urdst.by = tmp;

// eqn (62) of M&K
tmp = invsumd*(sqrtdl*urst.bz + sqrtdr*ulst.bz +
bxsig*sqrtdl*sqrtdr*((urst.mz*urst_d_inv) - (ulst.mz*ulst_d_inv)));
uldst.bz = urdst.bz = tmp;

// eqn (63) of M&K
tmp = spd[2]*bxi + (uldst.my*uldst.by + uldst.mz*uldst.bz)/uldst.d;
uldst.e = ulst.e - sqrtdl*bxsig*(vbstl - tmp);
urdst.e = urst.e + sqrtdr*bxsig*(vbstr - tmp);

//--- Step 6. Compute flux
uldst.d = spd[1] * (uldst.d - ulst.d);
Expand Down Expand Up @@ -344,6 +339,15 @@ void Hydro::RiemannSolver(const int k, const int j, const int il, const int iu,
flxi[IEN] = fl.e + ulst.e;
flxi[IBY] = fl.by + ulst.by;
flxi[IBZ] = fl.bz + ulst.bz;
} else if (spd[3] <= 0.0) {
// return Fr*
flxi[IDN] = fr.d + urst.d;
flxi[IVX] = fr.mx + urst.mx;
flxi[IVY] = fr.my + urst.my;
flxi[IVZ] = fr.mz + urst.mz;
flxi[IEN] = fr.e + urst.e;
flxi[IBY] = fr.by + urst.by;
flxi[IBZ] = fr.bz + urst.bz;
} else if (spd[2] >= 0.0) {
// return Fl**
flxi[IDN] = fl.d + ulst.d + uldst.d;
Expand All @@ -353,7 +357,7 @@ void Hydro::RiemannSolver(const int k, const int j, const int il, const int iu,
flxi[IEN] = fl.e + ulst.e + uldst.e;
flxi[IBY] = fl.by + ulst.by + uldst.by;
flxi[IBZ] = fl.bz + ulst.bz + uldst.bz;
} else if (spd[3] > 0.0) {
} else {
// return Fr**
flxi[IDN] = fr.d + urst.d + urdst.d;
flxi[IVX] = fr.mx + urst.mx + urdst.mx;
Expand All @@ -362,15 +366,6 @@ void Hydro::RiemannSolver(const int k, const int j, const int il, const int iu,
flxi[IEN] = fr.e + urst.e + urdst.e;
flxi[IBY] = fr.by + urst.by + urdst.by;
flxi[IBZ] = fr.bz + urst.bz + urdst.bz;
} else {
// return Fr*
flxi[IDN] = fr.d + urst.d;
flxi[IVX] = fr.mx + urst.mx;
flxi[IVY] = fr.my + urst.my;
flxi[IVZ] = fr.mz + urst.mz;
flxi[IEN] = fr.e + urst.e;
flxi[IBY] = fr.by + urst.by;
flxi[IBZ] = fr.bz + urst.bz;
}

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

0 comments on commit ee2b2c7

Please sign in to comment.