Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Fourth-order interpolation from fine to coarse level #2987

Merged
merged 1 commit into from
Oct 14, 2022
Merged

Fourth-order interpolation from fine to coarse level #2987

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
17 changes: 17 additions & 0 deletions Src/Base/AMReX_MultiFabUtil.H
View file
Open in desktop
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
Open in desktop
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();
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Does it make sense to use nlevels-1 here (and the following reduce_op.eval calls)? This is primarily a readability question (assuming I understood the length of mf)

Suggested change
auto const& a = mf.back()->const_arrays();
auto const& a = mf[nlevels-1]->const_arrays();

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Sorry. Too Late. I enabled auto-merge. Nevertheless, one could argue that mf.back() is better because one is less likely to make a mistake such as mf[nlevels].

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);
});
}
}
Comment on lines +1321 to +1373
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
#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);
});
}
}
#ifdef AMREX_USE_EB
if (has_eb) {
AMREX_ASSERT(mf[nlevels-1]->hasEBFabFactory());
auto const& f = dynamic_cast<EBFArrayBoxFactory const&>
(mf[nlevels-1]->Factory());
auto const& vfrac = f.getVolFrac();
auto const& va = vfrac.const_arrays();
reduce_op.eval(*mf[nlevels-1], 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[nlevels-1], 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[nlevels-1], 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[nlevels-1], 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);
}
}