Merge branch 'development' into fix_interp_for_ratio1

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();
         }
     }
 }