diff --git a/Source/hydro/Castro_ctu_hydro.cpp b/Source/hydro/Castro_ctu_hydro.cpp index d443bbe060..f4b8e9800e 100644 --- a/Source/hydro/Castro_ctu_hydro.cpp +++ b/Source/hydro/Castro_ctu_hydro.cpp @@ -219,11 +219,11 @@ Castro::construct_ctu_hydro_source(Real time, Real dt) // NOLINT(readability-co Array4 const U_old_arr = Sborder.array(mfi); - rho_inv.resize(qbx3, 1); + rho_inv.resize(qbx, 1); fab_size += rho_inv.nBytes(); Array4 const rho_inv_arr = rho_inv.array(); - amrex::ParallelFor(qbx3, + amrex::ParallelFor(qbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { rho_inv_arr(i,j,k) = 1.0 / U_old_arr(i,j,k,URHO); @@ -472,6 +472,8 @@ Castro::construct_ctu_hydro_source(Real time, Real dt) // NOLINT(readability-co shk_arr, 0, false); + enforce_reflect_states(xbx, 0, qxm_arr, qxp_arr); + #endif // 1-d @@ -587,6 +589,8 @@ Castro::construct_ctu_hydro_source(Real time, Real dt) // NOLINT(readability-co #endif + enforce_reflect_states(xbx, 0, ql_arr, qr_arr); + cmpflx_plus_godunov(xbx, ql_arr, qr_arr, flux0_arr, @@ -631,6 +635,8 @@ Castro::construct_ctu_hydro_source(Real time, Real dt) // NOLINT(readability-co #endif #endif + enforce_reflect_states(ybx, 1, ql_arr, qr_arr); + cmpflx_plus_godunov(ybx, ql_arr, qr_arr, flux1_arr, diff --git a/Source/hydro/Castro_mol.cpp b/Source/hydro/Castro_mol.cpp index 688b298a70..75d3982645 100644 --- a/Source/hydro/Castro_mol.cpp +++ b/Source/hydro/Castro_mol.cpp @@ -155,6 +155,20 @@ 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(); + + const auto dx = geom.CellSizeArray(); + + 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 +; amrex::ParallelFor(bx, NQ, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { @@ -165,7 +179,9 @@ Castro::mol_ppm_reconstruct(const Box& bx, Real sp; load_stencil(q_arr, idir, i, j, k, n, s); - ppm_reconstruct(s, flat, sm, sp); + ppm_reconstruct(s, i, j, k, idir, dx[idir], + lo_bc_test, hi_bc_test, is_axisymmetric, 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 53202100b8..352c0cde23 100644 --- a/Source/hydro/ppm.H +++ b/Source/hydro/ppm.H @@ -48,7 +48,7 @@ check_trace_source(amrex::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[nslp] the state to be reconstructed in zones i-2, i-1, i, i+1, i+2 +/// @param s Real[nslp] 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 @@ -56,67 +56,205 @@ check_trace_source(amrex::Array4 const& srcQ, const int idir, AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void ppm_reconstruct(const amrex::Real* s, + const int i, const int j, const int k, const int idir, + const amrex::Real dx, + const bool lo_bc_test, const bool hi_bc_test, const bool is_axisymmetric, + const amrex::GpuArray& domlo, const amrex::GpuArray& domhi, const amrex::Real flatn, amrex::Real& sm, amrex::Real& sp) { + // first we compute s_{i-1/2} -- the left interface value for zone i if (ppm_do_limiting) { - // Compute van Leer slopes + if (lo_bc_test && ((idir == 0 && i == domlo[0]) || + (idir == 1 && j == domlo[1]) || + (idir == 2 && k == domlo[2]))) { - amrex::Real dsl = 2.0_rt * (s[im1] - s[im2]); - amrex::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 - amrex::Real dsvl_l = 0.0_rt; - if (dsl*dsr > 0.0_rt) { - amrex::Real dsc = 0.5_rt * (s[i0] - s[im2]); - dsvl_l = std::copysign(1.0_rt, dsc) * amrex::min(std::abs(dsc), std::abs(dsl), std::abs(dsr)); - } + 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]); + } - dsl = 2.0_rt * (s[i0] - s[im1]); - dsr = 2.0_rt * (s[ip1] - s[i0]); + } else if (lo_bc_test && ((idir == 0 && i == domlo[0]+1) || + (idir == 1 && j == domlo[1]+1) || + (idir == 2 && k == domlo[2]+1))) { - amrex::Real dsvl_r = 0.0_rt; - if (dsl*dsr > 0.0_rt) { - amrex::Real dsc = 0.5_rt * (s[ip1] - s[im1]); - dsvl_r = std::copysign(1.0_rt, dsc) * amrex::min(std::abs(dsc), std::abs(dsl),std::abs(dsr)); - } + // 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 - // Interpolate s to edges + 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 = 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 - // 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])); + } else if (hi_bc_test && ((idir == 0 && i == domhi[0]) || + (idir == 1 && j == domhi[1]) || + (idir == 2 && k == domhi[2]))) { - // Compute van Leer slopes + // 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]); - dsl = 2.0_rt * (s[i0] - s[im1]); - dsr = 2.0_rt * (s[ip1] - s[i0]); + // Make sure sedge lies in between adjacent cell-centered values - dsvl_l = 0.0_rt; - if (dsl*dsr > 0.0_rt) { - amrex::Real dsc = 0.5_rt * (s[ip1] - s[im1]); - dsvl_l = std::copysign(1.0_rt, dsc) * amrex::min(std::abs(dsc), std::abs(dsl), std::abs(dsr)); - } + sm = amrex::max(sm, amrex::min(s[i0], s[im1])); + sm = amrex::min(sm, amrex::max(s[i0], s[im1])); + + } else { + + if (idir == 0 && is_axisymmetric) { + + const amrex::Real r_m12 = (static_cast(i) - 0.5_rt) * dx; + + amrex::Real dr4_term = std::pow(dx, 4) * (-12.0_rt * s[im2] + 60.0_rt * s[im1] + 60.0_rt * s[i0] - 12.0_rt * s[ip1]); + amrex::Real dr3_term = r_m12 * std::pow(dx, 3) * (-s[im2] - 27.0_rt * s[im1] + 27.0_rt * s[i0] + s[ip1]); + amrex::Real dr2_term = std::pow(r_m12, 2) * std::pow(dx, 2) * (30.0_rt * s[im2] - 210.0_rt * s[im1] - 210.0_rt * s[i0] + 30.0_rt * s[ip1]); + amrex::Real dr1_term = std::pow(r_m12, 3) * dx * (-s[im2] - 13.0_rt * s[im1] - 13.0_rt * s[i0] + s[ip1]); + amrex::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]); + + amrex::Real denom = 24.0_rt * (4.0_rt * std::pow(dx, 4) - 15.0_rt * std::pow(dx, 2) * std::pow(r_m12, 2) + 5.0_rt * std::pow(r_m12, 4)); + + sm = (dr4_term + dr3_term + dr2_term + dr1_term + dr0_term) / denom; + + } else { + + // Compute van Leer slopes - dsl = 2.0_rt * (s[ip1] - s[i0]); - dsr = 2.0_rt * (s[ip2] - s[ip1]); + amrex::Real dsl = 2.0_rt * (s[im1] - s[im2]); + amrex::Real dsr = 2.0_rt * (s[i0] - s[im1]); + + amrex::Real dsvl_l = 0.0_rt; + if (dsl*dsr > 0.0_rt) { + amrex::Real dsc = 0.5_rt * (s[i0] - s[im2]); + dsvl_l = std::copysign(1.0_rt, dsc) * amrex::min(std::abs(dsc), std::abs(dsl), std::abs(dsr)); + } + + dsl = 2.0_rt * (s[i0] - s[im1]); + dsr = 2.0_rt * (s[ip1] - s[i0]); + + amrex::Real dsvl_r = 0.0_rt; + if (dsl*dsr > 0.0_rt) { + amrex::Real dsc = 0.5_rt * (s[ip1] - s[im1]); + dsvl_r = std::copysign(1.0_rt, dsc) * amrex::min(std::abs(dsc), std::abs(dsl), std::abs(dsr)); + } + + // 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 = std::clamp(sm, std::min(s[i0], s[im1]), std::max(s[i0], s[im1])); - dsvl_r = 0.0_rt; - if (dsl*dsr > 0.0_rt) { - amrex::Real dsc = 0.5_rt * (s[ip2] - s[i0]); - dsvl_r = std::copysign(1.0_rt, dsc) * amrex::min(std::abs(dsc), std::abs(dsl), std::abs(dsr)); } - // 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, + + 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 + + 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]) || + (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 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]); + + } 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 { + + if (idir == 0 && is_axisymmetric) { + + const amrex::Real r_p12 = (static_cast(i) + 0.5_rt) * dx; + + amrex::Real dr4_term = std::pow(dx, 4) * (-12.0_rt * s[im1] + 60.0_rt * s[i0] + 60.0_rt * s[ip1] - 12.0_rt * s[ip2]); + amrex::Real dr3_term = r_p12 * std::pow(dx, 3) * (-s[im1] - 27.0_rt * s[i0] + 27.0_rt * s[ip1] + s[ip2]); + amrex::Real dr2_term = std::pow(r_p12, 2) * std::pow(dx, 2) * (30.0_rt * s[im1] - 210.0_rt * s[i0] - 210.0_rt * s[ip1] + 30.0_rt * s[ip2]); + amrex::Real dr1_term = std::pow(r_p12, 3) * dx * (-s[im1] - 13.0_rt * s[i0] - 13.0_rt * s[ip1] + s[ip2]); + amrex::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]); - sp = 0.5_rt * (s[ip1] + s[i0]) - (1.0_rt/6.0_rt) * (dsvl_r - dsvl_l); + amrex::Real denom = 24.0_rt * (4.0_rt * std::pow(dx, 4) - 15.0_rt * std::pow(dx, 2) * std::pow(r_p12, 2) + 5.0_rt * std::pow(r_p12, 4)); - // Make sure sedge lies in between adjacent cell-centered values + sp = (dr4_term + dr3_term + dr2_term + dr1_term + dr0_term) / denom; - sp = std::clamp(sp, std::min(s[ip1], s[i0]), std::max(s[ip1], s[i0])); + } else { + // Compute van Leer slopes + + amrex::Real dsl = 2.0_rt * (s[i0] - s[im1]); + amrex::Real dsr = 2.0_rt * (s[ip1] - s[i0]); + + amrex::Real dsvl_l = 0.0_rt; + if (dsl*dsr > 0.0_rt) { + amrex::Real dsc = 0.5_rt * (s[ip1] - s[im1]); + dsvl_l = std::copysign(1.0_rt, dsc) * amrex::min(std::abs(dsc), std::abs(dsl), std::abs(dsr)); + } + + dsl = 2.0_rt * (s[ip1] - s[i0]); + dsr = 2.0_rt * (s[ip2] - s[ip1]); + + amrex::Real dsvl_r = 0.0_rt; + if (dsl*dsr > 0.0_rt) { + amrex::Real dsc = 0.5_rt * (s[ip2] - s[i0]); + dsvl_r = std::copysign(1.0_rt, dsc) * amrex::min(std::abs(dsc), std::abs(dsl), std::abs(dsr)); + } + + // 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 = std::clamp(sp, std::min(s[ip1], s[i0]), std::max(s[ip1], s[i0])); + + } // Flatten the parabola @@ -127,14 +265,14 @@ ppm_reconstruct(const amrex::Real* s, // from Colella and Sekora (2008), not the original PPM paper. if ((sp - s[i0]) * (s[i0] - sm) <= 0.0_rt) { - sp = s[i0]; - sm = s[i0]; + sp = s[i0]; + sm = s[i0]; } else if (std::abs(sp - s[i0]) >= 2.0_rt * std::abs(sm - s[i0])) { - sp = 3.0_rt * s[i0] - 2.0_rt * sm; + sp = 3.0_rt * s[i0] - 2.0_rt * sm; } else if (std::abs(sm - s[i0]) >= 2.0_rt * std::abs(sp - s[i0])) { - sm = 3.0_rt * s[i0] - 2.0_rt * sp; + sm = 3.0_rt * s[i0] - 2.0_rt * sp; } } else { @@ -170,8 +308,12 @@ ppm_reconstruct(const amrex::Real* s, /// AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE bool -ppm_reconstruct_pslope(const amrex::Real* rho, const amrex::Real* p, const amrex::Real* src, const amrex::Real flatn, +ppm_reconstruct_pslope(const amrex::Real* rho, const amrex::Real* p, const amrex::Real* src, const amrex::Real 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 amrex::GpuArray& domlo, const amrex::GpuArray& domhi, + const amrex::Real flatn, amrex::Real& sm, amrex::Real& sp) { amrex::Real tp[nslp]; @@ -182,16 +324,20 @@ ppm_reconstruct_pslope(const amrex::Real* rho, const amrex::Real* p, const amrex amrex::Real pp1_hse = p0_hse + 0.25_rt*dx * (rho[i0] + rho[ip1]) * (src[i0] + src[ip1]); amrex::Real pp2_hse = pp1_hse + 0.25_rt*dx * (rho[ip1] + rho[ip2]) * (src[ip1] + src[ip2]); + amrex::Real pp3_hse = pp2_hse + 0.25_rt*dx * (rho[ip2] + rho[ip3]) * (src[ip2] + src[ip3]); amrex::Real pm1_hse = p0_hse - 0.25_rt*dx * (rho[i0] + rho[im1]) * (src[i0] + src[im1]); amrex::Real pm2_hse = pm1_hse - 0.25_rt*dx * (rho[im1] + rho[im2]) * (src[im1] + src[im2]); + amrex::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; bool do_hse = !ptest && rho[i0] >= pslope_cutoff_density; @@ -204,22 +350,26 @@ ppm_reconstruct_pslope(const amrex::Real* rho, const amrex::Real* p, const amrex 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, 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 @@ -262,7 +412,7 @@ ppm_int_profile(const amrex::Real sm, const amrex::Real sp, const amrex::Real sc // compute x-component of Ip and Im - Real s6 = 6.0_rt * sc - 3.0_rt * (sm + sp); + amrex::Real s6 = 6.0_rt * sc - 3.0_rt * (sm + sp); // Ip/m is the integral under the parabola for the extent // that a wave can travel over a timestep @@ -271,8 +421,8 @@ ppm_int_profile(const amrex::Real sm, const amrex::Real sp, const amrex::Real sc // Im integrates to the left edge of a cell // u-c wave - Real speed = u - c; - Real sigma = std::abs(speed) * dtdx; + amrex::Real speed = u - c; + amrex::Real sigma = std::abs(speed) * dtdx; // if speed == ZERO, then either branch is the same if (speed <= 0.0_rt) { @@ -323,15 +473,15 @@ ppm_int_profile(const amrex::Real sm, const amrex::Real sp, const amrex::Real sc /// AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void -ppm_int_profile_single(const Real sm, const Real sp, const Real sc, - const Real lam, const Real dtdx, - Real& Ip, Real& Im) { +ppm_int_profile_single(const amrex::Real sm, const amrex::Real sp, const amrex::Real sc, + const amrex::Real lam, const amrex::Real dtdx, + amrex::Real& Ip, amrex::Real& Im) { // Integrate the parabolic profile to the edge of the cell. // This is the MHD version. We come in already with the eigenvalues. // compute x-component of Ip and Im - Real s6 = 6.0_rt * sc - 3.0_rt * (sm + sp); + amrex::Real s6 = 6.0_rt * sc - 3.0_rt * (sm + sp); // Ip/m is the integral under the parabola for the extent // that a wave can travel over a timestep @@ -339,7 +489,7 @@ ppm_int_profile_single(const Real sm, const Real sp, const Real sc, // Ip integrates to the right edge of a cell // Im integrates to the left edge of a cell - Real sigma = std::abs(lam) * dtdx; + amrex::Real sigma = std::abs(lam) * dtdx; if (lam <= 0.0_rt) { Ip = sp; diff --git a/Source/hydro/reconstruction.H b/Source/hydro/reconstruction.H index 8f791a1ad5..72917cb8ab 100644 --- a/Source/hydro/reconstruction.H +++ b/Source/hydro/reconstruction.H @@ -5,11 +5,13 @@ namespace reconstruction { enum slope_indices { - im2 = 0, + im3 = 0, + im2, im1, i0, ip1, ip2, + ip3, nslp }; } @@ -23,25 +25,31 @@ load_stencil(amrex::Array4 const& q_arr, const int idir, using namespace reconstruction; 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); } @@ -58,25 +66,31 @@ load_passive_stencil(amrex::Array4 const& U_arr, 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); } diff --git a/Source/hydro/trace_ppm.cpp b/Source/hydro/trace_ppm.cpp index a824fbb56f..6c05d74612 100644 --- a/Source/hydro/trace_ppm.cpp +++ b/Source/hydro/trace_ppm.cpp @@ -58,6 +58,8 @@ Castro::trace_ppm(const Box& bx, const auto problo = geom.ProbLoArray(); const int coord = geom.Coord(); + bool is_axisymmetric = geom.Coord() == 1; + Real hdt = 0.5_rt * dt; #ifndef AMREX_USE_GPU @@ -153,6 +155,17 @@ Castro::trace_ppm(const Box& bx, Real lsmall_dens = small_dens; Real lsmall_pres = small_pres; + // 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_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 amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept @@ -200,14 +213,15 @@ Castro::trace_ppm(const Box& bx, Real sm; Real sp; - // reconstruct density Real Ip_rho[3]; 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, dL, 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, dtdL, Ip_rho, Im_rho); // reconstruct normal velocity @@ -218,7 +232,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, dL, lo_symm, hi_symm, is_axisymmetric, domlo, domhi, flat, sm, sp); ppm_int_profile_single(sm, sp, s[i0], un-cc, dtdL, Ip_un_0, Im_un_0); ppm_int_profile_single(sm, sp, s[i0], un+cc, dtdL, Ip_un_2, Im_un_2); @@ -244,7 +258,7 @@ 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); - in_hse = ppm_reconstruct_pslope(trho, s, src, flat, dL, sm, sp); + in_hse = ppm_reconstruct_pslope(trho, s, src, dL, i, j, k, idir, lo_symm, hi_symm, is_axisymmetric, domlo, domhi, flat, sm, sp); if (in_hse && castro::ppm_well_balanced) { // we are working with the perturbational pressure @@ -257,7 +271,7 @@ Castro::trace_ppm(const Box& bx, } } else { - ppm_reconstruct(s, flat, sm, sp); + ppm_reconstruct(s, i, j, k, idir, dL, lo_symm, hi_symm, is_axisymmetric, domlo, domhi, flat, sm, sp); ppm_int_profile(sm, sp, s[i0], un, cc, dtdL, Ip_p, Im_p); } @@ -267,7 +281,9 @@ 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, dL, 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, dtdL, Ip_rhoe, Im_rhoe); // reconstruct transverse velocities @@ -278,11 +294,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, dL, lo_symm, hi_symm, is_axisymmetric, domlo, domhi, flat, sm, sp); ppm_int_profile_single(sm, sp, s[i0], un, dtdL, 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, dL, lo_symm, hi_symm, is_axisymmetric, domlo, domhi, flat, sm, sp); ppm_int_profile_single(sm, sp, s[i0], un, dtdL, Ip_utt_1, Im_utt_1); // gamma_c @@ -293,7 +309,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, dL, lo_symm, hi_symm, is_axisymmetric, domlo, domhi, flat, sm, sp); ppm_int_profile_single(sm, sp, s[i0], un-cc, dtdL, Ip_gc_0, Im_gc_0); ppm_int_profile_single(sm, sp, s[i0], un+cc, dtdL, Ip_gc_2, Im_gc_2); @@ -323,7 +339,7 @@ Castro::trace_ppm(const Box& bx, add_geometric_rho_source(q_arr, dloga, i, j, k, QUN, s); } #endif - ppm_reconstruct(s, flat, sm, sp); + ppm_reconstruct(s, i, j, k, idir, dL, lo_symm, hi_symm, is_axisymmetric, domlo, domhi, flat, sm, sp); ppm_int_profile(sm, sp, s[i0], un, cc, dtdL, Ip_src_rho, Im_src_rho); } @@ -342,7 +358,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, dL, lo_symm, hi_symm, is_axisymmetric, domlo, domhi, flat, sm, sp); ppm_int_profile_single(sm, sp, s[i0], un-cc, dtdL, Ip_src_un_0, Im_src_un_0); ppm_int_profile_single(sm, sp, s[i0], un+cc, dtdL, Ip_src_un_2, Im_src_un_2); } @@ -369,7 +385,7 @@ Castro::trace_ppm(const Box& bx, add_geometric_p_source(q_arr, qaux_arr, dloga, i, j, k, QUN, s); } #endif - ppm_reconstruct(s, flat, sm, sp); + ppm_reconstruct(s, i, j, k, idir, dL, lo_symm, hi_symm, is_axisymmetric, domlo, domhi, flat, sm, sp); ppm_int_profile(sm, sp, s[i0], un, cc, dtdL, Ip_src_p, Im_src_p); } @@ -395,7 +411,7 @@ Castro::trace_ppm(const Box& bx, add_geometric_rhoe_source(q_arr, dloga, i, j, k, QUN, s); } #endif - ppm_reconstruct(s, flat, sm, sp); + ppm_reconstruct(s, i, j, k, idir, dL, lo_symm, hi_symm, is_axisymmetric, domlo, domhi, flat, sm, sp); ppm_int_profile(sm, sp, s[i0], un, cc, dtdL, Ip_src_rhoe, Im_src_rhoe); } @@ -412,7 +428,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, dL, lo_symm, hi_symm, is_axisymmetric, domlo, domhi, flat, sm, sp); ppm_int_profile_single(sm, sp, s[i0], un, dtdL, Ip_src_ut_1, Im_src_ut_1); } @@ -427,7 +443,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, dL, lo_symm, hi_symm, is_axisymmetric, domlo, domhi, flat, sm, sp); ppm_int_profile_single(sm, sp, s[i0], un, dtdL, Ip_src_utt_1, Im_src_utt_1); } @@ -445,7 +461,7 @@ Castro::trace_ppm(const Box& bx, const int n = qpassmap(ipassive); load_passive_stencil(U_arr, rho_inv_arr, idir, i, j, k, nc, s); - ppm_reconstruct(s, flat, sm, sp); + ppm_reconstruct(s, i, j, k, idir, dL, lo_symm, hi_symm, is_axisymmetric, domlo, domhi, flat, sm, sp); ppm_int_profile_single(sm, sp, s[i0], un, dtdL, Ip_passive, Im_passive); // Plus state on face i