From 3734079379bb6b2a3850d197241f6b2c3b3bfa7d Mon Sep 17 00:00:00 2001 From: Weiqun Zhang Date: Wed, 18 Sep 2024 22:49:24 -0500 Subject: [PATCH] FillPatchSingleLevel and FillPatchTwoLevels for ERF (#4158) These new functions do not take PhysBC functors. So it's the caller's responsibility to fill domain ghost cells before calling. --- Src/AmrCore/AMReX_FillPatchUtil.H | 73 ++++++++++++- Src/AmrCore/AMReX_FillPatchUtil_I.H | 161 ++++++++++++++++++---------- Src/Base/AMReX_PhysBCFunct.H | 50 +++++++++ 3 files changed, 229 insertions(+), 55 deletions(-) diff --git a/Src/AmrCore/AMReX_FillPatchUtil.H b/Src/AmrCore/AMReX_FillPatchUtil.H index 11a8eb560ec..b32c3e02947 100644 --- a/Src/AmrCore/AMReX_FillPatchUtil.H +++ b/Src/AmrCore/AMReX_FillPatchUtil.H @@ -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 std::enable_if_t::value> @@ -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& bcs); + const Vector& 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 + std::enable_if_t::value> + FillPatchSingleLevel (MF& mf, IntVect const& nghost, Real time, + const Vector& smf, IntVect const& snghost, + const Vector& 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 + std::enable_if_t::value> + FillPatchTwoLevels (MF& mf, IntVect const& nghost, + IntVect const& nghost_outside_domain, Real time, + const Vector& cmf, const Vector& ct, + const Vector& fmf, const Vector& ft, + int scomp, int dcomp, int ncomp, + const Geometry& cgeom, const Geometry& fgeom, + const IntVect& ratio, Interp* mapper, + const Vector& bcs, int bcscomp); #ifndef BL_NO_FORT enum InterpEM_t { InterpE, InterpB}; diff --git a/Src/AmrCore/AMReX_FillPatchUtil_I.H b/Src/AmrCore/AMReX_FillPatchUtil_I.H index fa26cc5ea9f..13e6ea89171 100644 --- a/Src/AmrCore/AMReX_FillPatchUtil_I.H +++ b/Src/AmrCore/AMReX_FillPatchUtil_I.H @@ -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) { + 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) @@ -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; @@ -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) { + 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); @@ -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 { @@ -504,10 +513,16 @@ namespace detail { auto solve_mask = make_mf_crse_mask(fpc, ncomp, mf.boxArray().ixType(), ratio); mf_set_domain_bndry(mf_crse_patch, cgeom); + if constexpr (std::is_same_v) { + 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) { + fbc.fp1_src_ghost = IntVect(0); + } FillPatchSingleLevel(mf_refined_patch, time, fmf, ft, scomp, 0, ncomp, fgeom, fbc, fbccomp); @@ -565,16 +580,29 @@ namespace detail { MF mf_crse_patch = make_mf_crse_patch(fpc, ncomp); mf_set_domain_bndry (mf_crse_patch, cgeom); + if constexpr (std::is_same_v) { + 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(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 + ) { + 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); @@ -583,6 +611,9 @@ namespace detail { } } + if constexpr(std::is_same_v) { + fbc.fp1_src_ghost = IntVect(0); + } FillPatchSingleLevel(mf, nghost, time, fmf, ft, scomp, dcomp, ncomp, fgeom, fbc, fbccomp); @@ -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); @@ -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) { + 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) { + 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); @@ -1075,6 +1120,8 @@ InterpFromCoarseLevel (Array 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; @@ -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& bcs) + const Vector& 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 +std::enable_if_t::value> +FillPatchSingleLevel (MF& mf, IntVect const& nghost, Real time, + const Vector& smf, IntVect const& snghost, + const Vector& 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 +std::enable_if_t::value> +FillPatchTwoLevels (MF& mf, IntVect const& nghost, + IntVect const& nghost_outside_domain, Real time, + const Vector& cmf, const Vector& ct, + const Vector& fmf, const Vector& ft, + int scomp, int dcomp, int ncomp, + const Geometry& cgeom, const Geometry& fgeom, + const IntVect& ratio, Interp* mapper, + const Vector& 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 diff --git a/Src/Base/AMReX_PhysBCFunct.H b/Src/Base/AMReX_PhysBCFunct.H index 699863f062a..c85c617a1f0 100644 --- a/Src/Base/AMReX_PhysBCFunct.H +++ b/Src/Base/AMReX_PhysBCFunct.H @@ -123,6 +123,56 @@ public: Real /*time*/, int /*bccomp*/) {} }; +class PhysBCFunctUseCoarseGhost +{ +public: + + PhysBCFunctUseCoarseGhost (IntVect const& a_fp1_src_ghost) + : fp1_src_ghost(a_fp1_src_ghost) {} + + template + PhysBCFunctUseCoarseGhost (MF const& cmf, IntVect const& a_nghost, + IntVect const& a_nghost_outside_domain, + IntVect const& ratio, Interp* mapper) + : nghost(a_nghost), + nghost_outside_domain(a_nghost_outside_domain), + cghost(cmf.nGrowVect()) + { + IndexType typ = cmf.ixType(); + + const auto& coarsener = mapper->BoxCoarsener(ratio); + Box tmp(-nghost, IntVect(32), typ); + Box tmp2 = coarsener.doit(tmp); + src_ghost = -tmp2.smallEnd(); + + tmp = Box(-nghost_outside_domain, IntVect(32), typ); + tmp2 = coarsener.doit(tmp); + src_ghost_outside_domain = -tmp2.smallEnd(); + + cghost.min(src_ghost); + AMREX_ALWAYS_ASSERT(cghost.allGE(src_ghost_outside_domain)); + } + + void operator() (MultiFab& /*mf*/, int /*dcomp*/, int /*ncomp*/, IntVect const& /*nghost*/, + Real /*time*/, int /*bccomp*/) {} + + IntVect nghost; // # of fine ghosts inside domain to be filled + + IntVect nghost_outside_domain; // # of fine ghosts outside domain to be filled + + // This is the number of coarse ghost cells needed to interpolate fine + // ghost cells inside the domain + IntVect src_ghost; + + // This is the number of coarse ghost cells needed to interpolate + // nghost_outside_domain fine ghost cells outside the domain + IntVect src_ghost_outside_domain; + + // This is the minimum number of ghost cells needed in coarse MF + IntVect cghost; + + IntVect fp1_src_ghost; // Used to pass information into FillPatchSingleLevel +}; template class PhysBCFunct