From 7710d239a22e199a68262ee237c49c68c4451b3d Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Mon, 18 Apr 2022 14:46:46 -0400 Subject: [PATCH 01/24] update an action to some later versions --- .github/workflows/gh-pages.yml | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/.github/workflows/gh-pages.yml b/.github/workflows/gh-pages.yml index f9b480b9b7..9a92fdab2b 100644 --- a/.github/workflows/gh-pages.yml +++ b/.github/workflows/gh-pages.yml @@ -8,18 +8,18 @@ on: jobs: deploy: - runs-on: ubuntu-18.04 + runs-on: ubuntu-latest steps: - - uses: actions/checkout@v2 + - uses: actions/checkout@v3 - - name: Install pandoc and doxygen + - name: Install pandoc and doxygen run: | sudo apt install pandoc doxygen - name: Setup Python - uses: actions/setup-python@v2 + uses: actions/setup-python@v3 with: - python-version: '3.8' + python-version: '3.9' - name: Upgrade pip run: | From a9517e7e3d3b9c581caf45f4f7037f4a0eca341a Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Fri, 22 Apr 2022 13:48:16 -0400 Subject: [PATCH 02/24] use one-sided stencils when necessary --- Source/hydro/ppm.H | 161 ++++++++++++++++++++++++++-------- Source/hydro/reconstruction.H | 18 ++-- 2 files changed, 135 insertions(+), 44 deletions(-) diff --git a/Source/hydro/ppm.H b/Source/hydro/ppm.H index 2509c7cc3b..ec046775c5 100644 --- a/Source/hydro/ppm.H +++ b/Source/hydro/ppm.H @@ -54,68 +54,151 @@ check_trace_source(Array4 const& srcQ, const int idir, AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void ppm_reconstruct(const Real* s, + const int i, const int j, const int k, const int idir, + const bool lo_bc_test, const bool hi_bc_test, + GpuArray const domlo, GpuArray const domhi, const Real flatn, Real& sm, Real& sp) { + // first we compute s_{i-1/2} -- the left interface value for zone i - // Compute van Leer slopes + if (lo_bc_test && ((idir == 0 && i == domlo[0]) || + (idir == 1 && j == domlo[1]) || + (idir == 2 && k == domlo[2]))) { - Real dsl = 2.0_rt * (s[im1] - s[im2]); - Real dsr = 2.0_rt * (s[i0] - s[im1]); + // use a stencil for when the current zone is on the left physical + // boundary. Then the left interface is on the physical boundary, + // MC Eq. 21 + sm = (1.0_rt/12.0_rt)*(25.0_rt*s[i0] - 23.0_rt*s[ip1] + 13.0_rt*s[ip2] - 3.0_rt*s[ip3]); - Real dsvl_l = 0.0_rt; - if (dsl*dsr > 0.0_rt) { - Real dsc = 0.5_rt * (s[i0] - s[im2]); - dsvl_l = std::copysign(1.0_rt, dsc) * amrex::min(std::abs(dsc),amrex::min(std::abs(dsl),std::abs(dsr))); - } + } else if (lo_bc_test && ((idir == 0 && i == domlo[0]+1) || + (idir == 1 && j == domlo[1]+1) || + (idir == 2 && k == domlo[2]+1))) { - dsl = 2.0_rt * (s[i0] - s[im1]); - dsr = 2.0_rt * (s[ip1] - s[i0]); + // use a stencil for when the current zone is one zone away from + // the left physical boundary, and then the left interface is one + // zone away from the boundary, MC Eq. 22 + sm = (1.0_rt/12.0_rt)*(3.0_rt*s[im1] + 13.0_rt*s[i0] - 5.0_rt*s[ip1] + s[ip2]); - Real dsvl_r = 0.0_rt; - if (dsl*dsr > 0.0_rt) { - Real dsc = 0.5_rt * (s[ip1] - s[im1]); - dsvl_r = std::copysign(1.0_rt, dsc) * amrex::min(std::abs(dsc),amrex::min(std::abs(dsl),std::abs(dsr))); - } + // Make sure sedge lies in between adjacent cell-centered values - // Interpolate s to edges + sm = amrex::max(sm, amrex::min(s[i0], s[im1])); + sm = amrex::min(sm, amrex::max(s[i0], s[im1])); - sm = 0.5_rt * (s[i0] + s[im1]) - (1.0_rt/6.0_rt) * (dsvl_r - dsvl_l); + } else if (hi_bc_test && ((idir == 0 && i == domhi[0]) || + (idir == 1 && j == domhi[1]) || + (idir == 2 && k == domhi[2]))) { - // Make sure sedge lies in between adjacent cell-centered values + // use a stencil for when the current zone is on the right physical boundary + // then the left interface is one zone away from the physical boundary + sm = (1.0_rt/12.0_rt)*(3.0_rt*s[i0] + 13.0_rt*s[im1] - 5.0_rt*s[im2] + s[im3]); - sm = amrex::max(sm, amrex::min(s[i0], s[im1])); - sm = amrex::min(sm, amrex::max(s[i0], s[im1])); + // Make sure sedge lies in between adjacent cell-centered values + sm = amrex::max(sm, amrex::min(s[i0], s[im1])); + sm = amrex::min(sm, amrex::max(s[i0], s[im1])); - // Compute van Leer slopes + } else { - dsl = 2.0_rt * (s[i0] - s[im1]); - dsr = 2.0_rt * (s[ip1] - s[i0]); + // Compute van Leer slopes - 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), amrex::min(std::abs(dsl),std::abs(dsr))); - } + Real dsl = 2.0_rt * (s[im1] - s[im2]); + Real dsr = 2.0_rt * (s[i0] - s[im1]); + + Real dsvl_l = 0.0_rt; + if (dsl*dsr > 0.0_rt) { + Real dsc = 0.5_rt * (s[i0] - s[im2]); + dsvl_l = std::copysign(1.0_rt, dsc) * amrex::min(std::abs(dsc),amrex::min(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; + if (dsl*dsr > 0.0_rt) { + Real dsc = 0.5_rt * (s[ip1] - s[im1]); + dsvl_r = std::copysign(1.0_rt, dsc) * amrex::min(std::abs(dsc),amrex::min(std::abs(dsl),std::abs(dsr))); + } - dsl = 2.0_rt * (s[ip1] - s[i0]); - dsr = 2.0_rt * (s[ip2] - s[ip1]); + // Interpolate s to edges + + sm = 0.5_rt * (s[i0] + s[im1]) - (1.0_rt/6.0_rt) * (dsvl_r - dsvl_l); + + // Make sure sedge lies in between adjacent cell-centered values + + sm = amrex::max(sm, amrex::min(s[i0], s[im1])); + sm = amrex::min(sm, amrex::max(s[i0], s[im1])); - 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), amrex::min(std::abs(dsl),std::abs(dsr))); } - // Interpolate s to edges + // now we compute s_{i+1/2} -- the right interface value for zone i + + if (lo_bc_test && ((idir == 0 && i == domlo[0]) || + (idir == 1 && j == domlo[1]) || + (idir == 2 && k == domlo[2]))) { + + // use a stencil for when the current zone is on the left physical + // boundary. Then the right interface is one zone away from the + // physical boundary, + sp = (1.0_rt/12.0_rt)*(3.0_rt*s[i0] + 13.0_rt*s[ip1] - 5.0_rt*s[ip2] + s[ip3]); + + // Make sure sedge lies in between adjacent cell-centered values + + sp = amrex::max(sp, amrex::min(s[ip1], s[i0])); + sp = amrex::min(sp, amrex::max(s[ip1], s[i0])); - sp = 0.5_rt * (s[ip1] + s[i0]) - (1.0_rt/6.0_rt) * (dsvl_r - dsvl_l); + } else if (hi_bc_test && ((idir == 0 && i == domhi[0]) || + (idir == 1 && j == domhi[1]) || + (idir == 2 && k == domhi[2]))) { - // Make sure sedge lies in between adjacent cell-centered values + // use a stencil for when the current zone is on the right physical boundary + // then the right interface is on the physical boundary + sp = (1.0_rt/12.0_rt)*(25.0_rt*s[i0] - 23.0_rt*s[im1] + 13.0_rt*s[im2] - 3.0_rt*s[im3]); - sp = amrex::max(sp, amrex::min(s[ip1], s[i0])); - sp = amrex::min(sp, amrex::max(s[ip1], s[i0])); + } else if (hi_bc_test && ((idir == 0 && i == domhi[0]-1) || + (idir == 1 && j == domhi[1]-1) || + (idir == 2 && k == domhi[2]-1))) { + // use a stencil for when the current zone is one zone away from + // the right physical boundary, and then the right interface is one + // zone away from the boundary + sp = (1.0_rt/12.0_rt)*(3.0_rt*s[ip1] + 13.0_rt*s[i0] - 5.0_rt*s[im1] + s[im2]); + + // Make sure sedge lies in between adjacent cell-centered values + + sp = amrex::max(sp, amrex::min(s[ip1], s[i0])); + sp = amrex::min(sp, amrex::max(s[ip1], s[i0])); + + } else { + + // Compute van Leer slopes + + Real dsl = 2.0_rt * (s[i0] - s[im1]); + Real dsr = 2.0_rt * (s[ip1] - s[i0]); + + Real 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), amrex::min(std::abs(dsl),std::abs(dsr))); + } + + dsl = 2.0_rt * (s[ip1] - s[i0]); + dsr = 2.0_rt * (s[ip2] - s[ip1]); + + Real dsvl_r = 0.0_rt; + 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), amrex::min(std::abs(dsl),std::abs(dsr))); + } + + // Interpolate s to edges + + sp = 0.5_rt * (s[ip1] + s[i0]) - (1.0_rt/6.0_rt) * (dsvl_r - dsvl_l); + + // Make sure sedge lies in between adjacent cell-centered values + + sp = amrex::max(sp, amrex::min(s[ip1], s[i0])); + sp = amrex::min(sp, amrex::max(s[ip1], s[i0])); + } // Flatten the parabola diff --git a/Source/hydro/reconstruction.H b/Source/hydro/reconstruction.H index 30a348c4fa..0f15f3b80b 100644 --- a/Source/hydro/reconstruction.H +++ b/Source/hydro/reconstruction.H @@ -1,11 +1,13 @@ #ifndef CASTRO_RECONSTRUCTION_H #define CASTRO_RECONSTRUCTION_H -constexpr int im2 = 0; -constexpr int im1 = 1; -constexpr int i0 = 2; -constexpr int ip1 = 3; -constexpr int ip2 = 4; +constexpr int im3 = 0; +constexpr int im2 = 1; +constexpr int im1 = 2; +constexpr int i0 = 3; +constexpr int ip1 = 4; +constexpr int ip2 = 5; +constexpr int ip3 = 6; AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void @@ -14,25 +16,31 @@ load_stencil(Array4 const& q_arr, const int idir, Real* s) { if (idir == 0) { + s[im3] = q_arr(i-3,j,k,ncomp); s[im2] = q_arr(i-2,j,k,ncomp); s[im1] = q_arr(i-1,j,k,ncomp); s[i0] = q_arr(i,j,k,ncomp); s[ip1] = q_arr(i+1,j,k,ncomp); s[ip2] = q_arr(i+2,j,k,ncomp); + s[ip3] = q_arr(i+3,j,k,ncomp); } else if (idir == 1) { + s[im3] = q_arr(i,j-3,k,ncomp); s[im2] = q_arr(i,j-2,k,ncomp); s[im1] = q_arr(i,j-1,k,ncomp); s[i0] = q_arr(i,j,k,ncomp); s[ip1] = q_arr(i,j+1,k,ncomp); s[ip2] = q_arr(i,j+2,k,ncomp); + s[ip3] = q_arr(i,j+3,k,ncomp); } else { + s[im3] = q_arr(i,j,k-3,ncomp); s[im2] = q_arr(i,j,k-2,ncomp); s[im1] = q_arr(i,j,k-1,ncomp); s[i0] = q_arr(i,j,k,ncomp); s[ip1] = q_arr(i,j,k+1,ncomp); s[ip2] = q_arr(i,j,k+2,ncomp); + s[ip3] = q_arr(i,j,k+3,ncomp); } From 9665212114dfafa066df4e2f62c1e65cc69d6442 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Sun, 24 Apr 2022 17:32:49 -0400 Subject: [PATCH 03/24] finish syncing --- Source/hydro/Castro_ctu.cpp | 146 ++++++++++++++++++++++++++++++++++++ Source/hydro/Castro_mol.cpp | 47 +++++------- Source/hydro/ppm.H | 2 +- Source/hydro/trace_ppm.cpp | 42 +++++++---- 4 files changed, 192 insertions(+), 45 deletions(-) diff --git a/Source/hydro/Castro_ctu.cpp b/Source/hydro/Castro_ctu.cpp index 79d29fa7da..2fa2950eaf 100644 --- a/Source/hydro/Castro_ctu.cpp +++ b/Source/hydro/Castro_ctu.cpp @@ -146,6 +146,152 @@ Castro::ctu_ppm_states(const Box& bx, const Box& vbx, #endif } + + // special care for reflecting BCs + const int* lo_bc = phys_bc.lo(); + const int* hi_bc = phys_bc.hi(); + + const auto domlo = geom.Domain().loVect3d(); + const auto domhi = geom.Domain().hiVect3d(); + + bool lo_bc_test = lo_bc[idir] == Symmetry; + bool hi_bc_test = hi_bc[idir] == Symmetry; + + // we have to do this after the loops above, since here we will + // consider interfaces, not zones + + if (idir == 0) { + if (lo_bc_test) { + + amrex::ParallelFor(bx, + [=] AMREX_GPU_HOST_DEVICE (int i, int j, int k) + { + + // reset the left state at domlo(0) if needed -- it is outside the domain + if (i == domlo[0]) { + for (int n = 0; n < NQ; n++) { + if (n == QU) { + qxm(i,j,k,QU) = -qxp(i,j,k,QU); + } else { + qxm(i,j,k,n) = qxp(i,j,k,n); + } + } + } + }); + + } + + if (hi_bc_test) { + + amrex::ParallelFor(bx, + [=] AMREX_GPU_HOST_DEVICE (int i, int j, int k) + { + + // reset the right state at domhi(0)+1 if needed -- it is outside the domain + if (i == domhi[0]+1) { + for (int n = 0; n < NQ; n++) { + if (n == QU) { + qxp(i,j,k,QU) = -qxm(i,j,k,QU); + } else { + qxp(i,j,k,n) = qxm(i,j,k,n); + } + } + } + }); + + } + +#if AMREX_SPACEDIM >= 2 + } else if (idir == 1) { + + if (lo_bc_test) { + + amrex::ParallelFor(bx, + [=] AMREX_GPU_HOST_DEVICE (int i, int j, int k) + { + + // reset the left state at domlo(0) if needed -- it is outside the domain + if (j == domlo[1]) { + for (int n = 0; n < NQ; n++) { + if (n == QV) { + qym(i,j,k,QV) = -qyp(i,j,k,QV); + } else { + qym(i,j,k,n) = qyp(i,j,k,n); + } + } + } + }); + + } + + if (hi_bc_test) { + + amrex::ParallelFor(bx, + [=] AMREX_GPU_HOST_DEVICE (int i, int j, int k) + { + + // reset the right state at domhi(0)+1 if needed -- it is outside the domain + if (j == domhi[1]+1) { + for (int n = 0; n < NQ; n++) { + if (n == QV) { + qyp(i,j,k,QV) = -qym(i,j,k,QV); + } else { + qyp(i,j,k,n) = qym(i,j,k,n); + } + } + } + }); + + } + +#endif +#if AMREX_SPACEDIM == 3 + } else { + + if (lo_bc_test) { + + amrex::ParallelFor(bx, + [=] AMREX_GPU_HOST_DEVICE (int i, int j, int k) + { + + // reset the left state at domlo(0) if needed -- it is outside the domain + if (k == domlo[2]) { + for (int n = 0; n < NQ; n++) { + if (n == QW) { + qzm(i,j,k,QW) = -qzp(i,j,k,QW); + } else { + qzm(i,j,k,n) = qzp(i,j,k,n); + } + } + } + }); + + } + + if (hi_bc_test) { + + amrex::ParallelFor(bx, + [=] AMREX_GPU_HOST_DEVICE (int i, int j, int k) + { + + // reset the right state at domhi(0)+1 if needed -- it is outside the domain + if (k == domhi[2]+1) { + for (int n = 0; n < NQ; n++) { + if (n == QW) { + qzp(i,j,k,QW) = -qzm(i,j,k,QW); + } else { + qzp(i,j,k,n) = qzm(i,j,k,n); + } + } + } + }); + + } +#endif + + } + + } } diff --git a/Source/hydro/Castro_mol.cpp b/Source/hydro/Castro_mol.cpp index e24d802aa8..fec951c3b7 100644 --- a/Source/hydro/Castro_mol.cpp +++ b/Source/hydro/Castro_mol.cpp @@ -48,7 +48,7 @@ Castro::mol_plm_reconstruct(const Box& bx, (idir == 1 && j == domhi[1]) || (idir == 2 && k == domhi[2])); - Real s[5]; + Real s[7]; Real flat = flatn_arr(i,j,k); if (idir == 0) { @@ -85,11 +85,11 @@ Castro::mol_plm_reconstruct(const Box& bx, [=] AMREX_GPU_HOST_DEVICE (int i, int j, int k) { - Real s[5]; + Real s[7]; Real flat = flatn_arr(i,j,k); - Real trho[5]; - Real src[5]; + Real trho[7]; + Real src[7]; bool lo_bc_test = lo_symm && ((idir == 0 && i == domlo[0]) || (idir == 1 && j == domlo[1]) || @@ -357,39 +357,30 @@ Castro::mol_ppm_reconstruct(const Box& bx, Array4 const& qm, Array4 const& qp) { + // special care for reflecting BCs + const int* lo_bc = phys_bc.lo(); + const int* hi_bc = phys_bc.hi(); + + const auto domlo = geom.Domain().loVect3d(); + const auto domhi = geom.Domain().hiVect3d(); + + bool lo_bc_test = lo_bc[idir] == Symmetry; + bool hi_bc_test = hi_bc[idir] == Symmetry; + amrex::ParallelFor(bx, NQ, [=] AMREX_GPU_HOST_DEVICE (int i, int j, int k, int n) { - Real s[5]; + Real s[7]; Real flat = flatn_arr(i,j,k); Real sm; Real sp; - if (idir == 0) { - s[im2] = q_arr(i-2,j,k,n); - s[im1] = q_arr(i-1,j,k,n); - s[i0] = q_arr(i,j,k,n); - s[ip1] = q_arr(i+1,j,k,n); - s[ip2] = q_arr(i+2,j,k,n); - - } else if (idir == 1) { - s[im2] = q_arr(i,j-2,k,n); - s[im1] = q_arr(i,j-1,k,n); - s[i0] = q_arr(i,j,k,n); - s[ip1] = q_arr(i,j+1,k,n); - s[ip2] = q_arr(i,j+2,k,n); - - } else { - s[im2] = q_arr(i,j,k-2,n); - s[im1] = q_arr(i,j,k-1,n); - s[i0] = q_arr(i,j,k,n); - s[ip1] = q_arr(i,j,k+1,n); - s[ip2] = q_arr(i,j,k+2,n); - - } + load_stencil(q_arr, idir, i, j, k, n, s); - ppm_reconstruct(s, flat, sm, sp); + ppm_reconstruct(s, i, j, k, idir, + lo_bc_test, hi_bc_test, domlo, domhi, + flat, sm, sp); if (idir == 0) { // right state at i-1/2 diff --git a/Source/hydro/ppm.H b/Source/hydro/ppm.H index ec046775c5..be6c6700ee 100644 --- a/Source/hydro/ppm.H +++ b/Source/hydro/ppm.H @@ -46,7 +46,7 @@ check_trace_source(Array4 const& srcQ, const int idir, /// Compute the coefficients of a parabolic reconstruction of the data in a zone. /// This uses the standard PPM limiters described in Colella & Woodward (1984) /// -/// @param s Real[5] the state to be reconstructed in zones i-2, i-1, i, i+1, i+2 +/// @param s Real[7] the state to be reconstructed in zones i-3, i-2, i-1, i, i+1, i+2, i+3 /// @param flatn The flattening coefficient /// @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 diff --git a/Source/hydro/trace_ppm.cpp b/Source/hydro/trace_ppm.cpp index 175d9d58de..942a733531 100644 --- a/Source/hydro/trace_ppm.cpp +++ b/Source/hydro/trace_ppm.cpp @@ -134,6 +134,16 @@ Castro::trace_ppm(const Box& bx, Real lsmall_dens = small_dens; Real lsmall_pres = small_pres; + // special care for the reflecting BCs + const int* lo_bc = phyc_bc.lo(); + const int* hi_bc = phys_bc.hi(); + + const auto domlo = geom.Domain().loVect3d(); + const auto domhi = geom.Domain().hiVect3d(); + + bool lo_bc_test = lo_bc[idir] == Symmetry; + bool hi_bc_test = hi_bc[idir] == Symmetry; + // Trace to left and right edges using upwind PPM amrex::ParallelFor(bx, [=] AMREX_GPU_HOST_DEVICE (int i, int j, int k) @@ -150,7 +160,7 @@ Castro::trace_ppm(const Box& bx, // do the parabolic reconstruction and compute the // integrals under the characteristic waves - Real s[5]; + Real s[7]; Real flat = flatn(i,j,k); Real sm; Real sp; @@ -162,7 +172,7 @@ Castro::trace_ppm(const Box& bx, Real Im_rho[3]; load_stencil(q_arr, idir, i, j, k, QRHO, s); - ppm_reconstruct(s, flat, sm, sp); + ppm_reconstruct(s, i, j, k, idir, lo_bc_test, hi_bc_test, domlo, domhi, flat, sm, sp); ppm_int_profile(sm, sp, s[i0], un, cc, dtdx, Ip_rho, Im_rho); // reconstruct normal velocity @@ -173,7 +183,7 @@ Castro::trace_ppm(const Box& bx, Real Im_un_2; load_stencil(q_arr, idir, i, j, k, QUN, s); - ppm_reconstruct(s, flat, sm, sp); + ppm_reconstruct(s, i, j, k, idir, lo_bc_test, hi_bc_test, domlo, domhi, flat, sm, sp); ppm_int_profile_single(sm, sp, s[i0], un-cc, dtdx, Ip_un_0, Im_un_0); ppm_int_profile_single(sm, sp, s[i0], un+cc, dtdx, Ip_un_2, Im_un_2); @@ -183,7 +193,7 @@ Castro::trace_ppm(const Box& bx, Real Im_p[3]; load_stencil(q_arr, idir, i, j, k, QPRES, s); - ppm_reconstruct(s, flat, sm, sp); + ppm_reconstruct(s, i, j, k, idir, lo_bc_test, hi_bc_test, domlo, domhi, flat, sm, sp); ppm_int_profile(sm, sp, s[i0], un, cc, dtdx, Ip_p, Im_p); // reconstruct rho e @@ -192,7 +202,7 @@ Castro::trace_ppm(const Box& bx, Real Im_rhoe[3]; load_stencil(q_arr, idir, i, j, k, QREINT, s); - ppm_reconstruct(s, flat, sm, sp); + ppm_reconstruct(s, i, j, k, idir, lo_bc_test, hi_bc_test, domlo, domhi, flat, sm, sp); ppm_int_profile(sm, sp, s[i0], un, cc, dtdx, Ip_rhoe, Im_rhoe); // reconstruct transverse velocities @@ -203,11 +213,11 @@ Castro::trace_ppm(const Box& bx, Real Im_utt_1; load_stencil(q_arr, idir, i, j, k, QUT, s); - ppm_reconstruct(s, flat, sm, sp); + ppm_reconstruct(s, i, j, k, idir, lo_bc_test, hi_bc_test, domlo, domhi, flat, sm, sp); ppm_int_profile_single(sm, sp, s[i0], un, dtdx, Ip_ut_1, Im_ut_1); load_stencil(q_arr, idir, i, j, k, QUTT, s); - ppm_reconstruct(s, flat, sm, sp); + ppm_reconstruct(s, i, j, k, idir, lo_bc_test, hi_bc_test, domlo, domhi, flat, sm, sp); ppm_int_profile_single(sm, sp, s[i0], un, dtdx, Ip_utt_1, Im_utt_1); // gamma_c @@ -218,7 +228,7 @@ Castro::trace_ppm(const Box& bx, Real Im_gc_2; load_stencil(qaux_arr, idir, i, j, k, QGAMC, s); - ppm_reconstruct(s, flat, sm, sp); + ppm_reconstruct(s, i, j, k, idir, lo_bc_test, hi_bc_test, domlo, domhi, flat, sm, sp); ppm_int_profile_single(sm, sp, s[i0], un-cc, dtdx, Ip_gc_0, Im_gc_0); ppm_int_profile_single(sm, sp, s[i0], un+cc, dtdx, Ip_gc_2, Im_gc_2); @@ -239,7 +249,7 @@ Castro::trace_ppm(const Box& bx, if (do_trace) { load_stencil(srcQ, idir, i, j, k, QRHO, s); - ppm_reconstruct(s, flat, sm, sp); + ppm_reconstruct(s, i, j, k, idir, lo_bc_test, hi_bc_test, domlo, domhi, flat, sm, sp); ppm_int_profile(sm, sp, s[i0], un, cc, dtdx, Ip_src_rho, Im_src_rho); } @@ -258,7 +268,7 @@ Castro::trace_ppm(const Box& bx, if (do_trace) { load_stencil(srcQ, idir, i, j, k, QUN, s); - ppm_reconstruct(s, flat, sm, sp); + ppm_reconstruct(s, i, j, k, idir, lo_bc_test, hi_bc_test, domlo, domhi, flat, sm, sp); ppm_int_profile_single(sm, sp, s[i0], un-cc, dtdx, Ip_src_un_0, Im_src_un_0); ppm_int_profile_single(sm, sp, s[i0], un+cc, dtdx, Ip_src_un_2, Im_src_un_2); } @@ -276,7 +286,7 @@ Castro::trace_ppm(const Box& bx, if (do_trace) { load_stencil(srcQ, idir, i, j, k, QPRES, s); - ppm_reconstruct(s, flat, sm, sp); + ppm_reconstruct(s, i, j, k, idir, lo_bc_test, hi_bc_test, domlo, domhi, flat, sm, sp); ppm_int_profile(sm, sp, s[i0], un, cc, dtdx, Ip_src_p, Im_src_p); } @@ -293,7 +303,7 @@ Castro::trace_ppm(const Box& bx, if (do_trace) { load_stencil(srcQ, idir, i, j, k, QREINT, s); - ppm_reconstruct(s, flat, sm, sp); + ppm_reconstruct(s, i, j, k, idir, lo_bc_test, hi_bc_test, domlo, domhi, flat, sm, sp); ppm_int_profile(sm, sp, s[i0], un, cc, dtdx, Ip_src_rhoe, Im_src_rhoe); } @@ -310,7 +320,7 @@ Castro::trace_ppm(const Box& bx, if (do_trace) { load_stencil(srcQ, idir, i, j, k, QUT, s); - ppm_reconstruct(s, flat, sm, sp); + ppm_reconstruct(s, i, j, k, idir, lo_bc_test, hi_bc_test, domlo, domhi, flat, sm, sp); ppm_int_profile_single(sm, sp, s[i0], un, dtdx, Ip_src_ut_1, Im_src_ut_1); } @@ -325,7 +335,7 @@ Castro::trace_ppm(const Box& bx, if (do_trace) { load_stencil(srcQ, idir, i, j, k, QUTT, s); - ppm_reconstruct(s, flat, sm, sp); + ppm_reconstruct(s, i, j, k, idir, lo_bc_test, hi_bc_test, domlo, domhi, flat, sm, sp); ppm_int_profile_single(sm, sp, s[i0], un, dtdx, Ip_src_utt_1, Im_src_utt_1); } @@ -347,13 +357,13 @@ Castro::trace_ppm(const Box& bx, load_stencil(q_arr, idir, i, j, k, n, s); - ppm_reconstruct(s, flat, sm, sp); + ppm_reconstruct(s, i, j, k, idir, lo_bc_test, hi_bc_test, domlo, domhi, flat, sm, sp); ppm_int_profile_single(sm, sp, s[i0], un, dtdx, Ip_passive, Im_passive); #ifdef PRIM_SPECIES_HAVE_SOURCES // if we turned this on, don't bother to check if it source is non-zero -- just trace load_stencil(srcQ, idir, i, j, k, n, s); - ppm_reconstruct(s, flat, sm, sp); + ppm_reconstruct(s, i, j, k, idir, lo_bc_test, hi_bc_test, domlo, domhi, flat, sm, sp); ppm_int_profile_single(sm, sp, s[i0], un, dtdx, Ip_src_passive, Im_src_passive); #endif From 51ac6c8058cd8a046dde77d9e86a7548e575b8d6 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Sun, 24 Apr 2022 17:35:53 -0400 Subject: [PATCH 04/24] fix compilation --- Source/hydro/trace_ppm.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Source/hydro/trace_ppm.cpp b/Source/hydro/trace_ppm.cpp index 942a733531..52e35d299e 100644 --- a/Source/hydro/trace_ppm.cpp +++ b/Source/hydro/trace_ppm.cpp @@ -135,7 +135,7 @@ Castro::trace_ppm(const Box& bx, Real lsmall_pres = small_pres; // special care for the reflecting BCs - const int* lo_bc = phyc_bc.lo(); + const int* lo_bc = phys_bc.lo(); const int* hi_bc = phys_bc.hi(); const auto domlo = geom.Domain().loVect3d(); From e22e9eb974a951c4a265c1fb8d3743ade7d1008e Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Sun, 24 Apr 2022 20:21:41 -0400 Subject: [PATCH 05/24] some diag --- Source/hydro/Castro_ctu.cpp | 17 +++++++++++++++++ 1 file changed, 17 insertions(+) diff --git a/Source/hydro/Castro_ctu.cpp b/Source/hydro/Castro_ctu.cpp index 2fa2950eaf..81712c1b57 100644 --- a/Source/hydro/Castro_ctu.cpp +++ b/Source/hydro/Castro_ctu.cpp @@ -157,9 +157,12 @@ Castro::ctu_ppm_states(const Box& bx, const Box& vbx, bool lo_bc_test = lo_bc[idir] == Symmetry; bool hi_bc_test = hi_bc[idir] == Symmetry; + const auto dx = geom.CellSizeArray(); + // we have to do this after the loops above, since here we will // consider interfaces, not zones + if (idir == 0) { if (lo_bc_test) { @@ -172,11 +175,25 @@ Castro::ctu_ppm_states(const Box& bx, const Box& vbx, for (int n = 0; n < NQ; n++) { if (n == QU) { qxm(i,j,k,QU) = -qxp(i,j,k,QU); +// } else if (n == QPRES) { +// // integrate from the right zone edge to this zone edge and reset both the left and right states here +// qxp(i,j,k,QPRES) = qxm(i+1,j,k,QPRES) - q_arr(i,j,k,QRHO) * srcQ(i,j,k,QU) * dx[0]; +// qxm(i,j,k,QPRES) = qxp(i,j,k,QPRES); + } else { qxm(i,j,k,n) = qxp(i,j,k,n); } } } + + // if (i == domlo[0]) { + // for (int ii = 0; ii < 4; ++ii) { + // // qxp q qxm + // std::cout << ii << " " << qxp(i+ii,j,k,QPRES) << " " << q_arr(i+ii,j,k,QPRES) << " " << qxm(i+ii+1,j,k,QPRES) << std::endl; + // } + // amrex::Error("stop"); + // } + }); } From baeb2dc8b92e4c20df9c8c0b3d5ae78343d9a5dc Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Tue, 26 Apr 2022 12:39:08 -0400 Subject: [PATCH 06/24] add some protection on the parabola values --- Source/hydro/trace_ppm.cpp | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/Source/hydro/trace_ppm.cpp b/Source/hydro/trace_ppm.cpp index 52e35d299e..eecb0a9035 100644 --- a/Source/hydro/trace_ppm.cpp +++ b/Source/hydro/trace_ppm.cpp @@ -173,6 +173,8 @@ Castro::trace_ppm(const Box& bx, load_stencil(q_arr, idir, i, j, k, QRHO, s); ppm_reconstruct(s, i, j, k, idir, lo_bc_test, hi_bc_test, domlo, domhi, flat, sm, sp); + sm = amrex::max(lsmall_dens, sm); + sp = amrex::max(lsmall_dens, sp); ppm_int_profile(sm, sp, s[i0], un, cc, dtdx, Ip_rho, Im_rho); // reconstruct normal velocity @@ -194,6 +196,8 @@ Castro::trace_ppm(const Box& bx, load_stencil(q_arr, idir, i, j, k, QPRES, s); ppm_reconstruct(s, i, j, k, idir, lo_bc_test, hi_bc_test, domlo, domhi, flat, sm, sp); + sm = amrex::max(lsmall_pres, sm); + sp = amrex::max(lsmall_pres, sp); ppm_int_profile(sm, sp, s[i0], un, cc, dtdx, Ip_p, Im_p); // reconstruct rho e @@ -203,6 +207,8 @@ Castro::trace_ppm(const Box& bx, load_stencil(q_arr, idir, i, j, k, QREINT, s); ppm_reconstruct(s, i, j, k, idir, lo_bc_test, hi_bc_test, domlo, domhi, flat, sm, sp); + sm = amrex::max(lsmall_pres, sm); + sp = amrex::max(lsmall_pres, sp); ppm_int_profile(sm, sp, s[i0], un, cc, dtdx, Ip_rhoe, Im_rhoe); // reconstruct transverse velocities From 83a5120662a05238864f98e9e70a80f865a15033 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Tue, 26 Apr 2022 13:35:57 -0400 Subject: [PATCH 07/24] add more reflect enforcement --- Source/hydro/Castro_ctu.cpp | 169 ++---------------------------- Source/hydro/Castro_ctu_hydro.cpp | 24 +++-- Source/hydro/Castro_hydro.H | 14 +++ Source/hydro/advection_util.cpp | 152 +++++++++++++++++++++++++++ 4 files changed, 186 insertions(+), 173 deletions(-) diff --git a/Source/hydro/Castro_ctu.cpp b/Source/hydro/Castro_ctu.cpp index 81712c1b57..3c3540051a 100644 --- a/Source/hydro/Castro_ctu.cpp +++ b/Source/hydro/Castro_ctu.cpp @@ -115,7 +115,7 @@ Castro::ctu_ppm_states(const Box& bx, const Box& vbx, for (int idir = 0; idir < AMREX_SPACEDIM; idir++) { if (idir == 0) { - trace_ppm(bx, + trace_ppm(bx, idir, q_arr, qaux_arr, srcQ, flatn, qxm, qxp, @@ -124,6 +124,8 @@ Castro::ctu_ppm_states(const Box& bx, const Box& vbx, #endif vbx, dt); + enforce_reflection_states(bx, idir, qxm, qxp); + #if AMREX_SPACEDIM >= 2 } else if (idir == 1) { trace_ppm(bx, @@ -134,6 +136,8 @@ Castro::ctu_ppm_states(const Box& bx, const Box& vbx, dloga, #endif vbx, dt); + + enforce_reflection_states(bx, idir, qym, qyp); #endif #if AMREX_SPACEDIM == 3 @@ -144,171 +148,10 @@ Castro::ctu_ppm_states(const Box& bx, const Box& vbx, qzm, qzp, vbx, dt); + enforce_reflection_states(bx, idir, qzm, qzp); #endif } - // special care for reflecting BCs - const int* lo_bc = phys_bc.lo(); - const int* hi_bc = phys_bc.hi(); - - const auto domlo = geom.Domain().loVect3d(); - const auto domhi = geom.Domain().hiVect3d(); - - bool lo_bc_test = lo_bc[idir] == Symmetry; - bool hi_bc_test = hi_bc[idir] == Symmetry; - - const auto dx = geom.CellSizeArray(); - - // we have to do this after the loops above, since here we will - // consider interfaces, not zones - - - if (idir == 0) { - if (lo_bc_test) { - - amrex::ParallelFor(bx, - [=] AMREX_GPU_HOST_DEVICE (int i, int j, int k) - { - - // reset the left state at domlo(0) if needed -- it is outside the domain - if (i == domlo[0]) { - for (int n = 0; n < NQ; n++) { - if (n == QU) { - qxm(i,j,k,QU) = -qxp(i,j,k,QU); -// } else if (n == QPRES) { -// // integrate from the right zone edge to this zone edge and reset both the left and right states here -// qxp(i,j,k,QPRES) = qxm(i+1,j,k,QPRES) - q_arr(i,j,k,QRHO) * srcQ(i,j,k,QU) * dx[0]; -// qxm(i,j,k,QPRES) = qxp(i,j,k,QPRES); - - } else { - qxm(i,j,k,n) = qxp(i,j,k,n); - } - } - } - - // if (i == domlo[0]) { - // for (int ii = 0; ii < 4; ++ii) { - // // qxp q qxm - // std::cout << ii << " " << qxp(i+ii,j,k,QPRES) << " " << q_arr(i+ii,j,k,QPRES) << " " << qxm(i+ii+1,j,k,QPRES) << std::endl; - // } - // amrex::Error("stop"); - // } - - }); - - } - - if (hi_bc_test) { - - amrex::ParallelFor(bx, - [=] AMREX_GPU_HOST_DEVICE (int i, int j, int k) - { - - // reset the right state at domhi(0)+1 if needed -- it is outside the domain - if (i == domhi[0]+1) { - for (int n = 0; n < NQ; n++) { - if (n == QU) { - qxp(i,j,k,QU) = -qxm(i,j,k,QU); - } else { - qxp(i,j,k,n) = qxm(i,j,k,n); - } - } - } - }); - - } - -#if AMREX_SPACEDIM >= 2 - } else if (idir == 1) { - - if (lo_bc_test) { - - amrex::ParallelFor(bx, - [=] AMREX_GPU_HOST_DEVICE (int i, int j, int k) - { - - // reset the left state at domlo(0) if needed -- it is outside the domain - if (j == domlo[1]) { - for (int n = 0; n < NQ; n++) { - if (n == QV) { - qym(i,j,k,QV) = -qyp(i,j,k,QV); - } else { - qym(i,j,k,n) = qyp(i,j,k,n); - } - } - } - }); - - } - - if (hi_bc_test) { - - amrex::ParallelFor(bx, - [=] AMREX_GPU_HOST_DEVICE (int i, int j, int k) - { - - // reset the right state at domhi(0)+1 if needed -- it is outside the domain - if (j == domhi[1]+1) { - for (int n = 0; n < NQ; n++) { - if (n == QV) { - qyp(i,j,k,QV) = -qym(i,j,k,QV); - } else { - qyp(i,j,k,n) = qym(i,j,k,n); - } - } - } - }); - - } - -#endif -#if AMREX_SPACEDIM == 3 - } else { - - if (lo_bc_test) { - - amrex::ParallelFor(bx, - [=] AMREX_GPU_HOST_DEVICE (int i, int j, int k) - { - - // reset the left state at domlo(0) if needed -- it is outside the domain - if (k == domlo[2]) { - for (int n = 0; n < NQ; n++) { - if (n == QW) { - qzm(i,j,k,QW) = -qzp(i,j,k,QW); - } else { - qzm(i,j,k,n) = qzp(i,j,k,n); - } - } - } - }); - - } - - if (hi_bc_test) { - - amrex::ParallelFor(bx, - [=] AMREX_GPU_HOST_DEVICE (int i, int j, int k) - { - - // reset the right state at domhi(0)+1 if needed -- it is outside the domain - if (k == domhi[2]+1) { - for (int n = 0; n < NQ; n++) { - if (n == QW) { - qzp(i,j,k,QW) = -qzm(i,j,k,QW); - } else { - qzp(i,j,k,n) = qzm(i,j,k,n); - } - } - } - }); - - } -#endif - - } - - } } diff --git a/Source/hydro/Castro_ctu_hydro.cpp b/Source/hydro/Castro_ctu_hydro.cpp index 7e190e35cf..b7c8fcaae1 100644 --- a/Source/hydro/Castro_ctu_hydro.cpp +++ b/Source/hydro/Castro_ctu_hydro.cpp @@ -617,9 +617,9 @@ Castro::construct_ctu_hydro_source(Real time, Real dt) vol_arr, hdt, hdtdy); - reset_edge_state_thermo(xbx, ql.array()); + reset_edge_state_thermo(xbx, ql_arr); - reset_edge_state_thermo(xbx, qr.array()); + reset_edge_state_thermo(xbx, qr_arr); // solve the final Riemann problem axross the x-interfaces @@ -629,6 +629,8 @@ Castro::construct_ctu_hydro_source(Real time, Real dt) #endif + enforce_reflection_states(xbx, 0, ql_arr, qr_arr); + cmpflx_plus_godunov(xbx, ql_arr, qr_arr, flux0_arr, @@ -659,9 +661,9 @@ Castro::construct_ctu_hydro_source(Real time, Real dt) vol_arr, hdt, hdtdx); - reset_edge_state_thermo(ybx, ql.array()); + reset_edge_state_thermo(ybx, ql_arr); - reset_edge_state_thermo(ybx, qr.array()); + reset_edge_state_thermo(ybx, qr_arr); // solve the final Riemann problem axross the y-interfaces @@ -672,6 +674,8 @@ Castro::construct_ctu_hydro_source(Real time, Real dt) #endif + enforce_reflection_states(ybx, 1, ql_arr, qr_arr); + cmpflx_plus_godunov(ybx, ql_arr, qr_arr, flux1_arr, @@ -991,9 +995,9 @@ Castro::construct_ctu_hydro_source(Real time, Real dt) qgdnvtmp2_arr, hdtdy, hdtdz); - reset_edge_state_thermo(xbx, ql.array()); + reset_edge_state_thermo(xbx, ql_arr); - reset_edge_state_thermo(xbx, qr.array()); + reset_edge_state_thermo(xbx, qr_arr); #ifdef SIMPLIFIED_SDC add_sdc_source_to_states(xbx, 0, dt, @@ -1069,9 +1073,9 @@ Castro::construct_ctu_hydro_source(Real time, Real dt) qgdnvtmp1_arr, hdtdx, hdtdz); - reset_edge_state_thermo(ybx, ql.array()); + reset_edge_state_thermo(ybx, ql_arr); - reset_edge_state_thermo(ybx, qr.array()); + reset_edge_state_thermo(ybx, qr_arr); #ifdef SIMPLIFIED_SDC add_sdc_source_to_states(ybx, 1, dt, @@ -1149,9 +1153,9 @@ Castro::construct_ctu_hydro_source(Real time, Real dt) qgdnvtmp2_arr, hdtdx, hdtdy); - reset_edge_state_thermo(zbx, ql.array()); + reset_edge_state_thermo(zbx, ql_arr); - reset_edge_state_thermo(zbx, qr.array()); + reset_edge_state_thermo(zbx, qr_arr); #ifdef SIMPLIFIED_SDC add_sdc_source_to_states(zbx, 2, dt, diff --git a/Source/hydro/Castro_hydro.H b/Source/hydro/Castro_hydro.H index c549e5a27d..36b9526977 100644 --- a/Source/hydro/Castro_hydro.H +++ b/Source/hydro/Castro_hydro.H @@ -118,6 +118,20 @@ amrex::Array4 const& q_arr, amrex::Array4 const& qaux_arr); +/// +/// for left / right interface states on a reflecting boundary, reflect the state +/// inside the domain to the outside +/// +/// @param bx the box to operate over +/// @param idir coordinate direction +/// @param qm left states on the interfaces +/// @param qp right states on the interfaces +/// + void enforce_reflection_states(const amrex::Box& bx, + const int idir, + amrex::Array4 const& qm, + amrex::Array4 const& qp); + /// /// compute the flattening coefficient. This is 0 if we are in a shock and /// 1 if we are not applying any flattening diff --git a/Source/hydro/advection_util.cpp b/Source/hydro/advection_util.cpp index addf019560..65d0860188 100644 --- a/Source/hydro/advection_util.cpp +++ b/Source/hydro/advection_util.cpp @@ -986,3 +986,155 @@ Castro::do_enforce_minimum_density(const Box& bx, } }); } + + + +void +Castro::enforce_reflection_states(const Box& bx, + const int idir, + Array4 const& qm, + Array4 const& qp) { + + const int* lo_bc = phys_bc.lo(); + const int* hi_bc = phys_bc.hi(); + + const auto domlo = geom.Domain().loVect3d(); + const auto domhi = geom.Domain().hiVect3d(); + + bool lo_bc_test = lo_bc[idir] == Symmetry; + bool hi_bc_test = hi_bc[idir] == Symmetry; + + const auto dx = geom.CellSizeArray(); + + if (idir == 0) { + if (lo_bc_test) { + + amrex::ParallelFor(bx, + [=] AMREX_GPU_HOST_DEVICE (int i, int j, int k) + { + + // reset the left state at domlo(0) if needed -- it is outside the domain + if (i == domlo[0]) { + for (int n = 0; n < NQ; n++) { + if (n == QU) { + qm(i,j,k,QU) = -qp(i,j,k,QU); + } else { + qm(i,j,k,n) = qp(i,j,k,n); + } + } + } + + }); + } + + if (hi_bc_test) { + + amrex::ParallelFor(bx, + [=] AMREX_GPU_HOST_DEVICE (int i, int j, int k) + { + + // reset the right state at domhi(0)+1 if needed -- it is outside the domain + if (i == domhi[0]+1) { + for (int n = 0; n < NQ; n++) { + if (n == QU) { + qp(i,j,k,QU) = -qm(i,j,k,QU); + } else { + qp(i,j,k,n) = qm(i,j,k,n); + } + } + } + }); + + } + +#if AMREX_SPACEDIM >= 2 + } else if (idir == 1) { + + if (lo_bc_test) { + + amrex::ParallelFor(bx, + [=] AMREX_GPU_HOST_DEVICE (int i, int j, int k) + { + + // reset the left state at domlo(0) if needed -- it is outside the domain + if (j == domlo[1]) { + for (int n = 0; n < NQ; n++) { + if (n == QV) { + qm(i,j,k,QV) = -qp(i,j,k,QV); + } else { + qm(i,j,k,n) = qp(i,j,k,n); + } + } + } + }); + + } + + if (hi_bc_test) { + + amrex::ParallelFor(bx, + [=] AMREX_GPU_HOST_DEVICE (int i, int j, int k) + { + + // reset the right state at domhi(0)+1 if needed -- it is outside the domain + if (j == domhi[1]+1) { + for (int n = 0; n < NQ; n++) { + if (n == QV) { + qp(i,j,k,QV) = -qm(i,j,k,QV); + } else { + qp(i,j,k,n) = qm(i,j,k,n); + } + } + } + }); + + } + +#endif +#if AMREX_SPACEDIM == 3 + } else { + + if (lo_bc_test) { + + amrex::ParallelFor(bx, + [=] AMREX_GPU_HOST_DEVICE (int i, int j, int k) + { + + // reset the left state at domlo(0) if needed -- it is outside the domain + if (k == domlo[2]) { + for (int n = 0; n < NQ; n++) { + if (n == QW) { + qm(i,j,k,QW) = -qp(i,j,k,QW); + } else { + qm(i,j,k,n) = qp(i,j,k,n); + } + } + } + }); + + } + + if (hi_bc_test) { + + amrex::ParallelFor(bx, + [=] AMREX_GPU_HOST_DEVICE (int i, int j, int k) + { + + // reset the right state at domhi(0)+1 if needed -- it is outside the domain + if (k == domhi[2]+1) { + for (int n = 0; n < NQ; n++) { + if (n == QW) { + qp(i,j,k,QW) = -qm(i,j,k,QW); + } else { + qp(i,j,k,n) = qm(i,j,k,n); + } + } + } + }); + + } +#endif + + } + +} From 7fae1f7335a67f27d2969f1f9e1dda939f7ab78c Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Thu, 19 May 2022 10:47:47 -0400 Subject: [PATCH 08/24] some merge fixes --- Source/hydro/Castro_ctu.cpp | 1 - Source/hydro/Castro_ctu_hydro.cpp | 4 ++-- Source/hydro/Castro_hydro.H | 13 ------------- 3 files changed, 2 insertions(+), 16 deletions(-) diff --git a/Source/hydro/Castro_ctu.cpp b/Source/hydro/Castro_ctu.cpp index d7dfb98383..5a75152322 100644 --- a/Source/hydro/Castro_ctu.cpp +++ b/Source/hydro/Castro_ctu.cpp @@ -151,7 +151,6 @@ Castro::ctu_ppm_states(const Box& bx, const Box& vbx, enforce_reflect_states(bx, 2, qzm, qzp); #endif } - } } diff --git a/Source/hydro/Castro_ctu_hydro.cpp b/Source/hydro/Castro_ctu_hydro.cpp index 156ed8125f..092a3d53fa 100644 --- a/Source/hydro/Castro_ctu_hydro.cpp +++ b/Source/hydro/Castro_ctu_hydro.cpp @@ -583,7 +583,7 @@ Castro::construct_ctu_hydro_source(Real time, Real dt) #endif - enforce_reflection_states(xbx, 0, ql_arr, qr_arr); + enforce_reflect_states(xbx, 0, ql_arr, qr_arr); cmpflx_plus_godunov(xbx, ql_arr, qr_arr, @@ -629,7 +629,7 @@ Castro::construct_ctu_hydro_source(Real time, Real dt) #endif #endif - enforce_reflection_states(ybx, 1, ql_arr, qr_arr); + enforce_reflect_states(ybx, 1, ql_arr, qr_arr); cmpflx_plus_godunov(ybx, ql_arr, qr_arr, diff --git a/Source/hydro/Castro_hydro.H b/Source/hydro/Castro_hydro.H index 72fd4282ce..60424017a1 100644 --- a/Source/hydro/Castro_hydro.H +++ b/Source/hydro/Castro_hydro.H @@ -118,19 +118,6 @@ amrex::Array4 const& q_arr, amrex::Array4 const& qaux_arr); -/// -/// for left / right interface states on a reflecting boundary, reflect the state -/// inside the domain to the outside -/// -/// @param bx the box to operate over -/// @param idir coordinate direction -/// @param qm left states on the interfaces -/// @param qp right states on the interfaces -/// - void enforce_reflection_states(const amrex::Box& bx, - const int idir, - amrex::Array4 const& qm, - amrex::Array4 const& qp); /// /// compute the flattening coefficient. This is 0 if we are in a shock and From 66e93504b08bd92097cce88b53bc78cb23eed3b4 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Thu, 19 May 2022 10:48:15 -0400 Subject: [PATCH 09/24] more merge fixes --- Source/hydro/Castro_hydro.H | 1 - 1 file changed, 1 deletion(-) diff --git a/Source/hydro/Castro_hydro.H b/Source/hydro/Castro_hydro.H index 60424017a1..fdf3348747 100644 --- a/Source/hydro/Castro_hydro.H +++ b/Source/hydro/Castro_hydro.H @@ -118,7 +118,6 @@ amrex::Array4 const& q_arr, amrex::Array4 const& qaux_arr); - /// /// compute the flattening coefficient. This is 0 if we are in a shock and /// 1 if we are not applying any flattening From 59c14244517decf32c346d341a95006c8086b73c Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Fri, 27 May 2022 12:02:14 -0400 Subject: [PATCH 10/24] more merging --- Source/hydro/ppm.H | 4 ++-- Source/hydro/trace_ppm.cpp | 10 ---------- 2 files changed, 2 insertions(+), 12 deletions(-) diff --git a/Source/hydro/ppm.H b/Source/hydro/ppm.H index cf6dc9263b..892d09d3ff 100644 --- a/Source/hydro/ppm.H +++ b/Source/hydro/ppm.H @@ -56,7 +56,7 @@ void ppm_reconstruct(const Real* s, const int i, const int j, const int k, const int idir, const bool lo_bc_test, const bool hi_bc_test, - GpuArray& const domlo, GpuArray& const domhi, + const GpuArray& domlo, const GpuArray& domhi, const Real flatn, Real& sm, Real& sp) { // first we compute s_{i-1/2} -- the left interface value for zone i @@ -244,7 +244,7 @@ void ppm_reconstruct_pressure(const Real* rho, const Real* p, const Real* src, const int i, const int j, const int k, const int idir, bool lo_bc_test, bool hi_bc_test, - GpuArray& const domlo, GpuArray& const domhi, + const GpuArray& domlo, const GpuArray& domhi, const Real dx, const Real flatn, Real& sm, Real& sp) { diff --git a/Source/hydro/trace_ppm.cpp b/Source/hydro/trace_ppm.cpp index 61ca2fdd9d..6533617d3c 100644 --- a/Source/hydro/trace_ppm.cpp +++ b/Source/hydro/trace_ppm.cpp @@ -143,16 +143,6 @@ Castro::trace_ppm(const Box& bx, Real lsmall_dens = small_dens; Real lsmall_pres = small_pres; - // special care for the reflecting BCs - const int* lo_bc = phys_bc.lo(); - const int* hi_bc = phys_bc.hi(); - - const auto domlo = geom.Domain().loVect3d(); - const auto domhi = geom.Domain().hiVect3d(); - - bool lo_bc_test = lo_bc[idir] == Symmetry; - bool hi_bc_test = hi_bc[idir] == Symmetry; - // Trace to left and right edges using upwind PPM amrex::ParallelFor(bx, [=] AMREX_GPU_HOST_DEVICE (int i, int j, int k) From cb56c001b0d7bf1899ae8be646f36ecefc99180a Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Fri, 27 May 2022 12:14:58 -0400 Subject: [PATCH 11/24] more fixing --- Source/hydro/trace_ppm.cpp | 35 ++++++++++++++--------------------- 1 file changed, 14 insertions(+), 21 deletions(-) diff --git a/Source/hydro/trace_ppm.cpp b/Source/hydro/trace_ppm.cpp index 6533617d3c..911af22aff 100644 --- a/Source/hydro/trace_ppm.cpp +++ b/Source/hydro/trace_ppm.cpp @@ -149,13 +149,6 @@ Castro::trace_ppm(const Box& bx, { - bool lo_bc_test = lo_symm && ((idir == 0 && i == domlo[0]) || - (idir == 1 && j == domlo[1]) || - (idir == 2 && k == domlo[2])); - - bool hi_bc_test = hi_symm && ((idir == 0 && i == domhi[0]) || - (idir == 1 && j == domhi[1]) || - (idir == 2 && k == domhi[2])); Real cc = qaux_arr(i,j,k,QC); @@ -202,7 +195,7 @@ Castro::trace_ppm(const Box& bx, Real Im_rho[3]; load_stencil(q_arr, idir, i, j, k, QRHO, s); - ppm_reconstruct(s, i, j, k, idir, lo_bc_test, hi_bc_test, domlo, domhi, flat, sm, sp); + ppm_reconstruct(s, i, j, k, idir, lo_symm, hi_symm, domlo, domhi, flat, sm, sp); sm = amrex::max(lsmall_dens, sm); sp = amrex::max(lsmall_dens, sp); ppm_int_profile(sm, sp, s[i0], un, cc, dtdx, Ip_rho, Im_rho); @@ -215,7 +208,7 @@ Castro::trace_ppm(const Box& bx, Real Im_un_2; load_stencil(q_arr, idir, i, j, k, QUN, s); - ppm_reconstruct(s, i, j, k, idir, lo_bc_test, hi_bc_test, domlo, domhi, flat, sm, sp); + ppm_reconstruct(s, i, j, k, idir, lo_symm, hi_symm, domlo, domhi, flat, sm, sp); ppm_int_profile_single(sm, sp, s[i0], un-cc, dtdx, Ip_un_0, Im_un_0); ppm_int_profile_single(sm, sp, s[i0], un+cc, dtdx, Ip_un_2, Im_un_2); @@ -233,7 +226,7 @@ Castro::trace_ppm(const Box& bx, load_stencil(srcQ, idir, i, j, k, QUN, src); ppm_reconstruct_pressure(trho, s, src, - i, j, k, idir, lo_bc_test, hi_bc_test, domlo, domhi, + i, j, k, idir, lo_symm, hi_symm, domlo, domhi, dx[idir], flat, sm, sp); sm = amrex::max(lsmall_pres, sm); sp = amrex::max(lsmall_pres, sp); @@ -246,7 +239,7 @@ Castro::trace_ppm(const Box& bx, Real Im_rhoe[3]; load_stencil(q_arr, idir, i, j, k, QREINT, s); - ppm_reconstruct(s, i, j, k, idir, lo_bc_test, hi_bc_test, domlo, domhi, flat, sm, sp); + ppm_reconstruct(s, i, j, k, idir, lo_symm, hi_symm, domlo, domhi, flat, sm, sp); sm = amrex::max(lsmall_pres, sm); sp = amrex::max(lsmall_pres, sp); ppm_int_profile(sm, sp, s[i0], un, cc, dtdx, Ip_rhoe, Im_rhoe); @@ -259,11 +252,11 @@ Castro::trace_ppm(const Box& bx, Real Im_utt_1; load_stencil(q_arr, idir, i, j, k, QUT, s); - ppm_reconstruct(s, i, j, k, idir, lo_bc_test, hi_bc_test, domlo, domhi, flat, sm, sp); + ppm_reconstruct(s, i, j, k, idir, lo_symm, hi_symm, domlo, domhi, flat, sm, sp); ppm_int_profile_single(sm, sp, s[i0], un, dtdx, Ip_ut_1, Im_ut_1); load_stencil(q_arr, idir, i, j, k, QUTT, s); - ppm_reconstruct(s, i, j, k, idir, lo_bc_test, hi_bc_test, domlo, domhi, flat, sm, sp); + ppm_reconstruct(s, i, j, k, idir, lo_symm, hi_symm, domlo, domhi, flat, sm, sp); ppm_int_profile_single(sm, sp, s[i0], un, dtdx, Ip_utt_1, Im_utt_1); // gamma_c @@ -274,7 +267,7 @@ Castro::trace_ppm(const Box& bx, Real Im_gc_2; load_stencil(qaux_arr, idir, i, j, k, QGAMC, s); - ppm_reconstruct(s, i, j, k, idir, lo_bc_test, hi_bc_test, domlo, domhi, flat, sm, sp); + ppm_reconstruct(s, i, j, k, idir, lo_symm, hi_symm, domlo, domhi, flat, sm, sp); ppm_int_profile_single(sm, sp, s[i0], un-cc, dtdx, Ip_gc_0, Im_gc_0); ppm_int_profile_single(sm, sp, s[i0], un+cc, dtdx, Ip_gc_2, Im_gc_2); @@ -295,7 +288,7 @@ Castro::trace_ppm(const Box& bx, if (do_trace) { load_stencil(srcQ, idir, i, j, k, QRHO, s); - ppm_reconstruct(s, i, j, k, idir, lo_bc_test, hi_bc_test, domlo, domhi, flat, sm, sp); + ppm_reconstruct(s, i, j, k, idir, lo_symm, hi_symm, domlo, domhi, flat, sm, sp); ppm_int_profile(sm, sp, s[i0], un, cc, dtdx, Ip_src_rho, Im_src_rho); } @@ -314,7 +307,7 @@ Castro::trace_ppm(const Box& bx, if (do_trace) { load_stencil(srcQ, idir, i, j, k, QUN, s); - ppm_reconstruct(s, i, j, k, idir, lo_bc_test, hi_bc_test, domlo, domhi, flat, sm, sp); + ppm_reconstruct(s, i, j, k, idir, lo_symm, hi_symm, domlo, domhi, flat, sm, sp); ppm_int_profile_single(sm, sp, s[i0], un-cc, dtdx, Ip_src_un_0, Im_src_un_0); ppm_int_profile_single(sm, sp, s[i0], un+cc, dtdx, Ip_src_un_2, Im_src_un_2); } @@ -332,7 +325,7 @@ Castro::trace_ppm(const Box& bx, if (do_trace) { load_stencil(srcQ, idir, i, j, k, QPRES, s); - ppm_reconstruct(s, i, j, k, idir, lo_bc_test, hi_bc_test, domlo, domhi, flat, sm, sp); + ppm_reconstruct(s, i, j, k, idir, lo_symm, hi_symm, domlo, domhi, flat, sm, sp); ppm_int_profile(sm, sp, s[i0], un, cc, dtdx, Ip_src_p, Im_src_p); } @@ -349,7 +342,7 @@ Castro::trace_ppm(const Box& bx, if (do_trace) { load_stencil(srcQ, idir, i, j, k, QREINT, s); - ppm_reconstruct(s, i, j, k, idir, lo_bc_test, hi_bc_test, domlo, domhi, flat, sm, sp); + ppm_reconstruct(s, i, j, k, idir, lo_symm, hi_symm, domlo, domhi, flat, sm, sp); ppm_int_profile(sm, sp, s[i0], un, cc, dtdx, Ip_src_rhoe, Im_src_rhoe); } @@ -366,7 +359,7 @@ Castro::trace_ppm(const Box& bx, if (do_trace) { load_stencil(srcQ, idir, i, j, k, QUT, s); - ppm_reconstruct(s, i, j, k, idir, lo_bc_test, hi_bc_test, domlo, domhi, flat, sm, sp); + ppm_reconstruct(s, i, j, k, idir, lo_symm, hi_symm, domlo, domhi, flat, sm, sp); ppm_int_profile_single(sm, sp, s[i0], un, dtdx, Ip_src_ut_1, Im_src_ut_1); } @@ -381,7 +374,7 @@ Castro::trace_ppm(const Box& bx, if (do_trace) { load_stencil(srcQ, idir, i, j, k, QUTT, s); - ppm_reconstruct(s, i, j, k, idir, lo_bc_test, hi_bc_test, domlo, domhi, flat, sm, sp); + ppm_reconstruct(s, i, j, k, idir, lo_symm, hi_symm, domlo, domhi, flat, sm, sp); ppm_int_profile_single(sm, sp, s[i0], un, dtdx, Ip_src_utt_1, Im_src_utt_1); } @@ -399,7 +392,7 @@ Castro::trace_ppm(const Box& bx, load_stencil(q_arr, idir, i, j, k, n, s); - ppm_reconstruct(s, i, j, k, idir, lo_bc_test, hi_bc_test, domlo, domhi, flat, sm, sp); + ppm_reconstruct(s, i, j, k, idir, lo_symm, hi_symm, domlo, domhi, flat, sm, sp); ppm_int_profile_single(sm, sp, s[i0], un, dtdx, Ip_passive, Im_passive); // Plus state on face i From 00e05ecc78c4b49a4a5cb06586f129e58a097cde Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Fri, 27 May 2022 12:23:54 -0400 Subject: [PATCH 12/24] add reflect --- Source/hydro/Castro_ctu_hydro.cpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/Source/hydro/Castro_ctu_hydro.cpp b/Source/hydro/Castro_ctu_hydro.cpp index 092a3d53fa..8b4760db50 100644 --- a/Source/hydro/Castro_ctu_hydro.cpp +++ b/Source/hydro/Castro_ctu_hydro.cpp @@ -460,6 +460,8 @@ Castro::construct_ctu_hydro_source(Real time, Real dt) shk_arr, 0, false); + enforce_reflect_states(xbx, 0, qxm_arr, qxp_arr); + #endif // 1-d From 9995481cd43380a3f4a4449140dc1b75691fd3b8 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Mon, 25 Jul 2022 11:23:30 -0400 Subject: [PATCH 13/24] fix one-sided + pslope --- Source/hydro/ppm.H | 12 ++++++++++-- 1 file changed, 10 insertions(+), 2 deletions(-) diff --git a/Source/hydro/ppm.H b/Source/hydro/ppm.H index 0925e7e210..575f2cc0f7 100644 --- a/Source/hydro/ppm.H +++ b/Source/hydro/ppm.H @@ -249,7 +249,7 @@ ppm_reconstruct_pslope(const Real* rho, const Real* p, const Real* src, const Real flatn, Real& sm, Real& sp) { - Real tp[5]; + Real tp[7]; // compute the hydrostatic pressure in each zone center starting with i0 @@ -257,16 +257,20 @@ ppm_reconstruct_pslope(const Real* rho, const Real* p, const Real* src, 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]); + Real pp3_hse = pp2_hse + 0.25_rt*dx * (rho[ip2] + rho[ip3]) * (src[ip2] + src[ip3]); 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]); + Real pm3_hse = pm2_hse - 0.25_rt*dx * (rho[im2] + rho[im3]) * (src[im2] + src[im3]); // this only makes sense if the pressures are positive bool ptest = p0_hse < 0.0 || pp1_hse < 0.0 || pp2_hse < 0.0 || + pp3_hse < 0.0 || pm1_hse < 0.0 || - pm2_hse < 0.0; + pm2_hse < 0.0 || + pm3_hse < 0.0; if (!ptest && rho[i0] >= pslope_cutoff_density) { @@ -277,19 +281,23 @@ ppm_reconstruct_pslope(const Real* rho, const Real* p, const Real* src, tp[ip1] = p[ip1] - pp1_hse; tp[ip2] = p[ip2] - pp2_hse; + tp[ip3] = p[ip3] - pp3_hse; tp[im1] = p[im1] - pm1_hse; tp[im2] = p[im2] - pm2_hse; + tp[im3] = p[im3] - pm3_hse; } else { // don't subtract off HSE + tp[im3] = p[im3]; tp[im2] = p[im2]; tp[im1] = p[im1]; tp[i0] = p[i0]; tp[ip1]= p[ip1]; tp[ip2] = p[ip2]; + tp[ip3] = p[ip3]; } ppm_reconstruct(tp, i, j, k, idir, lo_bc_test, hi_bc_test, domlo, domhi, flatn, sm, sp); From 0a4044b467c2e140e0251d9ad607d46edc8ad27f Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Mon, 25 Jul 2022 16:03:01 -0400 Subject: [PATCH 14/24] change centering of source terms --- Source/hydro/ppm.H | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/Source/hydro/ppm.H b/Source/hydro/ppm.H index 575f2cc0f7..3523111c68 100644 --- a/Source/hydro/ppm.H +++ b/Source/hydro/ppm.H @@ -255,13 +255,13 @@ ppm_reconstruct_pslope(const Real* rho, const Real* p, const Real* src, 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]); - Real pp3_hse = pp2_hse + 0.25_rt*dx * (rho[ip2] + rho[ip3]) * (src[ip2] + src[ip3]); + Real pp1_hse = p0_hse + 0.5_rt*dx * (rho[i0] * src[i0] + rho[ip1] * src[ip1]); + Real pp2_hse = pp1_hse + 0.5_rt*dx * (rho[ip1] * src[ip1] + rho[ip2] * src[ip2]); + Real pp3_hse = pp2_hse + 0.5_rt*dx * (rho[ip2] * src[ip2] + rho[ip3] * src[ip3]); - 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]); - Real pm3_hse = pm2_hse - 0.25_rt*dx * (rho[im2] + rho[im3]) * (src[im2] + src[im3]); + Real pm1_hse = p0_hse - 0.5_rt*dx * (rho[i0] * src[i0] + rho[im1] * src[im1]); + Real pm2_hse = pm1_hse - 0.5_rt*dx * (rho[im1] * src[im1] + rho[im2] * src[im2]); + Real pm3_hse = pm2_hse - 0.5_rt*dx * (rho[im2] * src[im2] + rho[im3] * src[im3]); // this only makes sense if the pressures are positive bool ptest = p0_hse < 0.0 || From 31d18554e11a4bcc9cf7ddbe0897d86e6b0de72f Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Tue, 26 Jul 2022 12:13:25 -0400 Subject: [PATCH 15/24] fix array sizes --- Source/hydro/trace_ppm.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Source/hydro/trace_ppm.cpp b/Source/hydro/trace_ppm.cpp index aef88ea3bc..bea52d9f74 100644 --- a/Source/hydro/trace_ppm.cpp +++ b/Source/hydro/trace_ppm.cpp @@ -224,8 +224,8 @@ Castro::trace_ppm(const Box& bx, load_stencil(q_arr, idir, i, j, k, QPRES, s); if (use_pslope) { - Real trho[5]; - Real src[5]; + Real trho[7]; + Real src[7]; load_stencil(q_arr, idir, i, j, k, QRHO, trho); load_stencil(srcQ, idir, i, j, k, QUN, src); From 4ad48f240507c9199548d560bc9f97723d5c050b Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Sun, 16 Apr 2023 09:59:14 -0400 Subject: [PATCH 16/24] sync --- 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 a157710cac..445cb9bbcf 100644 --- a/Source/hydro/Castro_ctu.cpp +++ b/Source/hydro/Castro_ctu.cpp @@ -116,7 +116,7 @@ Castro::ctu_ppm_states(const Box& bx, const Box& vbx, for (int idir = 0; idir < AMREX_SPACEDIM; idir++) { if (idir == 0) { - trace_ppm(bx, + trace_ppm(bx, idir, U_arr, rho_inv_arr, q_arr, qaux_arr, srcQ, qxm, qxp, From ff8fcf70ea53a827d242ce54c791c189d439e6ac Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Sun, 16 Apr 2023 11:15:02 -0400 Subject: [PATCH 17/24] fix --- Source/hydro/Castro_mol.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Source/hydro/Castro_mol.cpp b/Source/hydro/Castro_mol.cpp index 7401d8fc94..87750e8fb7 100644 --- a/Source/hydro/Castro_mol.cpp +++ b/Source/hydro/Castro_mol.cpp @@ -49,7 +49,7 @@ Castro::mol_plm_reconstruct(const Box& bx, (idir == 1 && j == domhi[1]) || (idir == 2 && k == domhi[2])); - Real s[7]; + Real s[nslp]; Real flat = flatn_arr(i,j,k); load_stencil(q_arr, idir, i, j, k, n, s); From 4e87524b09c089c1cbe876e4face20343df1e538 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Sun, 16 Apr 2023 14:12:11 -0400 Subject: [PATCH 18/24] implement the axisymmetric versions of the one-sided diffs --- Source/hydro/Castro_mol.cpp | 4 ++- Source/hydro/ppm.H | 58 +++++++++++++++++++++++-------------- Source/hydro/trace_ppm.cpp | 32 ++++++++++---------- 3 files changed, 57 insertions(+), 37 deletions(-) diff --git a/Source/hydro/Castro_mol.cpp b/Source/hydro/Castro_mol.cpp index 87750e8fb7..7f3747eece 100644 --- a/Source/hydro/Castro_mol.cpp +++ b/Source/hydro/Castro_mol.cpp @@ -162,6 +162,8 @@ Castro::mol_ppm_reconstruct(const Box& bx, bool lo_bc_test = lo_bc[idir] == Symmetry; bool hi_bc_test = hi_bc[idir] == Symmetry; + bool is_axisymmetric = geom.Coord() == 1 +; amrex::ParallelFor(bx, NQ, [=] AMREX_GPU_HOST_DEVICE (int i, int j, int k, int n) { @@ -173,7 +175,7 @@ Castro::mol_ppm_reconstruct(const Box& bx, load_stencil(q_arr, idir, i, j, k, n, s); ppm_reconstruct(s, i, j, k, idir, - lo_bc_test, hi_bc_test, domlo, domhi, + lo_bc_test, hi_bc_test, is_axisymmetric, domlo, domhi, flat, sm, sp); if (idir == 0) { diff --git a/Source/hydro/ppm.H b/Source/hydro/ppm.H index 99986da5f2..9a190a52a8 100644 --- a/Source/hydro/ppm.H +++ b/Source/hydro/ppm.H @@ -55,7 +55,7 @@ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void ppm_reconstruct(const Real* s, const int i, const int j, const int k, const int idir, - const bool lo_bc_test, const bool hi_bc_test, + const bool lo_bc_test, const bool hi_bc_test, const bool is_axisymmetric, const GpuArray& domlo, const GpuArray& domhi, const Real flatn, Real& sm, Real& sp) { @@ -65,37 +65,48 @@ ppm_reconstruct(const Real* s, (idir == 1 && j == domlo[1]) || (idir == 2 && k == domlo[2]))) { - // use a stencil for when the current zone is on the left physical - // boundary. Then the left interface is on the physical boundary, - // MC Eq. 21 - sm = (1.0_rt/12.0_rt)*(25.0_rt*s[i0] - 23.0_rt*s[ip1] + 13.0_rt*s[ip2] - 3.0_rt*s[ip3]); + // use a stencil for when the current zone is on the left physical + // boundary. Then the left interface is on the physical boundary, + // MC Eq. 21 + + if (idir == 0 && is_axisymmetric) { + sm = (1.0_rt/144.0_rt)*(415.0_rt*s[i0] - 483.0_rt*s[ip1] + 275.0_rt*s[ip2] - 63.0_rt*s[ip3]); + } else { + // Cartesian + sm = (1.0_rt/12.0_rt)*(25.0_rt*s[i0] - 23.0_rt*s[ip1] + 13.0_rt*s[ip2] - 3.0_rt*s[ip3]); + } } else if (lo_bc_test && ((idir == 0 && i == domlo[0]+1) || (idir == 1 && j == domlo[1]+1) || (idir == 2 && k == domlo[2]+1))) { - // use a stencil for when the current zone is one zone away from - // the left physical boundary, and then the left interface is one - // zone away from the boundary, MC Eq. 22 - sm = (1.0_rt/12.0_rt)*(3.0_rt*s[im1] + 13.0_rt*s[i0] - 5.0_rt*s[ip1] + s[ip2]); + // use a stencil for when the current zone is one zone away from + // the left physical boundary, and then the left interface is one + // zone away from the boundary, MC Eq. 22 - // Make sure sedge lies in between adjacent cell-centered values + if (idir == 0 && is_axisymmetric) { + sm = (1.0_rt/96.0_rt)*(37.0_rt*s[im1] + 87.0_rt*s[i0] - 35.0_rt*s[ip1] + 7.0_rt*s[ip2]); + } else { + sm = (1.0_rt/12.0_rt)*(3.0_rt*s[im1] + 13.0_rt*s[i0] - 5.0_rt*s[ip1] + s[ip2]); + } - sm = amrex::max(sm, amrex::min(s[i0], s[im1])); - sm = amrex::min(sm, amrex::max(s[i0], s[im1])); + // Make sure sedge lies in between adjacent cell-centered values + + sm = amrex::max(sm, amrex::min(s[i0], s[im1])); + sm = amrex::min(sm, amrex::max(s[i0], s[im1])); } else if (hi_bc_test && ((idir == 0 && i == domhi[0]) || (idir == 1 && j == domhi[1]) || (idir == 2 && k == domhi[2]))) { - // use a stencil for when the current zone is on the right physical boundary - // then the left interface is one zone away from the physical boundary - sm = (1.0_rt/12.0_rt)*(3.0_rt*s[i0] + 13.0_rt*s[im1] - 5.0_rt*s[im2] + s[im3]); + // use a stencil for when the current zone is on the right physical boundary + // then the left interface is one zone away from the physical boundary + sm = (1.0_rt/12.0_rt)*(3.0_rt*s[i0] + 13.0_rt*s[im1] - 5.0_rt*s[im2] + s[im3]); - // Make sure sedge lies in between adjacent cell-centered values + // Make sure sedge lies in between adjacent cell-centered values - sm = amrex::max(sm, amrex::min(s[i0], s[im1])); - sm = amrex::min(sm, amrex::max(s[i0], s[im1])); + sm = amrex::max(sm, amrex::min(s[i0], s[im1])); + sm = amrex::min(sm, amrex::max(s[i0], s[im1])); } else { @@ -139,7 +150,12 @@ ppm_reconstruct(const Real* s, // use a stencil for when the current zone is on the left physical // boundary. Then the right interface is one zone away from the // physical boundary, - sp = (1.0_rt/12.0_rt)*(3.0_rt*s[i0] + 13.0_rt*s[ip1] - 5.0_rt*s[ip2] + s[ip3]); + + if (idir == 0 && is_axisymmetric) { + sp = (1.0_rt/96.0_rt)*(37.0_rt*s[i0] + 87.0_rt*s[ip1] - 35.0_rt*s[ip2] + 7.0_rt*s[ip3]); + } else { + sp = (1.0_rt/12.0_rt)*(3.0_rt*s[i0] + 13.0_rt*s[ip1] - 5.0_rt*s[ip2] + s[ip3]); + } // Make sure sedge lies in between adjacent cell-centered values @@ -244,7 +260,7 @@ void ppm_reconstruct_pslope(const Real* rho, const Real* p, const Real* src, const Real dx, const int i, const int j, const int k, const int idir, - bool lo_bc_test, bool hi_bc_test, + bool lo_bc_test, bool hi_bc_test, const bool is_axisymmetric, const GpuArray& domlo, const GpuArray& domhi, const Real flatn, Real& sm, Real& sp) { @@ -300,7 +316,7 @@ ppm_reconstruct_pslope(const Real* rho, const Real* p, const Real* src, tp[ip3] = p[ip3]; } - ppm_reconstruct(tp, i, j, k, idir, lo_bc_test, hi_bc_test, domlo, domhi, flatn, sm, sp); + ppm_reconstruct(tp, i, j, k, idir, lo_bc_test, hi_bc_test, is_axisymmetric, domlo, domhi, flatn, sm, sp); // now correct sm and sp to be back to the full pressure by diff --git a/Source/hydro/trace_ppm.cpp b/Source/hydro/trace_ppm.cpp index 35ec1344d2..4fc7290120 100644 --- a/Source/hydro/trace_ppm.cpp +++ b/Source/hydro/trace_ppm.cpp @@ -57,6 +57,8 @@ Castro::trace_ppm(const Box& bx, const auto dx = geom.CellSizeArray(); + bool is_axisymmetric = geom.Coord() == 1; + Real hdt = 0.5_rt * dt; Real dtdx = dt / dx[idir]; @@ -204,7 +206,7 @@ Castro::trace_ppm(const Box& bx, Real Im_rho[3]; load_stencil(q_arr, idir, i, j, k, QRHO, s); - ppm_reconstruct(s, i, j, k, idir, lo_symm, hi_symm, domlo, domhi, flat, sm, sp); + ppm_reconstruct(s, i, j, k, idir, lo_symm, hi_symm, is_axisymmetric, domlo, domhi, flat, sm, sp); sm = amrex::max(lsmall_dens, sm); sp = amrex::max(lsmall_dens, sp); ppm_int_profile(sm, sp, s[i0], un, cc, dtdx, Ip_rho, Im_rho); @@ -217,7 +219,7 @@ Castro::trace_ppm(const Box& bx, Real Im_un_2; load_stencil(q_arr, idir, i, j, k, QUN, s); - ppm_reconstruct(s, i, j, k, idir, lo_symm, hi_symm, domlo, domhi, flat, sm, sp); + ppm_reconstruct(s, i, j, k, idir, lo_symm, hi_symm, is_axisymmetric, domlo, domhi, flat, sm, sp); ppm_int_profile_single(sm, sp, s[i0], un-cc, dtdx, Ip_un_0, Im_un_0); ppm_int_profile_single(sm, sp, s[i0], un+cc, dtdx, Ip_un_2, Im_un_2); @@ -235,10 +237,10 @@ Castro::trace_ppm(const Box& bx, 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, dx[idir], i, j, k, idir, lo_symm, hi_symm, domlo, domhi, flat, sm, sp); + ppm_reconstruct_pslope(trho, s, src, dx[idir], i, j, k, idir, lo_symm, hi_symm, is_axisymmetric, domlo, domhi, flat, sm, sp); } else { - ppm_reconstruct(s, i, j, k, idir, lo_symm, hi_symm, domlo, domhi, flat, sm, sp); + ppm_reconstruct(s, i, j, k, idir, lo_symm, hi_symm, is_axisymmetric, domlo, domhi, flat, sm, sp); } ppm_int_profile(sm, sp, s[i0], un, cc, dtdx, Ip_p, Im_p); @@ -248,7 +250,7 @@ Castro::trace_ppm(const Box& bx, Real Im_rhoe[3]; load_stencil(q_arr, idir, i, j, k, QREINT, s); - ppm_reconstruct(s, i, j, k, idir, lo_symm, hi_symm, domlo, domhi, flat, sm, sp); + ppm_reconstruct(s, i, j, k, idir, lo_symm, hi_symm, is_axisymmetric, domlo, domhi, flat, sm, sp); sm = amrex::max(lsmall_pres, sm); sp = amrex::max(lsmall_pres, sp); ppm_int_profile(sm, sp, s[i0], un, cc, dtdx, Ip_rhoe, Im_rhoe); @@ -261,11 +263,11 @@ Castro::trace_ppm(const Box& bx, Real Im_utt_1; load_stencil(q_arr, idir, i, j, k, QUT, s); - ppm_reconstruct(s, i, j, k, idir, lo_symm, hi_symm, domlo, domhi, flat, sm, sp); + ppm_reconstruct(s, i, j, k, idir, lo_symm, hi_symm, is_axisymmetric, domlo, domhi, flat, sm, sp); ppm_int_profile_single(sm, sp, s[i0], un, dtdx, Ip_ut_1, Im_ut_1); load_stencil(q_arr, idir, i, j, k, QUTT, s); - ppm_reconstruct(s, i, j, k, idir, lo_symm, hi_symm, domlo, domhi, flat, sm, sp); + ppm_reconstruct(s, i, j, k, idir, lo_symm, hi_symm, is_axisymmetric, domlo, domhi, flat, sm, sp); ppm_int_profile_single(sm, sp, s[i0], un, dtdx, Ip_utt_1, Im_utt_1); // gamma_c @@ -276,7 +278,7 @@ Castro::trace_ppm(const Box& bx, Real Im_gc_2; load_stencil(qaux_arr, idir, i, j, k, QGAMC, s); - ppm_reconstruct(s, i, j, k, idir, lo_symm, hi_symm, domlo, domhi, flat, sm, sp); + ppm_reconstruct(s, i, j, k, idir, lo_symm, hi_symm, is_axisymmetric, domlo, domhi, flat, sm, sp); ppm_int_profile_single(sm, sp, s[i0], un-cc, dtdx, Ip_gc_0, Im_gc_0); ppm_int_profile_single(sm, sp, s[i0], un+cc, dtdx, Ip_gc_2, Im_gc_2); @@ -297,7 +299,7 @@ Castro::trace_ppm(const Box& bx, if (do_trace) { load_stencil(srcQ, idir, i, j, k, QRHO, s); - ppm_reconstruct(s, i, j, k, idir, lo_symm, hi_symm, domlo, domhi, flat, sm, sp); + ppm_reconstruct(s, i, j, k, idir, lo_symm, hi_symm, is_axisymmetric, domlo, domhi, flat, sm, sp); ppm_int_profile(sm, sp, s[i0], un, cc, dtdx, Ip_src_rho, Im_src_rho); } @@ -316,7 +318,7 @@ Castro::trace_ppm(const Box& bx, if (do_trace) { load_stencil(srcQ, idir, i, j, k, QUN, s); - ppm_reconstruct(s, i, j, k, idir, lo_symm, hi_symm, domlo, domhi, flat, sm, sp); + ppm_reconstruct(s, i, j, k, idir, lo_symm, hi_symm, is_axisymmetric, domlo, domhi, flat, sm, sp); ppm_int_profile_single(sm, sp, s[i0], un-cc, dtdx, Ip_src_un_0, Im_src_un_0); ppm_int_profile_single(sm, sp, s[i0], un+cc, dtdx, Ip_src_un_2, Im_src_un_2); } @@ -334,7 +336,7 @@ Castro::trace_ppm(const Box& bx, if (do_trace) { load_stencil(srcQ, idir, i, j, k, QPRES, s); - ppm_reconstruct(s, i, j, k, idir, lo_symm, hi_symm, domlo, domhi, flat, sm, sp); + ppm_reconstruct(s, i, j, k, idir, lo_symm, hi_symm, is_axisymmetric, domlo, domhi, flat, sm, sp); ppm_int_profile(sm, sp, s[i0], un, cc, dtdx, Ip_src_p, Im_src_p); } @@ -351,7 +353,7 @@ Castro::trace_ppm(const Box& bx, if (do_trace) { load_stencil(srcQ, idir, i, j, k, QREINT, s); - ppm_reconstruct(s, i, j, k, idir, lo_symm, hi_symm, domlo, domhi, flat, sm, sp); + ppm_reconstruct(s, i, j, k, idir, lo_symm, hi_symm, is_axisymmetric, domlo, domhi, flat, sm, sp); ppm_int_profile(sm, sp, s[i0], un, cc, dtdx, Ip_src_rhoe, Im_src_rhoe); } @@ -368,7 +370,7 @@ Castro::trace_ppm(const Box& bx, if (do_trace) { load_stencil(srcQ, idir, i, j, k, QUT, s); - ppm_reconstruct(s, i, j, k, idir, lo_symm, hi_symm, domlo, domhi, flat, sm, sp); + ppm_reconstruct(s, i, j, k, idir, lo_symm, hi_symm, is_axisymmetric, domlo, domhi, flat, sm, sp); ppm_int_profile_single(sm, sp, s[i0], un, dtdx, Ip_src_ut_1, Im_src_ut_1); } @@ -383,7 +385,7 @@ Castro::trace_ppm(const Box& bx, if (do_trace) { load_stencil(srcQ, idir, i, j, k, QUTT, s); - ppm_reconstruct(s, i, j, k, idir, lo_symm, hi_symm, domlo, domhi, flat, sm, sp); + ppm_reconstruct(s, i, j, k, idir, lo_symm, hi_symm, is_axisymmetric, domlo, domhi, flat, sm, sp); ppm_int_profile_single(sm, sp, s[i0], un, dtdx, Ip_src_utt_1, Im_src_utt_1); } @@ -401,7 +403,7 @@ Castro::trace_ppm(const Box& bx, int n = qpassmap(ipassive); load_passive_stencil(U_arr, rho_inv_arr, idir, i, j, k, nc, s); - ppm_reconstruct(s, i, j, k, idir, lo_symm, hi_symm, domlo, domhi, flat, sm, sp); + ppm_reconstruct(s, i, j, k, idir, lo_symm, hi_symm, is_axisymmetric, domlo, domhi, flat, sm, sp); ppm_int_profile_single(sm, sp, s[i0], un, dtdx, Ip_passive, Im_passive); // Plus state on face i From 0a7c575575bef455be0c20b556108f176b4ff387 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Thu, 25 May 2023 20:22:09 -0400 Subject: [PATCH 19/24] more work on the one-sided stuff --- Source/hydro/Castro_mol.cpp | 4 +- Source/hydro/ppm.H | 128 +++++++++++++++++++++++------------- Source/hydro/trace_ppm.cpp | 30 ++++----- 3 files changed, 101 insertions(+), 61 deletions(-) diff --git a/Source/hydro/Castro_mol.cpp b/Source/hydro/Castro_mol.cpp index 7f3747eece..941def9a6e 100644 --- a/Source/hydro/Castro_mol.cpp +++ b/Source/hydro/Castro_mol.cpp @@ -159,6 +159,8 @@ Castro::mol_ppm_reconstruct(const Box& bx, const auto domlo = geom.Domain().loVect3d(); const auto domhi = geom.Domain().hiVect3d(); + const auto dx = geom.CellSizeArray(); + bool lo_bc_test = lo_bc[idir] == Symmetry; bool hi_bc_test = hi_bc[idir] == Symmetry; @@ -174,7 +176,7 @@ Castro::mol_ppm_reconstruct(const Box& bx, Real sp; load_stencil(q_arr, idir, i, j, k, n, s); - ppm_reconstruct(s, i, j, k, idir, + ppm_reconstruct(s, i, j, k, idir, dx, lo_bc_test, hi_bc_test, is_axisymmetric, domlo, domhi, flat, sm, sp); diff --git a/Source/hydro/ppm.H b/Source/hydro/ppm.H index 9a190a52a8..fea7ec624c 100644 --- a/Source/hydro/ppm.H +++ b/Source/hydro/ppm.H @@ -55,6 +55,7 @@ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void ppm_reconstruct(const Real* s, const int i, const int j, const int k, const int idir, + const GpuArray& dx, const bool lo_bc_test, const bool hi_bc_test, const bool is_axisymmetric, const GpuArray& domlo, const GpuArray& domhi, const Real flatn, Real& sm, Real& sp) { @@ -110,35 +111,53 @@ ppm_reconstruct(const Real* s, } else { - // Compute van Leer slopes + if (idir == 0 && is_axisymmetric) { - Real dsl = 2.0_rt * (s[im1] - s[im2]); - Real dsr = 2.0_rt * (s[i0] - s[im1]); + const Real dr = dx[0]; + const Real r_m12 = (static_cast(i) - 0.5_rt) * dr; - Real dsvl_l = 0.0_rt; - if (dsl*dsr > 0.0_rt) { - Real dsc = 0.5_rt * (s[i0] - s[im2]); - dsvl_l = std::copysign(1.0_rt, dsc) * amrex::min(std::abs(dsc),amrex::min(std::abs(dsl),std::abs(dsr))); - } + Real dr4_term = std::pow(dr, 4) * (-12.0_rt * s[im2] + 60.0_rt * s[im1] + 60.0_rt * s[i0] - 12.0_rt * s[ip1]); + Real dr3_term = r_m12 * std::pow(dr, 3) * (-s[im2] - 27.0_rt * s[im1] + 27.0_rt * s[i0] + s[ip1]); + Real dr2_term = std::pow(r_m12, 2) * std::pow(dr, 2) * (30.0_rt * s[im2] - 210.0_rt * s[im1] - 210.0_rt * s[i0] + 30.0_rt * s[ip1]); + Real dr1_term = std::pow(r_m12, 3) * dr * (-s[im2] - 13.0_rt * s[im1] - 13.0_rt * s[i0] + s[ip1]); + Real dr0_term = std::pow(r_m12, 4) * (-10.0_rt * s[im2] + 70.0_rt * s[im1] + 70.0_rt * s[i0] - 10.0_rt * s[ip1]); - dsl = 2.0_rt * (s[i0] - s[im1]); - dsr = 2.0_rt * (s[ip1] - s[i0]); + Real denom = 24.0_rt * (4.0_rt * std::pow(dr, 4) - 15.0_rt * std::pow(dr, 2) * std::pow(r_m12, 2) + 5.0_rt * std::pow(r_m12, 4)); - Real dsvl_r = 0.0_rt; - if (dsl*dsr > 0.0_rt) { - Real dsc = 0.5_rt * (s[ip1] - s[im1]); - dsvl_r = std::copysign(1.0_rt, dsc) * amrex::min(std::abs(dsc),amrex::min(std::abs(dsl),std::abs(dsr))); - } + sm = (dr4_term + dr3_term + dr2_term + dr1_term + dr0_term) / denom; + + } else { + + // Compute van Leer slopes + + Real dsl = 2.0_rt * (s[im1] - s[im2]); + Real dsr = 2.0_rt * (s[i0] - s[im1]); - // Interpolate s to edges + Real dsvl_l = 0.0_rt; + if (dsl*dsr > 0.0_rt) { + Real dsc = 0.5_rt * (s[i0] - s[im2]); + dsvl_l = std::copysign(1.0_rt, dsc) * amrex::min(std::abs(dsc),amrex::min(std::abs(dsl),std::abs(dsr))); + } - sm = 0.5_rt * (s[i0] + s[im1]) - (1.0_rt/6.0_rt) * (dsvl_r - dsvl_l); + dsl = 2.0_rt * (s[i0] - s[im1]); + dsr = 2.0_rt * (s[ip1] - s[i0]); + + Real dsvl_r = 0.0_rt; + if (dsl*dsr > 0.0_rt) { + Real dsc = 0.5_rt * (s[ip1] - s[im1]); + dsvl_r = std::copysign(1.0_rt, dsc) * amrex::min(std::abs(dsc),amrex::min(std::abs(dsl),std::abs(dsr))); + } + + // Interpolate s to edges + + sm = 0.5_rt * (s[i0] + s[im1]) - (1.0_rt/6.0_rt) * (dsvl_r - dsvl_l); + + } // Make sure sedge lies in between adjacent cell-centered values sm = amrex::max(sm, amrex::min(s[i0], s[im1])); sm = amrex::min(sm, amrex::max(s[i0], s[im1])); - } // now we compute s_{i+1/2} -- the right interface value for zone i @@ -186,29 +205,48 @@ ppm_reconstruct(const Real* s, } else { - // Compute van Leer slopes + if (idir == 0 && is_axisymmetric) { - Real dsl = 2.0_rt * (s[i0] - s[im1]); - Real dsr = 2.0_rt * (s[ip1] - s[i0]); + const Real dr = dx[0]; + const Real r_p12 = (static_cast(i) + 0.5_rt) * dr; - Real 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), amrex::min(std::abs(dsl),std::abs(dsr))); - } + Real dr4_term = std::pow(dr, 4) * (-12.0_rt * s[im1] + 60.0_rt * s[i0] + 60.0_rt * s[ip1] - 12.0_rt * s[ip2]); + Real dr3_term = r_p12 * std::pow(dr, 3) * (-s[im1] - 27.0_rt * s[i0] + 27.0_rt * s[ip1] + s[ip2]); + Real dr2_term = std::pow(r_p12, 2) * std::pow(dr, 2) * (30.0_rt * s[im1] - 210.0_rt * s[i0] - 210.0_rt * s[ip1] + 30.0_rt * s[ip2]); + Real dr1_term = std::pow(r_p12, 3) * dr * (-s[im1] - 13.0_rt * s[i0] - 13.0_rt * s[ip1] + s[ip2]); + Real dr0_term = std::pow(r_p12, 4) * (-10.0_rt * s[im1] + 70.0_rt * s[i0] + 70.0_rt * s[ip1] - 10.0_rt * s[ip2]); - dsl = 2.0_rt * (s[ip1] - s[i0]); - dsr = 2.0_rt * (s[ip2] - s[ip1]); + Real denom = 24.0_rt * (4.0_rt * std::pow(dr, 4) - 15.0_rt * std::pow(dr, 2) * std::pow(r_p12, 2) + 5.0_rt * std::pow(r_p12, 4)); - Real 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), amrex::min(std::abs(dsl),std::abs(dsr))); - } + sp = (dr4_term + dr3_term + dr2_term + dr1_term + dr0_term) / denom; + + } else { + + // Compute van Leer slopes - // Interpolate s to edges + Real dsl = 2.0_rt * (s[i0] - s[im1]); + Real dsr = 2.0_rt * (s[ip1] - s[i0]); - sp = 0.5_rt * (s[ip1] + s[i0]) - (1.0_rt/6.0_rt) * (dsvl_r - dsvl_l); + Real 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), amrex::min(std::abs(dsl),std::abs(dsr))); + } + + dsl = 2.0_rt * (s[ip1] - s[i0]); + dsr = 2.0_rt * (s[ip2] - s[ip1]); + + Real dsvl_r = 0.0_rt; + 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), amrex::min(std::abs(dsl),std::abs(dsr))); + } + + // Interpolate s to edges + + sp = 0.5_rt * (s[ip1] + s[i0]) - (1.0_rt/6.0_rt) * (dsvl_r - dsvl_l); + + } // Make sure sedge lies in between adjacent cell-centered values @@ -258,7 +296,7 @@ ppm_reconstruct(const Real* s, AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void ppm_reconstruct_pslope(const Real* rho, const Real* p, const Real* src, - const Real dx, + const GpuArray& dx, const int i, const int j, const int k, const int idir, bool lo_bc_test, bool hi_bc_test, const bool is_axisymmetric, const GpuArray& domlo, const GpuArray& domhi, @@ -271,13 +309,13 @@ ppm_reconstruct_pslope(const Real* rho, const Real* p, const Real* src, Real p0_hse = p[i0]; - Real pp1_hse = p0_hse + 0.5_rt*dx * (rho[i0] * src[i0] + rho[ip1] * src[ip1]); - Real pp2_hse = pp1_hse + 0.5_rt*dx * (rho[ip1] * src[ip1] + rho[ip2] * src[ip2]); - Real pp3_hse = pp2_hse + 0.5_rt*dx * (rho[ip2] * src[ip2] + rho[ip3] * src[ip3]); + Real pp1_hse = p0_hse + 0.5_rt * dx[idir] * (rho[i0] * src[i0] + rho[ip1] * src[ip1]); + Real pp2_hse = pp1_hse + 0.5_rt * dx[idir] * (rho[ip1] * src[ip1] + rho[ip2] * src[ip2]); + Real pp3_hse = pp2_hse + 0.5_rt * dx[idir] * (rho[ip2] * src[ip2] + rho[ip3] * src[ip3]); - Real pm1_hse = p0_hse - 0.5_rt*dx * (rho[i0] * src[i0] + rho[im1] * src[im1]); - Real pm2_hse = pm1_hse - 0.5_rt*dx * (rho[im1] * src[im1] + rho[im2] * src[im2]); - Real pm3_hse = pm2_hse - 0.5_rt*dx * (rho[im2] * src[im2] + rho[im3] * src[im3]); + Real pm1_hse = p0_hse - 0.5_rt * dx[idir] * (rho[i0] * src[i0] + rho[im1] * src[im1]); + Real pm2_hse = pm1_hse - 0.5_rt * dx[idir] * (rho[im1] * src[im1] + rho[im2] * src[im2]); + Real pm3_hse = pm2_hse - 0.5_rt * dx[idir] * (rho[im2] * src[im2] + rho[im3] * src[im3]); // this only makes sense if the pressures are positive bool ptest = p0_hse < 0.0 || @@ -316,15 +354,15 @@ ppm_reconstruct_pslope(const Real* rho, const Real* p, const Real* src, tp[ip3] = p[ip3]; } - ppm_reconstruct(tp, i, j, k, idir, lo_bc_test, hi_bc_test, is_axisymmetric, domlo, domhi, flatn, sm, sp); + ppm_reconstruct(tp, i, j, k, idir, dx, lo_bc_test, hi_bc_test, is_axisymmetric, domlo, domhi, flatn, sm, sp); // now correct sm and sp to be back to the full pressure by // adding in the hydrostatic pressure at the interface 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]; + sp += p[i0] + 0.5_rt * dx[idir] * rho[i0] * src[i0]; + sm += p[i0] - 0.5_rt * dx[idir] * rho[i0] * src[i0]; } } diff --git a/Source/hydro/trace_ppm.cpp b/Source/hydro/trace_ppm.cpp index 4623a0a529..e4845c0538 100644 --- a/Source/hydro/trace_ppm.cpp +++ b/Source/hydro/trace_ppm.cpp @@ -206,7 +206,7 @@ Castro::trace_ppm(const Box& bx, Real Im_rho[3]; load_stencil(q_arr, idir, i, j, k, QRHO, s); - ppm_reconstruct(s, i, j, k, idir, lo_symm, hi_symm, is_axisymmetric, domlo, domhi, flat, sm, sp); + ppm_reconstruct(s, i, j, k, idir, dx, lo_symm, hi_symm, is_axisymmetric, domlo, domhi, flat, sm, sp); sm = amrex::max(lsmall_dens, sm); sp = amrex::max(lsmall_dens, sp); ppm_int_profile(sm, sp, s[i0], un, cc, dtdx, Ip_rho, Im_rho); @@ -219,7 +219,7 @@ Castro::trace_ppm(const Box& bx, Real Im_un_2; load_stencil(q_arr, idir, i, j, k, QUN, s); - ppm_reconstruct(s, i, j, k, idir, lo_symm, hi_symm, is_axisymmetric, domlo, domhi, flat, sm, sp); + ppm_reconstruct(s, i, j, k, idir, dx, lo_symm, hi_symm, is_axisymmetric, domlo, domhi, flat, sm, sp); ppm_int_profile_single(sm, sp, s[i0], un-cc, dtdx, Ip_un_0, Im_un_0); ppm_int_profile_single(sm, sp, s[i0], un+cc, dtdx, Ip_un_2, Im_un_2); @@ -237,10 +237,10 @@ Castro::trace_ppm(const Box& bx, 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, dx[idir], i, j, k, idir, lo_symm, hi_symm, is_axisymmetric, domlo, domhi, flat, sm, sp); + ppm_reconstruct_pslope(trho, s, src, dx, i, j, k, idir, lo_symm, hi_symm, is_axisymmetric, domlo, domhi, flat, sm, sp); } else { - ppm_reconstruct(s, i, j, k, idir, lo_symm, hi_symm, is_axisymmetric, domlo, domhi, flat, sm, sp); + ppm_reconstruct(s, i, j, k, idir, dx, lo_symm, hi_symm, is_axisymmetric, domlo, domhi, flat, sm, sp); } ppm_int_profile(sm, sp, s[i0], un, cc, dtdx, Ip_p, Im_p); @@ -250,7 +250,7 @@ Castro::trace_ppm(const Box& bx, Real Im_rhoe[3]; load_stencil(q_arr, idir, i, j, k, QREINT, s); - ppm_reconstruct(s, i, j, k, idir, lo_symm, hi_symm, is_axisymmetric, domlo, domhi, flat, sm, sp); + ppm_reconstruct(s, i, j, k, idir, dx, lo_symm, hi_symm, is_axisymmetric, domlo, domhi, flat, sm, sp); sm = amrex::max(lsmall_pres, sm); sp = amrex::max(lsmall_pres, sp); ppm_int_profile(sm, sp, s[i0], un, cc, dtdx, Ip_rhoe, Im_rhoe); @@ -263,11 +263,11 @@ Castro::trace_ppm(const Box& bx, Real Im_utt_1; load_stencil(q_arr, idir, i, j, k, QUT, s); - ppm_reconstruct(s, i, j, k, idir, lo_symm, hi_symm, is_axisymmetric, domlo, domhi, flat, sm, sp); + ppm_reconstruct(s, i, j, k, idir, dx, lo_symm, hi_symm, is_axisymmetric, domlo, domhi, flat, sm, sp); ppm_int_profile_single(sm, sp, s[i0], un, dtdx, Ip_ut_1, Im_ut_1); load_stencil(q_arr, idir, i, j, k, QUTT, s); - ppm_reconstruct(s, i, j, k, idir, lo_symm, hi_symm, is_axisymmetric, domlo, domhi, flat, sm, sp); + ppm_reconstruct(s, i, j, k, idir, dx, lo_symm, hi_symm, is_axisymmetric, domlo, domhi, flat, sm, sp); ppm_int_profile_single(sm, sp, s[i0], un, dtdx, Ip_utt_1, Im_utt_1); // gamma_c @@ -278,7 +278,7 @@ Castro::trace_ppm(const Box& bx, Real Im_gc_2; load_stencil(qaux_arr, idir, i, j, k, QGAMC, s); - ppm_reconstruct(s, i, j, k, idir, lo_symm, hi_symm, is_axisymmetric, domlo, domhi, flat, sm, sp); + ppm_reconstruct(s, i, j, k, idir, dx, lo_symm, hi_symm, is_axisymmetric, domlo, domhi, flat, sm, sp); ppm_int_profile_single(sm, sp, s[i0], un-cc, dtdx, Ip_gc_0, Im_gc_0); ppm_int_profile_single(sm, sp, s[i0], un+cc, dtdx, Ip_gc_2, Im_gc_2); @@ -299,7 +299,7 @@ Castro::trace_ppm(const Box& bx, if (do_trace) { load_stencil(srcQ, idir, i, j, k, QRHO, s); - ppm_reconstruct(s, i, j, k, idir, lo_symm, hi_symm, is_axisymmetric, domlo, domhi, flat, sm, sp); + ppm_reconstruct(s, i, j, k, idir, dx, lo_symm, hi_symm, is_axisymmetric, domlo, domhi, flat, sm, sp); ppm_int_profile(sm, sp, s[i0], un, cc, dtdx, Ip_src_rho, Im_src_rho); } @@ -318,7 +318,7 @@ Castro::trace_ppm(const Box& bx, if (do_trace) { load_stencil(srcQ, idir, i, j, k, QUN, s); - ppm_reconstruct(s, i, j, k, idir, lo_symm, hi_symm, is_axisymmetric, domlo, domhi, flat, sm, sp); + ppm_reconstruct(s, i, j, k, idir, dx, lo_symm, hi_symm, is_axisymmetric, domlo, domhi, flat, sm, sp); ppm_int_profile_single(sm, sp, s[i0], un-cc, dtdx, Ip_src_un_0, Im_src_un_0); ppm_int_profile_single(sm, sp, s[i0], un+cc, dtdx, Ip_src_un_2, Im_src_un_2); } @@ -336,7 +336,7 @@ Castro::trace_ppm(const Box& bx, if (do_trace) { load_stencil(srcQ, idir, i, j, k, QPRES, s); - ppm_reconstruct(s, i, j, k, idir, lo_symm, hi_symm, is_axisymmetric, domlo, domhi, flat, sm, sp); + ppm_reconstruct(s, i, j, k, idir, dx, lo_symm, hi_symm, is_axisymmetric, domlo, domhi, flat, sm, sp); ppm_int_profile(sm, sp, s[i0], un, cc, dtdx, Ip_src_p, Im_src_p); } @@ -353,7 +353,7 @@ Castro::trace_ppm(const Box& bx, if (do_trace) { load_stencil(srcQ, idir, i, j, k, QREINT, s); - ppm_reconstruct(s, i, j, k, idir, lo_symm, hi_symm, is_axisymmetric, domlo, domhi, flat, sm, sp); + ppm_reconstruct(s, i, j, k, idir, dx, lo_symm, hi_symm, is_axisymmetric, domlo, domhi, flat, sm, sp); ppm_int_profile(sm, sp, s[i0], un, cc, dtdx, Ip_src_rhoe, Im_src_rhoe); } @@ -370,7 +370,7 @@ Castro::trace_ppm(const Box& bx, if (do_trace) { load_stencil(srcQ, idir, i, j, k, QUT, s); - ppm_reconstruct(s, i, j, k, idir, lo_symm, hi_symm, is_axisymmetric, domlo, domhi, flat, sm, sp); + ppm_reconstruct(s, i, j, k, idir, dx, lo_symm, hi_symm, is_axisymmetric, domlo, domhi, flat, sm, sp); ppm_int_profile_single(sm, sp, s[i0], un, dtdx, Ip_src_ut_1, Im_src_ut_1); } @@ -385,7 +385,7 @@ Castro::trace_ppm(const Box& bx, if (do_trace) { load_stencil(srcQ, idir, i, j, k, QUTT, s); - ppm_reconstruct(s, i, j, k, idir, lo_symm, hi_symm, is_axisymmetric, domlo, domhi, flat, sm, sp); + ppm_reconstruct(s, i, j, k, idir, dx, lo_symm, hi_symm, is_axisymmetric, domlo, domhi, flat, sm, sp); ppm_int_profile_single(sm, sp, s[i0], un, dtdx, Ip_src_utt_1, Im_src_utt_1); } @@ -403,7 +403,7 @@ Castro::trace_ppm(const Box& bx, int n = qpassmap(ipassive); load_passive_stencil(U_arr, rho_inv_arr, idir, i, j, k, nc, s); - ppm_reconstruct(s, i, j, k, idir, lo_symm, hi_symm, is_axisymmetric, domlo, domhi, flat, sm, sp); + ppm_reconstruct(s, i, j, k, idir, dx, lo_symm, hi_symm, is_axisymmetric, domlo, domhi, flat, sm, sp); ppm_int_profile_single(sm, sp, s[i0], un, dtdx, Ip_passive, Im_passive); // Plus state on face i From e6fd73b23430b26ef600e6aa044ee40b0c0aaf18 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Fri, 26 May 2023 08:25:48 -0400 Subject: [PATCH 20/24] fix 1d cylindrical Sedov --- Exec/hydro_tests/Sedov/inputs.1d.cyl | 3 --- Source/driver/Derive.cpp | 7 +++++++ 2 files changed, 7 insertions(+), 3 deletions(-) diff --git a/Exec/hydro_tests/Sedov/inputs.1d.cyl b/Exec/hydro_tests/Sedov/inputs.1d.cyl index 746269214c..7129c7338e 100644 --- a/Exec/hydro_tests/Sedov/inputs.1d.cyl +++ b/Exec/hydro_tests/Sedov/inputs.1d.cyl @@ -75,6 +75,3 @@ problem.nsub = 10 # EOS eos.eos_assume_neutral = 1 - -# FAB FORMAT -fab.format = ASCII diff --git a/Source/driver/Derive.cpp b/Source/driver/Derive.cpp index 8b66331caf..aaf3786f48 100644 --- a/Source/driver/Derive.cpp +++ b/Source/driver/Derive.cpp @@ -960,7 +960,9 @@ extern "C" const int coord_type = geomdata.Coord(); +#if AMREX_SPACEDIM == 2 auto problo = geomdata.ProbLoArray(); +#endif amrex::ParallelFor(bx, [=] AMREX_GPU_HOST_DEVICE (int i, int j, int k) @@ -1005,6 +1007,7 @@ extern "C" der(i,j,k,0) = std::sqrt(v1*v1 + v2*v2 + v3*v3); } else if (coord_type == 1) { +#if AMREX_SPACEDIM == 2 // 2-d axisymmetric -- the coordinate ordering is r, z, phi Real r = (static_cast(i) + 0.5_rt)*dx[0] + problo[0]; @@ -1030,6 +1033,10 @@ extern "C" der(i,j,k,0) = std::sqrt(vphi_z*vphi_z + (vr_z - vz_r)*(vr_z - vz_r) + (rvphi_r/r)*(rvphi_r/r)); +#else + // for 1-d axisymmetric, we just set vorticity to 0 + der(i,j,k,0) = 0.0_rt; +#endif } else if (coord_type == 2) { // 1-d spherical -- we don't really have a vorticity in this From 4d1a8251c7509309a5bccb2c892fc474a96fe1ff Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Fri, 26 May 2023 12:18:42 -0400 Subject: [PATCH 21/24] fix the area / volume factors for 1-d cylindrical we were neglecting this case completely --- Source/driver/Castro_util.H | 23 +++++++++++++++++++---- 1 file changed, 19 insertions(+), 4 deletions(-) diff --git a/Source/driver/Castro_util.H b/Source/driver/Castro_util.H index adf873f83a..a2cb5b1b6b 100644 --- a/Source/driver/Castro_util.H +++ b/Source/driver/Castro_util.H @@ -184,8 +184,16 @@ Real volume(const int& i, const int& j, const int& k, vol = dx[0]; - } - else { + } else if (coord == 1) { + + // 1-d cylindrical + + Real rl = geomdata.ProbLo()[0] + static_cast(i) * dx[0]; + Real rr = rl + dx[0]; + + vol = M_PI * dx[0] * (rr + rl); + + } else { // Spherical @@ -253,8 +261,15 @@ Real area(const int& i, const int& j, const int& k, a = 1.0_rt; - } - else { + } else if (coord == 1) { + + // 1-D cylindrical + + Real r = geomdata.ProbLo()[0] + static_cast(i) * dx[0]; + + a = 2.0_rt * M_PI * r; + + } else { // Spherical From c989b5d74e30f03ff565ad327a8f6886e09328a9 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Tue, 30 May 2023 14:37:52 -0400 Subject: [PATCH 22/24] fix composition --- Source/hydro/Castro_ctu_hydro.cpp | 4 ++-- Source/hydro/reconstruction.H | 6 ++++++ 2 files changed, 8 insertions(+), 2 deletions(-) diff --git a/Source/hydro/Castro_ctu_hydro.cpp b/Source/hydro/Castro_ctu_hydro.cpp index 499c10409d..22cc0dcbb0 100644 --- a/Source/hydro/Castro_ctu_hydro.cpp +++ b/Source/hydro/Castro_ctu_hydro.cpp @@ -205,12 +205,12 @@ Castro::construct_ctu_hydro_source(Real time, Real dt) Array4 const U_old_arr = Sborder.array(mfi); - rho_inv.resize(qbx3, 1); + rho_inv.resize(qbx, 1); Elixir elix_rho_inv = rho_inv.elixir(); fab_size += rho_inv.nBytes(); Array4 const rho_inv_arr = rho_inv.array(); - amrex::ParallelFor(qbx3, + amrex::ParallelFor(qbx, [=] AMREX_GPU_HOST_DEVICE (int i, int j, int k) { rho_inv_arr(i,j,k) = 1.0 / U_old_arr(i,j,k,URHO); diff --git a/Source/hydro/reconstruction.H b/Source/hydro/reconstruction.H index 1569e7a011..07b47989ad 100644 --- a/Source/hydro/reconstruction.H +++ b/Source/hydro/reconstruction.H @@ -62,25 +62,31 @@ load_passive_stencil(Array4 const& U_arr, Array4 const& using namespace reconstruction; if (idir == 0) { + s[im3] = U_arr(i-3,j,k,ncomp) * rho_inv_arr(i-3,j,k); s[im2] = U_arr(i-2,j,k,ncomp) * rho_inv_arr(i-2,j,k); s[im1] = U_arr(i-1,j,k,ncomp) * rho_inv_arr(i-1,j,k); s[i0] = U_arr(i,j,k,ncomp) * rho_inv_arr(i,j,k); s[ip1] = U_arr(i+1,j,k,ncomp) * rho_inv_arr(i+1,j,k); s[ip2] = U_arr(i+2,j,k,ncomp) * rho_inv_arr(i+2,j,k); + s[ip3] = U_arr(i+3,j,k,ncomp) * rho_inv_arr(i+3,j,k); } else if (idir == 1) { + s[im3] = U_arr(i,j-3,k,ncomp) * rho_inv_arr(i,j-3,k); s[im2] = U_arr(i,j-2,k,ncomp) * rho_inv_arr(i,j-2,k); s[im1] = U_arr(i,j-1,k,ncomp) * rho_inv_arr(i,j-1,k); s[i0] = U_arr(i,j,k,ncomp) * rho_inv_arr(i,j,k); s[ip1] = U_arr(i,j+1,k,ncomp) * rho_inv_arr(i,j+1,k); s[ip2] = U_arr(i,j+2,k,ncomp) * rho_inv_arr(i,j+2,k); + s[ip3] = U_arr(i,j+3,k,ncomp) * rho_inv_arr(i,j+3,k); } else { + s[im3] = U_arr(i,j,k-3,ncomp) * rho_inv_arr(i,j,k-3); s[im2] = U_arr(i,j,k-2,ncomp) * rho_inv_arr(i,j,k-2); s[im1] = U_arr(i,j,k-1,ncomp) * rho_inv_arr(i,j,k-1); s[i0] = U_arr(i,j,k,ncomp) * rho_inv_arr(i,j,k); s[ip1] = U_arr(i,j,k+1,ncomp) * rho_inv_arr(i,j,k+1); s[ip2] = U_arr(i,j,k+2,ncomp) * rho_inv_arr(i,j,k+2); + s[ip3] = U_arr(i,j,k+3,ncomp) * rho_inv_arr(i,j,k+3); } From d074a649dbcaa9bf618c2f4ef0a6980a4f203781 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Mon, 19 Aug 2024 16:43:52 -0400 Subject: [PATCH 23/24] fix build --- Source/hydro/Castro_mol.cpp | 6 +++--- Source/hydro/ppm.H | 16 ++++++++-------- Source/hydro/trace_ppm.cpp | 6 ++---- 3 files changed, 13 insertions(+), 15 deletions(-) diff --git a/Source/hydro/Castro_mol.cpp b/Source/hydro/Castro_mol.cpp index 7f8ec56896..adee235f93 100644 --- a/Source/hydro/Castro_mol.cpp +++ b/Source/hydro/Castro_mol.cpp @@ -163,9 +163,9 @@ Castro::mol_ppm_reconstruct(const Box& bx, const auto domhi = geom.Domain().hiVect3d(); const auto dx = geom.CellSizeArray(); - - bool lo_bc_test = lo_bc[idir] == Symmetry; - bool hi_bc_test = hi_bc[idir] == Symmetry; + + bool lo_bc_test = lo_bc[idir] == amrex::PhysBCType::symmetry; + bool hi_bc_test = hi_bc[idir] == amrex::PhysBCType::symmetry; bool is_axisymmetric = geom.Coord() == 1 ; diff --git a/Source/hydro/ppm.H b/Source/hydro/ppm.H index 688341c921..d34326ea6b 100644 --- a/Source/hydro/ppm.H +++ b/Source/hydro/ppm.H @@ -139,7 +139,7 @@ ppm_reconstruct(const Real* s, Real dsvl_l = 0.0_rt; if (dsl*dsr > 0.0_rt) { Real dsc = 0.5_rt * (s[i0] - s[im2]); - dsvl_l = std::copysign(1.0_rt, dsc) * amrex::min(std::abs(dsc),amrex::min(std::abs(dsl),std::abs(dsr))); + 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]); @@ -148,7 +148,7 @@ ppm_reconstruct(const Real* s, Real dsvl_r = 0.0_rt; if (dsl*dsr > 0.0_rt) { Real dsc = 0.5_rt * (s[ip1] - s[im1]); - dsvl_r = std::copysign(1.0_rt, dsc) * amrex::min(std::abs(dsc),amrex::min(std::abs(dsl),std::abs(dsr))); + dsvl_r = std::copysign(1.0_rt, dsc) * amrex::min(std::abs(dsc), std::abs(dsl), std::abs(dsr)); } // Interpolate s to edges @@ -159,8 +159,8 @@ ppm_reconstruct(const Real* s, // Make sure sedge lies in between adjacent cell-centered values - sm = amrex::max(sm, amrex::min(s[i0], s[im1])); - sm = amrex::min(sm, amrex::max(s[i0], s[im1])); + sm = std::clamp(sm, std::min(s[i0], s[im1]), std::max(s[i0], s[im1])); + } // now we compute s_{i+1/2} -- the right interface value for zone i @@ -233,7 +233,7 @@ ppm_reconstruct(const Real* s, Real 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), amrex::min(std::abs(dsl),std::abs(dsr))); + 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]); @@ -242,7 +242,7 @@ ppm_reconstruct(const Real* s, Real 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), amrex::min(std::abs(dsl),std::abs(dsr))); + dsvl_r = std::copysign(1.0_rt, dsc) * amrex::min(std::abs(dsc), std::abs(dsl), std::abs(dsr)); } // Interpolate s to edges @@ -253,8 +253,8 @@ ppm_reconstruct(const Real* s, // Make sure sedge lies in between adjacent cell-centered values - sp = amrex::max(sp, amrex::min(s[ip1], s[i0])); - sp = amrex::min(sp, amrex::max(s[ip1], s[i0])); + sp = std::clamp(sp, std::min(s[ip1], s[i0]), std::max(s[ip1], s[i0])); + } // Flatten the parabola diff --git a/Source/hydro/trace_ppm.cpp b/Source/hydro/trace_ppm.cpp index 28c0f1d512..2407755a6f 100644 --- a/Source/hydro/trace_ppm.cpp +++ b/Source/hydro/trace_ppm.cpp @@ -152,8 +152,8 @@ Castro::trace_ppm(const Box& bx, const auto domlo = geom.Domain().loVect3d(); const auto domhi = geom.Domain().hiVect3d(); - bool lo_symm = lo_bc[idir] == Symmetry; - bool hi_symm = hi_bc[idir] == Symmetry; + bool lo_symm = lo_bc[idir] == amrex::PhysBCType::symmetry; + bool hi_symm = hi_bc[idir] == amrex::PhysBCType::symmetry; // Trace to left and right edges using upwind PPM @@ -651,5 +651,3 @@ Castro::trace_ppm(const Box& bx, }); } - - From 147854be5910a5f357e299097d025bda1b3cc8d8 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Wed, 23 Oct 2024 20:36:47 -0400 Subject: [PATCH 24/24] fix compilation --- Source/hydro/Castro_mol.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Source/hydro/Castro_mol.cpp b/Source/hydro/Castro_mol.cpp index 71b030a1f5..75d3982645 100644 --- a/Source/hydro/Castro_mol.cpp +++ b/Source/hydro/Castro_mol.cpp @@ -179,7 +179,7 @@ Castro::mol_ppm_reconstruct(const Box& bx, Real sp; load_stencil(q_arr, idir, i, j, k, n, s); - ppm_reconstruct(s, i, j, k, idir, dx, + ppm_reconstruct(s, i, j, k, idir, dx[idir], lo_bc_test, hi_bc_test, is_axisymmetric, domlo, domhi, flat, sm, sp);