diff --git a/Src/AmrCore/AMReX_Interpolater.cpp b/Src/AmrCore/AMReX_Interpolater.cpp index 5ad0f440eff..80c00a24175 100644 --- a/Src/AmrCore/AMReX_Interpolater.cpp +++ b/Src/AmrCore/AMReX_Interpolater.cpp @@ -14,7 +14,8 @@ namespace amrex { * PCInterp, NodeBilinear, FaceLinear, CellConservativeLinear, and * CellBilinear are supported for all dimensions on cpu and gpu. * - * CellConservativeProtected only works in 2D and 3D on cpu and gpu. + * CellConservativeProtected only works in 2D and 3D on cpu and gpu + * and assumes that ratio > 1 in all directions * * CellQuadratic only works in 2D and 3D on cpu and gpu. * @@ -411,7 +412,11 @@ CellConservativeLinear::CoarseBox (const Box& fine, const IntVect& ratio) { Box crse = amrex::coarsen(fine,ratio); - crse.grow(1); + for (int dim = 0; dim < AMREX_SPACEDIM; dim++) { + if (ratio[dim] > 1) { + crse.grow(dim,1); + } + } return crse; } @@ -440,7 +445,7 @@ CellConservativeLinear::interp (const FArrayBox& crse, RunOn runon) { BL_PROFILE("CellConservativeLinear::interp()"); - BL_ASSERT(bcr.size() >= ncomp); + AMREX_ASSERT(bcr.size() >= ncomp); AMREX_ASSERT(fine.box().contains(fine_region)); @@ -454,7 +459,12 @@ CellConservativeLinear::interp (const FArrayBox& crse, Array4 const& finearr = fine.array(); const Box& crse_region = CoarseBox(fine_region,ratio); - const Box& cslope_bx = amrex::grow(crse_region,-1); + Box cslope_bx(crse_region); + for (int dim = 0; dim < AMREX_SPACEDIM; dim++) { + if (ratio[dim] > 1) { + cslope_bx.grow(dim,-1); + } + } FArrayBox ccfab(cslope_bx, ncomp*AMREX_SPACEDIM); Array4 const& tmp = ccfab.array(); @@ -478,7 +488,7 @@ CellConservativeLinear::interp (const FArrayBox& crse, AMREX_HOST_DEVICE_PARALLEL_FOR_3D_FLAG(runon, cslope_bx, i, j, k, { mf_cell_cons_lin_interp_llslope(i,j,k, tmp, crsearr, crse_comp, ncomp, - cdomain, bcrp); + cdomain, ratio, bcrp); }); } else { AMREX_HOST_DEVICE_PARALLEL_FOR_4D_FLAG(runon, cslope_bx, ncomp, i, j, k, n, @@ -504,7 +514,7 @@ CellConservativeLinear::interp (const FArrayBox& crse, AMREX_HOST_DEVICE_PARALLEL_FOR_3D_FLAG(runon, cslope_bx, i, j, k, { mf_cell_cons_lin_interp_llslope(i,j,k, tmp, crsearr, crse_comp, ncomp, - cdomain, bcrp); + cdomain, ratio, bcrp); }); } else { AMREX_HOST_DEVICE_PARALLEL_FOR_4D_FLAG(runon, cslope_bx, ncomp, i, j, k, n, @@ -528,7 +538,7 @@ CellConservativeLinear::interp (const FArrayBox& crse, AMREX_HOST_DEVICE_PARALLEL_FOR_3D_FLAG(runon, cslope_bx, i, j, k, { mf_cell_cons_lin_interp_llslope(i,j,k, tmp, crsearr, crse_comp, ncomp, - cdomain, bcrp); + cdomain, ratio, bcrp); }); } else { AMREX_HOST_DEVICE_PARALLEL_FOR_4D_FLAG(runon, cslope_bx, ncomp, i, j, k, n, @@ -586,7 +596,7 @@ CellQuadratic::interp (const FArrayBox& crse, #else BL_PROFILE("CellQuadratic::interp()"); - BL_ASSERT(bcr.size() >= ncomp); + AMREX_ASSERT(bcr.size() >= ncomp); // // Make box which is intersection of fine_region and domain of fine. @@ -595,7 +605,7 @@ CellQuadratic::interp (const FArrayBox& crse, // Make Box for slopes. Box cslope_bx = amrex::coarsen(target_fine_region,ratio); - BL_ASSERT(crse.box().contains(cslope_bx)); + AMREX_ASSERT(crse.box().contains(cslope_bx)); // Are we running on GPU? bool run_on_gpu = (runon == RunOn::Gpu && Gpu::inLaunchRegion()); @@ -739,6 +749,8 @@ CellConservativeProtected::protect (const FArrayBox& /*crse*/, Vector& /*bcr*/, RunOn runon) { + AMREX_ALWAYS_ASSERT(ratio.allGT(IntVect(1))); + #if (AMREX_SPACEDIM == 1) amrex::ignore_unused(fine,fine_state, ncomp,fine_region,ratio, @@ -838,13 +850,7 @@ CellConservativeQuartic::interp (const FArrayBox& crse, RunOn runon) { BL_PROFILE("CellConservativeQuartic::interp()"); - BL_ASSERT(ratio[0] == 2); -#if (AMREX_SPACEDIM >= 2) - BL_ASSERT(ratio[0] == ratio[1]); -#endif -#if (AMREX_SPACEDIM == 3) - BL_ASSERT(ratio[1] == ratio[2]); -#endif + AMREX_ASSERT(ratio == 2); amrex::ignore_unused(ratio); // @@ -867,16 +873,16 @@ Box FaceDivFree::CoarseBox (const Box& fine, int ratio) { - Box b = amrex::coarsen(fine,ratio).grow(1); - return b; + Box crse = amrex::coarsen(fine,ratio).grow(1); + return crse; } Box FaceDivFree::CoarseBox (const Box& fine, const IntVect& ratio) { - Box b = amrex::coarsen(fine,ratio).grow(1); - return b; + Box crse = amrex::coarsen(fine,ratio).grow(1); + return crse; } void @@ -894,7 +900,7 @@ FaceDivFree::interp (const FArrayBox& /*crse*/, int /*actual_state*/, RunOn /*runon*/) { - amrex::Abort("FaceDivFree does not work on a single MultiFab. Call 'interp_arr' instead."); + amrex::Abort("FaceDivFree does not work on a single FArrayBox. Call 'interp_arr' instead."); } void diff --git a/Src/AmrCore/AMReX_MFInterp_1D_C.H b/Src/AmrCore/AMReX_MFInterp_1D_C.H index d1ab2765d12..2d52ca6d57f 100644 --- a/Src/AmrCore/AMReX_MFInterp_1D_C.H +++ b/Src/AmrCore/AMReX_MFInterp_1D_C.H @@ -38,7 +38,7 @@ void mf_cell_cons_lin_interp_limit_minmax_llslope (int i, int, int, Array4 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mf_cell_cons_lin_interp_llslope (int i, int, int, Array4 const& slope, Array4 const& u, int scomp, int ncomp, - Box const& domain, BCRec const* bc) noexcept + Box const& domain, IntVect const& /*ratio*/, BCRec const* bc) noexcept { Real sfx = Real(1.0); diff --git a/Src/AmrCore/AMReX_MFInterp_2D_C.H b/Src/AmrCore/AMReX_MFInterp_2D_C.H index 4c6fa85768c..c51c37e4419 100644 --- a/Src/AmrCore/AMReX_MFInterp_2D_C.H +++ b/Src/AmrCore/AMReX_MFInterp_2D_C.H @@ -13,24 +13,38 @@ void mf_cell_cons_lin_interp_limit_minmax_llslope (int i, int j, int, Array4= Real(0.0)) ? amrex::min(std::abs(df),std::abs(db)) : Real(0.); - sx = std::copysign(Real(1.),dcx)*amrex::min(sx,std::abs(dcx)); - slope(i,j,0,ns ) = dcx; + if (ratio[0] > 1) { + dcx = mf_compute_slopes_x(i, j, 0, u, nu, domain, bc[ns]); + Real df = Real(2.0) * (u(i+1,j,0,nu) - u(i ,j,0,nu)); + Real db = Real(2.0) * (u(i ,j,0,nu) - u(i-1,j,0,nu)); + sx = (df*db >= Real(0.0)) ? amrex::min(std::abs(df),std::abs(db)) : Real(0.); + sx = std::copysign(Real(1.),dcx)*amrex::min(sx,std::abs(dcx)); + slope(i,j,0,ns ) = dcx; + } else { + slope(i,j,0,ns ) = Real(0.0); + } // y-direction - Real dcy = mf_compute_slopes_y(i, j, 0, u, nu, domain, bc[ns]); - df = Real(2.0) * (u(i,j+1,0,nu) - u(i,j ,0,nu)); - db = Real(2.0) * (u(i,j ,0,nu) - u(i,j-1,0,nu)); - Real sy = (df*db >= Real(0.0)) ? amrex::min(std::abs(df),std::abs(db)) : Real(0.); - sy = std::copysign(Real(1.),dcy)*amrex::min(sy,std::abs(dcy)); - slope(i,j,0,ns+ ncomp) = dcy; + if (ratio[1] > 1) { + dcy = mf_compute_slopes_y(i, j, 0, u, nu, domain, bc[ns]); + Real df = Real(2.0) * (u(i,j+1,0,nu) - u(i,j ,0,nu)); + Real db = Real(2.0) * (u(i,j ,0,nu) - u(i,j-1,0,nu)); + sy = (df*db >= Real(0.0)) ? amrex::min(std::abs(df),std::abs(db)) : Real(0.); + sy = std::copysign(Real(1.),dcy)*amrex::min(sy,std::abs(dcy)); + slope(i,j,0,ns+ ncomp) = dcy; + } else { + slope(i,j,0,ns+ ncomp) = Real(0.0); + } // adjust limited slopes to prevent new min/max for this component Real alpha = Real(1.0); @@ -39,8 +53,10 @@ void mf_cell_cons_lin_interp_limit_minmax_llslope (int i, int j, int, Array4 1 ? 1 : 0; + int jlim = ratio[1] > 1 ? 1 : 0; + for (int joff = -jlim; joff <= jlim; ++joff) { + for (int ioff = -ilim; ioff <= ilim; ++ioff) { umin = amrex::min(umin, u(i+ioff,j+joff,0,nu)); umax = amrex::max(umax, u(i+ioff,j+joff,0,nu)); }} @@ -73,7 +89,7 @@ void mf_cell_cons_lin_interp_limit_minmax_llslope (int i, int j, int, Array4 const& slope, Array4 const& u, int scomp, int ncomp, - Box const& domain, BCRec const* bc) noexcept + Box const& domain, IntVect const& ratio, BCRec const* bc) noexcept { Real sfx = Real(1.0); Real sfy = Real(1.0); @@ -82,26 +98,34 @@ void mf_cell_cons_lin_interp_llslope (int i, int j, int, Array4 const& slo int nu = ns + scomp; // x-direction - Real dc = mf_compute_slopes_x(i, j, 0, u, nu, domain, bc[ns]); - Real df = Real(2.0) * (u(i+1,j,0,nu) - u(i ,j,0,nu)); - Real db = Real(2.0) * (u(i ,j,0,nu) - u(i-1,j,0,nu)); - Real sx = (df*db >= Real(0.0)) ? amrex::min(std::abs(df),std::abs(db)) : Real(0.); - sx = std::copysign(Real(1.),dc)*amrex::min(sx,std::abs(dc)); - if (dc != Real(0.0)) { - sfx = amrex::min(sfx, sx / dc); + if (ratio[0] > 1) { + Real dc = mf_compute_slopes_x(i, j, 0, u, nu, domain, bc[ns]); + Real df = Real(2.0) * (u(i+1,j,0,nu) - u(i ,j,0,nu)); + Real db = Real(2.0) * (u(i ,j,0,nu) - u(i-1,j,0,nu)); + Real sx = (df*db >= Real(0.0)) ? amrex::min(std::abs(df),std::abs(db)) : Real(0.); + sx = std::copysign(Real(1.),dc)*amrex::min(sx,std::abs(dc)); + if (dc != Real(0.0)) { + sfx = amrex::min(sfx, sx / dc); + } + slope(i,j,0,ns ) = dc; + } else { + slope(i,j,0,ns ) = Real(0.0); } - slope(i,j,0,ns ) = dc; // y-direction - dc = mf_compute_slopes_y(i, j, 0, u, nu, domain, bc[ns]); - df = Real(2.0) * (u(i,j+1,0,nu) - u(i,j ,0,nu)); - db = Real(2.0) * (u(i,j ,0,nu) - u(i,j-1,0,nu)); - Real sy = (df*db >= Real(0.0)) ? amrex::min(std::abs(df),std::abs(db)) : Real(0.); - sy = std::copysign(Real(1.),dc)*amrex::min(sy,std::abs(dc)); - if (dc != Real(0.0)) { - sfy = amrex::min(sfy, sy / dc); + if (ratio[1] > 1) { + Real dc = mf_compute_slopes_y(i, j, 0, u, nu, domain, bc[ns]); + Real df = Real(2.0) * (u(i,j+1,0,nu) - u(i,j ,0,nu)); + Real db = Real(2.0) * (u(i,j ,0,nu) - u(i,j-1,0,nu)); + Real sy = (df*db >= Real(0.0)) ? amrex::min(std::abs(df),std::abs(db)) : Real(0.); + sy = std::copysign(Real(1.),dc)*amrex::min(sy,std::abs(dc)); + if (dc != Real(0.0)) { + sfy = amrex::min(sfy, sy / dc); + } + slope(i,j,0,ns+ ncomp) = dc; + } else { + slope(i,j,0,ns+ ncomp) = Real(0.0); } - slope(i,j,0,ns+ ncomp) = dc; } for (int ns = 0; ns < ncomp; ++ns) { @@ -118,19 +142,26 @@ void mf_cell_cons_lin_interp_mcslope (int i, int j, int /*k*/, int ns, Array4= Real(0.0)) ? amrex::min(std::abs(df),std::abs(db)) : Real(0.); - sx = std::copysign(Real(1.),dc)*amrex::min(sx,std::abs(dc)); + if (ratio[0] > 1) { + Real dc = mf_compute_slopes_x(i, j, 0, u, nu, domain, bc[ns]); + Real df = Real(2.0) * (u(i+1,j,0,nu) - u(i ,j,0,nu)); + Real db = Real(2.0) * (u(i ,j,0,nu) - u(i-1,j,0,nu)); + sx = (df*db >= Real(0.0)) ? amrex::min(std::abs(df),std::abs(db)) : Real(0.); + sx = std::copysign(Real(1.),dc)*amrex::min(sx,std::abs(dc)); + } // y-direction - dc = mf_compute_slopes_y(i, j, 0, u, nu, domain, bc[ns]); - df = Real(2.0) * (u(i,j+1,0,nu) - u(i,j ,0,nu)); - db = Real(2.0) * (u(i,j ,0,nu) - u(i,j-1,0,nu)); - Real sy = (df*db >= Real(0.0)) ? amrex::min(std::abs(df),std::abs(db)) : Real(0.); - sy = std::copysign(Real(1.),dc)*amrex::min(sy,std::abs(dc)); + if (ratio[1] > 1) { + Real dc = mf_compute_slopes_y(i, j, 0, u, nu, domain, bc[ns]); + Real df = Real(2.0) * (u(i,j+1,0,nu) - u(i,j ,0,nu)); + Real db = Real(2.0) * (u(i,j ,0,nu) - u(i,j-1,0,nu)); + sy = (df*db >= Real(0.0)) ? amrex::min(std::abs(df),std::abs(db)) : Real(0.); + sy = std::copysign(Real(1.),dc)*amrex::min(sy,std::abs(dc)); + } Real alpha = Real(1.0); if (sx != Real(0.0) || sy != Real(0.0)) { @@ -138,8 +169,10 @@ void mf_cell_cons_lin_interp_mcslope (int i, int j, int /*k*/, int ns, Array4 1 ? 1 : 0; + int jlim = ratio[1] > 1 ? 1 : 0; + for (int joff = -jlim; joff <= jlim; ++joff) { + for (int ioff = -ilim; ioff <= ilim; ++ioff) { umin = amrex::min(umin, u(i+ioff,j+joff,0,nu)); umax = amrex::max(umax, u(i+ioff,j+joff,0,nu)); }} @@ -177,19 +210,26 @@ void mf_cell_cons_lin_interp_mcslope_rz (int i, int j, int ns, Array4 cons { int nu = ns + scomp; + Real sx = Real(0.0); + Real sy = Real(0.0); + // x-direction - Real dc = mf_compute_slopes_x(i, j, 0, u, nu, domain, bc[ns]); - Real df = Real(2.0) * (u(i+1,j,0,nu) - u(i ,j,0,nu)); - Real db = Real(2.0) * (u(i ,j,0,nu) - u(i-1,j,0,nu)); - Real sx = (df*db >= Real(0.0)) ? amrex::min(std::abs(df),std::abs(db)) : Real(0.); - sx = std::copysign(Real(1.),dc)*amrex::min(sx,std::abs(dc)); + if (ratio[0] > 1) { + Real dc = mf_compute_slopes_x(i, j, 0, u, nu, domain, bc[ns]); + Real df = Real(2.0) * (u(i+1,j,0,nu) - u(i ,j,0,nu)); + Real db = Real(2.0) * (u(i ,j,0,nu) - u(i-1,j,0,nu)); + sx = (df*db >= Real(0.0)) ? amrex::min(std::abs(df),std::abs(db)) : Real(0.); + sx = std::copysign(Real(1.),dc)*amrex::min(sx,std::abs(dc)); + } // y-direction - dc = mf_compute_slopes_y(i, j, 0, u, nu, domain, bc[ns]); - df = Real(2.0) * (u(i,j+1,0,nu) - u(i,j ,0,nu)); - db = Real(2.0) * (u(i,j ,0,nu) - u(i,j-1,0,nu)); - Real sy = (df*db >= Real(0.0)) ? amrex::min(std::abs(df),std::abs(db)) : Real(0.); - sy = std::copysign(Real(1.),dc)*amrex::min(sy,std::abs(dc)); + if (ratio[1] > 1) { + Real dc = mf_compute_slopes_y(i, j, 0, u, nu, domain, bc[ns]); + Real df = Real(2.0) * (u(i,j+1,0,nu) - u(i,j ,0,nu)); + Real db = Real(2.0) * (u(i,j ,0,nu) - u(i,j-1,0,nu)); + sy = (df*db >= Real(0.0)) ? amrex::min(std::abs(df),std::abs(db)) : Real(0.); + sy = std::copysign(Real(1.),dc)*amrex::min(sy,std::abs(dc)); + } Real alpha = Real(1.0); if (sx != Real(0.0) || sy != Real(0.0)) { @@ -214,8 +254,10 @@ void mf_cell_cons_lin_interp_mcslope_rz (int i, int j, int ns, Array4 cons + std::abs(sy) * Real(ratio[1]-1)/Real(2*ratio[1]); Real umax = u(i,j,0,nu); Real umin = u(i,j,0,nu); - for (int joff = -1; joff <= 1; ++joff) { - for (int ioff = -1; ioff <= 1; ++ioff) { + int ilim = ratio[0] > 1 ? 1 : 0; + int jlim = ratio[1] > 1 ? 1 : 0; + for (int joff = -jlim; joff <= jlim; ++joff) { + for (int ioff = -ilim; ioff <= ilim; ++ioff) { umin = amrex::min(umin, u(i+ioff,j+joff,0,nu)); umax = amrex::max(umax, u(i+ioff,j+joff,0,nu)); }} diff --git a/Src/AmrCore/AMReX_MFInterp_3D_C.H b/Src/AmrCore/AMReX_MFInterp_3D_C.H index db18e1adb6e..4d2e92ae3ed 100644 --- a/Src/AmrCore/AMReX_MFInterp_3D_C.H +++ b/Src/AmrCore/AMReX_MFInterp_3D_C.H @@ -12,32 +12,52 @@ void mf_cell_cons_lin_interp_limit_minmax_llslope (int i, int j, int k, Array4= Real(0.0)) ? amrex::min(std::abs(df),std::abs(db)) : Real(0.); - sx = std::copysign(Real(1.),dcx)*amrex::min(sx,std::abs(dcx)); - slope(i,j,k,ns ) = dcx; // unlimited slope + if (ratio[0] > 1) { + dcx = mf_compute_slopes_x(i, j, k, u, nu, domain, bc[ns]); + Real df = Real(2.0) * (u(i+1,j,k,nu) - u(i ,j,k,nu)); + Real db = Real(2.0) * (u(i ,j,k,nu) - u(i-1,j,k,nu)); + sx = (df*db >= Real(0.0)) ? amrex::min(std::abs(df),std::abs(db)) : Real(0.); + sx = std::copysign(Real(1.),dcx)*amrex::min(sx,std::abs(dcx)); + slope(i,j,k,ns ) = dcx; // unlimited slope + } else { + slope(i,j,k,ns ) = Real(0.0); + } // y-direction - Real dcy = mf_compute_slopes_y(i, j, k, u, nu, domain, bc[ns]); - df = Real(2.0) * (u(i,j+1,k,nu) - u(i,j ,k,nu)); - db = Real(2.0) * (u(i,j ,k,nu) - u(i,j-1,k,nu)); - Real sy = (df*db >= Real(0.0)) ? amrex::min(std::abs(df),std::abs(db)) : Real(0.); - sy = std::copysign(Real(1.),dcy)*amrex::min(sy,std::abs(dcy)); - slope(i,j,k,ns+ ncomp) = dcy; // unlimited slope + if (ratio[1] > 1) { + dcy = mf_compute_slopes_y(i, j, k, u, nu, domain, bc[ns]); + Real df = Real(2.0) * (u(i,j+1,k,nu) - u(i,j ,k,nu)); + Real db = Real(2.0) * (u(i,j ,k,nu) - u(i,j-1,k,nu)); + sy = (df*db >= Real(0.0)) ? amrex::min(std::abs(df),std::abs(db)) : Real(0.); + sy = std::copysign(Real(1.),dcy)*amrex::min(sy,std::abs(dcy)); + slope(i,j,k,ns+ ncomp) = dcy; // unlimited slope + } else { + slope(i,j,k,ns+ ncomp) = Real(0.0); + } // z-direction - Real dcz = mf_compute_slopes_z(i, j, k, u, nu, domain, bc[ns]); - df = Real(2.0) * (u(i,j,k+1,nu) - u(i,j,k ,nu)); - db = Real(2.0) * (u(i,j,k ,nu) - u(i,j,k-1,nu)); - Real sz = (df*db >= Real(0.0)) ? amrex::min(std::abs(df),std::abs(db)) : Real(0.); - sz = std::copysign(Real(1.),dcz)*amrex::min(sz,std::abs(dcz)); - slope(i,j,k,ns+2*ncomp) = dcz; // unlimited slope + if (ratio[2] > 1) { + dcz = mf_compute_slopes_z(i, j, k, u, nu, domain, bc[ns]); + Real df = Real(2.0) * (u(i,j,k+1,nu) - u(i,j,k ,nu)); + Real db = Real(2.0) * (u(i,j,k ,nu) - u(i,j,k-1,nu)); + sz = (df*db >= Real(0.0)) ? amrex::min(std::abs(df),std::abs(db)) : Real(0.); + sz = std::copysign(Real(1.),dcz)*amrex::min(sz,std::abs(dcz)); + slope(i,j,k,ns+2*ncomp) = dcz; // unlimited slope + } else { + slope(i,j,k,ns+2*ncomp) = Real(0.0); + } // adjust limited slopes to prevent new min/max for this component Real alpha = 1.0; @@ -47,9 +67,12 @@ void mf_cell_cons_lin_interp_limit_minmax_llslope (int i, int j, int k, Array4 1 ? 1 : 0; + int jlim = ratio[1] > 1 ? 1 : 0; + int klim = ratio[2] > 1 ? 1 : 0; + for (int koff = -klim; koff <= klim; ++koff) { + for (int joff = -jlim; joff <= jlim; ++joff) { + for (int ioff = -ilim; ioff <= ilim; ++ioff) { umin = amrex::min(umin, u(i+ioff,j+joff,k+koff,nu)); umax = amrex::max(umax, u(i+ioff,j+joff,k+koff,nu)); }}} @@ -87,7 +110,7 @@ void mf_cell_cons_lin_interp_limit_minmax_llslope (int i, int j, int k, Array4 const& slope, Array4 const& u, int scomp, int ncomp, - Box const& domain, BCRec const* bc) noexcept + Box const& domain, IntVect const& ratio, BCRec const* bc) noexcept { Real sfx = Real(1.0); Real sfy = Real(1.0); @@ -97,37 +120,51 @@ void mf_cell_cons_lin_interp_llslope (int i, int j, int k, Array4 const& s int nu = ns + scomp; // x-direction - Real dc = mf_compute_slopes_x(i, j, k, u, nu, domain, bc[ns]); - Real df = Real(2.0) * (u(i+1,j,k,nu) - u(i ,j,k,nu)); - Real db = Real(2.0) * (u(i ,j,k,nu) - u(i-1,j,k,nu)); - Real sx = (df*db >= Real(0.0)) ? amrex::min(std::abs(df),std::abs(db)) : Real(0.); - sx = std::copysign(Real(1.),dc)*amrex::min(sx,std::abs(dc)); - if (dc != Real(0.0)) { - sfx = amrex::min(sfx, sx / dc); + if (ratio[0] > 1) { + Real dc = mf_compute_slopes_x(i, j, k, u, nu, domain, bc[ns]); + Real df = Real(2.0) * (u(i+1,j,k,nu) - u(i ,j,k,nu)); + Real db = Real(2.0) * (u(i ,j,k,nu) - u(i-1,j,k,nu)); + Real sx = (df*db >= Real(0.0)) ? amrex::min(std::abs(df),std::abs(db)) : Real(0.); + sx = std::copysign(Real(1.),dc)*amrex::min(sx,std::abs(dc)); + if (dc != Real(0.0)) { + sfx = amrex::min(sfx, sx / dc); + } + slope(i,j,k,ns ) = dc; + } else { + slope(i,j,k,ns ) = Real(0.0); } - slope(i,j,k,ns ) = dc; // y-direction - dc = mf_compute_slopes_y(i, j, k, u, nu, domain, bc[ns]); - df = Real(2.0) * (u(i,j+1,k,nu) - u(i,j ,k,nu)); - db = Real(2.0) * (u(i,j ,k,nu) - u(i,j-1,k,nu)); - Real sy = (df*db >= Real(0.0)) ? amrex::min(std::abs(df),std::abs(db)) : Real(0.); - sy = std::copysign(Real(1.),dc)*amrex::min(sy,std::abs(dc)); - if (dc != Real(0.0)) { - sfy = amrex::min(sfy, sy / dc); + if (ratio[1] > 1) { + Real dc = mf_compute_slopes_y(i, j, k, u, nu, domain, bc[ns]); + Real df = Real(2.0) * (u(i,j+1,k,nu) - u(i,j ,k,nu)); + Real db = Real(2.0) * (u(i,j ,k,nu) - u(i,j-1,k,nu)); + Real sy = (df*db >= Real(0.0)) ? amrex::min(std::abs(df),std::abs(db)) : Real(0.); + sy = std::copysign(Real(1.),dc)*amrex::min(sy,std::abs(dc)); + if (dc != Real(0.0)) { + sfy = amrex::min(sfy, sy / dc); + } + slope(i,j,k,ns+ ncomp) = dc; + } else { + slope(i,j,k,ns+ ncomp) = Real(0.0); } - slope(i,j,k,ns+ ncomp) = dc; + // z-direction - dc = mf_compute_slopes_z(i, j, k, u, nu, domain, bc[ns]); - df = Real(2.0) * (u(i,j,k+1,nu) - u(i,j,k ,nu)); - db = Real(2.0) * (u(i,j,k ,nu) - u(i,j,k-1,nu)); - Real sz = (df*db >= Real(0.0)) ? amrex::min(std::abs(df),std::abs(db)) : Real(0.); - sz = std::copysign(Real(1.),dc)*amrex::min(sz,std::abs(dc)); - if (dc != Real(0.0)) { - sfz = amrex::min(sfz, sz / dc); + if (ratio[2] > 1) { + Real dc = mf_compute_slopes_z(i, j, k, u, nu, domain, bc[ns]); + Real df = Real(2.0) * (u(i,j,k+1,nu) - u(i,j,k ,nu)); + Real db = Real(2.0) * (u(i,j,k ,nu) - u(i,j,k-1,nu)); + Real sz = (df*db >= Real(0.0)) ? amrex::min(std::abs(df),std::abs(db)) : Real(0.); + sz = std::copysign(Real(1.),dc)*amrex::min(sz,std::abs(dc)); + if (dc != Real(0.0)) { + sfz = amrex::min(sfz, sz / dc); + } + slope(i,j,k,ns+2*ncomp) = dc; + } else { + slope(i,j,k,ns+2*ncomp) = Real(0.0); } - slope(i,j,k,ns+2*ncomp) = dc; + } for (int ns = 0; ns < ncomp; ++ns) { @@ -145,26 +182,36 @@ void mf_cell_cons_lin_interp_mcslope (int i, int j, int k, int ns, Array4 { int nu = ns + scomp; + Real sx = Real(0.); + Real sy = Real(0.); + Real sz = Real(0.); + // x-direction - Real dc = mf_compute_slopes_x(i, j, k, u, nu, domain, bc[ns]); - Real df = Real(2.0) * (u(i+1,j,k,nu) - u(i ,j,k,nu)); - Real db = Real(2.0) * (u(i ,j,k,nu) - u(i-1,j,k,nu)); - Real sx = (df*db >= Real(0.0)) ? amrex::min(std::abs(df),std::abs(db)) : Real(0.); - sx = std::copysign(Real(1.),dc)*amrex::min(sx,std::abs(dc)); + if (ratio[0] > 1) { + Real dc = mf_compute_slopes_x(i, j, k, u, nu, domain, bc[ns]); + Real df = Real(2.0) * (u(i+1,j,k,nu) - u(i ,j,k,nu)); + Real db = Real(2.0) * (u(i ,j,k,nu) - u(i-1,j,k,nu)); + sx = (df*db >= Real(0.0)) ? amrex::min(std::abs(df),std::abs(db)) : Real(0.); + sx = std::copysign(Real(1.),dc)*amrex::min(sx,std::abs(dc)); + } // y-direction - dc = mf_compute_slopes_y(i, j, k, u, nu, domain, bc[ns]); - df = Real(2.0) * (u(i,j+1,k,nu) - u(i,j ,k,nu)); - db = Real(2.0) * (u(i,j ,k,nu) - u(i,j-1,k,nu)); - Real sy = (df*db >= Real(0.0)) ? amrex::min(std::abs(df),std::abs(db)) : Real(0.); - sy = std::copysign(Real(1.),dc)*amrex::min(sy,std::abs(dc)); + if (ratio[1] > 1) { + Real dc = mf_compute_slopes_y(i, j, k, u, nu, domain, bc[ns]); + Real df = Real(2.0) * (u(i,j+1,k,nu) - u(i,j ,k,nu)); + Real db = Real(2.0) * (u(i,j ,k,nu) - u(i,j-1,k,nu)); + sy = (df*db >= Real(0.0)) ? amrex::min(std::abs(df),std::abs(db)) : Real(0.); + sy = std::copysign(Real(1.),dc)*amrex::min(sy,std::abs(dc)); + } // z-direction - dc = mf_compute_slopes_z(i, j, k, u, nu, domain, bc[ns]); - df = Real(2.0) * (u(i,j,k+1,nu) - u(i,j,k ,nu)); - db = Real(2.0) * (u(i,j,k ,nu) - u(i,j,k-1,nu)); - Real sz = (df*db >= Real(0.0)) ? amrex::min(std::abs(df),std::abs(db)) : Real(0.); - sz = std::copysign(Real(1.),dc)*amrex::min(sz,std::abs(dc)); + if (ratio[2] > 1) { + Real dc = mf_compute_slopes_z(i, j, k, u, nu, domain, bc[ns]); + Real df = Real(2.0) * (u(i,j,k+1,nu) - u(i,j,k ,nu)); + Real db = Real(2.0) * (u(i,j,k ,nu) - u(i,j,k-1,nu)); + sz = (df*db >= Real(0.0)) ? amrex::min(std::abs(df),std::abs(db)) : Real(0.); + sz = std::copysign(Real(1.),dc)*amrex::min(sz,std::abs(dc)); + } Real alpha = 1.0; if (sx != Real(0.0) || sy != Real(0.0) || sz != Real(0.0)) { @@ -173,9 +220,12 @@ void mf_cell_cons_lin_interp_mcslope (int i, int j, int k, int ns, Array4 + std::abs(sz) * Real(ratio[2]-1)/Real(2*ratio[2]); Real umax = u(i,j,k,nu); Real umin = u(i,j,k,nu); - for (int koff = -1; koff <= 1; ++koff) { - for (int joff = -1; joff <= 1; ++joff) { - for (int ioff = -1; ioff <= 1; ++ioff) { + int ilim = ratio[0] > 1 ? 1 : 0; + int jlim = ratio[1] > 1 ? 1 : 0; + int klim = ratio[2] > 1 ? 1 : 0; + for (int koff = -klim; koff <= klim; ++koff) { + for (int joff = -jlim; joff <= jlim; ++joff) { + for (int ioff = -ilim; ioff <= ilim; ++ioff) { umin = amrex::min(umin, u(i+ioff,j+joff,k+koff,nu)); umax = amrex::max(umax, u(i+ioff,j+joff,k+koff,nu)); }}} diff --git a/Src/AmrCore/AMReX_MFInterpolater.cpp b/Src/AmrCore/AMReX_MFInterpolater.cpp index 12025c0747a..21e99e760f2 100644 --- a/Src/AmrCore/AMReX_MFInterpolater.cpp +++ b/Src/AmrCore/AMReX_MFInterpolater.cpp @@ -72,7 +72,11 @@ Box MFCellConsLinInterp::CoarseBox (const Box& fine, const IntVect& ratio) { Box crse = amrex::coarsen(fine,ratio); - crse.grow(1); + for (int dim = 0; dim < AMREX_SPACEDIM; dim++) { + if (ratio[dim] > 1) { + crse.grow(dim,1); + } + } return crse; } @@ -95,6 +99,11 @@ MFCellConsLinInterp::interp (MultiFab const& crsemf, int ccomp, MultiFab& finemf Box const& cdomain = cgeom.Domain(); + IntVect minus1; + for (int dim = 0; dim < AMREX_SPACEDIM; dim++) { + minus1[dim] = (ratio[dim] > 1) ? -1 : 0; + } + #ifdef AMREX_USE_GPU if (Gpu::inLaunchRegion()) { MultiFab crse_tmp(crsemf.boxArray(), crsemf.DistributionMap(), AMREX_SPACEDIM*nc, 0); @@ -112,14 +121,14 @@ MFCellConsLinInterp::interp (MultiFab const& crsemf, int ccomp, MultiFab& finemf Real drf = fgeom.CellSize(0); Real rlo = fgeom.Offset(0); if (do_linear_limiting) { - ParallelFor(crsemf, IntVect(-1), + ParallelFor(crsemf, minus1, [=] AMREX_GPU_DEVICE (int box_no, int i, int, int) noexcept { mf_cell_cons_lin_interp_llslope(i,0,0, tmp[box_no], crse[box_no], ccomp, nc, - cdomain, pbc); + cdomain, ratio, pbc); }); } else { - ParallelFor(crsemf, IntVect(-1), nc, + ParallelFor(crsemf, minus1, nc, [=] AMREX_GPU_DEVICE (int box_no, int i, int, int, int n) noexcept { mf_cell_cons_lin_interp_mcslope_sph(i, n, tmp[box_no], crse[box_no], ccomp, nc, @@ -141,14 +150,14 @@ MFCellConsLinInterp::interp (MultiFab const& crsemf, int ccomp, MultiFab& finemf Real drf = fgeom.CellSize(0); Real rlo = fgeom.Offset(0); if (do_linear_limiting) { - ParallelFor(crsemf, IntVect(-1), + ParallelFor(crsemf, minus1, [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int) noexcept { mf_cell_cons_lin_interp_llslope(i,j,0, tmp[box_no], crse[box_no], ccomp, nc, - cdomain, pbc); + cdomain, ratio, pbc); }); } else { - ParallelFor(crsemf, IntVect(-1), nc, + ParallelFor(crsemf, minus1, nc, [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int, int n) noexcept { mf_cell_cons_lin_interp_mcslope_rz(i,j,n, tmp[box_no], crse[box_no], ccomp, nc, @@ -168,14 +177,14 @@ MFCellConsLinInterp::interp (MultiFab const& crsemf, int ccomp, MultiFab& finemf #endif { if (do_linear_limiting) { - ParallelFor(crsemf, IntVect(-1), + ParallelFor(crsemf, minus1, [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k) noexcept { mf_cell_cons_lin_interp_llslope(i,j,k, tmp[box_no], crse[box_no], ccomp, nc, - cdomain, pbc); + cdomain, ratio, pbc); }); } else { - ParallelFor(crsemf, IntVect(-1), nc, + ParallelFor(crsemf, minus1, nc, [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k, int n) noexcept { mf_cell_cons_lin_interp_mcslope(i,j,k,n, tmp[box_no], crse[box_no], ccomp, nc, @@ -208,7 +217,7 @@ MFCellConsLinInterp::interp (MultiFab const& crsemf, int ccomp, MultiFab& finemf auto const& fine = finemf.array(mfi); auto const& crse = crsemf.const_array(mfi); - Box const& cbox = amrex::grow(crsemf[mfi].box(), -1); + Box const& cbox = amrex::grow(crsemf[mfi].box(), minus1); tmpfab.resize(cbox, AMREX_SPACEDIM*nc); auto const& tmp = tmpfab.array(); auto const& ctmp = tmpfab.const_array(); @@ -224,7 +233,7 @@ MFCellConsLinInterp::interp (MultiFab const& crsemf, int ccomp, MultiFab& finemf [&] (int i, int, int) noexcept { mf_cell_cons_lin_interp_llslope(i,0,0, tmp, crse, ccomp, nc, - cdomain, pbc); + cdomain, ratio, pbc); }); } else { amrex::LoopConcurrentOnCpu(cbox, nc, @@ -251,7 +260,7 @@ MFCellConsLinInterp::interp (MultiFab const& crsemf, int ccomp, MultiFab& finemf [&] (int i, int j, int) noexcept { mf_cell_cons_lin_interp_llslope(i,j,0, tmp, crse, ccomp, nc, - cdomain, pbc); + cdomain, ratio, pbc); }); } else { amrex::LoopConcurrentOnCpu(cbox, nc, @@ -276,7 +285,7 @@ MFCellConsLinInterp::interp (MultiFab const& crsemf, int ccomp, MultiFab& finemf [&] (int i, int j, int k) noexcept { mf_cell_cons_lin_interp_llslope(i,j,k, tmp, crse, ccomp, nc, - cdomain, pbc); + cdomain, ratio, pbc); }); } else { amrex::LoopConcurrentOnCpu(cbox, nc, @@ -304,7 +313,11 @@ Box MFCellConsLinMinmaxLimitInterp::CoarseBox (const Box& fine, const IntVect& ratio) { Box crse = amrex::coarsen(fine,ratio); - crse.grow(1); + for (int dim = 0; dim < AMREX_SPACEDIM; dim++) { + if (ratio[dim] > 1) { + crse.grow(dim,1); + } + } return crse; } @@ -327,6 +340,11 @@ MFCellConsLinMinmaxLimitInterp::interp (MultiFab const& crsemf, int ccomp, Multi Box const& cdomain = cgeom.Domain(); + IntVect minus1; + for (int dim = 0; dim < AMREX_SPACEDIM; dim++) { + minus1[dim] = (ratio[dim] > 1) ? -1 : 0; + } + #ifdef AMREX_USE_GPU if (Gpu::inLaunchRegion()) { MultiFab crse_tmp(crsemf.boxArray(), crsemf.DistributionMap(), AMREX_SPACEDIM*nc, 0); @@ -339,7 +357,7 @@ MFCellConsLinMinmaxLimitInterp::interp (MultiFab const& crsemf, int ccomp, Multi BCRec const* pbc = d_bc.data(); Gpu::copyAsync(Gpu::hostToDevice, bcs.begin()+bcomp, bcs.begin()+bcomp+nc, d_bc.begin()); - ParallelFor(crsemf, IntVect(-1), + ParallelFor(crsemf, minus1, [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k) noexcept { mf_cell_cons_lin_interp_limit_minmax_llslope(i,j,k, tmp[box_no], crse[box_no], ccomp, nc, @@ -370,7 +388,7 @@ MFCellConsLinMinmaxLimitInterp::interp (MultiFab const& crsemf, int ccomp, Multi auto const& fine = finemf.array(mfi); auto const& crse = crsemf.const_array(mfi); - Box const& cbox = amrex::grow(crsemf[mfi].box(), -1); + Box const& cbox = amrex::grow(crsemf[mfi].box(), minus1); tmpfab.resize(cbox, AMREX_SPACEDIM*nc); auto const& tmp = tmpfab.array(); auto const& ctmp = tmpfab.const_array();