Merge pull request erf-model#1276 from ewquon/grid_stretching

Add options to specify grid stretching

Docs/sphinx_doc/Inputs.rst

@@ -164,6 +164,61 @@ Examples of Usage
      level-1 grids will be created every 2 level-0 time steps, and new
      level-2 grids will be created every 2 level-1 time steps.
 
+
+Grid Stretching
+===============
+
+This automatically activates **erf.use_terrain**. By default, the
+problem-specific terrain is initialized to be flat at an elevation of z=0.
+These inputs are used to automatically generate the staggered z levels used to
+calculate the grid metric transformation. Alternatively, arbitrary z levels may
+be specified with the **erf.terrain_z_levels** parameter.
+
+.. _list-of-parameters-3:
+
+List of Parameters
+------------------
+
++-------------------------------+-----------------+-----------------+-------------+
+| Parameter                     | Definition      | Acceptable      | Default     |
+|                               |                 | Values          |             |
++===============================+=================+=================+=============+
+| **erf.grid_stretching_ratio** | scaling factor  | Real > 0        | 0 (no grid  |
+|                               | applied to      |                 | stretching) |
+|                               | delta z at each |                 |             |
+|                               | level           |                 |             |
++-------------------------------+-----------------+-----------------+-------------+
+| **erf.initial_dz**            | vertical grid   | Real > 0        | must be set |
+|                               | spacing for the |                 | if grid     |
+|                               | cell above the  |                 | stretching  |
+|                               | bottom surface  |                 | ratio is set|
++-------------------------------+-----------------+-----------------+-------------+
+| **erf.terrain_z_levels**      | nominal         | List of Real    | NONE        |
+|                               | staggered       |                 |             |
+|                               | z levels        |                 |             |
++-------------------------------+-----------------+-----------------+-------------+
+
+.. _notes-3:
+
+Notes
+-----
+
+- If both **erf.terrain_z_levels** and **erf.grid_stretching_ratio** are
+  specified, the simple grid stretching will be ignored.
+- The number of input **erf.terrain_z_levels** must be equal **amr.n_cell** in
+  the z direction + 1.
+
+.. _examples-of-usage-3:
+
+Examples of Usage
+-----------------
+
+-  **erf.grid_stretching_ratio** = 1.025
+
+-  | **erf.initial_dz** = 5.0
+   | the first cell center would be at z=2.5
+
+
 Regridding
 ==========
 

Exec/ABL/inputs_most_pbl

@@ -0,0 +1,65 @@
+# ------------------  INPUTS TO MAIN PROGRAM  -------------------
+max_step = 3000
+
+amrex.fpe_trap_invalid = 1
+
+fabarray.mfiter_tile_size = 1024 1024 1024
+
+# PROBLEM SIZE & GEOMETRY
+geometry.prob_extent =  1024    1024    1024
+amr.n_cell           =     4       4      64
+
+#erf.grid_stretching_ratio = 1.025
+#erf.initial_dz = 16.0
+
+geometry.is_periodic = 1 1 0
+
+# MOST BOUNDARY (DEFAULT IS ADIABATIC FOR THETA)
+zlo.type      = "Most"
+erf.most.z0   = 0.1
+erf.most.zref = 8.0
+
+zhi.type = "SlipWall"
+
+# TIME STEP CONTROL
+erf.fixed_dt           = 0.1  # fixed time step depending on grid resolution
+
+# DIAGNOSTICS & VERBOSITY
+erf.sum_interval   = 1       # timesteps between computing mass
+erf.v              = 1       # verbosity in ERF.cpp
+amr.v                = 1       # verbosity in Amr.cpp
+
+# REFINEMENT / REGRIDDING
+amr.max_level       = 0       # maximum level number allowed
+
+# CHECKPOINT FILES
+erf.check_file      = chk       # root name of checkpoint file
+erf.check_int       = -1        # number of timesteps between checkpoints
+
+# PLOTFILES
+erf.plot_file_1     = plt       # prefix of plotfile name
+erf.plot_int_1      = 100       # number of timesteps between plotfiles
+erf.plot_vars_1     = density rhoadv_0 x_velocity y_velocity z_velocity pressure temp theta
+
+# SOLVER CHOICE
+erf.alpha_T = 0.0
+erf.alpha_C = 1.0
+erf.use_gravity = false
+
+erf.molec_diff_type = "None"
+erf.les_type        = "Smagorinsky"
+erf.Cs              = 0.1
+
+erf.pbl_type        = "MYNN2.5"
+
+erf.init_type = "uniform"
+
+# PROBLEM PARAMETERS
+prob.rho_0 = 1.0
+prob.A_0 = 1.0
+
+prob.U_0 = 10.0
+prob.V_0 = 0.0
+prob.W_0 = 0.0
+prob.T_0 = 300.0
+

