diff --git a/Src/AmrCore/AMReX_FillPatchUtil.H b/Src/AmrCore/AMReX_FillPatchUtil.H index 11a8eb560e..b71102b58d 100644 --- a/Src/AmrCore/AMReX_FillPatchUtil.H +++ b/Src/AmrCore/AMReX_FillPatchUtil.H @@ -728,6 +728,75 @@ namespace amrex const IntVect& ratio, Interp* mapper, const Vector& bcs); + /** + * \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. + */ + 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); + #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 fa26cc5ea9..4192d75934 100644 --- a/Src/AmrCore/AMReX_FillPatchUtil_I.H +++ b/Src/AmrCore/AMReX_FillPatchUtil_I.H @@ -1075,6 +1075,8 @@ InterpFromCoarseLevel (Array const& mf, IntVect const& ngho const PreInterpHook& pre_interp, const PostInterpHook& post_interp) { + BL_PROFILE("InterpFromCoarseLevel"); + using FAB = typename MF::FABType::value_type; using iFAB = typename iMultiFab::FABType::value_type; @@ -1219,6 +1221,8 @@ InterpFromCoarseLevel (MF& mf, IntVect const& nghost, const IntVect& ratio, Interp* mapper, const Vector& bcs) { + BL_PROFILE("InterpFromCoarseLevel_nobc"); + const BoxArray& ba = mf.boxArray(); const DistributionMapping& dm = mf.DistributionMap(); @@ -1244,11 +1248,203 @@ InterpFromCoarseLevel (MF& mf, IntVect const& nghost, cgeom.periodicity()); Box fdomain_g = amrex::convert(fgeom.Domain(),typ); - fdomain_g.grow(nghost_outside_domain); + for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { + if (fgeom.isPeriodic(idim)) { + fdomain_g.grow(idim, nghost[idim]); + } else { + fdomain_g.grow(idim, nghost_outside_domain[idim]); + } + } FillPatchInterp(mf, dcomp, mf_crse_patch, 0, ncomp, nghost, cgeom, fgeom, fdomain_g, ratio, mapper, bcs, 0); } +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) +{ + BL_PROFILE("FillPatchSingleLevel_nobc"); + + AMREX_ASSERT(scomp+ncomp <= smf[0]->nComp()); + AMREX_ASSERT(dcomp+ncomp <= mf.nComp()); + AMREX_ASSERT(smf.size() == stime.size()); + AMREX_ASSERT(!smf.empty()); + AMREX_ASSERT(nghost.allLE(mf.nGrowVect())); + + 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, snghost, nghost, geom.periodicity()); + } + } + else if (smf.size() == 2) + { + BL_ASSERT(smf[0]->boxArray() == smf[1]->boxArray()); + MF raii; + MF * dmf; + int destcomp; + bool sameba; + if (mf.boxArray() == smf[0]->boxArray() && + mf.DistributionMap() == smf[0]->DistributionMap()) + { + dmf = &mf; + destcomp = dcomp; + sameba = true; + } else { + raii.define(smf[0]->boxArray(), smf[0]->DistributionMap(), ncomp, snghost, + MFInfo(), smf[0]->Factory()); + + dmf = &raii; + destcomp = 0; + sameba = false; + } + + if ((dmf != smf[0] && dmf != smf[1]) || scomp != dcomp) + { + IntVect interp_ghost = snghost; + 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.growntilebox(interp_ghost); + const Real t0 = stime[0]; + const Real t1 = stime[1]; + auto const sfab0 = smf[0]->array(mfi); + auto const sfab1 = smf[1]->array(mfi); + auto dfab = dmf->array(mfi); + + if (time == t0) + { + AMREX_HOST_DEVICE_PARALLEL_FOR_4D ( bx, ncomp, i, j, k, n, + { + dfab(i,j,k,n+destcomp) = sfab0(i,j,k,n+scomp); + }); + } + else if (time == t1) + { + AMREX_HOST_DEVICE_PARALLEL_FOR_4D ( bx, ncomp, i, j, k, n, + { + dfab(i,j,k,n+destcomp) = sfab1(i,j,k,n+scomp); + }); + } + else if (! amrex::almostEqual(t0,t1)) + { + Real alpha = (t1-time)/(t1-t0); + Real beta = (time-t0)/(t1-t0); + AMREX_HOST_DEVICE_PARALLEL_FOR_4D ( bx, ncomp, i, j, k, n, + { + dfab(i,j,k,n+destcomp) = alpha*sfab0(i,j,k,n+scomp) + + beta*sfab1(i,j,k,n+scomp); + }); + } + else + { + AMREX_HOST_DEVICE_PARALLEL_FOR_4D ( bx, ncomp, i, j, k, n, + { + dfab(i,j,k,n+destcomp) = sfab0(i,j,k,n+scomp); + }); + } + } + } + + if (sameba) + { + // Note that when sameba is true mf's BoxArray is nonoverlapping. + // So FillBoundary is safe. + mf.FillBoundary(dcomp, ncomp, nghost, geom.periodicity()); + } + else + { + mf.ParallelCopy(*dmf, 0, dcomp, ncomp, snghost, nghost, geom.periodicity()); + } + } + else { + amrex::Abort("FillPatchSingleLevel: high-order interpolation in time not implemented yet"); + } +} + +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) +{ + BL_PROFILE("FillPatchTwoLevels_nobc"); + + const IndexType& typ = mf.ixType(); + + AMREX_ALWAYS_ASSERT(mf.ixType.nodeCentered() || mf.cellCentered()); + + if (nghost.max() > 0 || mf.getBDKey() != fmf[0]->getBDKey()) + { + const InterpolaterBoxCoarsener& coarsener = mapper->BoxCoarsener(ratio); + + Box tmp(-nghost, IntVect(32), typ); + Box tmp2 = coarsener.doit(tmp); + // This is the number of coarse ghost cells needed to interpolate + // nghost fine ghost cells inside the domain + IntVect src_ghost = -tmp2.smallEnd(); + + // This is the number of coarse ghost cells needed to interpolate + // nghost_outside_domain fine ghost cells outside the domain + tmp = Box(-nghost_outside_domain, IntVect(32), typ); + tmp2 = coarsener.doit(tmp); + IntVect src_ghost_outside_domain = -tmp2.smallEnd(); + AMREX_ALWAYS_ASSERT(src_ghost_outside_domain.allGE(src_ghost)); + + IntVect cghost = cmf.nGrowVect(); + cghost.min(src_ghost); + // This is the minimum number of ghost cells needed in cmf. + AMREX_ALWAYS_ASSERT(cghost.allGE(src_ghost_outside_domain)); + + const FabArrayBase::FPinfo& fpc = FabArrayBase::TheFPinfo(*fmf[0], mf, + nghost, + coarsener, + fgeom, + cgeom, + nullptr); + + if ( ! fpc.ba_crse_patch.empty()) + { + MF mf_crse_patch(fpc.ba_crse_patch, fpc.dm_patch, ncomp, 0); + + FillPatchSingleLevel(mf_crse_patch, IntVect(0), time, cmf, cghost, + ct, scomp, 0, ncomp, cgeom); + + MF mf_fine_patch(fpc.ba_fine_patch, fpc.dm_patch, ncomp, 0); + + Box fdomain_g = amrex::convert(fgeom.Domain(),typ); + for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { + if (fgeom.isPeriodic(idim)) { + fdomain_g.grow(idim, nghost[idim]); + } else { + fdomain_g.grow(idim, nghost_outside_domain[idim]); + } + } + FillPatchInterp(mf_fine_patch, 0, mf_crse_patch, 0, + ncomp, IntVect(0), cgeom, fgeom, + fdomain_g, ratio, mapper, bcs, 0); + + mf.ParallelCopy(mf_fine_patch, 0, dcomp, ncomp, IntVect{0}, nghost); + } + } + + FillPatchSingleLevel(mf, nghost, time, fmf, IntVect(0), ft, + scomp, dcomp, ncomp, fgeom); +} + template std::enable_if_t::value> FillPatchNLevels (MF& mf, int level, const IntVect& nghost, Real time,