From 103db6ebe2b570910ac4dbd7d6611e59d80f1a0b Mon Sep 17 00:00:00 2001
From: Weiqun Zhang <WeiqunZhang@lbl.gov>
Date: Fri, 5 Aug 2022 15:25:33 -0700
Subject: [PATCH] EB: Add Fine Levels (#2881)

Add a new function EB2::addFineLevels() that can be used to add more fine
levels to the existing EB IndexSpace without changing the coarse levels.
This is useful for restarting with a larger amr.max_level.  The issue is we
build EB at the finest level first and then coarsen it to the coarse levels.
If the restart run has a different finest level, the EB on the coarse levels
could be different without using this new capability.
---
 Src/EB/AMReX_EB2.H                  | 10 ++++++++-
 Src/EB/AMReX_EB2.cpp                |  9 ++++++++
 Src/EB/AMReX_EB2_IndexSpaceI.H      | 32 +++++++++++++++++++++++++++--
 Src/EB/AMReX_EB2_IndexSpace_STL.H   |  1 +
 Src/EB/AMReX_EB2_IndexSpace_STL.cpp |  6 ++++++
 Tests/EB/CNS/Source/main.cpp        |  8 +++++++-
 6 files changed, 62 insertions(+), 4 deletions(-)

diff --git a/Src/EB/AMReX_EB2.H b/Src/EB/AMReX_EB2.H
index ff897276510..6a143bf2a9c 100644
--- a/Src/EB/AMReX_EB2.H
+++ b/Src/EB/AMReX_EB2.H
@@ -49,6 +49,7 @@ public:
     virtual const Level& getLevel (const Geometry & geom) const = 0;
     virtual const Geometry& getGeometry (const Box& domain) const = 0;
     virtual const Box& coarsestDomain () const = 0;
+    virtual void addFineLevels (int num_new_fine_levels) = 0;
 
 protected:
     static AMREX_EXPORT Vector<std::unique_ptr<IndexSpace> > m_instance;
@@ -80,16 +81,21 @@ public:
     virtual const Box& coarsestDomain () const final {
         return m_geom.back().Domain();
     }
+    virtual void addFineLevels (int num_new_fine_levels) final;
 
     using F = typename G::FunctionType;
 
 private:
 
+    G m_gshop;
+    bool m_build_coarse_level_by_coarsening;
+    bool m_extend_domain_face;
+    int m_num_coarsen_opt;
+
     Vector<GShopLevel<G> > m_gslevel;
     Vector<Geometry> m_geom;
     Vector<Box> m_domain;
     Vector<int> m_ngrow;
-    std::unique_ptr<F> m_impfunc;
 };
 
 #include <AMReX_EB2_IndexSpaceI.H>
@@ -125,6 +131,8 @@ void Build (const Geometry& geom,
 int maxCoarseningLevel (const Geometry& geom);
 int maxCoarseningLevel (IndexSpace const* ebis, const Geometry& geom);
 
+void addFineLevels (int num_new_fine_levels);
+
 }}
 
 #endif
diff --git a/Src/EB/AMReX_EB2.cpp b/Src/EB/AMReX_EB2.cpp
index 16f683cb019..4f2ad5bf873 100644
--- a/Src/EB/AMReX_EB2.cpp
+++ b/Src/EB/AMReX_EB2.cpp
@@ -230,6 +230,15 @@ Build (const Geometry& geom, int required_coarsening_level,
     }
 }
 
