From 3395e2ea13048ef9a6f2564aa8b0a2b790435bd0 Mon Sep 17 00:00:00 2001 From: Zhi Date: Wed, 11 Sep 2024 17:31:15 -0400 Subject: [PATCH 01/25] update mom_flux_has_p in Castro_util --- Source/driver/Castro_util.H | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/Source/driver/Castro_util.H b/Source/driver/Castro_util.H index e7b8751d12..b857e86a33 100644 --- a/Source/driver/Castro_util.H +++ b/Source/driver/Castro_util.H @@ -75,6 +75,14 @@ bool mom_flux_has_p (const int mom_dir, const int flux_dir, const int coord) return true; } + } else if (mom_dir == 1 && flux_dir == 1) { + + if (coord == CoordSys::spherical) { + return false; + } else { + return true + } + } else { return (mom_dir == flux_dir); } From 79562a84bcfce3ac58427002f55c967e24467dce Mon Sep 17 00:00:00 2001 From: Zhi Date: Wed, 11 Sep 2024 18:17:24 -0400 Subject: [PATCH 02/25] change places associated with mom_flux_has_p except castro_ctu/mol_hydro.cpp where ptheta needs to be taken care of later --- Source/hydro/Castro_ctu.cpp | 17 +++++++++++------ Source/hydro/Castro_hydro.H | 2 ++ Source/hydro/Castro_mol.cpp | 16 ++++++++++++---- Source/hydro/Castro_mol_hydro.cpp | 1 + 4 files changed, 26 insertions(+), 10 deletions(-) diff --git a/Source/hydro/Castro_ctu.cpp b/Source/hydro/Castro_ctu.cpp index 7998b7a107..7da960180d 100644 --- a/Source/hydro/Castro_ctu.cpp +++ b/Source/hydro/Castro_ctu.cpp @@ -42,7 +42,7 @@ Castro::consup_hydro(const Box& bx, #endif ) * volinv; - // Add the p div(u) source term to (rho e). + // Add the p div(u) source term to (rho e). if (n == UEINT) { Real pdu = (qx(i+1,j,k,GDPRES) + qx(i,j,k,GDPRES)) * @@ -65,13 +65,18 @@ Castro::consup_hydro(const Box& bx, U_new(i,j,k,n) = U_new(i,j,k,n) - dt * pdu; - } else if (n == UMX) { - // Add gradp term to momentum equation -- only for axisymmetric - // coords (and only for the radial flux). + } else if (n == UMX && !mom_flux_has_p(0, 0, geomdata.Coord())) { + // Add gradp term to momentum equation -- only for axisymmetric + // coords (and only for the radial flux). - if (!mom_flux_has_p(0, 0, geomdata.Coord())) { U_new(i,j,k,UMX) += - dt * (qx(i+1,j,k,GDPRES) - qx(i,j,k,GDPRES)) / geomdata.CellSize()[0]; - } + + } else if (n == UMY && !mom_flux_has_p(1, 1, geomdata.Coord())) { + // Add gradp term to polar(theta) momentum equation for Spherical 2D geometry + + Real r = geomdata.ProbLoArray()[0] + (static_cast(i) + 0.5_rt) * geomdata.Coord()[0]; + U_new(i,j,k,UMY) += - dt * (qy(i,j+1,k,GDPRES) - qy(i,j,k,GDPRES)) / (r * geomdata.CellSize()[0]); + } }); } diff --git a/Source/hydro/Castro_hydro.H b/Source/hydro/Castro_hydro.H index c70fcd64eb..aa8c9582c5 100644 --- a/Source/hydro/Castro_hydro.H +++ b/Source/hydro/Castro_hydro.H @@ -387,6 +387,7 @@ /// @param area1 area of y faces /// @param area2 area of z faces /// @param q0 Godunuv state on x faces +/// @param q1 Godunuv state on y faces /// @param vol cell volume /// void mol_consup(const amrex::Box& bx, @@ -412,6 +413,7 @@ #endif #if AMREX_SPACEDIM <= 2 amrex::Array4 const& q0, + amrex::Array4 const& q1, #endif amrex::Array4 const& vol); diff --git a/Source/hydro/Castro_mol.cpp b/Source/hydro/Castro_mol.cpp index fe2004a60c..49bbf9d2bc 100644 --- a/Source/hydro/Castro_mol.cpp +++ b/Source/hydro/Castro_mol.cpp @@ -226,6 +226,7 @@ Castro::mol_consup(const Box& bx, // NOLINT(readability-convert-member-function #endif #if AMREX_SPACEDIM <= 2 Array4 const& q0, + Array4 const& q1, #endif Array4 const& vol) { @@ -238,6 +239,7 @@ Castro::mol_consup(const Box& bx, // NOLINT(readability-convert-member-function #if AMREX_SPACEDIM <= 2 const auto dx = geom.CellSizeArray(); auto coord = geom.Coord(); + auto prob_lo = geom.ProbLoArray(); #endif amrex::ParallelFor(bx, NUM_STATE, @@ -258,12 +260,18 @@ Castro::mol_consup(const Box& bx, // NOLINT(readability-convert-member-function #endif #if AMREX_SPACEDIM <= 2 - if (n == UMX && do_hydro == 1) { - // Add gradp term to momentum equation -- only for axisymmetric - // coords (and only for the radial flux). + if (do_hydro == 1) { + if (n == UMX && !mom_flux_has_p(0, 0, coord)) { + // Add gradp term to radial momentum equation -- only for axisymmetric + // coords. - if (!mom_flux_has_p(0, 0, coord)) { update(i,j,k,UMX) -= (q0(i+1,j,k,GDPRES) - q0(i,j,k,GDPRES)) / dx[0]; + + } else if (n == UMY && !mom_flux_has_p(1, 1, coord)) { + // Add gradp term to polar(theta) momentum equation for Spherical 2D geometry + + Real r = prob_lo[0] + (static_cast(i) + 0.5_rt)*dx[0]; + update(i,j,k,UMY) -= (q1(i,j+1,k,GDPRES) - q1(i,j,k,GDPRES)) / (r * dx[1]); } } #endif diff --git a/Source/hydro/Castro_mol_hydro.cpp b/Source/hydro/Castro_mol_hydro.cpp index 9fee21334d..75899d67a1 100644 --- a/Source/hydro/Castro_mol_hydro.cpp +++ b/Source/hydro/Castro_mol_hydro.cpp @@ -620,6 +620,7 @@ Castro::construct_mol_hydro_source(Real time, Real dt, MultiFab& A_update) #endif #if AMREX_SPACEDIM <= 2 qe[0].array(), + qe[1].array(), #endif volume.array(mfi)); From ca04f9ec6d9397f8822ada8640518998ae17eaa8 Mon Sep 17 00:00:00 2001 From: Zhi Date: Wed, 11 Sep 2024 18:24:13 -0400 Subject: [PATCH 03/25] fix incorrect calculation --- Source/hydro/Castro_ctu.cpp | 4 ++-- Source/hydro/Castro_mol.cpp | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/Source/hydro/Castro_ctu.cpp b/Source/hydro/Castro_ctu.cpp index 7da960180d..c2f5e3fc40 100644 --- a/Source/hydro/Castro_ctu.cpp +++ b/Source/hydro/Castro_ctu.cpp @@ -74,8 +74,8 @@ Castro::consup_hydro(const Box& bx, } else if (n == UMY && !mom_flux_has_p(1, 1, geomdata.Coord())) { // Add gradp term to polar(theta) momentum equation for Spherical 2D geometry - Real r = geomdata.ProbLoArray()[0] + (static_cast(i) + 0.5_rt) * geomdata.Coord()[0]; - U_new(i,j,k,UMY) += - dt * (qy(i,j+1,k,GDPRES) - qy(i,j,k,GDPRES)) / (r * geomdata.CellSize()[0]); + Real r = geomdata.ProbLoArray()[0] + (static_cast(i) + 0.5_rt) * geomdata.Cellsize()[0]; + U_new(i,j,k,UMY) += - dt * (qy(i,j+1,k,GDPRES) - qy(i,j,k,GDPRES)) / (r * geomdata.CellSize()[1]); } }); diff --git a/Source/hydro/Castro_mol.cpp b/Source/hydro/Castro_mol.cpp index 49bbf9d2bc..688b298a70 100644 --- a/Source/hydro/Castro_mol.cpp +++ b/Source/hydro/Castro_mol.cpp @@ -270,7 +270,7 @@ Castro::mol_consup(const Box& bx, // NOLINT(readability-convert-member-function } else if (n == UMY && !mom_flux_has_p(1, 1, coord)) { // Add gradp term to polar(theta) momentum equation for Spherical 2D geometry - Real r = prob_lo[0] + (static_cast(i) + 0.5_rt)*dx[0]; + Real r = prob_lo[0] + (static_cast(i) + 0.5_rt) * dx[0]; update(i,j,k,UMY) -= (q1(i,j+1,k,GDPRES) - q1(i,j,k,GDPRES)) / (r * dx[1]); } } From 442b5d83937cb9f5a54785a0267e18eeb1d07e92 Mon Sep 17 00:00:00 2001 From: Zhi Date: Wed, 11 Sep 2024 18:30:13 -0400 Subject: [PATCH 04/25] remove redundant else and fix CoordSys::spherical -> CoordSys:SPHERICAL --- Source/driver/Castro_util.H | 8 +++----- 1 file changed, 3 insertions(+), 5 deletions(-) diff --git a/Source/driver/Castro_util.H b/Source/driver/Castro_util.H index b857e86a33..1b05c0cb84 100644 --- a/Source/driver/Castro_util.H +++ b/Source/driver/Castro_util.H @@ -71,17 +71,15 @@ bool mom_flux_has_p (const int mom_dir, const int flux_dir, const int coord) if (coord != CoordSys::cartesian) { return false; - } else { - return true; } + return true; } else if (mom_dir == 1 && flux_dir == 1) { - if (coord == CoordSys::spherical) { + if (coord == CoordSys::SPHERICAL) { return false; - } else { - return true } + return true } else { return (mom_dir == flux_dir); From 8606de8d7f0edd48e44447d5db7b08024b4fd38c Mon Sep 17 00:00:00 2001 From: Zhi Date: Wed, 11 Sep 2024 18:32:36 -0400 Subject: [PATCH 05/25] missing colon --- Source/driver/Castro_util.H | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Source/driver/Castro_util.H b/Source/driver/Castro_util.H index 1b05c0cb84..016215aa32 100644 --- a/Source/driver/Castro_util.H +++ b/Source/driver/Castro_util.H @@ -79,7 +79,7 @@ bool mom_flux_has_p (const int mom_dir, const int flux_dir, const int coord) if (coord == CoordSys::SPHERICAL) { return false; } - return true + return true; } else { return (mom_dir == flux_dir); From 3631380026584de0d83be596969927d0a5ddc50f Mon Sep 17 00:00:00 2001 From: Zhi Date: Wed, 11 Sep 2024 18:37:09 -0400 Subject: [PATCH 06/25] fix typo --- Source/hydro/Castro_ctu.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Source/hydro/Castro_ctu.cpp b/Source/hydro/Castro_ctu.cpp index c2f5e3fc40..7e69b4696b 100644 --- a/Source/hydro/Castro_ctu.cpp +++ b/Source/hydro/Castro_ctu.cpp @@ -74,7 +74,7 @@ Castro::consup_hydro(const Box& bx, } else if (n == UMY && !mom_flux_has_p(1, 1, geomdata.Coord())) { // Add gradp term to polar(theta) momentum equation for Spherical 2D geometry - Real r = geomdata.ProbLoArray()[0] + (static_cast(i) + 0.5_rt) * geomdata.Cellsize()[0]; + Real r = geom.ProbLoArray()[0] + (static_cast(i) + 0.5_rt) * geomdata.CellSize()[0]; U_new(i,j,k,UMY) += - dt * (qy(i,j+1,k,GDPRES) - qy(i,j,k,GDPRES)) / (r * geomdata.CellSize()[1]); } From f9df3a7ed3f1cf44cb7a110f3dff3dacf95e837a Mon Sep 17 00:00:00 2001 From: Zhi Date: Wed, 11 Sep 2024 18:41:09 -0400 Subject: [PATCH 07/25] add preprocessor to ensure we're in 2d or more --- Source/hydro/Castro_ctu.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/Source/hydro/Castro_ctu.cpp b/Source/hydro/Castro_ctu.cpp index 7e69b4696b..9948748880 100644 --- a/Source/hydro/Castro_ctu.cpp +++ b/Source/hydro/Castro_ctu.cpp @@ -71,12 +71,13 @@ Castro::consup_hydro(const Box& bx, U_new(i,j,k,UMX) += - dt * (qx(i+1,j,k,GDPRES) - qx(i,j,k,GDPRES)) / geomdata.CellSize()[0]; +#if AMREX_SPACEDIM >= 2 } else if (n == UMY && !mom_flux_has_p(1, 1, geomdata.Coord())) { // Add gradp term to polar(theta) momentum equation for Spherical 2D geometry Real r = geom.ProbLoArray()[0] + (static_cast(i) + 0.5_rt) * geomdata.CellSize()[0]; U_new(i,j,k,UMY) += - dt * (qy(i,j+1,k,GDPRES) - qy(i,j,k,GDPRES)) / (r * geomdata.CellSize()[1]); - +#endif } }); } From 727f49f75446341ab354be4217e5543216ba69ee Mon Sep 17 00:00:00 2001 From: Zhi Date: Wed, 11 Sep 2024 18:56:58 -0400 Subject: [PATCH 08/25] fix compilation --- Source/hydro/Castro_ctu.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/Source/hydro/Castro_ctu.cpp b/Source/hydro/Castro_ctu.cpp index 9948748880..c46860a3ea 100644 --- a/Source/hydro/Castro_ctu.cpp +++ b/Source/hydro/Castro_ctu.cpp @@ -69,14 +69,14 @@ Castro::consup_hydro(const Box& bx, // Add gradp term to momentum equation -- only for axisymmetric // coords (and only for the radial flux). - U_new(i,j,k,UMX) += - dt * (qx(i+1,j,k,GDPRES) - qx(i,j,k,GDPRES)) / geomdata.CellSize()[0]; + U_new(i,j,k,UMX) += - dt * (qx(i+1,j,k,GDPRES) - qx(i,j,k,GDPRES)) / geomdata.CellSize(1); #if AMREX_SPACEDIM >= 2 } else if (n == UMY && !mom_flux_has_p(1, 1, geomdata.Coord())) { // Add gradp term to polar(theta) momentum equation for Spherical 2D geometry - Real r = geom.ProbLoArray()[0] + (static_cast(i) + 0.5_rt) * geomdata.CellSize()[0]; - U_new(i,j,k,UMY) += - dt * (qy(i,j+1,k,GDPRES) - qy(i,j,k,GDPRES)) / (r * geomdata.CellSize()[1]); + Real r = geomdata.ProbLo(0) + (static_cast(i) + 0.5_rt) * geomdata.CellSize(0); + U_new(i,j,k,UMY) += - dt * (qy(i,j+1,k,GDPRES) - qy(i,j,k,GDPRES)) / (r * geomdata.CellSize(1)); #endif } }); From 0bd9707cfd7101dea0023fc63e91c5032e16185c Mon Sep 17 00:00:00 2001 From: Zhi Date: Thu, 12 Sep 2024 09:41:11 -0400 Subject: [PATCH 09/25] fix typo --- Source/hydro/Castro_ctu.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Source/hydro/Castro_ctu.cpp b/Source/hydro/Castro_ctu.cpp index c46860a3ea..6c222cdddf 100644 --- a/Source/hydro/Castro_ctu.cpp +++ b/Source/hydro/Castro_ctu.cpp @@ -69,7 +69,7 @@ Castro::consup_hydro(const Box& bx, // Add gradp term to momentum equation -- only for axisymmetric // coords (and only for the radial flux). - U_new(i,j,k,UMX) += - dt * (qx(i+1,j,k,GDPRES) - qx(i,j,k,GDPRES)) / geomdata.CellSize(1); + U_new(i,j,k,UMX) += - dt * (qx(i+1,j,k,GDPRES) - qx(i,j,k,GDPRES)) / geomdata.CellSize(0); #if AMREX_SPACEDIM >= 2 } else if (n == UMY && !mom_flux_has_p(1, 1, geomdata.Coord())) { From 1c1df5e9a64f4152825b8c7c1969a3836160fa62 Mon Sep 17 00:00:00 2001 From: Zhi Date: Thu, 12 Sep 2024 14:27:34 -0400 Subject: [PATCH 10/25] update P_theta --- Source/driver/Castro.H | 3 +++ Source/driver/Castro.cpp | 6 ++++++ 2 files changed, 9 insertions(+) diff --git a/Source/driver/Castro.H b/Source/driver/Castro.H index a3c6983c4a..0f3745160f 100644 --- a/Source/driver/Castro.H +++ b/Source/driver/Castro.H @@ -1301,6 +1301,9 @@ protected: #if (AMREX_SPACEDIM <= 2) amrex::MultiFab P_radial; #endif +#if (AMREX_SPACEDIM == 2) + amrex::MultiFab P_theta; +#endif #ifdef RADIATION amrex::Vector > rad_fluxes; #endif diff --git a/Source/driver/Castro.cpp b/Source/driver/Castro.cpp index d3fe48b979..1873e55e34 100644 --- a/Source/driver/Castro.cpp +++ b/Source/driver/Castro.cpp @@ -858,6 +858,12 @@ Castro::initMFs() } #endif +#if (AMREX_SPACEDIM == 2) + if (Geom().IsSPHERICAL()) { + P_theta.define(getEdgeBoxArray(1), dmap, 1, 0); + } +#endif + #ifdef RADIATION if (Radiation::rad_hydro_combined) { rad_fluxes.resize(AMREX_SPACEDIM); From 7bc26066ca3eefa51b99e9c035dc8645943feb78 Mon Sep 17 00:00:00 2001 From: Zhi Date: Mon, 16 Sep 2024 13:04:47 -0400 Subject: [PATCH 11/25] update P_theta flux --- Source/driver/Castro.cpp | 51 +++++++++++++++++++++++ Source/driver/Castro_advance.cpp | 6 +++ Source/driver/Castro_advance_ctu.cpp | 6 +++ Source/hydro/Castro_ctu_hydro.cpp | 54 +++++++++++++++++++----- Source/hydro/Castro_mol_hydro.cpp | 61 ++++++++++++++++++++++------ 5 files changed, 155 insertions(+), 23 deletions(-) diff --git a/Source/driver/Castro.cpp b/Source/driver/Castro.cpp index 1873e55e34..0639904605 100644 --- a/Source/driver/Castro.cpp +++ b/Source/driver/Castro.cpp @@ -2568,6 +2568,12 @@ Castro::FluxRegCrseInit() { } #endif +#if (AMREX_SPACEDIM == 2) + if (Geom().IsSPHERICAL()) { + fine_level.pres_reg.CrseInit(P_theta, 1, 0, 0, 1, pres_crse_scale); + } +#endif + #ifdef RADIATION if (Radiation::rad_hydro_combined) { for (int i = 0; i < AMREX_SPACEDIM; ++i) { @@ -2598,6 +2604,12 @@ Castro::FluxRegFineAdd() { } #endif +#if (AMREX_SPACEDIM == 2) + if (Geom().IsSPHERICAL()) { + getLevel(level).pres_reg.FineAdd(P_theta, 1, 0, 0, 1, pres_fine_scale); + } +#endif + #ifdef RADIATION if (Radiation::rad_hydro_combined) { for (int i = 0; i < AMREX_SPACEDIM; ++i) { @@ -2867,6 +2879,45 @@ Castro::reflux (int crse_level, int fine_level, bool in_post_timestep) } #endif +#if (AMREX_SPACEDIM == 2) + if (Geom().IsSPHERICAL()) { + + reg = &getLevel(lev).pres_reg; + + MultiFab dtheta(crse_lev.grids, crse_lev.dmap, 1, 0); + dtheta.setVal(crse_lev.geom.CellSize(1)); + + // Need to multiply by r for rdtheta + + reg->ClearInternalBorders(crse_lev.geom); + + reg->Reflux(crse_state, rdtheta, 1.0, 0, UMY, 1, crse_lev.geom); + + if (update_sources_after_reflux || !in_post_timestep) { + + MultiFab tmp_fluxes(crse_lev.P_theta.boxArray(), + crse_lev.P_theta.DistributionMap(), + crse_lev.P_theta.nComp(), crse_lev.P_theta.nGrow()); + + tmp_fluxes.setVal(0.0); + + for (OrientationIter fi; fi.isValid(); ++fi) + { + const FabSet& fs = (*reg)[fi()]; + if (fi().coordDir() == 0) { + fs.copyTo(tmp_fluxes, 0, 0, 0, tmp_fluxes.nComp()); + } + } + + MultiFab::Add(crse_lev.P_theta, tmp_fluxes, 0, 0, crse_lev.P_radial.nComp(), 0); + + } + + reg->setVal(0.0); + + } +#endif + #ifdef RADIATION // This follows the same logic as the pure hydro fluxes; see above for details. diff --git a/Source/driver/Castro_advance.cpp b/Source/driver/Castro_advance.cpp index 6f9f0b9d07..b3eb3cdf59 100644 --- a/Source/driver/Castro_advance.cpp +++ b/Source/driver/Castro_advance.cpp @@ -551,6 +551,12 @@ Castro::initialize_advance(Real time, Real dt, int amr_iteration) } #endif +#if (AMREX_SPACEDIM == 2) + if (Geom().IsSPHERICAL()) { + P_theta.setVal(0.0); + } +#endif + #ifdef RADIATION if (Radiation::rad_hydro_combined) { for (int dir = 0; dir < AMREX_SPACEDIM; ++dir) { diff --git a/Source/driver/Castro_advance_ctu.cpp b/Source/driver/Castro_advance_ctu.cpp index af387a92d2..d2fe4fa1da 100644 --- a/Source/driver/Castro_advance_ctu.cpp +++ b/Source/driver/Castro_advance_ctu.cpp @@ -204,6 +204,12 @@ Castro::retry_advance_ctu(Real dt, const advance_status& status) } #endif +#if (AMREX_SPACEDIM == 2) + if (Geom().IsSPHERICAL()) { + P_theta.setVal(0.0); + } +#endif + #ifdef RADIATION if (Radiation::rad_hydro_combined) { for (int dir = 0; dir < AMREX_SPACEDIM; ++dir) { diff --git a/Source/hydro/Castro_ctu_hydro.cpp b/Source/hydro/Castro_ctu_hydro.cpp index d443bbe060..0c5a0ec41e 100644 --- a/Source/hydro/Castro_ctu_hydro.cpp +++ b/Source/hydro/Castro_ctu_hydro.cpp @@ -179,6 +179,9 @@ Castro::construct_ctu_hydro_source(Real time, Real dt) // NOLINT(readability-co #if AMREX_SPACEDIM <= 2 FArrayBox pradial(The_Async_Arena()); #endif +#if AMREX_SPACEDIM == 2 + FArrayBox ptheta(The_Async_Arena()); +#endif #if AMREX_SPACEDIM == 3 FArrayBox qmyx(The_Async_Arena()), qpyx(The_Async_Arena()); FArrayBox qmzx(The_Async_Arena()), qpzx(The_Async_Arena()); @@ -444,6 +447,13 @@ Castro::construct_ctu_hydro_source(Real time, Real dt) // NOLINT(readability-co fab_size += pradial.nBytes(); #endif +#if AMREX_SPACEDIM == 2 + if (Geom().IsSPHERICAL()) { + ptheta.resize(ybx, 1); + } + fab_size += ptheta.nBytes(); +#endif + #ifdef SIMPLIFIED_SDC #ifdef REACTIONS Array4 const sdc_src_arr = SDC_react_source.array(mfi); @@ -1270,23 +1280,33 @@ Castro::construct_ctu_hydro_source(Real time, Real dt) // NOLINT(readability-co scale_rad_flux(nbx, rad_flux_arr, area_arr, dt); #endif - if (idir == 0) { #if AMREX_SPACEDIM <= 2 + // get the scaled radial pressure -- we need to treat this specially + + if (idir == 0 && !mom_flux_has_p(0, 0, coord)) { Array4 pradial_fab = pradial.array(); + + amrex::ParallelFor(nbx, + [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + pradial_fab(i,j,k) = qex_arr(i,j,k,GDPRES) * dt; + }); + } #endif - // get the scaled radial pressure -- we need to treat this specially -#if AMREX_SPACEDIM <= 2 - if (!mom_flux_has_p(0, 0, coord)) { - amrex::ParallelFor(nbx, - [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - pradial_fab(i,j,k) = qex_arr(i,j,k,GDPRES) * dt; - }); - } +#if AMREX_SPACEDIM == 2 + // get the scaled pressure in the theta direction -#endif + if (idir ==1 && !mom_flux_has_p(1, 1, coord)) { + Array4 ptheta_fab = ptheta.array(); + + amrex::ParallelFor(nbx, + [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + ptheta_fab(i,j,k) = qey_arr(i,j,k,GDPRES) * dt; + }); } +#endif // Store the fluxes from this advance. For simplified SDC integration we // only need to do this on the last iteration. @@ -1331,7 +1351,19 @@ Castro::construct_ctu_hydro_source(Real time, Real dt) // NOLINT(readability-co P_radial_fab(i,j,k,0) += pradial_fab(i,j,k,0); }); } +#endif +#if AMREX_SPACEDIM == 2 + if (idir == 1 && !mom_flux_has_p(1, 1, coord)) { + Array4 ptheta_fab = ptheta.array(); + Array4 P_theta_fab = P_theta.array(mfi); + + amrex::ParallelFor(mfi.nodaltilebox(0), + [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + P_theta_fab(i,j,k,0) += ptheta_fab(i,j,k,0); + }); + } #endif } // add_fluxes diff --git a/Source/hydro/Castro_mol_hydro.cpp b/Source/hydro/Castro_mol_hydro.cpp index 75899d67a1..930dabe5a1 100644 --- a/Source/hydro/Castro_mol_hydro.cpp +++ b/Source/hydro/Castro_mol_hydro.cpp @@ -84,6 +84,9 @@ Castro::construct_mol_hydro_source(Real time, Real dt, MultiFab& A_update) } #if AMREX_SPACEDIM <= 2 FArrayBox pradial(The_Async_Arena()); +#endif +#if AMREX_SPACEDIM == 2 + FArrayBox ptheta(The_Async_Arena()); #endif FArrayBox avis(The_Async_Arena()); @@ -657,6 +660,15 @@ Castro::construct_mol_hydro_source(Real time, Real dt, MultiFab& A_update) Array4 pradial_fab = pradial.array(); #endif + +#if AMREX_SPACEDIM == 2 + if (Geom().IsSPHERICAL()) { + ptheta.resize(ybx, 1); + } + + Array4 ptheta_fab = ptheta.array(); +#endif + #if AMREX_SPACEDIM == 1 Array4 const qex_arr = qe[0].array(); #endif @@ -674,23 +686,34 @@ Castro::construct_mol_hydro_source(Real time, Real dt, MultiFab& A_update) #endif flux_arr, area_arr, dt); - - if (idir == 0) { +#if AMREX_SPACEDIM <= 2 + if (idir == 0 && !mom_flux_has_p(0, 0, coord)) { // get the scaled radial pressure -- we need to treat this specially - Array4 const qex_fab = qe[idir].array(); - const int prescomp = GDPRES; + Array4 const qex_fab = qe[idir].array(); -#if AMREX_SPACEDIM <= 2 - if (!mom_flux_has_p(0, 0, coord)) { - amrex::ParallelFor(nbx, - [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - pradial_fab(i,j,k) = qex_fab(i,j,k,prescomp) * dt; - }); - } + amrex::ParallelFor(nbx, + [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + pradial_fab(i,j,k) = qex_fab(i,j,k,GDPRES) * dt; + }); + } #endif + +#if AMREX_SPACEDIM == 2 + if (idir == 1 && !mom_flux_has_p(1, 1, coord)) { + // get the scaled pressure in the theta direction + + Array4 const qey_fab = qe[idir].array(); + + amrex::ParallelFor(nbx, + [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + ptheta_fab(i,j,k) = qey_fab(i,j,k,GDPRES) * dt; + }); } +#endif + } @@ -729,6 +752,20 @@ Castro::construct_mol_hydro_source(Real time, Real dt, MultiFab& A_update) } #endif + +#if AMREX_SPACEDIM == 2 + if (Geom().IsSPHERICAL()) { + + Array4 P_theta_fab = P_theta.array(mfi); + const Real scale = stage_weight; + + AMREX_HOST_DEVICE_FOR_4D(mfi.nodaltilebox(0), 1, i, j, k, n, + { + P_theta_fab(i,j,k,0) += scale * ptheta_fab(i,j,k,0); + }); + + } +#endif } } // MFIter loop From 0ff40aed989f1fee0862849f0d50c9452d6cc4cb Mon Sep 17 00:00:00 2001 From: Zhi Date: Tue, 17 Sep 2024 14:38:05 -0400 Subject: [PATCH 12/25] change dtheta to rdtheta --- Source/driver/Castro.cpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/Source/driver/Castro.cpp b/Source/driver/Castro.cpp index 0639904605..8f7090a024 100644 --- a/Source/driver/Castro.cpp +++ b/Source/driver/Castro.cpp @@ -2884,8 +2884,8 @@ Castro::reflux (int crse_level, int fine_level, bool in_post_timestep) reg = &getLevel(lev).pres_reg; - MultiFab dtheta(crse_lev.grids, crse_lev.dmap, 1, 0); - dtheta.setVal(crse_lev.geom.CellSize(1)); + MultiFab rdtheta(crse_lev.grids, crse_lev.dmap, 1, 0); + rdtheta.setVal(crse_lev.geom.CellSize(1)); // Need to multiply by r for rdtheta @@ -2904,12 +2904,12 @@ Castro::reflux (int crse_level, int fine_level, bool in_post_timestep) for (OrientationIter fi; fi.isValid(); ++fi) { const FabSet& fs = (*reg)[fi()]; - if (fi().coordDir() == 0) { + if (fi().coordDir() == 1) { fs.copyTo(tmp_fluxes, 0, 0, 0, tmp_fluxes.nComp()); } } - MultiFab::Add(crse_lev.P_theta, tmp_fluxes, 0, 0, crse_lev.P_radial.nComp(), 0); + MultiFab::Add(crse_lev.P_theta, tmp_fluxes, 0, 0, crse_lev.P_theta.nComp(), 0); } From 579d949525c3532aa6a2a4a50d692841b9ceef96 Mon Sep 17 00:00:00 2001 From: Zhi Date: Tue, 17 Sep 2024 16:32:15 -0400 Subject: [PATCH 13/25] do rdtheta correctly --- Source/driver/Castro.cpp | 14 ++++++++++++-- 1 file changed, 12 insertions(+), 2 deletions(-) diff --git a/Source/driver/Castro.cpp b/Source/driver/Castro.cpp index 8f7090a024..c499274250 100644 --- a/Source/driver/Castro.cpp +++ b/Source/driver/Castro.cpp @@ -2885,9 +2885,19 @@ Castro::reflux (int crse_level, int fine_level, bool in_post_timestep) reg = &getLevel(lev).pres_reg; MultiFab rdtheta(crse_lev.grids, crse_lev.dmap, 1, 0); - rdtheta.setVal(crse_lev.geom.CellSize(1)); - // Need to multiply by r for rdtheta + auto problo = Geom().ProbLo(); + auto dr = crse_lev.geom.CellSize(0); + auto dtheta = crse_lev.geom.CellSize(1); + + auto const& ma = rdtheta.arrays(); + + amrex:ParallelFor(rdtheta, + [=] AMREX_GPU_DEVICE (int b, int i, int j, int k) noexcept + { + Real r = problo[0] + static_cast(i + 0.5_rt) * dr; + ma[b](i,j,k) = r * dtheta; + }); reg->ClearInternalBorders(crse_lev.geom); From 765969f63e46721043200c0b4c4977311d3b1d7d Mon Sep 17 00:00:00 2001 From: Zhi Date: Tue, 17 Sep 2024 16:43:25 -0400 Subject: [PATCH 14/25] fix compilation --- Source/driver/Castro.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Source/driver/Castro.cpp b/Source/driver/Castro.cpp index c499274250..439b4bdf63 100644 --- a/Source/driver/Castro.cpp +++ b/Source/driver/Castro.cpp @@ -2886,13 +2886,13 @@ Castro::reflux (int crse_level, int fine_level, bool in_post_timestep) MultiFab rdtheta(crse_lev.grids, crse_lev.dmap, 1, 0); - auto problo = Geom().ProbLo(); + const auto* problo = Geom().ProbLo(); auto dr = crse_lev.geom.CellSize(0); auto dtheta = crse_lev.geom.CellSize(1); auto const& ma = rdtheta.arrays(); - amrex:ParallelFor(rdtheta, + amrex::ParallelFor(rdtheta, [=] AMREX_GPU_DEVICE (int b, int i, int j, int k) noexcept { Real r = problo[0] + static_cast(i + 0.5_rt) * dr; From b280340edb3b0eebc332c77b2e42f7fec6e6502e Mon Sep 17 00:00:00 2001 From: zhi Date: Wed, 30 Oct 2024 18:36:29 -0400 Subject: [PATCH 15/25] Revert "implement full well-balancing in PPM (#2945)" This reverts commit db92cbae30e062c73c47bf8c686716b3526d09e0. --- Docs/source/hse.rst | 8 +- Exec/gravity_tests/hse_convergence/README.md | 11 +- .../hse_convergence/convergence_ppm.sh | 27 ----- .../hse_convergence/convergence_sdc.sh | 15 +-- .../ci-benchmarks/job_info_params.txt | 1 - Source/driver/_cpp_parameters | 4 - Source/hydro/ppm.H | 71 +++++------- Source/hydro/trace_ppm.cpp | 107 +++++++----------- 8 files changed, 86 insertions(+), 158 deletions(-) diff --git a/Docs/source/hse.rst b/Docs/source/hse.rst index de34b10b8d..a58acee2ce 100644 --- a/Docs/source/hse.rst +++ b/Docs/source/hse.rst @@ -74,12 +74,6 @@ 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 ------------------------- @@ -171,3 +165,5 @@ 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 d8d909421c..fab88d7bd7 100644 --- a/Exec/gravity_tests/hse_convergence/README.md +++ b/Exec/gravity_tests/hse_convergence/README.md @@ -31,13 +31,10 @@ To run this problem, use one of the convergence scripts: * convergence_sdc.sh : - 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 + 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 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 3dedc8b11c..27b4bb86ef 100755 --- a/Exec/gravity_tests/hse_convergence/convergence_ppm.sh +++ b/Exec/gravity_tests/hse_convergence/convergence_ppm.sh @@ -102,30 +102,3 @@ ${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 407c754759..038109dea1 100755 --- a/Exec/gravity_tests/hse_convergence/convergence_sdc.sh +++ b/Exec/gravity_tests/hse_convergence/convergence_sdc.sh @@ -67,15 +67,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 """ @@ -97,14 +96,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 """ @@ -126,6 +126,7 @@ 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/job_info_params.txt b/Exec/science/flame_wave/ci-benchmarks/job_info_params.txt index f3b20994e7..ff734d8ae0 100644 --- a/Exec/science/flame_wave/ci-benchmarks/job_info_params.txt +++ b/Exec/science/flame_wave/ci-benchmarks/job_info_params.txt @@ -64,7 +64,6 @@ 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/Source/driver/_cpp_parameters b/Source/driver/_cpp_parameters index 4983d32b86..f11b20719a 100644 --- a/Source/driver/_cpp_parameters +++ b/Source/driver/_cpp_parameters @@ -161,10 +161,6 @@ 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/ppm.H b/Source/hydro/ppm.H index 53202100b8..bf854da880 100644 --- a/Source/hydro/ppm.H +++ b/Source/hydro/ppm.H @@ -1,18 +1,16 @@ #include - -#include #include #ifndef PPM_H #define PPM_H -using namespace amrex::literals; +using namespace amrex; using namespace reconstruction; AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE int -check_trace_source(amrex::Array4 const& srcQ, const int idir, +check_trace_source(Array4 const& srcQ, const int idir, const int i, const int j, const int k, const int ncomp) { int do_trace = 0; @@ -55,29 +53,29 @@ check_trace_source(amrex::Array4 const& srcQ, const int idir, /// AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void -ppm_reconstruct(const amrex::Real* s, - const amrex::Real flatn, amrex::Real& sm, amrex::Real& sp) { +ppm_reconstruct(const Real* s, + const Real flatn, Real& sm, Real& sp) { if (ppm_do_limiting) { // Compute van Leer slopes - amrex::Real dsl = 2.0_rt * (s[im1] - s[im2]); - amrex::Real dsr = 2.0_rt * (s[i0] - s[im1]); + Real dsl = 2.0_rt * (s[im1] - s[im2]); + Real dsr = 2.0_rt * (s[i0] - s[im1]); - amrex::Real dsvl_l = 0.0_rt; + Real dsvl_l = 0.0_rt; if (dsl*dsr > 0.0_rt) { - amrex::Real dsc = 0.5_rt * (s[i0] - s[im2]); + 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]); - amrex::Real dsvl_r = 0.0_rt; + Real dsvl_r = 0.0_rt; if (dsl*dsr > 0.0_rt) { - amrex::Real dsc = 0.5_rt * (s[ip1] - s[im1]); + 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)); } @@ -96,8 +94,8 @@ ppm_reconstruct(const amrex::Real* s, dsvl_l = 0.0_rt; if (dsl*dsr > 0.0_rt) { - 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)); + 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]); @@ -105,8 +103,8 @@ ppm_reconstruct(const amrex::Real* s, dsvl_r = 0.0_rt; if (dsl*dsr > 0.0_rt) { - 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)); + 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 @@ -155,7 +153,7 @@ ppm_reconstruct(const amrex::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 @@ -166,25 +164,23 @@ ppm_reconstruct(const amrex::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 -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) { +void +ppm_reconstruct_pslope(const Real* rho, const Real* p, const Real* src, const Real flatn, + const Real dx, + Real& sm, Real& sp) { - amrex::Real tp[nslp]; + Real tp[nslp]; // compute the hydrostatic pressure in each zone center starting with i0 - amrex::Real p0_hse = p[i0]; + Real p0_hse = p[i0]; - 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 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 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]); + 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]); // this only makes sense if the pressures are positive bool ptest = p0_hse < 0.0 || @@ -194,9 +190,7 @@ ppm_reconstruct_pslope(const amrex::Real* rho, const amrex::Real* p, const amrex pm2_hse < 0.0; - bool do_hse = !ptest && rho[i0] >= pslope_cutoff_density; - - if (do_hse) { + if (!ptest && rho[i0] >= pslope_cutoff_density) { // subtract off the hydrostatic pressure @@ -224,16 +218,12 @@ ppm_reconstruct_pslope(const amrex::Real* rho, const amrex::Real* p, const amrex // 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 (do_hse && !castro::ppm_well_balanced) { + if (!ptest && rho[i0] >= pslope_cutoff_density) { sp += p[i0] + 0.5 * dx * rho[i0] * src[i0]; sm += p[i0] - 0.5 * dx * rho[i0] * src[i0]; } - return do_hse; - } @@ -254,14 +244,13 @@ ppm_reconstruct_pslope(const amrex::Real* rho, const amrex::Real* p, const amrex /// AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void -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) { +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) { // 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); // Ip/m is the integral under the parabola for the extent diff --git a/Source/hydro/trace_ppm.cpp b/Source/hydro/trace_ppm.cpp index a824fbb56f..a4fa02e6e5 100644 --- a/Source/hydro/trace_ppm.cpp +++ b/Source/hydro/trace_ppm.cpp @@ -229,37 +229,19 @@ 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 (castro::use_pslope || castro::ppm_well_balanced) { + if (use_pslope) { 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); - in_hse = ppm_reconstruct_pslope(trho, s, src, flat, dL, 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); - } + ppm_reconstruct_pslope(trho, s, src, flat, dL, sm, sp); } else { ppm_reconstruct(s, 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 @@ -477,10 +459,6 @@ 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]) || @@ -498,23 +476,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 + fac * hdt * Im_src_un_0; + Real un_ref = Im_un_0 + 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 = std::max(rho_ref, lsmall_dens); + rho_ref = amrex::max(rho_ref, lsmall_dens); Real rho_ref_inv = 1.0_rt/rho_ref; - p_ref = std::max(p_ref, lsmall_pres); + p_ref = amrex::max(p_ref, lsmall_pres); // For tracing - Real csq_ref = gam_g_ref * (p_m_hse + p_ref) * rho_ref_inv; + Real csq_ref = gam_g_ref*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_m_hse + p_ref) + rhoe_g_ref) * rho_ref_inv; + Real h_g_ref = (p_ref + rhoe_g_ref)*rho_ref_inv; // *m are the jumps carried by un-c // *p are the jumps carried by un+c @@ -524,16 +502,15 @@ 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 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 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]; + Real dup = un_ref - Im_un_2 - hdt*Im_src_un_2; + Real dptotp = p_ref - Im_p[2] - hdt*Im_src_p[2]; // {rho, u, p, (rho e)} eigensystem @@ -554,11 +531,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) = std::max(lsmall_dens, rho_ref + alphap + alpham + alpha0r); + qp(i,j,k,QRHO) = amrex::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) = std::max(castro::small_dens * castro::small_ener, + qp(i,j,k,QREINT) = amrex::max(castro::small_dens * castro::small_ener, rhoe_g_ref + (alphap + alpham)*h_g_ref + alpha0e_g); - qp(i,j,k,QPRES) = std::max(lsmall_pres, p_m_hse + p_ref + (alphap + alpham)*csq_ref); + qp(i,j,k,QPRES) = amrex::max(lsmall_pres, 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 @@ -567,8 +544,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; } @@ -581,35 +558,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 + fac * hdt * Ip_src_un_2; + Real un_ref = Ip_un_2 + 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 = std::max(rho_ref, lsmall_dens); + rho_ref = amrex::max(rho_ref, lsmall_dens); Real rho_ref_inv = 1.0_rt/rho_ref; - p_ref = std::max(p_ref, lsmall_pres); + p_ref = amrex::max(p_ref, lsmall_pres); // For tracing - Real csq_ref = gam_g_ref * (p_p_hse + p_ref) * rho_ref_inv; + Real csq_ref = gam_g_ref*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_p_hse + p_ref) + rhoe_g_ref) * rho_ref_inv; + Real h_g_ref = (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 - fac * hdt * Ip_src_un_0; - Real dptotm = p_ref - Ip_p[0] - hdt * Ip_src_p[0]; + 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 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 - fac * hdt * Ip_src_un_2; - Real dptotp = p_ref - Ip_p[2] - hdt * Ip_src_p[2]; + Real dup = un_ref - Ip_un_2 - hdt*Ip_src_un_2; + Real dptotp = p_ref - Ip_p[2] - hdt*Ip_src_p[2]; // {rho, u, p, (rho e)} eigensystem @@ -631,33 +608,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) = std::max(lsmall_dens, rho_ref + alphap + alpham + alpha0r); + qm(i+1,j,k,QRHO) = amrex::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) = std::max(castro::small_dens * castro::small_ener, + qm(i+1,j,k,QREINT) = amrex::max(castro::small_dens * castro::small_ener, rhoe_g_ref + (alphap + alpham)*h_g_ref + alpha0e_g); - qm(i+1,j,k,QPRES) = std::max(lsmall_pres, p_p_hse + p_ref + (alphap + alpham)*csq_ref); + qm(i+1,j,k,QPRES) = amrex::max(lsmall_pres, 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) = std::max(lsmall_dens, rho_ref + alphap + alpham + alpha0r); + qm(i,j+1,k,QRHO) = amrex::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) = 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); + 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); // 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) = std::max(lsmall_dens, rho_ref + alphap + alpham + alpha0r); + qm(i,j,k+1,QRHO) = amrex::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) = std::max(castro::small_dens * castro::small_ener, + qm(i,j,k+1,QREINT) = amrex::max(castro::small_dens * castro::small_ener, rhoe_g_ref + (alphap + alpham)*h_g_ref + alpha0e_g); - qm(i,j,k+1,QPRES) = std::max(lsmall_pres, p_p_hse + p_ref + (alphap + alpham)*csq_ref); + qm(i,j,k+1,QPRES) = amrex::max(lsmall_pres, p_ref + (alphap + alpham)*csq_ref); // transverse velocities qm(i,j,k+1,QUT) = Ip_ut_1 + hdt*Ip_src_ut_1; From 199ebeb091c359df51a12ed0c850b1660d58de10 Mon Sep 17 00:00:00 2001 From: zhi Date: Wed, 30 Oct 2024 18:53:45 -0400 Subject: [PATCH 16/25] Reapply "implement full well-balancing in PPM (#2945)" This reverts commit b280340edb3b0eebc332c77b2e42f7fec6e6502e. --- Docs/source/hse.rst | 8 +- Exec/gravity_tests/hse_convergence/README.md | 11 +- .../hse_convergence/convergence_ppm.sh | 27 +++++ .../hse_convergence/convergence_sdc.sh | 15 ++- .../ci-benchmarks/job_info_params.txt | 1 + Source/driver/_cpp_parameters | 4 + Source/hydro/ppm.H | 71 +++++++----- Source/hydro/trace_ppm.cpp | 107 +++++++++++------- 8 files changed, 158 insertions(+), 86 deletions(-) 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/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/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/ppm.H b/Source/hydro/ppm.H index bf854da880..53202100b8 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,13 +254,14 @@ 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. // 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 diff --git a/Source/hydro/trace_ppm.cpp b/Source/hydro/trace_ppm.cpp index a4fa02e6e5..a824fbb56f 100644 --- a/Source/hydro/trace_ppm.cpp +++ b/Source/hydro/trace_ppm.cpp @@ -229,19 +229,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, flat, dL, sm, sp); + in_hse = ppm_reconstruct_pslope(trho, s, src, flat, dL, 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, 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 @@ -459,6 +477,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]) || @@ -476,23 +498,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 @@ -502,15 +524,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 @@ -531,11 +554,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 @@ -544,8 +567,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; } @@ -558,35 +581,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 @@ -608,33 +631,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; From a9dfeaa0b3616363c88575667e7e0e89b8923168 Mon Sep 17 00:00:00 2001 From: zhi Date: Wed, 30 Oct 2024 21:49:33 -0400 Subject: [PATCH 17/25] add space --- Source/hydro/Castro_ctu_hydro.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Source/hydro/Castro_ctu_hydro.cpp b/Source/hydro/Castro_ctu_hydro.cpp index 0c5a0ec41e..4ad7a16a81 100644 --- a/Source/hydro/Castro_ctu_hydro.cpp +++ b/Source/hydro/Castro_ctu_hydro.cpp @@ -1297,7 +1297,7 @@ Castro::construct_ctu_hydro_source(Real time, Real dt) // NOLINT(readability-co #if AMREX_SPACEDIM == 2 // get the scaled pressure in the theta direction - if (idir ==1 && !mom_flux_has_p(1, 1, coord)) { + if (idir == 1 && !mom_flux_has_p(1, 1, coord)) { Array4 ptheta_fab = ptheta.array(); amrex::ParallelFor(nbx, From ebb142cb26f1c64845a606c5fecaafa1dace85a3 Mon Sep 17 00:00:00 2001 From: zhi Date: Tue, 12 Nov 2024 13:44:17 -0500 Subject: [PATCH 18/25] update pressure scaling with area and coarse volume --- Source/driver/Castro.cpp | 22 ++-------------------- Source/hydro/Castro_ctu_hydro.cpp | 4 ++-- Source/hydro/Castro_mol_hydro.cpp | 4 ++-- 3 files changed, 6 insertions(+), 24 deletions(-) diff --git a/Source/driver/Castro.cpp b/Source/driver/Castro.cpp index 448fbf31f0..184f930340 100644 --- a/Source/driver/Castro.cpp +++ b/Source/driver/Castro.cpp @@ -2875,12 +2875,9 @@ Castro::reflux (int crse_level, int fine_level, bool in_post_timestep) reg = &getLevel(lev).pres_reg; - MultiFab dr(crse_lev.grids, crse_lev.dmap, 1, 0); - dr.setVal(crse_lev.geom.CellSize(0)); - reg->ClearInternalBorders(crse_lev.geom); - reg->Reflux(crse_state, dr, 1.0, 0, UMX, 1, crse_lev.geom); + reg->Reflux(crse_state, crse_lev.volume, 1.0, 0, UMX, 1, crse_lev.geom); if (update_sources_after_reflux || !in_post_timestep) { @@ -2912,24 +2909,9 @@ Castro::reflux (int crse_level, int fine_level, bool in_post_timestep) reg = &getLevel(lev).pres_reg; - MultiFab rdtheta(crse_lev.grids, crse_lev.dmap, 1, 0); - - const auto* problo = Geom().ProbLo(); - auto dr = crse_lev.geom.CellSize(0); - auto dtheta = crse_lev.geom.CellSize(1); - - auto const& ma = rdtheta.arrays(); - - amrex::ParallelFor(rdtheta, - [=] AMREX_GPU_DEVICE (int b, int i, int j, int k) noexcept - { - Real r = problo[0] + static_cast(i + 0.5_rt) * dr; - ma[b](i,j,k) = r * dtheta; - }); - reg->ClearInternalBorders(crse_lev.geom); - reg->Reflux(crse_state, rdtheta, 1.0, 0, UMY, 1, crse_lev.geom); + reg->Reflux(crse_state, crse_lev.volume, 1.0, 0, UMY, 1, crse_lev.geom); if (update_sources_after_reflux || !in_post_timestep) { diff --git a/Source/hydro/Castro_ctu_hydro.cpp b/Source/hydro/Castro_ctu_hydro.cpp index 4ad7a16a81..b3b4befb76 100644 --- a/Source/hydro/Castro_ctu_hydro.cpp +++ b/Source/hydro/Castro_ctu_hydro.cpp @@ -1289,7 +1289,7 @@ Castro::construct_ctu_hydro_source(Real time, Real dt) // NOLINT(readability-co amrex::ParallelFor(nbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - pradial_fab(i,j,k) = qex_arr(i,j,k,GDPRES) * dt; + pradial_fab(i,j,k) = area_arr(i,j,k) * qex_arr(i,j,k,GDPRES) * dt; }); } #endif @@ -1303,7 +1303,7 @@ Castro::construct_ctu_hydro_source(Real time, Real dt) // NOLINT(readability-co amrex::ParallelFor(nbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - ptheta_fab(i,j,k) = qey_arr(i,j,k,GDPRES) * dt; + ptheta_fab(i,j,k) = area_arr(i,j,k) * qey_arr(i,j,k,GDPRES) * dt; }); } #endif diff --git a/Source/hydro/Castro_mol_hydro.cpp b/Source/hydro/Castro_mol_hydro.cpp index 930dabe5a1..0b9a3c518a 100644 --- a/Source/hydro/Castro_mol_hydro.cpp +++ b/Source/hydro/Castro_mol_hydro.cpp @@ -695,7 +695,7 @@ Castro::construct_mol_hydro_source(Real time, Real dt, MultiFab& A_update) amrex::ParallelFor(nbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - pradial_fab(i,j,k) = qex_fab(i,j,k,GDPRES) * dt; + pradial_fab(i,j,k) = area_arr(i,j,k) * qex_fab(i,j,k,GDPRES) * dt; }); } #endif @@ -709,7 +709,7 @@ Castro::construct_mol_hydro_source(Real time, Real dt, MultiFab& A_update) amrex::ParallelFor(nbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - ptheta_fab(i,j,k) = qey_fab(i,j,k,GDPRES) * dt; + ptheta_fab(i,j,k) = area_arr(i,j,k) * qey_fab(i,j,k,GDPRES) * dt; }); } #endif From 28a74b5aebecff78393fe8a764d7268c2ac1bb19 Mon Sep 17 00:00:00 2001 From: Zhi Date: Tue, 12 Nov 2024 15:00:20 -0500 Subject: [PATCH 19/25] update benchmark --- .../ci-benchmarks/wdmerger_collision_2D.out | 46 +++++++++---------- 1 file changed, 23 insertions(+), 23 deletions(-) diff --git a/Exec/science/wdmerger/ci-benchmarks/wdmerger_collision_2D.out b/Exec/science/wdmerger/ci-benchmarks/wdmerger_collision_2D.out index 54c364e26e..026c16a299 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.6940338039e-05 19441641.375 - xmom -5.4953770416e+14 1.3594264808e+14 - ymom -2.4933243003e+15 2.4933250525e+15 + density 8.6944011764e-05 19444371.155 + xmom -5.4917303392e+14 1.3509526313e+14 + ymom -2.4930150346e+15 2.4930157613e+15 zmom 0 0 - 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 + rho_E 7.4445974314e+11 5.0751963975e+24 + rho_e 7.0573031626e+11 5.0728530525e+24 + Temp 241944.6157 1402008261.1 + rho_He4 8.6944011764e-17 3.3044554838 + rho_C12 3.4777604705e-05 7776963.7255 + rho_O16 5.2166407058e-05 11666101.692 + rho_Ne20 8.6944011764e-17 167487.84509 + rho_Mg24 8.6944011764e-17 977.92441758 + rho_Si28 8.6944011764e-17 5.6766266908 + rho_S32 8.6944011764e-17 0.0001524753087 + rho_Ar36 8.6944011764e-17 1.9444373276e-05 + rho_Ca40 8.6944011764e-17 1.9444371178e-05 + rho_Ti44 8.6944011764e-17 1.9444371165e-05 + rho_Cr48 8.6944011764e-17 1.9444371164e-05 + rho_Fe52 8.6944011764e-17 1.9444371164e-05 + rho_Ni56 8.6944011764e-17 1.9444371164e-05 + phiGrav -5.8708220182e+17 -2.3375498709e+16 + grav_x -685053820.85 -51428.28215 + grav_y -739593455 739593428.31 grav_z 0 0 - rho_enuc -7.5385126121e+12 7.1503781318e+23 + rho_enuc -6.2415111706e+12 6.8835763053e+23 From 47af504c5c17d04102f398f1469c91bcaa0bec98 Mon Sep 17 00:00:00 2001 From: Zhi Date: Tue, 12 Nov 2024 15:16:14 -0500 Subject: [PATCH 20/25] Revert "update benchmark" This reverts commit 28a74b5aebecff78393fe8a764d7268c2ac1bb19. --- .../ci-benchmarks/wdmerger_collision_2D.out | 46 +++++++++---------- 1 file changed, 23 insertions(+), 23 deletions(-) diff --git a/Exec/science/wdmerger/ci-benchmarks/wdmerger_collision_2D.out b/Exec/science/wdmerger/ci-benchmarks/wdmerger_collision_2D.out index 026c16a299..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.6944011764e-05 19444371.155 - xmom -5.4917303392e+14 1.3509526313e+14 - ymom -2.4930150346e+15 2.4930157613e+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.4445974314e+11 5.0751963975e+24 - rho_e 7.0573031626e+11 5.0728530525e+24 - Temp 241944.6157 1402008261.1 - rho_He4 8.6944011764e-17 3.3044554838 - rho_C12 3.4777604705e-05 7776963.7255 - rho_O16 5.2166407058e-05 11666101.692 - rho_Ne20 8.6944011764e-17 167487.84509 - rho_Mg24 8.6944011764e-17 977.92441758 - rho_Si28 8.6944011764e-17 5.6766266908 - rho_S32 8.6944011764e-17 0.0001524753087 - rho_Ar36 8.6944011764e-17 1.9444373276e-05 - rho_Ca40 8.6944011764e-17 1.9444371178e-05 - rho_Ti44 8.6944011764e-17 1.9444371165e-05 - rho_Cr48 8.6944011764e-17 1.9444371164e-05 - rho_Fe52 8.6944011764e-17 1.9444371164e-05 - rho_Ni56 8.6944011764e-17 1.9444371164e-05 - phiGrav -5.8708220182e+17 -2.3375498709e+16 - grav_x -685053820.85 -51428.28215 - grav_y -739593455 739593428.31 + 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 -6.2415111706e+12 6.8835763053e+23 + rho_enuc -7.5385126121e+12 7.1503781318e+23 From 99464ed249c4b306fa08e4b700d9bd452a081129 Mon Sep 17 00:00:00 2001 From: Zhi Date: Tue, 12 Nov 2024 15:36:48 -0500 Subject: [PATCH 21/25] update benchmark --- .../ci-benchmarks/wdmerger_collision_2D.out | 46 +++++++++---------- 1 file changed, 23 insertions(+), 23 deletions(-) diff --git a/Exec/science/wdmerger/ci-benchmarks/wdmerger_collision_2D.out b/Exec/science/wdmerger/ci-benchmarks/wdmerger_collision_2D.out index 54c364e26e..026c16a299 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.6940338039e-05 19441641.375 - xmom -5.4953770416e+14 1.3594264808e+14 - ymom -2.4933243003e+15 2.4933250525e+15 + density 8.6944011764e-05 19444371.155 + xmom -5.4917303392e+14 1.3509526313e+14 + ymom -2.4930150346e+15 2.4930157613e+15 zmom 0 0 - 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 + rho_E 7.4445974314e+11 5.0751963975e+24 + rho_e 7.0573031626e+11 5.0728530525e+24 + Temp 241944.6157 1402008261.1 + rho_He4 8.6944011764e-17 3.3044554838 + rho_C12 3.4777604705e-05 7776963.7255 + rho_O16 5.2166407058e-05 11666101.692 + rho_Ne20 8.6944011764e-17 167487.84509 + rho_Mg24 8.6944011764e-17 977.92441758 + rho_Si28 8.6944011764e-17 5.6766266908 + rho_S32 8.6944011764e-17 0.0001524753087 + rho_Ar36 8.6944011764e-17 1.9444373276e-05 + rho_Ca40 8.6944011764e-17 1.9444371178e-05 + rho_Ti44 8.6944011764e-17 1.9444371165e-05 + rho_Cr48 8.6944011764e-17 1.9444371164e-05 + rho_Fe52 8.6944011764e-17 1.9444371164e-05 + rho_Ni56 8.6944011764e-17 1.9444371164e-05 + phiGrav -5.8708220182e+17 -2.3375498709e+16 + grav_x -685053820.85 -51428.28215 + grav_y -739593455 739593428.31 grav_z 0 0 - rho_enuc -7.5385126121e+12 7.1503781318e+23 + rho_enuc -6.2415111706e+12 6.8835763053e+23 From b11b943cd82b8d812d76a5d2beb860a029c85fa0 Mon Sep 17 00:00:00 2001 From: Zhi Date: Tue, 12 Nov 2024 15:55:08 -0500 Subject: [PATCH 22/25] updae benchmark --- .../ci-benchmarks/wdmerger_collision_2D.out | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/Exec/science/wdmerger/ci-benchmarks/wdmerger_collision_2D.out b/Exec/science/wdmerger/ci-benchmarks/wdmerger_collision_2D.out index 026c16a299..d951913be1 100644 --- a/Exec/science/wdmerger/ci-benchmarks/wdmerger_collision_2D.out +++ b/Exec/science/wdmerger/ci-benchmarks/wdmerger_collision_2D.out @@ -5,16 +5,16 @@ xmom -5.4917303392e+14 1.3509526313e+14 ymom -2.4930150346e+15 2.4930157613e+15 zmom 0 0 - rho_E 7.4445974314e+11 5.0751963975e+24 + rho_E 7.4445974314e+11 5.0751963974e+24 rho_e 7.0573031626e+11 5.0728530525e+24 Temp 241944.6157 1402008261.1 - rho_He4 8.6944011764e-17 3.3044554838 + rho_He4 8.6944011764e-17 3.3044554839 rho_C12 3.4777604705e-05 7776963.7255 rho_O16 5.2166407058e-05 11666101.692 - rho_Ne20 8.6944011764e-17 167487.84509 - rho_Mg24 8.6944011764e-17 977.92441758 - rho_Si28 8.6944011764e-17 5.6766266908 - rho_S32 8.6944011764e-17 0.0001524753087 + rho_Ne20 8.6944011764e-17 167487.84482 + rho_Mg24 8.6944011764e-17 977.92440606 + rho_Si28 8.6944011764e-17 5.6766266345 + rho_S32 8.6944011764e-17 0.00015247530548 rho_Ar36 8.6944011764e-17 1.9444373276e-05 rho_Ca40 8.6944011764e-17 1.9444371178e-05 rho_Ti44 8.6944011764e-17 1.9444371165e-05 @@ -25,5 +25,5 @@ grav_x -685053820.85 -51428.28215 grav_y -739593455 739593428.31 grav_z 0 0 - rho_enuc -6.2415111706e+12 6.8835763053e+23 + rho_enuc -9.0785617027e+12 6.8835762986e+23 From fc7cfce6a08ce3d09bd2cd0f4169046976774458 Mon Sep 17 00:00:00 2001 From: zhi Date: Thu, 14 Nov 2024 19:43:19 -0500 Subject: [PATCH 23/25] fix ptheta reflux --- Source/driver/Castro.cpp | 19 +++++++------------ 1 file changed, 7 insertions(+), 12 deletions(-) diff --git a/Source/driver/Castro.cpp b/Source/driver/Castro.cpp index 184f930340..b55e38d9ec 100644 --- a/Source/driver/Castro.cpp +++ b/Source/driver/Castro.cpp @@ -2875,9 +2875,9 @@ Castro::reflux (int crse_level, int fine_level, bool in_post_timestep) reg = &getLevel(lev).pres_reg; - reg->ClearInternalBorders(crse_lev.geom); + Reg->ClearInternalBorders(crse_lev.geom); - reg->Reflux(crse_state, crse_lev.volume, 1.0, 0, UMX, 1, crse_lev.geom); + reg->Reflux(crse_state, crse_lev.volume, 0, 1.0, 0, UMX, 1, crse_lev.geom); if (update_sources_after_reflux || !in_post_timestep) { @@ -2899,19 +2899,10 @@ Castro::reflux (int crse_level, int fine_level, bool in_post_timestep) } - reg->setVal(0.0); - - } -#endif - #if (AMREX_SPACEDIM == 2) if (Geom().IsSPHERICAL()) { - reg = &getLevel(lev).pres_reg; - - reg->ClearInternalBorders(crse_lev.geom); - - reg->Reflux(crse_state, crse_lev.volume, 1.0, 0, UMY, 1, crse_lev.geom); + reg->Reflux(crse_state, crse_lev.volume, 1, 1.0, 0, UMY, 1, crse_lev.geom); if (update_sources_after_reflux || !in_post_timestep) { @@ -2933,11 +2924,15 @@ Castro::reflux (int crse_level, int fine_level, bool in_post_timestep) } + } +#endif + reg->setVal(0.0); } #endif + #ifdef RADIATION // This follows the same logic as the pure hydro fluxes; see above for details. From 8bd8b60aefed17cd590b96e92db82cec2c88d548 Mon Sep 17 00:00:00 2001 From: zhi Date: Thu, 14 Nov 2024 19:50:22 -0500 Subject: [PATCH 24/25] fix typo --- Source/driver/Castro.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Source/driver/Castro.cpp b/Source/driver/Castro.cpp index b55e38d9ec..c29fdaa8f0 100644 --- a/Source/driver/Castro.cpp +++ b/Source/driver/Castro.cpp @@ -2875,7 +2875,7 @@ Castro::reflux (int crse_level, int fine_level, bool in_post_timestep) reg = &getLevel(lev).pres_reg; - Reg->ClearInternalBorders(crse_lev.geom); + reg->ClearInternalBorders(crse_lev.geom); reg->Reflux(crse_state, crse_lev.volume, 0, 1.0, 0, UMX, 1, crse_lev.geom); From 727f52e749be64aa01d231468ef39dddd566e369 Mon Sep 17 00:00:00 2001 From: zhi Date: Mon, 18 Nov 2024 13:13:21 -0500 Subject: [PATCH 25/25] fix compile --- Source/hydro/Castro_mol_hydro.cpp | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/Source/hydro/Castro_mol_hydro.cpp b/Source/hydro/Castro_mol_hydro.cpp index e1d06fc4b1..0772b71b79 100644 --- a/Source/hydro/Castro_mol_hydro.cpp +++ b/Source/hydro/Castro_mol_hydro.cpp @@ -657,16 +657,12 @@ Castro::construct_mol_hydro_source(Real time, Real dt, MultiFab& A_update) if (!Geom().IsCartesian()) { pradial.resize(xbx, 1); } - - Array4 pradial_fab = pradial.array(); #endif #if AMREX_SPACEDIM == 2 if (Geom().IsSPHERICAL()) { ptheta.resize(ybx, 1); } - - Array4 ptheta_fab = ptheta.array(); #endif for (int idir = 0; idir < AMREX_SPACEDIM; ++idir) { @@ -682,26 +678,28 @@ Castro::construct_mol_hydro_source(Real time, Real dt, MultiFab& A_update) // get the scaled radial pressure -- we need to treat this specially if (idir == 0 && !mom_flux_has_p(0, 0, coord)) { + Array4 pradial_fab = pradial.array(); Array4 const qex_arr = qe[idir].array(); amrex::ParallelFor(nbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - pradial_fab(i,j,k) = area_arr(i,j,k) * qex_fab(i,j,k,GDPRES) * dt; + pradial_fab(i,j,k) = area_arr(i,j,k) * qex_arr(i,j,k,GDPRES) * dt; }); } #endif #if AMREX_SPACEDIM == 2 - if (idir == 1 && !mom_flux_has_p(1, 1, coord)) { - // get the scaled pressure in the theta direction + // get the scaled pressure in the theta direction - Array4 const qey_fab = qe[idir].array(); + if (idir == 1 && !mom_flux_has_p(1, 1, coord)) { + Array4 ptheta_fab = ptheta.array(); + Array4 const qey_arr = qe[idir].array(); amrex::ParallelFor(nbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - ptheta_fab(i,j,k) = area_arr(i,j,k) * qey_fab(i,j,k,GDPRES) * dt; + ptheta_fab(i,j,k) = area_arr(i,j,k) * qey_arr(i,j,k,GDPRES) * dt; }); } #endif @@ -734,6 +732,7 @@ Castro::construct_mol_hydro_source(Real time, Real dt, MultiFab& A_update) #if AMREX_SPACEDIM <= 2 if (!Geom().IsCartesian()) { + Array4 pradial_fab = pradial.array(); Array4 P_radial_fab = P_radial.array(mfi); const Real scale = stage_weight; @@ -748,6 +747,7 @@ Castro::construct_mol_hydro_source(Real time, Real dt, MultiFab& A_update) #if AMREX_SPACEDIM == 2 if (Geom().IsSPHERICAL()) { + Array4 ptheta_fab = ptheta.array(); Array4 P_theta_fab = P_theta.array(mfi); const Real scale = stage_weight;