Skip to content

Commit

Permalink
FillPatchSingleLevel and FillPatchTwoLevels for ERF (AMReX-Codes#4158)
Browse files Browse the repository at this point in the history 
These new functions do not take PhysBC functors. So it's the caller's
responsibility to fill domain ghost cells before calling.
  • Loading branch information
@WeiqunZhang
WeiqunZhang authored Sep 19, 2024
1 parent 6c7668c commit 3734079
Show file tree
Hide file tree
Showing 3 changed files with 229 additions and 55 deletions.
73 changes: 72 additions & 1 deletion Src/AmrCore/AMReX_FillPatchUtil.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -718,6 +718,7 @@ namespace amrex
* \param ratio refinement ratio
* \param mapper spatial interpolater
* \param bcs boundar types for each component
* \param bcscomp starting component for bcs
*/
template <typename MF, typename Interp>
std::enable_if_t<IsFabArray<MF>::value>
Expand All @@ -726,7 +727,77 @@ namespace amrex
const MF& cmf, int scomp, int dcomp, int ncomp,
const Geometry& cgeom, const Geometry& fgeom,
const IntVect& ratio, Interp* mapper,
const Vector<BCRec>& bcs);
const Vector<BCRec>& bcs, int bcscomp);

/**
* \brief FillPatch with data from the current level
*
* In this version of FillPatchSingleLevel, it's the CALLER's
* responsibility to make sure that smf has `snghost` ghost cells
* already filled before calling this function. The destination
* MultiFab/FabArray is on the same AMR level as the source
* MultiFab/FabArray. If needed, interpolation in time is performed.
*
* \tparam MF the MultiFab/FabArray type
*
* \param mf destination MF
* \param nghost number of ghost cells of mf needed to be filled
* \param time time associated with mf
* \param smf source MFs
* \param snghost number of ghost cells in smf with valid data
* \param stime times associated smf
* \param scomp starting component of the source MFs
* \param dcomp starting component of the destination MF
* \param ncomp number of components
* \param geom Geometry for this level
*/
template <typename MF>
std::enable_if_t<IsFabArray<MF>::value>
FillPatchSingleLevel (MF& mf, IntVect const& nghost, Real time,
const Vector<MF*>& smf, IntVect const& snghost,
const Vector<Real>& stime, int scomp, int dcomp, int ncomp,
const Geometry& geom);

/**
* \brief FillPatch with data from the current level and the level below.
*
* In this version of FillPatchTwoLevels, it's the CALLER's
* responsibility to make sure all ghost cells of the coarse MF needed
* for interpolation are filled already before calling this
* function. It's assumed that the fine level MultiFab mf's BoxArray is
* coarsenable by the refinement ratio. There is no support for EB.
*
* \tparam MF the MultiFab/FabArray type
* \tparam Interp spatial interpolater
*
* \param mf destination MF on the fine level
* \param nghost number of ghost cells of mf inside domain needed to be filled
* \param nghost_outside_domain number of ghost cells of mf outside domain needed to be filled
* \param time time associated with mf
* \param cmf source MFs on the coarse level
* \param ct times associated cmf
* \param fmf source MFs on the fine level
* \param ft times associated fmf
* \param scomp starting component of the source MFs
* \param dcomp starting component of the destination MF
* \param ncomp number of components
* \param cgeom Geometry for the coarse level
* \param fgeom Geometry for the fine level
* \param ratio refinement ratio
* \param mapper spatial interpolater
* \param bcs boundary types for each component.
* \param bcscomp starting component for bcs
*/
template <typename MF, typename Interp>
std::enable_if_t<IsFabArray<MF>::value>
FillPatchTwoLevels (MF& mf, IntVect const& nghost,
IntVect const& nghost_outside_domain, Real time,
const Vector<MF*>& cmf, const Vector<Real>& ct,
const Vector<MF*>& fmf, const Vector<Real>& ft,
int scomp, int dcomp, int ncomp,
const Geometry& cgeom, const Geometry& fgeom,
const IntVect& ratio, Interp* mapper,
const Vector<BCRec>& bcs, int bcscomp);

