forest-->forest_drag and terrain-->terrain_drag (#2022)

Source/DataStructs/ERF_DataStruct.H

@@ -688,9 +688,10 @@ struct SolverChoice {
     amrex::Real windfarm_x_shift = -1.0;
     amrex::Real windfarm_y_shift = -1.0;
 
-    // Flag for valid canopy model
-    bool do_forest {false};
-    // Immersed Forcing
-    bool do_terrain {false};
+    // Use forest canopy model?
+    bool do_forest_drag {false};
+
+    // Use immersed forcing representation of terrain?
+    bool do_terrain_drag {false};
 };
 #endif

Source/ERF.H

@@ -1157,8 +1157,8 @@ private:
     std::unique_ptr<WriteBndryPlanes> m_w2d  = nullptr;
     std::unique_ptr<ReadBndryPlanes>  m_r2d  = nullptr;
     std::unique_ptr<ABLMost>          m_most = nullptr;
-    amrex::Vector<std::unique_ptr<ForestDrag>> m_forest;
-    amrex::Vector<std::unique_ptr<TerrainDrag>> m_terrain;
+    amrex::Vector<std::unique_ptr<ForestDrag>> m_forest_drag;
+    amrex::Vector<std::unique_ptr<TerrainDrag>> m_terrain_drag;
 
     //
     // Holds info for dynamically generated tagging criteria

Source/ERF.cpp

@@ -140,12 +140,12 @@ ERF::ERF_shared ()
     lsm_flux.resize(nlevs_max);
 
     // NOTE: size canopy model before readparams (if file exists, we construct)
-    m_forest.resize(nlevs_max);
-    for (int lev = 0; lev < max_level; ++lev) { m_forest[lev] = nullptr;}
+    m_forest_drag.resize(nlevs_max);
+    for (int lev = 0; lev < max_level; ++lev) { m_forest_drag[lev] = nullptr;}
 
-    // Immersed Forcing
-    m_terrain.resize(nlevs_max);
-    for (int lev = 0; lev < max_level; ++lev) { m_terrain[lev] = nullptr;}
+    // Immersed Forcing Representation of Terrain
+    m_terrain_drag.resize(nlevs_max);
+    for (int lev = 0; lev < max_level; ++lev) { m_terrain_drag[lev] = nullptr;}
 
 
     ReadParameters();
@@ -1628,19 +1628,19 @@ ERF::ReadParameters ()
 
         // Query the canopy model file name
         std::string forestfile;
-        solverChoice.do_forest = pp.query("forest_file", forestfile);
-        if (solverChoice.do_forest) {
+        solverChoice.do_forest_drag = pp.query("forest_file", forestfile);
+        if (solverChoice.do_forest_drag) {
             for (int lev = 0; lev <= max_level; ++lev) {
-                m_forest[lev] = std::make_unique<ForestDrag>(forestfile);
+                m_forest_drag[lev] = std::make_unique<ForestDrag>(forestfile);
             }
         }
 
         //Query the terrain file name
         std::string terrainfile;
-        solverChoice.do_terrain = pp.query("terrain_file", terrainfile);
-        if (solverChoice.do_terrain) {
+        solverChoice.do_terrain_drag = pp.query("terrain_file", terrainfile);
+        if (solverChoice.do_terrain_drag) {
             for (int lev = 0; lev <= max_level; ++lev) {
-                m_terrain[lev] = std::make_unique<TerrainDrag>(terrainfile);
+                m_terrain_drag[lev] = std::make_unique<TerrainDrag>(terrainfile);
             }
         }
 

Source/ERF_MakeNewLevel.cpp

@@ -141,9 +141,9 @@ void ERF::MakeNewLevelFromScratch (int lev, Real time, const BoxArray& ba_in,
     // ********************************************************************************************
     // Build the data structures for canopy model (depends upon z_phys)
     // ********************************************************************************************
-    if (solverChoice.do_forest) { m_forest[lev]->define_drag_field(ba, dm, geom[lev], z_phys_nd[lev].get()); }
+    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) { m_terrain[lev]->define_terrain_blank_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()); }
     //********************************************************************************************
     // Microphysics
     // *******************************************************************************************
@@ -221,9 +221,9 @@ ERF::MakeNewLevelFromCoarse (int lev, Real time, const BoxArray& ba,
     // ********************************************************************************************
     // Build the data structures for canopy model (depends upon z_phys)
     // ********************************************************************************************
-    if (solverChoice.do_forest) { m_forest[lev]->define_drag_field(ba, dm, geom[lev], z_phys_nd[lev].get()); }
+    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) { m_terrain[lev]->define_terrain_blank_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()); }
     //********************************************************************************************
     // Microphysics
     // *******************************************************************************************
@@ -351,9 +351,9 @@ ERF::RemakeLevel (int lev, Real time, const BoxArray& ba, const DistributionMapp
     // ********************************************************************************************
     // Build the data structures for canopy model (depends upon z_phys)
     // ********************************************************************************************
-    if (solverChoice.do_forest) { m_forest[lev]->define_drag_field(ba, dm, geom[lev], z_phys_nd[lev].get()); }
+    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) { m_terrain[lev]->define_terrain_blank_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 the physbcs objects (after initializing the terrain but before calling FillCoarsePatch
     // *****************************************************************************************************

Source/IO/ERF_Plotfile.cpp

@@ -491,7 +491,7 @@ ERF::WritePlotFile (int which, PlotFileType plotfile_type, Vector<std::string> p
 
         if (containerHasElement(plot_var_names, "terrain_IB_mask"))
         {
-            MultiFab* terrain_blank = m_terrain[lev]->get_terrain_blank_field();
+            MultiFab* terrain_blank = m_terrain_drag[lev]->get_terrain_blank_field();
 #ifdef _OPENMP
 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
 #endif

Source/SourceTerms/ERF_MakeMomSources.cpp

@@ -90,16 +90,16 @@ void make_mom_sources (int level,
     //    9. Forest canopy
     //   10. Immersed Forcing
     // *****************************************************************************
-    const bool l_use_ndiff    = solverChoice.use_NumDiff;
-    const bool use_terrain    = solverChoice.terrain_type != TerrainType::None;
-    const bool l_do_forest    = solverChoice.do_forest;
-    const bool l_do_terrain   = solverChoice.do_terrain;
+    const bool l_use_ndiff       = solverChoice.use_NumDiff;
+    const bool l_use_zphys       = (solverChoice.mesh_type != MeshType::ConstantDz);
+    const bool l_do_forest_drag  = solverChoice.do_forest_drag;
+    const bool l_do_terrain_drag = solverChoice.do_terrain_drag;
 
     // Check if terrain and immersed terrain clash
-    if(use_terrain && l_do_terrain){
-        amrex::Error(" Cannot use immersed forcing with terrain");
+    if(l_use_zphys && l_do_terrain_drag){
+        amrex::Error(" Cannot use immersed forcing for terrain with terrain-fitted coordinates");
     }
-    if(l_do_forest && l_do_terrain){
+    if(l_do_forest_drag && l_do_terrain_drag){
         amrex::Error(" Currently forest canopy cannot be used with immersed forcing");
     }
 
@@ -243,9 +243,9 @@ void make_mom_sources (int level,
         const Array4<const Real>& t_blank_arr = (terrain_blank) ? terrain_blank->const_array(mfi) :
                                                                Array4<const Real>{};
 
-        const Array4<const Real>& z_nd_arr = (use_terrain) ? z_phys_nd->const_array(mfi) :
+        const Array4<const Real>& z_nd_arr = (l_use_zphys) ? z_phys_nd->const_array(mfi) :
                                                              Array4<const Real>{};
-        const Array4<const Real>& z_cc_arr = (use_terrain) ? z_phys_cc->const_array(mfi) :
+        const Array4<const Real>& z_cc_arr = (l_use_zphys) ? z_phys_cc->const_array(mfi) :
                                                              Array4<const Real>{};
 
         // *****************************************************************************
@@ -505,7 +505,7 @@ void make_mom_sources (int level,
         // *****************************************************************************
         // 9. Add CANOPY source terms
         // *****************************************************************************
-        if (l_do_forest) {
+        if (l_do_forest_drag) {
             ParallelFor(tbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
             {
                 const Real ux = u(i, j, k);
@@ -543,7 +543,7 @@ void make_mom_sources (int level,
         // *****************************************************************************
         // 10. Add Immersed source terms
         // *****************************************************************************
-        if (l_do_terrain) {
+        if (l_do_terrain_drag) {
             const Real drag_coefficient=10.0/dz;
             const Real tiny = std::numeric_limits<amrex::Real>::epsilon();
             ParallelFor(tbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept

Source/TimeIntegration/ERF_TI_slow_rhs_fun.H

@@ -62,10 +62,10 @@
 
         // Canopy data for mom sources
         MultiFab* forest_drag = nullptr;
-        if (solverChoice.do_forest) { forest_drag = m_forest[level]->get_drag_field(); }
+        if (solverChoice.do_forest_drag) { forest_drag = m_forest_drag[level]->get_drag_field(); }
         // Immersed Forcing
         MultiFab* terrain_blank = nullptr;
-        if(solverChoice.do_terrain) { terrain_blank = m_terrain[level]->get_terrain_blank_field(); }
+        if(solverChoice.do_terrain_drag) { terrain_blank = m_terrain_drag[level]->get_terrain_blank_field(); }
 
         // Moving terrain
         if ( solverChoice.terrain_type == TerrainType::Moving )
@@ -435,10 +435,10 @@
 
         // Canopy data for mom sources
         MultiFab* forest_drag = nullptr;
-        if (solverChoice.do_forest) { forest_drag = m_forest[level]->get_drag_field(); }
+        if (solverChoice.do_forest_drag) { forest_drag = m_forest_drag[level]->get_drag_field(); }
         // Immersed Forcing
         MultiFab* terrain_blank = nullptr;
-        if(solverChoice.do_terrain) { terrain_blank = m_terrain[level]->get_terrain_blank_field(); }
+        if(solverChoice.do_terrain_drag) { terrain_blank = m_terrain_drag[level]->get_terrain_blank_field(); }
 
         make_sources(level, nrk, slow_dt, old_stage_time, S_data, S_prim, cc_src, z_phys_cc[level],
 #if defined(ERF_USE_RRTMGP)