Merge branch 'development' into pltfile_face_vels

CMake/BuildERFExe.cmake

@@ -169,6 +169,7 @@ function(build_erf_lib erf_lib_name)
        ${SRC_DIR}/SourceTerms/ERF_make_sources.cpp
        ${SRC_DIR}/SourceTerms/ERF_moist_set_rhs.cpp
        ${SRC_DIR}/SourceTerms/ERF_NumericalDiffusion.cpp
+       ${SRC_DIR}/SourceTerms/ERF_ForestDrag.cpp
        ${SRC_DIR}/TimeIntegration/ERF_ComputeTimestep.cpp
        ${SRC_DIR}/TimeIntegration/ERF_Advance.cpp
        ${SRC_DIR}/TimeIntegration/ERF_TimeStep.cpp

Docs/sphinx_doc/Inputs.rst

@@ -1096,6 +1096,9 @@ List of Parameters
 | **erf.input_sounding_file**         | Name(s) of the         | String(s)         | input_sounding_file |
 |                                     | input sounding file(s) |                   |                     |
 +-------------------------------------+------------------------+-------------------+---------------------+
+| **erf.forest_file**                 | Name(s) of the         | String            | None                |
+|                                     | canopy forest file     |                   |                     |
++-------------------------------------+------------------------+-------------------+---------------------+
 | **erf.input_sounding_time**         | Time(s) of the         | Real(s)           | false               |
 |                                     | input sounding file(s) |                   |                     |
 +-------------------------------------+------------------------+-------------------+---------------------+

Docs/sphinx_doc/LinearSolvers.rst

@@ -17,8 +17,9 @@ Multigrid can also be used when the union of grids at a level is not in itself r
 For simulations using grid stretching in the vertical but flat terrain, we must use the hybrid FFT solver in which
 we perform 2D transforms only in the lateral directions and couple the solution in the vertical direction with a tridiagonal solve.
 In both these cases we use a 7-point stencil.
+
 To solve the Poisson equation on terrain-fitted coordinates with general terrain,
 we rely on the FFT-preconditioned GMRES solver since the stencil effectively has variable coefficients and requires 19 points.
 
-   .. note::
-      **Currently only doubly periodic lateral boundary conditions are supported by the hybrid FFT, and therefore by the GMRES solver.  More general boundary conditions are a work in progress.**
+-   .. note::
+-      **The FFT solver / preconditioner can only be used when the union of grids at a level is itself rectangular.

Docs/sphinx_doc/index.rst

@@ -57,6 +57,7 @@ In addition to this documentation, there is API documentation for ERF generated
    theory/Forcings.rst
    Particles.rst
    theory/WindFarmModels.rst
+   theory/Forest.rst
    theory/UnitsAndConstants.rst
 
 .. toctree::

Docs/sphinx_doc/theory/Buoyancy.rst

@@ -7,11 +7,24 @@
 
 .. _sec:Buoyancy:
 
+ERF has several options for how to define the buoyancy force.
+
 Buoyancy
 =========
 