#ifndef BL_NO_FORT
enum InterpEM_t { InterpE, InterpB};
Expand Down
161 changes: 107 additions & 54 deletions Src/AmrCore/AMReX_FillPatchUtil_I.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -84,12 +84,17 @@ FillPatchSingleLevel (MF& mf, IntVect const& nghost, Real time,
AMREX_ASSERT(!smf.empty());
AMREX_ASSERT(nghost.allLE(mf.nGrowVect()));

IntVect src_ghost(0);
if constexpr (std::is_same_v<BC,PhysBCFunctUseCoarseGhost>) {
src_ghost = physbcf.fp1_src_ghost;
}

if (smf.size() == 1)
{
if (&mf == smf[0] && scomp == dcomp) {
mf.FillBoundary(dcomp, ncomp, nghost, geom.periodicity());
} else {
mf.ParallelCopy(*smf[0], scomp, dcomp, ncomp, IntVect{0}, nghost, geom.periodicity());
mf.ParallelCopy(*smf[0], scomp, dcomp, ncomp, src_ghost, nghost, geom.periodicity());
}
}
else if (smf.size() == 2)
Expand All @@ -106,7 +111,7 @@ FillPatchSingleLevel (MF& mf, IntVect const& nghost, Real time,
destcomp = dcomp;
sameba = true;
} else {
raii.define(smf[0]->boxArray(), smf[0]->DistributionMap(), ncomp, 0,
raii.define(smf[0]->boxArray(), smf[0]->DistributionMap(), ncomp, src_ghost,
MFInfo(), smf[0]->Factory());

dmf = &raii;
Expand All @@ -116,12 +121,19 @@ FillPatchSingleLevel (MF& mf, IntVect const& nghost, Real time,

if ((dmf != smf[0] && dmf != smf[1]) || scomp != dcomp)
{
IntVect interp_ghost(0);
if constexpr (std::is_same_v<BC,PhysBCFunctUseCoarseGhost>) {
interp_ghost = physbcf.fp1_src_ghost;
if (sameba) {
interp_ghost.min(nghost);
}
}
#ifdef AMREX_USE_OMP
#pragma omp parallel if (Gpu::notInLaunchRegion())
#endif
for (MFIter mfi(*dmf,TilingIfNotGPU()); mfi.isValid(); ++mfi)
{
const Box& bx = mfi.tilebox();
const Box& bx = mfi.growntilebox(interp_ghost);
const Real t0 = stime[0];
const Real t1 = stime[1];
auto const sfab0 = smf[0]->array(mfi);
Expand Down Expand Up @@ -170,10 +182,7 @@ FillPatchSingleLevel (MF& mf, IntVect const& nghost, Real time,
}
else
{
IntVect src_ngrow = IntVect::TheZeroVector();
IntVect dst_ngrow = nghost;

mf.ParallelCopy(*dmf, 0, dcomp, ncomp, src_ngrow, dst_ngrow, geom.periodicity());
mf.ParallelCopy(*dmf, 0, dcomp, ncomp, src_ghost, nghost, geom.periodicity());
}
}
else {
Expand Down Expand Up @@ -504,10 +513,16 @@ namespace detail {
auto solve_mask = make_mf_crse_mask<iMultiFab>(fpc, ncomp, mf.boxArray().ixType(), ratio);

mf_set_domain_bndry(mf_crse_patch, cgeom);
if constexpr (std::is_same_v<BC,PhysBCFunctUseCoarseGhost>) {
cbc.fp1_src_ghost = cbc.cghost;
}
FillPatchSingleLevel(mf_crse_patch, time, cmf, ct, scomp, 0, ncomp,
cgeom, cbc, cbccomp);

mf_set_domain_bndry(mf_refined_patch, fgeom);
if constexpr (std::is_same_v<BC,PhysBCFunctUseCoarseGhost>) {
fbc.fp1_src_ghost = IntVect(0);
}
FillPatchSingleLevel(mf_refined_patch, time, fmf, ft, scomp, 0, ncomp,
fgeom, fbc, fbccomp);

Expand Down Expand Up @@ -565,16 +580,29 @@ namespace detail {
MF mf_crse_patch = make_mf_crse_patch<MF>(fpc, ncomp);
mf_set_domain_bndry (mf_crse_patch, cgeom);

if constexpr (std::is_same_v<BC,PhysBCFunctUseCoarseGhost>) {
cbc.fp1_src_ghost = cbc.cghost;
}
FillPatchSingleLevel(mf_crse_patch, time, cmf, ct, scomp, 0, ncomp, cgeom, cbc, cbccomp);

MF mf_fine_patch = make_mf_fine_patch<MF>(fpc, ncomp);

detail::call_interp_hook(pre_interp, mf_crse_patch, 0, ncomp);

Box fdomain_g( amrex::convert(fgeom.Domain(),mf.ixType()) );
for (int i = 0; i < AMREX_SPACEDIM; ++i) {
if (fgeom.isPeriodic(i)) {
fdomain_g.grow(i, nghost[i]);
} else {
if constexpr (std::is_same_v
<BC, PhysBCFunctUseCoarseGhost>) {
fdomain_g.grow(i, fbc.nghost_outside_domain[i]);
}
}
}
FillPatchInterp(mf_fine_patch, 0, mf_crse_patch, 0,
ncomp, IntVect(0), cgeom, fgeom,
amrex::grow(amrex::convert(fgeom.Domain(),mf.ixType()),nghost),
ratio, mapper, bcs, bcscomp);
fdomain_g, ratio, mapper, bcs, bcscomp);

detail::call_interp_hook(post_interp, mf_fine_patch, 0, ncomp);

Expand All @@ -583,6 +611,9 @@ namespace detail {
}
}

if constexpr(std::is_same_v<BC, PhysBCFunctUseCoarseGhost>) {
fbc.fp1_src_ghost = IntVect(0);
}
FillPatchSingleLevel(mf, nghost, time, fmf, ft, scomp, dcomp, ncomp,
fgeom, fbc, fbccomp);

Expand Down Expand Up @@ -997,6 +1028,8 @@ InterpFromCoarseLevel (MF& mf, IntVect const& nghost, Real time,
const PreInterpHook& pre_interp,
const PostInterpHook& post_interp)
{
BL_PROFILE("InterpFromCoarseLevel");

using FAB = typename MF::FABType::value_type;

const InterpolaterBoxCoarsener& coarsener = mapper->BoxCoarsener(ratio);
Expand All @@ -1011,33 +1044,45 @@ InterpFromCoarseLevel (MF& mf, IntVect const& nghost, Real time,
Box fdomain_g( amrex::convert(fgeom.Domain(),mf.ixType()) );
for (int i = 0; i < AMREX_SPACEDIM; ++i) {
if (fgeom.isPeriodic(i)) {
fdomain_g.grow(i,nghost[i]);
fdomain_g.grow(i, nghost[i]);
} else {
if constexpr (std::is_same_v<BC, PhysBCFunctUseCoarseGhost>) {
fdomain_g.grow(i, fbc.nghost_outside_domain[i]);
}
}
}

BoxArray ba_crse_patch(ba.size());
{ // TODO: later we might want to cache this
for (int i = 0, N = ba.size(); i < N; ++i)
{
Box bx = amrex::convert(amrex::grow(ba[i],nghost), typ);
bx &= fdomain_g;
ba_crse_patch.set(i, coarsener.doit(bx));
MF mf_crse_patch;
IntVect send_ghost(0), recv_ghost(0);
if constexpr (std::is_same_v<BC, PhysBCFunctUseCoarseGhost>) {
mf_crse_patch.define(amrex::coarsen(ba,ratio), dm, ncomp, fbc.src_ghost);
send_ghost = fbc.cghost;
recv_ghost = fbc.src_ghost;
} else {
BoxArray ba_crse_patch(ba.size());
{ // TODO: later we might want to cache this
for (int i = 0, N = ba.size(); i < N; ++i)
{
Box bx = amrex::convert(amrex::grow(ba[i],nghost), typ);
bx &= fdomain_g;
ba_crse_patch.set(i, coarsener.doit(bx));
}
}
}

MF mf_crse_patch;
#ifdef AMREX_USE_EB
if (EB2::TopIndexSpaceIfPresent()) {
auto factory = makeEBFabFactory(cgeom, ba_crse_patch, dm, {0,0,0}, EBSupport::basic);
mf_crse_patch.define(ba_crse_patch, dm, ncomp, 0, MFInfo(), *factory);
} else
if (EB2::TopIndexSpaceIfPresent()) {
auto factory = makeEBFabFactory(cgeom, ba_crse_patch, dm, {0,0,0}, EBSupport::basic);
mf_crse_patch.define(ba_crse_patch, dm, ncomp, 0, MFInfo(), *factory);
} else
#endif
{
mf_crse_patch.define(ba_crse_patch, dm, ncomp, 0);
{
mf_crse_patch.define(ba_crse_patch, dm, ncomp, 0);
}
detail::mf_set_domain_bndry (mf_crse_patch, cgeom);
}
detail::mf_set_domain_bndry (mf_crse_patch, cgeom);

mf_crse_patch.ParallelCopy(cmf, scomp, 0, ncomp, cgeom.periodicity());
mf_crse_patch.ParallelCopy(cmf, scomp, 0, ncomp, send_ghost, recv_ghost,
cgeom.periodicity());

cbc(mf_crse_patch, 0, ncomp, mf_crse_patch.nGrowVect(), time, cbccomp);

Expand Down Expand Up @@ -1075,6 +1120,8 @@ InterpFromCoarseLevel (Array<MF*, AMREX_SPACEDIM> const& mf, IntVect const& ngho
const PreInterpHook& pre_interp,
const PostInterpHook& post_interp)
{
BL_PROFILE("InterpFromCoarseLevel(array)");

using FAB = typename MF::FABType::value_type;
using iFAB = typename iMultiFab::FABType::value_type;

Expand Down Expand Up @@ -1217,36 +1264,42 @@ InterpFromCoarseLevel (MF& mf, IntVect const& nghost,
const MF& cmf, int scomp, int dcomp, int ncomp,
const Geometry& cgeom, const Geometry& fgeom,
const IntVect& ratio, Interp* mapper,
const Vector<BCRec>& bcs)
const Vector<BCRec>& bcs, int bcscomp)
{
const BoxArray& ba = mf.boxArray();
const DistributionMapping& dm = mf.DistributionMap();

const IndexType& typ = ba.ixType();
BL_ASSERT(typ == cmf.boxArray().ixType());

const InterpolaterBoxCoarsener& coarsener = mapper->BoxCoarsener(ratio);
Box tmp(-nghost, IntVect(32), typ);
Box tmp2 = coarsener.doit(tmp);
IntVect src_ghost = -tmp2.smallEnd();

tmp = Box(-nghost_outside_domain, IntVect(32), typ);
tmp2 = coarsener.doit(tmp);
IntVect src_ghost_outside_domain = -tmp2.smallEnd();

IntVect cghost = cmf.nGrowVect();
cghost.min(src_ghost);

AMREX_ALWAYS_ASSERT(cghost.allGE(src_ghost_outside_domain));
PhysBCFunctUseCoarseGhost erfbc(cmf,nghost,nghost_outside_domain,ratio,mapper);
InterpFromCoarseLevel(mf, nghost, Real(0.0), cmf, scomp, dcomp, ncomp,
cgeom, fgeom, erfbc, 0, erfbc, 0, ratio, mapper,
bcs, bcscomp);
}

MF mf_crse_patch(amrex::coarsen(ba,ratio), dm, ncomp, src_ghost);
mf_crse_patch.ParallelCopy(cmf, scomp, 0, ncomp, cghost, src_ghost,
cgeom.periodicity());
template <typename MF>
std::enable_if_t<IsFabArray<MF>::value>
FillPatchSingleLevel (MF& mf, IntVect const& nghost, Real time,
const Vector<MF*>& smf, IntVect const& snghost,
const Vector<Real>& stime, int scomp, int dcomp, int ncomp,
const Geometry& geom)
{
PhysBCFunctUseCoarseGhost erfbc(snghost);
FillPatchSingleLevel(mf, nghost, time, smf, stime, scomp, dcomp, ncomp, geom,
erfbc, 0);
}

Box fdomain_g = amrex::convert(fgeom.Domain(),typ);
fdomain_g.grow(nghost_outside_domain);
FillPatchInterp(mf, dcomp, mf_crse_patch, 0, ncomp, nghost, cgeom, fgeom,
fdomain_g, ratio, mapper, bcs, 0);
template <typename MF, typename Interp>
std::enable_if_t<IsFabArray<MF>::value>
FillPatchTwoLevels (MF& mf, IntVect const& nghost,
IntVect const& nghost_outside_domain, Real time,
const Vector<MF*>& cmf, const Vector<Real>& ct,
const Vector<MF*>& fmf, const Vector<Real>& ft,
int scomp, int dcomp, int ncomp,
const Geometry& cgeom, const Geometry& fgeom,
const IntVect& ratio, Interp* mapper,
const Vector<BCRec>& bcs, int bcscomp)
{
PhysBCFunctUseCoarseGhost erfbc(*cmf[0], nghost, nghost_outside_domain, ratio,
mapper);
FillPatchTwoLevels(mf, nghost, time, cmf, ct, fmf, ft, scomp, dcomp, ncomp,
cgeom, fgeom, erfbc, 0, erfbc, 0, ratio, mapper,
bcs, bcscomp);
}

template <typename MF, typename BC, typename Interp>
Expand Down
Loading

0 comments on commit 3734079

Please sign in to comment.