From a5997b9cc51e7caebdc99b975ff3096da68ae5c6 Mon Sep 17 00:00:00 2001 From: Weiqun Zhang Date: Sun, 9 Jul 2023 12:39:08 -0700 Subject: [PATCH] New feature in EB: Add regular coarse levels This feature could be useful when the EB is so thin that it is only resolved at the finest level and the usual way of generate coarse levels fails with multi-cut or multi-value cells. When the EB is fully refined to the finest level and there is no subcycling, no EB information is really needed at coarse levels. amrex::TagCutCells has been updated to work with the new feature. --- Src/EB/AMReX_EB2.H | 4 + Src/EB/AMReX_EB2.cpp | 8 ++ Src/EB/AMReX_EB2_IndexSpaceI.H | 46 +++++++++++ Src/EB/AMReX_EB2_IndexSpace_STL.H | 1 + Src/EB/AMReX_EB2_IndexSpace_STL.cpp | 6 ++ Src/EB/AMReX_EB2_IndexSpace_chkpt_file.H | 1 + Src/EB/AMReX_EB2_IndexSpace_chkpt_file.cpp | 6 ++ Src/EB/AMReX_EB2_Level.H | 28 ++++++- Src/EB/AMReX_EB2_Level.cpp | 93 +++++++++++++++++++++- Src/EB/AMReX_EBAmrUtil.cpp | 43 +++++++--- Src/EB/AMReX_EBDataCollection.H | 5 ++ Src/EB/AMReX_EBDataCollection.cpp | 14 ++++ Src/EB/AMReX_EBFabFactory.H | 6 ++ Src/EB/AMReX_EBFabFactory.cpp | 6 ++ 14 files changed, 252 insertions(+), 15 deletions(-) diff --git a/Src/EB/AMReX_EB2.H b/Src/EB/AMReX_EB2.H index a73bf5bdd3d..f76fc8ddafe 100644 --- a/Src/EB/AMReX_EB2.H +++ b/Src/EB/AMReX_EB2.H @@ -56,6 +56,7 @@ public: [[nodiscard]] virtual const Geometry& getGeometry (const Box& domain) const = 0; [[nodiscard]] virtual const Box& coarsestDomain () const = 0; virtual void addFineLevels (int num_new_fine_levels) = 0; + virtual void addRegularCoarseLevels (int num_new_coarse_levels) = 0; protected: static AMREX_EXPORT Vector > m_instance; @@ -88,6 +89,7 @@ public: return m_geom.back().Domain(); } void addFineLevels (int num_new_fine_levels) final; + void addRegularCoarseLevels (int num_new_coarse_levels) final; using F = typename G::FunctionType; @@ -147,6 +149,8 @@ int maxCoarseningLevel (IndexSpace const* ebis, const Geometry& geom); void addFineLevels (int num_new_fine_levels); +void addRegularCoarseLevels (int num_new_coarse_levels); + } #endif diff --git a/Src/EB/AMReX_EB2.cpp b/Src/EB/AMReX_EB2.cpp index fab107dabb7..8786da3c124 100644 --- a/Src/EB/AMReX_EB2.cpp +++ b/Src/EB/AMReX_EB2.cpp @@ -240,6 +240,14 @@ void addFineLevels (int num_new_fine_levels) } } +void addRegularCoarseLevels (int num_new_coarse_levels) +{ + auto *p = const_cast(TopIndexSpace()); + if (p) { + p->addRegularCoarseLevels(num_new_coarse_levels); + } +} + void BuildFromChkptFile (std::string const& fname, const Geometry& geom, int required_coarsening_level, diff --git a/Src/EB/AMReX_EB2_IndexSpaceI.H b/Src/EB/AMReX_EB2_IndexSpaceI.H index e7db810b03b..e3b9960538c 100644 --- a/Src/EB/AMReX_EB2_IndexSpaceI.H +++ b/Src/EB/AMReX_EB2_IndexSpaceI.H @@ -106,3 +106,49 @@ IndexSpaceImp::addFineLevels (int num_new_fine_levels) m_domain.insert(m_domain.begin(), fine_isp.m_domain.begin(), fine_isp.m_domain.end()); m_ngrow.insert(m_ngrow.begin(), fine_isp.m_ngrow.begin(), fine_isp.m_ngrow.end()); } + +template +void +IndexSpaceImp::addRegularCoarseLevels (int num_new_coarse_levels) +{ + if (num_new_coarse_levels <= 0) { return; } + + auto nlevs_old = int(m_gslevel.size()); + int nlevs_new = nlevs_old + num_new_coarse_levels; + + Vector> new_gslevel; + new_gslevel.reserve(nlevs_new); + + Vector new_geom; + new_geom.reserve(nlevs_new); + + Vector new_domain; + new_domain.reserve(nlevs_new); + + Vector new_ngrow; + new_ngrow.reserve(nlevs_new); + + for (int ilev = 0; ilev < num_new_coarse_levels; ++ilev) { + int rr = 1 << (num_new_coarse_levels-ilev); // 2^(num_new_coarse_levels-ilev) + new_geom.push_back(amrex::coarsen(m_geom[0], rr)); + new_domain.push_back(new_geom.back().Domain()); + new_ngrow.push_back(m_ngrow[0]); + new_gslevel.push_back(GShopLevel::makeAllRegular(this, new_geom.back())); + } + + for (int ilev = 0; ilev < nlevs_old; ++ilev) { + new_gslevel.emplace_back(std::move(m_gslevel[ilev])); + new_geom.push_back (m_geom [ilev]); + new_domain.push_back(m_domain[ilev]); + new_ngrow.push_back (m_ngrow [ilev]); + } + + std::swap(new_gslevel, m_gslevel); + std::swap(new_geom , m_geom); + std::swap(new_domain , m_domain); + std::swap(new_ngrow , m_ngrow); + + for (int ilev = num_new_coarse_levels-1; ilev >= 0; --ilev) { + m_gslevel[ilev].buildCutCellMask(m_gslevel[ilev+1]); + } +} diff --git a/Src/EB/AMReX_EB2_IndexSpace_STL.H b/Src/EB/AMReX_EB2_IndexSpace_STL.H index 7ebaaeb7646..0c72d076ea2 100644 --- a/Src/EB/AMReX_EB2_IndexSpace_STL.H +++ b/Src/EB/AMReX_EB2_IndexSpace_STL.H @@ -34,6 +34,7 @@ public: return m_geom.back().Domain(); } void addFineLevels (int num_new_fine_levels) final; + void addRegularCoarseLevels (int num_new_coarse_levels) final; private: diff --git a/Src/EB/AMReX_EB2_IndexSpace_STL.cpp b/Src/EB/AMReX_EB2_IndexSpace_STL.cpp index 335b490ce48..662aaf14dd6 100644 --- a/Src/EB/AMReX_EB2_IndexSpace_STL.cpp +++ b/Src/EB/AMReX_EB2_IndexSpace_STL.cpp @@ -88,4 +88,10 @@ IndexSpaceSTL::addFineLevels (int /*num_new_fine_levels*/) amrex::Abort("IndexSpaceSTL::addFineLevels: todo"); } +void +IndexSpaceSTL::addRegularCoarseLevels (int /*num_new_coarse_levels*/) +{ + amrex::Abort("IndexSpaceSTL::addRegularCoarseLevels: todo"); +} + } diff --git a/Src/EB/AMReX_EB2_IndexSpace_chkpt_file.H b/Src/EB/AMReX_EB2_IndexSpace_chkpt_file.H index 2605118d3b2..cdac3efb7f2 100644 --- a/Src/EB/AMReX_EB2_IndexSpace_chkpt_file.H +++ b/Src/EB/AMReX_EB2_IndexSpace_chkpt_file.H @@ -33,6 +33,7 @@ public: return m_geom.back().Domain(); } void addFineLevels (int num_new_fine_levels) final; + void addRegularCoarseLevels (int num_new_coarse_levels) final; private: diff --git a/Src/EB/AMReX_EB2_IndexSpace_chkpt_file.cpp b/Src/EB/AMReX_EB2_IndexSpace_chkpt_file.cpp index 68688c5aab8..cd811d73688 100644 --- a/Src/EB/AMReX_EB2_IndexSpace_chkpt_file.cpp +++ b/Src/EB/AMReX_EB2_IndexSpace_chkpt_file.cpp @@ -83,4 +83,10 @@ IndexSpaceChkptFile::addFineLevels (int /*num_new_fine_levels*/) amrex::Abort("IndexSpaceChkptFile::addFineLevels: not supported"); } +void +IndexSpaceChkptFile::addRegularCoarseLevels (int /*num_new_coarse_levels*/) +{ + amrex::Abort("IndexSpaceChkptFile::addRegularCoarseLevels: not supported"); +} + } diff --git a/Src/EB/AMReX_EB2_Level.H b/Src/EB/AMReX_EB2_Level.H index d7f9a6f361c..daf55dbd7e7 100644 --- a/Src/EB/AMReX_EB2_Level.H +++ b/Src/EB/AMReX_EB2_Level.H @@ -5,6 +5,7 @@ #include #include #include +#include #include #include #include @@ -26,6 +27,7 @@ namespace amrex::EB2 { class IndexSpace; +template class GShopLevel; class Level { @@ -55,16 +57,20 @@ public: const DistributionMapping& DistributionMap () const noexcept { return m_dmap; } Level (IndexSpace const* is, const Geometry& geom) : m_geom(geom), m_parent(is) {} - void prepareForCoarsening (const Level& rhs, int max_grid_size, IntVect ngrow); + void prepareForCoarsening (const Level& rhs, int max_grid_size, IntVect const& ngrow); const Geometry& Geom () const noexcept { return m_geom; } IndexSpace const* getEBIndexSpace () const noexcept { return m_parent; } void write_to_chkpt_file (const std::string& fname, bool extend_domain_face, int max_grid_size) const; + bool hasEBInfo () const noexcept { return m_has_eb_info; } + void fillCutCellMask (iMultiFab& cutcellmask, const Geometry& geom) const; + // public: // for cuda int coarsenFromFine (Level& fineLevel, bool fill_boundary); void buildCellFlag (); + void buildCutCellMask (Level const& fine_level); protected: @@ -84,9 +90,21 @@ protected: Array m_areafrac; Array m_facecent; Array m_edgecent; + iMultiFab m_cutcellmask; bool m_allregular = false; bool m_ok = false; + bool m_has_eb_info = true; IndexSpace const* m_parent; + +private: + template friend class GShopLevel; + + // Need this function to work around a gcc bug. + void setRegularLevel () { + m_allregular = true; + m_ok = true; + m_has_eb_info = false; + } }; template @@ -101,6 +119,14 @@ public: GShopLevel (IndexSpace const* is, const Geometry& geom); void define_fine (G const& gshop, const Geometry& geom, int max_grid_size, int ngrow, bool extend_domain_face, int num_crse_opt); + + static GShopLevel + makeAllRegular(IndexSpace const* is, const Geometry& geom) + { + GShopLevel r(is, geom); + r.setRegularLevel(); + return r; + } }; template diff --git a/Src/EB/AMReX_EB2_Level.cpp b/Src/EB/AMReX_EB2_Level.cpp index 112bf9db044..106e42d1e43 100644 --- a/Src/EB/AMReX_EB2_Level.cpp +++ b/Src/EB/AMReX_EB2_Level.cpp @@ -7,7 +7,7 @@ namespace amrex::EB2 { void -Level::prepareForCoarsening (const Level& rhs, int max_grid_size, IntVect ngrow) +Level::prepareForCoarsening (const Level& rhs, int max_grid_size, IntVect const& ngrow) { BoxArray all_grids(amrex::grow(m_geom.Domain(),ngrow)); all_grids.maxSize(max_grid_size); @@ -917,6 +917,15 @@ Level::fillLevelSet (MultiFab& levelset, const Geometry& geom) const } } +void +Level::fillCutCellMask (iMultiFab& cutcellmask, const Geometry&) const +{ + if (!m_has_eb_info) { + cutcellmask.setVal(0); + cutcellmask.ParallelCopy(m_cutcellmask); + } +} + void Level::write_to_chkpt_file (const std::string& fname, bool extend_domain_face, int max_grid_size) const { @@ -927,4 +936,86 @@ Level::write_to_chkpt_file (const std::string& fname, bool extend_domain_face, i m_geom, m_ngrow, extend_domain_face, max_grid_size); } +void +Level::buildCutCellMask (Level const& fine_level) +{ + AMREX_ALWAYS_ASSERT(!m_has_eb_info); + + MFInfo mf_info; + mf_info.SetTag("EB2::Level"); + + m_dmap = fine_level.m_dmap; + const BoxArray& fine_grids = fine_level.m_grids; + if (fine_level.hasEBInfo()) + { + AMREX_ALWAYS_ASSERT(fine_grids.coarsenable(2)); + m_grids = amrex::coarsen(fine_grids,2); + m_cutcellmask.define(m_grids,m_dmap,1,0,mf_info); + + auto const& farrs = fine_level.m_cellflag.const_arrays(); + auto const& carrs = m_cutcellmask.arrays(); + ParallelFor(m_cutcellmask, + [=] AMREX_GPU_DEVICE (int bno, int i, int j, int k) + { + auto const& fa = farrs[bno]; + auto const& ca = carrs[bno]; +#if (AMREX_SPACEDIM == 3) + int k3d = 1; +#else + int k3d = 0; +#endif + bool cut = false; + for (int kk = k*2; kk <= k*2+k3d; ++kk) { + for (int jj = j*2; jj <= j*2+1 ; ++jj) { + for (int ii = i*2; ii <= i*2+1 ; ++ii) { + if (fa(ii,jj,kk).isSingleValued()) { cut = true; } + }}} + carrs[bno](i,j,k) = int(cut); + }); + Gpu::streamSynchronize(); + } + else + { + iMultiFab raii; + iMultiFab const* fine_mask; + if (fine_grids.coarsenable(2)) + { + fine_mask = &(fine_level.m_cutcellmask); + } + else + { + BoxList bl; + bl.reserve(fine_grids.size()); + for (Long ibox = 0; ibox < fine_grids.size(); ++ibox) { + bl.push_back(amrex::refine(amrex::coarsen(fine_grids[ibox],8),8)); + } + raii.define(BoxArray(std::move(bl)), fine_level.m_dmap, 1, 0); + raii.setVal(0); + raii.ParallelCopy(fine_level.m_cutcellmask); + fine_mask = &raii; + } + + m_grids = amrex::coarsen(fine_mask->boxArray(),2); + m_cutcellmask.define(m_grids,m_dmap,1,0,mf_info); + + auto const& farrs = fine_mask->const_arrays(); + auto const& carrs = m_cutcellmask.arrays(); + ParallelFor(m_cutcellmask, + [=] AMREX_GPU_DEVICE (int bno, int i, int j, int k) + { + auto const& fa = farrs[bno]; + auto const& ca = carrs[bno]; + constexpr int k3d = (AMREX_SPACEDIM == 3) ? 1 : 0; + bool cut = false; + for (int kk = k*2; kk <= k*2+k3d; ++kk) { + for (int jj = j*2; jj <= j*2+1 ; ++jj) { + for (int ii = i*2; ii <= i*2+1 ; ++ii) { + if (fa(ii,jj,kk)) { cut = true; } + }}} + carrs[bno](i,j,k) = int(cut); + }); + Gpu::streamSynchronize(); + } +} + } diff --git a/Src/EB/AMReX_EBAmrUtil.cpp b/Src/EB/AMReX_EBAmrUtil.cpp index 9bad74a4acc..3e756d40648 100644 --- a/Src/EB/AMReX_EBAmrUtil.cpp +++ b/Src/EB/AMReX_EBAmrUtil.cpp @@ -2,6 +2,7 @@ #include #include #include +#include #ifdef AMREX_USE_OMP #include @@ -16,28 +17,44 @@ TagCutCells (TagBoxArray& tags, const MultiFab& state) // const char clearval = TagBox::CLEAR; auto const& factory = dynamic_cast(state.Factory()); - auto const& flags = factory.getMultiEBCellFlagFab(); + + if (factory.hasEBInfo()) { + auto const& flags = factory.getMultiEBCellFlagFab(); #ifdef AMREX_USE_OMP #pragma omp parallel if (Gpu::notInLaunchRegion()) #endif - for (MFIter mfi(state, TilingIfNotGPU()); mfi.isValid(); ++mfi) - { - const Box& bx = mfi.tilebox(); + for (MFIter mfi(state, TilingIfNotGPU()); mfi.isValid(); ++mfi) + { + const Box& bx = mfi.tilebox(); - const auto& flag = flags[mfi]; + const auto& flag = flags[mfi]; - const FabType typ = flag.getType(bx); - if (typ != FabType::regular && typ != FabType::covered) - { - Array4 const& tagarr = tags.array(mfi); - Array4 const& flagarr = flag.const_array(); - AMREX_HOST_DEVICE_FOR_3D (bx, i, j, k, + const FabType typ = flag.getType(bx); + if (typ != FabType::regular && typ != FabType::covered) + { + Array4 const& tagarr = tags.array(mfi); + Array4 const& flagarr = flag.const_array(); + AMREX_HOST_DEVICE_FOR_3D (bx, i, j, k, + { + if (flagarr(i,j,k).isSingleValued()) { + tagarr(i,j,k) = tagval; + } + }); + } + } + } else { + auto* cutcell_mask = factory.getCutCellMask(); + if (cutcell_mask) { + auto const& ta = tags.arrays(); + auto const& ma = cutcell_mask->const_arrays(); + ParallelFor(state, [=] AMREX_GPU_DEVICE (int bno, int i, int j, int k) { - if (flagarr(i,j,k).isSingleValued()) { - tagarr(i,j,k) = tagval; + if (ma[bno](i,j,k)) { + ta[bno](i,j,k) = tagval; } }); + amrex::Gpu::streamSynchronize(); } } } diff --git a/Src/EB/AMReX_EBDataCollection.H b/Src/EB/AMReX_EBDataCollection.H index 13edcdb7243..17edaabef16 100644 --- a/Src/EB/AMReX_EBDataCollection.H +++ b/Src/EB/AMReX_EBDataCollection.H @@ -11,6 +11,7 @@ namespace amrex { template class FabArray; class MultiFab; +class iMultiFab; class MultiCutFab; namespace EB2 { class Level; } @@ -39,6 +40,7 @@ public: [[nodiscard]] Array getAreaFrac () const; [[nodiscard]] Array getFaceCent () const; [[nodiscard]] Array getEdgeCent () const; + [[nodiscard]] const iMultiFab* getCutCellMask () const; private: @@ -63,6 +65,9 @@ private: Array m_areafrac {{AMREX_D_DECL(nullptr, nullptr, nullptr)}}; Array m_facecent {{AMREX_D_DECL(nullptr, nullptr, nullptr)}}; Array m_edgecent {{AMREX_D_DECL(nullptr, nullptr, nullptr)}}; + + // for levels created by addRegularCoarseLevels only + iMultiFab* m_cutcellmask = nullptr; }; } diff --git a/Src/EB/AMReX_EBDataCollection.cpp b/Src/EB/AMReX_EBDataCollection.cpp index d38efc4af7f..78728c76907 100644 --- a/Src/EB/AMReX_EBDataCollection.cpp +++ b/Src/EB/AMReX_EBDataCollection.cpp @@ -1,6 +1,7 @@ #include #include +#include #include #include @@ -65,6 +66,12 @@ EBDataCollection::EBDataCollection (const EB2::Level& a_level, a_level.fillFaceCent(m_facecent, m_geom); a_level.fillEdgeCent(m_edgecent, m_geom); } + + if (! a_level.hasEBInfo()) { + m_cutcellmask = new iMultiFab(a_ba, a_dm, 1, 0, MFInfo(), + DefaultFabFactory()); + a_level.fillCutCellMask(*m_cutcellmask, m_geom); + } } EBDataCollection::~EBDataCollection () @@ -81,6 +88,7 @@ EBDataCollection::~EBDataCollection () delete m_facecent[idim]; delete m_edgecent[idim]; } + delete m_cutcellmask; } const FabArray& @@ -153,4 +161,10 @@ EBDataCollection::getBndryNormal () const return *m_bndrynorm; } +const iMultiFab* +EBDataCollection::getCutCellMask () const +{ + return m_cutcellmask; +} + } diff --git a/Src/EB/AMReX_EBFabFactory.H b/Src/EB/AMReX_EBFabFactory.H index aa8a3ece2d2..c9db5f39473 100644 --- a/Src/EB/AMReX_EBFabFactory.H +++ b/Src/EB/AMReX_EBFabFactory.H @@ -83,6 +83,12 @@ public: [[nodiscard]] const BoxArray& boxArray () const noexcept; [[nodiscard]] const Geometry& Geom () const noexcept { return m_geom; } + [[nodiscard]] bool hasEBInfo() const noexcept; + + //! Returns nullptr unless this level is built by EB2::addRegularCoarseLevels. + //! One should use getMultiEBCellFlagFab for normal levels. + [[nodiscard]] iMultiFab const* getCutCellMask () const noexcept { return m_ebdc->getCutCellMask(); } + private: EBSupport m_support; diff --git a/Src/EB/AMReX_EBFabFactory.cpp b/Src/EB/AMReX_EBFabFactory.cpp index 1fd928b3c0b..f5154578d6f 100644 --- a/Src/EB/AMReX_EBFabFactory.cpp +++ b/Src/EB/AMReX_EBFabFactory.cpp @@ -109,6 +109,12 @@ EBFArrayBoxFactory::boxArray () const noexcept return m_ebdc->getMultiEBCellFlagFab().boxArray(); } +bool +EBFArrayBoxFactory::hasEBInfo () const noexcept +{ + return m_parent->hasEBInfo(); +} + std::unique_ptr makeEBFabFactory (const Geometry& a_geom, const BoxArray& a_ba,