Merge branch 'development' into bomex

Docs/sphinx_doc/MOST.rst

@@ -149,13 +149,13 @@ In ERF, when the MOST boundary condition is applied, velocity and temperature in
 
 MOST Inputs
 ~~~~~~~~~~~~~~~~~~~
-To evaluate the fluxes with MOST, the surface rougness parameter :math:`z_{0}` must be specified. This quantity may be considered a constant or may be parameterized through the friction velocity :math:`u_{\star}`. ERF supports three methods for parameterizing the surface roughness: ``constant``, ``charnock``, and ``modified_charnock``. The latter two methods parameterize :math:`z_{0} = f(u_{\star})` and are described in `Jimenez & Dudhia, American Meteorological Society, 2018 <https://doi.org/10.1175/JAMC-D-17-0137.1>`_. The rougness calculation method may be specified with
+To evaluate the fluxes with MOST, the surface rougness parameter :math:`z_{0}` must be specified. This quantity may be considered a constant or may be parameterized through the friction velocity :math:`u_{\star}`. ERF supports four methods for parameterizing the surface roughness: ``constant``, ``charnock``, ``modified_charnock``, and ``wave_coupled``. The latter three methods parameterize :math:`z_{0} = f(u_{\star})` and are described in `Jimenez & Dudhia, American Meteorological Society, 2018 <https://doi.org/10.1175/JAMC-D-17-0137.1>`_ and `Warner et. al, Ocean Modelling, 2010 <https://doi.org/10.1016/j.ocemod.2010.07.010>`_. The rougness calculation method may be specified with
 
 ::
 
-   erf.most.roughness_type    = STRING    #Z_0 type (constant, charnock, modified_charnock)
+   erf.most.roughness_type    = STRING    #Z_0 type (constant, charnock, modified_charnock, wave_couples)
 
-If the ``charnock`` method is employed, the :math:`a` constant may be specified with ``erf.most.charnock_constant`` (defaults to 0.0185). If the ``modified_charnock`` method is employed, the depth :math:`d` may be specified with ``erf.most.modified_charnock_depth`` (defaults to 30 m).
+If the ``charnock`` method is employed, the :math:`a` constant may be specified with ``erf.most.charnock_constant`` (defaults to 0.0185). If the ``modified_charnock`` method is employed, the depth :math:`d` may be specified with ``erf.most.modified_charnock_depth`` (defaults to 30 m). If the ``wave_coupled`` method is employed, the user must provide wave height and mean wavelength data.
 
 When computing an average :math:`\overline{\phi}` for the MOST boundary, where :math:`\phi` denotes a generic variable, ERF supports a variety of approaches. Specifically, ``planar averages`` and ``local region averages`` may be computed with or without ``time averaging``. With each averaging methodology, the query point :math:`z` may be determined from the following procedures: specified vertical distance :math:`z_{ref}` from the bottom surface, specified :math:`k_{index}`, or (when employing terrain-fit coordinates) specified normal vector length :math:`z_{ref}`. The available inputs to the MOST boundary and their associated data types are
 

Docs/sphinx_doc/building.rst

@@ -263,7 +263,7 @@ For Perlmutter at NERSC, look at the general instructions for building ERF using
    module load PrgEnv-gnu
    module load cudatoolkit
 
-Then build ERF as, for example (specify your own path to the AMReX submodule in `ERF/Submodules/AMReX`):
+Then build ERF as, for example (specify your own path to the AMReX submodule in ``ERF/Submodules/AMReX``):
 
 ::
 
@@ -308,5 +308,47 @@ Finally, you can prepare your SLURM job script, using the following as a guide:
                ./ERF3d.gnu.MPI.CUDA.ex inputs_wrf_baseline max_step=100 ${GPU_AWARE_MPI}" \
              > test.out
 
-To submit your job script, do `sbatch [your job script]` and you can check its status by doing `squeue -u [your username]`.
+To submit your job script, do ``sbatch [your job script]`` and you can check its status by doing ``squeue -u [your username]``.
 
