Control EB generation level

Source/PeleLM.H

@@ -742,6 +742,8 @@ class PeleLM : public amrex::AmrCore {
 
    void makeEBGeometry();
 
+   int getRestartEBMaxLevel();
+
    void initialRedistribution();
 
    void initCoveredState();
@@ -1232,6 +1234,7 @@ class PeleLM : public amrex::AmrCore {
    int m_EB_refine_LevMin = -1;
    int m_EB_refine_LevAdapt = -1;
    int m_signDistNeeded = 0;
+   int m_EB_generate_max_level = -1;
    amrex::Real m_derefineEBBuffer = 3.0;
    std::unique_ptr<amrex::MultiFab> m_signedDist0 = nullptr;
    amrex::Vector<amrex::Real> coveredState_h;

Source/PeleLMEB.cpp

@@ -21,12 +21,32 @@ void PeleLM::makeEBGeometry()
     std::string geom_type;
     ppeb2.get("geom_type", geom_type);
 
+    // At what level should the EB be generated ?
+    // Default : max_level
+    int prev_max_lvl_eb = max_level;
+
+    // If restarting: use the level at which it was generated earlier
+    prev_max_lvl_eb  = (getRestartEBMaxLevel() > 0) ? getRestartEBMaxLevel()
+                                                    : prev_max_lvl_eb;
+
+    // Manual override if needed -> might incur issues with previously covered/uncovered
+    // cell flipping
+    ppeb2.query("max_level_generation",prev_max_lvl_eb);
+    AMREX_ALWAYS_ASSERT(prev_max_lvl_eb <= max_level);
+
+    // Stash away the level at which we ended up
+    m_EB_generate_max_level = prev_max_lvl_eb;
+
+    // Generate the EB data at prev_max_lvl_eb and create consistent coarse version from there.
     if (geom_type == "UserDefined") {
-        EBUserDefined(geom.back(), req_coarsening_level, max_coarsening_level);
+        EBUserDefined(geom[prev_max_lvl_eb], req_coarsening_level, max_coarsening_level);
     } else {
         // If geom_type is not an AMReX recognized type, it'll crash.
-        EB2::Build(geom.back(), req_coarsening_level, max_coarsening_level);
+        EB2::Build(geom[prev_max_lvl_eb], req_coarsening_level, max_coarsening_level);
     }
+
+    // Add finer level, might be inconsistent with the coarser level created above.
+    EB2::addFineLevels(max_level - prev_max_lvl_eb);
 }
 
 void PeleLM::redistributeAofS(int a_lev,
@@ -432,4 +452,44 @@ PeleLM::correct_vel_small_cells (Vector<MultiFab*> const& a_vel,
     }
 }
 
+int
+PeleLM::getRestartEBMaxLevel()
+{
+    if (m_restart_chkfile.empty()) {
+        return -1;
+    } else {
+        // Go and parse the line we need
+        std::string File(m_restart_chkfile + "/Header");
+
+        VisMF::IO_Buffer io_buffer(VisMF::GetIOBufferSize());
+
+        Vector<char> fileCharPtr;
+        ParallelDescriptor::ReadAndBcastFile(File, fileCharPtr);
+        std::string fileCharPtrString(fileCharPtr.dataPtr());
+        std::istringstream is(fileCharPtrString, std::istringstream::in);
+
+        std::string line, word;
+
+        // Title line
+        std::getline(is, line);
+
+        // Finest level
+        std::getline(is, line);
+
+        // Step count
+        std::getline(is, line);
+
+        // Either what we're looking for or m_cur_time
+        std::getline(is, line);
+
+        if (line.find('.') != std::string::npos) {
+            return -1;
+        } else {
+            int max_eb_rst = -1;
+            max_eb_rst = std::stoi(line);
+            return max_eb_rst;
+        }
+    }
+}
+
 #endif

Source/PeleLMPlot.cpp

@@ -338,9 +338,13 @@ void PeleLM::WriteHeader(const std::string& name, bool is_checkpoint) const
         }
 
         HeaderFile << finest_level << "\n";
-
-        // Time stepping controls
+
         HeaderFile << m_nstep << "\n";
+
+#ifdef AMREX_USE_EB
+        HeaderFile << m_EB_generate_max_level << "\n";
+#endif
+
         HeaderFile << m_cur_time << "\n";
         HeaderFile << m_dt << "\n";
         HeaderFile << m_prev_dt << "\n";
@@ -462,9 +466,26 @@ void PeleLM::ReadCheckPointFile()
    is >> m_nstep;
    GotoNextLine(is);
 
+#ifdef AMREX_USE_EB
+   // Finest level at which EB was generated
+   // actually used independently, so just skip ...
+   std::getline(is, line);
+
+   // ... but to be backward compatible, if we get a float,
+   // let's assume it's m_cur_time
+   if (line.find('.') != std::string::npos) {
+      m_cur_time = std::stod(line);
+   } else {
+      // Skip line and read current time
+      is >> m_cur_time;
+      GotoNextLine(is);
+   }
+#else 
+
    // Current time
    is >> m_cur_time;
    GotoNextLine(is);
+#endif
 
    // Time step size
    is >> m_dt;

Source/PeleLMSetup.cpp

@@ -124,7 +124,6 @@ void PeleLM::Setup() {
    prob_parm = new ProbParm{};
    prob_parm_d = (ProbParm*)The_Arena()->alloc(sizeof(ProbParm));
 
-
    // Problem parameters
    readProbParm();