From 95f1c2d4cd60c674660a50446f785ba5fa85fb90 Mon Sep 17 00:00:00 2001 From: Weiqun Zhang Date: Thu, 29 Sep 2022 11:04:58 -0700 Subject: [PATCH] FillPatcher class This adds a class FillPatcher for filling fine level data. It's not as general as the various FillPatch functions (e.g., FillPatchTwoLevels). However, it can reduce the amount of communication data. Suppose we use RK2 with subcycling and the refinement ratio is 2. For each step on level 0, there are two steps on level 1. With RK2, each fine step needs to call FillPatch twice. So the total number of FillPatch calls is 4 in the two fine steps. Using the free function, one ParallelCopy per FillPatch call is needed for copying coarse data for spatial interpolation. With the FillPatcher class, two ParallelCopy calls will be done to copy old and new coarse data. Then these data will be used in the four FillPatcher::fill calls. This new approach saves two ParallelCopy calls per coarse step for a two levels run. It could save more if the time stepping requires more substeps or the refinement ratio is higher. Note that many of our AMReX codes use a time stepping algorithm that needs only one FillPatch call per step. For those codes, this new approach will not save any communication for a refinement ratio of 2. However, it will save communication when the refinement ratio is 4. --- Src/Amr/AMReX_AmrLevel.H | 19 +- Src/Amr/AMReX_AmrLevel.cpp | 56 +++ Src/AmrCore/AMReX_FillPatchUtil.H | 11 +- Src/AmrCore/AMReX_FillPatcher.H | 331 ++++++++++++++++++ Src/AmrCore/CMakeLists.txt | 1 + Src/AmrCore/Make.package | 2 + Src/Base/AMReX_FArrayBox.H | 2 +- Src/Base/AMReX_Geometry.H | 4 + Src/Base/AMReX_Geometry.cpp | 20 +- .../Source/AdvancePhiAllLevels.cpp | 3 +- .../Source/AdvancePhiAtLevel.cpp | 3 +- .../Amr/Advection_AmrCore/Source/AmrCoreAdv.H | 11 +- .../Advection_AmrCore/Source/AmrCoreAdv.cpp | 64 +++- .../Source/Src_K/Make.package | 2 +- .../Advection_AmrLevel/Source/AmrLevelAdv.H | 2 +- .../Advection_AmrLevel/Source/AmrLevelAdv.cpp | 37 +- 16 files changed, 518 insertions(+), 50 deletions(-) create mode 100644 Src/AmrCore/AMReX_FillPatcher.H diff --git a/Src/Amr/AMReX_AmrLevel.H b/Src/Amr/AMReX_AmrLevel.H index 0aaf7fc2620..e4cd966d932 100644 --- a/Src/Amr/AMReX_AmrLevel.H +++ b/Src/Amr/AMReX_AmrLevel.H @@ -15,6 +15,7 @@ #include #include #include +#include #ifdef AMREX_USE_EB #include #endif @@ -243,12 +244,14 @@ public: Long countCells () const noexcept; //! Get the area not to tag. - const BoxArray& getAreaNotToTag() noexcept; - const Box& getAreaToTag() noexcept; + const BoxArray& getAreaNotToTag () noexcept; + const Box& getAreaToTag () noexcept; //! Construct the area not to tag. - void constructAreaNotToTag(); + void constructAreaNotToTag (); //! Set the area not to tag. - void setAreaNotToTag(BoxArray& ba) noexcept; + void setAreaNotToTag (BoxArray& ba) noexcept; + + void resetFillPatcher (); /** * \brief Error estimation for regridding. This is a pure virtual @@ -365,6 +368,10 @@ public: virtual void particle_redistribute (int /*lbase*/ = 0, bool /*a_init*/ = false) {;} #endif + // Fill with FillPatcher on level > 0 and AmrLevel::FillPatch on level 0. + void FillPatcherFill (amrex::MultiFab& mf, int nghost, amrex::Real time, + int state_index); + static void FillPatch (AmrLevel& amrlevel, MultiFab& leveldata, int boxGrow, @@ -425,7 +432,7 @@ protected: IntVect fine_ratio; // Refinement ratio to finer level. static DeriveList derive_lst; // List of derived quantities. static DescriptorList desc_lst; // List of state variables. - Vector state; // Array of state data. + Vector state; // Array of state data. BoxArray m_AreaNotToTag; //Area which shouldn't be tagged on this level. Box m_AreaToTag; //Area which is allowed to be tagged on this level. @@ -436,6 +443,8 @@ protected: std::unique_ptr > m_factory; + Vector>> m_fillpatcher; + private: mutable BoxArray edge_grids[AMREX_SPACEDIM]; // face-centered grids diff --git a/Src/Amr/AMReX_AmrLevel.cpp b/Src/Amr/AMReX_AmrLevel.cpp index a88489f9512..591d504d11e 100644 --- a/Src/Amr/AMReX_AmrLevel.cpp +++ b/Src/Amr/AMReX_AmrLevel.cpp @@ -102,6 +102,7 @@ AmrLevel::AmrLevel (Amr& papa, } state.resize(desc_lst.size()); + m_fillpatcher.resize(desc_lst.size()); #ifdef AMREX_USE_EB if (EB2::TopIndexSpaceIfPresent()) { @@ -451,6 +452,8 @@ AmrLevel::restart (Amr& papa, } } + m_fillpatcher.resize(ndesc); + if (parent->useFixedCoarseGrids()) constructAreaNotToTag(); post_step_regrid = 0; @@ -2096,6 +2099,59 @@ void AmrLevel::constructAreaNotToTag () } } +void +AmrLevel::resetFillPatcher () +{ + for (auto& fp : m_fillpatcher) { + fp.reset(); + } +} + +void +AmrLevel::FillPatcherFill (MultiFab& mf, int nghost, Real time, int state_index) +{ + if (level == 0) { + FillPatch(*this, mf, nghost, time, state_index, 0, mf.nComp()); + } else { + AmrLevel& fine_level = *this; + AmrLevel& crse_level = parent->getLevel(level-1); + const Geometry& geom_fine = fine_level.geom; + const Geometry& geom_crse = crse_level.geom; + + Vector smf_crse; + Vector stime_crse; + StateData& statedata_crse = crse_level.state[state_index]; + statedata_crse.getData(smf_crse,stime_crse,time); + StateDataPhysBCFunct physbcf_crse(statedata_crse,0,geom_crse); + + Vector smf_fine; + Vector stime_fine; + StateData& statedata_fine = fine_level.state[state_index]; + statedata_fine.getData(smf_fine,stime_fine,time); + StateDataPhysBCFunct physbcf_fine(statedata_fine,0,geom_fine); + + const StateDescriptor& desc = AmrLevel::desc_lst[state_index]; + + if (level > 1 &&!amrex::ProperlyNested(fine_level.crse_ratio, + parent->blockingFactor(fine_level.level), + nghost, mf.ixType(), desc.interp(0))) { + amrex::Abort("FillPatcherFill: Grids are not properly nested. Must increase blocking factor."); + } + + + auto& fillpatcher = m_fillpatcher[state_index]; + if (fillpatcher == nullptr) { + fillpatcher = std::make_unique> + (parent->boxArray(level), parent->DistributionMap(level), geom_fine, + parent->boxArray(level-1), parent->DistributionMap(level-1), geom_crse, + IntVect(nghost), desc.interp(0)); + } + + fillpatcher->fill(mf, time, smf_crse, stime_crse, smf_fine, stime_fine, + physbcf_crse, physbcf_fine, desc.getBCs()); + } +} + void AmrLevel::FillPatch (AmrLevel& amrlevel, MultiFab& leveldata, diff --git a/Src/AmrCore/AMReX_FillPatchUtil.H b/Src/AmrCore/AMReX_FillPatchUtil.H index 51a5f457391..495cbc180b6 100644 --- a/Src/AmrCore/AMReX_FillPatchUtil.H +++ b/Src/AmrCore/AMReX_FillPatchUtil.H @@ -28,12 +28,17 @@ namespace amrex { - template + template struct NullInterpHook { - void operator() (FAB& /*fab*/, const Box& /*bx*/, int /*icomp*/, int /*ncomp*/) const {} + template ::value,int> = 0> + void operator() (MFFAB& /*fab*/, const Box& /*bx*/, int /*icomp*/, int /*ncomp*/) const {} - void operator() (Array /*fab*/, const Box& /*bx*/, int /*icomp*/, int /*ncomp*/) const {} + template ::value,int> = 0> + void operator() (Array /*fab*/, const Box& /*bx*/, int /*icomp*/, int /*ncomp*/) const {} + + template ::value,int> = 0> + void operator() (MFFAB& /*mf*/, int /*icomp*/, int /*ncomp*/) const {} }; template diff --git a/Src/AmrCore/AMReX_FillPatcher.H b/Src/AmrCore/AMReX_FillPatcher.H new file mode 100644 index 00000000000..849c3f7db70 --- /dev/null +++ b/Src/AmrCore/AMReX_FillPatcher.H @@ -0,0 +1,331 @@ +#ifndef AMREX_FILLPATCHER_H_ +#define AMREX_FILLPATCHER_H_ +#include + +#include + +namespace amrex { + +/** + * \brief FillPatcher is for filling a fine level MultiFab/FabArray. + * + * This class is not as general as the FillPatchTwoLevels functions. It + * fills the fine ghost cells not overlapping any fine level valid cells + * with interpolation of the coarse data. Then it fills the fine ghost + * cells overlapping fine level valid cells with the fine level data. If + * the valid cells of the destination need to be filled, it will be done as + * well. Finally, it will fill the physical bounbary using the user + * provided functor. The `fill` member function can be used to do the + * operations just described. Alternatively, one can also use the + * `fillCoarseFineBounary` to fill the ghost cells at the coarse/fine + * boundary only. Then one can manually call FillBoundary to fill the other + * ghost cells, and use the physical BC functor to handle the physical + * boundeary. + * + * The communication of the coarse data needed for spatial interpolation is + * optimized at the cost of being error-prone. One must follow the + * following guidelines. + * + * (1) This class is for filling data during time stepping, not during + * regrid. The fine level data passed as input must have the same BoxArray + * and DistributionMapping as the destination. It's OK they are the same + * MultiFab. For AmrLevel based codes, AmrLevel::FillPatcherFill wil try to + * use FillPatcher if it can, and AmrLevel::Patch will use the fillpatch + * functions. + * + * (2) When to build? It is recommended that one uses `std::unique_ptr` to + * store the FillPatcher object, and build it only when it is needed and + * it's a nullptr. For AmrLevel based codes, the AmrLevel class will build + * it for you as needed when you call the AmrLevel::FillPatcherFill + * function. + * + * (3) When to destroy? Usually, we do time steppig on a coarse level + * first. Then we recursively do time stepping on fine levels. After the + * finer level finishes, we do reflux and averge the fine data down to the + * coarse level. After that we should destroy the FillPatcher object + * associated with these two levels, because the coarse data stored in the + * object has become outdated. For AmrCore based codes, you could use + * Tests/Amr/Advection_AmrCore as an example. For AmrLevel based codes, you + * should do this in the post_timestep virtual function (see + * Tests/Amr/Advection_AmrLevel for an example). + * + * (4) This class will try to fill all the components. Unlike the + * FillPatchTwoLevels funcitons, one cannot specify the number of + * components. If you do not need to fill all the components, you should + * create alias MulitFabs/FabArrays to work aroun it. + * + * (5) This only works for cell-centered and nodal data. + */ + +template +class FillPatcher +{ +public: + + /** + * \brief Constructor of FillPatcher + * + * \param fba fine level BoxArray + * \param fdm fine level DistributionMapping + * \param fgeom fine level DistributionMapping + * \param cba coarse level BoxArray + * \param cdm coarse level DistributionMapping + * \param cgeom coarse level DistributionMapping + * \param ngost max number of ghost cells to be filled at coarse/fine boundary + * \param interp for spatial interpolation + * \param eb_index_space optional argument for specifying EB IndexSpace + */ + FillPatcher (BoxArray const& fba, DistributionMapping const& fdm, + Geometry const& fgeom, + BoxArray const& cba, DistributionMapping const& cdm, + Geometry const& cgeom, + IntVect const& nghost, InterpBase* interp, +#ifdef AMREX_USE_EB + EB2::IndexSpace const* eb_index_space = EB2::TopIndexSpaceIfPresent()); +#else + EB2::IndexSpace const* eb_index_space = nullptr); +#endif + + /** + * \brief Function to fill data + * + * \param mf destination MultiFab/FabArray + * \param time time associated with the destination + * \param crse_data coarse level data + * \param crse_time time associated with the coarse data + * \param fine_data fine level data + * \param fine_time time associated with the fine data + * \param cbc for filling coarse level physical BC + * \param fbc for filling fine level physical BC + * \param bcs BCRec specifying physical boundary types + * \param nghost optional argument for specify number of ghost cell smaller + * than what's provided in the constructor + * \param pre_interp optional pre-interpolation hook for modifying the coarse data + * \param post_interp optional post-interpolation hook for modifying the fine data + */ + template , + typename PostInterpHook=NullInterpHook > + void fill (MF& mf, Real time, + Vector const& crse_data, Vector const& crse_time, + Vector const& fine_data, Vector const& fine_time, + BC& cbc, BC& fbc, Vector const& bcs, + IntVect const& nghost = IntVect(-1), + PreInterpHook const& pre_interp = {}, + PostInterpHook const& post_interp = {}); + + /** + * \brief Function to fill data at coarse/fine boundary only + * + * \param mf destination MultiFab/FabArray + * \param time time associated with the destination + * \param crse_data coarse level data + * \param crse_time time associated with the coarse data + * \param cbc for filling coarse level physical BC + * \param bcs BCRec specifying physical boundary types + * \param nghost optional argument for specify number of ghost cell smaller + * than what's provided in the constructor + * \param pre_interp optional pre-interpolation hook for modifying the coarse data + * \param post_interp optional post-interpolation hook for modifying the fine data + */ + template , + typename PostInterpHook=NullInterpHook > + void fillCoarseFineBoundary (MF& mf, Real time, + Vector const& crse_data, + Vector const& crse_time, + BC& cbc, Vector const& bcs, + IntVect const& nghost = IntVect(-1), + PreInterpHook const& pre_interp = {}, + PostInterpHook const& post_interp = {}); + +private: + + BoxArray m_fba; + BoxArray m_cba; + DistributionMapping m_fdm; + DistributionMapping m_cdm; + Geometry m_fgeom; + Geometry m_cgeom; + IntVect m_nghost; + InterpBase* m_interp; + EB2::IndexSpace const* m_eb_index_space = nullptr; + Vector>> m_cf_crse_data; + std::unique_ptr m_cf_fine_data; +}; + +template +FillPatcher::FillPatcher (BoxArray const& fba, DistributionMapping const& fdm, + Geometry const& fgeom, + BoxArray const& cba, DistributionMapping const& cdm, + Geometry const& cgeom, + IntVect const& nghost, InterpBase* interp, + EB2::IndexSpace const* eb_index_space) + : m_fba(fba), + m_cba(cba), + m_fdm(fdm), + m_cdm(cdm), + m_fgeom(fgeom), + m_cgeom(cgeom), + m_nghost(nghost), + m_interp(interp), + m_eb_index_space(eb_index_space) +{ + static_assert(IsFabArray::value, + "FillPatcher: MF must be FabArray type"); + AMREX_ALWAYS_ASSERT(m_fba.ixType().cellCentered() || m_fba.ixType().nodeCentered()); +} + +template +template +void +FillPatcher::fill (MF& mf, Real time, + Vector const& cmf, Vector const& ct, + Vector const& fmf, Vector const& ft, + BC& cbc, BC& fbc, Vector const& bcs, + IntVect const& a_nghost, + PreInterpHook const& pre_interp, + PostInterpHook const& post_interp) +{ + BL_PROFILE("FillPatcher::fill()"); + + IntVect nghost = (a_nghost == -1) ? m_nghost : a_nghost; + nghost.min(mf.nGrowVect()); + + AMREX_ASSERT(m_fba == fmf[0]->boxArray() && + m_fdm == fmf[0]->DistributionMap() && + mf.nComp() == fmf[0]->nComp()); + + fillCoarseFineBoundary(mf, time, cmf, ct, cbc, bcs, nghost, + pre_interp, post_interp); + + FillPatchSingleLevel(mf, nghost, time, fmf, ft, 0, 0, mf.nComp(), + m_fgeom, fbc, 0); +} + +template +template +void +FillPatcher::fillCoarseFineBoundary (MF& mf, Real time, + Vector const& cmf, + Vector const& ct, + BC& cbc, Vector const& bcs, + IntVect const& a_nghost, + PreInterpHook const& pre_interp, + PostInterpHook const& post_interp) +{ + BL_PROFILE("FillPatcher::fillCFB"); + + IntVect nghost = (a_nghost == -1) ? m_nghost : a_nghost; + nghost.min(mf.nGrowVect()); + AMREX_ALWAYS_ASSERT(nghost.allLE(m_nghost)); + + const int ncomp = mf.nComp(); + + AMREX_ASSERT(m_fba == mf.boxArray() && + m_fdm == mf.DistributionMap() && + m_cba == cmf[0]->boxArray() && + m_cdm == cmf[0]->DistributionMap() && + ncomp == cmf[0]->nComp()); + + IntVect ratio; + for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { + ratio[idim] = m_fgeom.Domain().length(idim) / m_cgeom.Domain().length(idim); + } + AMREX_ASSERT(m_fgeom.Domain() == amrex::refine(m_cgeom.Domain(),ratio)); + + const InterpolaterBoxCoarsener& coarsener = m_interp->BoxCoarsener(ratio); + const FabArrayBase::FPinfo& fpc = FabArrayBase::TheFPinfo(mf, mf, + m_nghost, + coarsener, + m_fgeom, + m_cgeom, + m_eb_index_space); + + if ( ! fpc.ba_crse_patch.empty()) + { + if (m_cf_fine_data == nullptr) { + m_cf_fine_data = std::make_unique + (make_mf_fine_patch(fpc, ncomp)); + } + + int ncmfs = cmf.size(); + for (int icmf = 0; icmf < ncmfs; ++icmf) { + Real t = ct[icmf]; + auto it = std::find_if(m_cf_crse_data.begin(), m_cf_crse_data.end(), + [=] (auto const& x) { + return amrex::almostEqual(x.first,t); + }); + + if (it == std::end(m_cf_crse_data)) { + MF mf_crse_patch = make_mf_crse_patch(fpc, ncomp); + mf_crse_patch.ParallelCopy(*cmf[icmf], m_cgeom.periodicity()); + + std::pair> tmp; + tmp.first = t; + tmp.second = std::make_unique(std::move(mf_crse_patch)); + m_cf_crse_data.push_back(std::move(tmp)); + } + } + + MF mf_crse_patch; + if (m_cf_crse_data.size() > 0 && + amrex::almostEqual(time, m_cf_crse_data[0].first)) + { + mf_crse_patch = MF(*m_cf_crse_data[0].second, amrex::make_alias, + 0, ncomp); + } + else if (m_cf_crse_data.size() > 1 && + amrex::almostEqual(time, m_cf_crse_data[1].first)) + { + mf_crse_patch = MF(*m_cf_crse_data[1].second, amrex::make_alias, + 0, ncomp); + } + else if (m_cf_crse_data.size() == 2) + { + mf_crse_patch = make_mf_crse_patch(fpc, ncomp); + int const ng_space_interp = 8; // Need to be big enough + Box domain = m_cgeom.growPeriodicDomain(ng_space_interp); + domain.convert(mf.ixType()); + Real t0 = m_cf_crse_data[0].first; + Real t1 = m_cf_crse_data[1].first; + Real alpha = (t1-time)/(t1-t0); + Real beta = (time-t0)/(t1-t0); + AMREX_ASSERT(alpha >= 0._rt && beta >= 0._rt); + auto const& a = mf_crse_patch.arrays(); + auto const& a0 = m_cf_crse_data[0].second->const_arrays(); + auto const& a1 = m_cf_crse_data[1].second->const_arrays(); + amrex::ParallelFor(mf_crse_patch, IntVect(0), ncomp, + [=] AMREX_GPU_DEVICE (int bi, int i, int j, int k, int n) noexcept + { + if (domain.contains(i,j,k)) { + a[bi](i,j,k,n) = alpha*a0[bi](i,j,k,n) + + beta*a1[bi](i,j,k,n); + } + }); + Gpu::streamSynchronize(); + } + else + { + amrex::Abort("FillPatcher: High order interpolation in time not supported. Or FillPatcher was not properly deleted."); + } + + cbc(mf_crse_patch, 0, ncomp, nghost, time, 0); + + pre_interp(mf_crse_patch, 0, ncomp); + + FillPatchInterp(*m_cf_fine_data, 0, mf_crse_patch, 0, + ncomp, IntVect(0), m_cgeom, m_fgeom, + amrex::grow(amrex::convert(m_fgeom.Domain(), + mf.ixType()),nghost), + ratio, m_interp, bcs, 0); + + post_interp(*m_cf_fine_data, 0, ncomp); + + mf.ParallelCopy(*m_cf_fine_data, 0, 0, ncomp, IntVect{0}, nghost); + } +} + +} + +#endif diff --git a/Src/AmrCore/CMakeLists.txt b/Src/AmrCore/CMakeLists.txt index f9ff24f243b..be7c87eee4f 100644 --- a/Src/AmrCore/CMakeLists.txt +++ b/Src/AmrCore/CMakeLists.txt @@ -12,6 +12,7 @@ target_sources(amrex AMReX_FluxRegister.cpp AMReX_FillPatchUtil.H AMReX_FillPatchUtil_I.H + AMReX_FillPatcher.H AMReX_FluxRegister.H AMReX_InterpBase.H AMReX_InterpBase.cpp diff --git a/Src/AmrCore/Make.package b/Src/AmrCore/Make.package index 5b3afa61ccb..df3c2e83d40 100644 --- a/Src/AmrCore/Make.package +++ b/Src/AmrCore/Make.package @@ -6,6 +6,8 @@ CEXE_sources += AMReX_AmrCore.cpp AMReX_Cluster.cpp AMReX_ErrorList.cpp AMReX_Fi AMReX_Interpolater.cpp AMReX_MFInterpolater.cpp AMReX_TagBox.cpp AMReX_AmrMesh.cpp \ AMReX_InterpBase.cpp +CEXE_headers += AMReX_FillPatcher.H + CEXE_headers += AMReX_Interp_C.H AMReX_Interp_$(DIM)D_C.H CEXE_headers += AMReX_MFInterp_C.H AMReX_MFInterp_$(DIM)D_C.H diff --git a/Src/Base/AMReX_FArrayBox.H b/Src/Base/AMReX_FArrayBox.H index 3d3cda3674b..b678986c0e9 100644 --- a/Src/Base/AMReX_FArrayBox.H +++ b/Src/Base/AMReX_FArrayBox.H @@ -272,7 +272,7 @@ public: virtual ~FArrayBox () noexcept override {} FArrayBox (FArrayBox&& rhs) noexcept = default; - FArrayBox& operator= (FArrayBox&&) = default; + FArrayBox& operator= (FArrayBox&&) noexcept = default; FArrayBox (const FArrayBox&) = delete; FArrayBox& operator= (const FArrayBox&) = delete; diff --git a/Src/Base/AMReX_Geometry.H b/Src/Base/AMReX_Geometry.H index 0e0a49f540e..890ec2e0f7e 100644 --- a/Src/Base/AMReX_Geometry.H +++ b/Src/Base/AMReX_Geometry.H @@ -420,9 +420,13 @@ public: const Box& src, Vector& out) const noexcept; + //! Return domain box with non-periodic directions grown by ngrow. + Box growNonPeriodicDomain (IntVect const& ngrow) const noexcept; //! Return domain box with non-periodic directions grown by ngrow. Box growNonPeriodicDomain (int ngrow) const noexcept; //! Return domain box with periodic directions grown by ngrow. + Box growPeriodicDomain (IntVect const& ngrow) const noexcept; + //! Return domain box with periodic directions grown by ngrow. Box growPeriodicDomain (int ngrow) const noexcept; //! Set periodicity flags and return the old flags. diff --git a/Src/Base/AMReX_Geometry.cpp b/Src/Base/AMReX_Geometry.cpp index 2f80f2eb947..235c7bb7674 100644 --- a/Src/Base/AMReX_Geometry.cpp +++ b/Src/Base/AMReX_Geometry.cpp @@ -473,29 +473,41 @@ Geometry::periodicShift (const Box& target, } Box -Geometry::growNonPeriodicDomain (int ngrow) const noexcept +Geometry::growNonPeriodicDomain (IntVect const& ngrow) const noexcept { Box b = Domain(); for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { if (!isPeriodic(idim)) { - b.grow(idim,ngrow); + b.grow(idim,ngrow[idim]); } } return b; } Box -Geometry::growPeriodicDomain (int ngrow) const noexcept +Geometry::growPeriodicDomain (IntVect const& ngrow) const noexcept { Box b = Domain(); for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { if (isPeriodic(idim)) { - b.grow(idim,ngrow); + b.grow(idim,ngrow[idim]); } } return b; } +Box +Geometry::growNonPeriodicDomain (int ngrow) const noexcept +{ + return growNonPeriodicDomain(IntVect(ngrow)); +} + +Box +Geometry::growPeriodicDomain (int ngrow) const noexcept +{ + return growPeriodicDomain(IntVect(ngrow)); +} + void Geometry::computeRoundoffDomain () { diff --git a/Tests/Amr/Advection_AmrCore/Source/AdvancePhiAllLevels.cpp b/Tests/Amr/Advection_AmrCore/Source/AdvancePhiAllLevels.cpp index b5e48e6e409..4f97cbf3184 100644 --- a/Tests/Amr/Advection_AmrCore/Source/AdvancePhiAllLevels.cpp +++ b/Tests/Amr/Advection_AmrCore/Source/AdvancePhiAllLevels.cpp @@ -35,7 +35,8 @@ AmrCoreAdv::AdvancePhiAllLevels (Real time, Real dt_lev, int /*iteration*/) // State with ghost cells MultiFab Sborder(grids[lev], dmap[lev], phi_new[lev].nComp(), num_grow); - FillPatch(lev, time, Sborder, 0, Sborder.nComp()); + FillPatch(lev, time, Sborder, 0, Sborder.nComp(), + FillPatchType::fillpatch_function); #ifdef AMREX_USE_OMP #pragma omp parallel if (Gpu::notInLaunchRegion()) diff --git a/Tests/Amr/Advection_AmrCore/Source/AdvancePhiAtLevel.cpp b/Tests/Amr/Advection_AmrCore/Source/AdvancePhiAtLevel.cpp index 3ddd055eda0..7a5e1abbaa7 100644 --- a/Tests/Amr/Advection_AmrCore/Source/AdvancePhiAtLevel.cpp +++ b/Tests/Amr/Advection_AmrCore/Source/AdvancePhiAtLevel.cpp @@ -33,7 +33,8 @@ AmrCoreAdv::AdvancePhiAtLevel (int lev, Real time, Real dt_lev, int /*iteration* // State with ghost cells MultiFab Sborder(grids[lev], dmap[lev], S_new.nComp(), num_grow); - FillPatch(lev, time, Sborder, 0, Sborder.nComp()); + FillPatch(lev, time, Sborder, 0, Sborder.nComp(), + FillPatchType::fillpatch_class); #ifdef AMREX_USE_OMP #pragma omp parallel if (Gpu::notInLaunchRegion()) diff --git a/Tests/Amr/Advection_AmrCore/Source/AmrCoreAdv.H b/Tests/Amr/Advection_AmrCore/Source/AmrCoreAdv.H index e330d30e740..1b6832d8663 100644 --- a/Tests/Amr/Advection_AmrCore/Source/AmrCoreAdv.H +++ b/Tests/Amr/Advection_AmrCore/Source/AmrCoreAdv.H @@ -10,6 +10,7 @@ #include #include #include +#include #ifdef AMREX_USE_OMP # include @@ -98,15 +99,18 @@ private: // more flexible version of AverageDown() that lets you average down across multiple levels void AverageDownTo (int crse_lev); + enum class FillPatchType { fillpatch_class, fillpatch_function }; + // compute a new multifab by coping in phi from valid region and filling ghost cells // works for single level and 2-level cases (fill fine grid ghost by interpolating from coarse) - void FillPatch (int lev, amrex::Real time, amrex::MultiFab& mf, int icomp, int ncomp); + void FillPatch (int lev, amrex::Real time, amrex::MultiFab& mf, int icomp, + int ncomp, FillPatchType fptype); // fill an entire multifab by interpolating from the coarser level // this comes into play when a new level of refinement appears void FillCoarsePatch (int lev, amrex::Real time, amrex::MultiFab& mf, int icomp, int ncomp); - // utility to copy in data from phi_old and/or phi_new into another multifab + // Pack pointers to phi_old and/or phi_new and associated times. void GetData (int lev, amrex::Real time, amrex::Vector& data, amrex::Vector& datatime); @@ -165,6 +169,9 @@ private: // used in the reflux operation amrex::Vector > flux_reg; + // This is for fillpatch during timestepping, but not for regridding. + amrex::Vector>> fillpatcher; + // Velocity on all faces at all levels amrex::Vector< amrex::Array > facevel; diff --git a/Tests/Amr/Advection_AmrCore/Source/AmrCoreAdv.cpp b/Tests/Amr/Advection_AmrCore/Source/AmrCoreAdv.cpp index 62c9dc7417e..d3482264350 100644 --- a/Tests/Amr/Advection_AmrCore/Source/AmrCoreAdv.cpp +++ b/Tests/Amr/Advection_AmrCore/Source/AmrCoreAdv.cpp @@ -2,7 +2,6 @@ #include #include #include -#include #include #include #include @@ -95,6 +94,10 @@ AmrCoreAdv::AmrCoreAdv () // with the lev/lev-1 interface (and has grid spacing associated with lev-1) // therefore flux_reg[0] is never actually used in the reflux operation flux_reg.resize(nlevs_max+1); + + // fillpatcher[lev] is for filling data on level lev using the data on + // lev-1 and lev. + fillpatcher.resize(nlevs_max+1); } AmrCoreAdv::~AmrCoreAdv () @@ -230,7 +233,8 @@ AmrCoreAdv::RemakeLevel (int lev, Real time, const BoxArray& ba, MultiFab new_state(ba, dm, ncomp, ng); MultiFab old_state(ba, dm, ncomp, ng); - FillPatch(lev, time, new_state, 0, ncomp); + // Must use fillpatch_function + FillPatch(lev, time, new_state, 0, ncomp, FillPatchType::fillpatch_function); std::swap(new_state, phi_new[lev]); std::swap(old_state, phi_old[lev]); @@ -257,6 +261,7 @@ AmrCoreAdv::ClearLevel (int lev) phi_new[lev].clear(); phi_old[lev].clear(); flux_reg[lev].reset(nullptr); + fillpatcher[lev].reset(nullptr); } // Make a new level from scratch using provided BoxArray and DistributionMapping. @@ -418,7 +423,8 @@ AmrCoreAdv::AverageDownTo (int crse_lev) // compute a new multifab by coping in phi from valid region and filling ghost cells // works for single level and 2-level cases (fill fine grid ghost by interpolating from coarse) void -AmrCoreAdv::FillPatch (int lev, Real time, MultiFab& mf, int icomp, int ncomp) +AmrCoreAdv::FillPatch (int lev, Real time, MultiFab& mf, int icomp, int ncomp, + FillPatchType fptype) { if (lev == 0) { @@ -450,16 +456,32 @@ AmrCoreAdv::FillPatch (int lev, Real time, MultiFab& mf, int icomp, int ncomp) Interpolater* mapper = &cell_cons_interp; + if (fptype == FillPatchType::fillpatch_class) { + if (fillpatcher[lev] == nullptr) { + fillpatcher[lev] = std::make_unique> + (boxArray(lev ), DistributionMap(lev ), Geom(lev ), + boxArray(lev-1), DistributionMap(lev-1), Geom(lev-1), + mf.nGrowVect(), mapper); + } + } + if(Gpu::inLaunchRegion()) { GpuBndryFuncFab gpu_bndry_func(AmrCoreFill{}); PhysBCFunct > cphysbc(geom[lev-1],bcs,gpu_bndry_func); PhysBCFunct > fphysbc(geom[lev],bcs,gpu_bndry_func); - amrex::FillPatchTwoLevels(mf, time, cmf, ctime, fmf, ftime, - 0, icomp, ncomp, geom[lev-1], geom[lev], - cphysbc, 0, fphysbc, 0, refRatio(lev-1), - mapper, bcs, 0); + if (fptype == FillPatchType::fillpatch_class) { + AMREX_ASSERT(icomp == 0); // Otherwise we will have to make some + // alias multifabs. + fillpatcher[lev]->fill(mf, time, cmf, ctime, fmf, ftime, + cphysbc, fphysbc, bcs); + } else { + amrex::FillPatchTwoLevels(mf, time, cmf, ctime, fmf, ftime, + 0, icomp, ncomp, geom[lev-1], geom[lev], + cphysbc, 0, fphysbc, 0, refRatio(lev-1), + mapper, bcs, 0); + } } else { @@ -467,10 +489,17 @@ AmrCoreAdv::FillPatch (int lev, Real time, MultiFab& mf, int icomp, int ncomp) PhysBCFunct cphysbc(geom[lev-1],bcs,bndry_func); PhysBCFunct fphysbc(geom[lev],bcs,bndry_func); - amrex::FillPatchTwoLevels(mf, time, cmf, ctime, fmf, ftime, - 0, icomp, ncomp, geom[lev-1], geom[lev], - cphysbc, 0, fphysbc, 0, refRatio(lev-1), - mapper, bcs, 0); + if (fptype == FillPatchType::fillpatch_class) { + AMREX_ASSERT(icomp == 0); // Otherwise we will have to make some + // alias multifabs. + fillpatcher[lev]->fill(mf, time, cmf, ctime, fmf, ftime, + cphysbc, fphysbc, bcs); + } else { + amrex::FillPatchTwoLevels(mf, time, cmf, ctime, fmf, ftime, + 0, icomp, ncomp, geom[lev-1], geom[lev], + cphysbc, 0, fphysbc, 0, refRatio(lev-1), + mapper, bcs, 0); + } } } } @@ -513,21 +542,18 @@ AmrCoreAdv::FillCoarsePatch (int lev, Real time, MultiFab& mf, int icomp, int nc } } -// utility to copy in data from phi_old and/or phi_new into another multifab void AmrCoreAdv::GetData (int lev, Real time, Vector& data, Vector& datatime) { data.clear(); datatime.clear(); - const Real teps = (t_new[lev] - t_old[lev]) * 1.e-3; - - if (time > t_new[lev] - teps && time < t_new[lev] + teps) + if (amrex::almostEqual(time, t_new[lev])) { data.push_back(&phi_new[lev]); datatime.push_back(t_new[lev]); } - else if (time > t_old[lev] - teps && time < t_old[lev] + teps) + else if (amrex::almostEqual(time, t_old[lev])) { data.push_back(&phi_old[lev]); datatime.push_back(t_old[lev]); @@ -631,6 +657,8 @@ AmrCoreAdv::timeStepWithSubcycling (int lev, Real time, int iteration) } AverageDownTo(lev); // average lev+1 down to lev + + fillpatcher[lev+1].reset(); // Because the data on lev have changed. } @@ -694,6 +722,10 @@ AmrCoreAdv::timeStepNoSubcycling (Real time, int iteration) // Make sure the coarser levels are consistent with the finer levels AverageDown (); + for (auto& fp : fillpatcher) { + fp.reset(); // Because the data have changed. + } + for (int lev = 0; lev <= finest_level; lev++) ++istep[lev]; diff --git a/Tests/Amr/Advection_AmrCore/Source/Src_K/Make.package b/Tests/Amr/Advection_AmrCore/Source/Src_K/Make.package index e98f493727c..5254ff6f63f 100644 --- a/Tests/Amr/Advection_AmrCore/Source/Src_K/Make.package +++ b/Tests/Amr/Advection_AmrCore/Source/Src_K/Make.package @@ -1,3 +1,3 @@ CEXE_headers += Adv_K.H -CEXE_headers += compute_flux_K_$(DIM).H +CEXE_headers += compute_flux_$(DIM)D_K.H CEXE_headers += slope_K.H diff --git a/Tests/Amr/Advection_AmrLevel/Source/AmrLevelAdv.H b/Tests/Amr/Advection_AmrLevel/Source/AmrLevelAdv.H index 1e5bacbc497..faf56357e29 100644 --- a/Tests/Amr/Advection_AmrLevel/Source/AmrLevelAdv.H +++ b/Tests/Amr/Advection_AmrLevel/Source/AmrLevelAdv.H @@ -231,7 +231,7 @@ protected: /* * The data. */ - amrex::FluxRegister* flux_reg; + std::unique_ptr flux_reg; /* * Static data members. diff --git a/Tests/Amr/Advection_AmrLevel/Source/AmrLevelAdv.cpp b/Tests/Amr/Advection_AmrLevel/Source/AmrLevelAdv.cpp index db69749a85f..5c2cb99a3ac 100644 --- a/Tests/Amr/Advection_AmrLevel/Source/AmrLevelAdv.cpp +++ b/Tests/Amr/Advection_AmrLevel/Source/AmrLevelAdv.cpp @@ -36,7 +36,6 @@ int AmrLevelAdv::do_tracers = 0; */ AmrLevelAdv::AmrLevelAdv () { - flux_reg = 0; } /** @@ -51,9 +50,9 @@ AmrLevelAdv::AmrLevelAdv (Amr& papa, : AmrLevel(papa,lev,level_geom,bl,dm,time) { - flux_reg = 0; - if (level > 0 && do_reflux) - flux_reg = new FluxRegister(grids,dmap,crse_ratio,level,NUM_STATE); + if (level > 0 && do_reflux) { + flux_reg = std::make_unique(grids,dmap,crse_ratio,level,NUM_STATE); + } } /** @@ -61,7 +60,6 @@ AmrLevelAdv::AmrLevelAdv (Amr& papa, */ AmrLevelAdv::~AmrLevelAdv () { - delete flux_reg; } /** @@ -74,9 +72,9 @@ AmrLevelAdv::restart (Amr& papa, { AmrLevel::restart(papa,is,bReadSpecial); - BL_ASSERT(flux_reg == 0); - if (level > 0 && do_reflux) - flux_reg = new FluxRegister(grids,dmap,crse_ratio,level,NUM_STATE); + if (level > 0 && do_reflux) { + flux_reg = std::make_unique(grids,dmap,crse_ratio,level,NUM_STATE); + } } /** @@ -88,11 +86,11 @@ AmrLevelAdv::checkPoint (const std::string& dir, VisMF::How how, bool dump_old) { - AmrLevel::checkPoint(dir, os, how, dump_old); + AmrLevel::checkPoint(dir, os, how, dump_old); #ifdef AMREX_PARTICLES - if (do_tracers && level == 0) { - TracerPC->WritePlotFile(dir, "Tracer"); - } + if (do_tracers && level == 0) { + TracerPC->WritePlotFile(dir, "Tracer"); + } #endif } @@ -285,7 +283,8 @@ AmrLevelAdv::advance (Real time, // State with ghost cells MultiFab Sborder(grids, dmap, NUM_STATE, NUM_GROW); - FillPatch(*this, Sborder, NUM_GROW, time, Phi_Type, 0, NUM_STATE); + // We use FillPatcher to do fillpatch here if we can + FillPatcherFill(Sborder, NUM_GROW, time, Phi_Type); // MF to hold the mac velocity MultiFab Umac[BL_SPACEDIM]; @@ -601,11 +600,19 @@ AmrLevelAdv::post_timestep (int iteration) // int finest_level = parent->finestLevel(); - if (do_reflux && level < finest_level) + if (do_reflux && level < finest_level) { reflux(); + } - if (level < finest_level) + if (level < finest_level) { avgDown(); + } + + if (level < finest_level) { + // fillpatcher on level+1 needs to be reset because data on this + // level have changed. + getLevel(level+1).resetFillPatcher(); + } #ifdef AMREX_PARTICLES if (TracerPC)