-ERF has three options for how to define the buoyancy force.  Even in the absence of moisture these
-expressions are not equivalent.
+When using the anelastic formulation, the expression for buoyancy used by ERF is
+
+.. math::
+     \mathbf{B} = -\rho^\prime \mathbf{g} \approx -\rho_0 \mathbf{g} ( \frac{\theta_v^\prime}{\overline{\theta_0}}
+
+where we define
+
+.. math::
+     \theta_v = \theta_d (1 - (1 - \frac{R_v}{R_d}) (q_v + q_c) - \frac{R_v}{R_d} q_c)
+
+
+When using the compressible formulation, the following choices are available.
+
 
 Density of the mixture
 -----------------------
@@ -52,56 +65,10 @@ for the background state. For eg., a usual scenario is that of a background stat
 the total background density :math:`\rho_0 = \rho_{d_0}(1 + q_{v_0})`, where :math:`\rho_{d_0}` and :math:`q_{v_0}` are the background dry density and vapor mixing ratio respectively.
 As a check, we observe that :math:`\rho^\prime_0 = 0`, which means that the background state is not buoyant.
 
-Type 2
-------
-
-The second option for the buoyancy force is
-
-.. math::
-   \mathbf{B} = -\rho_0 \mathbf{g} ( 0.61 q_v^\prime - q_c^\prime - q_i^\prime - q_p^\prime
-                  + \frac{T^\prime}{\bar{T}} (1.0 + 0.61 \bar{q_v} - \bar{q_i} - \bar{q_c} - \bar{q_p}) )
-
-To derive this expression, we define :math:`T_v = T (1 + 0.61 q_v − q_c − q_i - q_p)`, then we can write
-
-.. math::
-    p = \rho (R_d q_d + R_v q_v) T = \rho R_d T (1 + 0.61 q_v − q_c − q_i - q_p ) = \rho R_d T_v
-
-
-Starting from :math:`p = \rho R_d T_v` and neglecting :math:`\frac{p^\prime}{\bar{p}}`, we now write
-
-.. math::
-   \frac{\rho^\prime}{\overline{\rho}} = -\frac{T_v^\prime}{\overline{T_v}}
-
-and define
-
-.. math::
-
-   T_v^\prime = T_v - \overline{T_v} \approx \overline{T} [ 0.61 q_v^\prime - (q_c^\prime + q_i^\prime + q_p^\prime)] +
-               (T - \overline{T}) [1+ 0.61 \bar{q_v} - \bar{q_c} - \bar{q_i} - \bar{q_p} ] .
-
-where we have retained only first order terms in perturbational quantities.
-
-Then
-
-.. math::
-
-   \mathbf{B} = \rho^\prime \mathbf{g} = -\overline{\rho} \frac{\overline{T}}{\overline{T_v}} \mathbf{g} [ 0.61 q_v^\prime - q_c^\prime - q_i^\prime - q_p^\prime ) + \frac{T^\prime}{\overline{T_v}} (1.0 + 0.61 \bar{q_v} - \bar{q_i} - \bar{q_c} - \bar{q_p}) ]
-
-where the overbar represents a horizontal average of the current state and the perturbation is defined relative to that average.
-
-Again keeping only the first order terms in the mass mixing ratios, we can simplify this to
-
-.. math::
-   \mathbf{B} = \rho^\prime \mathbf{g} = -\rho_0 \mathbf{g} [ 0.61 q_v^\prime - q_c^\prime + q_i^\prime + q_p^\prime
-                  + \frac{T^\prime}{\overline{T}} (1.0 + 0.61 \bar{q_v} - \bar{q_i} - \bar{q_c} - \bar{q_p}) ]
-
-We note that this reduces to Type 3 if the horizontal averages of the moisture terms are all zero.
-
 Type 3
 ------
 
