From 17a192ffe2db6847b3ebb41cad0d4fbcc905a4cd Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Tue, 25 Jun 2024 14:10:37 -0400 Subject: [PATCH 1/5] move the HLL solvers into their own header --- Source/hydro/HLL_solvers.H | 592 +++++++++++++++++++++++++++++++++ Source/hydro/riemann.H | 140 -------- Source/hydro/riemann.cpp | 3 +- Source/hydro/riemann_solvers.H | 440 ------------------------ 4 files changed, 593 insertions(+), 582 deletions(-) create mode 100644 Source/hydro/HLL_solvers.H diff --git a/Source/hydro/HLL_solvers.H b/Source/hydro/HLL_solvers.H new file mode 100644 index 0000000000..0acf814735 --- /dev/null +++ b/Source/hydro/HLL_solvers.H @@ -0,0 +1,592 @@ +#ifndef HLL_SOLVERS_H +#define HLL_SOLVERS_H + +#include + +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE +void +HLLC_state(const int idir, const Real S_k, const Real S_c, + const Real* qstate, Real* U) { + + Real u_k = 0.0; + if (idir == 0) { + u_k = qstate[QU]; + } else if (idir == 1) { + u_k = qstate[QV]; + } else if (idir == 2) { + u_k = qstate[QW]; + } + + Real hllc_factor = qstate[QRHO]*(S_k - u_k)/(S_k - S_c); + U[URHO] = hllc_factor; + + if (idir == 0) { + U[UMX] = hllc_factor*S_c; + U[UMY] = hllc_factor*qstate[QV]; + U[UMZ] = hllc_factor*qstate[QW]; + + } else if (idir == 1) { + U[UMX] = hllc_factor*qstate[QU]; + U[UMY] = hllc_factor*S_c; + U[UMZ] = hllc_factor*qstate[QW]; + + } else { + U[UMX] = hllc_factor*qstate[QU]; + U[UMY] = hllc_factor*qstate[QV]; + U[UMZ] = hllc_factor*S_c; + } + + U[UEDEN] = hllc_factor*(qstate[QREINT]/qstate[QRHO] + + 0.5_rt*(qstate[QU]*qstate[QU] + qstate[QV]*qstate[QV] + qstate[QW]*qstate[QW]) + + (S_c - u_k)*(S_c + qstate[QPRES]/(qstate[QRHO]*(S_k - u_k)))); + U[UEINT] = hllc_factor*qstate[QREINT]/qstate[QRHO]; + + U[UTEMP] = 0.0; // we don't evolve T + +#ifdef SHOCK_VAR + U[USHK] = 0.0; +#endif + +#ifdef NSE_NET + U[UMUP] = 0.0; + U[UMUN] = 0.0; +#endif + + for (int ipassive = 0; ipassive < npassive; ipassive++) { + const int n = upassmap(ipassive); + const int nqs = qpassmap(ipassive); + U[n] = hllc_factor*qstate[nqs]; + } +} + + +/// +/// given a conserved state, compute the flux in direction idir +/// +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE +void +compute_flux(const int idir, const Real bnd_fac, const int coord, + const Real* U, const Real p, + Real* F) { + + Real u_flx = U[UMX+idir]/U[URHO]; + + if (bnd_fac == 0) { + u_flx = 0.0; + } + + F[URHO] = U[URHO]*u_flx; + + F[UMX] = U[UMX]*u_flx; + F[UMY] = U[UMY]*u_flx; + F[UMZ] = U[UMZ]*u_flx; + + auto mom_check = mom_flux_has_p(idir, idir, coord); + + if (mom_check) { + // we do not include the pressure term in any non-Cartesian + // coordinate directions + F[UMX+idir] = F[UMX+idir] + p; + } + + F[UEINT] = U[UEINT]*u_flx; + F[UEDEN] = (U[UEDEN] + p)*u_flx; + + F[UTEMP] = 0.0; + +#ifdef SHOCK_VAR + F[USHK] = 0.0; +#endif + +#ifdef NSE_NET + F[UMUP] = 0.0; + F[UMUN] = 0.0; +#endif + + for (int ipassive=0; ipassive < npassive; ipassive++) { + const int n = upassmap(ipassive); + F[n] = U[n]*u_flx; + } +} + +/// +/// convert the primitive variable state to the conserved state +/// +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE +void +cons_state(const Real* qstate, Real* U) { + + U[URHO] = qstate[QRHO]; + + // since we advect all 3 velocity components regardless of dimension, this + // will be general + U[UMX] = qstate[QRHO]*qstate[QU]; + U[UMY] = qstate[QRHO]*qstate[QV]; + U[UMZ] = qstate[QRHO]*qstate[QW]; + + U[UEDEN] = qstate[QREINT] + 0.5_rt*qstate[QRHO]*(qstate[QU]*qstate[QU] + qstate[QV]*qstate[QV] + qstate[QW]*qstate[QW]); + U[UEINT] = qstate[QREINT]; + + // we don't care about T here, but initialize it to make NaN + // checking happy + U[UTEMP] = 0.0; + +#ifdef SHOCK_VAR + U[USHK] = 0.0; +#endif + +#ifdef NSE_NET + U[UMUP] = 0.0; + U[UMUN] = 0.0; +#endif + + for (int ipassive = 0; ipassive < npassive; ipassive++) { + const int n = upassmap(ipassive); + const int nqs = qpassmap(ipassive); + U[n] = qstate[QRHO]*qstate[nqs]; + } +} + + +/// +/// A simple HLL Riemann solver for pure hydrodynamics. This takes just a +/// single interface's data and returns the HLL flux +/// +/// @param ql the left interface state +/// @param qr the right interface state +/// @param cl sound speed on the left interface +/// @param cr sound speed on the right interface +/// @param idir coordinate direction for the solve (0 = x, 1 = y, 2 = z) +/// @param coord geometry type (0 = Cartesian, 1 = axisymmetric, 2 = spherical) +/// @param f the HLL fluxes +/// +AMREX_GPU_HOST_DEVICE AMREX_INLINE +void +HLL(const Real* ql, const Real* qr, + const Real cl, const Real cr, + const int idir, const int coord, + Real* flux_hll) { + + // This is the HLLE solver. We should apply it to zone averages + // (not reconstructed states) at an interface in the presence of + // shocks to avoid the odd-even decoupling / carbuncle phenomenon. + // + // See: Einfeldt, B. et al. 1991, JCP, 92, 273 + // Einfeldt, B. 1988, SIAM J NA, 25, 294 + + + constexpr Real small_hll = 1.e-10_rt; + + int ivel, ivelt, iveltt; + int imom, imomt, imomtt; + + if (idir == 0) { + ivel = QU; + ivelt = QV; + iveltt = QW; + + imom = UMX; + imomt = UMY; + imomtt = UMZ; + + } else if (idir == 1) { + ivel = QV; + ivelt = QU; + iveltt = QW; + + imom = UMY; + imomt = UMX; + imomtt = UMZ; + + } else { + ivel = QW; + ivelt = QU; + iveltt = QV; + + imom = UMZ; + imomt = UMX; + imomtt = UMY; + } + + Real rhol_sqrt = std::sqrt(ql[QRHO]); + Real rhor_sqrt = std::sqrt(qr[QRHO]); + + Real rhod = 1.0_rt/(rhol_sqrt + rhor_sqrt); + + + // compute the average sound speed. This uses an approximation from + // E88, eq. 5.6, 5.7 that assumes gamma falls between 1 + // and 5/3 + Real cavg = std::sqrt( (rhol_sqrt*cl*cl + rhor_sqrt*cr*cr)*rhod + + 0.5_rt*rhol_sqrt*rhor_sqrt*rhod*rhod*std::pow(qr[ivel] - ql[ivel], 2)); + + + // Roe eigenvalues (E91, eq. 5.3b) + Real uavg = (rhol_sqrt*ql[ivel] + rhor_sqrt*qr[ivel])*rhod; + + Real a1 = uavg - cavg; + Real a4 = uavg + cavg; + + + // signal speeds (E91, eq. 4.5) + Real bl = amrex::min(a1, ql[ivel] - cl); + Real br = amrex::max(a4, qr[ivel] + cr); + + Real bm = amrex::min(0.0_rt, bl); + Real bp = amrex::max(0.0_rt, br); + + Real bd = bp - bm; + + if (std::abs(bd) < small_hll*amrex::max(std::abs(bm), std::abs(bp))) { + return; + } + + // we'll overwrite the passed in flux with the HLL flux + + bd = 1.0_rt/bd; + + // compute the fluxes according to E91, eq. 4.4b -- note that the + // min/max above picks the correct flux if we are not in the star + // region + + // density flux + Real fl_tmp = ql[QRHO]*ql[ivel]; + Real fr_tmp = qr[QRHO]*qr[ivel]; + + flux_hll[URHO] = (bp*fl_tmp - bm*fr_tmp)*bd + bp*bm*bd*(qr[QRHO] - ql[QRHO]); + + // normal momentum flux. Note for 1-d and 2-d non cartesian + // r-coordinate, we leave off the pressure term and handle that + // separately in the update, to accommodate different geometries + fl_tmp = ql[QRHO]*ql[ivel]*ql[ivel]; + fr_tmp = qr[QRHO]*qr[ivel]*qr[ivel]; + if (mom_flux_has_p(idir, idir, coord)) { + fl_tmp = fl_tmp + ql[QPRES]; + fr_tmp = fr_tmp + qr[QPRES]; + } + + flux_hll[imom] = (bp*fl_tmp - bm*fr_tmp)*bd + bp*bm*bd*(qr[QRHO]*qr[ivel] - ql[QRHO]*ql[ivel]); + + // transverse momentum flux + fl_tmp = ql[QRHO]*ql[ivel]*ql[ivelt]; + fr_tmp = qr[QRHO]*qr[ivel]*qr[ivelt]; + + flux_hll[imomt] = (bp*fl_tmp - bm*fr_tmp)*bd + bp*bm*bd*(qr[QRHO]*qr[ivelt] - ql[QRHO]*ql[ivelt]); + + + fl_tmp = ql[QRHO]*ql[ivel]*ql[iveltt]; + fr_tmp = qr[QRHO]*qr[ivel]*qr[iveltt]; + + flux_hll[imomtt] = (bp*fl_tmp - bm*fr_tmp)*bd + bp*bm*bd*(qr[QRHO]*qr[iveltt] - ql[QRHO]*ql[iveltt]); + + // total energy flux + Real rhoEl = ql[QREINT] + 0.5_rt*ql[QRHO]*(ql[ivel]*ql[ivel] + ql[ivelt]*ql[ivelt] + ql[iveltt]*ql[iveltt]); + fl_tmp = ql[ivel]*(rhoEl + ql[QPRES]); + + Real rhoEr = qr[QREINT] + 0.5_rt*qr[QRHO]*(qr[ivel]*qr[ivel] + qr[ivelt]*qr[ivelt] + qr[iveltt]*qr[iveltt]); + fr_tmp = qr[ivel]*(rhoEr + qr[QPRES]); + + flux_hll[UEDEN] = (bp*fl_tmp - bm*fr_tmp)*bd + bp*bm*bd*(rhoEr - rhoEl); + + + // eint flux + fl_tmp = ql[QREINT]*ql[ivel]; + fr_tmp = qr[QREINT]*qr[ivel]; + + flux_hll[UEINT] = (bp*fl_tmp - bm*fr_tmp)*bd + bp*bm*bd*(qr[QREINT] - ql[QREINT]); + + + // passively-advected scalar fluxes + for (int ipassive = 0; ipassive < npassive; ipassive++) { + const int n = upassmap(ipassive); + const int nqs = qpassmap(ipassive); + + fl_tmp = ql[QRHO]*ql[nqs]*ql[ivel]; + fr_tmp = qr[QRHO]*qr[nqs]*qr[ivel]; + + flux_hll[n] = (bp*fl_tmp - bm*fr_tmp)*bd + bp*bm*bd*(qr[QRHO]*qr[nqs] - ql[QRHO]*ql[nqs]); + } +} + + +/// +/// A HLLC Riemann solver for pure hydrodynamics +/// +/// @param ql the left interface state +/// @param qr the right interface state +/// @param qaux_arr the auxiliary state +/// @param uflx the flux through the interface +/// @param qint an approximate Godunov state on the interface +/// @param idir coordinate direction for the solve (0 = x, 1 = y, 2 = z) +/// +AMREX_GPU_HOST_DEVICE AMREX_INLINE +void +HLLC(const int i, const int j, const int k, const int idir, + Array4 const& ql, + Array4 const& qr, + Array4 const& qaux_arr, + Array4 const& uflx, + Array4 const& qgdnv, const bool store_full_state, + const GeometryData& geom, + const bool special_bnd_lo, const bool special_bnd_hi, + GpuArray const& domlo, GpuArray const& domhi) { + + // this is an implementation of the HLLC solver described in Toro's + // book. it uses the simplest estimate of the wave speeds, since + // those should work for a general EOS. We also initially do the + // CGF Riemann construction to get pstar and ustar, since we'll + // need to know the pressure and velocity on the interface for the + // pdV term in the internal energy update. + + const Real small = 1.e-8_rt; + + int iu; + int sx, sy, sz; + + if (idir == 0) { + iu = QU; + sx = 1; + sy = 0; + sz = 0; + + } else if (idir == 1) { + iu = QV; + sx = 0; + sy = 1; + sz = 0; + + } else { + iu = QW; + sx = 0; + sy = 0; + sz = 1; + + } + + + int coord = geom.Coord(); + + // deal with hard walls + Real bnd_fac = 1.0_rt; + + if (idir == 0) { + if ((i == domlo[0] && special_bnd_lo) || + (i == domhi[0]+1 && special_bnd_hi)) { + bnd_fac = 0.0_rt; + } + + } else if (idir == 1) { + if ((j == domlo[1] && special_bnd_lo) || + (j == domhi[1]+1 && special_bnd_hi)) { + bnd_fac = 0.0_rt; + } + + } else { + if ((k == domlo[2] && special_bnd_lo) || + (k == domhi[2]+1 && special_bnd_hi)) { + bnd_fac = 0.0_rt; + } + } + + + Real rl = amrex::max(ql(i,j,k,QRHO), small_dens); + + // pick left velocities based on direction + Real ul = ql(i,j,k,iu); + + Real pl = amrex::max(ql(i,j,k,QPRES), small_pres); + + Real rr = amrex::max(qr(i,j,k,QRHO), small_dens); + + // pick right velocities based on direction + Real ur = qr(i,j,k,iu); + + Real pr = amrex::max(qr(i,j,k,QPRES), small_pres); + + // now we essentially do the CGF solver to get p and u on the + // interface, but we won't use these in any flux construction. + Real csmall = amrex::max(small, amrex::max(small * qaux_arr(i,j,k,QC), + small * qaux_arr(i-sx,j-sy,k-sz,QC))); + Real cavg = 0.5_rt*(qaux_arr(i,j,k,QC) + qaux_arr(i-sx,j-sy,k-sz,QC)); + + Real gamcl = qaux_arr(i-sx,j-sy,k-sz,QGAMC); + Real gamcr = qaux_arr(i,j,k,QGAMC); + +#ifdef TRUE_SDC + if (use_reconstructed_gamma1 == 1) { + gamcl = ql(i,j,k,QGC); + gamcr = qr(i,j,k,QGC); + } +#endif + + Real wsmall = small_dens*csmall; + Real wl = amrex::max(wsmall, std::sqrt(std::abs(gamcl*pl*rl))); + Real wr = amrex::max(wsmall, std::sqrt(std::abs(gamcr*pr*rr))); + + Real wwinv = 1.0_rt/(wl + wr); + Real pstar = ((wr*pl + wl*pr) + wl*wr*(ul - ur))*wwinv; + Real ustar = ((wl*ul + wr*ur) + (pl - pr))*wwinv; + + pstar = amrex::max(pstar, small_pres); + // for symmetry preservation, if ustar is really small, then we + // set it to zero + if (std::abs(ustar) < riemann_constants::smallu*0.5_rt*(std::abs(ul) + std::abs(ur))){ + ustar = 0.0_rt; + } + + Real ro; + Real uo; + Real po; + Real gamco; + + if (ustar > 0.0_rt) { + ro = rl; + uo = ul; + po = pl; + gamco = gamcl; + + } else if (ustar < 0.0_rt) { + ro = rr; + uo = ur; + po = pr; + gamco = gamcr; + + } else { + ro = 0.5_rt*(rl + rr); + uo = 0.5_rt*(ul + ur); + po = 0.5_rt*(pl + pr); + gamco = 0.5_rt*(gamcl + gamcr); + } + + ro = amrex::max(small_dens, ro); + + Real roinv = 1.0_rt/ro; + Real co = std::sqrt(std::abs(gamco*po*roinv)); + co = amrex::max(csmall, co); + Real co2inv = 1.0_rt/(co*co); + + Real rstar = ro + (pstar - po)*co2inv; + rstar = amrex::max(small_dens, rstar); + + Real cstar = std::sqrt(std::abs(gamco*pstar/rstar)); + cstar = max(cstar, csmall); + + Real sgnm = std::copysign(1.0_rt, ustar); + Real spout = co - sgnm*uo; + Real spin = cstar - sgnm*ustar; + Real ushock = 0.5_rt*(spin + spout); + + if (pstar-po > 0.0_rt) { + spin = ushock; + spout = ushock; + } + + Real scr = spout-spin; + if (spout-spin == 0.0_rt) { + scr = small*cavg; + } + + Real frac = (1.0_rt + (spout + spin)/scr)*0.5_rt; + frac = amrex::max(0.0_rt, amrex::min(1.0_rt, frac)); + + Real qint[NQ] = {0.0}; + + qint[QRHO] = frac*rstar + (1.0_rt - frac)*ro; + qint[iu] = frac*ustar + (1.0_rt - frac)*uo; + qint[QPRES] = frac*pstar + (1.0_rt - frac)*po; + + + // now we do the HLLC construction + + // use the simplest estimates of the wave speeds + Real S_l = amrex::min(ul - std::sqrt(gamcl*pl/rl), ur - std::sqrt(gamcr*pr/rr)); + Real S_r = amrex::max(ul + std::sqrt(gamcl*pl/rl), ur + std::sqrt(gamcr*pr/rr)); + + // estimate of the contact speed -- this is Toro Eq. 10.8 + Real S_c = (pr - pl + rl*ul*(S_l - ul) - rr*ur*(S_r - ur))/ + (rl*(S_l - ul) - rr*(S_r - ur)); + + Real q_zone[NQ]; + Real U_state[NUM_STATE]; + Real U_hllc_state[NUM_STATE]; + Real F_state[NUM_STATE]; + + if (S_r <= 0.0_rt) { + // R region + for (int n = 0; n < NQ; n++) { + q_zone[n] = qr(i,j,k,n); + } + cons_state(q_zone, U_state); + compute_flux(idir, bnd_fac, coord, + U_state, pr, F_state); + + } else if (S_r > 0.0_rt && S_c <= 0.0_rt) { + // R* region + for (int n = 0; n < NQ; n++) { + q_zone[n] = qr(i,j,k,n); + } + cons_state(q_zone, U_state); + compute_flux(idir, bnd_fac, coord, + U_state, pr, F_state); + HLLC_state(idir, S_r, S_c, q_zone, U_hllc_state); + + // correct the flux + for (int n = 0; n < NUM_STATE; n++) { + F_state[n] = F_state[n] + S_r*(U_hllc_state[n] - U_state[n]); + } + + } else if (S_c > 0.0_rt && S_l < 0.0_rt) { + // L* region + for (int n = 0; n < NQ; n++) { + q_zone[n] = ql(i,j,k,n); + } + cons_state(q_zone, U_state); + compute_flux(idir, bnd_fac, coord, + U_state, pl, F_state); + HLLC_state(idir, S_l, S_c, q_zone, U_hllc_state); + + // correct the flux + for (int n = 0; n < NUM_STATE; n++) { + F_state[n] = F_state[n] + S_l*(U_hllc_state[n] - U_state[n]); + } + + } else { + // L region + for (int n = 0; n < NQ; n++) { + q_zone[n] = ql(i,j,k,n); + } + cons_state(q_zone, U_state); + compute_flux(idir, bnd_fac, coord, + U_state, pl, F_state); + } + + for (int n = 0; n < NUM_STATE; n++) { + uflx(i,j,k,n) = F_state[n]; + } + + + // now store the state + + if (store_full_state) { + + for (int n = 0; n < NQ; n++) { + qgdnv(i,j,k,n) = qint[n]; + } + + } else { + +#ifdef HYBRID_MOMENTUM + qgdnv(i,j,k,GDRHO) = qint[QRHO]; +#endif + qgdnv(i,j,k,GDU) = qint[QU]; + qgdnv(i,j,k,GDV) = qint[QV]; + qgdnv(i,j,k,GDW) = qint[QW]; + qgdnv(i,j,k,GDPRES) = qint[QPRES]; + + } + + +} + + +#endif diff --git a/Source/hydro/riemann.H b/Source/hydro/riemann.H index c62c7a1647..47166f060c 100644 --- a/Source/hydro/riemann.H +++ b/Source/hydro/riemann.H @@ -375,144 +375,4 @@ pstar_bisection(Real& pstar_lo, Real& pstar_hi, pstar = pstar_c; } - -AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -void -cons_state(const Real* qstate, Real* U) { - - U[URHO] = qstate[QRHO]; - - // since we advect all 3 velocity components regardless of dimension, this - // will be general - U[UMX] = qstate[QRHO]*qstate[QU]; - U[UMY] = qstate[QRHO]*qstate[QV]; - U[UMZ] = qstate[QRHO]*qstate[QW]; - - U[UEDEN] = qstate[QREINT] + 0.5_rt*qstate[QRHO]*(qstate[QU]*qstate[QU] + qstate[QV]*qstate[QV] + qstate[QW]*qstate[QW]); - U[UEINT] = qstate[QREINT]; - - // we don't care about T here, but initialize it to make NaN - // checking happy - U[UTEMP] = 0.0; - -#ifdef SHOCK_VAR - U[USHK] = 0.0; -#endif - -#ifdef NSE_NET - U[UMUP] = 0.0; - U[UMUN] = 0.0; -#endif - - for (int ipassive = 0; ipassive < npassive; ipassive++) { - const int n = upassmap(ipassive); - const int nqs = qpassmap(ipassive); - U[n] = qstate[QRHO]*qstate[nqs]; - } -} - -AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -void -HLLC_state(const int idir, const Real S_k, const Real S_c, - const Real* qstate, Real* U) { - - Real u_k = 0.0; - if (idir == 0) { - u_k = qstate[QU]; - } else if (idir == 1) { - u_k = qstate[QV]; - } else if (idir == 2) { - u_k = qstate[QW]; - } - - Real hllc_factor = qstate[QRHO]*(S_k - u_k)/(S_k - S_c); - U[URHO] = hllc_factor; - - if (idir == 0) { - U[UMX] = hllc_factor*S_c; - U[UMY] = hllc_factor*qstate[QV]; - U[UMZ] = hllc_factor*qstate[QW]; - - } else if (idir == 1) { - U[UMX] = hllc_factor*qstate[QU]; - U[UMY] = hllc_factor*S_c; - U[UMZ] = hllc_factor*qstate[QW]; - - } else { - U[UMX] = hllc_factor*qstate[QU]; - U[UMY] = hllc_factor*qstate[QV]; - U[UMZ] = hllc_factor*S_c; - } - - U[UEDEN] = hllc_factor*(qstate[QREINT]/qstate[QRHO] + - 0.5_rt*(qstate[QU]*qstate[QU] + qstate[QV]*qstate[QV] + qstate[QW]*qstate[QW]) + - (S_c - u_k)*(S_c + qstate[QPRES]/(qstate[QRHO]*(S_k - u_k)))); - U[UEINT] = hllc_factor*qstate[QREINT]/qstate[QRHO]; - - U[UTEMP] = 0.0; // we don't evolve T - -#ifdef SHOCK_VAR - U[USHK] = 0.0; -#endif - -#ifdef NSE_NET - U[UMUP] = 0.0; - U[UMUN] = 0.0; -#endif - - for (int ipassive = 0; ipassive < npassive; ipassive++) { - const int n = upassmap(ipassive); - const int nqs = qpassmap(ipassive); - U[n] = hllc_factor*qstate[nqs]; - } -} - -AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -void -compute_flux(const int idir, const Real bnd_fac, const int coord, - const Real* U, const Real p, - Real* F) { - - // given a conserved state, compute the flux in direction idir - - Real u_flx = U[UMX+idir]/U[URHO]; - - if (bnd_fac == 0) { - u_flx = 0.0; - } - - F[URHO] = U[URHO]*u_flx; - - F[UMX] = U[UMX]*u_flx; - F[UMY] = U[UMY]*u_flx; - F[UMZ] = U[UMZ]*u_flx; - - auto mom_check = mom_flux_has_p(idir, idir, coord); - - if (mom_check) { - // we do not include the pressure term in any non-Cartesian - // coordinate directions - F[UMX+idir] = F[UMX+idir] + p; - } - - F[UEINT] = U[UEINT]*u_flx; - F[UEDEN] = (U[UEDEN] + p)*u_flx; - - F[UTEMP] = 0.0; - -#ifdef SHOCK_VAR - F[USHK] = 0.0; -#endif - -#ifdef NSE_NET - F[UMUP] = 0.0; - F[UMUN] = 0.0; -#endif - - for (int ipassive=0; ipassive < npassive; ipassive++) { - const int n = upassmap(ipassive); - F[n] = U[n]*u_flx; - } -} - #endif diff --git a/Source/hydro/riemann.cpp b/Source/hydro/riemann.cpp index db8641a52b..f674b954a8 100644 --- a/Source/hydro/riemann.cpp +++ b/Source/hydro/riemann.cpp @@ -1,6 +1,7 @@ #include #include +#include #ifdef RADIATION #include @@ -193,5 +194,3 @@ Castro::cmpflx_plus_godunov(const Box& bx, }); } - - diff --git a/Source/hydro/riemann_solvers.H b/Source/hydro/riemann_solvers.H index e83f0488e7..51a804e380 100644 --- a/Source/hydro/riemann_solvers.H +++ b/Source/hydro/riemann_solvers.H @@ -819,446 +819,6 @@ riemannus(const RiemannState& ql, const RiemannState& qr, const RiemannAux& raux } -/// -/// A simple HLL Riemann solver for pure hydrodynamics. This takes just a -/// single interface's data and returns the HLL flux -/// -/// @param ql the left interface state -/// @param qr the right interface state -/// @param cl sound speed on the left interface -/// @param cr sound speed on the right interface -/// @param idir coordinate direction for the solve (0 = x, 1 = y, 2 = z) -/// @param coord geometry type (0 = Cartesian, 1 = axisymmetric, 2 = spherical) -/// @param f the HLL fluxes -/// -AMREX_GPU_HOST_DEVICE AMREX_INLINE -void -HLL(const Real* ql, const Real* qr, - const Real cl, const Real cr, - const int idir, const int coord, - Real* flux_hll) { - - // This is the HLLE solver. We should apply it to zone averages - // (not reconstructed states) at an interface in the presence of - // shocks to avoid the odd-even decoupling / carbuncle phenomenon. - // - // See: Einfeldt, B. et al. 1991, JCP, 92, 273 - // Einfeldt, B. 1988, SIAM J NA, 25, 294 - - - constexpr Real small_hll = 1.e-10_rt; - - int ivel, ivelt, iveltt; - int imom, imomt, imomtt; - - if (idir == 0) { - ivel = QU; - ivelt = QV; - iveltt = QW; - - imom = UMX; - imomt = UMY; - imomtt = UMZ; - - } else if (idir == 1) { - ivel = QV; - ivelt = QU; - iveltt = QW; - - imom = UMY; - imomt = UMX; - imomtt = UMZ; - - } else { - ivel = QW; - ivelt = QU; - iveltt = QV; - - imom = UMZ; - imomt = UMX; - imomtt = UMY; - } - - Real rhol_sqrt = std::sqrt(ql[QRHO]); - Real rhor_sqrt = std::sqrt(qr[QRHO]); - - Real rhod = 1.0_rt/(rhol_sqrt + rhor_sqrt); - - - // compute the average sound speed. This uses an approximation from - // E88, eq. 5.6, 5.7 that assumes gamma falls between 1 - // and 5/3 - Real cavg = std::sqrt( (rhol_sqrt*cl*cl + rhor_sqrt*cr*cr)*rhod + - 0.5_rt*rhol_sqrt*rhor_sqrt*rhod*rhod*std::pow(qr[ivel] - ql[ivel], 2)); - - - // Roe eigenvalues (E91, eq. 5.3b) - Real uavg = (rhol_sqrt*ql[ivel] + rhor_sqrt*qr[ivel])*rhod; - - Real a1 = uavg - cavg; - Real a4 = uavg + cavg; - - - // signal speeds (E91, eq. 4.5) - Real bl = amrex::min(a1, ql[ivel] - cl); - Real br = amrex::max(a4, qr[ivel] + cr); - - Real bm = amrex::min(0.0_rt, bl); - Real bp = amrex::max(0.0_rt, br); - - Real bd = bp - bm; - - if (std::abs(bd) < small_hll*amrex::max(std::abs(bm), std::abs(bp))) { - return; - } - - // we'll overwrite the passed in flux with the HLL flux - - bd = 1.0_rt/bd; - - // compute the fluxes according to E91, eq. 4.4b -- note that the - // min/max above picks the correct flux if we are not in the star - // region - - // density flux - Real fl_tmp = ql[QRHO]*ql[ivel]; - Real fr_tmp = qr[QRHO]*qr[ivel]; - - flux_hll[URHO] = (bp*fl_tmp - bm*fr_tmp)*bd + bp*bm*bd*(qr[QRHO] - ql[QRHO]); - - // normal momentum flux. Note for 1-d and 2-d non cartesian - // r-coordinate, we leave off the pressure term and handle that - // separately in the update, to accommodate different geometries - fl_tmp = ql[QRHO]*ql[ivel]*ql[ivel]; - fr_tmp = qr[QRHO]*qr[ivel]*qr[ivel]; - if (mom_flux_has_p(idir, idir, coord)) { - fl_tmp = fl_tmp + ql[QPRES]; - fr_tmp = fr_tmp + qr[QPRES]; - } - - flux_hll[imom] = (bp*fl_tmp - bm*fr_tmp)*bd + bp*bm*bd*(qr[QRHO]*qr[ivel] - ql[QRHO]*ql[ivel]); - - // transverse momentum flux - fl_tmp = ql[QRHO]*ql[ivel]*ql[ivelt]; - fr_tmp = qr[QRHO]*qr[ivel]*qr[ivelt]; - - flux_hll[imomt] = (bp*fl_tmp - bm*fr_tmp)*bd + bp*bm*bd*(qr[QRHO]*qr[ivelt] - ql[QRHO]*ql[ivelt]); - - - fl_tmp = ql[QRHO]*ql[ivel]*ql[iveltt]; - fr_tmp = qr[QRHO]*qr[ivel]*qr[iveltt]; - - flux_hll[imomtt] = (bp*fl_tmp - bm*fr_tmp)*bd + bp*bm*bd*(qr[QRHO]*qr[iveltt] - ql[QRHO]*ql[iveltt]); - - // total energy flux - Real rhoEl = ql[QREINT] + 0.5_rt*ql[QRHO]*(ql[ivel]*ql[ivel] + ql[ivelt]*ql[ivelt] + ql[iveltt]*ql[iveltt]); - fl_tmp = ql[ivel]*(rhoEl + ql[QPRES]); - - Real rhoEr = qr[QREINT] + 0.5_rt*qr[QRHO]*(qr[ivel]*qr[ivel] + qr[ivelt]*qr[ivelt] + qr[iveltt]*qr[iveltt]); - fr_tmp = qr[ivel]*(rhoEr + qr[QPRES]); - - flux_hll[UEDEN] = (bp*fl_tmp - bm*fr_tmp)*bd + bp*bm*bd*(rhoEr - rhoEl); - - - // eint flux - fl_tmp = ql[QREINT]*ql[ivel]; - fr_tmp = qr[QREINT]*qr[ivel]; - - flux_hll[UEINT] = (bp*fl_tmp - bm*fr_tmp)*bd + bp*bm*bd*(qr[QREINT] - ql[QREINT]); - - - // passively-advected scalar fluxes - for (int ipassive = 0; ipassive < npassive; ipassive++) { - const int n = upassmap(ipassive); - const int nqs = qpassmap(ipassive); - - fl_tmp = ql[QRHO]*ql[nqs]*ql[ivel]; - fr_tmp = qr[QRHO]*qr[nqs]*qr[ivel]; - - flux_hll[n] = (bp*fl_tmp - bm*fr_tmp)*bd + bp*bm*bd*(qr[QRHO]*qr[nqs] - ql[QRHO]*ql[nqs]); - } -} - - -/// -/// A HLLC Riemann solver for pure hydrodynamics -/// -/// @param ql the left interface state -/// @param qr the right interface state -/// @param qaux_arr the auxiliary state -/// @param uflx the flux through the interface -/// @param qint an approximate Godunov state on the interface -/// @param idir coordinate direction for the solve (0 = x, 1 = y, 2 = z) -/// -AMREX_GPU_HOST_DEVICE AMREX_INLINE -void -HLLC(const int i, const int j, const int k, const int idir, - Array4 const& ql, - Array4 const& qr, - Array4 const& qaux_arr, - Array4 const& uflx, - Array4 const& qgdnv, const bool store_full_state, - const GeometryData& geom, - const bool special_bnd_lo, const bool special_bnd_hi, - GpuArray const& domlo, GpuArray const& domhi) { - - // this is an implementation of the HLLC solver described in Toro's - // book. it uses the simplest estimate of the wave speeds, since - // those should work for a general EOS. We also initially do the - // CGF Riemann construction to get pstar and ustar, since we'll - // need to know the pressure and velocity on the interface for the - // pdV term in the internal energy update. - - const Real small = 1.e-8_rt; - - int iu; - int sx, sy, sz; - - if (idir == 0) { - iu = QU; - sx = 1; - sy = 0; - sz = 0; - - } else if (idir == 1) { - iu = QV; - sx = 0; - sy = 1; - sz = 0; - - } else { - iu = QW; - sx = 0; - sy = 0; - sz = 1; - - } - - - int coord = geom.Coord(); - - // deal with hard walls - Real bnd_fac = 1.0_rt; - - if (idir == 0) { - if ((i == domlo[0] && special_bnd_lo) || - (i == domhi[0]+1 && special_bnd_hi)) { - bnd_fac = 0.0_rt; - } - - } else if (idir == 1) { - if ((j == domlo[1] && special_bnd_lo) || - (j == domhi[1]+1 && special_bnd_hi)) { - bnd_fac = 0.0_rt; - } - - } else { - if ((k == domlo[2] && special_bnd_lo) || - (k == domhi[2]+1 && special_bnd_hi)) { - bnd_fac = 0.0_rt; - } - } - - - Real rl = amrex::max(ql(i,j,k,QRHO), small_dens); - - // pick left velocities based on direction - Real ul = ql(i,j,k,iu); - - Real pl = amrex::max(ql(i,j,k,QPRES), small_pres); - - Real rr = amrex::max(qr(i,j,k,QRHO), small_dens); - - // pick right velocities based on direction - Real ur = qr(i,j,k,iu); - - Real pr = amrex::max(qr(i,j,k,QPRES), small_pres); - - // now we essentially do the CGF solver to get p and u on the - // interface, but we won't use these in any flux construction. - Real csmall = amrex::max(small, amrex::max(small * qaux_arr(i,j,k,QC), - small * qaux_arr(i-sx,j-sy,k-sz,QC))); - Real cavg = 0.5_rt*(qaux_arr(i,j,k,QC) + qaux_arr(i-sx,j-sy,k-sz,QC)); - - Real gamcl = qaux_arr(i-sx,j-sy,k-sz,QGAMC); - Real gamcr = qaux_arr(i,j,k,QGAMC); - -#ifdef TRUE_SDC - if (use_reconstructed_gamma1 == 1) { - gamcl = ql(i,j,k,QGC); - gamcr = qr(i,j,k,QGC); - } -#endif - - Real wsmall = small_dens*csmall; - Real wl = amrex::max(wsmall, std::sqrt(std::abs(gamcl*pl*rl))); - Real wr = amrex::max(wsmall, std::sqrt(std::abs(gamcr*pr*rr))); - - Real wwinv = 1.0_rt/(wl + wr); - Real pstar = ((wr*pl + wl*pr) + wl*wr*(ul - ur))*wwinv; - Real ustar = ((wl*ul + wr*ur) + (pl - pr))*wwinv; - - pstar = amrex::max(pstar, small_pres); - // for symmetry preservation, if ustar is really small, then we - // set it to zero - if (std::abs(ustar) < riemann_constants::smallu*0.5_rt*(std::abs(ul) + std::abs(ur))){ - ustar = 0.0_rt; - } - - Real ro; - Real uo; - Real po; - Real gamco; - - if (ustar > 0.0_rt) { - ro = rl; - uo = ul; - po = pl; - gamco = gamcl; - - } else if (ustar < 0.0_rt) { - ro = rr; - uo = ur; - po = pr; - gamco = gamcr; - - } else { - ro = 0.5_rt*(rl + rr); - uo = 0.5_rt*(ul + ur); - po = 0.5_rt*(pl + pr); - gamco = 0.5_rt*(gamcl + gamcr); - } - - ro = amrex::max(small_dens, ro); - - Real roinv = 1.0_rt/ro; - Real co = std::sqrt(std::abs(gamco*po*roinv)); - co = amrex::max(csmall, co); - Real co2inv = 1.0_rt/(co*co); - - Real rstar = ro + (pstar - po)*co2inv; - rstar = amrex::max(small_dens, rstar); - - Real cstar = std::sqrt(std::abs(gamco*pstar/rstar)); - cstar = max(cstar, csmall); - - Real sgnm = std::copysign(1.0_rt, ustar); - Real spout = co - sgnm*uo; - Real spin = cstar - sgnm*ustar; - Real ushock = 0.5_rt*(spin + spout); - - if (pstar-po > 0.0_rt) { - spin = ushock; - spout = ushock; - } - - Real scr = spout-spin; - if (spout-spin == 0.0_rt) { - scr = small*cavg; - } - - Real frac = (1.0_rt + (spout + spin)/scr)*0.5_rt; - frac = amrex::max(0.0_rt, amrex::min(1.0_rt, frac)); - - Real qint[NQ] = {0.0}; - - qint[QRHO] = frac*rstar + (1.0_rt - frac)*ro; - qint[iu] = frac*ustar + (1.0_rt - frac)*uo; - qint[QPRES] = frac*pstar + (1.0_rt - frac)*po; - - - // now we do the HLLC construction - - // use the simplest estimates of the wave speeds - Real S_l = amrex::min(ul - std::sqrt(gamcl*pl/rl), ur - std::sqrt(gamcr*pr/rr)); - Real S_r = amrex::max(ul + std::sqrt(gamcl*pl/rl), ur + std::sqrt(gamcr*pr/rr)); - - // estimate of the contact speed -- this is Toro Eq. 10.8 - Real S_c = (pr - pl + rl*ul*(S_l - ul) - rr*ur*(S_r - ur))/ - (rl*(S_l - ul) - rr*(S_r - ur)); - - Real q_zone[NQ]; - Real U_state[NUM_STATE]; - Real U_hllc_state[NUM_STATE]; - Real F_state[NUM_STATE]; - - if (S_r <= 0.0_rt) { - // R region - for (int n = 0; n < NQ; n++) { - q_zone[n] = qr(i,j,k,n); - } - cons_state(q_zone, U_state); - compute_flux(idir, bnd_fac, coord, - U_state, pr, F_state); - - } else if (S_r > 0.0_rt && S_c <= 0.0_rt) { - // R* region - for (int n = 0; n < NQ; n++) { - q_zone[n] = qr(i,j,k,n); - } - cons_state(q_zone, U_state); - compute_flux(idir, bnd_fac, coord, - U_state, pr, F_state); - HLLC_state(idir, S_r, S_c, q_zone, U_hllc_state); - - // correct the flux - for (int n = 0; n < NUM_STATE; n++) { - F_state[n] = F_state[n] + S_r*(U_hllc_state[n] - U_state[n]); - } - - } else if (S_c > 0.0_rt && S_l < 0.0_rt) { - // L* region - for (int n = 0; n < NQ; n++) { - q_zone[n] = ql(i,j,k,n); - } - cons_state(q_zone, U_state); - compute_flux(idir, bnd_fac, coord, - U_state, pl, F_state); - HLLC_state(idir, S_l, S_c, q_zone, U_hllc_state); - - // correct the flux - for (int n = 0; n < NUM_STATE; n++) { - F_state[n] = F_state[n] + S_l*(U_hllc_state[n] - U_state[n]); - } - - } else { - // L region - for (int n = 0; n < NQ; n++) { - q_zone[n] = ql(i,j,k,n); - } - cons_state(q_zone, U_state); - compute_flux(idir, bnd_fac, coord, - U_state, pl, F_state); - } - - for (int n = 0; n < NUM_STATE; n++) { - uflx(i,j,k,n) = F_state[n]; - } - - - // now store the state - - if (store_full_state) { - - for (int n = 0; n < NQ; n++) { - qgdnv(i,j,k,n) = qint[n]; - } - - } else { - -#ifdef HYBRID_MOMENTUM - qgdnv(i,j,k,GDRHO) = qint[QRHO]; -#endif - qgdnv(i,j,k,GDU) = qint[QU]; - qgdnv(i,j,k,GDV) = qint[QV]; - qgdnv(i,j,k,GDW) = qint[QW]; - qgdnv(i,j,k,GDPRES) = qint[QPRES]; - - } - - -} - AMREX_GPU_HOST_DEVICE AMREX_INLINE From a3a13c7b93a1d5f8d1b96181952a947b9839693a Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Tue, 25 Jun 2024 14:51:00 -0400 Subject: [PATCH 2/5] split out the 2shock stuff --- Source/hydro/Make.package | 1 + Source/hydro/riemann.H | 177 +----- Source/hydro/riemann_2shock_solvers.H | 751 ++++++++++++++++++++++++++ Source/hydro/riemann_solvers.H | 607 +-------------------- 4 files changed, 778 insertions(+), 758 deletions(-) create mode 100644 Source/hydro/riemann_2shock_solvers.H diff --git a/Source/hydro/Make.package b/Source/hydro/Make.package index c867ef538d..d176be31e3 100644 --- a/Source/hydro/Make.package +++ b/Source/hydro/Make.package @@ -28,6 +28,7 @@ endif CEXE_headers += ppm.H CEXE_sources += riemann.cpp CEXE_headers += riemann_solvers.H +CEXE_headers += riemann_2shock_solvers.H CEXE_sources += riemann_util.cpp CEXE_headers += riemann.H CEXE_headers += slope.H diff --git a/Source/hydro/riemann.H b/Source/hydro/riemann.H index 47166f060c..c8be48f47a 100644 --- a/Source/hydro/riemann.H +++ b/Source/hydro/riemann.H @@ -1,31 +1,33 @@ #ifndef CASTRO_RIEMANN_H #define CASTRO_RIEMANN_H -using namespace amrex; +#include + +using namespace amrex::literals; namespace riemann_constants { - const Real smlp1 = 1.e-10_rt; - const Real small = 1.e-8_rt; - const Real smallu = 1.e-12_rt; + const amrex::Real smlp1 = 1.e-10_rt; + const amrex::Real small = 1.e-8_rt; + const amrex::Real smallu = 1.e-12_rt; } struct RiemannState { - Real rho; - Real p; - Real rhoe; - Real gamc; - Real un; - Real ut; - Real utt; + amrex::Real rho; + amrex::Real p; + amrex::Real rhoe; + amrex::Real gamc; + amrex::Real un; + amrex::Real ut; + amrex::Real utt; #ifdef RADIATION - Real lam[NGROUPS]; - Real er[NGROUPS]; - Real p_g; - Real rhoe_g; - Real gamcg; + amrex::Real lam[NGROUPS]; + amrex::Real er[NGROUPS]; + amrex::Real p_g; + amrex::Real rhoe_g; + amrex::Real gamcg; #endif }; @@ -33,9 +35,9 @@ struct RiemannState struct RiemannAux { - Real csmall; - Real cavg; - Real bnd_fac; + amrex::Real csmall; + amrex::Real cavg; + amrex::Real bnd_fac; }; @@ -66,12 +68,12 @@ std::ostream& operator<< (std::ostream& o, RiemannAux const& ra) AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void load_input_states(const int i, const int j, const int k, const int idir, - Array4 const& qleft_arr, - Array4 const& qright_arr, - Array4 const& qaux_arr, + amrex::Array4 const& qleft_arr, + amrex::Array4 const& qright_arr, + amrex::Array4 const& qaux_arr, RiemannState& ql, RiemannState& qr, RiemannAux& raux) { - const Real small = 1.e-8_rt; + const amrex::Real small = 1.e-8_rt; // fill the states directly from inputs @@ -245,134 +247,5 @@ load_input_states(const int i, const int j, const int k, const int idir, } -AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -void -wsqge(const Real p, const Real v, - const Real gam, const Real gdot, Real& gstar, - const Real gmin, const Real gmax, const Real csq, - const Real pstar, Real& wsq) { - - // compute the lagrangian wave speeds -- this is the approximate - // version for the Colella & Glaz algorithm - - - // First predict a value of game across the shock - - // CG Eq. 31 - gstar = (pstar-p)*gdot/(pstar+p) + gam; - gstar = amrex::max(gmin, amrex::min(gmax, gstar)); - - // Now use that predicted value of game with the R-H jump conditions - // to compute the wave speed. - - // this is CG Eq. 34 - Real alpha = pstar - (gstar - 1.0_rt)*p/(gam - 1.0_rt); - if (alpha == 0.0_rt) { - alpha = riemann_constants::smlp1*(pstar + p); - } - - Real beta = pstar + 0.5_rt*(gstar - 1.0_rt)*(pstar+p); - - wsq = (pstar-p)*beta/(v*alpha); - - if (std::abs(pstar - p) < riemann_constants::smlp1*(pstar + p)) { - wsq = csq; - } - wsq = amrex::max(wsq, (0.5_rt * (gam - 1.0_rt)/gam)*csq); -} - - -AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -void -pstar_bisection(Real& pstar_lo, Real& pstar_hi, - const Real ul, const Real pl, const Real taul, - const Real gamel, const Real clsql, - const Real ur, const Real pr, const Real taur, - const Real gamer, const Real clsqr, - const Real gdot, const Real gmin, const Real gmax, - const int lcg_maxiter, const Real lcg_tol, - Real& pstar, Real& gamstar, - bool& converged, GpuArray& pstar_hist_extra) { - - // we want to zero - // f(p*) = u*_l(p*) - u*_r(p*) - // we'll do bisection - // - // this version is for the approximate Colella & Glaz - // version - - - // lo bounds - Real wlsq = 0.0; - wsqge(pl, taul, gamel, gdot, - gamstar, gmin, gmax, clsql, pstar_lo, wlsq); - - Real wrsq = 0.0; - wsqge(pr, taur, gamer, gdot, - gamstar, gmin, gmax, clsqr, pstar_lo, wrsq); - - Real wl = 1.0_rt / std::sqrt(wlsq); - Real wr = 1.0_rt / std::sqrt(wrsq); - - Real ustar_l = ul - (pstar_lo - pstar)*wl; - Real ustar_r = ur + (pstar_lo - pstar)*wr; - - Real f_lo = ustar_l - ustar_r; - - // hi bounds - wsqge(pl, taul, gamel, gdot, - gamstar, gmin, gmax, clsql, pstar_hi, wlsq); - - wsqge(pr, taur, gamer, gdot, - gamstar, gmin, gmax, clsqr, pstar_hi, wrsq); - - // wl = 1.0_rt / std::sqrt(wlsq); - // wr = 1.0_rt / std::sqrt(wrsq); - - //ustar_l = ul - (pstar_hi - pstar)*wl; - //ustar_r = ur + (pstar_hi - pstar)*wr; - - //Real f_hi = ustar_l - ustar_r; - - // bisection - converged = false; - Real pstar_c = 0.0; - - for (int iter = 0; iter < PSTAR_BISECT_FACTOR*lcg_maxiter; iter++) { - - pstar_c = 0.5_rt * (pstar_lo + pstar_hi); - pstar_hist_extra[iter] = pstar_c; - - wsqge(pl, taul, gamel, gdot, - gamstar, gmin, gmax, clsql, pstar_c, wlsq); - - wsqge(pr, taur, gamer, gdot, - gamstar, gmin, gmax, clsqr, pstar_c, wrsq); - - wl = 1.0_rt / std::sqrt(wlsq); - wr = 1.0_rt / std::sqrt(wrsq); - - ustar_l = ul - (pstar_c - pl)*wl; - ustar_r = ur - (pstar_c - pr)*wr; - - Real f_c = ustar_l - ustar_r; - - if ( 0.5_rt * std::abs(pstar_lo - pstar_hi) < lcg_tol * pstar_c ) { - converged = true; - break; - } - - if (f_lo * f_c < 0.0_rt) { - // root is in the left half - pstar_hi = pstar_c; - //f_hi = f_c; - } else { - pstar_lo = pstar_c; - f_lo = f_c; - } - } - - pstar = pstar_c; -} #endif diff --git a/Source/hydro/riemann_2shock_solvers.H b/Source/hydro/riemann_2shock_solvers.H new file mode 100644 index 0000000000..a1064545c9 --- /dev/null +++ b/Source/hydro/riemann_2shock_solvers.H @@ -0,0 +1,751 @@ +#ifndef RIEMANN_2SHOCK_SOLVERS_H +#define RIEMANN_2SHOCK_SOLVERS_H + +#include + +using namespace amrex::literals; + +#include +#ifdef RADIATION +#include +#include +#endif + + +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE +void +wsqge(const amrex::Real p, const amrex::Real v, + const amrex::Real gam, const amrex::Real gdot, amrex::Real& gstar, + const amrex::Real gmin, const amrex::Real gmax, const amrex::Real csq, + const amrex::Real pstar, amrex::Real& wsq) { + + // compute the lagrangian wave speeds -- this is the approximate + // version for the Colella & Glaz algorithm + + + // First predict a value of game across the shock + + // CG Eq. 31 + gstar = (pstar-p)*gdot/(pstar+p) + gam; + gstar = amrex::max(gmin, amrex::min(gmax, gstar)); + + // Now use that predicted value of game with the R-H jump conditions + // to compute the wave speed. + + // this is CG Eq. 34 + amrex::Real alpha = pstar - (gstar - 1.0_rt)*p/(gam - 1.0_rt); + if (alpha == 0.0_rt) { + alpha = riemann_constants::smlp1*(pstar + p); + } + + amrex::Real beta = pstar + 0.5_rt*(gstar - 1.0_rt)*(pstar+p); + + wsq = (pstar-p)*beta/(v*alpha); + + if (std::abs(pstar - p) < riemann_constants::smlp1*(pstar + p)) { + wsq = csq; + } + wsq = amrex::max(wsq, (0.5_rt * (gam - 1.0_rt)/gam)*csq); +} + + +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE +void +pstar_bisection(amrex::Real& pstar_lo, amrex::Real& pstar_hi, + const amrex::Real ul, const amrex::Real pl, const amrex::Real taul, + const amrex::Real gamel, const amrex::Real clsql, + const amrex::Real ur, const amrex::Real pr, const amrex::Real taur, + const amrex::Real gamer, const amrex::Real clsqr, + const amrex::Real gdot, const amrex::Real gmin, const amrex::Real gmax, + const int lcg_maxiter, const amrex::Real lcg_tol, + amrex::Real& pstar, amrex::Real& gamstar, + bool& converged, amrex::GpuArray& pstar_hist_extra) { + + // we want to zero + // f(p*) = u*_l(p*) - u*_r(p*) + // we'll do bisection + // + // this version is for the approximate Colella & Glaz + // version + + + // lo bounds + amrex::Real wlsq = 0.0; + wsqge(pl, taul, gamel, gdot, + gamstar, gmin, gmax, clsql, pstar_lo, wlsq); + + amrex::Real wrsq = 0.0; + wsqge(pr, taur, gamer, gdot, + gamstar, gmin, gmax, clsqr, pstar_lo, wrsq); + + amrex::Real wl = 1.0_rt / std::sqrt(wlsq); + amrex::Real wr = 1.0_rt / std::sqrt(wrsq); + + amrex::Real ustar_l = ul - (pstar_lo - pstar)*wl; + amrex::Real ustar_r = ur + (pstar_lo - pstar)*wr; + + amrex::Real f_lo = ustar_l - ustar_r; + + // hi bounds + wsqge(pl, taul, gamel, gdot, + gamstar, gmin, gmax, clsql, pstar_hi, wlsq); + + wsqge(pr, taur, gamer, gdot, + gamstar, gmin, gmax, clsqr, pstar_hi, wrsq); + + // wl = 1.0_rt / std::sqrt(wlsq); + // wr = 1.0_rt / std::sqrt(wrsq); + + //ustar_l = ul - (pstar_hi - pstar)*wl; + //ustar_r = ur + (pstar_hi - pstar)*wr; + + //amrex::Real f_hi = ustar_l - ustar_r; + + // bisection + converged = false; + amrex::Real pstar_c = 0.0; + + for (int iter = 0; iter < PSTAR_BISECT_FACTOR*lcg_maxiter; iter++) { + + pstar_c = 0.5_rt * (pstar_lo + pstar_hi); + pstar_hist_extra[iter] = pstar_c; + + wsqge(pl, taul, gamel, gdot, + gamstar, gmin, gmax, clsql, pstar_c, wlsq); + + wsqge(pr, taur, gamer, gdot, + gamstar, gmin, gmax, clsqr, pstar_c, wrsq); + + wl = 1.0_rt / std::sqrt(wlsq); + wr = 1.0_rt / std::sqrt(wrsq); + + ustar_l = ul - (pstar_c - pl)*wl; + ustar_r = ur - (pstar_c - pr)*wr; + + amrex::Real f_c = ustar_l - ustar_r; + + if ( 0.5_rt * std::abs(pstar_lo - pstar_hi) < lcg_tol * pstar_c ) { + converged = true; + break; + } + + if (f_lo * f_c < 0.0_rt) { + // root is in the left half + pstar_hi = pstar_c; + //f_hi = f_c; + } else { + pstar_lo = pstar_c; + f_lo = f_c; + } + } + + pstar = pstar_c; +} + + +/// +/// The Colella-Glaz Riemann solver for pure hydrodynamics. This is a +/// two shock approximate state Riemann solver. +/// +/// @param bx the box to operate over +/// @param ql the left interface state +/// @param qr the right interface state +/// @param qaux_arr the auxiliary state +/// @param qint the full Godunov state on the interface +/// @param idir coordinate direction for the solve (0 = x, 1 = y, 2 = z) +/// +AMREX_GPU_HOST_DEVICE AMREX_INLINE +void +riemanncg(const RiemannState& ql, const RiemannState& qr, const RiemannAux& raux, + RiemannState& qint) { + + // this implements the approximate Riemann solver of Colella & Glaz + // (1985) + // + + constexpr amrex::Real weakwv = 1.e-3_rt; + +#ifndef AMREX_USE_GPU + amrex::GpuArray pstar_hist; +#endif + + + // common quantities + amrex::Real taul = 1.0_rt / ql.rho; + amrex::Real taur = 1.0_rt / qr.rho; + + // lagrangian sound speeds + amrex::Real clsql = ql.gamc * ql.p * ql.rho; + amrex::Real clsqr = qr.gamc * qr.p * qr.rho; + + // Note: in the original Colella & Glaz paper, they predicted + // gamma_e to the interfaces using a special (non-hyperbolic) + // evolution equation. In Castro, we instead bring (rho e) + // to the edges, so we construct the necessary gamma_e here from + // what we have on the interfaces. + amrex::Real gamel = ql.p / ql.rhoe + 1.0_rt; + amrex::Real gamer = qr.p / qr.rhoe + 1.0_rt; + + // these should consider a wider average of the cell-centered + // gammas + amrex::Real gmin = amrex::min(amrex::min(gamel, gamer), 1.0_rt); + amrex::Real gmax = amrex::max(amrex::max(gamel, gamer), 2.0_rt); + + amrex::Real game_bar = 0.5_rt*(gamel + gamer); + amrex::Real gamc_bar = 0.5_rt*(ql.gamc + qr.gamc); + + amrex::Real gdot = 2.0_rt*(1.0_rt - game_bar/gamc_bar)*(game_bar - 1.0_rt); + + amrex::Real wsmall = small_dens * raux.csmall; + amrex::Real wl = amrex::max(wsmall, std::sqrt(std::abs(clsql))); + amrex::Real wr = amrex::max(wsmall, std::sqrt(std::abs(clsqr))); + + // make an initial guess for pstar -- this is a two-shock + // approximation + //pstar = ((wr*pl + wl*pr) + wl*wr*(ul - ur))/(wl + wr) + amrex::Real pstar = ql.p + ( (qr.p - ql.p) - wr*(qr.un - ql.un) ) * wl / (wl + wr); + pstar = amrex::max(pstar, small_pres); + + // get the shock speeds -- this computes W_s from CG Eq. 34 + amrex::Real gamstar = 0.0; + + amrex::Real wlsq = 0.0; + wsqge(ql.p, taul, gamel, gdot, gamstar, + gmin, gmax, clsql, pstar, wlsq); + + amrex::Real wrsq = 0.0; + wsqge(qr.p, taur, gamer, gdot, gamstar, + gmin, gmax, clsqr, pstar, wrsq); + + amrex::Real pstar_old = pstar; + + wl = std::sqrt(wlsq); + wr = std::sqrt(wrsq); + + // R-H jump conditions give ustar across each wave -- these + // should be equal when we are done iterating. Our notation + // here is a little funny, comparing to CG, ustar_l = u*_L and + // ustar_r = u*_R. + amrex::Real ustar_l = ql.un - (pstar - ql.p) / wl; + amrex::Real ustar_r = qr.un + (pstar - qr.p) / wr; + + // revise our pstar guess + // pstar = ((wr*pl + wl*pr) + wl*wr*(ul - ur))/(wl + wr) + pstar = ql.p + ( (qr.p - ql.p) - wr * (qr.un - ql.un) ) * wl / (wl + wr); + pstar = amrex::max(pstar, small_pres); + + // secant iteration + bool converged = false; + + int iter = 0; + while ((iter < cg_maxiter && !converged) || iter < 2) { + + wsqge(ql.p, taul, gamel, gdot, gamstar, + gmin, gmax, clsql, pstar, wlsq); + + wsqge(qr.p, taur, gamer, gdot, gamstar, + gmin, gmax, clsqr, pstar, wrsq); + + + // NOTE: these are really the inverses of the wave speeds! + wl = 1.0_rt / std::sqrt(wlsq); + wr = 1.0_rt / std::sqrt(wrsq); + + amrex::Real ustar_r_old = ustar_r; + amrex::Real ustar_l_old = ustar_l; + + ustar_r = qr.un - (qr.p - pstar) * wr; + ustar_l = ql.un + (ql.p - pstar) * wl; + + amrex::Real dpditer = std::abs(pstar_old - pstar); + + // Here we are going to do the Secant iteration version in + // CG. Note that what we call zp and zm here are not + // actually the Z_p = |dp*/du*_p| defined in CG, by rather + // simply |du*_p| (or something that looks like dp/Z!). + amrex::Real zp = std::abs(ustar_l - ustar_l_old); + if (zp - weakwv * raux.cavg <= 0.0_rt) { + zp = dpditer * wl; + } + + amrex::Real zm = std::abs(ustar_r - ustar_r_old); + if (zm - weakwv * raux.cavg <= 0.0_rt) { + zm = dpditer * wr; + } + + // the new pstar is found via CG Eq. 18 + amrex::Real denom = dpditer / amrex::max(zp + zm, riemann_constants::small * raux.cavg); + pstar_old = pstar; + pstar = pstar - denom*(ustar_r - ustar_l); + pstar = amrex::max(pstar, small_pres); + + amrex::Real err = std::abs(pstar - pstar_old); + if (err < cg_tol*pstar) { + converged = true; + } + +#ifndef AMREX_USE_GPU + pstar_hist[iter] = pstar; +#endif + + iter++; + } + + // If we failed to converge using the secant iteration, we + // can either stop here; or, revert to the original + // two-shock estimate for pstar; or do a bisection root + // find using the bounds established by the most recent + // iterations. + + if (!converged) { + + if (cg_blend == 0) { + +#ifndef AMREX_USE_GPU + std::cout << "pstar history: " << std::endl; + for (int iter_l=0; iter_l < cg_maxiter; iter_l++) { + std::cout << iter_l << " " << pstar_hist[iter_l] << std::endl; + } + + std::cout << std::endl; + std::cout << "left state: " << std::endl << ql << std::endl; + std::cout << "right state: " << std::endl << qr << std::endl; + std::cout << "aux information: " << std::endl << raux << std::endl; + + amrex::Error("ERROR: non-convergence in the Riemann solver"); +#endif + + } else if (cg_blend == 1) { + + pstar = ql.p + ( (qr.p - ql.p) - wr * (qr.un - ql.un) ) * wl / (wl + wr); + + } else if (cg_blend == 2) { + + // we don't store the history if we are in CUDA, so + // we can't do this +#ifndef AMREX_USE_GPU + // first try to find a reasonable bounds + amrex::Real pstarl = 1.e200; + amrex::Real pstaru = -1.e200; + for (int n = cg_maxiter-6; n < cg_maxiter; n++) { + pstarl = amrex::min(pstarl, pstar_hist[n]); + pstaru = amrex::max(pstaru, pstar_hist[n]); + } + + pstarl = amrex::max(pstarl, small_pres); + pstaru = amrex::max(pstaru, small_pres); + + amrex::GpuArray pstar_hist_extra; + + pstar_bisection(pstarl, pstaru, + ql.un, ql.p, taul, gamel, clsql, + qr.un, qr.p, taur, gamer, clsqr, + gdot, gmin, gmax, + cg_maxiter, cg_tol, + pstar, gamstar, converged, pstar_hist_extra); + + if (!converged) { + + std::cout << "pstar history: " << std::endl; + for (int iter_l = 0; iter_l < cg_maxiter; iter_l++) { + std::cout << iter_l << " " << pstar_hist[iter_l] << std::endl; + } + std::cout << "pstar extra history: " << std::endl; + for (int iter_l = 0; iter_l < PSTAR_BISECT_FACTOR*cg_maxiter; iter_l++) { + std::cout << iter_l << " " << pstar_hist_extra[iter_l] << std::endl; + } + + std::cout << std::endl; + std::cout << "left state: " << std::endl << ql << std::endl; + std::cout << "right state: " << std::endl << qr << std::endl; + std::cout << "aux information: " << std::endl << raux << std::endl; + + amrex::Error("ERROR: non-convergence in the Riemann solver"); + } + +#endif + } else { + +#ifndef AMREX_USE_GPU + amrex::Error("ERROR: unrecognized cg_blend option."); +#endif + } + + } + + // we converged! construct the single ustar for the region + // between the left and right waves, using the updated wave speeds + ustar_r = qr.un - (qr.p - pstar) * wr; // careful -- here wl, wr are 1/W + ustar_l = ql.un + (ql.p - pstar) * wl; + + amrex::Real ustar = 0.5_rt * (ustar_l + ustar_r); + + // for symmetry preservation, if ustar is really small, then we + // set it to zero + if (std::abs(ustar) < riemann_constants::smallu * 0.5_rt * (std::abs(ql.un) + std::abs(qr.un))) { + ustar = 0.0_rt; + } + + // sample the solution -- here we look first at the direction + // that the contact is moving. This tells us if we need to + // worry about the L/L* states or the R*/R states. + amrex::Real ro; + amrex::Real uo; + amrex::Real po; + amrex::Real tauo; + amrex::Real gamco; + amrex::Real gameo; + + if (ustar > 0.0_rt) { + //ro = ql.rho; + uo = ql.un; + po = ql.p; + tauo = taul; + gamco = ql.gamc; + gameo = gamel; + + } else if (ustar < 0.0_rt) { + //ro = qr.rho; + uo = qr.un; + po = qr.p; + tauo = taur; + gamco = qr.gamc; + gameo = gamer; + + } else { + //ro = 0.5_rt * (ql.rho + qr.rho); + uo = 0.5_rt * (ql.un + qr.un); + po = 0.5_rt * (ql.p + qr.p); + tauo = 0.5_rt * (taul + taur); + gamco = 0.5_rt * (ql.gamc + qr.gamc); + gameo = 0.5_rt * (gamel + gamer); + } + + // use tau = 1/rho as the independent variable here + ro = amrex::max(small_dens, 1.0_rt/tauo); + tauo = 1.0_rt/ro; + + amrex::Real co = std::sqrt(std::abs(gamco*po*tauo)); + co = amrex::max(raux.csmall, co); + amrex::Real clsq = std::pow(co*ro, 2); + + // now that we know which state (left or right) we need to worry + // about, get the value of gamstar and wosq across the wave we + // are dealing with. + amrex::Real wosq = 0.0; + wsqge(po, tauo, gameo, gdot, gamstar, + gmin, gmax, clsq, pstar, wosq); + + amrex::Real sgnm = std::copysign(1.0_rt, ustar); + + amrex::Real wo = std::sqrt(wosq); + amrex::Real dpjmp = pstar - po; + + // is this max really necessary? + //rstar=max(ONE-ro*dpjmp/wosq, (gameo-ONE)/(gameo+ONE)) + amrex::Real rstar = 1.0_rt - ro*dpjmp/wosq; + rstar = ro/rstar; + rstar = amrex::max(small_dens, rstar); + + amrex::Real cstar = std::sqrt(std::abs(gamco * pstar / rstar)); + cstar = amrex::max(cstar, raux.csmall); + + amrex::Real spout = co - sgnm*uo; + amrex::Real spin = cstar - sgnm*ustar; + + //ushock = 0.5_rt*(spin + spout) + amrex::Real ushock = wo*tauo - sgnm*uo; + + if (pstar - po >= 0.0_rt) { + spin = ushock; + spout = ushock; + } + + amrex::Real frac = 0.5_rt*(1.0_rt + (spin + spout)/amrex::max(amrex::max(spout-spin, spin+spout), + riemann_constants::small * raux.cavg)); + + // the transverse velocity states only depend on the + // direction that the contact moves + if (ustar > 0.0_rt) { + qint.ut = ql.ut; + qint.utt = ql.utt; + } else if (ustar < 0.0_rt) { + qint.ut = qr.ut; + qint.utt = qr.utt; + } else { + qint.ut = 0.5_rt * (ql.ut + qr.ut); + qint.utt = 0.5_rt * (ql.utt + qr.utt); + } + + // linearly interpolate between the star and normal state -- this covers the + // case where we are inside the rarefaction fan. + qint.rho = frac*rstar + (1.0_rt - frac)*ro; + qint.un = frac*ustar + (1.0_rt - frac)*uo; + qint.p = frac*pstar + (1.0_rt - frac)*po; + amrex::Real game_int = frac*gamstar + (1.0_rt-frac)*gameo; + + // now handle the cases where instead we are fully in the + // star or fully in the original (l/r) state + if (spout < 0.0_rt) { + qint.rho = ro; + qint.un = uo; + qint.p = po; + game_int = gameo; + } + + if (spin >= 0.0_rt) { + qint.rho = rstar; + qint.un = ustar; + qint.p = pstar; + game_int = gamstar; + } + + qint.p = amrex::max(qint.p, small_pres); + + qint.un = qint.un * raux.bnd_fac; + + // compute the total energy from the internal, p/(gamma - 1), and the kinetic + qint.rhoe = qint.p / (game_int - 1.0_rt); + + // we'll do the passive scalars separately + +} + + +/// +/// The Colella-Glaz-Ferguson Riemann solver for hydrodynamics and +/// radiation hydrodynamics. This is a two shock approximate state +/// Riemann solver. +/// +/// @param bx the box to operate over +/// @param ql the left interface state +/// @param qr the right interface state +/// @param qaux_arr the auxiliary state +/// @param qint the full Godunov state on the interface +/// @param lambda_int the radiation flux limiter on the interface +/// +AMREX_GPU_HOST_DEVICE AMREX_INLINE +void +riemannus(const RiemannState& ql, const RiemannState& qr, const RiemannAux& raux, + RiemannState& qint) { + + // Colella, Glaz, and Ferguson solver + // + // this is a 2-shock solver that uses a very simple approximation for the + // star state, and carries an auxiliary jump condition for (rho e) to + // deal with a real gas + + + // estimate the star state: pstar, ustar + + amrex::Real wsmall = small_dens * raux.csmall; + + // this is Castro I: Eq. 33 + + amrex::Real wl = amrex::max(wsmall, std::sqrt(std::abs(ql.gamc * ql.p * ql.rho))); + amrex::Real wr = amrex::max(wsmall, std::sqrt(std::abs(qr.gamc * qr.p * qr.rho))); + + amrex::Real wwinv = 1.0_rt/(wl + wr); + amrex::Real pstar = ((wr * ql.p + wl * qr.p) + wl * wr * (ql.un - qr.un)) * wwinv; + amrex::Real ustar = ((wl * ql.un + wr * qr.un) + (ql.p - qr.p)) * wwinv; + + pstar = amrex::max(pstar, small_pres); + + // for symmetry preservation, if ustar is really small, then we + // set it to zero + + if (std::abs(ustar) < riemann_constants::smallu * 0.5_rt * (std::abs(ql.un) + std::abs(qr.un))) { + ustar = 0.0_rt; + } + + // look at the contact to determine which region we are in + + // this just determines which of the left or right states is still + // in play. We still need to look at the other wave to determine + // if the star state or this state is on the interface. + amrex::Real sgnm = std::copysign(1.0_rt, ustar); + if (ustar == 0.0_rt) { + sgnm = 0.0_rt; + } + + amrex::Real fp = 0.5_rt*(1.0_rt + sgnm); + amrex::Real fm = 0.5_rt*(1.0_rt - sgnm); + + amrex::Real ro = fp * ql.rho + fm * qr.rho; + amrex::Real uo = fp * ql.un + fm * qr.un; + amrex::Real po = fp * ql.p + fm * qr.p; + amrex::Real reo = fp * ql.rhoe + fm * qr.rhoe; + amrex::Real gamco = fp * ql.gamc + fm * qr.gamc; +#ifdef RADIATION + for (int g = 0; g < NGROUPS; g++) { + qint.lam[g] = fp * ql.lam[g] + fm * qr.lam[g]; + } + + if (ustar == 0) { + // harmonic average + for (int g = 0; g < NGROUPS; g++) { + qint.lam[g] = 2.0_rt * (ql.lam[g] * qr.lam[g]) / (ql.lam[g] + qr.lam[g] + 1.e-50_rt); + } + } + + amrex::Real po_g = fp * ql.p_g + fm * qr.p_g; + amrex::Real reo_r[NGROUPS]; + amrex::Real po_r[NGROUPS]; + for (int g = 0; g < NGROUPS; g++) { + reo_r[g] = fp * ql.er[g] + fm * qr.er[g]; + po_r[g] = qint.lam[g] * reo_r[g]; + } + amrex::Real reo_g = fp * ql.rhoe_g + fm * qr.rhoe_g; + amrex::Real gamco_g = fp * ql.gamcg + fm * qr.gamcg; +#endif + + ro = amrex::max(small_dens, ro); + + amrex::Real roinv = 1.0_rt / ro; + + amrex::Real co = std::sqrt(std::abs(gamco * po * roinv)); + co = amrex::max(raux.csmall, co); + amrex::Real co2inv = 1.0_rt / (co*co); + + // we can already deal with the transverse velocities -- they + // only jump across the contact + + qint.ut = fp * ql.ut + fm * qr.ut; + qint.utt = fp * ql.utt + fm * qr.utt; + + // compute the rest of the star state + + amrex::Real drho = (pstar - po)*co2inv; + amrex::Real rstar = ro + drho; + rstar = amrex::max(small_dens, rstar); + +#ifdef RADIATION + amrex::Real estar_g = reo_g + drho*(reo_g + po_g)*roinv; + + amrex::Real co_g = std::sqrt(std::abs(gamco_g*po_g*roinv)); + co_g = amrex::max(raux.csmall, co_g); + + amrex::Real pstar_g = po_g + drho*co_g*co_g; + pstar_g = amrex::max(pstar_g, small_pres); + + amrex::Real estar_r[NGROUPS]; + for (int g = 0; g < NGROUPS; g++) { + estar_r[g] = reo_r[g] + drho*(reo_r[g] + po_r[g])*roinv; + } +#else + amrex::Real entho = (reo + po)*roinv*co2inv; + amrex::Real estar = reo + (pstar - po)*entho; +#endif + + amrex::Real cstar = std::sqrt(std::abs(gamco*pstar/rstar)); + cstar = amrex::max(cstar, raux.csmall); + + // finish sampling the solution + + // look at the remaining wave to determine if the star state or the + // 'o' state above is on the interface + + // the values of u +/- c on either side of the non-contact wave + amrex::Real spout = co - sgnm*uo; + amrex::Real spin = cstar - sgnm*ustar; + + // a simple estimate of the shock speed + amrex::Real ushock = 0.5_rt*(spin + spout); + + if (pstar-po > 0.0_rt) { + spin = ushock; + spout = ushock; + } + + amrex::Real scr = spout - spin; + if (spout-spin == 0.0_rt) { + scr = riemann_constants::small * raux.cavg; + } + + // interpolate for the case that we are in a rarefaction + amrex::Real frac = (1.0_rt + (spout + spin)/scr)*0.5_rt; + frac = amrex::max(0.0_rt, amrex::min(1.0_rt, frac)); + + qint.rho = frac*rstar + (1.0_rt - frac)*ro; + qint.un = frac*ustar + (1.0_rt - frac)*uo; + +#ifdef RADIATION + amrex::Real pgdnv_t = frac*pstar + (1.0_rt - frac)*po; + amrex::Real pgdnv_g = frac*pstar_g + (1.0_rt - frac)*po_g; + amrex::Real regdnv_g = frac*estar_g + (1.0_rt - frac)*reo_g; + amrex::Real regdnv_r[NGROUPS]; + for (int g = 0; g < NGROUPS; g++) { + regdnv_r[g] = frac*estar_r[g] + (1.0_rt - frac)*reo_r[g]; + } +#else + qint.p = frac*pstar + (1.0_rt - frac)*po; + amrex::Real regdnv = frac*estar + (1.0_rt - frac)*reo; +#endif + + // as it stands now, we set things assuming that the rarefaction + // spans the interface. We overwrite that here depending on the + // wave speeds + + // look at the speeds on either side of the remaining wave + // to determine which region we are in + if (spout < 0.0_rt) { + // the l or r state is on the interface + qint.rho = ro; + qint.un = uo; +#ifdef RADIATION + pgdnv_t = po; + pgdnv_g = po_g; + regdnv_g = reo_g; + for (int g = 0; g < NGROUPS; g++) { + regdnv_r[g] = reo_r[g]; + } +#else + qint.p = po; + regdnv = reo; +#endif + } + + if (spin >= 0.0_rt) { + // the star state is on the interface + qint.rho = rstar; + qint.un = ustar; +#ifdef RADIATION + pgdnv_t = pstar; + pgdnv_g = pstar_g; + regdnv_g = estar_g; + for (int g = 0; g < NGROUPS; g++) { + regdnv_r[g] = estar_r[g]; + } +#else + qint.p = pstar; + regdnv = estar; +#endif + } + +#ifdef RADIATION + for (int g = 0; g < NGROUPS; g++) { + qint.er[g] = amrex::max(regdnv_r[g], 0.0_rt); + } + + qint.p_g = pgdnv_g; + qint.p = pgdnv_t; + qint.rhoe_g = regdnv_g; + + qint.rhoe = regdnv_g; + for (int g = 0; g < NGROUPS; g++) { + qint.rhoe += regdnv_r[g]; + } + +#else + qint.p = amrex::max(qint.p, small_pres); + qint.rhoe = regdnv; +#endif + + // Enforce that fluxes through a symmetry plane or wall are hard zero. + qint.un = qint.un * raux.bnd_fac; + + // we'll do the passive scalars separately + +} + +#endif diff --git a/Source/hydro/riemann_solvers.H b/Source/hydro/riemann_solvers.H index 51a804e380..ded2d5a169 100644 --- a/Source/hydro/riemann_solvers.H +++ b/Source/hydro/riemann_solvers.H @@ -3,6 +3,7 @@ #include #include +#include #ifdef HYBRID_MOMENTUM #include #endif @@ -214,612 +215,6 @@ compute_flux_q(const int i, const int j, const int k, const int idir, } -/// -/// The Colella-Glaz Riemann solver for pure hydrodynamics. This is a -/// two shock approximate state Riemann solver. -/// -/// @param bx the box to operate over -/// @param ql the left interface state -/// @param qr the right interface state -/// @param qaux_arr the auxiliary state -/// @param qint the full Godunov state on the interface -/// @param idir coordinate direction for the solve (0 = x, 1 = y, 2 = z) -/// -AMREX_GPU_HOST_DEVICE AMREX_INLINE -void -riemanncg(const RiemannState& ql, const RiemannState& qr, const RiemannAux& raux, - RiemannState& qint) { - - // this implements the approximate Riemann solver of Colella & Glaz - // (1985) - // - - constexpr Real weakwv = 1.e-3_rt; - -#ifndef AMREX_USE_GPU - GpuArray pstar_hist; -#endif - - - // common quantities - Real taul = 1.0_rt / ql.rho; - Real taur = 1.0_rt / qr.rho; - - // lagrangian sound speeds - Real clsql = ql.gamc * ql.p * ql.rho; - Real clsqr = qr.gamc * qr.p * qr.rho; - - // Note: in the original Colella & Glaz paper, they predicted - // gamma_e to the interfaces using a special (non-hyperbolic) - // evolution equation. In Castro, we instead bring (rho e) - // to the edges, so we construct the necessary gamma_e here from - // what we have on the interfaces. - Real gamel = ql.p / ql.rhoe + 1.0_rt; - Real gamer = qr.p / qr.rhoe + 1.0_rt; - - // these should consider a wider average of the cell-centered - // gammas - Real gmin = amrex::min(amrex::min(gamel, gamer), 1.0_rt); - Real gmax = amrex::max(amrex::max(gamel, gamer), 2.0_rt); - - Real game_bar = 0.5_rt*(gamel + gamer); - Real gamc_bar = 0.5_rt*(ql.gamc + qr.gamc); - - Real gdot = 2.0_rt*(1.0_rt - game_bar/gamc_bar)*(game_bar - 1.0_rt); - - Real wsmall = small_dens * raux.csmall; - Real wl = amrex::max(wsmall, std::sqrt(std::abs(clsql))); - Real wr = amrex::max(wsmall, std::sqrt(std::abs(clsqr))); - - // make an initial guess for pstar -- this is a two-shock - // approximation - //pstar = ((wr*pl + wl*pr) + wl*wr*(ul - ur))/(wl + wr) - Real pstar = ql.p + ( (qr.p - ql.p) - wr*(qr.un - ql.un) ) * wl / (wl + wr); - pstar = amrex::max(pstar, small_pres); - - // get the shock speeds -- this computes W_s from CG Eq. 34 - Real gamstar = 0.0; - - Real wlsq = 0.0; - wsqge(ql.p, taul, gamel, gdot, gamstar, - gmin, gmax, clsql, pstar, wlsq); - - Real wrsq = 0.0; - wsqge(qr.p, taur, gamer, gdot, gamstar, - gmin, gmax, clsqr, pstar, wrsq); - - Real pstar_old = pstar; - - wl = std::sqrt(wlsq); - wr = std::sqrt(wrsq); - - // R-H jump conditions give ustar across each wave -- these - // should be equal when we are done iterating. Our notation - // here is a little funny, comparing to CG, ustar_l = u*_L and - // ustar_r = u*_R. - Real ustar_l = ql.un - (pstar - ql.p) / wl; - Real ustar_r = qr.un + (pstar - qr.p) / wr; - - // revise our pstar guess - // pstar = ((wr*pl + wl*pr) + wl*wr*(ul - ur))/(wl + wr) - pstar = ql.p + ( (qr.p - ql.p) - wr * (qr.un - ql.un) ) * wl / (wl + wr); - pstar = amrex::max(pstar, small_pres); - - // secant iteration - bool converged = false; - - int iter = 0; - while ((iter < cg_maxiter && !converged) || iter < 2) { - - wsqge(ql.p, taul, gamel, gdot, gamstar, - gmin, gmax, clsql, pstar, wlsq); - - wsqge(qr.p, taur, gamer, gdot, gamstar, - gmin, gmax, clsqr, pstar, wrsq); - - - // NOTE: these are really the inverses of the wave speeds! - wl = 1.0_rt / std::sqrt(wlsq); - wr = 1.0_rt / std::sqrt(wrsq); - - Real ustar_r_old = ustar_r; - Real ustar_l_old = ustar_l; - - ustar_r = qr.un - (qr.p - pstar) * wr; - ustar_l = ql.un + (ql.p - pstar) * wl; - - Real dpditer = std::abs(pstar_old - pstar); - - // Here we are going to do the Secant iteration version in - // CG. Note that what we call zp and zm here are not - // actually the Z_p = |dp*/du*_p| defined in CG, by rather - // simply |du*_p| (or something that looks like dp/Z!). - Real zp = std::abs(ustar_l - ustar_l_old); - if (zp - weakwv * raux.cavg <= 0.0_rt) { - zp = dpditer * wl; - } - - Real zm = std::abs(ustar_r - ustar_r_old); - if (zm - weakwv * raux.cavg <= 0.0_rt) { - zm = dpditer * wr; - } - - // the new pstar is found via CG Eq. 18 - Real denom = dpditer / amrex::max(zp + zm, riemann_constants::small * raux.cavg); - pstar_old = pstar; - pstar = pstar - denom*(ustar_r - ustar_l); - pstar = amrex::max(pstar, small_pres); - - Real err = std::abs(pstar - pstar_old); - if (err < cg_tol*pstar) { - converged = true; - } - -#ifndef AMREX_USE_GPU - pstar_hist[iter] = pstar; -#endif - - iter++; - } - - // If we failed to converge using the secant iteration, we - // can either stop here; or, revert to the original - // two-shock estimate for pstar; or do a bisection root - // find using the bounds established by the most recent - // iterations. - - if (!converged) { - - if (cg_blend == 0) { - -#ifndef AMREX_USE_GPU - std::cout << "pstar history: " << std::endl; - for (int iter_l=0; iter_l < cg_maxiter; iter_l++) { - std::cout << iter_l << " " << pstar_hist[iter_l] << std::endl; - } - - std::cout << std::endl; - std::cout << "left state: " << std::endl << ql << std::endl; - std::cout << "right state: " << std::endl << qr << std::endl; - std::cout << "aux information: " << std::endl << raux << std::endl; - - amrex::Error("ERROR: non-convergence in the Riemann solver"); -#endif - - } else if (cg_blend == 1) { - - pstar = ql.p + ( (qr.p - ql.p) - wr * (qr.un - ql.un) ) * wl / (wl + wr); - - } else if (cg_blend == 2) { - - // we don't store the history if we are in CUDA, so - // we can't do this -#ifndef AMREX_USE_GPU - // first try to find a reasonable bounds - Real pstarl = 1.e200; - Real pstaru = -1.e200; - for (int n = cg_maxiter-6; n < cg_maxiter; n++) { - pstarl = amrex::min(pstarl, pstar_hist[n]); - pstaru = amrex::max(pstaru, pstar_hist[n]); - } - - pstarl = amrex::max(pstarl, small_pres); - pstaru = amrex::max(pstaru, small_pres); - - GpuArray pstar_hist_extra; - - pstar_bisection(pstarl, pstaru, - ql.un, ql.p, taul, gamel, clsql, - qr.un, qr.p, taur, gamer, clsqr, - gdot, gmin, gmax, - cg_maxiter, cg_tol, - pstar, gamstar, converged, pstar_hist_extra); - - if (!converged) { - - std::cout << "pstar history: " << std::endl; - for (int iter_l = 0; iter_l < cg_maxiter; iter_l++) { - std::cout << iter_l << " " << pstar_hist[iter_l] << std::endl; - } - std::cout << "pstar extra history: " << std::endl; - for (int iter_l = 0; iter_l < PSTAR_BISECT_FACTOR*cg_maxiter; iter_l++) { - std::cout << iter_l << " " << pstar_hist_extra[iter_l] << std::endl; - } - - std::cout << std::endl; - std::cout << "left state: " << std::endl << ql << std::endl; - std::cout << "right state: " << std::endl << qr << std::endl; - std::cout << "aux information: " << std::endl << raux << std::endl; - - amrex::Error("ERROR: non-convergence in the Riemann solver"); - } - -#endif - } else { - -#ifndef AMREX_USE_GPU - amrex::Error("ERROR: unrecognized cg_blend option."); -#endif - } - - } - - // we converged! construct the single ustar for the region - // between the left and right waves, using the updated wave speeds - ustar_r = qr.un - (qr.p - pstar) * wr; // careful -- here wl, wr are 1/W - ustar_l = ql.un + (ql.p - pstar) * wl; - - Real ustar = 0.5_rt * (ustar_l + ustar_r); - - // for symmetry preservation, if ustar is really small, then we - // set it to zero - if (std::abs(ustar) < riemann_constants::smallu * 0.5_rt * (std::abs(ql.un) + std::abs(qr.un))) { - ustar = 0.0_rt; - } - - // sample the solution -- here we look first at the direction - // that the contact is moving. This tells us if we need to - // worry about the L/L* states or the R*/R states. - Real ro; - Real uo; - Real po; - Real tauo; - Real gamco; - Real gameo; - - if (ustar > 0.0_rt) { - //ro = ql.rho; - uo = ql.un; - po = ql.p; - tauo = taul; - gamco = ql.gamc; - gameo = gamel; - - } else if (ustar < 0.0_rt) { - //ro = qr.rho; - uo = qr.un; - po = qr.p; - tauo = taur; - gamco = qr.gamc; - gameo = gamer; - - } else { - //ro = 0.5_rt * (ql.rho + qr.rho); - uo = 0.5_rt * (ql.un + qr.un); - po = 0.5_rt * (ql.p + qr.p); - tauo = 0.5_rt * (taul + taur); - gamco = 0.5_rt * (ql.gamc + qr.gamc); - gameo = 0.5_rt * (gamel + gamer); - } - - // use tau = 1/rho as the independent variable here - ro = amrex::max(small_dens, 1.0_rt/tauo); - tauo = 1.0_rt/ro; - - Real co = std::sqrt(std::abs(gamco*po*tauo)); - co = amrex::max(raux.csmall, co); - Real clsq = std::pow(co*ro, 2); - - // now that we know which state (left or right) we need to worry - // about, get the value of gamstar and wosq across the wave we - // are dealing with. - Real wosq = 0.0; - wsqge(po, tauo, gameo, gdot, gamstar, - gmin, gmax, clsq, pstar, wosq); - - Real sgnm = std::copysign(1.0_rt, ustar); - - Real wo = std::sqrt(wosq); - Real dpjmp = pstar - po; - - // is this max really necessary? - //rstar=max(ONE-ro*dpjmp/wosq, (gameo-ONE)/(gameo+ONE)) - Real rstar = 1.0_rt - ro*dpjmp/wosq; - rstar = ro/rstar; - rstar = amrex::max(small_dens, rstar); - - Real cstar = std::sqrt(std::abs(gamco * pstar / rstar)); - cstar = amrex::max(cstar, raux.csmall); - - Real spout = co - sgnm*uo; - Real spin = cstar - sgnm*ustar; - - //ushock = 0.5_rt*(spin + spout) - Real ushock = wo*tauo - sgnm*uo; - - if (pstar - po >= 0.0_rt) { - spin = ushock; - spout = ushock; - } - - Real frac = 0.5_rt*(1.0_rt + (spin + spout)/amrex::max(amrex::max(spout-spin, spin+spout), - riemann_constants::small * raux.cavg)); - - // the transverse velocity states only depend on the - // direction that the contact moves - if (ustar > 0.0_rt) { - qint.ut = ql.ut; - qint.utt = ql.utt; - } else if (ustar < 0.0_rt) { - qint.ut = qr.ut; - qint.utt = qr.utt; - } else { - qint.ut = 0.5_rt * (ql.ut + qr.ut); - qint.utt = 0.5_rt * (ql.utt + qr.utt); - } - - // linearly interpolate between the star and normal state -- this covers the - // case where we are inside the rarefaction fan. - qint.rho = frac*rstar + (1.0_rt - frac)*ro; - qint.un = frac*ustar + (1.0_rt - frac)*uo; - qint.p = frac*pstar + (1.0_rt - frac)*po; - Real game_int = frac*gamstar + (1.0_rt-frac)*gameo; - - // now handle the cases where instead we are fully in the - // star or fully in the original (l/r) state - if (spout < 0.0_rt) { - qint.rho = ro; - qint.un = uo; - qint.p = po; - game_int = gameo; - } - - if (spin >= 0.0_rt) { - qint.rho = rstar; - qint.un = ustar; - qint.p = pstar; - game_int = gamstar; - } - - qint.p = amrex::max(qint.p, small_pres); - - qint.un = qint.un * raux.bnd_fac; - - // compute the total energy from the internal, p/(gamma - 1), and the kinetic - qint.rhoe = qint.p / (game_int - 1.0_rt); - - // we'll do the passive scalars separately - -} - - -/// -/// The Colella-Glaz-Ferguson Riemann solver for hydrodynamics and -/// radiation hydrodynamics. This is a two shock approximate state -/// Riemann solver. -/// -/// @param bx the box to operate over -/// @param ql the left interface state -/// @param qr the right interface state -/// @param qaux_arr the auxiliary state -/// @param qint the full Godunov state on the interface -/// @param lambda_int the radiation flux limiter on the interface -/// -AMREX_GPU_HOST_DEVICE AMREX_INLINE -void -riemannus(const RiemannState& ql, const RiemannState& qr, const RiemannAux& raux, - RiemannState& qint) { - - // Colella, Glaz, and Ferguson solver - // - // this is a 2-shock solver that uses a very simple approximation for the - // star state, and carries an auxiliary jump condition for (rho e) to - // deal with a real gas - - - // estimate the star state: pstar, ustar - - Real wsmall = small_dens * raux.csmall; - - // this is Castro I: Eq. 33 - - Real wl = amrex::max(wsmall, std::sqrt(std::abs(ql.gamc * ql.p * ql.rho))); - Real wr = amrex::max(wsmall, std::sqrt(std::abs(qr.gamc * qr.p * qr.rho))); - - Real wwinv = 1.0_rt/(wl + wr); - Real pstar = ((wr * ql.p + wl * qr.p) + wl * wr * (ql.un - qr.un)) * wwinv; - Real ustar = ((wl * ql.un + wr * qr.un) + (ql.p - qr.p)) * wwinv; - - pstar = amrex::max(pstar, small_pres); - - // for symmetry preservation, if ustar is really small, then we - // set it to zero - - if (std::abs(ustar) < riemann_constants::smallu * 0.5_rt * (std::abs(ql.un) + std::abs(qr.un))) { - ustar = 0.0_rt; - } - - // look at the contact to determine which region we are in - - // this just determines which of the left or right states is still - // in play. We still need to look at the other wave to determine - // if the star state or this state is on the interface. - Real sgnm = std::copysign(1.0_rt, ustar); - if (ustar == 0.0_rt) { - sgnm = 0.0_rt; - } - - Real fp = 0.5_rt*(1.0_rt + sgnm); - Real fm = 0.5_rt*(1.0_rt - sgnm); - - Real ro = fp * ql.rho + fm * qr.rho; - Real uo = fp * ql.un + fm * qr.un; - Real po = fp * ql.p + fm * qr.p; - Real reo = fp * ql.rhoe + fm * qr.rhoe; - Real gamco = fp * ql.gamc + fm * qr.gamc; -#ifdef RADIATION - for (int g = 0; g < NGROUPS; g++) { - qint.lam[g] = fp * ql.lam[g] + fm * qr.lam[g]; - } - - if (ustar == 0) { - // harmonic average - for (int g = 0; g < NGROUPS; g++) { - qint.lam[g] = 2.0_rt * (ql.lam[g] * qr.lam[g]) / (ql.lam[g] + qr.lam[g] + 1.e-50_rt); - } - } - - Real po_g = fp * ql.p_g + fm * qr.p_g; - Real reo_r[NGROUPS]; - Real po_r[NGROUPS]; - for (int g = 0; g < NGROUPS; g++) { - reo_r[g] = fp * ql.er[g] + fm * qr.er[g]; - po_r[g] = qint.lam[g] * reo_r[g]; - } - Real reo_g = fp * ql.rhoe_g + fm * qr.rhoe_g; - Real gamco_g = fp * ql.gamcg + fm * qr.gamcg; -#endif - - ro = amrex::max(small_dens, ro); - - Real roinv = 1.0_rt / ro; - - Real co = std::sqrt(std::abs(gamco * po * roinv)); - co = amrex::max(raux.csmall, co); - Real co2inv = 1.0_rt / (co*co); - - // we can already deal with the transverse velocities -- they - // only jump across the contact - - qint.ut = fp * ql.ut + fm * qr.ut; - qint.utt = fp * ql.utt + fm * qr.utt; - - // compute the rest of the star state - - Real drho = (pstar - po)*co2inv; - Real rstar = ro + drho; - rstar = amrex::max(small_dens, rstar); - -#ifdef RADIATION - Real estar_g = reo_g + drho*(reo_g + po_g)*roinv; - - Real co_g = std::sqrt(std::abs(gamco_g*po_g*roinv)); - co_g = amrex::max(raux.csmall, co_g); - - Real pstar_g = po_g + drho*co_g*co_g; - pstar_g = amrex::max(pstar_g, small_pres); - - Real estar_r[NGROUPS]; - for (int g = 0; g < NGROUPS; g++) { - estar_r[g] = reo_r[g] + drho*(reo_r[g] + po_r[g])*roinv; - } -#else - Real entho = (reo + po)*roinv*co2inv; - Real estar = reo + (pstar - po)*entho; -#endif - - Real cstar = std::sqrt(std::abs(gamco*pstar/rstar)); - cstar = amrex::max(cstar, raux.csmall); - - // finish sampling the solution - - // look at the remaining wave to determine if the star state or the - // 'o' state above is on the interface - - // the values of u +/- c on either side of the non-contact wave - Real spout = co - sgnm*uo; - Real spin = cstar - sgnm*ustar; - - // a simple estimate of the shock speed - Real ushock = 0.5_rt*(spin + spout); - - if (pstar-po > 0.0_rt) { - spin = ushock; - spout = ushock; - } - - Real scr = spout - spin; - if (spout-spin == 0.0_rt) { - scr = riemann_constants::small * raux.cavg; - } - - // interpolate for the case that we are in a rarefaction - Real frac = (1.0_rt + (spout + spin)/scr)*0.5_rt; - frac = amrex::max(0.0_rt, amrex::min(1.0_rt, frac)); - - qint.rho = frac*rstar + (1.0_rt - frac)*ro; - qint.un = frac*ustar + (1.0_rt - frac)*uo; - -#ifdef RADIATION - Real pgdnv_t = frac*pstar + (1.0_rt - frac)*po; - Real pgdnv_g = frac*pstar_g + (1.0_rt - frac)*po_g; - Real regdnv_g = frac*estar_g + (1.0_rt - frac)*reo_g; - Real regdnv_r[NGROUPS]; - for (int g = 0; g < NGROUPS; g++) { - regdnv_r[g] = frac*estar_r[g] + (1.0_rt - frac)*reo_r[g]; - } -#else - qint.p = frac*pstar + (1.0_rt - frac)*po; - Real regdnv = frac*estar + (1.0_rt - frac)*reo; -#endif - - // as it stands now, we set things assuming that the rarefaction - // spans the interface. We overwrite that here depending on the - // wave speeds - - // look at the speeds on either side of the remaining wave - // to determine which region we are in - if (spout < 0.0_rt) { - // the l or r state is on the interface - qint.rho = ro; - qint.un = uo; -#ifdef RADIATION - pgdnv_t = po; - pgdnv_g = po_g; - regdnv_g = reo_g; - for (int g = 0; g < NGROUPS; g++) { - regdnv_r[g] = reo_r[g]; - } -#else - qint.p = po; - regdnv = reo; -#endif - } - - if (spin >= 0.0_rt) { - // the star state is on the interface - qint.rho = rstar; - qint.un = ustar; -#ifdef RADIATION - pgdnv_t = pstar; - pgdnv_g = pstar_g; - regdnv_g = estar_g; - for (int g = 0; g < NGROUPS; g++) { - regdnv_r[g] = estar_r[g]; - } -#else - qint.p = pstar; - regdnv = estar; -#endif - } - -#ifdef RADIATION - for (int g = 0; g < NGROUPS; g++) { - qint.er[g] = amrex::max(regdnv_r[g], 0.0_rt); - } - - qint.p_g = pgdnv_g; - qint.p = pgdnv_t; - qint.rhoe_g = regdnv_g; - - qint.rhoe = regdnv_g; - for (int g = 0; g < NGROUPS; g++) { - qint.rhoe += regdnv_r[g]; - } - -#else - qint.p = amrex::max(qint.p, small_pres); - qint.rhoe = regdnv; -#endif - - // Enforce that fluxes through a symmetry plane or wall are hard zero. - qint.un = qint.un * raux.bnd_fac; - - // we'll do the passive scalars separately - -} - - AMREX_GPU_HOST_DEVICE AMREX_INLINE void From d61bf01fdd8e01c993c0fc670ba5747cae9df944 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Tue, 25 Jun 2024 14:51:34 -0400 Subject: [PATCH 3/5] update Make.package --- Source/hydro/Make.package | 1 + 1 file changed, 1 insertion(+) diff --git a/Source/hydro/Make.package b/Source/hydro/Make.package index c867ef538d..fdff393206 100644 --- a/Source/hydro/Make.package +++ b/Source/hydro/Make.package @@ -27,6 +27,7 @@ endif CEXE_headers += ppm.H CEXE_sources += riemann.cpp +CEXE_headers += HLL_solvers.H CEXE_headers += riemann_solvers.H CEXE_sources += riemann_util.cpp CEXE_headers += riemann.H From 2115bdc6b72de0d8ebc75b0a0ec30c3381049db9 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Tue, 25 Jun 2024 19:49:27 -0400 Subject: [PATCH 4/5] add namespace --- Source/hydro/HLL_solvers.H | 952 +++++++++++++++++++------------------ Source/hydro/riemann.cpp | 22 +- 2 files changed, 488 insertions(+), 486 deletions(-) diff --git a/Source/hydro/HLL_solvers.H b/Source/hydro/HLL_solvers.H index 0acf814735..67f7a0302a 100644 --- a/Source/hydro/HLL_solvers.H +++ b/Source/hydro/HLL_solvers.H @@ -3,590 +3,592 @@ #include -AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -void -HLLC_state(const int idir, const Real S_k, const Real S_c, - const Real* qstate, Real* U) { - - Real u_k = 0.0; - if (idir == 0) { - u_k = qstate[QU]; - } else if (idir == 1) { - u_k = qstate[QV]; - } else if (idir == 2) { - u_k = qstate[QW]; - } - - Real hllc_factor = qstate[QRHO]*(S_k - u_k)/(S_k - S_c); - U[URHO] = hllc_factor; - - if (idir == 0) { - U[UMX] = hllc_factor*S_c; - U[UMY] = hllc_factor*qstate[QV]; - U[UMZ] = hllc_factor*qstate[QW]; - - } else if (idir == 1) { - U[UMX] = hllc_factor*qstate[QU]; - U[UMY] = hllc_factor*S_c; - U[UMZ] = hllc_factor*qstate[QW]; - - } else { - U[UMX] = hllc_factor*qstate[QU]; - U[UMY] = hllc_factor*qstate[QV]; - U[UMZ] = hllc_factor*S_c; - } - - U[UEDEN] = hllc_factor*(qstate[QREINT]/qstate[QRHO] + - 0.5_rt*(qstate[QU]*qstate[QU] + qstate[QV]*qstate[QV] + qstate[QW]*qstate[QW]) + - (S_c - u_k)*(S_c + qstate[QPRES]/(qstate[QRHO]*(S_k - u_k)))); - U[UEINT] = hllc_factor*qstate[QREINT]/qstate[QRHO]; - - U[UTEMP] = 0.0; // we don't evolve T - -#ifdef SHOCK_VAR - U[USHK] = 0.0; -#endif - -#ifdef NSE_NET - U[UMUP] = 0.0; - U[UMUN] = 0.0; -#endif - - for (int ipassive = 0; ipassive < npassive; ipassive++) { - const int n = upassmap(ipassive); - const int nqs = qpassmap(ipassive); - U[n] = hllc_factor*qstate[nqs]; - } -} - - -/// -/// given a conserved state, compute the flux in direction idir -/// -AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -void -compute_flux(const int idir, const Real bnd_fac, const int coord, - const Real* U, const Real p, - Real* F) { - - Real u_flx = U[UMX+idir]/U[URHO]; - - if (bnd_fac == 0) { - u_flx = 0.0; - } - - F[URHO] = U[URHO]*u_flx; - - F[UMX] = U[UMX]*u_flx; - F[UMY] = U[UMY]*u_flx; - F[UMZ] = U[UMZ]*u_flx; - - auto mom_check = mom_flux_has_p(idir, idir, coord); - - if (mom_check) { - // we do not include the pressure term in any non-Cartesian - // coordinate directions - F[UMX+idir] = F[UMX+idir] + p; - } - - F[UEINT] = U[UEINT]*u_flx; - F[UEDEN] = (U[UEDEN] + p)*u_flx; +namespace HLL { + + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE + void + HLLC_state(const int idir, const Real S_k, const Real S_c, + const Real* qstate, Real* U) { + + Real u_k = 0.0; + if (idir == 0) { + u_k = qstate[QU]; + } else if (idir == 1) { + u_k = qstate[QV]; + } else if (idir == 2) { + u_k = qstate[QW]; + } + + Real hllc_factor = qstate[QRHO]*(S_k - u_k)/(S_k - S_c); + U[URHO] = hllc_factor; + + if (idir == 0) { + U[UMX] = hllc_factor*S_c; + U[UMY] = hllc_factor*qstate[QV]; + U[UMZ] = hllc_factor*qstate[QW]; + + } else if (idir == 1) { + U[UMX] = hllc_factor*qstate[QU]; + U[UMY] = hllc_factor*S_c; + U[UMZ] = hllc_factor*qstate[QW]; + + } else { + U[UMX] = hllc_factor*qstate[QU]; + U[UMY] = hllc_factor*qstate[QV]; + U[UMZ] = hllc_factor*S_c; + } + + U[UEDEN] = hllc_factor*(qstate[QREINT]/qstate[QRHO] + + 0.5_rt*(qstate[QU]*qstate[QU] + qstate[QV]*qstate[QV] + qstate[QW]*qstate[QW]) + + (S_c - u_k)*(S_c + qstate[QPRES]/(qstate[QRHO]*(S_k - u_k)))); + U[UEINT] = hllc_factor*qstate[QREINT]/qstate[QRHO]; + + U[UTEMP] = 0.0; // we don't evolve T + + #ifdef SHOCK_VAR + U[USHK] = 0.0; + #endif + + #ifdef NSE_NET + U[UMUP] = 0.0; + U[UMUN] = 0.0; + #endif + + for (int ipassive = 0; ipassive < npassive; ipassive++) { + const int n = upassmap(ipassive); + const int nqs = qpassmap(ipassive); + U[n] = hllc_factor*qstate[nqs]; + } + } - F[UTEMP] = 0.0; -#ifdef SHOCK_VAR - F[USHK] = 0.0; -#endif + /// + /// given a conserved state, compute the flux in direction idir + /// + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE + void + compute_flux(const int idir, const Real bnd_fac, const int coord, + const Real* U, const Real p, + Real* F) { -#ifdef NSE_NET - F[UMUP] = 0.0; - F[UMUN] = 0.0; -#endif + Real u_flx = U[UMX+idir]/U[URHO]; - for (int ipassive=0; ipassive < npassive; ipassive++) { - const int n = upassmap(ipassive); - F[n] = U[n]*u_flx; - } -} + if (bnd_fac == 0) { + u_flx = 0.0; + } -/// -/// convert the primitive variable state to the conserved state -/// -AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -void -cons_state(const Real* qstate, Real* U) { + F[URHO] = U[URHO]*u_flx; - U[URHO] = qstate[QRHO]; + F[UMX] = U[UMX]*u_flx; + F[UMY] = U[UMY]*u_flx; + F[UMZ] = U[UMZ]*u_flx; - // since we advect all 3 velocity components regardless of dimension, this - // will be general - U[UMX] = qstate[QRHO]*qstate[QU]; - U[UMY] = qstate[QRHO]*qstate[QV]; - U[UMZ] = qstate[QRHO]*qstate[QW]; + auto mom_check = mom_flux_has_p(idir, idir, coord); - U[UEDEN] = qstate[QREINT] + 0.5_rt*qstate[QRHO]*(qstate[QU]*qstate[QU] + qstate[QV]*qstate[QV] + qstate[QW]*qstate[QW]); - U[UEINT] = qstate[QREINT]; + if (mom_check) { + // we do not include the pressure term in any non-Cartesian + // coordinate directions + F[UMX+idir] = F[UMX+idir] + p; + } - // we don't care about T here, but initialize it to make NaN - // checking happy - U[UTEMP] = 0.0; + F[UEINT] = U[UEINT]*u_flx; + F[UEDEN] = (U[UEDEN] + p)*u_flx; -#ifdef SHOCK_VAR - U[USHK] = 0.0; -#endif - -#ifdef NSE_NET - U[UMUP] = 0.0; - U[UMUN] = 0.0; -#endif + F[UTEMP] = 0.0; - for (int ipassive = 0; ipassive < npassive; ipassive++) { - const int n = upassmap(ipassive); - const int nqs = qpassmap(ipassive); - U[n] = qstate[QRHO]*qstate[nqs]; - } -} + #ifdef SHOCK_VAR + F[USHK] = 0.0; + #endif + #ifdef NSE_NET + F[UMUP] = 0.0; + F[UMUN] = 0.0; + #endif -/// -/// A simple HLL Riemann solver for pure hydrodynamics. This takes just a -/// single interface's data and returns the HLL flux -/// -/// @param ql the left interface state -/// @param qr the right interface state -/// @param cl sound speed on the left interface -/// @param cr sound speed on the right interface -/// @param idir coordinate direction for the solve (0 = x, 1 = y, 2 = z) -/// @param coord geometry type (0 = Cartesian, 1 = axisymmetric, 2 = spherical) -/// @param f the HLL fluxes -/// -AMREX_GPU_HOST_DEVICE AMREX_INLINE -void -HLL(const Real* ql, const Real* qr, - const Real cl, const Real cr, - const int idir, const int coord, - Real* flux_hll) { + for (int ipassive=0; ipassive < npassive; ipassive++) { + const int n = upassmap(ipassive); + F[n] = U[n]*u_flx; + } + } - // This is the HLLE solver. We should apply it to zone averages - // (not reconstructed states) at an interface in the presence of - // shocks to avoid the odd-even decoupling / carbuncle phenomenon. - // - // See: Einfeldt, B. et al. 1991, JCP, 92, 273 - // Einfeldt, B. 1988, SIAM J NA, 25, 294 + /// + /// convert the primitive variable state to the conserved state + /// + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE + void + cons_state(const Real* qstate, Real* U) { + + U[URHO] = qstate[QRHO]; + + // since we advect all 3 velocity components regardless of dimension, this + // will be general + U[UMX] = qstate[QRHO]*qstate[QU]; + U[UMY] = qstate[QRHO]*qstate[QV]; + U[UMZ] = qstate[QRHO]*qstate[QW]; + + U[UEDEN] = qstate[QREINT] + 0.5_rt*qstate[QRHO]*(qstate[QU]*qstate[QU] + qstate[QV]*qstate[QV] + qstate[QW]*qstate[QW]); + U[UEINT] = qstate[QREINT]; + + // we don't care about T here, but initialize it to make NaN + // checking happy + U[UTEMP] = 0.0; + + #ifdef SHOCK_VAR + U[USHK] = 0.0; + #endif + + #ifdef NSE_NET + U[UMUP] = 0.0; + U[UMUN] = 0.0; + #endif + + for (int ipassive = 0; ipassive < npassive; ipassive++) { + const int n = upassmap(ipassive); + const int nqs = qpassmap(ipassive); + U[n] = qstate[QRHO]*qstate[nqs]; + } + } - constexpr Real small_hll = 1.e-10_rt; + /// + /// A simple HLL Riemann solver for pure hydrodynamics. This takes just a + /// single interface's data and returns the HLL flux + /// + /// @param ql the left interface state + /// @param qr the right interface state + /// @param cl sound speed on the left interface + /// @param cr sound speed on the right interface + /// @param idir coordinate direction for the solve (0 = x, 1 = y, 2 = z) + /// @param coord geometry type (0 = Cartesian, 1 = axisymmetric, 2 = spherical) + /// @param f the HLL fluxes + /// + AMREX_GPU_HOST_DEVICE AMREX_INLINE + void + HLL(const Real* ql, const Real* qr, + const Real cl, const Real cr, + const int idir, const int coord, + Real* flux_hll) { - int ivel, ivelt, iveltt; - int imom, imomt, imomtt; + // This is the HLLE solver. We should apply it to zone averages + // (not reconstructed states) at an interface in the presence of + // shocks to avoid the odd-even decoupling / carbuncle phenomenon. + // + // See: Einfeldt, B. et al. 1991, JCP, 92, 273 + // Einfeldt, B. 1988, SIAM J NA, 25, 294 - if (idir == 0) { - ivel = QU; - ivelt = QV; - iveltt = QW; - imom = UMX; - imomt = UMY; - imomtt = UMZ; + constexpr Real small_hll = 1.e-10_rt; - } else if (idir == 1) { - ivel = QV; - ivelt = QU; - iveltt = QW; + int ivel, ivelt, iveltt; + int imom, imomt, imomtt; - imom = UMY; - imomt = UMX; - imomtt = UMZ; + if (idir == 0) { + ivel = QU; + ivelt = QV; + iveltt = QW; - } else { - ivel = QW; - ivelt = QU; - iveltt = QV; + imom = UMX; + imomt = UMY; + imomtt = UMZ; - imom = UMZ; - imomt = UMX; - imomtt = UMY; - } + } else if (idir == 1) { + ivel = QV; + ivelt = QU; + iveltt = QW; - Real rhol_sqrt = std::sqrt(ql[QRHO]); - Real rhor_sqrt = std::sqrt(qr[QRHO]); + imom = UMY; + imomt = UMX; + imomtt = UMZ; - Real rhod = 1.0_rt/(rhol_sqrt + rhor_sqrt); + } else { + ivel = QW; + ivelt = QU; + iveltt = QV; + imom = UMZ; + imomt = UMX; + imomtt = UMY; + } - // compute the average sound speed. This uses an approximation from - // E88, eq. 5.6, 5.7 that assumes gamma falls between 1 - // and 5/3 - Real cavg = std::sqrt( (rhol_sqrt*cl*cl + rhor_sqrt*cr*cr)*rhod + - 0.5_rt*rhol_sqrt*rhor_sqrt*rhod*rhod*std::pow(qr[ivel] - ql[ivel], 2)); + Real rhol_sqrt = std::sqrt(ql[QRHO]); + Real rhor_sqrt = std::sqrt(qr[QRHO]); + Real rhod = 1.0_rt/(rhol_sqrt + rhor_sqrt); - // Roe eigenvalues (E91, eq. 5.3b) - Real uavg = (rhol_sqrt*ql[ivel] + rhor_sqrt*qr[ivel])*rhod; - Real a1 = uavg - cavg; - Real a4 = uavg + cavg; + // compute the average sound speed. This uses an approximation from + // E88, eq. 5.6, 5.7 that assumes gamma falls between 1 + // and 5/3 + Real cavg = std::sqrt( (rhol_sqrt*cl*cl + rhor_sqrt*cr*cr)*rhod + + 0.5_rt*rhol_sqrt*rhor_sqrt*rhod*rhod*std::pow(qr[ivel] - ql[ivel], 2)); - // signal speeds (E91, eq. 4.5) - Real bl = amrex::min(a1, ql[ivel] - cl); - Real br = amrex::max(a4, qr[ivel] + cr); + // Roe eigenvalues (E91, eq. 5.3b) + Real uavg = (rhol_sqrt*ql[ivel] + rhor_sqrt*qr[ivel])*rhod; - Real bm = amrex::min(0.0_rt, bl); - Real bp = amrex::max(0.0_rt, br); + Real a1 = uavg - cavg; + Real a4 = uavg + cavg; - Real bd = bp - bm; - if (std::abs(bd) < small_hll*amrex::max(std::abs(bm), std::abs(bp))) { - return; - } + // signal speeds (E91, eq. 4.5) + Real bl = amrex::min(a1, ql[ivel] - cl); + Real br = amrex::max(a4, qr[ivel] + cr); - // we'll overwrite the passed in flux with the HLL flux + Real bm = amrex::min(0.0_rt, bl); + Real bp = amrex::max(0.0_rt, br); - bd = 1.0_rt/bd; + Real bd = bp - bm; - // compute the fluxes according to E91, eq. 4.4b -- note that the - // min/max above picks the correct flux if we are not in the star - // region + if (std::abs(bd) < small_hll*amrex::max(std::abs(bm), std::abs(bp))) { + return; + } - // density flux - Real fl_tmp = ql[QRHO]*ql[ivel]; - Real fr_tmp = qr[QRHO]*qr[ivel]; + // we'll overwrite the passed in flux with the HLL flux - flux_hll[URHO] = (bp*fl_tmp - bm*fr_tmp)*bd + bp*bm*bd*(qr[QRHO] - ql[QRHO]); + bd = 1.0_rt/bd; - // normal momentum flux. Note for 1-d and 2-d non cartesian - // r-coordinate, we leave off the pressure term and handle that - // separately in the update, to accommodate different geometries - fl_tmp = ql[QRHO]*ql[ivel]*ql[ivel]; - fr_tmp = qr[QRHO]*qr[ivel]*qr[ivel]; - if (mom_flux_has_p(idir, idir, coord)) { - fl_tmp = fl_tmp + ql[QPRES]; - fr_tmp = fr_tmp + qr[QPRES]; - } + // compute the fluxes according to E91, eq. 4.4b -- note that the + // min/max above picks the correct flux if we are not in the star + // region - flux_hll[imom] = (bp*fl_tmp - bm*fr_tmp)*bd + bp*bm*bd*(qr[QRHO]*qr[ivel] - ql[QRHO]*ql[ivel]); + // density flux + Real fl_tmp = ql[QRHO]*ql[ivel]; + Real fr_tmp = qr[QRHO]*qr[ivel]; - // transverse momentum flux - fl_tmp = ql[QRHO]*ql[ivel]*ql[ivelt]; - fr_tmp = qr[QRHO]*qr[ivel]*qr[ivelt]; + flux_hll[URHO] = (bp*fl_tmp - bm*fr_tmp)*bd + bp*bm*bd*(qr[QRHO] - ql[QRHO]); - flux_hll[imomt] = (bp*fl_tmp - bm*fr_tmp)*bd + bp*bm*bd*(qr[QRHO]*qr[ivelt] - ql[QRHO]*ql[ivelt]); + // normal momentum flux. Note for 1-d and 2-d non cartesian + // r-coordinate, we leave off the pressure term and handle that + // separately in the update, to accommodate different geometries + fl_tmp = ql[QRHO]*ql[ivel]*ql[ivel]; + fr_tmp = qr[QRHO]*qr[ivel]*qr[ivel]; + if (mom_flux_has_p(idir, idir, coord)) { + fl_tmp = fl_tmp + ql[QPRES]; + fr_tmp = fr_tmp + qr[QPRES]; + } + flux_hll[imom] = (bp*fl_tmp - bm*fr_tmp)*bd + bp*bm*bd*(qr[QRHO]*qr[ivel] - ql[QRHO]*ql[ivel]); - fl_tmp = ql[QRHO]*ql[ivel]*ql[iveltt]; - fr_tmp = qr[QRHO]*qr[ivel]*qr[iveltt]; + // transverse momentum flux + fl_tmp = ql[QRHO]*ql[ivel]*ql[ivelt]; + fr_tmp = qr[QRHO]*qr[ivel]*qr[ivelt]; - flux_hll[imomtt] = (bp*fl_tmp - bm*fr_tmp)*bd + bp*bm*bd*(qr[QRHO]*qr[iveltt] - ql[QRHO]*ql[iveltt]); + flux_hll[imomt] = (bp*fl_tmp - bm*fr_tmp)*bd + bp*bm*bd*(qr[QRHO]*qr[ivelt] - ql[QRHO]*ql[ivelt]); - // total energy flux - Real rhoEl = ql[QREINT] + 0.5_rt*ql[QRHO]*(ql[ivel]*ql[ivel] + ql[ivelt]*ql[ivelt] + ql[iveltt]*ql[iveltt]); - fl_tmp = ql[ivel]*(rhoEl + ql[QPRES]); - Real rhoEr = qr[QREINT] + 0.5_rt*qr[QRHO]*(qr[ivel]*qr[ivel] + qr[ivelt]*qr[ivelt] + qr[iveltt]*qr[iveltt]); - fr_tmp = qr[ivel]*(rhoEr + qr[QPRES]); + fl_tmp = ql[QRHO]*ql[ivel]*ql[iveltt]; + fr_tmp = qr[QRHO]*qr[ivel]*qr[iveltt]; - flux_hll[UEDEN] = (bp*fl_tmp - bm*fr_tmp)*bd + bp*bm*bd*(rhoEr - rhoEl); + flux_hll[imomtt] = (bp*fl_tmp - bm*fr_tmp)*bd + bp*bm*bd*(qr[QRHO]*qr[iveltt] - ql[QRHO]*ql[iveltt]); + // total energy flux + Real rhoEl = ql[QREINT] + 0.5_rt*ql[QRHO]*(ql[ivel]*ql[ivel] + ql[ivelt]*ql[ivelt] + ql[iveltt]*ql[iveltt]); + fl_tmp = ql[ivel]*(rhoEl + ql[QPRES]); - // eint flux - fl_tmp = ql[QREINT]*ql[ivel]; - fr_tmp = qr[QREINT]*qr[ivel]; + Real rhoEr = qr[QREINT] + 0.5_rt*qr[QRHO]*(qr[ivel]*qr[ivel] + qr[ivelt]*qr[ivelt] + qr[iveltt]*qr[iveltt]); + fr_tmp = qr[ivel]*(rhoEr + qr[QPRES]); - flux_hll[UEINT] = (bp*fl_tmp - bm*fr_tmp)*bd + bp*bm*bd*(qr[QREINT] - ql[QREINT]); + flux_hll[UEDEN] = (bp*fl_tmp - bm*fr_tmp)*bd + bp*bm*bd*(rhoEr - rhoEl); - // passively-advected scalar fluxes - for (int ipassive = 0; ipassive < npassive; ipassive++) { - const int n = upassmap(ipassive); - const int nqs = qpassmap(ipassive); + // eint flux + fl_tmp = ql[QREINT]*ql[ivel]; + fr_tmp = qr[QREINT]*qr[ivel]; - fl_tmp = ql[QRHO]*ql[nqs]*ql[ivel]; - fr_tmp = qr[QRHO]*qr[nqs]*qr[ivel]; + flux_hll[UEINT] = (bp*fl_tmp - bm*fr_tmp)*bd + bp*bm*bd*(qr[QREINT] - ql[QREINT]); - flux_hll[n] = (bp*fl_tmp - bm*fr_tmp)*bd + bp*bm*bd*(qr[QRHO]*qr[nqs] - ql[QRHO]*ql[nqs]); - } -} + // passively-advected scalar fluxes + for (int ipassive = 0; ipassive < npassive; ipassive++) { + const int n = upassmap(ipassive); + const int nqs = qpassmap(ipassive); -/// -/// A HLLC Riemann solver for pure hydrodynamics -/// -/// @param ql the left interface state -/// @param qr the right interface state -/// @param qaux_arr the auxiliary state -/// @param uflx the flux through the interface -/// @param qint an approximate Godunov state on the interface -/// @param idir coordinate direction for the solve (0 = x, 1 = y, 2 = z) -/// -AMREX_GPU_HOST_DEVICE AMREX_INLINE -void -HLLC(const int i, const int j, const int k, const int idir, - Array4 const& ql, - Array4 const& qr, - Array4 const& qaux_arr, - Array4 const& uflx, - Array4 const& qgdnv, const bool store_full_state, - const GeometryData& geom, - const bool special_bnd_lo, const bool special_bnd_hi, - GpuArray const& domlo, GpuArray const& domhi) { - - // this is an implementation of the HLLC solver described in Toro's - // book. it uses the simplest estimate of the wave speeds, since - // those should work for a general EOS. We also initially do the - // CGF Riemann construction to get pstar and ustar, since we'll - // need to know the pressure and velocity on the interface for the - // pdV term in the internal energy update. - - const Real small = 1.e-8_rt; - - int iu; - int sx, sy, sz; - - if (idir == 0) { - iu = QU; - sx = 1; - sy = 0; - sz = 0; - - } else if (idir == 1) { - iu = QV; - sx = 0; - sy = 1; - sz = 0; - - } else { - iu = QW; - sx = 0; - sy = 0; - sz = 1; + fl_tmp = ql[QRHO]*ql[nqs]*ql[ivel]; + fr_tmp = qr[QRHO]*qr[nqs]*qr[ivel]; + flux_hll[n] = (bp*fl_tmp - bm*fr_tmp)*bd + bp*bm*bd*(qr[QRHO]*qr[nqs] - ql[QRHO]*ql[nqs]); + } } - int coord = geom.Coord(); - - // deal with hard walls - Real bnd_fac = 1.0_rt; - - if (idir == 0) { - if ((i == domlo[0] && special_bnd_lo) || - (i == domhi[0]+1 && special_bnd_hi)) { - bnd_fac = 0.0_rt; - } - - } else if (idir == 1) { - if ((j == domlo[1] && special_bnd_lo) || - (j == domhi[1]+1 && special_bnd_hi)) { - bnd_fac = 0.0_rt; - } + /// + /// A HLLC Riemann solver for pure hydrodynamics + /// + /// @param ql the left interface state + /// @param qr the right interface state + /// @param qaux_arr the auxiliary state + /// @param uflx the flux through the interface + /// @param qint an approximate Godunov state on the interface + /// @param idir coordinate direction for the solve (0 = x, 1 = y, 2 = z) + /// + AMREX_GPU_HOST_DEVICE AMREX_INLINE + void + HLLC(const int i, const int j, const int k, const int idir, + Array4 const& ql, + Array4 const& qr, + Array4 const& qaux_arr, + Array4 const& uflx, + Array4 const& qgdnv, const bool store_full_state, + const GeometryData& geom, + const bool special_bnd_lo, const bool special_bnd_hi, + GpuArray const& domlo, GpuArray const& domhi) { + + // this is an implementation of the HLLC solver described in Toro's + // book. it uses the simplest estimate of the wave speeds, since + // those should work for a general EOS. We also initially do the + // CGF Riemann construction to get pstar and ustar, since we'll + // need to know the pressure and velocity on the interface for the + // pdV term in the internal energy update. + + const Real small = 1.e-8_rt; + + int iu; + int sx, sy, sz; + + if (idir == 0) { + iu = QU; + sx = 1; + sy = 0; + sz = 0; + + } else if (idir == 1) { + iu = QV; + sx = 0; + sy = 1; + sz = 0; + + } else { + iu = QW; + sx = 0; + sy = 0; + sz = 1; - } else { - if ((k == domlo[2] && special_bnd_lo) || - (k == domhi[2]+1 && special_bnd_hi)) { - bnd_fac = 0.0_rt; } - } - - - Real rl = amrex::max(ql(i,j,k,QRHO), small_dens); - - // pick left velocities based on direction - Real ul = ql(i,j,k,iu); - - Real pl = amrex::max(ql(i,j,k,QPRES), small_pres); - Real rr = amrex::max(qr(i,j,k,QRHO), small_dens); - // pick right velocities based on direction - Real ur = qr(i,j,k,iu); + int coord = geom.Coord(); - Real pr = amrex::max(qr(i,j,k,QPRES), small_pres); + // deal with hard walls + Real bnd_fac = 1.0_rt; - // now we essentially do the CGF solver to get p and u on the - // interface, but we won't use these in any flux construction. - Real csmall = amrex::max(small, amrex::max(small * qaux_arr(i,j,k,QC), - small * qaux_arr(i-sx,j-sy,k-sz,QC))); - Real cavg = 0.5_rt*(qaux_arr(i,j,k,QC) + qaux_arr(i-sx,j-sy,k-sz,QC)); + if (idir == 0) { + if ((i == domlo[0] && special_bnd_lo) || + (i == domhi[0]+1 && special_bnd_hi)) { + bnd_fac = 0.0_rt; + } - Real gamcl = qaux_arr(i-sx,j-sy,k-sz,QGAMC); - Real gamcr = qaux_arr(i,j,k,QGAMC); + } else if (idir == 1) { + if ((j == domlo[1] && special_bnd_lo) || + (j == domhi[1]+1 && special_bnd_hi)) { + bnd_fac = 0.0_rt; + } -#ifdef TRUE_SDC - if (use_reconstructed_gamma1 == 1) { - gamcl = ql(i,j,k,QGC); - gamcr = qr(i,j,k,QGC); - } -#endif + } else { + if ((k == domlo[2] && special_bnd_lo) || + (k == domhi[2]+1 && special_bnd_hi)) { + bnd_fac = 0.0_rt; + } + } - Real wsmall = small_dens*csmall; - Real wl = amrex::max(wsmall, std::sqrt(std::abs(gamcl*pl*rl))); - Real wr = amrex::max(wsmall, std::sqrt(std::abs(gamcr*pr*rr))); - Real wwinv = 1.0_rt/(wl + wr); - Real pstar = ((wr*pl + wl*pr) + wl*wr*(ul - ur))*wwinv; - Real ustar = ((wl*ul + wr*ur) + (pl - pr))*wwinv; + Real rl = amrex::max(ql(i,j,k,QRHO), small_dens); - pstar = amrex::max(pstar, small_pres); - // for symmetry preservation, if ustar is really small, then we - // set it to zero - if (std::abs(ustar) < riemann_constants::smallu*0.5_rt*(std::abs(ul) + std::abs(ur))){ - ustar = 0.0_rt; - } + // pick left velocities based on direction + Real ul = ql(i,j,k,iu); - Real ro; - Real uo; - Real po; - Real gamco; - - if (ustar > 0.0_rt) { - ro = rl; - uo = ul; - po = pl; - gamco = gamcl; - - } else if (ustar < 0.0_rt) { - ro = rr; - uo = ur; - po = pr; - gamco = gamcr; - - } else { - ro = 0.5_rt*(rl + rr); - uo = 0.5_rt*(ul + ur); - po = 0.5_rt*(pl + pr); - gamco = 0.5_rt*(gamcl + gamcr); - } + Real pl = amrex::max(ql(i,j,k,QPRES), small_pres); - ro = amrex::max(small_dens, ro); + Real rr = amrex::max(qr(i,j,k,QRHO), small_dens); - Real roinv = 1.0_rt/ro; - Real co = std::sqrt(std::abs(gamco*po*roinv)); - co = amrex::max(csmall, co); - Real co2inv = 1.0_rt/(co*co); + // pick right velocities based on direction + Real ur = qr(i,j,k,iu); - Real rstar = ro + (pstar - po)*co2inv; - rstar = amrex::max(small_dens, rstar); + Real pr = amrex::max(qr(i,j,k,QPRES), small_pres); - Real cstar = std::sqrt(std::abs(gamco*pstar/rstar)); - cstar = max(cstar, csmall); + // now we essentially do the CGF solver to get p and u on the + // interface, but we won't use these in any flux construction. + Real csmall = amrex::max(small, amrex::max(small * qaux_arr(i,j,k,QC), + small * qaux_arr(i-sx,j-sy,k-sz,QC))); + Real cavg = 0.5_rt*(qaux_arr(i,j,k,QC) + qaux_arr(i-sx,j-sy,k-sz,QC)); - Real sgnm = std::copysign(1.0_rt, ustar); - Real spout = co - sgnm*uo; - Real spin = cstar - sgnm*ustar; - Real ushock = 0.5_rt*(spin + spout); + Real gamcl = qaux_arr(i-sx,j-sy,k-sz,QGAMC); + Real gamcr = qaux_arr(i,j,k,QGAMC); - if (pstar-po > 0.0_rt) { - spin = ushock; - spout = ushock; - } + #ifdef TRUE_SDC + if (use_reconstructed_gamma1 == 1) { + gamcl = ql(i,j,k,QGC); + gamcr = qr(i,j,k,QGC); + } + #endif - Real scr = spout-spin; - if (spout-spin == 0.0_rt) { - scr = small*cavg; - } + Real wsmall = small_dens*csmall; + Real wl = amrex::max(wsmall, std::sqrt(std::abs(gamcl*pl*rl))); + Real wr = amrex::max(wsmall, std::sqrt(std::abs(gamcr*pr*rr))); - Real frac = (1.0_rt + (spout + spin)/scr)*0.5_rt; - frac = amrex::max(0.0_rt, amrex::min(1.0_rt, frac)); + Real wwinv = 1.0_rt/(wl + wr); + Real pstar = ((wr*pl + wl*pr) + wl*wr*(ul - ur))*wwinv; + Real ustar = ((wl*ul + wr*ur) + (pl - pr))*wwinv; - Real qint[NQ] = {0.0}; + pstar = amrex::max(pstar, small_pres); + // for symmetry preservation, if ustar is really small, then we + // set it to zero + if (std::abs(ustar) < riemann_constants::smallu*0.5_rt*(std::abs(ul) + std::abs(ur))){ + ustar = 0.0_rt; + } - qint[QRHO] = frac*rstar + (1.0_rt - frac)*ro; - qint[iu] = frac*ustar + (1.0_rt - frac)*uo; - qint[QPRES] = frac*pstar + (1.0_rt - frac)*po; + Real ro; + Real uo; + Real po; + Real gamco; + + if (ustar > 0.0_rt) { + ro = rl; + uo = ul; + po = pl; + gamco = gamcl; + + } else if (ustar < 0.0_rt) { + ro = rr; + uo = ur; + po = pr; + gamco = gamcr; + + } else { + ro = 0.5_rt*(rl + rr); + uo = 0.5_rt*(ul + ur); + po = 0.5_rt*(pl + pr); + gamco = 0.5_rt*(gamcl + gamcr); + } + ro = amrex::max(small_dens, ro); - // now we do the HLLC construction + Real roinv = 1.0_rt/ro; + Real co = std::sqrt(std::abs(gamco*po*roinv)); + co = amrex::max(csmall, co); + Real co2inv = 1.0_rt/(co*co); - // use the simplest estimates of the wave speeds - Real S_l = amrex::min(ul - std::sqrt(gamcl*pl/rl), ur - std::sqrt(gamcr*pr/rr)); - Real S_r = amrex::max(ul + std::sqrt(gamcl*pl/rl), ur + std::sqrt(gamcr*pr/rr)); + Real rstar = ro + (pstar - po)*co2inv; + rstar = amrex::max(small_dens, rstar); - // estimate of the contact speed -- this is Toro Eq. 10.8 - Real S_c = (pr - pl + rl*ul*(S_l - ul) - rr*ur*(S_r - ur))/ - (rl*(S_l - ul) - rr*(S_r - ur)); + Real cstar = std::sqrt(std::abs(gamco*pstar/rstar)); + cstar = max(cstar, csmall); - Real q_zone[NQ]; - Real U_state[NUM_STATE]; - Real U_hllc_state[NUM_STATE]; - Real F_state[NUM_STATE]; + Real sgnm = std::copysign(1.0_rt, ustar); + Real spout = co - sgnm*uo; + Real spin = cstar - sgnm*ustar; + Real ushock = 0.5_rt*(spin + spout); - if (S_r <= 0.0_rt) { - // R region - for (int n = 0; n < NQ; n++) { - q_zone[n] = qr(i,j,k,n); - } - cons_state(q_zone, U_state); - compute_flux(idir, bnd_fac, coord, - U_state, pr, F_state); - - } else if (S_r > 0.0_rt && S_c <= 0.0_rt) { - // R* region - for (int n = 0; n < NQ; n++) { - q_zone[n] = qr(i,j,k,n); + if (pstar-po > 0.0_rt) { + spin = ushock; + spout = ushock; } - cons_state(q_zone, U_state); - compute_flux(idir, bnd_fac, coord, - U_state, pr, F_state); - HLLC_state(idir, S_r, S_c, q_zone, U_hllc_state); - // correct the flux - for (int n = 0; n < NUM_STATE; n++) { - F_state[n] = F_state[n] + S_r*(U_hllc_state[n] - U_state[n]); + Real scr = spout-spin; + if (spout-spin == 0.0_rt) { + scr = small*cavg; } - } else if (S_c > 0.0_rt && S_l < 0.0_rt) { - // L* region - for (int n = 0; n < NQ; n++) { - q_zone[n] = ql(i,j,k,n); + Real frac = (1.0_rt + (spout + spin)/scr)*0.5_rt; + frac = amrex::max(0.0_rt, amrex::min(1.0_rt, frac)); + + Real qint[NQ] = {0.0}; + + qint[QRHO] = frac*rstar + (1.0_rt - frac)*ro; + qint[iu] = frac*ustar + (1.0_rt - frac)*uo; + qint[QPRES] = frac*pstar + (1.0_rt - frac)*po; + + + // now we do the HLLC construction + + // use the simplest estimates of the wave speeds + Real S_l = amrex::min(ul - std::sqrt(gamcl*pl/rl), ur - std::sqrt(gamcr*pr/rr)); + Real S_r = amrex::max(ul + std::sqrt(gamcl*pl/rl), ur + std::sqrt(gamcr*pr/rr)); + + // estimate of the contact speed -- this is Toro Eq. 10.8 + Real S_c = (pr - pl + rl*ul*(S_l - ul) - rr*ur*(S_r - ur))/ + (rl*(S_l - ul) - rr*(S_r - ur)); + + Real q_zone[NQ]; + Real U_state[NUM_STATE]; + Real U_hllc_state[NUM_STATE]; + Real F_state[NUM_STATE]; + + if (S_r <= 0.0_rt) { + // R region + for (int n = 0; n < NQ; n++) { + q_zone[n] = qr(i,j,k,n); + } + cons_state(q_zone, U_state); + compute_flux(idir, bnd_fac, coord, + U_state, pr, F_state); + + } else if (S_r > 0.0_rt && S_c <= 0.0_rt) { + // R* region + for (int n = 0; n < NQ; n++) { + q_zone[n] = qr(i,j,k,n); + } + cons_state(q_zone, U_state); + compute_flux(idir, bnd_fac, coord, + U_state, pr, F_state); + HLLC_state(idir, S_r, S_c, q_zone, U_hllc_state); + + // correct the flux + for (int n = 0; n < NUM_STATE; n++) { + F_state[n] = F_state[n] + S_r*(U_hllc_state[n] - U_state[n]); + } + + } else if (S_c > 0.0_rt && S_l < 0.0_rt) { + // L* region + for (int n = 0; n < NQ; n++) { + q_zone[n] = ql(i,j,k,n); + } + cons_state(q_zone, U_state); + compute_flux(idir, bnd_fac, coord, + U_state, pl, F_state); + HLLC_state(idir, S_l, S_c, q_zone, U_hllc_state); + + // correct the flux + for (int n = 0; n < NUM_STATE; n++) { + F_state[n] = F_state[n] + S_l*(U_hllc_state[n] - U_state[n]); + } + + } else { + // L region + for (int n = 0; n < NQ; n++) { + q_zone[n] = ql(i,j,k,n); + } + cons_state(q_zone, U_state); + compute_flux(idir, bnd_fac, coord, + U_state, pl, F_state); } - cons_state(q_zone, U_state); - compute_flux(idir, bnd_fac, coord, - U_state, pl, F_state); - HLLC_state(idir, S_l, S_c, q_zone, U_hllc_state); - // correct the flux for (int n = 0; n < NUM_STATE; n++) { - F_state[n] = F_state[n] + S_l*(U_hllc_state[n] - U_state[n]); + uflx(i,j,k,n) = F_state[n]; } - } else { - // L region - for (int n = 0; n < NQ; n++) { - q_zone[n] = ql(i,j,k,n); - } - cons_state(q_zone, U_state); - compute_flux(idir, bnd_fac, coord, - U_state, pl, F_state); - } - for (int n = 0; n < NUM_STATE; n++) { - uflx(i,j,k,n) = F_state[n]; - } + // now store the state + if (store_full_state) { - // now store the state + for (int n = 0; n < NQ; n++) { + qgdnv(i,j,k,n) = qint[n]; + } - if (store_full_state) { + } else { - for (int n = 0; n < NQ; n++) { - qgdnv(i,j,k,n) = qint[n]; - } + #ifdef HYBRID_MOMENTUM + qgdnv(i,j,k,GDRHO) = qint[QRHO]; + #endif + qgdnv(i,j,k,GDU) = qint[QU]; + qgdnv(i,j,k,GDV) = qint[QV]; + qgdnv(i,j,k,GDW) = qint[QW]; + qgdnv(i,j,k,GDPRES) = qint[QPRES]; - } else { + } -#ifdef HYBRID_MOMENTUM - qgdnv(i,j,k,GDRHO) = qint[QRHO]; -#endif - qgdnv(i,j,k,GDU) = qint[QU]; - qgdnv(i,j,k,GDV) = qint[QV]; - qgdnv(i,j,k,GDW) = qint[QW]; - qgdnv(i,j,k,GDPRES) = qint[QPRES]; } - - } - #endif diff --git a/Source/hydro/riemann.cpp b/Source/hydro/riemann.cpp index f674b954a8..dc6eb0df98 100644 --- a/Source/hydro/riemann.cpp +++ b/Source/hydro/riemann.cpp @@ -122,14 +122,14 @@ Castro::cmpflx_plus_godunov(const Box& bx, } else if (riemann_solver == 2) { // HLLC - HLLC(i, j, k, idir, - qm, qp, - qaux_arr, - flx, - qgdnv, store_full_state, - geomdata, - special_bnd_lo, special_bnd_hi, - domlo, domhi); + HLL::HLLC(i, j, k, idir, + qm, qp, + qaux_arr, + flx, + qgdnv, store_full_state, + geomdata, + special_bnd_lo, special_bnd_hi, + domlo, domhi); #ifndef AMREX_USE_GPU } else { @@ -182,9 +182,9 @@ Castro::cmpflx_plus_godunov(const Box& bx, flx_zone[n] = flx(i,j,k,n); } - HLL(ql_zone, qr_zone, cl, cr, - idir, coord, - flx_zone); + HLL::HLL(ql_zone, qr_zone, cl, cr, + idir, coord, + flx_zone); for (int n = 0; n < NUM_STATE; n++) { flx(i,j,k,n) = flx_zone[n]; From 8fdcbf54bef10543e0d397747b9b6abfcdef9a87 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Thu, 27 Jun 2024 09:45:29 -0400 Subject: [PATCH 5/5] add namespace --- Source/hydro/riemann_2shock_solvers.H | 1185 +++++++++++++------------ Source/hydro/riemann_solvers.H | 6 +- 2 files changed, 595 insertions(+), 596 deletions(-) diff --git a/Source/hydro/riemann_2shock_solvers.H b/Source/hydro/riemann_2shock_solvers.H index a1064545c9..4f1485eef8 100644 --- a/Source/hydro/riemann_2shock_solvers.H +++ b/Source/hydro/riemann_2shock_solvers.H @@ -11,741 +11,742 @@ using namespace amrex::literals; #include #endif +namespace TwoShock { -AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -void -wsqge(const amrex::Real p, const amrex::Real v, - const amrex::Real gam, const amrex::Real gdot, amrex::Real& gstar, - const amrex::Real gmin, const amrex::Real gmax, const amrex::Real csq, - const amrex::Real pstar, amrex::Real& wsq) { + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE + void + wsqge(const amrex::Real p, const amrex::Real v, + const amrex::Real gam, const amrex::Real gdot, amrex::Real& gstar, + const amrex::Real gmin, const amrex::Real gmax, const amrex::Real csq, + const amrex::Real pstar, amrex::Real& wsq) { - // compute the lagrangian wave speeds -- this is the approximate - // version for the Colella & Glaz algorithm + // compute the lagrangian wave speeds -- this is the approximate + // version for the Colella & Glaz algorithm - // First predict a value of game across the shock + // First predict a value of game across the shock - // CG Eq. 31 - gstar = (pstar-p)*gdot/(pstar+p) + gam; - gstar = amrex::max(gmin, amrex::min(gmax, gstar)); + // CG Eq. 31 + gstar = (pstar-p)*gdot/(pstar+p) + gam; + gstar = amrex::max(gmin, amrex::min(gmax, gstar)); - // Now use that predicted value of game with the R-H jump conditions - // to compute the wave speed. + // Now use that predicted value of game with the R-H jump conditions + // to compute the wave speed. - // this is CG Eq. 34 - amrex::Real alpha = pstar - (gstar - 1.0_rt)*p/(gam - 1.0_rt); - if (alpha == 0.0_rt) { - alpha = riemann_constants::smlp1*(pstar + p); - } + // this is CG Eq. 34 + amrex::Real alpha = pstar - (gstar - 1.0_rt)*p/(gam - 1.0_rt); + if (alpha == 0.0_rt) { + alpha = riemann_constants::smlp1*(pstar + p); + } - amrex::Real beta = pstar + 0.5_rt*(gstar - 1.0_rt)*(pstar+p); + amrex::Real beta = pstar + 0.5_rt*(gstar - 1.0_rt)*(pstar+p); - wsq = (pstar-p)*beta/(v*alpha); + wsq = (pstar-p)*beta/(v*alpha); - if (std::abs(pstar - p) < riemann_constants::smlp1*(pstar + p)) { - wsq = csq; - } - wsq = amrex::max(wsq, (0.5_rt * (gam - 1.0_rt)/gam)*csq); -} + if (std::abs(pstar - p) < riemann_constants::smlp1*(pstar + p)) { + wsq = csq; + } + wsq = amrex::max(wsq, (0.5_rt * (gam - 1.0_rt)/gam)*csq); + } -AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -void -pstar_bisection(amrex::Real& pstar_lo, amrex::Real& pstar_hi, - const amrex::Real ul, const amrex::Real pl, const amrex::Real taul, - const amrex::Real gamel, const amrex::Real clsql, - const amrex::Real ur, const amrex::Real pr, const amrex::Real taur, - const amrex::Real gamer, const amrex::Real clsqr, - const amrex::Real gdot, const amrex::Real gmin, const amrex::Real gmax, - const int lcg_maxiter, const amrex::Real lcg_tol, - amrex::Real& pstar, amrex::Real& gamstar, - bool& converged, amrex::GpuArray& pstar_hist_extra) { + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE + void + pstar_bisection(amrex::Real& pstar_lo, amrex::Real& pstar_hi, + const amrex::Real ul, const amrex::Real pl, const amrex::Real taul, + const amrex::Real gamel, const amrex::Real clsql, + const amrex::Real ur, const amrex::Real pr, const amrex::Real taur, + const amrex::Real gamer, const amrex::Real clsqr, + const amrex::Real gdot, const amrex::Real gmin, const amrex::Real gmax, + const int lcg_maxiter, const amrex::Real lcg_tol, + amrex::Real& pstar, amrex::Real& gamstar, + bool& converged, amrex::GpuArray& pstar_hist_extra) { - // we want to zero - // f(p*) = u*_l(p*) - u*_r(p*) - // we'll do bisection - // - // this version is for the approximate Colella & Glaz - // version + // we want to zero + // f(p*) = u*_l(p*) - u*_r(p*) + // we'll do bisection + // + // this version is for the approximate Colella & Glaz + // version - // lo bounds - amrex::Real wlsq = 0.0; - wsqge(pl, taul, gamel, gdot, - gamstar, gmin, gmax, clsql, pstar_lo, wlsq); + // lo bounds + amrex::Real wlsq = 0.0; + wsqge(pl, taul, gamel, gdot, + gamstar, gmin, gmax, clsql, pstar_lo, wlsq); - amrex::Real wrsq = 0.0; - wsqge(pr, taur, gamer, gdot, - gamstar, gmin, gmax, clsqr, pstar_lo, wrsq); + amrex::Real wrsq = 0.0; + wsqge(pr, taur, gamer, gdot, + gamstar, gmin, gmax, clsqr, pstar_lo, wrsq); - amrex::Real wl = 1.0_rt / std::sqrt(wlsq); - amrex::Real wr = 1.0_rt / std::sqrt(wrsq); + amrex::Real wl = 1.0_rt / std::sqrt(wlsq); + amrex::Real wr = 1.0_rt / std::sqrt(wrsq); - amrex::Real ustar_l = ul - (pstar_lo - pstar)*wl; - amrex::Real ustar_r = ur + (pstar_lo - pstar)*wr; + amrex::Real ustar_l = ul - (pstar_lo - pstar)*wl; + amrex::Real ustar_r = ur + (pstar_lo - pstar)*wr; - amrex::Real f_lo = ustar_l - ustar_r; + amrex::Real f_lo = ustar_l - ustar_r; - // hi bounds - wsqge(pl, taul, gamel, gdot, - gamstar, gmin, gmax, clsql, pstar_hi, wlsq); + // hi bounds + wsqge(pl, taul, gamel, gdot, + gamstar, gmin, gmax, clsql, pstar_hi, wlsq); - wsqge(pr, taur, gamer, gdot, - gamstar, gmin, gmax, clsqr, pstar_hi, wrsq); + wsqge(pr, taur, gamer, gdot, + gamstar, gmin, gmax, clsqr, pstar_hi, wrsq); - // wl = 1.0_rt / std::sqrt(wlsq); - // wr = 1.0_rt / std::sqrt(wrsq); + // wl = 1.0_rt / std::sqrt(wlsq); + // wr = 1.0_rt / std::sqrt(wrsq); - //ustar_l = ul - (pstar_hi - pstar)*wl; - //ustar_r = ur + (pstar_hi - pstar)*wr; + //ustar_l = ul - (pstar_hi - pstar)*wl; + //ustar_r = ur + (pstar_hi - pstar)*wr; - //amrex::Real f_hi = ustar_l - ustar_r; + //amrex::Real f_hi = ustar_l - ustar_r; - // bisection - converged = false; - amrex::Real pstar_c = 0.0; + // bisection + converged = false; + amrex::Real pstar_c = 0.0; - for (int iter = 0; iter < PSTAR_BISECT_FACTOR*lcg_maxiter; iter++) { + for (int iter = 0; iter < PSTAR_BISECT_FACTOR*lcg_maxiter; iter++) { - pstar_c = 0.5_rt * (pstar_lo + pstar_hi); - pstar_hist_extra[iter] = pstar_c; + pstar_c = 0.5_rt * (pstar_lo + pstar_hi); + pstar_hist_extra[iter] = pstar_c; - wsqge(pl, taul, gamel, gdot, - gamstar, gmin, gmax, clsql, pstar_c, wlsq); + wsqge(pl, taul, gamel, gdot, + gamstar, gmin, gmax, clsql, pstar_c, wlsq); - wsqge(pr, taur, gamer, gdot, - gamstar, gmin, gmax, clsqr, pstar_c, wrsq); + wsqge(pr, taur, gamer, gdot, + gamstar, gmin, gmax, clsqr, pstar_c, wrsq); - wl = 1.0_rt / std::sqrt(wlsq); - wr = 1.0_rt / std::sqrt(wrsq); + wl = 1.0_rt / std::sqrt(wlsq); + wr = 1.0_rt / std::sqrt(wrsq); - ustar_l = ul - (pstar_c - pl)*wl; - ustar_r = ur - (pstar_c - pr)*wr; + ustar_l = ul - (pstar_c - pl)*wl; + ustar_r = ur - (pstar_c - pr)*wr; - amrex::Real f_c = ustar_l - ustar_r; + amrex::Real f_c = ustar_l - ustar_r; - if ( 0.5_rt * std::abs(pstar_lo - pstar_hi) < lcg_tol * pstar_c ) { - converged = true; - break; - } + if ( 0.5_rt * std::abs(pstar_lo - pstar_hi) < lcg_tol * pstar_c ) { + converged = true; + break; + } + + if (f_lo * f_c < 0.0_rt) { + // root is in the left half + pstar_hi = pstar_c; + //f_hi = f_c; + } else { + pstar_lo = pstar_c; + f_lo = f_c; + } + } - if (f_lo * f_c < 0.0_rt) { - // root is in the left half - pstar_hi = pstar_c; - //f_hi = f_c; - } else { - pstar_lo = pstar_c; - f_lo = f_c; + pstar = pstar_c; } - } - pstar = pstar_c; -} + /// + /// The Colella-Glaz Riemann solver for pure hydrodynamics. This is a + /// two shock approximate state Riemann solver. + /// + /// @param bx the box to operate over + /// @param ql the left interface state + /// @param qr the right interface state + /// @param qaux_arr the auxiliary state + /// @param qint the full Godunov state on the interface + /// @param idir coordinate direction for the solve (0 = x, 1 = y, 2 = z) + /// + AMREX_GPU_HOST_DEVICE AMREX_INLINE + void + riemanncg(const RiemannState& ql, const RiemannState& qr, const RiemannAux& raux, + RiemannState& qint) { + + // this implements the approximate Riemann solver of Colella & Glaz + // (1985) + // + + constexpr amrex::Real weakwv = 1.e-3_rt; + + #ifndef AMREX_USE_GPU + amrex::GpuArray pstar_hist; + #endif + + + // common quantities + amrex::Real taul = 1.0_rt / ql.rho; + amrex::Real taur = 1.0_rt / qr.rho; + + // lagrangian sound speeds + amrex::Real clsql = ql.gamc * ql.p * ql.rho; + amrex::Real clsqr = qr.gamc * qr.p * qr.rho; + + // Note: in the original Colella & Glaz paper, they predicted + // gamma_e to the interfaces using a special (non-hyperbolic) + // evolution equation. In Castro, we instead bring (rho e) + // to the edges, so we construct the necessary gamma_e here from + // what we have on the interfaces. + amrex::Real gamel = ql.p / ql.rhoe + 1.0_rt; + amrex::Real gamer = qr.p / qr.rhoe + 1.0_rt; + + // these should consider a wider average of the cell-centered + // gammas + amrex::Real gmin = amrex::min(amrex::min(gamel, gamer), 1.0_rt); + amrex::Real gmax = amrex::max(amrex::max(gamel, gamer), 2.0_rt); + + amrex::Real game_bar = 0.5_rt*(gamel + gamer); + amrex::Real gamc_bar = 0.5_rt*(ql.gamc + qr.gamc); + + amrex::Real gdot = 2.0_rt*(1.0_rt - game_bar/gamc_bar)*(game_bar - 1.0_rt); + + amrex::Real wsmall = small_dens * raux.csmall; + amrex::Real wl = amrex::max(wsmall, std::sqrt(std::abs(clsql))); + amrex::Real wr = amrex::max(wsmall, std::sqrt(std::abs(clsqr))); + + // make an initial guess for pstar -- this is a two-shock + // approximation + //pstar = ((wr*pl + wl*pr) + wl*wr*(ul - ur))/(wl + wr) + amrex::Real pstar = ql.p + ( (qr.p - ql.p) - wr*(qr.un - ql.un) ) * wl / (wl + wr); + pstar = amrex::max(pstar, small_pres); -/// -/// The Colella-Glaz Riemann solver for pure hydrodynamics. This is a -/// two shock approximate state Riemann solver. -/// -/// @param bx the box to operate over -/// @param ql the left interface state -/// @param qr the right interface state -/// @param qaux_arr the auxiliary state -/// @param qint the full Godunov state on the interface -/// @param idir coordinate direction for the solve (0 = x, 1 = y, 2 = z) -/// -AMREX_GPU_HOST_DEVICE AMREX_INLINE -void -riemanncg(const RiemannState& ql, const RiemannState& qr, const RiemannAux& raux, - RiemannState& qint) { - - // this implements the approximate Riemann solver of Colella & Glaz - // (1985) - // - - constexpr amrex::Real weakwv = 1.e-3_rt; - -#ifndef AMREX_USE_GPU - amrex::GpuArray pstar_hist; -#endif + // get the shock speeds -- this computes W_s from CG Eq. 34 + amrex::Real gamstar = 0.0; + amrex::Real wlsq = 0.0; + wsqge(ql.p, taul, gamel, gdot, gamstar, + gmin, gmax, clsql, pstar, wlsq); - // common quantities - amrex::Real taul = 1.0_rt / ql.rho; - amrex::Real taur = 1.0_rt / qr.rho; + amrex::Real wrsq = 0.0; + wsqge(qr.p, taur, gamer, gdot, gamstar, + gmin, gmax, clsqr, pstar, wrsq); - // lagrangian sound speeds - amrex::Real clsql = ql.gamc * ql.p * ql.rho; - amrex::Real clsqr = qr.gamc * qr.p * qr.rho; + amrex::Real pstar_old = pstar; - // Note: in the original Colella & Glaz paper, they predicted - // gamma_e to the interfaces using a special (non-hyperbolic) - // evolution equation. In Castro, we instead bring (rho e) - // to the edges, so we construct the necessary gamma_e here from - // what we have on the interfaces. - amrex::Real gamel = ql.p / ql.rhoe + 1.0_rt; - amrex::Real gamer = qr.p / qr.rhoe + 1.0_rt; + wl = std::sqrt(wlsq); + wr = std::sqrt(wrsq); - // these should consider a wider average of the cell-centered - // gammas - amrex::Real gmin = amrex::min(amrex::min(gamel, gamer), 1.0_rt); - amrex::Real gmax = amrex::max(amrex::max(gamel, gamer), 2.0_rt); + // R-H jump conditions give ustar across each wave -- these + // should be equal when we are done iterating. Our notation + // here is a little funny, comparing to CG, ustar_l = u*_L and + // ustar_r = u*_R. + amrex::Real ustar_l = ql.un - (pstar - ql.p) / wl; + amrex::Real ustar_r = qr.un + (pstar - qr.p) / wr; - amrex::Real game_bar = 0.5_rt*(gamel + gamer); - amrex::Real gamc_bar = 0.5_rt*(ql.gamc + qr.gamc); + // revise our pstar guess + // pstar = ((wr*pl + wl*pr) + wl*wr*(ul - ur))/(wl + wr) + pstar = ql.p + ( (qr.p - ql.p) - wr * (qr.un - ql.un) ) * wl / (wl + wr); + pstar = amrex::max(pstar, small_pres); - amrex::Real gdot = 2.0_rt*(1.0_rt - game_bar/gamc_bar)*(game_bar - 1.0_rt); + // secant iteration + bool converged = false; - amrex::Real wsmall = small_dens * raux.csmall; - amrex::Real wl = amrex::max(wsmall, std::sqrt(std::abs(clsql))); - amrex::Real wr = amrex::max(wsmall, std::sqrt(std::abs(clsqr))); + int iter = 0; + while ((iter < cg_maxiter && !converged) || iter < 2) { - // make an initial guess for pstar -- this is a two-shock - // approximation - //pstar = ((wr*pl + wl*pr) + wl*wr*(ul - ur))/(wl + wr) - amrex::Real pstar = ql.p + ( (qr.p - ql.p) - wr*(qr.un - ql.un) ) * wl / (wl + wr); - pstar = amrex::max(pstar, small_pres); + wsqge(ql.p, taul, gamel, gdot, gamstar, + gmin, gmax, clsql, pstar, wlsq); - // get the shock speeds -- this computes W_s from CG Eq. 34 - amrex::Real gamstar = 0.0; + wsqge(qr.p, taur, gamer, gdot, gamstar, + gmin, gmax, clsqr, pstar, wrsq); - amrex::Real wlsq = 0.0; - wsqge(ql.p, taul, gamel, gdot, gamstar, - gmin, gmax, clsql, pstar, wlsq); - amrex::Real wrsq = 0.0; - wsqge(qr.p, taur, gamer, gdot, gamstar, - gmin, gmax, clsqr, pstar, wrsq); + // NOTE: these are really the inverses of the wave speeds! + wl = 1.0_rt / std::sqrt(wlsq); + wr = 1.0_rt / std::sqrt(wrsq); - amrex::Real pstar_old = pstar; + amrex::Real ustar_r_old = ustar_r; + amrex::Real ustar_l_old = ustar_l; - wl = std::sqrt(wlsq); - wr = std::sqrt(wrsq); + ustar_r = qr.un - (qr.p - pstar) * wr; + ustar_l = ql.un + (ql.p - pstar) * wl; - // R-H jump conditions give ustar across each wave -- these - // should be equal when we are done iterating. Our notation - // here is a little funny, comparing to CG, ustar_l = u*_L and - // ustar_r = u*_R. - amrex::Real ustar_l = ql.un - (pstar - ql.p) / wl; - amrex::Real ustar_r = qr.un + (pstar - qr.p) / wr; + amrex::Real dpditer = std::abs(pstar_old - pstar); - // revise our pstar guess - // pstar = ((wr*pl + wl*pr) + wl*wr*(ul - ur))/(wl + wr) - pstar = ql.p + ( (qr.p - ql.p) - wr * (qr.un - ql.un) ) * wl / (wl + wr); - pstar = amrex::max(pstar, small_pres); + // Here we are going to do the Secant iteration version in + // CG. Note that what we call zp and zm here are not + // actually the Z_p = |dp*/du*_p| defined in CG, by rather + // simply |du*_p| (or something that looks like dp/Z!). + amrex::Real zp = std::abs(ustar_l - ustar_l_old); + if (zp - weakwv * raux.cavg <= 0.0_rt) { + zp = dpditer * wl; + } - // secant iteration - bool converged = false; + amrex::Real zm = std::abs(ustar_r - ustar_r_old); + if (zm - weakwv * raux.cavg <= 0.0_rt) { + zm = dpditer * wr; + } - int iter = 0; - while ((iter < cg_maxiter && !converged) || iter < 2) { + // the new pstar is found via CG Eq. 18 + amrex::Real denom = dpditer / amrex::max(zp + zm, riemann_constants::small * raux.cavg); + pstar_old = pstar; + pstar = pstar - denom*(ustar_r - ustar_l); + pstar = amrex::max(pstar, small_pres); - wsqge(ql.p, taul, gamel, gdot, gamstar, - gmin, gmax, clsql, pstar, wlsq); + amrex::Real err = std::abs(pstar - pstar_old); + if (err < cg_tol*pstar) { + converged = true; + } - wsqge(qr.p, taur, gamer, gdot, gamstar, - gmin, gmax, clsqr, pstar, wrsq); + #ifndef AMREX_USE_GPU + pstar_hist[iter] = pstar; + #endif + iter++; + } - // NOTE: these are really the inverses of the wave speeds! - wl = 1.0_rt / std::sqrt(wlsq); - wr = 1.0_rt / std::sqrt(wrsq); + // If we failed to converge using the secant iteration, we + // can either stop here; or, revert to the original + // two-shock estimate for pstar; or do a bisection root + // find using the bounds established by the most recent + // iterations. - amrex::Real ustar_r_old = ustar_r; - amrex::Real ustar_l_old = ustar_l; + if (!converged) { - ustar_r = qr.un - (qr.p - pstar) * wr; - ustar_l = ql.un + (ql.p - pstar) * wl; + if (cg_blend == 0) { - amrex::Real dpditer = std::abs(pstar_old - pstar); + #ifndef AMREX_USE_GPU + std::cout << "pstar history: " << std::endl; + for (int iter_l=0; iter_l < cg_maxiter; iter_l++) { + std::cout << iter_l << " " << pstar_hist[iter_l] << std::endl; + } - // Here we are going to do the Secant iteration version in - // CG. Note that what we call zp and zm here are not - // actually the Z_p = |dp*/du*_p| defined in CG, by rather - // simply |du*_p| (or something that looks like dp/Z!). - amrex::Real zp = std::abs(ustar_l - ustar_l_old); - if (zp - weakwv * raux.cavg <= 0.0_rt) { - zp = dpditer * wl; - } + std::cout << std::endl; + std::cout << "left state: " << std::endl << ql << std::endl; + std::cout << "right state: " << std::endl << qr << std::endl; + std::cout << "aux information: " << std::endl << raux << std::endl; - amrex::Real zm = std::abs(ustar_r - ustar_r_old); - if (zm - weakwv * raux.cavg <= 0.0_rt) { - zm = dpditer * wr; - } + amrex::Error("ERROR: non-convergence in the Riemann solver"); + #endif - // the new pstar is found via CG Eq. 18 - amrex::Real denom = dpditer / amrex::max(zp + zm, riemann_constants::small * raux.cavg); - pstar_old = pstar; - pstar = pstar - denom*(ustar_r - ustar_l); - pstar = amrex::max(pstar, small_pres); + } else if (cg_blend == 1) { - amrex::Real err = std::abs(pstar - pstar_old); - if (err < cg_tol*pstar) { - converged = true; - } + pstar = ql.p + ( (qr.p - ql.p) - wr * (qr.un - ql.un) ) * wl / (wl + wr); -#ifndef AMREX_USE_GPU - pstar_hist[iter] = pstar; -#endif + } else if (cg_blend == 2) { + + // we don't store the history if we are in CUDA, so + // we can't do this + #ifndef AMREX_USE_GPU + // first try to find a reasonable bounds + amrex::Real pstarl = 1.e200; + amrex::Real pstaru = -1.e200; + for (int n = cg_maxiter-6; n < cg_maxiter; n++) { + pstarl = amrex::min(pstarl, pstar_hist[n]); + pstaru = amrex::max(pstaru, pstar_hist[n]); + } + + pstarl = amrex::max(pstarl, small_pres); + pstaru = amrex::max(pstaru, small_pres); + + amrex::GpuArray pstar_hist_extra; - iter++; - } + pstar_bisection(pstarl, pstaru, + ql.un, ql.p, taul, gamel, clsql, + qr.un, qr.p, taur, gamer, clsqr, + gdot, gmin, gmax, + cg_maxiter, cg_tol, + pstar, gamstar, converged, pstar_hist_extra); - // If we failed to converge using the secant iteration, we - // can either stop here; or, revert to the original - // two-shock estimate for pstar; or do a bisection root - // find using the bounds established by the most recent - // iterations. + if (!converged) { - if (!converged) { + std::cout << "pstar history: " << std::endl; + for (int iter_l = 0; iter_l < cg_maxiter; iter_l++) { + std::cout << iter_l << " " << pstar_hist[iter_l] << std::endl; + } + std::cout << "pstar extra history: " << std::endl; + for (int iter_l = 0; iter_l < PSTAR_BISECT_FACTOR*cg_maxiter; iter_l++) { + std::cout << iter_l << " " << pstar_hist_extra[iter_l] << std::endl; + } - if (cg_blend == 0) { + std::cout << std::endl; + std::cout << "left state: " << std::endl << ql << std::endl; + std::cout << "right state: " << std::endl << qr << std::endl; + std::cout << "aux information: " << std::endl << raux << std::endl; -#ifndef AMREX_USE_GPU - std::cout << "pstar history: " << std::endl; - for (int iter_l=0; iter_l < cg_maxiter; iter_l++) { - std::cout << iter_l << " " << pstar_hist[iter_l] << std::endl; + amrex::Error("ERROR: non-convergence in the Riemann solver"); + } + + #endif + } else { + + #ifndef AMREX_USE_GPU + amrex::Error("ERROR: unrecognized cg_blend option."); + #endif } - std::cout << std::endl; - std::cout << "left state: " << std::endl << ql << std::endl; - std::cout << "right state: " << std::endl << qr << std::endl; - std::cout << "aux information: " << std::endl << raux << std::endl; + } - amrex::Error("ERROR: non-convergence in the Riemann solver"); -#endif + // we converged! construct the single ustar for the region + // between the left and right waves, using the updated wave speeds + ustar_r = qr.un - (qr.p - pstar) * wr; // careful -- here wl, wr are 1/W + ustar_l = ql.un + (ql.p - pstar) * wl; - } else if (cg_blend == 1) { + amrex::Real ustar = 0.5_rt * (ustar_l + ustar_r); - pstar = ql.p + ( (qr.p - ql.p) - wr * (qr.un - ql.un) ) * wl / (wl + wr); + // for symmetry preservation, if ustar is really small, then we + // set it to zero + if (std::abs(ustar) < riemann_constants::smallu * 0.5_rt * (std::abs(ql.un) + std::abs(qr.un))) { + ustar = 0.0_rt; + } - } else if (cg_blend == 2) { + // sample the solution -- here we look first at the direction + // that the contact is moving. This tells us if we need to + // worry about the L/L* states or the R*/R states. + amrex::Real ro; + amrex::Real uo; + amrex::Real po; + amrex::Real tauo; + amrex::Real gamco; + amrex::Real gameo; + + if (ustar > 0.0_rt) { + //ro = ql.rho; + uo = ql.un; + po = ql.p; + tauo = taul; + gamco = ql.gamc; + gameo = gamel; + + } else if (ustar < 0.0_rt) { + //ro = qr.rho; + uo = qr.un; + po = qr.p; + tauo = taur; + gamco = qr.gamc; + gameo = gamer; - // we don't store the history if we are in CUDA, so - // we can't do this -#ifndef AMREX_USE_GPU - // first try to find a reasonable bounds - amrex::Real pstarl = 1.e200; - amrex::Real pstaru = -1.e200; - for (int n = cg_maxiter-6; n < cg_maxiter; n++) { - pstarl = amrex::min(pstarl, pstar_hist[n]); - pstaru = amrex::max(pstaru, pstar_hist[n]); - } + } else { + //ro = 0.5_rt * (ql.rho + qr.rho); + uo = 0.5_rt * (ql.un + qr.un); + po = 0.5_rt * (ql.p + qr.p); + tauo = 0.5_rt * (taul + taur); + gamco = 0.5_rt * (ql.gamc + qr.gamc); + gameo = 0.5_rt * (gamel + gamer); + } - pstarl = amrex::max(pstarl, small_pres); - pstaru = amrex::max(pstaru, small_pres); + // use tau = 1/rho as the independent variable here + ro = amrex::max(small_dens, 1.0_rt/tauo); + tauo = 1.0_rt/ro; - amrex::GpuArray pstar_hist_extra; + amrex::Real co = std::sqrt(std::abs(gamco*po*tauo)); + co = amrex::max(raux.csmall, co); + amrex::Real clsq = std::pow(co*ro, 2); - pstar_bisection(pstarl, pstaru, - ql.un, ql.p, taul, gamel, clsql, - qr.un, qr.p, taur, gamer, clsqr, - gdot, gmin, gmax, - cg_maxiter, cg_tol, - pstar, gamstar, converged, pstar_hist_extra); + // now that we know which state (left or right) we need to worry + // about, get the value of gamstar and wosq across the wave we + // are dealing with. + amrex::Real wosq = 0.0; + wsqge(po, tauo, gameo, gdot, gamstar, + gmin, gmax, clsq, pstar, wosq); - if (!converged) { + amrex::Real sgnm = std::copysign(1.0_rt, ustar); - std::cout << "pstar history: " << std::endl; - for (int iter_l = 0; iter_l < cg_maxiter; iter_l++) { - std::cout << iter_l << " " << pstar_hist[iter_l] << std::endl; - } - std::cout << "pstar extra history: " << std::endl; - for (int iter_l = 0; iter_l < PSTAR_BISECT_FACTOR*cg_maxiter; iter_l++) { - std::cout << iter_l << " " << pstar_hist_extra[iter_l] << std::endl; - } + amrex::Real wo = std::sqrt(wosq); + amrex::Real dpjmp = pstar - po; - std::cout << std::endl; - std::cout << "left state: " << std::endl << ql << std::endl; - std::cout << "right state: " << std::endl << qr << std::endl; - std::cout << "aux information: " << std::endl << raux << std::endl; + // is this max really necessary? + //rstar=max(ONE-ro*dpjmp/wosq, (gameo-ONE)/(gameo+ONE)) + amrex::Real rstar = 1.0_rt - ro*dpjmp/wosq; + rstar = ro/rstar; + rstar = amrex::max(small_dens, rstar); - amrex::Error("ERROR: non-convergence in the Riemann solver"); - } + amrex::Real cstar = std::sqrt(std::abs(gamco * pstar / rstar)); + cstar = amrex::max(cstar, raux.csmall); -#endif + amrex::Real spout = co - sgnm*uo; + amrex::Real spin = cstar - sgnm*ustar; + + //ushock = 0.5_rt*(spin + spout) + amrex::Real ushock = wo*tauo - sgnm*uo; + + if (pstar - po >= 0.0_rt) { + spin = ushock; + spout = ushock; + } + + amrex::Real frac = 0.5_rt*(1.0_rt + (spin + spout)/amrex::max(amrex::max(spout-spin, spin+spout), + riemann_constants::small * raux.cavg)); + + // the transverse velocity states only depend on the + // direction that the contact moves + if (ustar > 0.0_rt) { + qint.ut = ql.ut; + qint.utt = ql.utt; + } else if (ustar < 0.0_rt) { + qint.ut = qr.ut; + qint.utt = qr.utt; } else { + qint.ut = 0.5_rt * (ql.ut + qr.ut); + qint.utt = 0.5_rt * (ql.utt + qr.utt); + } -#ifndef AMREX_USE_GPU - amrex::Error("ERROR: unrecognized cg_blend option."); -#endif + // linearly interpolate between the star and normal state -- this covers the + // case where we are inside the rarefaction fan. + qint.rho = frac*rstar + (1.0_rt - frac)*ro; + qint.un = frac*ustar + (1.0_rt - frac)*uo; + qint.p = frac*pstar + (1.0_rt - frac)*po; + amrex::Real game_int = frac*gamstar + (1.0_rt-frac)*gameo; + + // now handle the cases where instead we are fully in the + // star or fully in the original (l/r) state + if (spout < 0.0_rt) { + qint.rho = ro; + qint.un = uo; + qint.p = po; + game_int = gameo; } - } - - // we converged! construct the single ustar for the region - // between the left and right waves, using the updated wave speeds - ustar_r = qr.un - (qr.p - pstar) * wr; // careful -- here wl, wr are 1/W - ustar_l = ql.un + (ql.p - pstar) * wl; - - amrex::Real ustar = 0.5_rt * (ustar_l + ustar_r); - - // for symmetry preservation, if ustar is really small, then we - // set it to zero - if (std::abs(ustar) < riemann_constants::smallu * 0.5_rt * (std::abs(ql.un) + std::abs(qr.un))) { - ustar = 0.0_rt; - } - - // sample the solution -- here we look first at the direction - // that the contact is moving. This tells us if we need to - // worry about the L/L* states or the R*/R states. - amrex::Real ro; - amrex::Real uo; - amrex::Real po; - amrex::Real tauo; - amrex::Real gamco; - amrex::Real gameo; - - if (ustar > 0.0_rt) { - //ro = ql.rho; - uo = ql.un; - po = ql.p; - tauo = taul; - gamco = ql.gamc; - gameo = gamel; - - } else if (ustar < 0.0_rt) { - //ro = qr.rho; - uo = qr.un; - po = qr.p; - tauo = taur; - gamco = qr.gamc; - gameo = gamer; - - } else { - //ro = 0.5_rt * (ql.rho + qr.rho); - uo = 0.5_rt * (ql.un + qr.un); - po = 0.5_rt * (ql.p + qr.p); - tauo = 0.5_rt * (taul + taur); - gamco = 0.5_rt * (ql.gamc + qr.gamc); - gameo = 0.5_rt * (gamel + gamer); - } - - // use tau = 1/rho as the independent variable here - ro = amrex::max(small_dens, 1.0_rt/tauo); - tauo = 1.0_rt/ro; - - amrex::Real co = std::sqrt(std::abs(gamco*po*tauo)); - co = amrex::max(raux.csmall, co); - amrex::Real clsq = std::pow(co*ro, 2); - - // now that we know which state (left or right) we need to worry - // about, get the value of gamstar and wosq across the wave we - // are dealing with. - amrex::Real wosq = 0.0; - wsqge(po, tauo, gameo, gdot, gamstar, - gmin, gmax, clsq, pstar, wosq); - - amrex::Real sgnm = std::copysign(1.0_rt, ustar); - - amrex::Real wo = std::sqrt(wosq); - amrex::Real dpjmp = pstar - po; - - // is this max really necessary? - //rstar=max(ONE-ro*dpjmp/wosq, (gameo-ONE)/(gameo+ONE)) - amrex::Real rstar = 1.0_rt - ro*dpjmp/wosq; - rstar = ro/rstar; - rstar = amrex::max(small_dens, rstar); - - amrex::Real cstar = std::sqrt(std::abs(gamco * pstar / rstar)); - cstar = amrex::max(cstar, raux.csmall); - - amrex::Real spout = co - sgnm*uo; - amrex::Real spin = cstar - sgnm*ustar; - - //ushock = 0.5_rt*(spin + spout) - amrex::Real ushock = wo*tauo - sgnm*uo; - - if (pstar - po >= 0.0_rt) { - spin = ushock; - spout = ushock; - } - - amrex::Real frac = 0.5_rt*(1.0_rt + (spin + spout)/amrex::max(amrex::max(spout-spin, spin+spout), - riemann_constants::small * raux.cavg)); - - // the transverse velocity states only depend on the - // direction that the contact moves - if (ustar > 0.0_rt) { - qint.ut = ql.ut; - qint.utt = ql.utt; - } else if (ustar < 0.0_rt) { - qint.ut = qr.ut; - qint.utt = qr.utt; - } else { - qint.ut = 0.5_rt * (ql.ut + qr.ut); - qint.utt = 0.5_rt * (ql.utt + qr.utt); - } - - // linearly interpolate between the star and normal state -- this covers the - // case where we are inside the rarefaction fan. - qint.rho = frac*rstar + (1.0_rt - frac)*ro; - qint.un = frac*ustar + (1.0_rt - frac)*uo; - qint.p = frac*pstar + (1.0_rt - frac)*po; - amrex::Real game_int = frac*gamstar + (1.0_rt-frac)*gameo; - - // now handle the cases where instead we are fully in the - // star or fully in the original (l/r) state - if (spout < 0.0_rt) { - qint.rho = ro; - qint.un = uo; - qint.p = po; - game_int = gameo; - } - - if (spin >= 0.0_rt) { - qint.rho = rstar; - qint.un = ustar; - qint.p = pstar; - game_int = gamstar; - } - - qint.p = amrex::max(qint.p, small_pres); - - qint.un = qint.un * raux.bnd_fac; - - // compute the total energy from the internal, p/(gamma - 1), and the kinetic - qint.rhoe = qint.p / (game_int - 1.0_rt); - - // we'll do the passive scalars separately + if (spin >= 0.0_rt) { + qint.rho = rstar; + qint.un = ustar; + qint.p = pstar; + game_int = gamstar; + } -} + qint.p = amrex::max(qint.p, small_pres); + qint.un = qint.un * raux.bnd_fac; -/// -/// The Colella-Glaz-Ferguson Riemann solver for hydrodynamics and -/// radiation hydrodynamics. This is a two shock approximate state -/// Riemann solver. -/// -/// @param bx the box to operate over -/// @param ql the left interface state -/// @param qr the right interface state -/// @param qaux_arr the auxiliary state -/// @param qint the full Godunov state on the interface -/// @param lambda_int the radiation flux limiter on the interface -/// -AMREX_GPU_HOST_DEVICE AMREX_INLINE -void -riemannus(const RiemannState& ql, const RiemannState& qr, const RiemannAux& raux, - RiemannState& qint) { + // compute the total energy from the internal, p/(gamma - 1), and the kinetic + qint.rhoe = qint.p / (game_int - 1.0_rt); - // Colella, Glaz, and Ferguson solver - // - // this is a 2-shock solver that uses a very simple approximation for the - // star state, and carries an auxiliary jump condition for (rho e) to - // deal with a real gas + // we'll do the passive scalars separately + } - // estimate the star state: pstar, ustar - amrex::Real wsmall = small_dens * raux.csmall; + /// + /// The Colella-Glaz-Ferguson Riemann solver for hydrodynamics and + /// radiation hydrodynamics. This is a two shock approximate state + /// Riemann solver. + /// + /// @param bx the box to operate over + /// @param ql the left interface state + /// @param qr the right interface state + /// @param qaux_arr the auxiliary state + /// @param qint the full Godunov state on the interface + /// @param lambda_int the radiation flux limiter on the interface + /// + AMREX_GPU_HOST_DEVICE AMREX_INLINE + void + riemannus(const RiemannState& ql, const RiemannState& qr, const RiemannAux& raux, + RiemannState& qint) { - // this is Castro I: Eq. 33 + // Colella, Glaz, and Ferguson solver + // + // this is a 2-shock solver that uses a very simple approximation for the + // star state, and carries an auxiliary jump condition for (rho e) to + // deal with a real gas - amrex::Real wl = amrex::max(wsmall, std::sqrt(std::abs(ql.gamc * ql.p * ql.rho))); - amrex::Real wr = amrex::max(wsmall, std::sqrt(std::abs(qr.gamc * qr.p * qr.rho))); - amrex::Real wwinv = 1.0_rt/(wl + wr); - amrex::Real pstar = ((wr * ql.p + wl * qr.p) + wl * wr * (ql.un - qr.un)) * wwinv; - amrex::Real ustar = ((wl * ql.un + wr * qr.un) + (ql.p - qr.p)) * wwinv; + // estimate the star state: pstar, ustar - pstar = amrex::max(pstar, small_pres); + amrex::Real wsmall = small_dens * raux.csmall; - // for symmetry preservation, if ustar is really small, then we - // set it to zero + // this is Castro I: Eq. 33 - if (std::abs(ustar) < riemann_constants::smallu * 0.5_rt * (std::abs(ql.un) + std::abs(qr.un))) { - ustar = 0.0_rt; - } + amrex::Real wl = amrex::max(wsmall, std::sqrt(std::abs(ql.gamc * ql.p * ql.rho))); + amrex::Real wr = amrex::max(wsmall, std::sqrt(std::abs(qr.gamc * qr.p * qr.rho))); - // look at the contact to determine which region we are in + amrex::Real wwinv = 1.0_rt/(wl + wr); + amrex::Real pstar = ((wr * ql.p + wl * qr.p) + wl * wr * (ql.un - qr.un)) * wwinv; + amrex::Real ustar = ((wl * ql.un + wr * qr.un) + (ql.p - qr.p)) * wwinv; - // this just determines which of the left or right states is still - // in play. We still need to look at the other wave to determine - // if the star state or this state is on the interface. - amrex::Real sgnm = std::copysign(1.0_rt, ustar); - if (ustar == 0.0_rt) { - sgnm = 0.0_rt; - } + pstar = amrex::max(pstar, small_pres); - amrex::Real fp = 0.5_rt*(1.0_rt + sgnm); - amrex::Real fm = 0.5_rt*(1.0_rt - sgnm); + // for symmetry preservation, if ustar is really small, then we + // set it to zero - amrex::Real ro = fp * ql.rho + fm * qr.rho; - amrex::Real uo = fp * ql.un + fm * qr.un; - amrex::Real po = fp * ql.p + fm * qr.p; - amrex::Real reo = fp * ql.rhoe + fm * qr.rhoe; - amrex::Real gamco = fp * ql.gamc + fm * qr.gamc; -#ifdef RADIATION - for (int g = 0; g < NGROUPS; g++) { - qint.lam[g] = fp * ql.lam[g] + fm * qr.lam[g]; - } + if (std::abs(ustar) < riemann_constants::smallu * 0.5_rt * (std::abs(ql.un) + std::abs(qr.un))) { + ustar = 0.0_rt; + } + + // look at the contact to determine which region we are in - if (ustar == 0) { - // harmonic average + // this just determines which of the left or right states is still + // in play. We still need to look at the other wave to determine + // if the star state or this state is on the interface. + amrex::Real sgnm = std::copysign(1.0_rt, ustar); + if (ustar == 0.0_rt) { + sgnm = 0.0_rt; + } + + amrex::Real fp = 0.5_rt*(1.0_rt + sgnm); + amrex::Real fm = 0.5_rt*(1.0_rt - sgnm); + + amrex::Real ro = fp * ql.rho + fm * qr.rho; + amrex::Real uo = fp * ql.un + fm * qr.un; + amrex::Real po = fp * ql.p + fm * qr.p; + amrex::Real reo = fp * ql.rhoe + fm * qr.rhoe; + amrex::Real gamco = fp * ql.gamc + fm * qr.gamc; + #ifdef RADIATION for (int g = 0; g < NGROUPS; g++) { - qint.lam[g] = 2.0_rt * (ql.lam[g] * qr.lam[g]) / (ql.lam[g] + qr.lam[g] + 1.e-50_rt); + qint.lam[g] = fp * ql.lam[g] + fm * qr.lam[g]; } - } - - amrex::Real po_g = fp * ql.p_g + fm * qr.p_g; - amrex::Real reo_r[NGROUPS]; - amrex::Real po_r[NGROUPS]; - for (int g = 0; g < NGROUPS; g++) { - reo_r[g] = fp * ql.er[g] + fm * qr.er[g]; - po_r[g] = qint.lam[g] * reo_r[g]; - } - amrex::Real reo_g = fp * ql.rhoe_g + fm * qr.rhoe_g; - amrex::Real gamco_g = fp * ql.gamcg + fm * qr.gamcg; -#endif - ro = amrex::max(small_dens, ro); + if (ustar == 0) { + // harmonic average + for (int g = 0; g < NGROUPS; g++) { + qint.lam[g] = 2.0_rt * (ql.lam[g] * qr.lam[g]) / (ql.lam[g] + qr.lam[g] + 1.e-50_rt); + } + } - amrex::Real roinv = 1.0_rt / ro; + amrex::Real po_g = fp * ql.p_g + fm * qr.p_g; + amrex::Real reo_r[NGROUPS]; + amrex::Real po_r[NGROUPS]; + for (int g = 0; g < NGROUPS; g++) { + reo_r[g] = fp * ql.er[g] + fm * qr.er[g]; + po_r[g] = qint.lam[g] * reo_r[g]; + } + amrex::Real reo_g = fp * ql.rhoe_g + fm * qr.rhoe_g; + amrex::Real gamco_g = fp * ql.gamcg + fm * qr.gamcg; + #endif - amrex::Real co = std::sqrt(std::abs(gamco * po * roinv)); - co = amrex::max(raux.csmall, co); - amrex::Real co2inv = 1.0_rt / (co*co); + ro = amrex::max(small_dens, ro); - // we can already deal with the transverse velocities -- they - // only jump across the contact + amrex::Real roinv = 1.0_rt / ro; - qint.ut = fp * ql.ut + fm * qr.ut; - qint.utt = fp * ql.utt + fm * qr.utt; + amrex::Real co = std::sqrt(std::abs(gamco * po * roinv)); + co = amrex::max(raux.csmall, co); + amrex::Real co2inv = 1.0_rt / (co*co); - // compute the rest of the star state + // we can already deal with the transverse velocities -- they + // only jump across the contact - amrex::Real drho = (pstar - po)*co2inv; - amrex::Real rstar = ro + drho; - rstar = amrex::max(small_dens, rstar); + qint.ut = fp * ql.ut + fm * qr.ut; + qint.utt = fp * ql.utt + fm * qr.utt; -#ifdef RADIATION - amrex::Real estar_g = reo_g + drho*(reo_g + po_g)*roinv; + // compute the rest of the star state - amrex::Real co_g = std::sqrt(std::abs(gamco_g*po_g*roinv)); - co_g = amrex::max(raux.csmall, co_g); + amrex::Real drho = (pstar - po)*co2inv; + amrex::Real rstar = ro + drho; + rstar = amrex::max(small_dens, rstar); - amrex::Real pstar_g = po_g + drho*co_g*co_g; - pstar_g = amrex::max(pstar_g, small_pres); + #ifdef RADIATION + amrex::Real estar_g = reo_g + drho*(reo_g + po_g)*roinv; - amrex::Real estar_r[NGROUPS]; - for (int g = 0; g < NGROUPS; g++) { - estar_r[g] = reo_r[g] + drho*(reo_r[g] + po_r[g])*roinv; - } -#else - amrex::Real entho = (reo + po)*roinv*co2inv; - amrex::Real estar = reo + (pstar - po)*entho; -#endif + amrex::Real co_g = std::sqrt(std::abs(gamco_g*po_g*roinv)); + co_g = amrex::max(raux.csmall, co_g); - amrex::Real cstar = std::sqrt(std::abs(gamco*pstar/rstar)); - cstar = amrex::max(cstar, raux.csmall); + amrex::Real pstar_g = po_g + drho*co_g*co_g; + pstar_g = amrex::max(pstar_g, small_pres); - // finish sampling the solution + amrex::Real estar_r[NGROUPS]; + for (int g = 0; g < NGROUPS; g++) { + estar_r[g] = reo_r[g] + drho*(reo_r[g] + po_r[g])*roinv; + } + #else + amrex::Real entho = (reo + po)*roinv*co2inv; + amrex::Real estar = reo + (pstar - po)*entho; + #endif - // look at the remaining wave to determine if the star state or the - // 'o' state above is on the interface + amrex::Real cstar = std::sqrt(std::abs(gamco*pstar/rstar)); + cstar = amrex::max(cstar, raux.csmall); - // the values of u +/- c on either side of the non-contact wave - amrex::Real spout = co - sgnm*uo; - amrex::Real spin = cstar - sgnm*ustar; + // finish sampling the solution - // a simple estimate of the shock speed - amrex::Real ushock = 0.5_rt*(spin + spout); + // look at the remaining wave to determine if the star state or the + // 'o' state above is on the interface - if (pstar-po > 0.0_rt) { - spin = ushock; - spout = ushock; - } + // the values of u +/- c on either side of the non-contact wave + amrex::Real spout = co - sgnm*uo; + amrex::Real spin = cstar - sgnm*ustar; - amrex::Real scr = spout - spin; - if (spout-spin == 0.0_rt) { - scr = riemann_constants::small * raux.cavg; - } + // a simple estimate of the shock speed + amrex::Real ushock = 0.5_rt*(spin + spout); - // interpolate for the case that we are in a rarefaction - amrex::Real frac = (1.0_rt + (spout + spin)/scr)*0.5_rt; - frac = amrex::max(0.0_rt, amrex::min(1.0_rt, frac)); + if (pstar-po > 0.0_rt) { + spin = ushock; + spout = ushock; + } - qint.rho = frac*rstar + (1.0_rt - frac)*ro; - qint.un = frac*ustar + (1.0_rt - frac)*uo; + amrex::Real scr = spout - spin; + if (spout-spin == 0.0_rt) { + scr = riemann_constants::small * raux.cavg; + } -#ifdef RADIATION - amrex::Real pgdnv_t = frac*pstar + (1.0_rt - frac)*po; - amrex::Real pgdnv_g = frac*pstar_g + (1.0_rt - frac)*po_g; - amrex::Real regdnv_g = frac*estar_g + (1.0_rt - frac)*reo_g; - amrex::Real regdnv_r[NGROUPS]; - for (int g = 0; g < NGROUPS; g++) { - regdnv_r[g] = frac*estar_r[g] + (1.0_rt - frac)*reo_r[g]; - } -#else - qint.p = frac*pstar + (1.0_rt - frac)*po; - amrex::Real regdnv = frac*estar + (1.0_rt - frac)*reo; -#endif + // interpolate for the case that we are in a rarefaction + amrex::Real frac = (1.0_rt + (spout + spin)/scr)*0.5_rt; + frac = amrex::max(0.0_rt, amrex::min(1.0_rt, frac)); - // as it stands now, we set things assuming that the rarefaction - // spans the interface. We overwrite that here depending on the - // wave speeds + qint.rho = frac*rstar + (1.0_rt - frac)*ro; + qint.un = frac*ustar + (1.0_rt - frac)*uo; - // look at the speeds on either side of the remaining wave - // to determine which region we are in - if (spout < 0.0_rt) { - // the l or r state is on the interface - qint.rho = ro; - qint.un = uo; -#ifdef RADIATION - pgdnv_t = po; - pgdnv_g = po_g; - regdnv_g = reo_g; + #ifdef RADIATION + amrex::Real pgdnv_t = frac*pstar + (1.0_rt - frac)*po; + amrex::Real pgdnv_g = frac*pstar_g + (1.0_rt - frac)*po_g; + amrex::Real regdnv_g = frac*estar_g + (1.0_rt - frac)*reo_g; + amrex::Real regdnv_r[NGROUPS]; for (int g = 0; g < NGROUPS; g++) { - regdnv_r[g] = reo_r[g]; + regdnv_r[g] = frac*estar_r[g] + (1.0_rt - frac)*reo_r[g]; + } + #else + qint.p = frac*pstar + (1.0_rt - frac)*po; + amrex::Real regdnv = frac*estar + (1.0_rt - frac)*reo; + #endif + + // as it stands now, we set things assuming that the rarefaction + // spans the interface. We overwrite that here depending on the + // wave speeds + + // look at the speeds on either side of the remaining wave + // to determine which region we are in + if (spout < 0.0_rt) { + // the l or r state is on the interface + qint.rho = ro; + qint.un = uo; + #ifdef RADIATION + pgdnv_t = po; + pgdnv_g = po_g; + regdnv_g = reo_g; + for (int g = 0; g < NGROUPS; g++) { + regdnv_r[g] = reo_r[g]; + } + #else + qint.p = po; + regdnv = reo; + #endif } -#else - qint.p = po; - regdnv = reo; -#endif - } - if (spin >= 0.0_rt) { - // the star state is on the interface - qint.rho = rstar; - qint.un = ustar; -#ifdef RADIATION - pgdnv_t = pstar; - pgdnv_g = pstar_g; - regdnv_g = estar_g; + if (spin >= 0.0_rt) { + // the star state is on the interface + qint.rho = rstar; + qint.un = ustar; + #ifdef RADIATION + pgdnv_t = pstar; + pgdnv_g = pstar_g; + regdnv_g = estar_g; + for (int g = 0; g < NGROUPS; g++) { + regdnv_r[g] = estar_r[g]; + } + #else + qint.p = pstar; + regdnv = estar; + #endif + } + + #ifdef RADIATION for (int g = 0; g < NGROUPS; g++) { - regdnv_r[g] = estar_r[g]; + qint.er[g] = amrex::max(regdnv_r[g], 0.0_rt); } -#else - qint.p = pstar; - regdnv = estar; -#endif - } -#ifdef RADIATION - for (int g = 0; g < NGROUPS; g++) { - qint.er[g] = amrex::max(regdnv_r[g], 0.0_rt); - } - - qint.p_g = pgdnv_g; - qint.p = pgdnv_t; - qint.rhoe_g = regdnv_g; - - qint.rhoe = regdnv_g; - for (int g = 0; g < NGROUPS; g++) { - qint.rhoe += regdnv_r[g]; - } - -#else - qint.p = amrex::max(qint.p, small_pres); - qint.rhoe = regdnv; -#endif + qint.p_g = pgdnv_g; + qint.p = pgdnv_t; + qint.rhoe_g = regdnv_g; - // Enforce that fluxes through a symmetry plane or wall are hard zero. - qint.un = qint.un * raux.bnd_fac; + qint.rhoe = regdnv_g; + for (int g = 0; g < NGROUPS; g++) { + qint.rhoe += regdnv_r[g]; + } - // we'll do the passive scalars separately + #else + qint.p = amrex::max(qint.p, small_pres); + qint.rhoe = regdnv; + #endif -} + // Enforce that fluxes through a symmetry plane or wall are hard zero. + qint.un = qint.un * raux.bnd_fac; + // we'll do the passive scalars separately + + } +} #endif diff --git a/Source/hydro/riemann_solvers.H b/Source/hydro/riemann_solvers.H index ded2d5a169..b3c3509dbf 100644 --- a/Source/hydro/riemann_solvers.H +++ b/Source/hydro/riemann_solvers.H @@ -321,15 +321,13 @@ riemann_state(const int i, const int j, const int k, const int idir, if (riemann_solver == 0) { // Colella, Glaz, & Ferguson solver - riemannus(ql, qr, raux, - qint); + TwoShock::riemannus(ql, qr, raux, qint); } else if (riemann_solver == 1) { // Colella & Glaz solver #ifndef RADIATION - riemanncg(ql, qr, raux, - qint); + TwoShock::riemanncg(ql, qr, raux, qint); #endif #ifndef AMREX_USE_GPU