From 885e9c0cb6f3a51a61a1b36776dc0deb9bc4b7a1 Mon Sep 17 00:00:00 2001 From: Zhi Chen <62574124+zhichen3@users.noreply.github.com> Date: Mon, 18 Nov 2024 21:16:50 -0500 Subject: [PATCH 1/4] fix indexing in geometric source term in ppm tracing for spherical coord (#2995) Supposed to load stencil in j instead of i when the normal direction is theta for spherical coord. --- Source/hydro/reconstruction.H | 69 +++++++++++++++++++++++++++-------- 1 file changed, 54 insertions(+), 15 deletions(-) diff --git a/Source/hydro/reconstruction.H b/Source/hydro/reconstruction.H index 8f791a1ad5..94a0f1b71f 100644 --- a/Source/hydro/reconstruction.H +++ b/Source/hydro/reconstruction.H @@ -103,11 +103,24 @@ add_geometric_rho_source(amrex::Array4 const& q_arr, // ncomp should be QU for idir == 0 and QV for idir == 1 - s[im2] += -dloga(i-2,j,k) * q_arr(i-2,j,k,QRHO) * q_arr(i-2,j,k,ncomp); - s[im1] += -dloga(i-1,j,k) * q_arr(i-1,j,k,QRHO) * q_arr(i-1,j,k,ncomp); - s[i0] += -dloga(i,j,k) * q_arr(i,j,k,QRHO) * q_arr(i,j,k,ncomp); - s[ip1] += -dloga(i+1,j,k) * q_arr(i+1,j,k,QRHO) * q_arr(i+1,j,k,ncomp); - s[ip2] += -dloga(i+2,j,k) * q_arr(i+2,j,k,QRHO) * q_arr(i+2,j,k,ncomp); + if (ncomp == QU) { + + s[im2] += -dloga(i-2,j,k) * q_arr(i-2,j,k,QRHO) * q_arr(i-2,j,k,ncomp); + s[im1] += -dloga(i-1,j,k) * q_arr(i-1,j,k,QRHO) * q_arr(i-1,j,k,ncomp); + s[i0] += -dloga(i,j,k) * q_arr(i,j,k,QRHO) * q_arr(i,j,k,ncomp); + s[ip1] += -dloga(i+1,j,k) * q_arr(i+1,j,k,QRHO) * q_arr(i+1,j,k,ncomp); + s[ip2] += -dloga(i+2,j,k) * q_arr(i+2,j,k,QRHO) * q_arr(i+2,j,k,ncomp); + + } else if (ncomp == QV) { + + s[im2] += -dloga(i,j-2,k) * q_arr(i,j-2,k,QRHO) * q_arr(i,j-2,k,ncomp); + s[im1] += -dloga(i,j-1,k) * q_arr(i,j-1,k,QRHO) * q_arr(i,j-1,k,ncomp); + s[i0] += -dloga(i,j,k) * q_arr(i,j,k,QRHO) * q_arr(i,j,k,ncomp); + s[ip1] += -dloga(i,j+1,k) * q_arr(i,j+1,k,QRHO) * q_arr(i,j+1,k,ncomp); + s[ip2] += -dloga(i,j+2,k) * q_arr(i,j+2,k,QRHO) * q_arr(i,j+2,k,ncomp); + + } + } AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE @@ -131,11 +144,24 @@ add_geometric_rhoe_source(amrex::Array4 const& q_arr, // ncomp should be QU for idir == 0 and QV for idir == 1 - s[im2] += -dloga(i-2,j,k) * (q_arr(i-2,j,k,QREINT) + q_arr(i-2,j,k,QPRES)) * q_arr(i-2,j,k,ncomp); - s[im1] += -dloga(i-1,j,k) * (q_arr(i-1,j,k,QREINT) + q_arr(i-1,j,k,QPRES)) * q_arr(i-1,j,k,ncomp); - s[i0] += -dloga(i,j,k) * (q_arr(i,j,k,QREINT) + q_arr(i,j,k,QPRES)) * q_arr(i,j,k,ncomp); - s[ip1] += -dloga(i+1,j,k) * (q_arr(i+1,j,k,QREINT) + q_arr(i+1,j,k,QPRES)) * q_arr(i+1,j,k,ncomp); - s[ip2] += -dloga(i+2,j,k) * (q_arr(i+2,j,k,QREINT) + q_arr(i+2,j,k,QPRES)) * q_arr(i+2,j,k,ncomp); + if (ncomp == QU) { + + s[im2] += -dloga(i-2,j,k) * (q_arr(i-2,j,k,QREINT) + q_arr(i-2,j,k,QPRES)) * q_arr(i-2,j,k,ncomp); + s[im1] += -dloga(i-1,j,k) * (q_arr(i-1,j,k,QREINT) + q_arr(i-1,j,k,QPRES)) * q_arr(i-1,j,k,ncomp); + s[i0] += -dloga(i,j,k) * (q_arr(i,j,k,QREINT) + q_arr(i,j,k,QPRES)) * q_arr(i,j,k,ncomp); + s[ip1] += -dloga(i+1,j,k) * (q_arr(i+1,j,k,QREINT) + q_arr(i+1,j,k,QPRES)) * q_arr(i+1,j,k,ncomp); + s[ip2] += -dloga(i+2,j,k) * (q_arr(i+2,j,k,QREINT) + q_arr(i+2,j,k,QPRES)) * q_arr(i+2,j,k,ncomp); + + } else if (ncomp == QV) { + + s[im2] += -dloga(i,j-2,k) * (q_arr(i,j-2,k,QREINT) + q_arr(i,j-2,k,QPRES)) * q_arr(i,j-2,k,ncomp); + s[im1] += -dloga(i,j-1,k) * (q_arr(i,j-1,k,QREINT) + q_arr(i,j-1,k,QPRES)) * q_arr(i,j-1,k,ncomp); + s[i0] += -dloga(i,j,k) * (q_arr(i,j,k,QREINT) + q_arr(i,j,k,QPRES)) * q_arr(i,j,k,ncomp); + s[ip1] += -dloga(i,j+1,k) * (q_arr(i,j+1,k,QREINT) + q_arr(i,j+1,k,QPRES)) * q_arr(i,j+1,k,ncomp); + s[ip2] += -dloga(i,j+2,k) * (q_arr(i,j+2,k,QREINT) + q_arr(i,j+2,k,QPRES)) * q_arr(i,j+2,k,ncomp); + + } + } AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE @@ -160,11 +186,24 @@ add_geometric_p_source(amrex::Array4 const& q_arr, // ncomp should be QU for idir == 0 and QV for idir == 1 - s[im2] += -dloga(i-2,j,k) * q_arr(i-2,j,k,QPRES) * qaux_arr(i-2,j,k,QGAMC) * q_arr(i-2,j,k,ncomp); - s[im1] += -dloga(i-1,j,k) * q_arr(i-1,j,k,QPRES) * qaux_arr(i-1,j,k,QGAMC) * q_arr(i-1,j,k,ncomp); - s[i0] += -dloga(i,j,k) * q_arr(i,j,k,QPRES) * qaux_arr(i,j,k,QGAMC) * q_arr(i,j,k,ncomp); - s[ip1] += -dloga(i+1,j,k) * q_arr(i+1,j,k,QPRES) * qaux_arr(i+1,j,k,QGAMC) * q_arr(i+1,j,k,ncomp); - s[ip2] += -dloga(i+2,j,k) * q_arr(i+2,j,k,QPRES) * qaux_arr(i+2,j,k,QGAMC) * q_arr(i+2,j,k,ncomp); + + if (ncomp == QU) { + + s[im2] += -dloga(i-2,j,k) * q_arr(i-2,j,k,QPRES) * qaux_arr(i-2,j,k,QGAMC) * q_arr(i-2,j,k,ncomp); + s[im1] += -dloga(i-1,j,k) * q_arr(i-1,j,k,QPRES) * qaux_arr(i-1,j,k,QGAMC) * q_arr(i-1,j,k,ncomp); + s[i0] += -dloga(i,j,k) * q_arr(i,j,k,QPRES) * qaux_arr(i,j,k,QGAMC) * q_arr(i,j,k,ncomp); + s[ip1] += -dloga(i+1,j,k) * q_arr(i+1,j,k,QPRES) * qaux_arr(i+1,j,k,QGAMC) * q_arr(i+1,j,k,ncomp); + s[ip2] += -dloga(i+2,j,k) * q_arr(i+2,j,k,QPRES) * qaux_arr(i+2,j,k,QGAMC) * q_arr(i+2,j,k,ncomp); + + } else if (ncomp == QV) { + + s[im2] += -dloga(i,j-2,k) * q_arr(i,j-2,k,QPRES) * qaux_arr(i,j-2,k,QGAMC) * q_arr(i,j-2,k,ncomp); + s[im1] += -dloga(i,j-1,k) * q_arr(i,j-1,k,QPRES) * qaux_arr(i,j-1,k,QGAMC) * q_arr(i,j-1,k,ncomp); + s[i0] += -dloga(i,j,k) * q_arr(i,j,k,QPRES) * qaux_arr(i,j,k,QGAMC) * q_arr(i,j,k,ncomp); + s[ip1] += -dloga(i,j+1,k) * q_arr(i,j+1,k,QPRES) * qaux_arr(i,j+1,k,QGAMC) * q_arr(i,j+1,k,ncomp); + s[ip2] += -dloga(i,j+2,k) * q_arr(i,j+2,k,QPRES) * qaux_arr(i,j+2,k,QGAMC) * q_arr(i,j+2,k,ncomp); + } + } From 1a22b17c56aa9eda31c0298d5af48e7dacf5b053 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Tue, 19 Nov 2024 20:41:28 -0500 Subject: [PATCH 2/4] fix field name in xrb_spherical plotting script (#2996) --- Exec/science/xrb_spherical/analysis/slice.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/Exec/science/xrb_spherical/analysis/slice.py b/Exec/science/xrb_spherical/analysis/slice.py index 7e55d3848f..b0ec8e40f0 100755 --- a/Exec/science/xrb_spherical/analysis/slice.py +++ b/Exec/science/xrb_spherical/analysis/slice.py @@ -66,12 +66,12 @@ def slice(fname:str, field:str, if field in ["x_velocity", "y_velocity", "z_velocity"]: sp.set_cmap(field, "coolwarm") elif field == "Temp": - sp.set_zlim(f, 5.e7, 2.5e9) - sp.set_cmap(f, "magma_r") + sp.set_zlim(field. 5.e7, 2.5e9) + sp.set_cmap(field, "magma_r") elif field == "enuc": - sp.set_zlim(f, 1.e18, 1.e20) + sp.set_zlim(field, 1.e18, 1.e20) elif field == "density": - sp.set_zlim(f, 1.e-3, 5.e8) + sp.set_zlim(field, 1.e-3, 5.e8) # sp.annotate_text((0.05, 0.05), f"{currentTime.in_cgs():8.5f} s") sp.save(f"{ds}_{loc}") From 8b3196a4b85f5aa5f4e2d19d72bbf40e2ce3e48c Mon Sep 17 00:00:00 2001 From: Zhi Chen <62574124+zhichen3@users.noreply.github.com> Date: Thu, 21 Nov 2024 07:36:03 -0500 Subject: [PATCH 3/4] make coverity happy (#2998) the original intention is to not execute the loop when spacedim=2, but coverity complains that its dead code. --- Source/driver/Castro.cpp | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/Source/driver/Castro.cpp b/Source/driver/Castro.cpp index 0133c68bfe..a3ceec8f8a 100644 --- a/Source/driver/Castro.cpp +++ b/Source/driver/Castro.cpp @@ -851,12 +851,13 @@ Castro::buildMetrics () dLogArea[dir].clear(); geom.GetDLogA(dLogArea[dir],grids,dmap,dir,NUM_GROW); } - for (int dir = AMREX_SPACEDIM; dir < 2; dir++) - { - dLogArea[dir].clear(); - dLogArea[dir].define(grids, dmap, 1, 0); - dLogArea[dir].setVal(0.0); - } + +#if (AMREX_SPACEDIM == 1) + dLogArea[1].clear(); + dLogArea[1].define(grids, dmap, 1, 0); + dLogArea[1].setVal(0.0); +#endif + #endif wall_time_start = 0.0; From f02444dc649462316a1b626f95d0d7198a0355f4 Mon Sep 17 00:00:00 2001 From: Zhi Chen <62574124+zhichen3@users.noreply.github.com> Date: Thu, 21 Nov 2024 10:51:41 -0500 Subject: [PATCH 4/4] add 2d spherical plm support (#2997) add 2d spherical plm support. Mainly adds the source terms. --- Source/hydro/trace_plm.cpp | 55 +++++++++++++++++++++++++++----------- 1 file changed, 39 insertions(+), 16 deletions(-) diff --git a/Source/hydro/trace_plm.cpp b/Source/hydro/trace_plm.cpp index 8b075b47eb..20e75ba02c 100644 --- a/Source/hydro/trace_plm.cpp +++ b/Source/hydro/trace_plm.cpp @@ -33,6 +33,8 @@ Castro::trace_plm(const Box& bx, const int idir, // vbx is the valid box (no ghost cells) const auto dx = geom.CellSizeArray(); + const auto problo = geom.ProbLoArray(); + const int coord = geom.Coord(); const int* lo_bc = phys_bc.lo(); const int* hi_bc = phys_bc.hi(); @@ -43,8 +45,6 @@ Castro::trace_plm(const Box& bx, const int idir, const auto domlo = geom.Domain().loVect3d(); const auto domhi = geom.Domain().hiVect3d(); - const Real dtdx = dt/dx[idir]; - auto vlo = vbx.loVect3d(); auto vhi = vbx.hiVect3d(); @@ -100,6 +100,16 @@ Castro::trace_plm(const Box& bx, const int idir, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { + Real dtdL = dt / dx[idir]; + Real dL = dx[idir]; + + // Want dt/(rdtheta) instead of dt/dtheta for 2d Spherical + if (coord == 2 && idir == 1) { + Real r = problo[0] + static_cast(i + 0.5_rt) * dx[0]; + dL = r * dx[1]; + dtdL = dt / dL; + } + bool lo_bc_test = lo_symm && ((idir == 0 && i == domlo[0]) || (idir == 1 && j == domlo[1]) || (idir == 2 && k == domlo[2])); @@ -166,7 +176,7 @@ Castro::trace_plm(const Box& bx, const int idir, load_stencil(srcQ, idir, i, j, k, QUN, src); Real dp = dq[IEIGN_P]; - pslope(trho, s, src, flat, lo_bc_test, hi_bc_test, dx[idir], dp); + pslope(trho, s, src, flat, lo_bc_test, hi_bc_test, dL, dp); dq[IEIGN_P] = dp; } @@ -185,7 +195,7 @@ Castro::trace_plm(const Box& bx, const int idir, // construct the right state on the i interface - Real ref_fac = 0.5_rt*(1.0_rt + dtdx*amrex::min(e[0], 0.0_rt)); + Real ref_fac = 0.5_rt*(1.0_rt + dtdL*amrex::min(e[0], 0.0_rt)); Real rho_ref = rho - ref_fac*dq[IEIGN_RHO]; Real un_ref = un - ref_fac*dq[IEIGN_UN]; Real ut_ref = ut - ref_fac*dq[IEIGN_UT]; @@ -195,8 +205,8 @@ Castro::trace_plm(const Box& bx, const int idir, // this is -(1/2) ( 1 + dt/dx lambda) (l . dq) r Real trace_fac0 = 0.0_rt; // FOURTH*dtdx*(e(1) - e(1))*(1.0_rt - sign(1.0_rt, e(1))) - Real trace_fac1 = 0.25_rt*dtdx*(e[0] - e[1])*(1.0_rt - std::copysign(1.0_rt, e[1])); - Real trace_fac2 = 0.25_rt*dtdx*(e[0] - e[2])*(1.0_rt - std::copysign(1.0_rt, e[2])); + Real trace_fac1 = 0.25_rt*dtdL*(e[0] - e[1])*(1.0_rt - std::copysign(1.0_rt, e[1])); + Real trace_fac2 = 0.25_rt*dtdL*(e[0] - e[2])*(1.0_rt - std::copysign(1.0_rt, e[2])); Real apright = trace_fac2*alphap; Real amright = trace_fac0*alpham; @@ -230,7 +240,7 @@ Castro::trace_plm(const Box& bx, const int idir, // now construct the left state on the i+1 interface - ref_fac = 0.5_rt*(1.0_rt - dtdx*amrex::max(e[2], 0.0_rt)); + ref_fac = 0.5_rt*(1.0_rt - dtdL*amrex::max(e[2], 0.0_rt)); rho_ref = rho + ref_fac*dq[IEIGN_RHO]; un_ref = un + ref_fac*dq[IEIGN_UN]; ut_ref = ut + ref_fac*dq[IEIGN_UT]; @@ -238,8 +248,8 @@ Castro::trace_plm(const Box& bx, const int idir, p_ref = p + ref_fac*dq[IEIGN_P]; rhoe_ref = rhoe + ref_fac*dq[IEIGN_RE]; - trace_fac0 = 0.25_rt*dtdx*(e[2] - e[0])*(1.0_rt + std::copysign(1.0_rt, e[0])); - trace_fac1 = 0.25_rt*dtdx*(e[2] - e[1])*(1.0_rt + std::copysign(1.0_rt, e[1])); + trace_fac0 = 0.25_rt*dtdL*(e[2] - e[0])*(1.0_rt + std::copysign(1.0_rt, e[0])); + trace_fac1 = 0.25_rt*dtdL*(e[2] - e[1])*(1.0_rt + std::copysign(1.0_rt, e[1])); trace_fac2 = 0.0_rt; // FOURTH*dtdx*(e(3) - e(3))*(1.0_rt + sign(1.0_rt, e(3))) Real apleft = trace_fac2*alphap; @@ -301,29 +311,42 @@ Castro::trace_plm(const Box& bx, const int idir, } #if (AMREX_SPACEDIM < 3) - // geometry source terms -- these only apply to the x-states - if (idir == 0 && dloga(i,j,k) != 0.0_rt) { - Real courn = dtdx*(cc + std::abs(un)); + // geometry source terms + // these only apply to the x-states for cylindrical and spherical + // or y-states for spherical + + if (dloga(i,j,k) != 0.0_rt) { + Real courn = dtdL*(cc + std::abs(un)); Real eta = (1.0_rt-courn)/(cc*dt*std::abs(dloga(i,j,k))); Real dlogatmp = amrex::min(eta, 1.0_rt)*dloga(i,j,k); Real sourcr = -0.5_rt*dt*rho*dlogatmp*un; Real sourcp = sourcr*csq; Real source = sourcp*enth; - if (i <= vhi[0]) { + if (idir == 0 && i <= vhi[0]) { qm(i+1,j,k,QRHO) += sourcr; qm(i+1,j,k,QRHO) = amrex::max(qm(i+1,j,k,QRHO), lsmall_dens); qm(i+1,j,k,QPRES) += sourcp; qm(i+1,j,k,QREINT) += source; } - if (i >= vlo[0]) { + if (idir == 1 && j <= vhi[1]) { + qm(i,j+1,k,QRHO) += sourcr; + qm(i,j+1,k,QRHO) = amrex::max(qm(i,j+1,k,QRHO), lsmall_dens); + qm(i,j+1,k,QPRES) += sourcp; + qm(i,j+1,k,QREINT) += source; + } + + if ((idir == 0 && i >= vlo[0]) || + (idir == 1 && j >= vlo[1])) { qp(i,j,k,QRHO) += sourcr; qp(i,j,k,QRHO) = amrex::max(qp(i,j,k,QRHO), lsmall_dens); qp(i,j,k,QPRES) += sourcp; qp(i,j,k,QREINT) += source; } + } + #endif for (int ipassive = 0; ipassive < npassive; ipassive++) { @@ -340,12 +363,12 @@ Castro::trace_plm(const Box& bx, const int idir, (idir == 1 && j >= vlo[1]) || (idir == 2 && k >= vlo[2])) { - Real spzero = un >= 0.0_rt ? -1.0_rt : un*dtdx; + Real spzero = un >= 0.0_rt ? -1.0_rt : un*dtdL; qp(i,j,k,n) = s[i0] + 0.5_rt*(-1.0_rt - spzero)*dX; } // Left state - Real spzero = un >= 0.0_rt ? un*dtdx : 1.0_rt; + Real spzero = un >= 0.0_rt ? un*dtdL : 1.0_rt; Real acmpleft = 0.5_rt*(1.0_rt - spzero )*dX; if (idir == 0 && i <= vhi[0]) {