+void addFineLevels (int num_new_fine_levels)
+{
+    BL_PROFILE("EB2::addFineLevels()");
+    auto p = const_cast<IndexSpace*>(TopIndexSpace());
+    if (p) {
+        p->addFineLevels(num_new_fine_levels);
+    }
+}
+
 namespace {
 static int comp_max_crse_level (Box cdomain, const Box& domain)
 {
diff --git a/Src/EB/AMReX_EB2_IndexSpaceI.H b/Src/EB/AMReX_EB2_IndexSpaceI.H
index cdad6b31f6f..e7db810b03b 100644
--- a/Src/EB/AMReX_EB2_IndexSpaceI.H
+++ b/Src/EB/AMReX_EB2_IndexSpaceI.H
@@ -5,6 +5,10 @@ IndexSpaceImp<G>::IndexSpaceImp (const G& gshop, const Geometry& geom,
                                  int max_coarsening_level,
                                  int ngrow, bool build_coarse_level_by_coarsening,
                                  bool extend_domain_face, int num_coarsen_opt)
+    : m_gshop(gshop),
+      m_build_coarse_level_by_coarsening(build_coarse_level_by_coarsening),
+      m_extend_domain_face(extend_domain_face),
+      m_num_coarsen_opt(num_coarsen_opt)
 {
     // build finest level (i.e., level 0) first
     AMREX_ALWAYS_ASSERT(required_coarsening_level >= 0 && required_coarsening_level <= 30);
@@ -56,8 +60,6 @@ IndexSpaceImp<G>::IndexSpaceImp (const G& gshop, const Geometry& geom,
         m_domain.push_back(cdomain);
         m_ngrow.push_back(ng);
     }
-
-    m_impfunc = std::make_unique<F>(gshop.GetImpFunc());
 }
 
 
@@ -78,3 +80,29 @@ IndexSpaceImp<G>::getGeometry (const Box& dom) const
     int i = std::distance(m_domain.begin(), it);
     return m_geom[i];
 }
+
+template <typename G>
+void
+IndexSpaceImp<G>::addFineLevels (int num_new_fine_levels)
+{
+    if (num_new_fine_levels <= 0) { return; }
+
+    if (m_num_coarsen_opt > 0) {
+        m_num_coarsen_opt += num_new_fine_levels;
+    }
+
+    IndexSpaceImp<G> fine_isp(m_gshop, amrex::refine(m_geom[0], 1<<num_new_fine_levels),
+                              num_new_fine_levels-1, num_new_fine_levels-1,
+                              m_ngrow[0], m_build_coarse_level_by_coarsening,
+                              m_extend_domain_face, m_num_coarsen_opt);
+
+    fine_isp.m_gslevel.reserve(m_domain.size()+num_new_fine_levels);
+    for (int i = 0; i < m_domain.size(); ++i) {
+        fine_isp.m_gslevel.emplace_back(std::move(m_gslevel[i]));
+    }
+    std::swap(fine_isp.m_gslevel, m_gslevel);
+
+    m_geom.insert(m_geom.begin(), fine_isp.m_geom.begin(), fine_isp.m_geom.end());
+    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());
+}
diff --git a/Src/EB/AMReX_EB2_IndexSpace_STL.H b/Src/EB/AMReX_EB2_IndexSpace_STL.H
index 4c4bb61ca35..a5e13ba0591 100644
--- a/Src/EB/AMReX_EB2_IndexSpace_STL.H
+++ b/Src/EB/AMReX_EB2_IndexSpace_STL.H
@@ -33,6 +33,7 @@ public:
     virtual const Box& coarsestDomain () const final {
         return m_geom.back().Domain();
     }
+    virtual void addFineLevels (int num_new_fine_levels) final;
 
 private:
 
diff --git a/Src/EB/AMReX_EB2_IndexSpace_STL.cpp b/Src/EB/AMReX_EB2_IndexSpace_STL.cpp
index 2c2c77e4c0e..1a81cfe72e9 100644
--- a/Src/EB/AMReX_EB2_IndexSpace_STL.cpp
+++ b/Src/EB/AMReX_EB2_IndexSpace_STL.cpp
@@ -82,4 +82,10 @@ IndexSpaceSTL::getGeometry (const Box& dom) const
     return m_geom[i];
 }
 
+void
+IndexSpaceSTL::addFineLevels (int /*num_new_fine_levels*/)
+{
+    amrex::Abort("IndexSpaceSTL::addFineLevels: todo");
+}
+
 }}
diff --git a/Tests/EB/CNS/Source/main.cpp b/Tests/EB/CNS/Source/main.cpp
index aa851c47956..393431f8f79 100644
--- a/Tests/EB/CNS/Source/main.cpp
+++ b/Tests/EB/CNS/Source/main.cpp
@@ -53,7 +53,13 @@ int main (int argc, char* argv[])
         AmrLevel::SetEBSupportLevel(EBSupport::full);
         AmrLevel::SetEBMaxGrowCells(CNS::numGrow(),4,2);
 
-        initialize_EB2(amr.Geom(amr.maxLevel()), amr.maxLevel(), amr.maxLevel());
+        int max_eb_level = amr.maxLevel();
+        ParmParse pp("amr");
+        pp.query("max_eb_level", max_eb_level);
+        initialize_EB2(amr.Geom(max_eb_level), max_eb_level, max_eb_level);
+        if (max_eb_level < amr.maxLevel()) {
+            EB2::addFineLevels(amr.maxLevel() - max_eb_level);
+        }
 
         amr.init(strt_time,stop_time);