From 3395e2ea13048ef9a6f2564aa8b0a2b790435bd0 Mon Sep 17 00:00:00 2001 From: Zhi Date: Wed, 11 Sep 2024 17:31:15 -0400 Subject: [PATCH 1/9] 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 2/9] 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 3/9] 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 4/9] 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 5/9] 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 6/9] 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 7/9] 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 8/9] 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 9/9] 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())) {