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

Jump to bottom

FillPatchSingleLevel and FillPatchTwoLevels for ERF #4158

Merged
merged 6 commits into from
Sep 19, 2024
Merged

FillPatchSingleLevel and FillPatchTwoLevels for ERF #4158

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
5cd8398
FillPatchSingleLevel and FillPatchTwoLevels for ERF
WeiqunZhang Sep 17, 2024
cc713b0
Fix compliation
WeiqunZhang Sep 18, 2024
1d4d15d
Refactor InterpFromCoarseLevel
WeiqunZhang Sep 18, 2024
c5ffc29
bcscomp
WeiqunZhang Sep 18, 2024
6bd4b72
Face data
WeiqunZhang Sep 18, 2024
d3783fb
fix warning
WeiqunZhang Sep 18, 2024
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
73 changes: 72 additions & 1 deletion Src/AmrCore/AMReX_FillPatchUtil.H
View file
Open in desktop
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
Open in desktop
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
Loading