enable incompressible to work with multilevel but only OneWay coupling (

#1612) * enable incompressible to work with multilevel but only OneWay coupling * code cleanup in ERF_slow_rhs_pre and ERF_slow_rhs_inc * add file to cmake * more changes to sync inc with pre --------- Co-authored-by: Aaron M. Lattanzi <[email protected]>
erf-model · May 10, 2024 · 084af37 · 084af37
1 parent 17472ad
commit 084af37
Show file tree

Hide file tree

Showing 17 changed files with 777 additions and 948 deletions.
diff --git a/CMake/BuildERFExe.cmake b/CMake/BuildERFExe.cmake
@@ -165,6 +165,7 @@ function(build_erf_lib erf_lib_name)
        ${SRC_DIR}/TimeIntegration/ERF_advance_lsm.cpp
        ${SRC_DIR}/TimeIntegration/ERF_advance_radiation.cpp
        ${SRC_DIR}/TimeIntegration/ERF_make_fast_coeffs.cpp
+       ${SRC_DIR}/TimeIntegration/ERF_make_tau_terms.cpp
        ${SRC_DIR}/TimeIntegration/ERF_slow_rhs_pre.cpp
        ${SRC_DIR}/TimeIntegration/ERF_slow_rhs_post.cpp
        ${SRC_DIR}/TimeIntegration/ERF_fast_rhs_N.cpp

diff --git a/Exec/Make.ERF b/Exec/Make.ERF
@@ -177,6 +177,7 @@ endif
 
 ifeq ($(USE_POISSON_SOLVE), TRUE)
   DEFINES += -DERF_USE_POISSON_SOLVE
+  USERSuffix += .INC
 endif
 
 ifeq ($(USE_PARTICLES), TRUE)

diff --git a/Source/ERF.H b/Source/ERF.H
@@ -123,17 +123,16 @@ public:
 
 #ifdef ERF_USE_POISSON_SOLVE
     // Project the velocities to be divergence-free
-    void project_velocities (                                        amrex::Vector<amrex::MultiFab >& vars);
-    void project_velocities (int lev_min, int lev_max, amrex::Vector<amrex::Vector<amrex::MultiFab>>& vars);
+    void project_velocities (int lev, amrex::Real dt, amrex::Vector<amrex::MultiFab >& vars, amrex::MultiFab& p);
+
+    // Project the velocities to be divergence-free with a thin body
+    void project_velocities_tb (int lev, amrex::Real dt, amrex::Vector<amrex::MultiFab >& vars, amrex::MultiFab& p);
 
     // Define the projection bc's based on the domain bc types
     amrex::Array<amrex::LinOpBCType,AMREX_SPACEDIM>
       get_projection_bc (amrex::Orientation::Side side) const noexcept;
     bool projection_has_dirichlet (amrex::Array<amrex::LinOpBCType,AMREX_SPACEDIM> bcs) const;
 
-    // Project the velocities to be divergence-free with a thin body
-    void project_velocities_tb (              amrex::Vector<amrex::MultiFab >& vars, const amrex::Real dt);
-    void project_velocities_tb (amrex::Vector<amrex::Vector<amrex::MultiFab>>& vars, const amrex::Real dt);
 #endif
 
     // Init (NOT restart or regrid)
@@ -605,6 +604,10 @@ private:
 #endif
     amrex::Vector<std::unique_ptr<MRISplitIntegrator<amrex::Vector<amrex::MultiFab> > > > mri_integrator_mem;
 
+#ifdef ERF_USE_POISSON_SOLVE
+    amrex::Vector<amrex::MultiFab> pp_inc;
+#endif
+
     // Vector over levels of routines to impose physical boundary conditions
     amrex::Vector<std::unique_ptr<ERFPhysBCFunct_cons>> physbcs_cons;
     amrex::Vector<std::unique_ptr<ERFPhysBCFunct_u>> physbcs_u;
@@ -840,6 +843,9 @@ private:
 #endif
 
     static int verbose;
+#ifdef ERF_USE_POISSON_SOLVE
+    static int mg_verbose;
+#endif
 
     // Diagnostic output interval
     static int sum_interval;

diff --git a/Source/ERF.cpp b/Source/ERF.cpp
@@ -38,6 +38,9 @@ int  ERF::fixed_mri_dt_ratio = 0;
 
 // Dictate verbosity in screen output
 int ERF::verbose       = 0;
+#ifdef ERF_USE_POISSON_SOLVE
+int ERF::mg_verbose    = 0;
+#endif
 
 // Frequency of diagnostic output
 int  ERF::sum_interval  = -1;
@@ -182,6 +185,11 @@ ERF::ERF ()
     vars_new.resize(nlevs_max);
     vars_old.resize(nlevs_max);
 
+#ifdef ERF_USE_POISSON_SOLVE
+    pp_inc.resize(nlevs_max);
+#endif
+
+
     rU_new.resize(nlevs_max);
     rV_new.resize(nlevs_max);
     rW_new.resize(nlevs_max);
@@ -819,7 +827,12 @@ ERF::InitData ()
         // Note -- this projection is only defined for no terrain
         if (solverChoice.project_initial_velocity) {
             AMREX_ALWAYS_ASSERT(solverChoice.use_terrain == 0);
-            project_velocities(0, finest_level, vars_new);
+            Real dummy_dt = 1.0;
+            for (int lev = 0; lev <= finest_level; ++lev)
+            {
+                project_velocities(lev, dummy_dt, vars_new[lev], pp_inc[lev]);
+                pp_inc[lev].setVal(0.);
+            }
         }
     }
 #endif
