diff --git a/Docs/source/hse.rst b/Docs/source/hse.rst index a58acee2ce..de34b10b8d 100644 --- a/Docs/source/hse.rst +++ b/Docs/source/hse.rst @@ -74,6 +74,12 @@ then done on the parabola, again working with :math:`\tilde{p}`. Finally, the parabola values are updated to include the hydrostatic pressure. +.. index:: castro.ppm_well_balanced + +We can do better with PPM, and only use the perturbational pressure, +$\tilde{p}$, in the characteristic tracing and then add back the +hydrostatic pressure to the interface afterwards. This is done via +``castro.ppm_well_balanced=1``. Fully fourth-order method ------------------------- @@ -165,5 +171,3 @@ test the different HSE approaches. This sets up a 1-d X-ray burst atmosphere (based on the ``flame_wave`` setup). Richardson extrapolation can be used to measure the convergence rate (or just look at how the peak velocity changes). - - diff --git a/Exec/gravity_tests/hse_convergence/README.md b/Exec/gravity_tests/hse_convergence/README.md index fab88d7bd7..d8d909421c 100644 --- a/Exec/gravity_tests/hse_convergence/README.md +++ b/Exec/gravity_tests/hse_convergence/README.md @@ -31,10 +31,13 @@ To run this problem, use one of the convergence scripts: * convergence_sdc.sh : - this uses the TRUE_SDC integration, first with SDC-2 + PLM and - reflecting BCs, the SDC-2 + PPM and reflecting BCs, then the same - but HSE BCs, and finally SDC-4 + reflect + this uses the `TRUE_SDC` integration, with the following variations: + 1. SDC-2 + PLM and reflecting BCs + 2. SDC-2 + PPM and reflecting BCs + 3. SDC-2 + PLM with HSE BCs + 4. SDC-2 + PPM with HSE BCs + 5. SDC-4 + reflect These tests show that the PLM + reflect (which uses the - well-balanced use_pslope) and the SDC-4 + reflect give the lowest + well-balanced `use_pslope`) and the SDC-4 + reflect give the lowest errors and expected (or better) convergence. diff --git a/Exec/gravity_tests/hse_convergence/convergence_ppm.sh b/Exec/gravity_tests/hse_convergence/convergence_ppm.sh index 27b4bb86ef..3dedc8b11c 100755 --- a/Exec/gravity_tests/hse_convergence/convergence_ppm.sh +++ b/Exec/gravity_tests/hse_convergence/convergence_ppm.sh @@ -102,3 +102,30 @@ ${EXEC} inputs.ppm.512 ${RUNPARAMS} >& 512.out pfile=`ls -t | grep -i hse_512_plt | head -1` fextrema.gnu.ex -v magvel ${pfile} | grep -i magvel >> ${ofile} + +## ppm + reflect + well-balanced + +ofile=ppm-reflect-wellbalanced.converge.out + +RUNPARAMS=" +castro.lo_bc=3 +castro.hi_bc=3 +castro.ppm_well_balanced=1 +""" + +${EXEC} inputs.ppm.64 ${RUNPARAMS} >& 64.out +pfile=`ls -t | grep -i hse_64_plt | head -1` +fextrema.gnu.ex -v magvel ${pfile} | grep -i magvel > ${ofile} + +${EXEC} inputs.ppm.128 ${RUNPARAMS} >& 128.out +pfile=`ls -t | grep -i hse_128_plt | head -1` +fextrema.gnu.ex -v magvel ${pfile} | grep -i magvel >> ${ofile} + +${EXEC} inputs.ppm.256 ${RUNPARAMS} >& 256.out +pfile=`ls -t | grep -i hse_256_plt | head -1` +fextrema.gnu.ex -v magvel ${pfile} | grep -i magvel >> ${ofile} + +${EXEC} inputs.ppm.512 ${RUNPARAMS} >& 512.out +pfile=`ls -t | grep -i hse_512_plt | head -1` +fextrema.gnu.ex -v magvel ${pfile} | grep -i magvel >> ${ofile} + diff --git a/Exec/gravity_tests/hse_convergence/convergence_sdc.sh b/Exec/gravity_tests/hse_convergence/convergence_sdc.sh index 038109dea1..407c754759 100755 --- a/Exec/gravity_tests/hse_convergence/convergence_sdc.sh +++ b/Exec/gravity_tests/hse_convergence/convergence_sdc.sh @@ -67,14 +67,15 @@ pfile=`ls -t | grep -i hse_512_plt | head -1` fextrema.gnu.ex -v magvel ${pfile} | grep -i magvel >> ${ofile} -## sdc-2 + ppm +## sdc-2 + HSE -ofile=sdc2-ppm.converge.out +ofile=sdc2.converge.out RUNPARAMS=" castro.time_integration_method=2 castro.sdc_order=2 -castro.ppm_type=1 +castro.ppm_type=0 +castro.use_pslope=1 castro.limit_fourth_order=1 castro.use_reconstructed_gamma1=1 """ @@ -96,15 +97,14 @@ pfile=`ls -t | grep -i hse_512_plt | head -1` fextrema.gnu.ex -v magvel ${pfile} | grep -i magvel >> ${ofile} -## sdc-2 + HSE +## sdc-2 + ppm -ofile=sdc2.converge.out +ofile=sdc2-ppm.converge.out RUNPARAMS=" castro.time_integration_method=2 castro.sdc_order=2 -castro.ppm_type=0 -castro.use_pslope=1 +castro.ppm_type=1 castro.limit_fourth_order=1 castro.use_reconstructed_gamma1=1 """ @@ -126,7 +126,6 @@ pfile=`ls -t | grep -i hse_512_plt | head -1` fextrema.gnu.ex -v magvel ${pfile} | grep -i magvel >> ${ofile} - ## sdc-4 + reflect ofile=sdc4-reflect.converge.out diff --git a/Exec/science/flame_wave/ci-benchmarks/grid_diag.out b/Exec/science/flame_wave/ci-benchmarks/grid_diag.out index 3bfc78ca85..5cea0e9518 100644 --- a/Exec/science/flame_wave/ci-benchmarks/grid_diag.out +++ b/Exec/science/flame_wave/ci-benchmarks/grid_diag.out @@ -1,5 +1,5 @@ # COLUMN 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 # TIMESTEP TIME MASS XMOM YMOM ZMOM ANG. MOM. X ANG. MOM. Y ANG. MOM. Z KIN. ENERGY INT. ENERGY GAS ENERGY GRAV. ENERGY TOTAL ENERGY CENTER OF MASS X-LOC CENTER OF MASS Y-LOC CENTER OF MASS Z-LOC CENTER OF MASS X-VEL CENTER OF MASS Y-VEL CENTER OF MASS Z-VEL MAXIMUM TEMPERATURE MAXIMUM DENSITY MAXIMUM T_S / T_E 0 0.0000000000000000e+00 1.4700685690736437e+20 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 2.9163021935541997e+37 2.9163021935541997e+37 0.0000000000000000e+00 2.9163021935541997e+37 2.7306641645125415e+04 6.8087525965685256e+02 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 1.3928678114307070e+09 3.0715027784567930e+07 0.0000000000000000e+00 - 1 5.9806047036494849e-08 1.4700700110495903e+20 3.3564221217849212e+23 7.2876252590226133e+23 1.2612378556916341e+20 3.1318391966955516e+23 -2.7141686347841742e+24 1.9066097493127676e+28 7.8605324800324878e+29 2.9163081863790706e+37 2.9163082649843738e+37 0.0000000000000000e+00 2.9163082649843738e+37 2.7306641717925955e+04 6.8087482283157317e+02 0.0000000000000000e+00 2.2831716153358752e+03 4.9573321027203638e+03 8.5794407491595881e-01 1.3929561431310942e+09 3.0715326643212341e+07 3.3873763041849706e-04 - 2 1.2559269877663918e-07 1.4700717750115462e+20 7.0484895445039834e+23 1.5304335700823335e+24 5.5619061962888426e+20 1.3810900989697918e+24 -1.1969287484132617e+25 4.0039912979259806e+28 3.0853275723280981e+30 2.9163145904135894e+37 2.9163148989463411e+37 0.0000000000000000e+00 2.9163148989463411e+37 2.7306641954244522e+04 6.8087461718998338e+02 0.0000000000000000e+00 4.7946567401095926e+03 1.0410604407871944e+04 3.7834249257966728e+00 1.3933944939675827e+09 3.0715562157563787e+07 3.3857586854903492e-04 + 1 5.9806047036494849e-08 1.4700700110495903e+20 3.3564221215590610e+23 7.2876252589700228e+23 1.2612378556594561e+20 3.1318391966786556e+23 -2.7141686347838510e+24 1.9066097493139366e+28 7.8605324800584047e+29 2.9163081863790706e+37 2.9163082649843738e+37 0.0000000000000000e+00 2.9163082649843738e+37 2.7306641717925955e+04 6.8087482283158556e+02 0.0000000000000000e+00 2.2831716151822361e+03 4.9573321026845897e+03 8.5794407489407010e-01 1.3929529135508094e+09 3.0715326539046615e+07 3.3873763041849706e-04 + 2 1.2559269877663918e-07 1.4700717750115462e+20 7.0484895413073359e+23 1.5304335700221490e+24 5.5619061951241526e+20 1.3810900989051759e+24 -1.1969287484120623e+25 4.0039912979436722e+28 3.0853275724059276e+30 2.9163145904135894e+37 2.9163148989463416e+37 0.0000000000000000e+00 2.9163148989463416e+37 2.7306641954244522e+04 6.8087461719003727e+02 0.0000000000000000e+00 4.7946567379351090e+03 1.0410604407462544e+04 3.7834249250044056e+00 1.3933796339641171e+09 3.0715561813912578e+07 3.3857586854903492e-04 diff --git a/Exec/science/flame_wave/ci-benchmarks/job_info_params.txt b/Exec/science/flame_wave/ci-benchmarks/job_info_params.txt index ff734d8ae0..f3b20994e7 100644 --- a/Exec/science/flame_wave/ci-benchmarks/job_info_params.txt +++ b/Exec/science/flame_wave/ci-benchmarks/job_info_params.txt @@ -64,6 +64,7 @@ castro.dual_energy_eta1 = 1 castro.dual_energy_eta2 = 0.0001 castro.use_pslope = 0 + castro.ppm_well_balanced = 0 castro.pslope_cutoff_density = -1e+20 castro.limit_fluxes_on_small_dens = 0 castro.speed_limit = 0 diff --git a/Exec/science/flame_wave/ci-benchmarks/species_diag.out b/Exec/science/flame_wave/ci-benchmarks/species_diag.out index f3c2dcfe7a..f2e2da45d1 100644 --- a/Exec/science/flame_wave/ci-benchmarks/species_diag.out +++ b/Exec/science/flame_wave/ci-benchmarks/species_diag.out @@ -1,5 +1,5 @@ # COLUMN 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 # TIMESTEP TIME Mass He4 Mass C12 Mass O16 Mass Ne20 Mass Mg24 Mass Si28 Mass S32 Mass Ar36 Mass Ca40 Mass Ti44 Mass Cr48 Mass Fe52 Mass Ni56 0 0.0000000000000000e+00 2.8142378112400895e-15 7.3932235330395157e-24 7.3932235330395157e-24 7.3932235330395157e-24 7.3932235330395157e-24 7.3932235330395157e-24 7.3932235330395157e-24 7.3932235330395157e-24 7.3932235330395157e-24 7.3932235330395157e-24 7.3932235330395157e-24 7.3932235330395157e-24 7.1117997526547418e-14 - 1 5.9806047036494849e-08 2.8142323856670446e-15 5.4329394398418374e-21 7.4194443455948811e-24 7.3934514829832856e-24 7.3935541949306949e-24 7.3933935989353088e-24 7.3932321468740865e-24 7.3932312019180345e-24 7.3932308328287333e-24 7.3932307869392985e-24 7.3932307850741519e-24 7.3932307849933044e-24 7.1118070045957103e-14 - 2 1.2559269877663918e-07 2.8142264165381050e-15 1.1401993450709496e-20 7.4929528148974785e-24 7.3944694290388374e-24 7.3939305042586412e-24 7.3935838894001731e-24 7.3932425308122493e-24 7.3932405326301655e-24 7.3932397568503098e-24 7.3932396603620498e-24 7.3932396564405566e-24 7.3932396562705492e-24 7.1118158758588142e-14 + 1 5.9806047036494849e-08 2.8142323856670446e-15 5.4329394398523593e-21 7.4194443428333686e-24 7.3934514828234477e-24 7.3935541949253214e-24 7.3933935989348298e-24 7.3932321468740835e-24 7.3932312019180345e-24 7.3932308328287333e-24 7.3932307869392985e-24 7.3932307850741519e-24 7.3932307849933044e-24 7.1118070045957103e-14 + 2 1.2559269877663918e-07 2.8142264165381050e-15 1.1401993450700379e-20 7.4929526498089992e-24 7.3944694218302123e-24 7.3939305039811099e-24 7.3935838893792728e-24 7.3932425308121303e-24 7.3932405326301640e-24 7.3932397568503098e-24 7.3932396603620498e-24 7.3932396564405580e-24 7.3932396562705492e-24 7.1118158758588142e-14 diff --git a/Exec/science/wdmerger/ci-benchmarks/wdmerger_collision_2D.out b/Exec/science/wdmerger/ci-benchmarks/wdmerger_collision_2D.out index 50abd61af7..54c364e26e 100644 --- a/Exec/science/wdmerger/ci-benchmarks/wdmerger_collision_2D.out +++ b/Exec/science/wdmerger/ci-benchmarks/wdmerger_collision_2D.out @@ -1,29 +1,29 @@ plotfile = plt00086 time = 1.25 variables minimum value maximum value - density 8.6936468335e-05 19568007.309 - xmom -5.4959005475e+14 1.3559838981e+14 - ymom -2.5540206053e+15 2.5540214496e+15 + density 8.6940338039e-05 19441641.375 + xmom -5.4953770416e+14 1.3594264808e+14 + ymom -2.4933243003e+15 2.4933250525e+15 zmom 0 0 - rho_E 7.4987846833e+11 5.0676543192e+24 - rho_e 7.1083104678e+11 5.0648015049e+24 - Temp 242292.44287 1409581841.1 - rho_He4 8.6936468335e-17 3.5973411848 - rho_C12 3.4774587333e-05 7826946.398 - rho_O16 5.2161881e-05 11740633.889 - rho_Ne20 8.6936468335e-17 181819.49413 - rho_Mg24 8.6936468335e-17 1190.7712386 - rho_Si28 8.6936468335e-17 6.6817951193 - rho_S32 8.6936468335e-17 0.00019451843176 - rho_Ar36 8.6936468335e-17 1.9568007618e-05 - rho_Ca40 8.6936468335e-17 1.9568007341e-05 - rho_Ti44 8.6936468335e-17 1.9568007318e-05 - rho_Cr48 8.6936468335e-17 1.9568007318e-05 - rho_Fe52 8.6936468335e-17 1.9568007318e-05 - rho_Ni56 8.6936468335e-17 1.9568007318e-05 - phiGrav -5.8709462562e+17 -2.3375498549e+16 - grav_x -685026429.13 -51428.265677 - grav_y -739654246.49 739654206.24 + rho_E 7.4973602186e+11 5.0768248379e+24 + rho_e 7.1068648973e+11 5.0744783673e+24 + Temp 242282.60874 1404450633.1 + rho_He4 8.6940338039e-17 3.398107373 + rho_C12 3.4776135215e-05 7775850.9371 + rho_O16 5.2164202823e-05 11664450.012 + rho_Ne20 8.6940338039e-17 172485.537 + rho_Mg24 8.6940338039e-17 1043.054267 + rho_Si28 8.6940338039e-17 5.9869391361 + rho_S32 8.6940338039e-17 0.00016459247232 + rho_Ar36 8.6940338039e-17 1.9441643669e-05 + rho_Ca40 8.6940338039e-17 1.9441641397e-05 + rho_Ti44 8.6940338039e-17 1.9441641384e-05 + rho_Cr48 8.6940338039e-17 1.9441641384e-05 + rho_Fe52 8.6940338039e-17 1.9441641384e-05 + rho_Ni56 8.6940338039e-17 1.9441641384e-05 + phiGrav -5.870743119e+17 -2.337549858e+16 + grav_x -685044085.4 -51428.268861 + grav_y -739591083.78 739591039.26 grav_z 0 0 - rho_enuc -4.7815621457e+12 7.6360058391e+23 + rho_enuc -7.5385126121e+12 7.1503781318e+23 diff --git a/Source/driver/_cpp_parameters b/Source/driver/_cpp_parameters index f11b20719a..4983d32b86 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 characteristic +# 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/advection_util.cpp b/Source/hydro/advection_util.cpp index 27d3befa3c..515778c4d1 100644 --- a/Source/hydro/advection_util.cpp +++ b/Source/hydro/advection_util.cpp @@ -270,7 +270,6 @@ Castro::divu(const Box& bx, if (i == 0) { ux = 0.0_rt; - vy = 0.0_rt; // is this part correct? } else { Real rl = (i - 0.5_rt) * dx[0] + problo[0]; Real rr = (i + 0.5_rt) * dx[0] + problo[0]; @@ -282,13 +281,13 @@ Castro::divu(const Box& bx, // Take 1/r d/dr(r*u) ux = (rr * ur - rl * ul) * dxinv / rc; + } - // These are transverse averages in the x-direction - Real vb = 0.5_rt * (q_arr(i,j-1,k,QV) + q_arr(i-1,j-1,k,QV)); - Real vt = 0.5_rt * (q_arr(i,j,k,QV) + q_arr(i-1,j,k,QV)); + // These are transverse averages in the x-direction + Real vb = 0.5_rt * (q_arr(i,j-1,k,QV) + q_arr(i-1,j-1,k,QV)); + Real vt = 0.5_rt * (q_arr(i,j,k,QV) + q_arr(i-1,j,k,QV)); - vy = (vt - vb) * dyinv; - } + vy = (vt - vb) * dyinv; } else { diff --git a/Source/hydro/ppm.H b/Source/hydro/ppm.H index 2b7fea274d..352c0cde23 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,18 +55,17 @@ check_trace_source(Array4 const& srcQ, const int idir, /// AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void -ppm_reconstruct(const Real* s, +ppm_reconstruct(const amrex::Real* s, const int i, const int j, const int k, const int idir, - const Real dx, + const amrex::Real dx, const bool lo_bc_test, const bool hi_bc_test, const bool is_axisymmetric, - const GpuArray& domlo, const GpuArray& domhi, - const Real flatn, Real& sm, Real& sp) { + const amrex::GpuArray& domlo, const amrex::GpuArray& domhi, + const amrex::Real flatn, amrex::Real& sm, amrex::Real& sp) { // first we compute s_{i-1/2} -- the left interface value for zone i if (ppm_do_limiting) { - if (lo_bc_test && ((idir == 0 && i == domlo[0]) || (idir == 1 && j == domlo[1]) || (idir == 2 && k == domlo[2]))) { @@ -76,7 +77,7 @@ ppm_reconstruct(const Real* s, if (idir == 0 && is_axisymmetric) { sm = (1.0_rt/144.0_rt)*(415.0_rt*s[i0] - 483.0_rt*s[ip1] + 275.0_rt*s[ip2] - 63.0_rt*s[ip3]); } else { - // Cartesian + // Cartesian sm = (1.0_rt/12.0_rt)*(25.0_rt*s[i0] - 23.0_rt*s[ip1] + 13.0_rt*s[ip2] - 3.0_rt*s[ip3]); } @@ -116,15 +117,15 @@ ppm_reconstruct(const Real* s, if (idir == 0 && is_axisymmetric) { - const Real r_m12 = (static_cast(i) - 0.5_rt) * dx; + const amrex::Real r_m12 = (static_cast(i) - 0.5_rt) * dx; - Real dr4_term = std::pow(dx, 4) * (-12.0_rt * s[im2] + 60.0_rt * s[im1] + 60.0_rt * s[i0] - 12.0_rt * s[ip1]); - Real dr3_term = r_m12 * std::pow(dx, 3) * (-s[im2] - 27.0_rt * s[im1] + 27.0_rt * s[i0] + s[ip1]); - Real dr2_term = std::pow(r_m12, 2) * std::pow(dx, 2) * (30.0_rt * s[im2] - 210.0_rt * s[im1] - 210.0_rt * s[i0] + 30.0_rt * s[ip1]); - Real dr1_term = std::pow(r_m12, 3) * dx * (-s[im2] - 13.0_rt * s[im1] - 13.0_rt * s[i0] + s[ip1]); - Real dr0_term = std::pow(r_m12, 4) * (-10.0_rt * s[im2] + 70.0_rt * s[im1] + 70.0_rt * s[i0] - 10.0_rt * s[ip1]); + amrex::Real dr4_term = std::pow(dx, 4) * (-12.0_rt * s[im2] + 60.0_rt * s[im1] + 60.0_rt * s[i0] - 12.0_rt * s[ip1]); + amrex::Real dr3_term = r_m12 * std::pow(dx, 3) * (-s[im2] - 27.0_rt * s[im1] + 27.0_rt * s[i0] + s[ip1]); + amrex::Real dr2_term = std::pow(r_m12, 2) * std::pow(dx, 2) * (30.0_rt * s[im2] - 210.0_rt * s[im1] - 210.0_rt * s[i0] + 30.0_rt * s[ip1]); + amrex::Real dr1_term = std::pow(r_m12, 3) * dx * (-s[im2] - 13.0_rt * s[im1] - 13.0_rt * s[i0] + s[ip1]); + amrex::Real dr0_term = std::pow(r_m12, 4) * (-10.0_rt * s[im2] + 70.0_rt * s[im1] + 70.0_rt * s[i0] - 10.0_rt * s[ip1]); - Real denom = 24.0_rt * (4.0_rt * std::pow(dx, 4) - 15.0_rt * std::pow(dx, 2) * std::pow(r_m12, 2) + 5.0_rt * std::pow(r_m12, 4)); + amrex::Real denom = 24.0_rt * (4.0_rt * std::pow(dx, 4) - 15.0_rt * std::pow(dx, 2) * std::pow(r_m12, 2) + 5.0_rt * std::pow(r_m12, 4)); sm = (dr4_term + dr3_term + dr2_term + dr1_term + dr0_term) / denom; @@ -132,21 +133,21 @@ ppm_reconstruct(const Real* s, // 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)); } @@ -209,15 +210,15 @@ ppm_reconstruct(const Real* s, if (idir == 0 && is_axisymmetric) { - const Real r_p12 = (static_cast(i) + 0.5_rt) * dx; + const amrex::Real r_p12 = (static_cast(i) + 0.5_rt) * dx; - Real dr4_term = std::pow(dx, 4) * (-12.0_rt * s[im1] + 60.0_rt * s[i0] + 60.0_rt * s[ip1] - 12.0_rt * s[ip2]); - Real dr3_term = r_p12 * std::pow(dx, 3) * (-s[im1] - 27.0_rt * s[i0] + 27.0_rt * s[ip1] + s[ip2]); - Real dr2_term = std::pow(r_p12, 2) * std::pow(dx, 2) * (30.0_rt * s[im1] - 210.0_rt * s[i0] - 210.0_rt * s[ip1] + 30.0_rt * s[ip2]); - Real dr1_term = std::pow(r_p12, 3) * dx * (-s[im1] - 13.0_rt * s[i0] - 13.0_rt * s[ip1] + s[ip2]); - Real dr0_term = std::pow(r_p12, 4) * (-10.0_rt * s[im1] + 70.0_rt * s[i0] + 70.0_rt * s[ip1] - 10.0_rt * s[ip2]); + amrex::Real dr4_term = std::pow(dx, 4) * (-12.0_rt * s[im1] + 60.0_rt * s[i0] + 60.0_rt * s[ip1] - 12.0_rt * s[ip2]); + amrex::Real dr3_term = r_p12 * std::pow(dx, 3) * (-s[im1] - 27.0_rt * s[i0] + 27.0_rt * s[ip1] + s[ip2]); + amrex::Real dr2_term = std::pow(r_p12, 2) * std::pow(dx, 2) * (30.0_rt * s[im1] - 210.0_rt * s[i0] - 210.0_rt * s[ip1] + 30.0_rt * s[ip2]); + amrex::Real dr1_term = std::pow(r_p12, 3) * dx * (-s[im1] - 13.0_rt * s[i0] - 13.0_rt * s[ip1] + s[ip2]); + amrex::Real dr0_term = std::pow(r_p12, 4) * (-10.0_rt * s[im1] + 70.0_rt * s[i0] + 70.0_rt * s[ip1] - 10.0_rt * s[ip2]); - Real denom = 24.0_rt * (4.0_rt * std::pow(dx, 4) - 15.0_rt * std::pow(dx, 2) * std::pow(r_p12, 2) + 5.0_rt * std::pow(r_p12, 4)); + amrex::Real denom = 24.0_rt * (4.0_rt * std::pow(dx, 4) - 15.0_rt * std::pow(dx, 2) * std::pow(r_p12, 2) + 5.0_rt * std::pow(r_p12, 4)); sp = (dr4_term + dr3_term + dr2_term + dr1_term + dr0_term) / denom; @@ -225,21 +226,21 @@ ppm_reconstruct(const Real* s, // Compute van Leer slopes - Real dsl = 2.0_rt * (s[i0] - s[im1]); - Real dsr = 2.0_rt * (s[ip1] - s[i0]); + amrex::Real dsl = 2.0_rt * (s[i0] - s[im1]); + amrex::Real dsr = 2.0_rt * (s[ip1] - s[i0]); - 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[ip1] - s[im1]); + 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]); dsr = 2.0_rt * (s[ip2] - s[ip1]); - 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[ip2] - s[i0]); + 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)); } @@ -292,7 +293,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 @@ -303,29 +304,31 @@ 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 dx, +bool +ppm_reconstruct_pslope(const amrex::Real* rho, const amrex::Real* p, const amrex::Real* src, + const amrex::Real dx, const int i, const int j, const int k, const int idir, bool lo_bc_test, bool hi_bc_test, const bool is_axisymmetric, - const GpuArray& domlo, const GpuArray& domhi, - const Real flatn, - Real& sm, Real& sp) { + const amrex::GpuArray& domlo, const amrex::GpuArray& domhi, + const amrex::Real flatn, + 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.5_rt * dx * (rho[i0] * src[i0] + rho[ip1] * src[ip1]); - Real pp2_hse = pp1_hse + 0.5_rt * dx * (rho[ip1] * src[ip1] + rho[ip2] * src[ip2]); - Real pp3_hse = pp2_hse + 0.5_rt * dx * (rho[ip2] * src[ip2] + rho[ip3] * src[ip3]); + 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]); + amrex::Real pp3_hse = pp2_hse + 0.25_rt*dx * (rho[ip2] + rho[ip3]) * (src[ip2] + src[ip3]); - Real pm1_hse = p0_hse - 0.5_rt * dx * (rho[i0] * src[i0] + rho[im1] * src[im1]); - Real pm2_hse = pm1_hse - 0.5_rt * dx * (rho[im1] * src[im1] + rho[im2] * src[im2]); - Real pm3_hse = pm2_hse - 0.5_rt * dx * (rho[im2] * src[im2] + rho[im3] * src[im3]); + 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]); + amrex::Real pm3_hse = pm2_hse - 0.25_rt*dx * (rho[im2] + rho[im3]) * (src[im2] + src[im3]); // this only makes sense if the pressures are positive bool ptest = p0_hse < 0.0 || @@ -337,7 +340,9 @@ ppm_reconstruct_pslope(const Real* rho, const Real* p, const Real* src, pm3_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 @@ -369,12 +374,16 @@ ppm_reconstruct_pslope(const Real* rho, const Real* p, const Real* src, // 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) { - sp += p[i0] + 0.5_rt * dx * rho[i0] * src[i0]; - sm += p[i0] - 0.5_rt * dx * rho[i0] * src[i0]; + 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; + } @@ -395,14 +404,15 @@ ppm_reconstruct_pslope(const Real* rho, const Real* p, const Real* src, /// 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. // compute x-component of Ip and Im - Real s6 = 6.0_rt * sc - 3.0_rt * (sm + sp); + + amrex::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 @@ -411,8 +421,8 @@ ppm_int_profile(const Real sm, const Real sp, const Real sc, // Im integrates to the left edge of a cell // u-c wave - Real speed = u - c; - Real sigma = std::abs(speed) * dtdx; + amrex::Real speed = u - c; + amrex::Real sigma = std::abs(speed) * dtdx; // if speed == ZERO, then either branch is the same if (speed <= 0.0_rt) { @@ -463,15 +473,15 @@ ppm_int_profile(const Real sm, const Real sp, const Real sc, /// AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void -ppm_int_profile_single(const Real sm, const Real sp, const Real sc, - const Real lam, const Real dtdx, - Real& Ip, Real& Im) { +ppm_int_profile_single(const amrex::Real sm, const amrex::Real sp, const amrex::Real sc, + const amrex::Real lam, const amrex::Real dtdx, + amrex::Real& Ip, amrex::Real& Im) { // Integrate the parabolic profile to the edge of the cell. // This is the MHD version. We come in already with the eigenvalues. // compute x-component of Ip and Im - Real s6 = 6.0_rt * sc - 3.0_rt * (sm + sp); + amrex::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 @@ -479,7 +489,7 @@ ppm_int_profile_single(const Real sm, const Real sp, const Real sc, // Ip integrates to the right edge of a cell // Im integrates to the left edge of a cell - Real sigma = std::abs(lam) * dtdx; + amrex::Real sigma = std::abs(lam) * dtdx; if (lam <= 0.0_rt) { Ip = sp; diff --git a/Source/hydro/trace_ppm.cpp b/Source/hydro/trace_ppm.cpp index 9ef75ff4bd..6c05d74612 100644 --- a/Source/hydro/trace_ppm.cpp +++ b/Source/hydro/trace_ppm.cpp @@ -243,19 +243,37 @@ Castro::trace_ppm(const Box& bx, load_stencil(q_arr, idir, i, j, k, QPRES, s); - if (use_pslope) { + 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 (castro::use_pslope || castro::ppm_well_balanced) { Real trho[nslp]; Real src[nslp]; 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, dL, i, j, k, idir, lo_symm, hi_symm, is_axisymmetric, domlo, domhi, flat, sm, sp); + in_hse = ppm_reconstruct_pslope(trho, s, src, dL, i, j, k, idir, lo_symm, hi_symm, is_axisymmetric, domlo, domhi, flat, 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, dtdL, Ip_p, Im_p); + p_m_hse = s[i0] - 0.5_rt * dL * trho[i0] * src[i0]; + p_p_hse = s[i0] + 0.5_rt * dL * trho[i0] * src[i0]; + + } else { + ppm_int_profile(sm, sp, s[i0], un, cc, dtdL, Ip_p, Im_p); + } } else { ppm_reconstruct(s, i, j, k, idir, dL, lo_symm, hi_symm, is_axisymmetric, domlo, domhi, flat, sm, sp); + ppm_int_profile(sm, sp, s[i0], un, cc, dtdL, Ip_p, Im_p); } - ppm_int_profile(sm, sp, s[i0], un, cc, dtdL, Ip_p, Im_p); // reconstruct rho e @@ -475,6 +493,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]) || @@ -492,23 +514,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 @@ -518,15 +540,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 dup = un_ref - Im_un_2 - hdt*Im_src_un_2; - Real dptotp = p_ref - Im_p[2] - hdt*Im_src_p[2]; + 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 - fac * hdt * Im_src_un_2; + Real dptotp = p_ref - Im_p[2] - hdt * Im_src_p[2]; // {rho, u, p, (rho e)} eigensystem @@ -547,11 +570,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 @@ -560,8 +583,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; } @@ -574,35 +597,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 @@ -624,33 +647,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;