Source/DataStructs/DataStruct.H

@@ -37,6 +37,18 @@ struct SolverChoice {
         // Do we have terrain (or grid stretching)?
         pp.query("use_terrain", use_terrain);
 
+        pp.query("grid_stretching_ratio", grid_stretching_ratio);
+        AMREX_ASSERT_WITH_MESSAGE((grid_stretching_ratio >= 0.),
+                                  "The grid stretching ratio must be greater than 0");
+        if (grid_stretching_ratio > 0) {
+            if (!use_terrain) {
+                amrex::Print() << "Turning terrain on to enable grid stretching" << std::endl;
+                use_terrain = true;
+            }
+            pp.query("zsurface", zsurf);
+            pp.get("initial_dz", dz0);
+        }
+
         // Do we set map scale factors to 0.5 instead of 1 for testing?
         pp.query("test_mapfactor", test_mapfactor);
 
@@ -266,6 +278,11 @@ struct SolverChoice {
     amrex::Real c_p = Cp_d; // specific heat at constant pressure for dry air [J/(kg-K)]
     amrex::Real rdOcp;
 
+    // Staggered z levels for vertical grid stretching
+    amrex::Real grid_stretching_ratio   = 0;
+    amrex::Real zsurf                   = 0.0;
+    amrex::Real dz0;
+
 #if defined(ERF_USE_WARM_NO_PRECIP)
     amrex::Real tau_cond = 1.0; // Default time of 1 sec -- this is somewhat arbitray
 #endif

Source/ERF.H

@@ -550,6 +550,8 @@ private:
     amrex::Vector<std::unique_ptr<amrex::MultiFab>> SFS_hfx1_lev, SFS_hfx2_lev, SFS_hfx3_lev;
     amrex::Vector<std::unique_ptr<amrex::MultiFab>> SFS_diss_lev;
 
+    // Terrain / grid stretching
+    amrex::Vector<amrex::Real> zlevels_stag;
     amrex::Vector<std::unique_ptr<amrex::MultiFab>> z_phys_nd;
     amrex::Vector<std::unique_ptr<amrex::MultiFab>> z_phys_cc;
     amrex::Vector<std::unique_ptr<amrex::MultiFab>>   detJ_cc;

Source/ERF.cpp

@@ -122,6 +122,14 @@ ERF::ERF ()
 
     // Geometry on all levels has been defined already.
 
+    if (solverChoice.use_terrain) {
+        init_zlevels(zlevels_stag,
+                     geom[0],
+                     solverChoice.grid_stretching_ratio,
+                     solverChoice.zsurf,
+                     solverChoice.dz0);
+    }
+
     // No valid BoxArray and DistributionMapping have been defined.
     // But the arrays for them have been resized.
 

Source/ERF_make_new_level.cpp

@@ -410,7 +410,7 @@ ERF::update_terrain_arrays (int lev, Real time)
     if (solverChoice.use_terrain) {
         if (init_type != "real" && init_type != "metgrid") {
             prob->init_custom_terrain(geom[lev],*z_phys_nd[lev],time);
-            init_terrain_grid(geom[lev],*z_phys_nd[lev]);
+            init_terrain_grid(geom[lev],*z_phys_nd[lev],zlevels_stag);
         }
         if (lev>0 && (init_type == "real" || init_type == "metgrid")) {
             PhysBCFunctNoOp empty_bc;

Source/Initialization/ERF_init_from_metgrid.cpp

@@ -103,7 +103,7 @@ ERF::init_from_metgrid(int lev)
     } // mf
 
     // This defines all the z(i,j,k) values given z(i,j,0) from above.
-    init_terrain_grid(geom[lev], *z_phys);
+    init_terrain_grid(geom[lev], *z_phys, zlevels_stag);
 
     // This makes the Jacobian.
     make_J  (geom[lev],*z_phys,  *detJ_cc[lev]);

Source/TimeIntegration/TI_fast_rhs_fun.H

@@ -29,13 +29,13 @@ auto fast_rhs_fun = [&](int fast_step, int n_sub,
             // Make "old" fast geom -- store in z_phys_nd for convenience
             if (verbose) Print() << "Making geometry at start of substep time: " << old_substep_time << std::endl;
             prob->init_custom_terrain(fine_geom,*z_phys_nd[level],old_substep_time);
-            init_terrain_grid  (fine_geom,*z_phys_nd[level]);
+            init_terrain_grid  (fine_geom,*z_phys_nd[level], zlevels_stag);
             make_J             (fine_geom,*z_phys_nd[level], *detJ_cc[level]);
 
             // Make "new" fast geom
             if (verbose) Print() << "Making geometry for end of substep time :" << new_substep_time << std::endl;
             prob->init_custom_terrain(fine_geom,*z_phys_nd_new[level],new_substep_time);
-            init_terrain_grid  (fine_geom,*z_phys_nd_new[level]);
+            init_terrain_grid  (fine_geom,*z_phys_nd_new[level], zlevels_stag);
             make_J             (fine_geom,*z_phys_nd_new[level], *detJ_cc_new[level]);
 
             Real inv_dt   = 1./dtau;

Source/TimeIntegration/TI_slow_rhs_fun.H

@@ -29,17 +29,17 @@
 
             if (verbose) Print() << "Re-making old geometry at old time   : " << old_step_time << std::endl;
             prob->init_custom_terrain(fine_geom,*z_phys_nd[level],old_step_time);
-            init_terrain_grid  (fine_geom,*z_phys_nd[level]);
+            init_terrain_grid  (fine_geom,*z_phys_nd[level], zlevels_stag);
             make_J             (fine_geom,*z_phys_nd[level], *detJ_cc[level]);
 
             if (verbose) Print() << "Making src geometry at old_stage_time:  " << old_stage_time << std::endl;
             prob->init_custom_terrain(fine_geom,*z_phys_nd_src[level],old_stage_time);
-            init_terrain_grid  (fine_geom,*z_phys_nd_src[level]);
+            init_terrain_grid  (fine_geom,*z_phys_nd_src[level], zlevels_stag);
             make_J             (fine_geom,*z_phys_nd_src[level], *detJ_cc_src[level]);
 
             if (verbose) Print() << "Making new geometry at new_stage_time: " << new_stage_time << std::endl;
             prob->init_custom_terrain(fine_geom,*z_phys_nd_new[level],new_stage_time);
-            init_terrain_grid  (fine_geom,*z_phys_nd_new[level]);
+            init_terrain_grid  (fine_geom,*z_phys_nd_new[level], zlevels_stag);
             make_J             (fine_geom,*z_phys_nd_new[level], *detJ_cc_new[level]);
 
             Real inv_dt  = 1./slow_dt;

Source/Utils/TerrainMetrics.H

@@ -10,8 +10,16 @@
  * Utility routines for constructing terrain metric terms
  */
 
-// Declare function for ERF.cpp
-void init_terrain_grid( const amrex::Geometry& geom, amrex::MultiFab& z_phys_nd);
+// Declare functions for ERF.cpp
+void init_zlevels (amrex::Vector<amrex::Real>& zlevels_stag,
+                   const amrex::Geometry& geom,
+                   const amrex::Real grid_stretching_ratio,
+                   const amrex::Real zsurf,
+                   const amrex::Real dz0);
+
+void init_terrain_grid (const amrex::Geometry& geom,
+                        amrex::MultiFab& z_phys_nd,
+                        amrex::Vector<amrex::Real> const& z_levels_h);
 
 //*****************************************************************************************
 // Compute terrain metric terms at cell-center

Source/Utils/TerrainMetrics.cpp

@@ -5,6 +5,59 @@
 
 using namespace amrex;
 
+void
+init_zlevels (amrex::Vector<amrex::Real>& zlevels_stag,
+              const amrex::Geometry& geom,
+              const amrex::Real grid_stretching_ratio,
+              const amrex::Real zsurf,
+              const amrex::Real dz0)
+{
+    auto dx = geom.CellSizeArray();
+    const amrex::Box& domain = geom.Domain();
+    int nz = domain.length(2)+1; // staggered
+    zlevels_stag.resize(nz);
+
+    if (grid_stretching_ratio == 0) {
+        // This is the default for z_levels
+        for (int k = 0; k < nz; k++)
+        {
+            zlevels_stag[k] = k * dx[2];
+        }
+    } else {
+        // Create stretched grid based on initial dz and stretching ratio
+        zlevels_stag[0] = zsurf;
+        amrex::Real dz = dz0;
+        amrex::Print() << "Stretched grid levels: " << zsurf;
+        for (int k = 1; k < nz; k++)
+        {
+            zlevels_stag[k] = zlevels_stag[k-1] + dz;
+            amrex::Print() << " " << zlevels_stag[k];
+            dz *= grid_stretching_ratio;
+        }
+        amrex::Print() << std::endl;
+    }
+
+    // Try reading in terrain_z_levels, which allows arbitrarily spaced grid
+    // levels to be specified and will take precedence over the
+    // grid_stretching_ratio parameter
+    ParmParse pp("erf");
+    int n_zlevels = pp.countval("terrain_z_levels");
+    if (n_zlevels > 0)
+    {
+        if (n_zlevels != nz) {
+            amrex::Print() << "You supplied " << n_zlevels << " staggered terrain_z_levels " << std::endl;
+            amrex::Print() << "but n_cell+1 in the z-direction is " <<  nz << std::endl;
+            amrex::Abort("You must specify a z_level for every value of k");
+        }
+
+        if (grid_stretching_ratio > 0) {
+            amrex::Print() << "Note: Found terrain_z_levels, ignoring grid_stretching_ratio" << std::endl;
+        }
+
+        pp.queryarr("terrain_z_levels", zlevels_stag, 0, nz);
+    }
+}
+
 /**
  * Computation of the terrain grid from BTF, STF, or Sullivan TF model
  *
@@ -14,7 +67,7 @@ using namespace amrex;
  */
 
 void
-init_terrain_grid (const Geometry& geom, MultiFab& z_phys_nd)
+init_terrain_grid (const Geometry& geom, MultiFab& z_phys_nd, amrex::Vector<Real> const& z_levels_h)
 {
   auto dx = geom.CellSizeArray();
   auto ProbHiArr = geom.ProbHiArray();
@@ -24,35 +77,15 @@ init_terrain_grid (const Geometry& geom, MultiFab& z_phys_nd)
   int domlo_x = domain.smallEnd(0); int domhi_x = domain.bigEnd(0) + 1;
   int domlo_y = domain.smallEnd(1); int domhi_y = domain.bigEnd(1) + 1;
   int domlo_z = domain.smallEnd(2); int domhi_z = domain.bigEnd(2) + 1;
+  int nz = domain.length(2)+1;
 
   // User-selected method from inputs file (BTF default)
   ParmParse pp("erf");
   int terrain_smoothing = 0;
   pp.query("terrain_smoothing", terrain_smoothing);
 
   // "custom" refers to terrain analytically specified, such as WitchOfAgnesi
-  std::string terrain_type = "custom";
-
-  int nz = domain.length(2)+1;
-  amrex::Vector<Real> z_levels_h;
-  z_levels_h.resize(nz);
-
-  // This is the default for z_levels
-  for (int k = 0; k < nz; k++)
-  {
-      z_levels_h[k] = k * dx[2];
-  }
-
-  // But we can read them in from the inputs file as an alternative
-  int n_zlevels = pp.countval("terrain_z_levels");
-
-  pp.queryarr("terrain_z_levels", z_levels_h, 0, nz);
-
-  if (n_zlevels > 0 && n_zlevels != nz) {
-      amrex::Print() << "You supplied " << n_zlevels << " z_levels " << std::endl;
-      amrex::Print() << "but n_cell in the z-direction is " <<  nz << std::endl;
-      amrex::Abort("You must specify a z_level for every value of k");
-  }
+  //std::string terrain_type = "custom";
 
    amrex::Gpu::DeviceVector<Real> z_levels_d;
    z_levels_d.resize(nz);

Source/prob_common.H

@@ -107,16 +107,19 @@ public:
     /**
      * Function to perform custom initialization of terrain
      *
+     * Note: Terrain functionality can also be used to provide grid stretching.
+     *
      * @param[in] geom container for geometric information
      * @param[out] z_phys_nd height coordinate at nodes
      * @param[in] time current time
     */
     virtual void
     init_custom_terrain (const amrex::Geometry& /*geom*/,
-                         amrex::MultiFab& /*z_phys_nd*/,
+                         amrex::MultiFab& z_phys_nd,
                          const amrex::Real& /*time*/)
     {
-        amrex::Error("Should never call init_custom_terrain for "+name()+" problem");
+        amrex::Print() << "Initialized flat terrain at z=0" << std::endl;
+        z_phys_nd.setVal(0.0);
     }
 
 #ifdef ERF_USE_TERRAIN_VELOCITY