Skip to content

Commit

Permalink
generalize interpolaters to not need or use extra values in a directi… (
Browse files Browse the repository at this point in the history 
#3415)

…on for which the refinement ratio is 1

The previous version of the interpolation routines in AmrCore would grow
the "CoarseBox" in all directions even if one of the directions had
refinement ratio 1. This PR changes that behavior so the CoarseBox is
not grown, and slopes are not computed, in any direction for which the
refinement ratio is 1.
  • Loading branch information
@asalmgren
asalmgren committed Jul 16, 2023
1 parent 5a90c26 commit 01b750d
Show file tree
Hide file tree
Showing 5 changed files with 274 additions and 158 deletions.
48 changes: 27 additions & 21 deletions Src/AmrCore/AMReX_Interpolater.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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.
*
Expand Down Expand Up @@ -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;
}

Expand Down Expand Up @@ -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));

Expand All @@ -454,7 +459,12 @@ CellConservativeLinear::interp (const FArrayBox& crse,
Array4<Real> 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<Real> const& tmp = ccfab.array();
Expand All @@ -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,
Expand All @@ -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,
Expand All @@ -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,
Expand Down Expand Up @@ -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.
Expand All @@ -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());
Expand Down Expand Up @@ -739,6 +749,8 @@ CellConservativeProtected::protect (const FArrayBox& /*crse*/,
Vector<BCRec>& /*bcr*/,
RunOn runon)
{
AMREX_ALWAYS_ASSERT(ratio.allGT(IntVect(1)));

#if (AMREX_SPACEDIM == 1)
amrex::ignore_unused(fine,fine_state,
ncomp,fine_region,ratio,
Expand Down Expand Up @@ -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);

//
Expand All @@ -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
Expand All @@ -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
Expand Down
2 changes: 1 addition & 1 deletion Src/AmrCore/AMReX_MFInterp_1D_C.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -38,7 +38,7 @@ void mf_cell_cons_lin_interp_limit_minmax_llslope (int i, int, int, Array4<Real>
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
void mf_cell_cons_lin_interp_llslope (int i, int, int, Array4<Real> const& slope,
Array4<Real const> 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);

Expand Down
152 changes: 97 additions & 55 deletions Src/AmrCore/AMReX_MFInterp_2D_C.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -13,24 +13,38 @@ void mf_cell_cons_lin_interp_limit_minmax_llslope (int i, int j, int, Array4<Rea
Real sfx = Real(1.0);
Real sfy = Real(1.0);

Real sx = Real(0.0);
Real sy = Real(0.0);

Real dcx = Real(0.0);
Real dcy = Real(0.0);

for (int ns = 0; ns < ncomp; ++ns) {
int nu = ns + scomp;

// x-direction
Real 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));
Real 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;
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);
Expand All @@ -39,8 +53,10 @@ void mf_cell_cons_lin_interp_limit_minmax_llslope (int i, int j, int, Array4<Rea
+ 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));
}}
Expand Down Expand Up @@ -73,7 +89,7 @@ void mf_cell_cons_lin_interp_limit_minmax_llslope (int i, int j, int, Array4<Rea
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
void mf_cell_cons_lin_interp_llslope (int i, int j, int, Array4<Real> const& slope,
Array4<Real const> 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);
Expand All @@ -82,26 +98,34 @@ void mf_cell_cons_lin_interp_llslope (int i, int j, int, Array4<Real> 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) {
Expand All @@ -118,28 +142,37 @@ void mf_cell_cons_lin_interp_mcslope (int i, int j, int /*k*/, int ns, Array4<Re
{
int nu = ns + scomp;

Real sx = Real(0.);
Real sy = Real(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)) {
Real dumax = std::abs(sx) * Real(ratio[0]-1)/Real(2*ratio[0])
+ 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));
}}
Expand Down Expand Up @@ -177,19 +210,26 @@ void mf_cell_cons_lin_interp_mcslope_rz (int i, int j, int ns, Array4<Real> 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)) {
Expand All @@ -214,8 +254,10 @@ void mf_cell_cons_lin_interp_mcslope_rz (int i, int j, int ns, Array4<Real> 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));
}}
Expand Down
Loading

0 comments on commit 01b750d

Please sign in to comment.