Skip to content

Commit

Permalink
Fourth-order interpolation from fine to coarse level (#2987)
Browse files Browse the repository at this point in the history 
For fourth-order finite-difference methods with data at cell centers, we
cannot use the usual averageDown function to overwrite coarse level data
with fine data. We actually need to do interpolation.
  • Loading branch information
@WeiqunZhang
WeiqunZhang authored Oct 14, 2022
1 parent 975b830 commit c841ae8
Show file tree
Hide file tree
Showing 2 changed files with 218 additions and 114 deletions.
17 changes: 17 additions & 0 deletions Src/Base/AMReX_MultiFabUtil.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -243,6 +243,23 @@ namespace amrex
Vector<Geometry> const& geom,
Vector<IntVect> const& ratio,
bool local = false);

/**
* \brief Fourth-order interpolation from fine to coarse level.
*
* This is for high-order "average-down" of finite-difference data. If
* ghost cell data are used, it's the caller's responsibility to fill
* the ghost cells before calling this function.
*
* \param cmf coarse data
* \param scomp starting component
* \param ncomp number of component
* \param fmf fine data
* \param ratio refinement ratio.
*/
void FourthOrderInterpFromFineToCoarse (MultiFab& cmf, int scomp, int ncomp,
MultiFab const& fmf,
IntVect const& ratio);
}

namespace amrex {
Expand Down
315 changes: 201 additions & 114 deletions Src/Base/AMReX_MultiFabUtil.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1227,157 +1227,244 @@ namespace amrex
return hv;
}

