Improve various docstrings (#1171)

* Update references to floris.simulation and add docstring to get_defaults

* Improve docstrings for user-level functionality

* Update args docstring

* Formatting

* Update docstring

* Fix WETO software hyperlink

* Fix documentation for cut planes and add documentation for visualizations

* Add Katic reference

* Switch from r to r_squared to be explicit and avoid confusion

* Another reference to Katic

* Add notes regarding CH initial exponent

* Spelling fix

* Square image, update email addresses

* Update todo note

Author: misi9170

docs/dev_guide.md

@@ -174,7 +174,7 @@ git commit -m "Update so and so"
 
 In order to maintain a level of confidence in the software, FLORIS is expected
 to maintain a reasonable level of test coverage. To that end, unit
-tests for a the low-level code in the `floris.simulation` package are included.
+tests for a the low-level code in the `floris.core` package are included.
 
 The full testing suite can by executed by running the command ``pytest`` from
 the highest directory in the repository. A testing-only class is included
@@ -404,7 +404,7 @@ def function(
 ```
 
 Some models require a special grid and/or solver, and that mapping happens in
-[floris.simulation.Floris](https://github.com/NatLabRockies/floris/blob/main/floris/simulation/floris.py#L145).
+[floris.core.core.Core](https://github.com/NatLabRockies/floris/blob/main/floris/core/core.py).
 Generally, a specific kind of solver requires one or a number of specific grid-types.
 For example, `full_flow_sequential_solver` requires either `FlowFieldGrid` or
 `FlowFieldPlanarGrid`.

docs/docs_image.png

@@ -146,7 +146,7 @@ @Article{annoni2018analysis
 
 @article{crespo1996turbulence,
   title={Turbulence characteristics in wind-turbine wakes},
-  author={Crespo, A and Hernandez, J},
+  author={Crespo, A. and Hernández, J.},
   journal={Journal of wind engineering and industrial aerodynamics},
   volume={61},
   number={1},
@@ -369,3 +369,36 @@ @techreport{zahle_IEA22MW_2024
 	author = {Zahle, Frederik and Barlas, Athanasios and Loenbaek, Kenneth and Bortolotti, Pietro and Zalkind, Daniel and Wang, Lu and Labuschagne, Casper and Sethuraman, Latha and Barter, Garrett},
 	year = {2024},
 }
+
+@article{fleming_sr_2022,
+	title = {Serial-Refine Method for Fast Wake-Steering Yaw Optimization},
+	volume = {2265},
+	issn = {1742-6588, 1742-6596},
+	url = {https://iopscience.iop.org/article/10.1088/1742-6596/2265/3/032109},
+	doi = {10.1088/1742-6596/2265/3/032109},
+	number = {3},
+	journal = {Journal of Physics: Conference Series (TORQUE)},
+	author = {Fleming, Paul A. and Stanley, Andrew P. J. and Bay, Christopher J. and King, Jennifer and Simley, Eric and Doekemeijer, Bart M. and Mudafort, Rafael},
+	year = {2022},
+	pages = {032109},
+}
+
+@inproceedings{katic_sos_1986,
+	address = {Rome, Italy},
+	title = {A simple model for cluster efficiency},
+	volume = {1},
+	author = {Katic, I and Højstrup, J and Jensen, N O},
+	year = {1986},
+	pages = {407--410},
+}
+
+@article{zehtabiyan_rezaie_CH_2023,
+	title = {A short note on turbulence characteristics in wind-turbine wakes},
+	volume = {240},
+	issn = {0167-6105},
+	doi = {10.1016/j.jweia.2023.105504},
+	journal = {Journal of Wind Engineering and Industrial Aerodynamics},
+	author = {Zehtabiyan-Rezaie, Navid and Abkar, Mahdi},
+	year = {2023},
+	pages = {105504},
+}

docs/references.bib

@@ -213,7 +213,9 @@
     "\n",
     "CrespoHernandez is a wake-turbulence model that is used to compute additional variability introduced\n",
     "to the flow field by operation of a wind turbine. Implementation of the model follows the original\n",
-    "formulation and limitations outlined in {cite:t}`crespo1996turbulence`."
+    "formulation and limitations outlined in {cite:t}`crespo1996turbulence`.\n",
+    "\n",
+    "The default parameter values used in FLORIS for CrespoHernandez differ from those reported in {cite:t}`crespo1996turbulence` following subsequent calibration. However, {cite:t}`zehtabiyan_rezaie_CH_2023` argue that the sign of certain parameters are not physically consistent (and also misreported in subsequent literature). See the `CrespoHernandez` class docstring for more details."
    ]
   },
   {
@@ -347,8 +349,7 @@
    "source": [
     "### Sum of Squares Freestream Superposition (SOSFS)\n",
     "\n",
-    "This model combines the wakes via a sum of squares of the new wake to add and the existing\n",
-    "flow field."
+    "This model combines the wakes via a sum of squares of the new wake to add and the existing flow field. For more information, refer to :cite:`katic_sos_1986`."
    ]
   },
   {

docs/wake_models.ipynb

@@ -1,7 +1,8 @@
-"""Example: Check turbine power curves
+"""Example: Reference turbines
 
-For each turbine in the turbine library, make a small figure showing that its power
-curve and power loss to yaw are reasonable and reasonably smooth
+For each reference wind turbine in the turbine library, make a small figure
+showing its power and thrust coefficient curves and demonstrate its power loss
+to yaw.
 """
 
 

examples/examples_turbine/001_reference_turbines.py

@@ -10,6 +10,8 @@ class SOSFS(BaseModel):
     """
     SOSFS uses sum of squares freestream superposition to combine the
     wake velocity deficits to the base flow field.
+
+    For more information, refer to :cite:`katic_sos_1986`.
     """
 
     def prepare_function(self) -> dict:

floris/core/wake_combination/sosfs.py

@@ -297,29 +297,31 @@ def wake_added_yaw(
     # top vortex
     # NOTE: this is the top of the grid, not the top of the rotor
     zT = z_i - (HH + D / 2) + NUM_EPS  # distance from the top of the grid
-    rT = ne.evaluate("yLocs ** 2 + zT ** 2")  # TODO: This is (-) in the paper
+    # NOTE: This is (-) in the paper, but (+) is consistent with the
+    # Martínez-Tossas et al. (2019) source.
+    rT_squared = ne.evaluate("yLocs ** 2 + zT ** 2")
     # This looks like spanwise decay;
     # it defines the vortex profile in the spanwise directions
-    core_shape = ne.evaluate("1 - exp(-rT / (eps ** 2))")
-    v_top = ne.evaluate("(Gamma_top * zT) / (2 * pi * rT) * core_shape")
+    core_shape = ne.evaluate("1 - exp(-rT_squared / (eps ** 2))")
+    v_top = ne.evaluate("(Gamma_top * zT) / (2 * pi * rT_squared) * core_shape")
     v_top = np.mean( v_top, axis=(2,3) )
     # w_top = (-1 * Gamma_top * yLocs) / (2 * pi * rT) * core_shape * decay
 
     # bottom vortex
     zB = z_i - (HH - D / 2) + NUM_EPS
-    rB = ne.evaluate("yLocs ** 2 + zB ** 2")
-    core_shape = ne.evaluate("1 - exp(-rB / (eps ** 2))")
-    v_bottom = ne.evaluate("(Gamma_bottom * zB) / (2 * pi * rB) * core_shape")
+    rB_squared = ne.evaluate("yLocs ** 2 + zB ** 2")
+    core_shape = ne.evaluate("1 - exp(-rB_squared / (eps ** 2))")
+    v_bottom = ne.evaluate("(Gamma_bottom * zB) / (2 * pi * rB_squared) * core_shape")
     v_bottom = np.mean( v_bottom, axis=(2,3) )
     # w_bottom = (-1 * Gamma_bottom * yLocs) / (2 * pi * rB) * core_shape * decay
 
     # wake rotation vortex
     zC = z_i - HH + NUM_EPS
-    rC = ne.evaluate("yLocs ** 2 + zC ** 2")
-    core_shape = ne.evaluate("1 - exp(-rC / (eps ** 2))")
-    v_core = ne.evaluate("(Gamma_wake_rotation * zC) / (2 * pi * rC) * core_shape")
+    rC_squared = ne.evaluate("yLocs ** 2 + zC ** 2")
+    core_shape = ne.evaluate("1 - exp(-rC_squared / (eps ** 2))")
+    v_core = ne.evaluate("(Gamma_wake_rotation * zC) / (2 * pi * rC_squared) * core_shape")
     v_core = np.mean( v_core, axis=(2,3) )
-    # w_core = (-1 * Gamma_wake_rotation * yLocs) / (2 * pi * rC) * core_shape * decay
+    # w_core = (-1 * Gamma_wake_rotation * yLocs) / (2 * pi * rC_squared) * core_shape * decay
 
     # Cap the effective yaw values between -45 and 45 degrees
     val = 2 * (avg_v - v_core) / (v_top + v_bottom)
@@ -398,51 +400,59 @@ def calculate_transverse_velocity(
 
     # top vortex
     zT = z - (HH + D / 2) + NUM_EPS
-    rT = ne.evaluate("yLocs ** 2 + zT ** 2")  # TODO: This is - in the paper
+    # NOTE: This is (-) in the paper, but (+) is consistent with the
+    # Martínez-Tossas et al. (2019) source.
+    rT_squared = ne.evaluate("yLocs ** 2 + zT ** 2")
     # This looks like spanwise decay;
     # it defines the vortex profile in the spanwise directions
-    core_shape = ne.evaluate("1 - exp(-rT / (eps ** 2))")
-    V1 = ne.evaluate("(Gamma_top * zT) / (2 * pi * rT) * core_shape * decay")
-    W1 = ne.evaluate("(-1 * Gamma_top * yLocs) / (2 * pi * rT) * core_shape * decay")
+    core_shape = ne.evaluate("1 - exp(-rT_squared / (eps ** 2))")
+    V1 = ne.evaluate("(Gamma_top * zT) / (2 * pi * rT_squared) * core_shape * decay")
+    W1 = ne.evaluate("(-1 * Gamma_top * yLocs) / (2 * pi * rT_squared) * core_shape * decay")
 
     # bottom vortex
     zB = z - (HH - D / 2) + NUM_EPS
-    rB = ne.evaluate("yLocs ** 2 + zB ** 2")
-    core_shape = ne.evaluate("1 - exp(-rB / (eps ** 2))")
-    V2 = ne.evaluate("(Gamma_bottom * zB) / (2 * pi * rB) * core_shape * decay")
-    W2 = ne.evaluate("(-1 * Gamma_bottom * yLocs) / (2 * pi * rB) * core_shape * decay")
+    rB_squared = ne.evaluate("yLocs ** 2 + zB ** 2")
+    core_shape = ne.evaluate("1 - exp(-rB_squared / (eps ** 2))")
+    V2 = ne.evaluate("(Gamma_bottom * zB) / (2 * pi * rB_squared) * core_shape * decay")
+    W2 = ne.evaluate("(-1 * Gamma_bottom * yLocs) / (2 * pi * rB_squared) * core_shape * decay")
 
     # wake rotation vortex
     zC = z - HH + NUM_EPS
-    rC = ne.evaluate("yLocs ** 2 + zC ** 2")
-    core_shape = ne.evaluate("1 - exp(-rC / (eps ** 2))")
-    V5 = ne.evaluate("(Gamma_wake_rotation * zC) / (2 * pi * rC) * core_shape * decay")
-    W5 = ne.evaluate("(-1 * Gamma_wake_rotation * yLocs) / (2 * pi * rC) * core_shape * decay")
+    rC_squared = ne.evaluate("yLocs ** 2 + zC ** 2")
+    core_shape = ne.evaluate("1 - exp(-rC_squared / (eps ** 2))")
+    V5 = ne.evaluate("(Gamma_wake_rotation * zC) / (2 * pi * rC_squared) * core_shape * decay")
+    W5 = ne.evaluate(
+        "(-1 * Gamma_wake_rotation * yLocs) / (2 * pi * rC_squared) * core_shape * decay"
+    )
 
     ### Boundary condition - ground mirror vortex
 
     # top vortex - ground
     zTb = z + (HH + D / 2) + NUM_EPS
-    rTb = ne.evaluate("yLocs ** 2 + zTb ** 2")
+    rTb_squared = ne.evaluate("yLocs ** 2 + zTb ** 2")
     # This looks like spanwise decay;
     # it defines the vortex profile in the spanwise directions
-    core_shape = ne.evaluate("1 - exp(-rTb / (eps ** 2))")
-    V3 = ne.evaluate("(-1 * Gamma_top * zTb) / (2 * pi * rTb) * core_shape * decay")
-    W3 = ne.evaluate("(Gamma_top * yLocs) / (2 * pi * rTb) * core_shape * decay")
+    core_shape = ne.evaluate("1 - exp(-rTb_squared / (eps ** 2))")
+    V3 = ne.evaluate("(-1 * Gamma_top * zTb) / (2 * pi * rTb_squared) * core_shape * decay")
+    W3 = ne.evaluate("(Gamma_top * yLocs) / (2 * pi * rTb_squared) * core_shape * decay")
 
     # bottom vortex - ground
     zBb = z + (HH - D / 2) + NUM_EPS
-    rBb = ne.evaluate("yLocs ** 2 + zBb ** 2")
-    core_shape = ne.evaluate("1 - exp(-rBb / (eps ** 2))")
-    V4 = ne.evaluate("(-1 * Gamma_bottom * zBb) / (2 * pi * rBb) * core_shape * decay")
-    W4 = ne.evaluate("(Gamma_bottom * yLocs) / (2 * pi * rBb) * core_shape * decay")
+    rBb_squared = ne.evaluate("yLocs ** 2 + zBb ** 2")
+    core_shape = ne.evaluate("1 - exp(-rBb_squared / (eps ** 2))")
+    V4 = ne.evaluate("(-1 * Gamma_bottom * zBb) / (2 * pi * rBb_squared) * core_shape * decay")
+    W4 = ne.evaluate("(Gamma_bottom * yLocs) / (2 * pi * rBb_squared) * core_shape * decay")
 
     # wake rotation vortex - ground effect
     zCb = z + HH + NUM_EPS
-    rCb = ne.evaluate("yLocs ** 2 + zCb ** 2")
-    core_shape = ne.evaluate("1 - exp(-rCb / (eps ** 2))")
-    V6 = ne.evaluate("(-1 * Gamma_wake_rotation * zCb) / (2 * pi * rCb) * core_shape * decay")
-    W6 = ne.evaluate("(Gamma_wake_rotation * yLocs) / (2 * pi * rCb) * core_shape * decay")
+    rCb_squared = ne.evaluate("yLocs ** 2 + zCb ** 2")
+    core_shape = ne.evaluate("1 - exp(-rCb_squared / (eps ** 2))")
+    V6 = ne.evaluate(
+        "(-1 * Gamma_wake_rotation * zCb) / (2 * pi * rCb_squared) * core_shape * decay"
+    )
+    W6 = ne.evaluate(
+        "(Gamma_wake_rotation * yLocs) / (2 * pi * rCb_squared) * core_shape * decay"
+    )
 
     # total spanwise velocity
     V = V1 + V2 + V3 + V4 + V5 + V6

floris/core/wake_deflection/gauss.py

@@ -23,13 +23,29 @@ class CrespoHernandez(BaseModel):
     turbine. Implementation of the model follows the original formulation and
     limitations outlined in :cite:`cht-crespo1996turbulence`.
 
+    Note: The values for default parameters provided here differ from those in
+    :cite:`cht-crespo1996turbulence. Following their recommendations, the
+    default parameters would instead be:
+        - initial: -0.0325*
+        - constant: 0.73
+        - ai: 0.8325
+        - downstream: -0.32
+    * The "initial" parameter is given as -0.0325 in :cite:`cht-crespo1996turbulence`,
+    but the negative exponent is not clear in the scans of the paper found on the internet,
+    and several subsequent paper cite the exponent as positive (0.0325). This discrepancy
+    is noted in :cite:`zehtabiyan_rezaie_CH_2023`. Moreover, :cite:`zehtabiyan_rezaie_CH_2023`
+    argues that positive values for this exponent are not representative of the physical
+    phenomena occurring. For more details, see https://github.com/NREL/floris/issues/773.
+    Nonetheless, the default value here is set to 0.1 for consistency with previous
+    FLORIS versions. The default value may be updated in a future release.
+
     Args:
         parameter_dictionary (dict): Model-specific parameters.
             Default values are used when a parameter is not included
             in `parameter_dictionary`. Possible key-value pairs include:
 
-            -   **initial** (*float*): The initial ambient turbulence
-                intensity, expressed as a decimal fraction.
+            -   **initial** (*float*): The exponent on the initial ambient
+                turbulence intensity.
             -   **constant** (*float*): The constant used to scale the
                 wake-added turbulence intensity.
             -   **ai** (*float*): The axial induction factor exponent used

floris/core/wake_turbulence/crespo_hernandez.py

@@ -132,7 +132,7 @@ def function(
             sigma_z *= (x >= xR)
             sigma_z += np.ones_like(sigma_z) * (x < xR) * 0.5 * rotor_diameter_i
 
-            r, C = rC(
+            r_squared, C = rC(
                 wind_veer,
                 sigma_y,
                 sigma_z,
@@ -146,7 +146,7 @@ def function(
                 rotor_diameter_i,
             )
 
-            near_wake_deficit = gaussian_function(C, r, 1, np.sqrt(0.5))
+            near_wake_deficit = gaussian_function(C, r_squared, 1, np.sqrt(0.5))
             near_wake_deficit *= near_wake_mask
 
             velocity_deficit += near_wake_deficit
@@ -160,7 +160,7 @@ def function(
             sigma_y = (ky * (x - x0) + sigma_y0) * far_wake_mask + sigma_y0 * (x < x0)
             sigma_z = (kz * (x - x0) + sigma_z0) * far_wake_mask + sigma_z0 * (x < x0)
 
-            r, C = rC(
+            r_squared, C = rC(
                 wind_veer,
                 sigma_y,
                 sigma_z,
@@ -174,7 +174,7 @@ def function(
                 rotor_diameter_i,
             )
 
-            far_wake_deficit = gaussian_function(C, r, 1, np.sqrt(0.5))
+            far_wake_deficit = gaussian_function(C, r_squared, 1, np.sqrt(0.5))
             far_wake_deficit *= far_wake_mask
 
             velocity_deficit += far_wake_deficit
@@ -189,7 +189,7 @@ def rC(wind_veer, sigma_y, sigma_z, y, y_i, delta, z, HH, Ct, yaw, D):
     # a = cosd(wind_veer) ** 2 / (2 * sigma_y ** 2) + sind(wind_veer) ** 2 / (2 * sigma_z ** 2)
     # b = -sind(2 * wind_veer) / (4 * sigma_y ** 2) + sind(2 * wind_veer) / (4 * sigma_z ** 2)
     # c = sind(wind_veer) ** 2 / (2 * sigma_y ** 2) + cosd(wind_veer) ** 2 / (2 * sigma_z ** 2)
-    # r = (
+    # r_squared = (
     #     a * (y - y_i - delta) ** 2
     #     - 2 * b * (y - y_i - delta) * (z - HH)
     #     + c * (z - HH) ** 2
@@ -204,7 +204,7 @@ def rC(wind_veer, sigma_y, sigma_z, y, y_i, delta, z, HH, Ct, yaw, D):
     # c = sind(wind_veer) ** 2 / (twox_sigmay_2) + cosd(wind_veer) ** 2 / (twox_sigmaz_2)
     # delta_y = y - y_i - delta
     # delta_z = z - HH
-    # r = (a * (delta_y ** 2) - 2 * b * (delta_y) * (delta_z) + c * (delta_z ** 2))
+    # r_squared = (a * (delta_y ** 2) - 2 * b * (delta_y) * (delta_z) + c * (delta_z ** 2))
     # C = 1 - np.sqrt(np.clip(1 - (Ct * cosd(yaw) / (8.0 * sigma_y * sigma_z / (D * D))), 0.0, 1.0))
 
     ## Numexpr
@@ -218,12 +218,12 @@ def rC(wind_veer, sigma_y, sigma_z, y, y_i, delta, z, HH, Ct, yaw, D):
     c = ne.evaluate(
         "sin(wind_veer) ** 2 / (2 * sigma_y ** 2) + cos(wind_veer) ** 2 / (2 * sigma_z ** 2)"
     )
-    r = ne.evaluate(
+    r_squared = ne.evaluate(
         "a * ((y - y_i - delta) ** 2) - 2 * b * (y - y_i - delta) * (z - HH) + c * ((z - HH) ** 2)"
     )
     d = np.clip(1 - (Ct * cosd(yaw) / ( 8.0 * sigma_y * sigma_z / (D * D) )), 0.0, 1.0)
     C = ne.evaluate("1 - sqrt(d)")
-    return r, C
+    return r_squared, C
 
 
 def mask_upstream_wake(mesh_y_rotated, x_coord_rotated, y_coord_rotated, turbine_yaw):
@@ -232,6 +232,6 @@ def mask_upstream_wake(mesh_y_rotated, x_coord_rotated, y_coord_rotated, turbine
     return xR, yR
 
 
-def gaussian_function(C, r, n, sigma):
-    result = ne.evaluate("C * exp(-1 * r ** n / (2 * sigma ** 2))")
+def gaussian_function(C, r_squared, n, sigma):
+    result = ne.evaluate("C * exp(-1 * r_squared ** n / (2 * sigma ** 2))")
     return result