+
+Kestrel (NREL)
+~~~~~~~~~~~~~~
+
+The `Kestrel <https://nrel.github.io/HPC/Documentation/Systems/Kestrel/>`_ cluster is an HPE Cray machine
+composed primarily of CPU compute nodes with 104 core
+Intel Xeon Sapphire Rapids nodes. It also contains a GPU partition with 4 Nvidia H100 GPUs per node.
+
+As with Perlmutter, the GNU Make build system is preferred. To compile and run on CPUs, the default modules
+loaded when logging into Kestrel can be used. If you are unsure about your environment, you can reset to
+the default modules: ::
+
+  module restore
+
+Then, build ERF using the cray compilers (if wishing to use other compilers, you can swap the ``PrgEnv-cray`` module
+for another module as appropriate, see Kestrel user documentation for more details): ::
+
+  make realclean; make -j COMP=cray
+
+For compiling and running on GPUs, the following commands can be used to set up your environment: ::
+
+  module restore; module load PrgEnv-gnu/8.5.0; module load cray-libsci/23.05.1.4; module load cmake; module load cuda/12.3; module load cray-mpich/8.1.28; module load craype/2.7.30;
+
+And then compile: ::
+
+  make realclean; make -j COMP=gnu USE_CUDA=TRUE
+
+When running on Kestrel, GPU node hours are charged allocation units (AUs) at 10 times the rate of CPU node hours.
+For ERF, the performance running on a Kestrel GPU node with 4 GPUs is typically 10-20x running on a CPU node
+with 96-104 MPI ranks per node, so the performance gain from on on GPUs is likely worth the higher charge
+rate for node hours, in addition to providing faster time to solution. However, for smaller problem sizes,
+or problems distributed across too many nodes (resulting in fewer than around 1 million cells/GPU),
+the compute capability of the GPUs may be unsaturated and the performance gain from running on GPUs
+may not justify the higher AU charge. The trade-off is problem dependent, so users may wish to assess
+performance for their particular case and objectives in terms of wall time, AUs used, etc to determine the
+optimal strategy if running large jobs.
+
+Another note about using Kestrel is that partial node allocations are possible, which means the full memory
+available on each node may not be assigned by default. In general, using the ``--exclusive`` flag when
+requesting nodes through the slurm scheduler, which will allocate entire nodes exlcusively for your request,
+is recommended. Otherwise, memory intensive operations such as CUDA compilation may fail. You can alternatively
+request a particular amount of memory with the ``--mem=XXX`` or ``--mem-per-cpu=XXX`` slurm inputs.

Docs/sphinx_doc/theory/WindFarmModels.rst

@@ -4,7 +4,9 @@ Wind farm models
 Introduction
 -------------
 
-ERF supports models for wind farm parametrization in which the effects of wind turbines are represented by imposing a momentum sink on the mean flow and/or turbulent kinetic energy (TKE). Currently only the Fitch model (`Fitch et al. 2012`_) is supported.
+ERF supports models for wind farm parametrization in which the effects of wind turbines are represented by imposing a momentum sink on the mean flow and/or turbulent kinetic energy (TKE). Currently only the Fitch model (`Fitch et al. 2012`_) and Explicit Wake Parametrization (EWP) model (`Volker et al. 2015`_) are supported.
+
+.. _Fitch model:
 
 Fitch model
 ------------
@@ -17,26 +19,27 @@ The Fitch model for wind farms introduced in `Fitch et al. 2012`_  models the ef
 
     \frac{\partial u_{ijk}}{\partial t} &= \frac{u_{ijk}}{|V|_{ijk}}\frac{\partial |V|_{ijk}}{\partial t} \\
     \frac{\partial v_{ijk}}{\partial t} &= \frac{v_{ijk}}{|V|_{ijk}}\frac{\partial |V|_{ijk}}{\partial t} \\
-    \frac{\partial \text{TKE}_{ijk}}{\partial t} &=  \frac{0.5N^{ij}C_{TKE}(|V|_{ijk})|V|_{ijk}^3A_{ijk}}{z_{k+1}-z_k}
+    \frac{\partial \text{TKE}_{ijk}}{\partial t} &=  \frac{0.5N_{ij}C_{TKE}(|V|_{ijk})|V|_{ijk}^3A_{ijk}}{\Delta x \Delta y (z_{k+1}-z_k)}
 
 where
 
 .. math::
 