Real volumeWeightedSum (Vector<MultiFab const*> const& mf, int icomp,
Vector<Geometry> const& geom,
Vector<IntVect> const& ratio,
bool local)
{
ReduceOps<ReduceOpSum> reduce_op;
ReduceData<Real> reduce_data(reduce_op);
Real volumeWeightedSum (Vector<MultiFab const*> const& mf, int icomp,
Vector<Geometry> const& geom,
Vector<IntVect> const& ratio,
bool local)
{
ReduceOps<ReduceOpSum> reduce_op;
ReduceData<Real> reduce_data(reduce_op);

#ifdef AMREX_USE_EB
bool has_eb = !(mf[0]->isAllRegular());
bool has_eb = !(mf[0]->isAllRegular());
#endif

int nlevels = mf.size();
for (int ilev = 0; ilev < nlevels-1; ++ilev) {
iMultiFab mask = makeFineMask(*mf[ilev], *mf[ilev+1], IntVect(0),
ratio[ilev],Periodicity::NonPeriodic(),
0, 1);
auto const& m = mask.const_arrays();
auto const& a = mf[ilev]->const_arrays();
auto const dx = geom[ilev].CellSizeArray();
Real dv = AMREX_D_TERM(dx[0],*dx[1],*dx[2]);
#ifdef AMREX_USE_EB
if (has_eb) {
AMREX_ASSERT(mf[ilev]->hasEBFabFactory());
auto const& f = dynamic_cast<EBFArrayBoxFactory const&>
(mf[ilev]->Factory());
auto const& vfrac = f.getVolFrac();
auto const& va = vfrac.const_arrays();
reduce_op.eval(*mf[ilev], IntVect(0), reduce_data,
[=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k) noexcept
-> Real
{
return m[box_no](i,j,k) ? Real(0.)
: dv*a[box_no](i,j,k,icomp)*va[box_no](i,j,k);
});
} else
#endif
{
#if (AMREX_SPACEDIM == 1)
if (geom[ilev].IsSPHERICAL()) {
const auto rlo = geom[ilev].CellSize(0);
reduce_op.eval(*mf[ilev], IntVect(0), reduce_data,
[=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k)
noexcept -> Real
{
if (m[box_no](i,j,k)) {
return Real(0.);
} else {
constexpr Real pi = Real(3.1415926535897932);
Real ri = rlo + dx[0]*i;
Real ro = ri + dx[0];
return Real(4./3.)*pi*(ro-ri)*(ro*ro+ro*ri+ri*ri)
* a[box_no](i,j,k,icomp);
}
});
} else
#elif (AMREX_SPACEDIM == 2)
if (geom[ilev].IsRZ()) {
const auto rlo = geom[ilev].CellSize(0);
reduce_op.eval(*mf[ilev], IntVect(0), reduce_data,
[=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k)
noexcept -> Real
{
if (m[box_no](i,j,k)) {
return Real(0.);
} else {
Real ri = rlo + dx[0]*i;
Real ro = ri + dx[0];
constexpr Real pi = Real(3.1415926535897932);
return pi*dx[1]*dx[0]*(ro+ri)
* a[box_no](i,j,k,icomp);
}
});
} else
#endif
{
reduce_op.eval(*mf[ilev], IntVect(0), reduce_data,
[=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k)
noexcept -> Real
{
return m[box_no](i,j,k) ? Real(0.)
: dv*a[box_no](i,j,k,icomp);
});
}
}
Gpu::streamSynchronize();
}

auto const& a = mf.back()->const_arrays();
auto const dx = geom[nlevels-1].CellSizeArray();
int nlevels = mf.size();
for (int ilev = 0; ilev < nlevels-1; ++ilev) {
iMultiFab mask = makeFineMask(*mf[ilev], *mf[ilev+1], IntVect(0),
ratio[ilev],Periodicity::NonPeriodic(),
0, 1);
auto const& m = mask.const_arrays();
auto const& a = mf[ilev]->const_arrays();
auto const dx = geom[ilev].CellSizeArray();
Real dv = AMREX_D_TERM(dx[0],*dx[1],*dx[2]);
#ifdef AMREX_USE_EB
if (has_eb) {
AMREX_ASSERT(mf.back()->hasEBFabFactory());
AMREX_ASSERT(mf[ilev]->hasEBFabFactory());
auto const& f = dynamic_cast<EBFArrayBoxFactory const&>
(mf.back()->Factory());
(mf[ilev]->Factory());
auto const& vfrac = f.getVolFrac();
auto const& va = vfrac.const_arrays();
reduce_op.eval(*mf.back(), IntVect(0), reduce_data,
reduce_op.eval(*mf[ilev], IntVect(0), reduce_data,
[=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k) noexcept
-> Real
{
return dv*a[box_no](i,j,k,icomp)*va[box_no](i,j,k);
return m[box_no](i,j,k) ? Real(0.)
: dv*a[box_no](i,j,k,icomp)*va[box_no](i,j,k);
});
} else
#endif
{
#if (AMREX_SPACEDIM == 1)
if (geom[nlevels-1].IsSPHERICAL()) {
const auto rlo = geom[nlevels-1].CellSize(0);
reduce_op.eval(*mf.back(), IntVect(0), reduce_data,
if (geom[ilev].IsSPHERICAL()) {
const auto rlo = geom[ilev].CellSize(0);
reduce_op.eval(*mf[ilev], IntVect(0), reduce_data,
[=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k)
noexcept -> Real
{
constexpr Real pi = Real(3.1415926535897932);
Real ri = rlo + dx[0]*i;
Real ro = ri + dx[0];
return Real(4./3.)*pi*(ro-ri)*(ro*ro+ro*ri+ri*ri)
* a[box_no](i,j,k,icomp);
if (m[box_no](i,j,k)) {
return Real(0.);
} else {
constexpr Real pi = Real(3.1415926535897932);
Real ri = rlo + dx[0]*i;
Real ro = ri + dx[0];
return Real(4./3.)*pi*(ro-ri)*(ro*ro+ro*ri+ri*ri)
* a[box_no](i,j,k,icomp);
}
});
} else
#elif (AMREX_SPACEDIM == 2)
if (geom[nlevels-1].IsRZ()) {
const auto rlo = geom[nlevels-1].CellSize(0);
reduce_op.eval(*mf.back(), IntVect(0), reduce_data,
if (geom[ilev].IsRZ()) {
const auto rlo = geom[ilev].CellSize(0);
reduce_op.eval(*mf[ilev], IntVect(0), reduce_data,
[=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k)
noexcept -> Real
{
Real ri = rlo + dx[0]*i;
Real ro = ri + dx[0];
constexpr Real pi = Real(3.1415926535897932);
return pi*dx[1]*dx[0]*(ro+ri)
* a[box_no](i,j,k,icomp);
if (m[box_no](i,j,k)) {
return Real(0.);
} else {
Real ri = rlo + dx[0]*i;
Real ro = ri + dx[0];
constexpr Real pi = Real(3.1415926535897932);
return pi*dx[1]*dx[0]*(ro+ri)
* a[box_no](i,j,k,icomp);
}
});
} else
#endif
{
reduce_op.eval(*mf.back(), IntVect(0), reduce_data,
[=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k) noexcept
reduce_op.eval(*mf[ilev], IntVect(0), reduce_data,
[=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k)
noexcept -> Real
{
return dv*a[box_no](i,j,k,icomp);
return m[box_no](i,j,k) ? Real(0.)
: dv*a[box_no](i,j,k,icomp);
});
}
}
Gpu::streamSynchronize();
}

auto const& a = mf.back()->const_arrays();
auto const dx = geom[nlevels-1].CellSizeArray();
Real dv = AMREX_D_TERM(dx[0],*dx[1],*dx[2]);
#ifdef AMREX_USE_EB
if (has_eb) {
AMREX_ASSERT(mf.back()->hasEBFabFactory());
auto const& f = dynamic_cast<EBFArrayBoxFactory const&>
(mf.back()->Factory());
auto const& vfrac = f.getVolFrac();
auto const& va = vfrac.const_arrays();
reduce_op.eval(*mf.back(), IntVect(0), reduce_data,
[=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k) noexcept
-> Real
{
return dv*a[box_no](i,j,k,icomp)*va[box_no](i,j,k);
});
} else
#endif
{
#if (AMREX_SPACEDIM == 1)
if (geom[nlevels-1].IsSPHERICAL()) {
const auto rlo = geom[nlevels-1].CellSize(0);
reduce_op.eval(*mf.back(), IntVect(0), reduce_data,
[=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k)
noexcept -> Real
{
constexpr Real pi = Real(3.1415926535897932);
Real ri = rlo + dx[0]*i;
Real ro = ri + dx[0];
return Real(4./3.)*pi*(ro-ri)*(ro*ro+ro*ri+ri*ri)
* a[box_no](i,j,k,icomp);
});
} else
#elif (AMREX_SPACEDIM == 2)
if (geom[nlevels-1].IsRZ()) {
const auto rlo = geom[nlevels-1].CellSize(0);
reduce_op.eval(*mf.back(), IntVect(0), reduce_data,
[=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k)
noexcept -> Real
{
Real ri = rlo + dx[0]*i;
Real ro = ri + dx[0];
constexpr Real pi = Real(3.1415926535897932);
return pi*dx[1]*dx[0]*(ro+ri)
* a[box_no](i,j,k,icomp);
});
} else
#endif
{
reduce_op.eval(*mf.back(), IntVect(0), reduce_data,
[=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k) noexcept
{
return dv*a[box_no](i,j,k,icomp);
});
}
}

auto const& hv = reduce_data.value(reduce_op);
Real r = amrex::get<0>(hv);

if (!local) {
ParallelAllReduce::Sum(r, ParallelContext::CommunicatorSub());
}
return r;
}

void FourthOrderInterpFromFineToCoarse (MultiFab& cmf, int scomp, int ncomp,
MultiFab const& fmf,
IntVect const& ratio)
{
AMREX_ASSERT(AMREX_D_TERM( (ratio[0] == 2 || ratio[0] == 4),
&& (ratio[1] == 2 || ratio[1] == 4),
&& (ratio[2] == 2 || ratio[2] == 4)));

MultiFab tmp(amrex::coarsen(fmf.boxArray(), ratio), fmf.DistributionMap(),
ncomp, 0);

#ifdef AMREX_USE_OMP
#pragma omp parallel if (Gpu::notInLaunchRegion())
#endif
{
#if (AMREX_SPACEDIM > 1)
FArrayBox xtmp;
#if (AMREX_SPACEDIM > 2)
FArrayBox ytmp;
#endif
#endif
for (MFIter mfi(tmp,TilingIfNotGPU()); mfi.isValid(); ++mfi) {
Box const& bx = mfi.tilebox();
auto const& fa = fmf.const_array(mfi,scomp);

auto const& hv = reduce_data.value(reduce_op);
Real r = amrex::get<0>(hv);
Box xbx = bx;
#if (AMREX_SPACEDIM == 1)
auto const& xa = tmp.array(mfi);
#else
xbx.refine(IntVect(AMREX_D_DECL(1,ratio[1],ratio[2])));
if (ratio[1] == 2) { xbx.grow(1,1); }
#if (AMREX_SPACEDIM == 3)
if (ratio[2] == 2) { xbx.grow(2,1); }
#endif
xtmp.resize(xbx,ncomp);
Elixir eli = xtmp.elixir();
auto const& xa = xtmp.array();
#endif
AMREX_HOST_DEVICE_PARALLEL_FOR_4D(xbx, ncomp, i, j, k, n,
{
int ii = 2*i;
xa(i,j,k,n) = Real(1./16)*(Real(9.)*(fa(ii ,j,k,n) +
fa(ii+1,j,k,n))
- fa(ii-1,j,k,n)
- fa(ii+2,j,k,n));
});

if (!local) {
ParallelAllReduce::Sum(r, ParallelContext::CommunicatorSub());
#if (AMREX_SPACEDIM > 1)
Box ybx = bx;
auto const& xca = xtmp.const_array();
#if (AMREX_SPACEDIM == 2)
auto const& ya = tmp.array(mfi);
#else
ybx.refine(IntVect(AMREX_D_DECL(1,1,ratio[2])));
if (ratio[2] == 2) { ybx.grow(2,1); }
ytmp.resize(ybx,ncomp);
eli.append(ytmp.elixir());
auto const& ya = ytmp.array();
#endif
AMREX_HOST_DEVICE_PARALLEL_FOR_4D(ybx, ncomp, i, j, k, n,
{
int jj = 2*j;
ya(i,j,k,n) = Real(1./16)*(Real(9.)*(xca(i,jj ,k,n) +
xca(i,jj+1,k,n))
- xca(i,jj-1,k,n)
- xca(i,jj+2,k,n));
});

#if (AMREX_SPACEDIM == 3)
auto const& yca = ytmp.const_array();
auto const& ca = tmp.array(mfi);
AMREX_HOST_DEVICE_PARALLEL_FOR_4D(bx, ncomp, i, j, k, n,
{
int kk = 2*k;
ca(i,j,k,n) = Real(1./16)*(Real(9.)*(yca(i,j,kk ,n) +
yca(i,j,kk+1,n))
- yca(i,j,kk-1,n)
- yca(i,j,kk+2,n));
});
#endif
#endif
}
return r;
}

cmf.ParallelCopy(tmp, 0, scomp, ncomp);
}
}

0 comments on commit c841ae8

Please sign in to comment.