-The third formulation of the buoyancy term assumes that the horizontal averages of the moisture quantities are negligible,
-which removes the need to compute horizontal averages of these quantities.   This reduces the Type 2 expression to the following:
+This formulation of the buoyancy term assumes that the horizontal averages of the moisture quantities are negligible.
 
 .. math::
      \mathbf{B} = \rho^\prime \mathbf{g} \approx -\rho_0 \mathbf{g} ( \frac{T^\prime}{\overline{T}}
@@ -120,7 +87,7 @@ This expression for buoyancy is from `khairoutdinov2003cloud`_ and `bryan2002ben
 .. math::
 
     \begin{equation}
-    \mathbf{B} = \rho'\mathbf{g} \approx -\rho\Bigg(\frac{T'}{T} + 0.61 q_v' - q_c - q_p - \frac{p'}{p}\Bigg)\mathbf{g},
+    \mathbf{B} = \rho'\mathbf{g} \approx -\rho \mathbf{g} \Bigg(\frac{T'}{T} + 0.61 q_v' - q_c - q_p - \frac{p'}{p}\Bigg)
     \end{equation}
 
 The derivation follows. The total density is given by :math:`\rho = \rho_d(1 + q_v + q_c + q_p)`, which can be written as

Docs/sphinx_doc/theory/Forest.rst

@@ -0,0 +1,31 @@
+Forest Model
+--------------
+
+The forest model provides an option to include the drag from forested regions to be included in the momentum equation. The
+drag force is calculated as follows:
+
+.. math::
+
+   F_i= - C_d L(x,y,z) U_i | U_i |
+
+
+Here :math:`C_d` is the coefficient of drag for the forested region and :math:`L(x,y,z)` is the leaf area density (LAD) for the
+forested region. A three-dimensional model for the LAD is usually unavailable and is also cumbersome to use if there are thousands
+of trees. Two different models are available as an alternative:
+
+.. math::
+   L=\frac{LAI}{h}
+
+.. math::
+   L(z)=L_m \left(\frac{h - z_m}{h - z}\right)^n  exp\left[n \left(1 -\frac{h - z_m}{h - z}\right )\right]
+
+Here :math:`LAI` is the leaf area index and is available from measurements, :math:`h` is the height of the tree, :math:`z_m` is the location
+of the maximum LAD, :math:`L_m` is the maximum value of LAD at :math:`z_m` and :math:`n` is a model constant with values  6 (below :math:`z_m`) and 0.5
+(above :math:`z_m`), respectively. :math:`L_m` is computed by integrating the following equation (see `Lalic and Mihailovic (2004)
+<https://doi.org/10.1175/1520-0450(2004)043<0641:AERDLD>2.0.CO;2>`_):
+
+.. math::
+   LAI = \int_{0}^{h} L(z) dz
+
+The simplified model with uniform LAD is recommended for forested regions with no knowledge of the individual trees. LAI values can be used from
+climate model look-up tables for different regions around the world if no local remote sensing data is available.

Exec/ABL/erf_forest_def

@@ -0,0 +1,4 @@
+1  1024 512 45 200  0.2 6 0.8
+1  1024 1024 35 200  0.2 6 0.8
+1  1024 1224 75 200  0.2 6 0.8 
+2  1024 1524 120 200 0.2 10 0.8

Exec/ABL/inputs_canopy

@@ -0,0 +1,64 @@
+# ------------------  INPUTS TO MAIN PROGRAM  -------------------
+max_step = 50
+
+amrex.fpe_trap_invalid = 1
+
+fabarray.mfiter_tile_size = 1024 1024 1024
+
+# PROBLEM SIZE & GEOMETRY
+geometry.prob_extent =  2048     2048    1024
+amr.n_cell           =    64       64      32
+
+geometry.is_periodic = 1 1 0
+
+zlo.type = "SlipWall"
+zhi.type = "SlipWall"
+
+# TIME STEP CONTROL
+erf.fixed_dt = 2.0  # fixed time step depending on grid resolution
+erf.cfl      = 0.5
+
+# 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       = 50         # number of timesteps between checkpoints
+
+# PLOTFILES
+erf.plot_file_1     = plt        # prefix of plotfile name
+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"
+erf.les_type        = "Smagorinsky"
+erf.Cs              = 0.1
+
+erf.init_type = "uniform"
+
+#erf.forest_file = erf_forest_def
+
+# 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
+
+# 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.W_0_Pert_Mag = 0.0

Source/BoundaryConditions/ERF_ABLMost.H

@@ -431,15 +431,15 @@ public:
     get_t_surf (const int& lev) { return t_surf[lev].get(); }
 
     amrex::Real
-    get_zref () {return m_ma.get_zref();}
+    get_zref () { return m_ma.get_zref(); }
 
     const amrex::FArrayBox*
-    get_z0 (const int& lev) {return &z_0[lev];}
+    get_z0 (const int& lev) { return &z_0[lev]; }
 
     bool have_variable_sea_roughness() { return m_var_z0; }
 
     const amrex::iMultiFab*
-    get_lmask(const int& lev) {return m_lmask_lev[lev][0];}
+    get_lmask(const int& lev) { return m_lmask_lev[lev][0]; }
 
     int
     lmask_min_reduce (amrex::iMultiFab& lmask, const int& nghost)

Source/DataStructs/ERF_DataStruct.H

@@ -126,14 +126,15 @@ struct SolverChoice {
             if (moisture_type == MoistureType::Kessler_NoRain ||
                 moisture_type == MoistureType::SAM ||
                 moisture_type == MoistureType::SAM_NoIce ||
-                moisture_type == MoistureType::SAM_NoPrecip_NoIce) {
-                buoyancy_type = 1; // asserted in make buoyancy
+                moisture_type == MoistureType::SAM_NoPrecip_NoIce)
+            {
+                buoyancy_type = 1; // uses Rhoprime
             } else {
                 buoyancy_type = 2; // uses Tprime
             }
         }
 
-        // Which expression (1,2 or 3) to use for buoyancy
+        // Which expression (1,2/3 or 4) to use for buoyancy
         pp.query("buoyancy_type", buoyancy_type);
 
         // What type of land surface model to use
@@ -196,7 +197,7 @@ struct SolverChoice {
         if (any_anelastic == 1) {
             constant_density = true;
             project_initial_velocity = true;
-            buoyancy_type = 2; // uses Tprime
+            buoyancy_type = 3; // (This isn't actually used when anelastic is set)
         } else {
             pp.query("project_initial_velocity", project_initial_velocity);
 
@@ -644,5 +645,8 @@ struct SolverChoice {
     amrex::Real turb_disk_angle = -1.0;
     amrex::Real windfarm_x_shift = -1.0;
     amrex::Real windfarm_y_shift = -1.0;
+
+    // Flag for valid canopy model
+    bool do_forest {false};
 };
 #endif

Source/ERF.H

@@ -45,6 +45,7 @@
 #include <ERF_PhysBCFunct.H>
 #include <ERF_FillPatcher.H>
 #include <ERF_SampleData.H>
+#include <ERF_ForestDrag.H>
 
 #ifdef ERF_USE_PARTICLES
 #include "ERF_ParticleData.H"
@@ -1124,6 +1125,7 @@ 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;
 
     //
     // Holds info for dynamically generated tagging criteria

Source/ERF.cpp

@@ -30,7 +30,7 @@ Real ERF::change_max    =  1.1;
 int  ERF::fixed_mri_dt_ratio = 0;
 
 // Dictate verbosity in screen output
-int ERF::verbose        = 0;
+int  ERF::verbose       = 0;
 int  ERF::mg_verbose    = 0;
 bool ERF::use_fft       = false;
 
@@ -129,6 +129,10 @@ ERF::ERF_shared ()
     lsm_data.resize(nlevs_max);
     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; }
+
     ReadParameters();
     initializeMicrophysics(nlevs_max);
 
@@ -300,6 +304,7 @@ ERF::ERF_shared ()
         Hwave_onegrid[lev] = nullptr;
         Lwave_onegrid[lev] = nullptr;
     }
+
     // Theta prim for MOST
     Theta_prim.resize(nlevs_max);
 
@@ -1590,6 +1595,15 @@ ERF::ReadParameters ()
         pp.query("cf_width", cf_width);
         pp.query("cf_set_width", cf_set_width);
 
+        // Query the canopy model file name
+        std::string forestfile;
+        solverChoice.do_forest = pp.query("forest_file", forestfile);
+        if (solverChoice.do_forest) {
+            for (int lev = 0; lev <= max_level; ++lev) {
+                m_forest[lev] = std::make_unique<ForestDrag>(forestfile);
+            }
+        }
+
         // AmrMesh iterate on grids?
         bool iterate(true);
         pp_amr.query("iterate_grids",iterate);

Source/ERF_make_new_level.cpp

@@ -131,6 +131,11 @@ 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()); }
+
     //********************************************************************************************
     // Microphysics
     // *******************************************************************************************
@@ -205,6 +210,11 @@ 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()); }
+
     //********************************************************************************************
     // Microphysics
     // *******************************************************************************************
@@ -329,6 +339,11 @@ 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()); }
+
     // *****************************************************************************************************
     // Create the physbcs objects (after initializing the terrain but before calling FillCoarsePatch
     // *****************************************************************************************************