New feature in EB: Add regular coarse levels (#3414)

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

@@ -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<std::unique_ptr<IndexSpace> > 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

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<IndexSpace*>(TopIndexSpace());
+ if (p) {
+ p->addRegularCoarseLevels(num_new_coarse_levels);
+ }
+}
+
 void
 BuildFromChkptFile (std::string const& fname,
  const Geometry& geom, int required_coarsening_level,

Src/EB/AMReX_EB2_IndexSpaceI.H

@@ -106,3 +106,49 @@ IndexSpaceImp<G>::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 <typename G>
+void
+IndexSpaceImp<G>::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<GShopLevel<G>> new_gslevel;
+ new_gslevel.reserve(nlevs_new);
+
+ Vector<Geometry> new_geom;
+ new_geom.reserve(nlevs_new);
+
+ Vector<Box> new_domain;
+ new_domain.reserve(nlevs_new);
+
+ Vector<int> 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<G>::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]);
+ }
+}

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:
 

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");
+}
+
 }

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:
 

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");
+}
+
 }

Src/EB/AMReX_EB2_Level.H

@@ -5,6 +5,7 @@
 #include <AMReX_ParmParse.H>
 #include <AMReX_Geometry.H>
 #include <AMReX_MultiFab.H>
+#include <AMReX_iMultiFab.H>
 #include <AMReX_LayoutData.H>
 #include <AMReX_VisMF.H>
 #include <AMReX_Array.H>
@@ -26,6 +27,7 @@
 namespace amrex::EB2 {
 
 class IndexSpace;
+template <typename G> 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<MultiFab,AMREX_SPACEDIM> m_areafrac;
  Array<MultiFab,AMREX_SPACEDIM> m_facecent;
  Array<MultiFab,AMREX_SPACEDIM> 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 <typename G> 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 <typename G>
@@ -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<G>
+ makeAllRegular(IndexSpace const* is, const Geometry& geom)
+ {
+ GShopLevel<G> r(is, geom);
+ r.setRegularLevel();
+ return r;
+ }
 };
 
 template <typename G>

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,85 @@ 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];
+#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;
+ int nboxes = int(fine_grids.size());
+ bl.reserve(nboxes);
+ for (int ibox = 0; ibox < nboxes; ++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];
+ 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();
+ }
+}
+
 }

Src/EB/AMReX_EBAmrUtil.cpp

@@ -2,6 +2,7 @@
 #include <AMReX_EBAmrUtil.H>
 #include <AMReX_EBFArrayBox.H>
 #include <AMReX_EBCellFlag.H>
+#include <AMReX_iMultiFab.H>
 
 #ifdef AMREX_USE_OMP
 #include <omp.h>
@@ -16,28 +17,44 @@ TagCutCells (TagBoxArray& tags, const MultiFab& state)
 // const char clearval = TagBox::CLEAR;
 
  auto const& factory = dynamic_cast<EBFArrayBoxFactory const&>(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<char> const& tagarr = tags.array(mfi);
- Array4<EBCellFlag const> 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<char> const& tagarr = tags.array(mfi);
+ Array4<EBCellFlag const> 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 const* 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();
  }
  }
 }