From 3395e2ea13048ef9a6f2564aa8b0a2b790435bd0 Mon Sep 17 00:00:00 2001 From: Zhi Date: Wed, 11 Sep 2024 17:31:15 -0400 Subject: [PATCH 01/12] 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/12] 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/12] 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/12] 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/12] 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/12] 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/12] 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/12] 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/12] 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 eefc1680f8af0d8fd6bf239fcbbe1446354c6c32 Mon Sep 17 00:00:00 2001 From: Zhi Date: Wed, 18 Sep 2024 19:05:57 -0400 Subject: [PATCH 10/12] include pressure flux during transverse correction term --- Source/hydro/trans.cpp | 16 ++++++++++++++-- 1 file changed, 14 insertions(+), 2 deletions(-) diff --git a/Source/hydro/trans.cpp b/Source/hydro/trans.cpp index acbc4947df..bd12b522ec 100644 --- a/Source/hydro/trans.cpp +++ b/Source/hydro/trans.cpp @@ -292,9 +292,17 @@ Castro::actual_trans_single(const Box& bx, // NOLINT(readability-convert-member // // in cylindrical coords -- note that the p term is not // in a divergence for UMX in the x-direction, so there - // are no area factors. For this geometry, we do not + // are no area factors. + // + // Similarly for Spherical2D geometry, we have: + // d(rho u)/dt + d(rho u v)/(rdtheta) = -1/r^2 d(r^2 rho u u)/dr - dp/dr + // d(rho v)/dt + d(rho u v)/dr = -1/(r sin(theta)) d(sin(theta) rho v v)/dtheta - 1/r dp/dtheta + // + // For these non-cartesian geometries, we do not // include p in our definition of the flux in the - // x-direction, for we need to fix this now. + // x-direction for Cylindrical2D or both x- and y-direction for spherical 2D + // So we need to fix this now. + Real runewn = run - hdt * (area_t(ir,jr,kr) * flux_t(ir,jr,kr,UMX) - area_t(il,jl,kl) * flux_t(il,jl,kl,UMX)) * volinv; if (idir_t == 0 && !mom_flux_has_p(0, idir_t, coord)) { @@ -302,6 +310,10 @@ Castro::actual_trans_single(const Box& bx, // NOLINT(readability-convert-member } Real rvnewn = rvn - hdt * (area_t(ir,jr,kr) * flux_t(ir,jr,kr,UMY) - area_t(il,jl,kl) * flux_t(il,jl,kl,UMY)) * volinv; + if (idir_t == 1 && !mom_flux_has_p(1, idir_t, coord)) { + Real r = geom.ProbLo(0) + static_cast(il + 0.5_rt) * geom.CellSize(0); + rvnewn = rvnewn - cdtdy / r * (pgp - pgm); + } Real rwnewn = rwn - hdt * (area_t(ir,jr,kr) * flux_t(ir,jr,kr,UMZ) - area_t(il,jl,kl) * flux_t(il,jl,kl,UMZ)) * volinv; Real renewn = ren - hdt * (area_t(ir,jr,kr) * flux_t(ir,jr,kr,UEDEN) - From 791db2f52f88d4c2782fd7e68ca226ead6e3f879 Mon Sep 17 00:00:00 2001 From: Zhi Date: Wed, 18 Sep 2024 19:11:41 -0400 Subject: [PATCH 11/12] fix typo --- Source/hydro/trans.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Source/hydro/trans.cpp b/Source/hydro/trans.cpp index bd12b522ec..a21ae49208 100644 --- a/Source/hydro/trans.cpp +++ b/Source/hydro/trans.cpp @@ -312,7 +312,7 @@ Castro::actual_trans_single(const Box& bx, // NOLINT(readability-convert-member area_t(il,jl,kl) * flux_t(il,jl,kl,UMY)) * volinv; if (idir_t == 1 && !mom_flux_has_p(1, idir_t, coord)) { Real r = geom.ProbLo(0) + static_cast(il + 0.5_rt) * geom.CellSize(0); - rvnewn = rvnewn - cdtdy / r * (pgp - pgm); + rvnewn = rvnewn - cdtdx / r * (pgp - pgm); } Real rwnewn = rwn - hdt * (area_t(ir,jr,kr) * flux_t(ir,jr,kr,UMZ) - area_t(il,jl,kl) * flux_t(il,jl,kl,UMZ)) * volinv; From c4ebf03ddd0928438385819c9f4b63b5dd926c96 Mon Sep 17 00:00:00 2001 From: Zhi Date: Wed, 18 Sep 2024 20:25:15 -0400 Subject: [PATCH 12/12] fix gpu compilation --- Source/hydro/trans.cpp | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/Source/hydro/trans.cpp b/Source/hydro/trans.cpp index a21ae49208..18dd50c2e0 100644 --- a/Source/hydro/trans.cpp +++ b/Source/hydro/trans.cpp @@ -101,6 +101,8 @@ Castro::actual_trans_single(const Box& bx, // NOLINT(readability-convert-member #if AMREX_SPACEDIM == 2 int coord = geom.Coord(); + const auto dx = geom.CellSizeArray(); + const auto problo = geom.ProbLoArray(); #endif bool reset_density = transverse_reset_density; @@ -311,7 +313,7 @@ Castro::actual_trans_single(const Box& bx, // NOLINT(readability-convert-member Real rvnewn = rvn - hdt * (area_t(ir,jr,kr) * flux_t(ir,jr,kr,UMY) - area_t(il,jl,kl) * flux_t(il,jl,kl,UMY)) * volinv; if (idir_t == 1 && !mom_flux_has_p(1, idir_t, coord)) { - Real r = geom.ProbLo(0) + static_cast(il + 0.5_rt) * geom.CellSize(0); + Real r = problo[0] + static_cast(il + 0.5_rt) * dx[0]; rvnewn = rvnewn - cdtdx / r * (pgp - pgm); } Real rwnewn = rwn - hdt * (area_t(ir,jr,kr) * flux_t(ir,jr,kr,UMZ) -