-    \frac{\partial |V|_{ijk}}{\partial t} = \frac{0.5N^{ij}C_T(|V|_{ijk})|V|_{ijk}^2A_{ijk}}{z_{k+1}-z_k}
+    \frac{\partial |V|_{ijk}}{\partial t} = \frac{0.5N_{ij}C_T(|V|_{ijk})|V|_{ijk}^2A_{ijk}}{\Delta x \Delta y (z_{k+1}-z_k)}
+
+where `u` and `v` are horizontal components of velocity, `|V|` is the velocity magnitude, :math:`N_{ij}` is the number of turbines in cell :math:`(i,j)` (see Fig. :numref:`fig:WindFarm`), :math:`C_T` is the coefficient of thrust of the turbines, :math:`C_{TKE}` is the fraction of energy converted to turbulent kinetic energy -- both of which are functions of the velocity magnitude, and :math:`A_{ijk}` is the area intersected by the swept area of the turbine between levels :math:`z=z_k` and :math:`z= z_{k+1}`, and is explained in the next section.
 
-where `u` and `v` are horizontal components of velocity, `|V|` is the velocity magnitude, :math:`N^{ij}` is the number of turbines in cell :math:`(i,j)`, :math:`C_T` is the coefficient of thrust of the turbines, :math:`C_{TKE}` is the fraction of energy converted to turbulent kinetic energy -- both of which are functions of the velocity magnitude, and :math:`A_{ijk}` is the area intersected by the swept area of the turbine between levels :math:`z=z_k` and :math:`z= z_{k+1}`, and is explained in the next section.
 
 Intersected area :math:`A_{ijk}`
-_________________________________
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-Consider :math:`A_k^{k+1}` -- the area intersected by the swept area of the wind turbine between :math:`z=z_k` and :math:`z = z_{k+1}`. We have (see Figs. `2` and `3` below)
+Consider :math:`A_k^{k+1}` -- the area intersected by the swept area of the wind turbine between :math:`z=z_k` and :math:`z = z_{k+1}`. We have (see Figs. :numref:`fig:WindTurbine_Fitch` and :numref:`fig:Fitch_Aijk` below)
 
 .. math::
 
     A_k = \frac{\pi R^2}{2} - A_{ks}
 
-where :math:`A_{ks}` is the area of the segment of the circle as shown in Fig. `3`. We have from geometry, :math:`d_k = \min(|z_k - z_c|,R)` is the perpendicular distance of the center of the turbine to :math:`z = z_k`, where :math:`z_c` is the height of the center of the turbine from the ground. The area of the segment is
+where :math:`A_{ks}` is the area of the segment of the circle as shown in Fig. :numref:`fig:Fitch_Aijk`. We have from geometry, :math:`d_k = \min(|z_k - h|,R)` is the perpendicular distance of the center of the turbine to :math:`z = z_k`, where :math:`h` is the height of the center of the turbine from the ground. The area of the segment is
 
 .. math::
 
@@ -50,32 +53,164 @@ Hence, we have the intersected area :math:`A_{ijk}\equiv A_k^{k+1}` as
 
     A_k^{k+1} =
     \begin{cases}
-        |A_k - A_{k+1}| & \text{if } (z_k - z_c)(z_{k+1}-z_c) > 0 \\
-        |A_k + A_{k+1}| & \text{if } (z_k - z_c)(z_{k+1}-z_c) \le 0 \\
+        |A_k - A_{k+1}| & \text{if } (z_k - h)(z_{k+1}-h) > 0 \\
+        |A_k + A_{k+1}| & \text{if } (z_k - h)(z_{k+1}-h) \le 0 \\
     \end{cases}
 
 An example of the Fitch model is in ``Exec/Fitch``
 
-.. 1:
+.. _fig:WindFarm:
 
 .. figure:: ../figures/WindFarm_Fitch.png
    :width: 300
    :align: center
 
    Horizontal view of the wind farm in the Fitch model -- shows a wind farm in cell `(i,j)` with 5 wind turbines. The turbines are discretized only in the vertical direction.
 
-.. 2:
+.. _fig:WindTurbine_Fitch:
 
 .. figure:: ../figures/WindTurbine_Fitch.png
-   :width: 300
+   :width: 400
    :align: center
 
    The vertical discretization of the wind turbine in the Fitch model.
 
