From b60f575838b59d939e967b254c1dcd419e187fdb Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Tue, 20 Aug 2024 13:29:23 -0400 Subject: [PATCH] implement full well-balancing in PPM --- Source/driver/_cpp_parameters | 4 ++ Source/hydro/ppm.H | 77 +++++++++++++---------- Source/hydro/trace_ppm.cpp | 111 ++++++++++++++++++++-------------- 3 files changed, 114 insertions(+), 78 deletions(-) diff --git a/Source/driver/_cpp_parameters b/Source/driver/_cpp_parameters index 44f78e6ebf..8a2eb146d2 100644 --- a/Source/driver/_cpp_parameters +++ b/Source/driver/_cpp_parameters @@ -161,6 +161,10 @@ dual_energy_eta2 Real 1.0e-4 # does well with HSE use_pslope bool 0 +# for PPM, do we only use the perturbational pressure in the charactersitic +# tracing? This is more indepth than the simple `use_pslope` approach. +ppm_well_balanced bool 0 + # if we are using pslope, below what density to we turn off the well-balanced # reconstruction? pslope_cutoff_density Real -1.e20 diff --git a/Source/hydro/ppm.H b/Source/hydro/ppm.H index bf854da880..fea87bc29e 100644 --- a/Source/hydro/ppm.H +++ b/Source/hydro/ppm.H @@ -1,16 +1,18 @@ #include + +#include #include #ifndef PPM_H #define PPM_H -using namespace amrex; +using namespace amrex::literals; using namespace reconstruction; AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE int -check_trace_source(Array4 const& srcQ, const int idir, +check_trace_source(amrex::Array4 const& srcQ, const int idir, const int i, const int j, const int k, const int ncomp) { int do_trace = 0; @@ -53,29 +55,29 @@ check_trace_source(Array4 const& srcQ, const int idir, /// AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void -ppm_reconstruct(const Real* s, - const Real flatn, Real& sm, Real& sp) { +ppm_reconstruct(const amrex::Real* s, + const amrex::Real flatn, amrex::Real& sm, amrex::Real& sp) { if (ppm_do_limiting) { // Compute van Leer slopes - Real dsl = 2.0_rt * (s[im1] - s[im2]); - Real dsr = 2.0_rt * (s[i0] - s[im1]); + amrex::Real dsl = 2.0_rt * (s[im1] - s[im2]); + amrex::Real dsr = 2.0_rt * (s[i0] - s[im1]); - Real dsvl_l = 0.0_rt; + amrex::Real dsvl_l = 0.0_rt; if (dsl*dsr > 0.0_rt) { - Real dsc = 0.5_rt * (s[i0] - s[im2]); + amrex::Real dsc = 0.5_rt * (s[i0] - s[im2]); dsvl_l = std::copysign(1.0_rt, dsc) * amrex::min(std::abs(dsc), std::abs(dsl), std::abs(dsr)); } dsl = 2.0_rt * (s[i0] - s[im1]); dsr = 2.0_rt * (s[ip1] - s[i0]); - Real dsvl_r = 0.0_rt; + amrex::Real dsvl_r = 0.0_rt; if (dsl*dsr > 0.0_rt) { - Real dsc = 0.5_rt * (s[ip1] - s[im1]); + amrex::Real dsc = 0.5_rt * (s[ip1] - s[im1]); dsvl_r = std::copysign(1.0_rt, dsc) * amrex::min(std::abs(dsc), std::abs(dsl),std::abs(dsr)); } @@ -94,8 +96,8 @@ ppm_reconstruct(const Real* s, dsvl_l = 0.0_rt; if (dsl*dsr > 0.0_rt) { - Real dsc = 0.5_rt * (s[ip1] - s[im1]); - dsvl_l = std::copysign(1.0_rt, dsc) * amrex::min(std::abs(dsc), std::abs(dsl), std::abs(dsr)); + amrex::Real dsc = 0.5_rt * (s[ip1] - s[im1]); + dsvl_l = std::copysign(1.0_rt, dsc) * amrex::min(std::abs(dsc), std::abs(dsl), std::abs(dsr)); } dsl = 2.0_rt * (s[ip1] - s[i0]); @@ -103,8 +105,8 @@ ppm_reconstruct(const Real* s, dsvl_r = 0.0_rt; if (dsl*dsr > 0.0_rt) { - Real dsc = 0.5_rt * (s[ip2] - s[i0]); - dsvl_r = std::copysign(1.0_rt, dsc) * amrex::min(std::abs(dsc), std::abs(dsl), std::abs(dsr)); + amrex::Real dsc = 0.5_rt * (s[ip2] - s[i0]); + dsvl_r = std::copysign(1.0_rt, dsc) * amrex::min(std::abs(dsc), std::abs(dsl), std::abs(dsr)); } // Interpolate s to edges @@ -153,7 +155,7 @@ ppm_reconstruct(const Real* s, /// the gravitational force by only limiting on the wave-generating pressure. /// This uses the standard PPM limiters described in Colella & Woodward (1984) /// -/// @param rho Real[nslp] giving the density in zones i-2, i-1, i, i+1, i+2 +/// @param rho Real[nslp] giving the density in zones i-2, i-1, i, i+1, i+2 /// @param p Real[nslp] giving the pressure in zones i-2, i-1, i, i+1, i+2 /// @param src Real[nslp] the source in the velocity equation (e.g. g) in zones /// i-2, i-1, i, i+2, i+2 @@ -164,23 +166,25 @@ ppm_reconstruct(const Real* s, /// @param sm The value of the parabola on the left edge of the zone /// @param sp The value of the parabola on the right edge of the zone /// +/// @return a bool indicating whether HSE correction was done +/// AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -void -ppm_reconstruct_pslope(const Real* rho, const Real* p, const Real* src, const Real flatn, - const Real dx, - Real& sm, Real& sp) { +bool +ppm_reconstruct_pslope(const amrex::Real* rho, const amrex::Real* p, const amrex::Real* src, const amrex::Real flatn, + const amrex::Real dx, + amrex::Real& sm, amrex::Real& sp) { - Real tp[nslp]; + amrex::Real tp[nslp]; // compute the hydrostatic pressure in each zone center starting with i0 - Real p0_hse = p[i0]; + amrex::Real p0_hse = p[i0]; - Real pp1_hse = p0_hse + 0.25_rt*dx * (rho[i0] + rho[ip1]) * (src[i0] + src[ip1]); - Real pp2_hse = pp1_hse + 0.25_rt*dx * (rho[ip1] + rho[ip2]) * (src[ip1] + src[ip2]); + amrex::Real pp1_hse = p0_hse + 0.25_rt*dx * (rho[i0] + rho[ip1]) * (src[i0] + src[ip1]); + amrex::Real pp2_hse = pp1_hse + 0.25_rt*dx * (rho[ip1] + rho[ip2]) * (src[ip1] + src[ip2]); - Real pm1_hse = p0_hse - 0.25_rt*dx * (rho[i0] + rho[im1]) * (src[i0] + src[im1]); - Real pm2_hse = pm1_hse - 0.25_rt*dx * (rho[im1] + rho[im2]) * (src[im1] + src[im2]); + amrex::Real pm1_hse = p0_hse - 0.25_rt*dx * (rho[i0] + rho[im1]) * (src[i0] + src[im1]); + amrex::Real pm2_hse = pm1_hse - 0.25_rt*dx * (rho[im1] + rho[im2]) * (src[im1] + src[im2]); // this only makes sense if the pressures are positive bool ptest = p0_hse < 0.0 || @@ -190,7 +194,9 @@ ppm_reconstruct_pslope(const Real* rho, const Real* p, const Real* src, const Re pm2_hse < 0.0; - if (!ptest && rho[i0] >= pslope_cutoff_density) { + bool do_hse = !ptest && rho[i0] >= pslope_cutoff_density; + + if (do_hse) { // subtract off the hydrostatic pressure @@ -218,12 +224,16 @@ ppm_reconstruct_pslope(const Real* rho, const Real* p, const Real* src, const Re // now correct sm and sp to be back to the full pressure by // adding in the hydrostatic pressure at the interface + // if we are doing the full well-balanced method, then we + // don't do this until after the characteristic tracing - if (!ptest && rho[i0] >= pslope_cutoff_density) { + if (do_hse && !castro::ppm_well_balanced) { sp += p[i0] + 0.5 * dx * rho[i0] * src[i0]; sm += p[i0] - 0.5 * dx * rho[i0] * src[i0]; } + return do_hse; + } @@ -244,14 +254,15 @@ ppm_reconstruct_pslope(const Real* rho, const Real* p, const Real* src, const Re /// AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void -ppm_int_profile(const Real sm, const Real sp, const Real sc, - const Real u, const Real c, const Real dtdx, - Real* Ip, Real* Im) { +ppm_int_profile(const amrex::Real sm, const amrex::Real sp, const amrex::Real sc, + const amrex::Real u, const amrex::Real c, const amrex::Real dtdx, + amrex::Real* Ip, amrex::Real* Im) { - // Integrate the parabolic profile to the edge of the cell. + // Integrate the parabolic profile to the edge of the cell. - // compute x-component of Ip and Im - Real s6 = 6.0_rt * sc - 3.0_rt * (sm + sp); + // compute x-component of Ip and Im + + Real s6 = 6.0_rt * sc - 3.0_rt * (sm + sp); // Ip/m is the integral under the parabola for the extent // that a wave can travel over a timestep diff --git a/Source/hydro/trace_ppm.cpp b/Source/hydro/trace_ppm.cpp index 2595b9fe65..2c9bd52027 100644 --- a/Source/hydro/trace_ppm.cpp +++ b/Source/hydro/trace_ppm.cpp @@ -217,6 +217,14 @@ Castro::trace_ppm(const Box& bx, load_stencil(q_arr, idir, i, j, k, QPRES, s); + bool in_hse{}; + + // HSE pressure on interfaces -- needed if we are dealing with + // perturbation pressure as the parabolic reconstruction + + amrex::Real p_m_hse{}; + amrex::Real p_p_hse{}; + if (use_pslope) { Real trho[nslp]; Real src[nslp]; @@ -224,12 +232,22 @@ Castro::trace_ppm(const Box& bx, load_stencil(q_arr, idir, i, j, k, QRHO, trho); load_stencil(srcQ, idir, i, j, k, QUN, src); - ppm_reconstruct_pslope(trho, s, src, flat, dx[idir], sm, sp); + in_hse = ppm_reconstruct_pslope(trho, s, src, flat, dx[idir], sm, sp); + + if (in_hse && castro::ppm_well_balanced) { + // we are working with the perturbational pressure + ppm_int_profile(sm, sp, 0.0_rt, un, cc, dtdx, Ip_p, Im_p); + p_m_hse = s[i0] - 0.5_rt * dx[idir] * trho[i0] * src[i0]; + p_p_hse = s[i0] + 0.5_rt * dx[idir] * trho[i0] * src[i0]; + + } else { + ppm_int_profile(sm, sp, s[i0], un, cc, dtdx, Ip_p, Im_p); + } } else { ppm_reconstruct(s, flat, sm, sp); + ppm_int_profile(sm, sp, s[i0], un, cc, dtdx, Ip_p, Im_p); } - ppm_int_profile(sm, sp, s[i0], un, cc, dtdx, Ip_p, Im_p); // reconstruct rho e @@ -420,6 +438,10 @@ Castro::trace_ppm(const Box& bx, } + // for well-balanced, the velocity sources should not be added + + amrex::Real fac = (castro::ppm_well_balanced && in_hse) ? 0.0_rt : 1.0_rt; + // plus state on face i if ((idir == 0 && i >= vlo[0]) || @@ -437,23 +459,23 @@ Castro::trace_ppm(const Box& bx, // the dt * g term will be the only non-zero contribution). // We ignore the effect of the source term for gamma. Real rho_ref = Im_rho[0] + hdt * Im_src_rho[0]; - Real un_ref = Im_un_0 + hdt * Im_src_un_0; + Real un_ref = Im_un_0 + fac * hdt * Im_src_un_0; Real p_ref = Im_p[0] + hdt * Im_src_p[0]; Real rhoe_g_ref = Im_rhoe[0] + hdt * Im_src_rhoe[0]; Real gam_g_ref = Im_gc_0; - rho_ref = amrex::max(rho_ref, lsmall_dens); + rho_ref = std::max(rho_ref, lsmall_dens); Real rho_ref_inv = 1.0_rt/rho_ref; - p_ref = amrex::max(p_ref, lsmall_pres); + p_ref = std::max(p_ref, lsmall_pres); // For tracing - Real csq_ref = gam_g_ref*p_ref*rho_ref_inv; + Real csq_ref = gam_g_ref * (p_m_hse + p_ref) * rho_ref_inv; Real cc_ref = std::sqrt(csq_ref); Real cc_ref_inv = 1.0_rt/cc_ref; - Real h_g_ref = (p_ref + rhoe_g_ref)*rho_ref_inv; + Real h_g_ref = ((p_m_hse + p_ref) + rhoe_g_ref) * rho_ref_inv; // *m are the jumps carried by un-c // *p are the jumps carried by un+c @@ -463,15 +485,16 @@ Castro::trace_ppm(const Box& bx, // we also add the sources here so they participate in the tracing - Real dum = un_ref - Im_un_0 - hdt*Im_src_un_0; - Real dptotm = p_ref - Im_p[0] - hdt*Im_src_p[0]; - Real drho = rho_ref - Im_rho[1] - hdt*Im_src_rho[1]; - Real dptot = p_ref - Im_p[1] - hdt*Im_src_p[1]; - Real drhoe_g = rhoe_g_ref - Im_rhoe[1] - hdt*Im_src_rhoe[1]; + Real dum = un_ref - Im_un_0 - fac * hdt * Im_src_un_0; + Real dptotm = p_ref - Im_p[0] - hdt * Im_src_p[0]; + + Real drho = rho_ref - Im_rho[1] - hdt * Im_src_rho[1]; + Real dptot = p_ref - Im_p[1] - hdt * Im_src_p[1]; + Real drhoe_g = rhoe_g_ref - Im_rhoe[1] - hdt * Im_src_rhoe[1]; - Real dup = un_ref - Im_un_2 - hdt*Im_src_un_2; - Real dptotp = p_ref - Im_p[2] - hdt*Im_src_p[2]; + Real dup = un_ref - Im_un_2 - fac * hdt * Im_src_un_2; + Real dptotp = p_ref - Im_p[2] - hdt * Im_src_p[2]; // {rho, u, p, (rho e)} eigensystem @@ -492,11 +515,11 @@ Castro::trace_ppm(const Box& bx, // The final interface states are just // q_s = q_ref - sum(l . dq) r // note that the a{mpz}right as defined above have the minus already - qp(i,j,k,QRHO) = amrex::max(lsmall_dens, rho_ref + alphap + alpham + alpha0r); + qp(i,j,k,QRHO) = std::max(lsmall_dens, rho_ref + alphap + alpham + alpha0r); qp(i,j,k,QUN) = un_ref + (alphap - alpham)*cc_ref*rho_ref_inv; - qp(i,j,k,QREINT) = amrex::max(castro::small_dens * castro::small_ener, + qp(i,j,k,QREINT) = std::max(castro::small_dens * castro::small_ener, rhoe_g_ref + (alphap + alpham)*h_g_ref + alpha0e_g); - qp(i,j,k,QPRES) = amrex::max(lsmall_pres, p_ref + (alphap + alpham)*csq_ref); + qp(i,j,k,QPRES) = std::max(lsmall_pres, p_m_hse + p_ref + (alphap + alpham)*csq_ref); // Transverse velocities -- there's no projection here, so we // don't need a reference state. We only care about the state @@ -505,8 +528,8 @@ Castro::trace_ppm(const Box& bx, // Recall that I already takes the limit of the parabola // in the event that the wave is not moving toward the // interface - qp(i,j,k,QUT) = Im_ut_1 + hdt*Im_src_ut_1; - qp(i,j,k,QUTT) = Im_utt_1 + hdt*Im_src_utt_1; + qp(i,j,k,QUT) = Im_ut_1 + hdt * Im_src_ut_1; + qp(i,j,k,QUTT) = Im_utt_1 + hdt * Im_src_utt_1; } @@ -519,35 +542,35 @@ Castro::trace_ppm(const Box& bx, // Set the reference state // This will be the fastest moving state to the right Real rho_ref = Ip_rho[2] + hdt * Ip_src_rho[2]; - Real un_ref = Ip_un_2 + hdt * Ip_src_un_2; + Real un_ref = Ip_un_2 + fac * hdt * Ip_src_un_2; Real p_ref = Ip_p[2] + hdt * Ip_src_p[2]; Real rhoe_g_ref = Ip_rhoe[2] + hdt * Ip_src_rhoe[2]; Real gam_g_ref = Ip_gc_2; - rho_ref = amrex::max(rho_ref, lsmall_dens); + rho_ref = std::max(rho_ref, lsmall_dens); Real rho_ref_inv = 1.0_rt/rho_ref; - p_ref = amrex::max(p_ref, lsmall_pres); + p_ref = std::max(p_ref, lsmall_pres); // For tracing - Real csq_ref = gam_g_ref*p_ref*rho_ref_inv; + Real csq_ref = gam_g_ref * (p_p_hse + p_ref) * rho_ref_inv; Real cc_ref = std::sqrt(csq_ref); Real cc_ref_inv = 1.0_rt/cc_ref; - Real h_g_ref = (p_ref + rhoe_g_ref)*rho_ref_inv; + Real h_g_ref = ((p_p_hse + p_ref) + rhoe_g_ref) * rho_ref_inv; // *m are the jumps carried by u-c // *p are the jumps carried by u+c - Real dum = un_ref - Ip_un_0 - hdt*Ip_src_un_0; - Real dptotm = p_ref - Ip_p[0] - hdt*Ip_src_p[0]; + Real dum = un_ref - Ip_un_0 - fac * hdt * Ip_src_un_0; + Real dptotm = p_ref - Ip_p[0] - hdt * Ip_src_p[0]; - Real drho = rho_ref - Ip_rho[1] - hdt*Ip_src_rho[1]; - Real dptot = p_ref - Ip_p[1] - hdt*Ip_src_p[1]; - Real drhoe_g = rhoe_g_ref - Ip_rhoe[1] - hdt*Ip_src_rhoe[1]; + Real drho = rho_ref - Ip_rho[1] - hdt * Ip_src_rho[1]; + Real dptot = p_ref - Ip_p[1] - hdt * Ip_src_p[1]; + Real drhoe_g = rhoe_g_ref - Ip_rhoe[1] - hdt * Ip_src_rhoe[1]; - Real dup = un_ref - Ip_un_2 - hdt*Ip_src_un_2; - Real dptotp = p_ref - Ip_p[2] - hdt*Ip_src_p[2]; + Real dup = un_ref - Ip_un_2 - fac * hdt * Ip_src_un_2; + Real dptotp = p_ref - Ip_p[2] - hdt * Ip_src_p[2]; // {rho, u, p, (rho e)} eigensystem @@ -569,33 +592,33 @@ Castro::trace_ppm(const Box& bx, // q_s = q_ref - sum (l . dq) r // note that the a{mpz}left as defined above have the minus already if (idir == 0) { - qm(i+1,j,k,QRHO) = amrex::max(lsmall_dens, rho_ref + alphap + alpham + alpha0r); + qm(i+1,j,k,QRHO) = std::max(lsmall_dens, rho_ref + alphap + alpham + alpha0r); qm(i+1,j,k,QUN) = un_ref + (alphap - alpham)*cc_ref*rho_ref_inv; - qm(i+1,j,k,QREINT) = amrex::max(castro::small_dens * castro::small_ener, + qm(i+1,j,k,QREINT) = std::max(castro::small_dens * castro::small_ener, rhoe_g_ref + (alphap + alpham)*h_g_ref + alpha0e_g); - qm(i+1,j,k,QPRES) = amrex::max(lsmall_pres, p_ref + (alphap + alpham)*csq_ref); + qm(i+1,j,k,QPRES) = std::max(lsmall_pres, p_p_hse + p_ref + (alphap + alpham)*csq_ref); // transverse velocities qm(i+1,j,k,QUT) = Ip_ut_1 + hdt*Ip_src_ut_1; qm(i+1,j,k,QUTT) = Ip_utt_1 + hdt*Ip_src_utt_1; } else if (idir == 1) { - qm(i,j+1,k,QRHO) = amrex::max(lsmall_dens, rho_ref + alphap + alpham + alpha0r); + qm(i,j+1,k,QRHO) = std::max(lsmall_dens, rho_ref + alphap + alpham + alpha0r); qm(i,j+1,k,QUN) = un_ref + (alphap - alpham)*cc_ref*rho_ref_inv; - qm(i,j+1,k,QREINT) = amrex::max(castro::small_dens * castro::small_ener, - rhoe_g_ref + (alphap + alpham)*h_g_ref + alpha0e_g); - qm(i,j+1,k,QPRES) = amrex::max(lsmall_pres, p_ref + (alphap + alpham)*csq_ref); + qm(i,j+1,k,QREINT) = std::max(castro::small_dens * castro::small_ener, + rhoe_g_ref + (alphap + alpham)*h_g_ref + alpha0e_g); + qm(i,j+1,k,QPRES) = std::max(lsmall_pres, p_p_hse + p_ref + (alphap + alpham)*csq_ref); // transverse velocities qm(i,j+1,k,QUT) = Ip_ut_1 + hdt*Ip_src_ut_1; qm(i,j+1,k,QUTT) = Ip_utt_1 + hdt*Ip_src_utt_1; } else if (idir == 2) { - qm(i,j,k+1,QRHO) = amrex::max(lsmall_dens, rho_ref + alphap + alpham + alpha0r); + qm(i,j,k+1,QRHO) = std::max(lsmall_dens, rho_ref + alphap + alpham + alpha0r); qm(i,j,k+1,QUN) = un_ref + (alphap - alpham)*cc_ref*rho_ref_inv; - qm(i,j,k+1,QREINT) = amrex::max(castro::small_dens * castro::small_ener, + qm(i,j,k+1,QREINT) = std::max(castro::small_dens * castro::small_ener, rhoe_g_ref + (alphap + alpham)*h_g_ref + alpha0e_g); - qm(i,j,k+1,QPRES) = amrex::max(lsmall_pres, p_ref + (alphap + alpham)*csq_ref); + qm(i,j,k+1,QPRES) = std::max(lsmall_pres, p_p_hse + p_ref + (alphap + alpham)*csq_ref); // transverse velocities qm(i,j,k+1,QUT) = Ip_ut_1 + hdt*Ip_src_ut_1; @@ -619,14 +642,14 @@ Castro::trace_ppm(const Box& bx, if (i <= vhi[0]) { qm(i+1,j,k,QRHO) = qm(i+1,j,k,QRHO) + sourcr; - qm(i+1,j,k,QRHO) = amrex::max(qm(i+1,j,k,QRHO), lsmall_dens); + qm(i+1,j,k,QRHO) = std::max(qm(i+1,j,k,QRHO), lsmall_dens); qm(i+1,j,k,QPRES) = qm(i+1,j,k,QPRES) + sourcp; qm(i+1,j,k,QREINT) = qm(i+1,j,k,QREINT) + source; } if (i >= vlo[0]) { qp(i,j,k,QRHO) = qp(i,j,k,QRHO) + sourcr; - qp(i,j,k,QRHO) = amrex::max(qp(i,j,k,QRHO), lsmall_dens); + qp(i,j,k,QRHO) = std::max(qp(i,j,k,QRHO), lsmall_dens); qp(i,j,k,QPRES) = qp(i,j,k,QPRES) + sourcp; qp(i,j,k,QREINT) = qp(i,j,k,QREINT) + source; } @@ -635,5 +658,3 @@ Castro::trace_ppm(const Box& bx, }); } - -