@@ -1228,6 +1241,9 @@ ERF::ReadParameters ()
 
         // Verbosity
         pp.query("v", verbose);
+#ifdef ERF_USE_POISSON_SOLVE
+        pp.query("mg_v", mg_verbose);
+#endif
 
         // Frequency of diagnostic output
         pp.query("sum_interval", sum_interval);

diff --git a/Source/ERF_make_new_arrays.cpp b/Source/ERF_make_new_arrays.cpp
@@ -141,6 +141,11 @@ ERF::init_stuff (int lev, const BoxArray& ba, const DistributionMapping& dm,
     lev_new[Vars::zvel].define(convert(ba, IntVect(0,0,1)), dm, 1, IntVect(ngrow_vels,ngrow_vels,0));
     lev_old[Vars::zvel].define(convert(ba, IntVect(0,0,1)), dm, 1, IntVect(ngrow_vels,ngrow_vels,0));
 
+#ifdef ERF_USE_POISSON_SOLVE
+    pp_inc[lev].define(ba, dm, 1, 1);
+    pp_inc[lev].setVal(0.0);
+#endif
+
     // ********************************************************************************************
     // These are just used for scratch in the time integrator but we might as well define them here
     // ********************************************************************************************

diff --git a/Source/ERF_make_new_level.cpp b/Source/ERF_make_new_level.cpp
@@ -499,6 +499,10 @@ ERF::ClearLevel (int lev)
     rW_new[lev].clear();
     rW_old[lev].clear();
 
+#ifdef ERF_USE_POISSON_SOLVE
+    pp_inc[lev].clear();
+#endif
+
     // Clears the integrator memory
     mri_integrator_mem[lev].reset();
 

diff --git a/Source/IO/Plotfile.cpp b/Source/IO/Plotfile.cpp
@@ -80,16 +80,7 @@ ERF::setPlotVariables (const std::string& pp_plot_var_names, Vector<std::string>
     }
 #endif
 
-#ifdef ERF_USE_EB
-    if (containerHasElement(plot_var_names, "volfrac")) {
-        tmp_plot_names.push_back("volfrac");
-    }
-#endif
-
     plot_var_names = tmp_plot_names;
-
-    for (int i=0; i<plot_var_names.size(); i++) {
-    }
 }
 
 void
@@ -1040,6 +1031,13 @@ ERF::WritePlotFile (int which, Vector<std::string> plot_var_names)
         }
 #endif
 
+#ifdef ERF_USE_POISSON_SOLVE
+        if (containerHasElement(plot_var_names, "pp_inc")) {
+            MultiFab::Copy(mf[lev], pp_inc[lev], 0, mf_comp, 1, 0);
+            mf_comp += 1;
+        }
+#endif
+
 #ifdef ERF_COMPUTE_ERROR
         // Next, check for error in velocities and if desired, output them -- note we output none or all, not just some
         if (containerHasElement(plot_var_names, "xvel_err") ||