diff --git a/Src/EB/AMReX_EB2_Level.H b/Src/EB/AMReX_EB2_Level.H index daf55dbd7e7..6c6f75f784b 100644 --- a/Src/EB/AMReX_EB2_Level.H +++ b/Src/EB/AMReX_EB2_Level.H @@ -62,6 +62,8 @@ public: const Geometry& Geom () const noexcept { return m_geom; } IndexSpace const* getEBIndexSpace () const noexcept { return m_parent; } + IntVect const& nGrowVect () const noexcept { return m_ngrow; } + 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; } diff --git a/Src/EB/AMReX_EBDataCollection.H b/Src/EB/AMReX_EBDataCollection.H index 17edaabef16..6a8556dc9a5 100644 --- a/Src/EB/AMReX_EBDataCollection.H +++ b/Src/EB/AMReX_EBDataCollection.H @@ -42,6 +42,9 @@ public: [[nodiscard]] Array getEdgeCent () const; [[nodiscard]] const iMultiFab* getCutCellMask () const; + // public for cuda + void extendDataOutsideDomain (IntVect const& level_ng); + private: Vector m_ngrow; diff --git a/Src/EB/AMReX_EBDataCollection.cpp b/Src/EB/AMReX_EBDataCollection.cpp index 78728c76907..3a232ab5ab9 100644 --- a/Src/EB/AMReX_EBDataCollection.cpp +++ b/Src/EB/AMReX_EBDataCollection.cpp @@ -5,6 +5,7 @@ #include #include +#include #include namespace amrex { @@ -13,7 +14,7 @@ EBDataCollection::EBDataCollection (const EB2::Level& a_level, const Geometry& a_geom, const BoxArray& a_ba_in, const DistributionMapping& a_dm, - Vector a_ngrow, EBSupport a_support) + Vector a_ngrow, EBSupport a_support) : m_ngrow(std::move(a_ngrow)), m_support(a_support), m_geom(a_geom) @@ -33,6 +34,8 @@ EBDataCollection::EBDataCollection (const EB2::Level& a_level, if (m_support >= EBSupport::volume) { + AMREX_ALWAYS_ASSERT(m_ngrow[1] <= m_ngrow[0]); + m_volfrac = new MultiFab(a_ba, a_dm, 1, m_ngrow[1], MFInfo(), FArrayBoxFactory()); a_level.fillVolFrac(*m_volfrac, m_geom); @@ -42,6 +45,8 @@ EBDataCollection::EBDataCollection (const EB2::Level& a_level, if (m_support == EBSupport::full) { + AMREX_ALWAYS_ASSERT(m_ngrow[2] <= m_ngrow[0]); + const int ng = m_ngrow[2]; m_bndrycent = new MultiCutFab(a_ba, a_dm, AMREX_SPACEDIM, ng, *m_cellflags); @@ -72,6 +77,147 @@ EBDataCollection::EBDataCollection (const EB2::Level& a_level, DefaultFabFactory()); a_level.fillCutCellMask(*m_cutcellmask, m_geom); } + + extendDataOutsideDomain(a_level.nGrowVect()); +} + +void EBDataCollection::extendDataOutsideDomain (IntVect const& level_ng) +{ + if (m_cellflags == nullptr) { return; } + + int const ngrow = m_ngrow[0]; + Box const& data_domain = amrex::grow(m_geom.Domain(), ngrow); + + Box level_domain = m_geom.Domain(); + for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { + if (m_geom.isPeriodic(idim)) { + level_domain.grow(idim, ngrow); + } else { + level_domain.grow(idim, level_ng[idim]); + } + } + + if (level_domain.contains(data_domain)) { return; } + + Box const& level_nodal_domain = amrex::surroundingNodes(level_domain); + Array lev_ap_domain + {AMREX_D_DECL(amrex::surroundingNodes(level_domain,0), + amrex::surroundingNodes(level_domain,1), + amrex::surroundingNodes(level_domain,2))}; + +#ifdef AMREX_USE_OMP +#pragma omp parallel if (Gpu::notInLaunchRegion()) +#endif + for (MFIter mfi(*m_cellflags); mfi.isValid(); ++mfi) { + Box const& bx = mfi.fabbox(); + if (! level_domain.contains(bx)) { + Box const& nbx = amrex::surroundingNodes(bx); + auto const& ls_a = m_levelset->array(mfi); + auto const& flag_a = m_cellflags->array(mfi); + Array4 vfrc_a; + if (m_volfrac) { + vfrc_a = m_volfrac->array(mfi); + } + amrex::ParallelFor(nbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + if (! level_nodal_domain.contains(i,j,k)) { + int ii = std::clamp(i, level_nodal_domain.smallEnd(0), + level_nodal_domain.bigEnd(0)); + int jj = std::clamp(j, level_nodal_domain.smallEnd(1), + level_nodal_domain.bigEnd(1)); +#if (AMREX_SPACEDIM > 2) + int kk = std::clamp(k, level_nodal_domain.smallEnd(2), + level_nodal_domain.bigEnd(2)); +#else + int kk = 0; +#endif + ls_a(i,j,k) = ls_a(ii,jj,kk); + } + }); + amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + if (! level_domain.contains(i,j,k)) { + EBCellFlag flag; + int ii = std::clamp(i, level_domain.smallEnd(0), + level_domain.bigEnd(0)); + int jj = std::clamp(j, level_domain.smallEnd(1), + level_domain.bigEnd(1)); +#if (AMREX_SPACEDIM > 2) + int kk = std::clamp(k, level_domain.smallEnd(2), + level_domain.bigEnd(2)); +#else + int kk = 0; +#endif + if (flag_a(ii,jj,kk).isCovered()) { + flag.setCovered(); + } else if (flag_a(ii,jj,kk).isRegular()) { + flag.setRegular(); + } else { + int ncov = 0; +#if (AMREX_SPACEDIM > 2) + for (int ko = k; ko <= k+1; ++ko) { +#else + for (int ko = k; ko <= k ; ++ko) { +#endif + for (int jo = j; jo <= j+1; ++jo) { + for (int io = i; io <= i+1; ++io) { + if (ls_a(io,jo,ko) >= Real(0.0)) { ++ncov; } + }}} + constexpr int nnodes = (AMREX_SPACEDIM == 2) ? 4 : 8; + if (ncov == nnodes) { + flag.setCovered(); + } else if (ncov == 0) { + flag.setRegular(); + } else { + flag.setSingleValued(); + } + } + if (flag.isCovered()) { + flag_a(i,j,k).setCovered(); + if (vfrc_a && vfrc_a.contains(i,j,k)) { + vfrc_a(i,j,k) = Real(0.0); + } + } else if (flag.isRegular()) { + flag_a(i,j,k).setRegular(); + if (vfrc_a && vfrc_a.contains(i,j,k)) { + vfrc_a(i,j,k) = Real(1.0); + } + } else { + flag_a(i,j,k).setSingleValued(); + if (vfrc_a && vfrc_a.contains(i,j,k)) { + vfrc_a(i,j,k) = vfrc_a(ii,jj,kk); + } + } + } + }); + + for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { + if (m_areafrac[idim]) { + auto const& ap = m_areafrac[idim]->array(mfi); + Box apbx = Box(ap); + if (apbx.smallEnd(idim) == nbx.smallEnd(idim)) { + apbx.growLo(idim,-1); + } + if (apbx.bigEnd(idim) == nbx.bigEnd(idim)) { + apbx.growHi(idim,-1); + } + auto lev_apidim_domain = lev_ap_domain[idim]; + Dim3 off = IntVect::TheDimensionVector(idim).dim3(); + amrex::ParallelFor(apbx, + [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + if (! lev_apidim_domain.contains(i,j,k)) { + if (flag_a(i-off.x,j-off.y,k-off.z).isCovered() || + flag_a(i,j,k).isCovered()) + { + ap(i,j,k) = Real(0.0); + } + } + }); + } + } + } + } } EBDataCollection::~EBDataCollection ()