Cleanup inputs, input files, docs (#2025)

* Remove QKE from list of plotfile outputs * Add note about eddy viscosity * Add eddy viscosity (nut) output * Cleanup comments in input files * Remove deprecated rayleigh option * Update numerical diffusion input param, enable num diff for momenta * Rename erf.dynamicViscosity to erf.dynamic_viscosity for consistent input naming * Strip out erf.rho0_trans = 1 (default) from input files * Clarify comments * Put realistic Deardorff coeffs (follow Moeng 1984) * Cleanup prob params for initial fields * Remove redundant erf.KE_0 initialiation; can already accomplish this with prob.KE_0 * Remove some bogus inputs * Removed specified alpha for molec_diff_type=None, only read diffusivities if not molec_diff_type==None * Update best practices page * Add 2D to best practices page * Cleanup/clarify screen output * Merge oops * Forgot to update test inputs * More fixes to test files * Add prob.rhoKE_0 init option for ABL, other fixes to not break testing * Comment out unused vars --------- Co-authored-by: Aaron M. Lattanzi <[email protected]>
erf-model · Dec 18, 2024 · 336c0b3 · 336c0b3
1 parent f445446
commit 336c0b3
Show file tree

Hide file tree

Showing 169 changed files with 664 additions and 812 deletions.
diff --git a/Docs/sphinx_doc/BestPractices.rst b/Docs/sphinx_doc/BestPractices.rst
@@ -45,41 +45,61 @@ Large-Eddy Simulations
 * Time Integration
 
   - Split timestepping offers some computational cost savings but still does
-    not allow you to run with an incompressible/anelastic time-step size.
-  - The acoustic CFL should be less than 0.5, with 4--6 fast timesteps
-    (substeps) according to WRF best practices.
+    not allow you to run with an incompressible/anelastic time-step size in
+    general.
+  - The acoustic CFL should conservatively be less than or equal to 0.5, with
+    4--6 fast timesteps (substeps) according to WRF best practices. If not
+    explicitly specified (through ``erf.fixed_mri_dt_ratio`` or
+    ``erf.fixed_fast_dt``, the number of substeps in ERF is chosen based on the
+    same algorithm as WRF. If the user follows the recommendation that
+    dt [s] ~ 6 dx [km],
+    then 4 substeps will be used, giving an effective CFL of approximately 0.5.
+    This meets the stability criteria from Wicker & Skamarock 2002 that, for a
+    5th-order scheme, the CFL be less than 1.42/sqrt(3) = 0.820.
 
     .. code-block:: python
 
        erf.fixed_dt           = 0.06  # slow timestep
 
        # These are equivalent and result in a fixed fast timestep size
-       #   if dx=10, speed of sound ~ 350 m/s
-       erf.fixed_fast_dt      = 0.01  # ==> CFL~0.35
-       #erf.fixed_mri_dt_ratio = 6
-
-       # Alternatively, let ERF chose the fast timestep
-       #erf.cfl                = 0.5
-
-  - We note that ERF LESs with up to 10 fast timesteps have successfully been
+       #   if dx=10, speed of sound ~ 300 m/s
+       erf.fixed_mri_dt_ratio = 4
+       #   or
+       #erf.fixed_fast_dt      = 0.015  # ==> CFL~0.45
+       #   or, let ERF chose the fast timestep
+       #erf.substepping_cfl    = 0.5
+
+  - Following the WRF guidelines for dt is conservative. More aggressive time
+    integration--i.e., larger time steps with more substeps--is possible. We
+    note that ERF LESs with 10 or more fast timesteps have successfully been
     run but your mileage may vary.
 
 
 Single-Column Model
 -------------------
 
-* Currently, ERF does not have the ability to run a true single-column model
-  (SCM). The grid size in the lateral directions must have a minimum number of
-  cells. This will give comparable results, e.g.:
+* An SCM is set up with a single cell in the lateral directions:
 
   .. code-block:: python
 
      geometry.prob_extent = 400  400  400
-     amr.n_cell           =   2    2   64
+     amr.n_cell           =   1    1   64
      geometry.is_periodic =   1    1    0
 
-  When set up this way, the solution is not sensitive to horizontal problem
-  extent.
-
 * An SCM was successfully run with third-order advection in the horizontal and
   vertical.
+
+
+2-D Cases
+---------
+
+* A 2-D planar domain can be configured as follows:
+
+  .. code-block:: python
+
+     geometry.prob_extent = 10000  100  1000
+     amr.n_cell           =   100    1    20
+     geometry.is_periodic =     0    0     0
+
+     ylo.type = "SlipWall"
+     yhi.type = "SlipWall"
diff --git a/Docs/sphinx_doc/Inputs.rst b/Docs/sphinx_doc/Inputs.rst
@@ -906,7 +906,7 @@ List of Parameters
 |                                  | viscosity and      | "Constant", or      |              |
 |                                  | diffusivity?       | "ConstantAlpha"     |              |
 +----------------------------------+--------------------+---------------------+--------------+
-| **erf.dynamicViscosity**         | Viscous coeff. if  | Real                | 0.0          |
+| **erf.dynamic_viscosity**        | Viscous coeff. if  | Real                | 0.0          |
 |                                  | DNS                |                     |              |
 +----------------------------------+--------------------+---------------------+--------------+
 | **erf.Cs**                       | Constant           | Real                | 0.0          |
@@ -918,13 +918,11 @@ List of Parameters
 | **erf.Sc_t**                     | Turbulent Schmidt  | Real                | 1.0          |
 |                                  | Number             |                     |              |
 +----------------------------------+--------------------+---------------------+--------------+
-| **erf.use_NumDiff**              | Use 6th order      | "true",             | "false"      |
-|                                  | numerical diffusion| "false"             |              |
-|                                  |                    |                     |              |
-+----------------------------------+--------------------+---------------------+--------------+
-| **erf.NumDiffCoeff**             | Coefficient for    | Real                | 0.0          |
-|                                  | 6th order          | [0.0,  1.0]         |              |
-|                                  | numerical diffusion|                     |              |
+| **erf.num_diff_coeff**           | Coefficient for    | Real                | 0.0          |
+|                                  | 6th-order          | [0.0,  1.0]         |              |
+|                                  | numerical          |                     |              |
+|                                  | diffusion, set to 0|                     |              |
+|                                  | to disable         |                     |              |
 +----------------------------------+--------------------+---------------------+--------------+
 
 Note: in the equations for the evolution of momentum, potential temperature and advected scalars, the

diff --git a/Docs/sphinx_doc/Plotfiles.rst b/Docs/sphinx_doc/Plotfiles.rst
@@ -187,19 +187,13 @@ Output Options
 |                             |                  |
 |                             |                  |
 +-----------------------------+------------------+
-| **KE**                      | Kinetic energy   |
-|                             |                  |
-|                             |                  |
-+-----------------------------+------------------+
-| **QKE**                     | Turbulent        |
+| **KE**                      | SGS turbulent    |
 |                             | kinetic energy   |
-|                             | * 2              |
-+-----------------------------+------------------+
-| **rhoKE**                   | Density * KE     |
-|                             |                  |
+|                             | (from Deardorff  |
+|                             |  or MYNN)        |
 |                             |                  |
 +-----------------------------+------------------+
-| **rhoQKE**                  | Density * QKE    |
+| **rhoKE**                   | Density * KE     |
 |                             |                  |
 |                             |                  |
 +-----------------------------+------------------+
@@ -263,13 +257,22 @@ Output Options
 |                             | mass points      |
 |                             |                  |
 +-----------------------------+------------------+
+| **nut**                     | Eddy viscosity,  |
+|                             | nu_t             |
++-----------------------------+------------------+
 | **Kmv**                     | Vertical         |
 |                             | Eddy Diffusivity |
 |                             | of Momentum      |
 +-----------------------------+------------------+
 | **Kmh**                     | Horizontal       |
 |                             | Eddy Diffusivity |
 |                             | of Momentum      |
+|                             | (Note: For LES,  |
+|                             | this is the      |
+|                             | _dynamic_ eddy   |
+|                             | viscosity, mu_t  |
+|                             | = rho * nu_t     |
+|                             | and Kmh==Kmv)    |
 +-----------------------------+------------------+
 | **Khv**                     | Vertical         |
 |                             | Eddy Diffusivity |

diff --git a/Exec/ABL/ERF_Prob.H b/Exec/ABL/ERF_Prob.H
@@ -11,7 +11,8 @@ struct ProbParm : ProbParmDefaults {
   amrex::Real rho_0 = 0.0;
   amrex::Real T_0   = 0.0;
   amrex::Real A_0   = 1.0;
-  amrex::Real KE_0  = 0.1;
+  amrex::Real KE_0  = 0.1;  // initialize to calculated rho_hse times input KE
+  amrex::Real rhoKE_0 = -1; // alternatively, initialize to specified (rho*KE)
 
   amrex::Real KE_decay_height = -1;
   amrex::Real KE_decay_order  =  1;

diff --git a/Exec/ABL/ERF_Prob.cpp b/Exec/ABL/ERF_Prob.cpp
@@ -17,6 +17,7 @@ Problem::Problem(const amrex::Real* problo, const amrex::Real* probhi)
   pp.query("T_0", parms.T_0);
   pp.query("A_0", parms.A_0);
   pp.query("KE_0", parms.KE_0);
+  pp.query("rhoKE_0", parms.rhoKE_0);
   pp.query("KE_decay_height", parms.KE_decay_height);
   pp.query("KE_decay_order", parms.KE_decay_order);
 
@@ -126,7 +127,11 @@ Problem::init_custom_pert(
     // Set an initial value for SGS KE
     if (state_pert.nComp() > RhoKE_comp) {
         // Deardorff
-        state_pert(i, j, k, RhoKE_comp) = r_hse(i,j,k) * parms_d.KE_0;
+        if (parms_d.rhoKE_0 > 0) {
+            state_pert(i, j, k, RhoKE_comp) = parms_d.rhoKE_0;
+        } else {
+            state_pert(i, j, k, RhoKE_comp) = r_hse(i,j,k) * parms_d.KE_0;
+        }
         if (parms_d.KE_decay_height > 0) {
             // scale initial SGS kinetic energy with height
             state_pert(i, j, k, RhoKE_comp) *= max(

diff --git a/Exec/ABL/inputs.read b/Exec/ABL/inputs.read
@@ -45,8 +45,6 @@ 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"

diff --git a/Exec/ABL/inputs.write b/Exec/ABL/inputs.write
@@ -38,8 +38,6 @@ erf.plot_int_1      = 10         # 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"

diff --git a/Exec/ABL/inputs_DataSampler b/Exec/ABL/inputs_DataSampler
@@ -48,8 +48,6 @@ erf.plot_int_1      = 10         # 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"

diff --git a/Exec/ABL/inputs_GABLS1_deardorff b/Exec/ABL/inputs_GABLS1_deardorff
@@ -76,7 +76,6 @@ erf.les_type = "Deardorff"
 erf.Ck       = 0.1  # Coefficient in Moeng1984, Eqn. 19
 erf.Ce       = 0.7  # Note: Ce_lcoeff = C_e - 1.9*C_k
 erf.Ce_wall  = 3.9  # To account for "wall effects"
-erf.sigma_k  = 0.5  # TKE diffusion coeff = 2*Km = Km * inv_sigma_k, Moeng1984, Eqn. 15
 
 erf.rayleigh_damp_U = false
 erf.rayleigh_damp_V = false

diff --git a/Exec/ABL/inputs_canopy b/Exec/ABL/inputs_canopy
@@ -36,8 +36,6 @@ erf.plot_int_1      = 50        # 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"

diff --git a/Exec/ABL/inputs_deardorff b/Exec/ABL/inputs_deardorff
@@ -35,30 +35,27 @@ erf.plot_int_1        = 10        # number of timesteps between plotfiles
 erf.plot_vars_1       = density rhoKE 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 = "Deardorff"
-erf.Ck       = 0.1
-erf.sigma_k  = 1.0
-erf.Ce       = 0.1
+erf.Ck       = 0.1  # Coefficient in Moeng1984, Eqn. 19
+erf.Ce       = 0.7  # Note: Ce_lcoeff = C_e - 1.9*C_k
+erf.Ce_wall  = 3.9  # To account for "wall effects"
 
 erf.init_type = "uniform"
-erf.KE_0 = 0.1 # for Deardorff
 
-# PROBLEM PARAMETERS
+# PROBLEM PARAMETERS -- set these for uniform init only
 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
+prob.A_0   = 1.0 # advected scalar
+prob.U_0   = 10.0
+prob.V_0   = 0.0
+prob.W_0   = 0.0
+prob.T_0   = 300.0
+prob.KE_0  = 0.1
 
 # Higher values of perturbations lead to instability
 # Instability seems to be coming from BC
 prob.U_0_Pert_Mag = 0.08
-prob.V_0_Pert_Mag = 0.08 #
+prob.V_0_Pert_Mag = 0.08
 prob.W_0_Pert_Mag = 0.0
diff --git a/Exec/ABL/inputs_deardorff_msf b/Exec/ABL/inputs_deardorff_msf
@@ -38,31 +38,27 @@ erf.plot_int_1        = 1        # number of timesteps between plotfiles
 erf.plot_vars_1       = density rhoKE rhoadv_0 x_velocity y_velocity z_velocity pressure temp theta scalar
 
 # SOLVER CHOICE
-erf.alpha_T = 1.0
-erf.alpha_C = 1.0
 erf.use_gravity = false
 
-erf.molec_diff_type = "Constant"
+erf.molec_diff_type = "None"
 erf.les_type = "Deardorff"
-erf.rho0_trans = 1.0 
-erf.dynamicViscosity = 75.0 
-erf.Ck       = 0.1
-erf.sigma_k  = 1.0
-erf.Ce       = 0.1
+erf.Ck       = 0.1  # Coefficient in Moeng1984, Eqn. 19
+erf.Ce       = 0.7  # Note: Ce_lcoeff = C_e - 1.9*C_k
+erf.Ce_wall  = 3.9  # To account for "wall effects"
 
 erf.init_type = "uniform"
-erf.KE_0 = 0.1 # for Deardorff
 
-# PROBLEM PARAMETERS
+# PROBLEM PARAMETERS -- set these for uniform init only
 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.A_0   = 1.0 # advected scalar
+prob.U_0   = 10.0
+prob.V_0   = 0.0
+prob.W_0   = 0.0
+prob.T_0   = 300.0
+prob.KE_0  = 0.1
 
 # Higher values of perturbations lead to instability
 # Instability seems to be coming from BC
 prob.U_0_Pert_Mag = 0.08
-prob.V_0_Pert_Mag = 0.08 #
+prob.V_0_Pert_Mag = 0.08
 prob.W_0_Pert_Mag = 0.0
diff --git a/Exec/ABL/inputs_deardorff_no_msf b/Exec/ABL/inputs_deardorff_no_msf
@@ -37,31 +37,27 @@ erf.plot_int_1        = 1         # number of timesteps between plotfiles
 erf.plot_vars_1       = density rhoKE rhoadv_0 x_velocity y_velocity z_velocity pressure temp theta scalar
 
 # SOLVER CHOICE
-erf.alpha_T = 1.0
-erf.alpha_C = 1.0
 erf.use_gravity = false
 
-erf.molec_diff_type = "Constant"
+erf.molec_diff_type = "None"
 erf.les_type = "Deardorff"
-erf.rho0_trans = 1.0 
-erf.dynamicViscosity = 75.0 
-erf.Ck       = 0.1
-erf.sigma_k  = 1.0
-erf.Ce       = 0.1
+erf.Ck       = 0.1  # Coefficient in Moeng1984, Eqn. 19
+erf.Ce       = 0.7  # Note: Ce_lcoeff = C_e - 1.9*C_k
+erf.Ce_wall  = 3.9  # To account for "wall effects"
 
 erf.init_type = "uniform"
-erf.KE_0 = 0.1 # for Deardorff
 
-# PROBLEM PARAMETERS
+# PROBLEM PARAMETERS -- set these for uniform init only
 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.A_0   = 1.0 # advected scalar
+prob.U_0   = 10.0
+prob.V_0   = 0.0
+prob.W_0   = 0.0
+prob.T_0   = 300.0
+prob.KE_0  = 0.1
 
 # Higher values of perturbations lead to instability
 # Instability seems to be coming from BC
 prob.U_0_Pert_Mag = 0.08
-prob.V_0_Pert_Mag = 0.08 #
+prob.V_0_Pert_Mag = 0.08
 prob.W_0_Pert_Mag = 0.0
diff --git a/Exec/ABL/inputs_most b/Exec/ABL/inputs_most
@@ -39,8 +39,6 @@ erf.plot_int_1      = 10        # 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"

diff --git a/Exec/ABL/inputs_most_bomex b/Exec/ABL/inputs_most_bomex
@@ -43,8 +43,6 @@ erf.plot_int_1      = 10        # number of timesteps between plotfiles
 erf.plot_vars_1     = density rhoadv_0 x_velocity y_velocity z_velocity pressure temp theta qv qc qp
 
 # SOLVER CHOICE
-erf.alpha_T = 0.0
-erf.alpha_C = 1.0
 erf.use_gravity = false
 
 erf.molec_diff_type = "None"

diff --git a/Exec/ABL/inputs_most_pbl b/Exec/ABL/inputs_most_pbl
@@ -42,8 +42,6 @@ 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"