RANS Part 1: Poisson wall distance (#2031)

* Rename init_immersed_body to avoid confusion * Create wall dist array * Add RANS options, begin implementing walldist calc, add plotfile output Poisson solve and dist calc implemented for uniform dz for now * Fix boxarray, force 3D solve, call setLevelBC * Add masking to walldist poisson solve * Simplify code, use nodal solver instead * Update mask for nodal solver * Update walldist calc to properly handle nodal values Don't use poisson dist for thin bodies; directly calculate wall dist in trivial case of flat bottom boundary * Add clarifying print statements * Fix direct calc * GPU fix * Remove debug print, additional GPU fixes --------- Co-authored-by: Aaron M. Lattanzi <[email protected]>
erf-model · Dec 19, 2024 · e63a233 · e63a233
1 parent c864f79
commit e63a233
Show file tree

Hide file tree

Showing 12 changed files with 375 additions and 22 deletions.
diff --git a/CMake/BuildERFExe.cmake b/CMake/BuildERFExe.cmake
@@ -190,6 +190,7 @@ function(build_erf_lib erf_lib_name)
        ${SRC_DIR}/Utils/ERF_MomentumToVelocity.cpp
        ${SRC_DIR}/LinearSolvers/ERF_PoissonSolve.cpp
        ${SRC_DIR}/LinearSolvers/ERF_PoissonSolve_tb.cpp
+       ${SRC_DIR}/LinearSolvers/ERF_PoissonWallDist.cpp
        ${SRC_DIR}/LinearSolvers/ERF_ComputeDivergence.cpp 
        ${SRC_DIR}/LinearSolvers/ERF_SolveWithGMRES.cpp
        ${SRC_DIR}/LinearSolvers/ERF_SolveWithMLMG.cpp

diff --git a/CMake/SetAmrexOptions.cmake b/CMake/SetAmrexOptions.cmake
@@ -29,7 +29,7 @@ set(AMReX_FORTRAN OFF)
 
 set(AMReX_LINEAR_SOLVERS ON)
 set(AMReX_LINEAR_SOLVERS_EM OFF)
-set(AMReX_LINEAR_SOLVERS_INCFLO OFF)
+set(AMReX_LINEAR_SOLVERS_INCFLO ON)
 
 set(AMReX_PARTICLES OFF)
 if(ERF_ENABLE_PARTICLES)

diff --git a/Source/DataStructs/ERF_AdvStruct.H b/Source/DataStructs/ERF_AdvStruct.H
@@ -184,6 +184,7 @@ struct AdvChoice {
         pp.queryarr("zero_xflux_faces", zero_xflux);
         pp.queryarr("zero_yflux_faces", zero_yflux);
         pp.queryarr("zero_zflux_faces", zero_zflux);
+        have_zero_flux_faces = ((zero_xflux.size() > 0) || (zero_yflux.size() > 0) || (zero_zflux.size() > 0));
     }
 
     void display()
@@ -295,9 +296,10 @@ struct AdvChoice {
     amrex::Real moistscal_horiz_upw_frac = 1.0;
     amrex::Real moistscal_vert_upw_frac  = 1.0;
 
-    // Thin immersed body
+    // Thin immersed bodies
     amrex::Vector<amrex::IntVect> zero_xflux;
     amrex::Vector<amrex::IntVect> zero_yflux;
     amrex::Vector<amrex::IntVect> zero_zflux;
+    bool have_zero_flux_faces = false;
 };
 #endif
diff --git a/Source/DataStructs/ERF_TurbStruct.H b/Source/DataStructs/ERF_TurbStruct.H
@@ -5,6 +5,8 @@
 
 AMREX_ENUM(LESType, None, Smagorinsky, Deardorff);
 
+AMREX_ENUM(RANSType, None, kEqn);
+
 AMREX_ENUM(PBLType, None, MYNN25, YSU);
 
 template <typename T>
@@ -36,6 +38,10 @@ struct TurbChoice {
         std::string les_type_string = "None";
         query_one_or_per_level(pp, "les_type", les_type, lev, max_level);
 
+        // Which RANS closure?
+        std::string rans_type_string = "None";
+        query_one_or_per_level(pp, "rans_type", rans_type, lev, max_level);
+
         // Which PBL Closure
         static std::string pbl_type_string = "None";
         query_one_or_per_level(pp, "pbl_type", pbl_type, lev, max_level);
@@ -188,6 +194,9 @@ struct TurbChoice {
 
     amrex::Real theta_ref = 300.0;
 
+    // RANS type
+    RANSType rans_type;
+
     // PBL model
     PBLType pbl_type;
 

diff --git a/Source/ERF.H b/Source/ERF.H
@@ -177,6 +177,9 @@ public:
     void project_velocities_tb (int lev, amrex::Real dt, amrex::Vector<amrex::MultiFab >& vars,
                                 amrex::MultiFab& p);
 
+    // Calculate wall distance by solving a Poisson equation
+    void poisson_wall_dist (int lev);
+
 #ifdef ERF_USE_FFT
     void solve_with_fft (int lev, amrex::MultiFab& rhs, amrex::MultiFab& p,
                          amrex::Array<amrex::MultiFab,AMREX_SPACEDIM>& fluxes);
@@ -404,7 +407,7 @@ public:
 
     void init_from_hse (int lev);
 
-    void init_immersed_body (int lev, const amrex::BoxArray& ba, const amrex::DistributionMapping& dm);
+    void init_thin_body (int lev, const amrex::BoxArray& ba, const amrex::DistributionMapping& dm);
 
 #ifdef ERF_USE_MULTIBLOCK
     // constructor used when ERF is created by a multiblock driver
@@ -853,6 +856,9 @@ private:
 
     amrex::Vector<std::unique_ptr<amrex::MultiFab>> z_t_rk;
 
+    // Wall distance function
+    amrex::Vector<std::unique_ptr<amrex::MultiFab>>  walldist;
+
     // Map scale factors
     amrex::Vector<std::unique_ptr<amrex::MultiFab>> mapfac_m;
     amrex::Vector<std::unique_ptr<amrex::MultiFab>> mapfac_u;
@@ -980,8 +986,10 @@ private:
                                                     "Kmv","Kmh",
                                                     // eddy diffusivity of heat
                                                     "Khv","Khh",
-                                                    // mynn pbl lengthscale
+                                                    // turbulence lengthscale
                                                     "Lturb",
+                                                    // wall distance
+                                                    "walldist",
                                                     // moisture vars
                                                     "qt", "qv", "qc", "qi", "qp", "qrain", "qsnow", "qgraup", "qsat",
                                                     "rain_accum", "snow_accum", "graup_accum",

diff --git a/Source/ERF.cpp b/Source/ERF.cpp
@@ -288,6 +288,9 @@ ERF::ERF_shared ()
 
     z_t_rk.resize(nlevs_max);
 
+    // Wall distance
+    walldist.resize(nlevs_max);
+
     // Mapfactors
     mapfac_m.resize(nlevs_max);
     mapfac_u.resize(nlevs_max);
@@ -669,9 +672,7 @@ ERF::InitData_post ()
             AverageDown();
         }
 
-        if ((solverChoice.advChoice.zero_xflux.size() > 0) ||
-            (solverChoice.advChoice.zero_yflux.size() > 0) ||
-            (solverChoice.advChoice.zero_zflux.size() > 0))
+        if (solverChoice.advChoice.have_zero_flux_faces)
         {
             AMREX_ALWAYS_ASSERT_WITH_MESSAGE(finest_level == 0,
                 "Thin immersed body with refinement not currently supported.");

diff --git a/Source/ERF_MakeNewArrays.cpp b/Source/ERF_MakeNewArrays.cpp
@@ -87,16 +87,26 @@ ERF::init_stuff (int lev, const BoxArray& ba, const DistributionMapping& dm,
                z_t_rk[lev] = nullptr;
     }
 
-   // We use these area arrays regardless of terrain, EB or none of the above
-   detJ_cc[lev] = std::make_unique<MultiFab>(ba,dm,1,1);
-        ax[lev] = std::make_unique<MultiFab>(convert(ba,IntVect(1,0,0)),dm,1,1);
-        ay[lev] = std::make_unique<MultiFab>(convert(ba,IntVect(0,1,0)),dm,1,1);
-        az[lev] = std::make_unique<MultiFab>(convert(ba,IntVect(0,0,1)),dm,1,1);
+    // We use these area arrays regardless of terrain, EB or none of the above
+    detJ_cc[lev] = std::make_unique<MultiFab>(ba,dm,1,1);
+         ax[lev] = std::make_unique<MultiFab>(convert(ba,IntVect(1,0,0)),dm,1,1);
+         ay[lev] = std::make_unique<MultiFab>(convert(ba,IntVect(0,1,0)),dm,1,1);
+         az[lev] = std::make_unique<MultiFab>(convert(ba,IntVect(0,0,1)),dm,1,1);
 
-   detJ_cc[lev]->setVal(1.0);
-        ax[lev]->setVal(1.0);
-        ay[lev]->setVal(1.0);
-        az[lev]->setVal(1.0);
+    detJ_cc[lev]->setVal(1.0);
+         ax[lev]->setVal(1.0);
+         ay[lev]->setVal(1.0);
+         az[lev]->setVal(1.0);
+
+    // ********************************************************************************************
+    // Create wall distance array for RANS modeling
+    // ********************************************************************************************
+    if (solverChoice.turbChoice[lev].rans_type != RANSType::None) {
+        walldist[lev] = std::make_unique<MultiFab>(ba,dm,1,1);
+        walldist[lev]->setVal(1e23);
+    } else {
+        walldist[lev] = nullptr;
+    }
 
     // ********************************************************************************************
     // These are the persistent containers for the old and new data
@@ -237,7 +247,7 @@ ERF::init_stuff (int lev, const BoxArray& ba, const DistributionMapping& dm,
 
 #if defined(ERF_USE_WINDFARM)
     //*********************************************************
-    // Variables for Ftich model for windfarm parametrization
+    // Variables for Fitch model for windfarm parametrization
     //*********************************************************
     if (solverChoice.windfarm_type == WindFarmType::Fitch){
         vars_windfarm[lev].define(ba, dm, 5, ngrow_state); // V, dVabsdt, dudt, dvdt, dTKEdt

diff --git a/Source/ERF_MakeNewLevel.cpp b/Source/ERF_MakeNewLevel.cpp
@@ -101,7 +101,7 @@ void ERF::MakeNewLevelFromScratch (int lev, Real time, const BoxArray& ba_in,
     // ********************************************************************************************
     // Thin immersed body
     // *******************************************************************************************
-    init_immersed_body(lev, ba, dm);
+    init_thin_body(lev, ba, dm);
 
     // ********************************************************************************************
     // Initialize the integrator class
@@ -144,6 +144,14 @@ void ERF::MakeNewLevelFromScratch (int lev, Real time, const BoxArray& ba_in,
     if (solverChoice.do_forest_drag) { m_forest_drag[lev]->define_drag_field(ba, dm, geom[lev], z_phys_nd[lev].get()); }
 
     if (solverChoice.do_terrain_drag) { m_terrain_drag[lev]->define_terrain_blank_field(ba, dm, geom[lev], z_phys_nd[lev].get()); }
+
+    //********************************************************************************************
+    // Create wall distance field for RANS model (depends upon z_phys)
+    // *******************************************************************************************
+    if (solverChoice.turbChoice[lev].rans_type != RANSType::None) {
+        poisson_wall_dist(lev);
+    }
+
     //********************************************************************************************
     // Microphysics
     // *******************************************************************************************
@@ -496,7 +504,7 @@ ERF::ClearLevel (int lev)
 }
 
 void
-ERF::init_immersed_body (int lev, const BoxArray& ba, const DistributionMapping& dm)
+ERF::init_thin_body (int lev, const BoxArray& ba, const DistributionMapping& dm)
 {
     //********************************************************************************************
     // Thin immersed body
@@ -544,7 +552,7 @@ ERF::init_immersed_body (int lev, const BoxArray& ba, const DistributionMapping&
     }
 
     if (solverChoice.advChoice.zero_yflux.size() > 0) {
-        amrex::Print() << "Setting up thin interface boundary for "
+        amrex::Print() << "Setting up thin immersed body for "
             << solverChoice.advChoice.zero_yflux.size() << " yfaces" << std::endl;
         BoxArray ba_yf(ba);
         ba_yf.surroundingNodes(1);
@@ -576,7 +584,7 @@ ERF::init_immersed_body (int lev, const BoxArray& ba, const DistributionMapping&
     }
 
     if (solverChoice.advChoice.zero_zflux.size() > 0) {
-        amrex::Print() << "Setting up thin interface boundary for "
+        amrex::Print() << "Setting up thin immersed body for "
             << solverChoice.advChoice.zero_zflux.size() << " zfaces" << std::endl;
         BoxArray ba_zf(ba);
         ba_zf.surroundingNodes(2);

diff --git a/Source/IO/ERF_Plotfile.cpp b/Source/IO/ERF_Plotfile.cpp
@@ -1045,6 +1045,10 @@ ERF::WritePlotFile (int which, PlotFileType plotfile_type, Vector<std::string> p
             MultiFab::Copy(mf[lev],*eddyDiffs_lev[lev],EddyDiff::Turb_lengthscale,mf_comp,1,0);
             mf_comp ++;
         }
+        if (containerHasElement(plot_var_names, "walldist")) {
+            MultiFab::Copy(mf[lev],*walldist[lev],0,mf_comp,1,0);
+            mf_comp ++;
+        }
 
         // TODO: The size of the q variables can vary with different
         //       moisture models. Therefore, certain components may