-.. 3:
+.. _fig:Fitch_Aijk:
 
 .. figure:: ../figures/FitchModel_A_ijk.png
    :width: 400
    :align: center
 
    The area terminology in the Fitch model. The circle represents the area swept by the turbine blades.
+
+
+.. _explicit-wake-parametrization-ewp-model:
+
+Explicit Wake Parametrization (EWP) model
+-----------------------------------------
+
+The Explicit Wake Parametrization (EWP) model [`Volker et al. 2015`_] is very similar to the Fitch model, and models the effect of wind farms as source terms in the governing equations for the horizontal components of momentum (i.e., :math:`x` and :math:`y` momentum) and the turbulent kinetic energy (TKE). At a given cell :math:`(i,j,k)`, the source terms in the governing equations are:
+
+.. math::
+    \frac{\partial u_{ijk}}{\partial t} = -\sqrt{\frac{\pi}{8}}\frac{N_{ij}c_tr_0^2\overline{u}_0^2}{\Delta x \Delta y \sigma_e}
+    \exp\left\{-\frac{1}{2}\left(\frac{z-h}{\sigma_e}\right)^2\right\}\cos(\phi(k))
+
+.. math::
+    \frac{\partial v_{ijk}}{\partial t} = -\sqrt{\frac{\pi}{8}}\frac{N_{ij}c_tr_0^2\overline{u}_0^2}{\Delta x \Delta y \sigma_e}
+    \exp\left\{-\frac{1}{2}\left(\frac{z-h}{\sigma_e}\right)^2\right\}\sin(\phi(k))
+
+.. math::
+    \frac{\partial \text{TKE}_{ijk}}{\partial t} = -N_{ij}\rho A_rc_t\langle \overline{u}_{i,h}\overline{u'^2}_{i,h}\rangle
+
+with
+
+.. math::
+    \sigma_e = \frac{\overline{u}_0}{3KL}\left[\left(\frac{2KL}{\overline{u}_0} + \sigma_0^2\right)^{\frac{3}{2}} - \sigma_0^3\right]
+
+where :math:`N_{ij}` is the number of turbines in cell :math:`(i,j)`, :math:`c_t` is the thrust coefficient, :math:`r_0` is the rotor radius, :math:`\overline{u}_0` is the mean advection velocity at hub height, :math:`h` is the hub height, :math:`\sigma_0 \approx 1.7 r_0` [`Volker et al. 2017`_] is a length scale that accounts for near-wake expansion, :math:`L` is the downstream distance that the wake travels within the cell approximated as a fraction of the cell size, :math:`K` is the turbulence eddy diffusivity (:math:`m^2/s`), :math:`\Delta x` and :math:`\Delta y` are the mesh sizes in the horizontal directions, and :math:`\phi(k)` is the wind direction with respect to the x-axis. :math:`\overline{u}_{i,h}` and :math:`\overline{u'}_{i,h}` are the mean and the fluctuating values of the velocity components (subscript :math:`i` is the direction index) at the hub height :math:`h`, :math:`A_r = \pi r^2` is the swept area of the rotor and :math:`\rho` is the density of air.
+
+The EWP model does not have a concept of intersected area by the turbine rotor like the Fitch model (see :ref:`Fitch model`). The exponential factor in the source terms for the velocities models the effect of the rotor for the momentum sinks (see Fig. :numref:`fig:WindTurbine_EWP`), and the turbulent kinetic energy source term only depends on the density, hub-height mean velocities and fluctuations, and the total swept area of the rotor :math:`A_r`.
+
+.. _fig:WindTurbine_EWP:
+
+.. figure:: ../figures/WindTurbine_EWP.png
+   :width: 400
+   :align: center
+
+   In the EWP model, the exponential factor in the source terms for the velocities models the effect of the rotor for the momentum sinks unlike the Fitch model which computes the
+   intersected area (see Fig. :numref:`fig:WindTurbine_Fitch`).
+
+.. _`Volker et al. 2015`: https://gmd.copernicus.org/articles/8/3715/2015/
+.. _`Volker et al. 2017`: https://iopscience.iop.org/article/10.1088/1748-9326/aa5d86
+
+
+.. _Inputs:
+
+Inputs for wind farm parametrization models
+------------------------------------------------------------
+
+Fitch, EWP
+~~~~~~~~~~~
+
+The following are the inputs required for simulations with Fitch, EWP models.
+
+.. code-block:: cpp
+
+    // The parametrization model - Fitch, EWP
+    erf.windfarm_type = "Fitch"
+
+    // How are the turbine locations specified? - using latitude-longitude
+    // format or x-y format? lat_lon or x_y
+    erf.windfarm_loc_type = "lat_lon"
+
+    // If using lat_lon, then the latitude and longitude of
+    // the lower bottom corner of the domain has to be specified
+    // The following means 35 deg N, 100 deg W (note the negative sign)
+    erf.latitude_lo      =   35.0
+    erf.longitude_lo     = -100.0
+
+    // Table having the wind turbine locations
+    erf.windfarm_loc_table = "windturbines_1WT.txt"
+
+    // Table having the specifications of the wind turbines. All turbines are assumed to
+    // have the same specifications
+    erf.windfarm_spec_table = "wind-turbine_1WT.tbl"
+
+1. ``erf.windfarm_type`` has to be one of the supported models - ``Fitch``, ``EWP``.
+2. ``erf.windfarm_loc_type`` is a variable to specify how the wind turbine locations in the wind farm is specified. If using the latitude and longitude of the turbine location, this has to be ``lat_lon`` or if using x and y co-ordinates to specify the turbine locations, this input is ``x_y``.
+
+   - If using ``lat_lon`` format, ``erf.latitude_lo`` and ``erf.longitude_lo`` specifies the latitude and longitude of the lower bottom corner of the domain box.  ie. if the domain box is specified as
+
+     .. code-block:: cpp
+
+        geometry.prob_lo     = -25000.   0.  -10000
+        geometry.prob_hi     =  25000. 10000. 10000.
+
+     then ``(erf.latitude_lo, erf.longitude_lo)`` corresponds to ``(-25000, 0, -10000)``.
+
+   - If using ``x_y`` format, there is no need to specify the ``erf.latitude_lo`` and ``erf.longitude_lo``.
+
+3. The ``erf.windfarm_loc_table`` specifies the locations of the wind turbines in the wind farm.
+
+   - For the latitude-longitude format, an example is as below. Each line specifies the latitude and longitude of the wind turbine location. The third entry simply has to be always 1 (WRF requires the third entry to be always 1, so maintaining same format here). The first entry means that the turbine is located at ``35.7828 deg N, 99.0168 deg W`` (note the negative sign in the entry corresponding to West).
+
+     .. code-block:: cpp
+
+        35.7828828829 -99.0168 1
+        35.8219219219 -99.0168 1
+        35.860960961 -99.0168 1
+        35.9 -99.0168 1
+        35.939039039 -99.0168 1
+        35.9780780781 -99.0168 1
+        36.0171171171 -99.0168 1
+        35.7828828829 -98.9705333333 1
+
+   - For the x-y format, an example is as below. Each line specifies the x and y co-ordinates of the wind turbine location in metres
+
+     .. code-block:: cpp
+
+        89264.99080053 91233.3333309002
+        89259.1966417755 95566.6666710007
+        89254.1277665419 99900.0000000001
+        89249.7842982733 104233.333329
+        89246.1663427532 108566.6666691
+        89243.2739881117 112899.9999981
+        93458.6633652711 86900.0000019001
+        93450.4377452458 91233.3333309002
+        93442.9032518779 95566.6666710007
+
+4. The ``erf.windfarm_spec_table`` gives the specifications of the wind turbines in the wind farm. All wind turbines are assumed to have the same specifications. Here is a sample specifications table.
+
+.. code-block:: cpp
+
+    4
+    119.0 178.0 0.130 2.0
+    9   0.805    50.0
+    10   0.805    50.0
+    11   0.805    50.0
+    12   0.805    50.0
+
+The first line is the number of pairs of entries for the power curve and thrust coefficient (there are 4 entries in this table which are in the last four lines of the table).
+The second line gives the height in meters of the turbine hub, the diameter in
+meters of the rotor, the standing thrust coefficient, and the nominal power of the turbine in MW.
+The remaining lines (four in this case) contain the three values of: wind speed (m/s), thrust coefficient, and power production in kW.