Nacelle model (#3)

* Added body force term and changed filters from \tilde to \widetilde

* Modified momentum forcing source term to include density averaging.

* Added documentation for ABL source term

* Removed comment for debugging

* Updated ablForcingEdge norm gold standard for latest ABL forcing subroutine

* Added Shreyas' comments on the pull request.

* Removed commented code for density

* Added braces to fix scope of if-statements

* Fixed infenitesimally thins slice sentece in documentation

* Changed order of variable declarations to avoid compiler warning

* Setting up for a different spreading function for the nacelle point for actuator line simulations of wind turbines

* Added Nacelle Model

* Adding input of air density to calculate force on Nacelle

* Nacelle Model incorporated.

* Nacelle model is only active if cd is grater than zero

* Assign standard epsilon value to hub when nacelle not active

* Updated reg test for nrel5MW by adding the nacelle and gold standard file

* Added documentation file for the turbine models

* Changed the gold standard to work with latest version of OpenFAST, added turbine documentation file, and fixed minor formatting in actuator line file.

* Resolving existing conflicts in windEnergy.rst.

* Removing nrel5MWactuatorLine.norm file.

docs/references/references.bib

@@ -305,3 +305,13 @@ @techreport{FuegoTheoryManual:2016
 	Number = {SAND2016-10163},
 	Title = {SIERRA Low Mach Module: Fuego Theory Manual - Version 4.42},
 	Year = {2016}}
+
+@phdthesis{Martinez-Tossas2017,
+  author       = {Luis A Mart{\'i}nez-Tossas}, 
+  title        = {Large Eddy Simulations and Theoretical Analysis of Wind
+                  Turbine Aerodynamics Using an Actuator Line Model},
+  school       = {Johns Hopkins University},
+  year         = 2017,
+  address      = {Baltimore, MD USA},
+  month        = {July}
+}

docs/source/theory/windEnergy.rst

@@ -179,18 +179,18 @@ described here. The Navier-Stokes equations are written as
 
 .. math::
    :label: fav-mom-nalu
-           
+
    {\bf F}_i (\rho^{n+1}, u_i^{n+1}, P^{n+1}) - \int \left . \frac{\partial \rho u_i}{\partial t} \right |^{n+1} {\rm d}V &= 0, \\
    {\bf F}_i (\rho^{n+1}, u_i^{n+1}, P^{n+1}) - \frac{ (\gamma_1 \rho^{n+1} {u_i}^{n+1} + \gamma_2 \rho^n {u_i}^{n} + \gamma_3 \rho^n {u_i}^{n-1})}{\Delta t} \Delta V &=0,
 
 where
 
 .. math::
-           
+
    {\bf F}_i (\rho^{n+1} u_i^{n+1}) &= - \int \rho^{n+1} u_i^{n+1} u_j^{n+1} n_j {\rm d}S  + \int \tau_{ij}^{n+1} n_j {\rm d}S - \int P^{n+1} n_i {\rm d}S - \int \left(\rho^{n+1} - \rho_{\circ} \right) g_i {\rm d}V, \\
    &= - \int u_i^{n+1} \dot{m}^{n+1}  + \int \tau_{ij}^{n+1} n_j {\rm d}S  - \int P^{n+1} n_i {\rm d}S - \int \left(\rho^{n+1} - \rho_{\circ} \right) g_i {\rm d}V. \\
-   
-   
+
+
 and :math:`\gamma_i` are factors for BDF2 time discretization scheme (see
 :numref:`theory_time_discretization`). As is typical of incompressible flow
 solvers, the mass flow rate through the sub-control surfaces is tracked
@@ -208,7 +208,7 @@ solution to the Navier-Stokes equations at time step :math:`n+1` as
 
 .. math::
    :label: fav-mom-nalu-newton
-           
+
    \mathbf{A}_{ij} \; \delta u_{j} &= {\bf F}_i^{*} - \frac{ (\gamma_1 \rho^{*} {u_i}^{*} + \gamma_2 \rho^n {u_i}^{n} + \gamma_3 \rho^n {u_i}^{n-1})}{\Delta t} \Delta V, \\
    \textrm{where } \delta u_{j} &= u_i^{**} - u_i^*, \\
    \mathbf{A}_{ij} &= \left ( \frac{ \gamma_1 \rho^{*}}{\Delta t} \Delta V \delta_{ij} - \left . \frac{\partial F_i}{\partial u_j} \right |^{*} \right ), \\
@@ -219,15 +219,15 @@ After each Newton or *outer* iteration, :math:`\phi^{**}` is a better approximat
 
 .. math::
    :label: linearize-f-phi-star
-           
-   {\bf F}_i^* = \left . \frac{\partial F_i}{\partial u_j} \right |^{*} u_j^{*} - \int P^{*} n_i {\rm d}S - \int \left(\rho^{*} - \rho_{\circ} \right) g_i {\rm d}V 
-   
+
+   {\bf F}_i^* = \left . \frac{\partial F_i}{\partial u_j} \right |^{*} u_j^{*} - \int P^{*} n_i {\rm d}S - \int \left(\rho^{*} - \rho_{\circ} \right) g_i {\rm d}V
+
 Applying Eq. :eq:`linearize-f-phi-star` to Eq. :eq:`fav-mom-nalu-newton`, we get the
 linearized momentum predictor equation solved in Nalu-Wind.
 
-.. math::   
+.. math::
    :label: fav-mom-nalu-linearize-f
-      
+
    {\bf A}_{ij} \; \delta u_j &= \left . \frac{\partial F_i}{\partial u_j} \right |^{*} u_j^{*} - \int P^{*} n_i {\rm d}S - \int \left(\rho^{*} - \rho_{\circ} \right) g_i {\rm d}V  \\
    & \quad \quad  - \frac{ (\gamma_1 \rho^{*} {u_i}^{*} + \gamma_2 \rho^{n} {u_i}^{n} + \gamma_3 \rho^{n-1} {u_i}^{n-1})}{\Delta t} \Delta V \\
    {\bf A}_{ij} \; \delta u_j &= \left (\frac{ \gamma_1 \rho^{*}}{\Delta t} \Delta V \delta_{ij} - \left . \frac{\partial F_i}{\partial u_j} \right |^{*} \right ) {u_j}^{*} - \int P^{*} n_i {\rm d}S - \int \left(\rho^{*} - \rho_{\circ} \right) g_i {\rm d}V \\
@@ -237,7 +237,7 @@ linearized momentum predictor equation solved in Nalu-Wind.
 
 :math:`u_i^{**}` will not satisfy the continuity equation. A correction step is
 performed at the end of each outer iteration to make :math:`u_i^{**}`
-satisfy the continuity equation as 
+satisfy the continuity equation as
 
 .. math::
 
@@ -250,14 +250,14 @@ equation to be satisfied along with the splitting and stabilization errors is
 
 .. math::
    :label: eq-continuity
-           
+
    {\bf D } \rho u^{**} = b + \left ({\bf L_1} - {\bf D} \tau_3 {\bf G} \right ) \Delta P^{**} + \left ({\bf L_2} - {\bf D} \tau_2 {\bf G} \right ) P^{*}
 
 where :math:`b` contains any source terms when the velocity field is not
 divergence free and the other terms are the errors due to pressure stabilization
 as shown by Domino :cite:`Domino:2006`. The final pressure Poisson equation
 solved to enforce continuity at each outer iteration is
-   
+
 .. math::
    :label: eq-pressure
 
@@ -267,7 +267,7 @@ solved to enforce continuity at each outer iteration is
    b + {\bf L_1} \Delta P^{**} &= {\bf D} (\rho \widehat{u}) - \left ({\bf L_2} - {\bf D} \tau_2 {\bf G} \right ) P^{*} \\
    -{\bf L_1} \Delta P^{**} &= {\bf D} \rho \widehat{u} + {\bf D} \tau_2 {\bf G} P^{*} - {\bf L_2} P^{*} \\
    -{\bf L_1} \Delta P^{**} &= - {\bf D} \rho \widehat{u} - {\bf D} \tau_2 {\bf G} P^{*} + {\bf L_2} P^{*} + b
-  
+
 Thus, the final set of equations solved at each outer iteration is
 
 .. math::
@@ -371,9 +371,9 @@ given by
 
 .. math::
    :label: force-integral
-           
+
    f_i(x,y,z) = \int_0^L g\left(\vec{r}\left(l\right)\right) F'_i\left(l\right) \: \textrm{d} l.
-   
+
 
 Here, the projection function’s position vector is a function of
 position on the actuator line. The part of the line nearest to the point in
@@ -392,9 +392,9 @@ This is convenient because the integral given in Equation
 
 .. math::
    :label: force-summation
-           
+
    f_i(x,y,z) = \sum_{k=0}^N g(\vec{r}^k) F_i^k.
-   
+
 
 This summation well approximates the integral given in Equation
 :eq:`force-integral` so long as the ratio of actuator element size to
@@ -414,6 +414,50 @@ module to OpenFAST by supplying the velocity field information.
 Nalu-Wind -- OpenFAST Coupling Algorithm
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
+A nacelle model is implemented using a Gaussian drag body force. The model
+implements a drag force in a direction opposite to velocity field at the center
+of the Gaussian.
+The width of the Gaussian kernel is determined using the reference
+area and drag coefficient of the nacelle as shown by Martinez-Tossas.
+:cite:`Martinez-Tossas2017`
+
+.. math::
+
+   \epsilon_d = \sqrt{2 \, c_d \, A/ \pi}
+
+
+where :math:`c_d` is the drag coefficient,
+and
+:math:`A` is the reference area.
+This value of :math:`\epsilon_d` ensures that the momentum thickness of the
+generated wake is of the right size.
+The velocity sampled at the center of the Gaussian is corrected
+to obtain the upstream velocity.
+
+.. math::
+
+   {\widetilde{u}_i}_\infty = \frac{1}{1 -
+                               c_d \, A / (4 \, \pi \, \epsilon_d^2)}
+                               {\widetilde{u}_i}_c
+
+
+where :math:`\widetilde{u_i}_c` is the velocity at the center of the Gaussian
+and :math:`{\widetilde{u}_i}_\infty` is the free-stream velocity used to compute
+the drag force.
+The drag body force is then
+
+.. math::
+
+   f_d(x,y,z) = \frac{1}{2} \, c_d \, A\, {\widetilde{u}_i}_\infty^2 \,
+                \frac{1}{\pi^{3/2} \epsilon_d^3} e^{-|\vec{r}|^2/\epsilon_d^2}
+
+
+where :math:`\vec{r}` is the position vector between the fluid point of
+interest and the center of the Gaussian force.
+
+Nalu -- OpenFAST Coupling Algorithm
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
 The actuator line implementation allows for flexible blades that are not
 necessarily straight (prebend and sweep). The current implementation requires a
 fixed time step when coupled to OpenFAST, but allows the time step in Nalu-Wind to be

docs/source/user/nalu_run/nalu_inp.rst

@@ -634,7 +634,7 @@ specifies the component of the gravity vector, defined in
 :inpfile:`solution_options.gravity`, that should be used in the
 definition of the Monin-Obukhov length scale calculation.  The
 entry :inpfile:`reference_temperature` is the reference temperature
-used in calculation of the Monin-Obukhov length scale. 
+used in calculation of the Monin-Obukhov length scale.
 
 When there is mesh motion involved the wall boundary must specify a user
 function to determine relative velocity at the surface.
@@ -1134,6 +1134,8 @@ Actuator
 .. inpfile:: actuator.turb_id
 
    A unique turbine id for each turbine
+
+.. include:: ./turbine_modeling.rst
 
 
 Turbulence averaging
@@ -1200,7 +1202,7 @@ Turbulence averaging
      "Moving window" average where the window size is set to to the time
      filter width. The contribution of any quantity before the moving window
      towards the average value reduces exponentially with every time step.
-   
+
 .. inpfile:: turbulence_averaging.time_filter_interval
 
    Number indicating the time filter size over which to calculate the
@@ -1239,12 +1241,12 @@ Turbulence averaging
 
 .. inpfile:: turbulence_averaging.specifications.compute_resolved_stress
 
-   A boolean flag indicating whether the average resolved stress is 
+   A boolean flag indicating whether the average resolved stress is
    computed as :math:`< \bar\rho \widetilde{u_i} \widetilde{u_j} >`.
    The default value is ``no``. When this option is turned on, the Favre
    average of the resolved velocity, :math:`< \bar{\rho} \widetilde{u_j} >`, is
    computed as well.
-   
+
 .. inpfile:: turbulence_averaging.specifications.compute_temperature_resolved_flux
 
    A boolean flag indicating whether the average resolved temperature flux is
@@ -1261,16 +1263,16 @@ Turbulence averaging
    by the turbulence model is used. The sub-filter scale kinetic energy is used
    to determine the isotropic component of the sub-filter stress. As described
    in the section :ref:`supp_eqn_set_mom_cons`, the Yoshizawa model is used to
-   compute the sub-filter kinetic energy when it is not transported. 
-   
+   compute the sub-filter kinetic energy when it is not transported.
+
 .. inpfile:: turbulence_averaging.specifications.compute_temperature_sfs_flux
 
    A boolean flag indicating whether the average sub-filter scale flux of
    temperature is computed. The default value is ``no``. The sub-filter scale
    stress model is assumed to be of an eddy diffusivity type and the turbulent
    diffusivity computed by the turbulence model is used along with a constant
    turbulent Prandtl number obtained from the Realm.
-   
+
 .. inpfile:: turbulence_averaging.specifications.compute_favre_stress
 
    A boolean flag indicating whether the Favre stress is computed. The

docs/source/user/nalu_run/turbine_modeling.rst

@@ -0,0 +1,156 @@
+Actuator Turbine Model
+``````````````````````
+
+.. inpfile:: actuator
+
+   ``actuator`` subsection defines the inputs for actuator line simulations. A
+   sample section is shown below for running actuator line simulations
+   coupled to OpenFAST with two turbines.
+
+.. code-block:: yaml
+
+     actuator:
+       type: ActLineFAST
+       search_method: boost_rtree
+       search_target_part: Unspecified-2-HEX
+
+       n_turbines_glob: 2
+       dry_run:  False
+       debug:    False
+       t_start: 0.0
+       simStart: init # init/trueRestart/restartDriverInitFAST
+       t_max:    5.0
+       n_every_checkpoint: 100
+
+       Turbine0:
+         procNo: 0
+         num_force_pts_blade: 50
+         num_force_pts_tower: 20
+         nacelle_cd: 1.0
+         nacelle_area: 1.0
+         air_density: 1.225
+         epsilon: [ 5.0, 5.0, 5.0 ]
+         turbine_base_pos: [ 0.0, 0.0, -90.0 ]
+         turbine_hub_pos: [ 0.0, 0.0, 0.0 ]
+         restart_filename: ""
+         FAST_input_filename: "Test01.fst"
+         turb_id:  1
+         turbine_name: machine_zero
+
+       Turbine1:
+         procNo: 0
+         num_force_pts_blade: 50
+         num_force_pts_tower: 20
+         nacelle_cd: 1.0
+         nacelle_area: 1.0
+         air_density: 1.225
+         epsilon: [ 5.0, 5.0, 5.0 ]
+         turbine_base_pos: [ 250.0, 0.0, -90.0 ]
+         turbine_hub_pos: [ 250.0, 0.0, 0.0 ]
+         restart_filename: ""
+         FAST_input_filename: "Test02.fst"
+         turb_id:  2
+         turbine_name: machine_one
+
+
+.. inpfile:: actuator.type
+
+   Type of actuator source. Options are ``ActLineFAST`` and ``ActLinePointDrag``. Only ``ActLineFAST`` is documented here.
+
+.. inpfile:: actuator.search_method
+
+   String specifying the type of search method used to identify the nodes within the search radius of the actuator points. Options are ``boost_rtree`` and ``stk_kdtree``. The default is ``stk_kdtree`` when the ``search_type`` is not specified.
+
+.. inpfile:: search_target_part
+
+   String or an array of strings specifying the parts of the mesh to be searched to identify the nodes near the actuator points.
+
+.. inpfile:: actuator.n_turbines_glob
+
+   Total number of turbines in the simulation. The input file must contain a number of turbine specific sections (`Turbine0`, `Turbine1`, ..., `Turbine(n-1)`) that is consistent with `nTurbinesGlob`.
+
+.. inpfile:: actuator.debug
+
+   Enable debug outputs if set to true
+
+.. inpfile:: actuator.dry_run
+
+   The simulation will not run if dryRun is set to true. However, the simulation will read the input files, allocate turbines to processors and prepare to run the individual turbine instances. This flag is useful to test the setup of the simulation before running it.
+
+.. inpfile:: actuator.simStart
+
+   Flag indicating whether the simulation starts from scratch or restart. ``simStart`` takes on one of three values:
+
+   * ``init`` - Use this option when starting a simulation from `t=0s`.
+   * ``trueRestart`` - While OpenFAST allows for restart of a turbine simulation, external components like the Bladed style controller may not. Use this option when all components of the simulation are known to restart.
+   * ``restartDriverInitFAST`` - When the ``restartDriverInitFAST`` option is selected, the individual turbine models start from `t=0s` and run up to the specified restart time using the inflow data stored at the actuator nodes from a hdf5 file. The C++ API stores the inflow data at the actuator nodes in a hdf5 file at every OpenFAST time step and then reads it back when using this restart option. This restart option is especially useful when the glue code is a CFD solver.
+
+.. inpfile:: actuator.t_start
+
+   Start time of the simulation
+
+.. inpfile:: actuator.t_end
+
+   End time of the simulation. ``t_end`` <= ``t_max``
+
+.. inpfile:: actuator.t_max
+
+   Max time of the simulation
+
+
+.. note::
+
+   ``t_max`` can only be set when OpenFAST is running from `t=0s` and ``simStart`` is ``init``. ``t_max`` can not be changed on a restart. OpenFAST will not be able to run beyond ``t_max``. Choose ``t_max`` to be large enough to accomodate any possible future extensions of runs. One can change ``t_start`` and ``t_end`` to start and stop the simulation any number of times as long as ``t_end`` <= ``t_max``.
+
+.. inpfile:: actuator.dt_fast
+
+   Time step for OpenFAST. All turbines should have the same time step.
+
+.. inpfile:: actuator.n_every_checkpoint
+
+   Restart files will be written every so many time steps
+
+**Turbine specific input options**
+
+.. inpfile:: actuator.turbine_base_pos
+
+   The position of the turbine base for actuator-line simulations
+
+.. inpfile:: actuator.num_force_pts_blade
+
+   The number of actuator points along each blade for actuator-line simulations
+
+.. inpfile:: actuator.num_force_pts_tower
+
+   The number of actuator points along the tower for actuator-line simulations.
+
+.. inpfile:: actuator.nacelle_cd
+
+   The drag coefficient for the nacelle. If this is set to zero, or not
+   defined, the code will not implement the nacelle model.
+
+.. inpfile:: actuator.nacelle_area
+
+   The reference area for the nacelle. This is only used if the nacelle
+   model is used.
+
+.. inpfile:: actuator.air_density
+
+   The air density. This is only used to compute the nacelle force. It should
+   match the density being used in both Nalu and OpenFAST.
+
+.. inpfile:: actuator.epsilon
+
+   The spreading width :math:`\epsilon` in the Gaussian spreading function in the `[chordwise, spanwise, chord normal]` coordinate system to spread the forces from the actuator point to the nodes. Nalu currently only supports an isotropic Gaussian spreading function and uses only the value in the first component along the `chordwise` direction.
+
+.. inpfile:: actuator.restart_filename
+
+   The checkpoint file for this turbine when restarting a simulation
+
+.. inpfile:: actuator.FAST_input_filename
+
+   The FAST input file for this turbine
+
+.. inpfile:: actuator.turb_id
+
+   A unique turbine id for each turbine