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<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
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<IndexSpace*>(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<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]);
+    }
+}
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 <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>
diff --git a/Src/EB/AMReX_EB2_Level.cpp b/Src/EB/AMReX_EB2_Level.cpp
index 112bf9db044..9df1423f6f1 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,84 @@ 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;
+            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];
+            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..fd0f80f419d 100644
--- a/Src/EB/AMReX_EBAmrUtil.cpp
+++ b/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();
         }
     }
 }
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 T> class FabArray;
 class MultiFab;
+class iMultiFab;
 class MultiCutFab;
 namespace EB2 { class Level; }
 
@@ -39,6 +40,7 @@ public:
     [[nodiscard]] Array<const MultiCutFab*, AMREX_SPACEDIM> getAreaFrac () const;
     [[nodiscard]] Array<const MultiCutFab*, AMREX_SPACEDIM> getFaceCent () const;
     [[nodiscard]] Array<const MultiCutFab*, AMREX_SPACEDIM> getEdgeCent () const;
+    [[nodiscard]] const iMultiFab* getCutCellMask () const;
 
 private:
 
@@ -63,6 +65,9 @@ private:
     Array<MultiCutFab*,AMREX_SPACEDIM> m_areafrac {{AMREX_D_DECL(nullptr, nullptr, nullptr)}};
     Array<MultiCutFab*,AMREX_SPACEDIM> m_facecent {{AMREX_D_DECL(nullptr, nullptr, nullptr)}};
     Array<MultiCutFab*,AMREX_SPACEDIM> 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 <AMReX_EBDataCollection.H>
 #include <AMReX_MultiFab.H>
+#include <AMReX_iMultiFab.H>
 #include <AMReX_MultiCutFab.H>
 
 #include <AMReX_EB2_Level.H>
@@ -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<IArrayBox>());
+        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<EBCellFlagFab>&
@@ -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<EBFArrayBoxFactory>
 makeEBFabFactory (const Geometry& a_geom,
                   const BoxArray& a_ba,