Add utility to build turbine yaml from absolute power, thrust curves (#…

…729)

* Structure of file built out.

* Main functionality built out.

* Test built.

* passes tests.

* isort, ruff.

* Issues with 3.8, 3.9 tests.

* isort issue when removing import.

* newline?

* Reverting; __future__ import to handle testing on 3.8, 3.9

* Now defaults to providing dictionary back; option to print to yaml.

* Test updated; ruff, isort.

* Improving documentation and exposition of function arguments.'

* Update comments, syntax, and typos

* Ensure path independence in unit tests

---------

Co-authored-by: Rafael M Mudafort <[email protected]>

docs/examples.md

@@ -172,6 +172,12 @@ This example builds on example 30. Specifically, fictional data for varying Cp/C
 wave period, Ts, and wave height, Hs, is used to show the difference in power performance for
 different wave heights.
 
+### 32_specify_turbine_power_curve.py
+
+This example demonstrates how to generate a turbine dictionary or yaml input file based on
+a specified power and thrust curve. The power and thrust curves may be specified as power
+and thrust coefficients or as absolute values.
+
 ## Optimization
 
 These examples demonstrate use of the optimization routines

examples/32_specify_turbine_power_curve.py

@@ -0,0 +1,78 @@
+# Copyright 2021 NREL
+
+# Licensed under the Apache License, Version 2.0 (the "License"); you may not
+# use this file except in compliance with the License. You may obtain a copy of
+# the License at http://www.apache.org/licenses/LICENSE-2.0
+
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS, WITHOUT
+# WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the
+# License for the specific language governing permissions and limitations under
+# the License.
+
+# See https://floris.readthedocs.io for documentation
+
+
+import matplotlib.pyplot as plt
+import numpy as np
+
+import floris.tools.visualization as wakeviz
+from floris.tools import FlorisInterface
+from floris.turbine_library.turbine_utilities import build_turbine_dict
+
+
+"""
+This example demonstrates how to specify a turbine model based on a power
+and thrust curve for the wind turbine, as well as possible physical parameters
+(which default to the parameters of the NREL 5MW reference turbine).
+
+Note that it is also possible to have a .yaml created, if the file_path
+argument to build_turbine_dict is set.
+"""
+
+# Generate an example turbine power and thrust curve for use in the FLORIS model
+turbine_data_dict = {
+    "wind_speed":[0, 2, 4, 6, 8, 10, 12, 14, 16, 18, 20],
+    "power_absolute":[0, 30, 200, 500, 1000, 2000, 4000, 4000, 4000, 4000, 4000],
+    "thrust_coefficient":[0, 0.9, 0.9, 0.8, 0.7, 0.6, 0.5, 0.4, 0.3, 0.25, 0.2]
+}
+
+turbine_dict = build_turbine_dict(
+    turbine_data_dict,
+    "example_turbine",
+    file_path=None,
+    generator_efficiency=1,
+    hub_height=90,
+    pP=1.88,
+    pT=1.88,
+    rotor_diameter=126,
+    TSR=8,
+    air_density=1.225,
+    ref_tilt_cp_ct=5
+)
+
+fi = FlorisInterface("inputs/gch.yaml")
+wind_speeds = np.linspace(1, 15, 100)
+# Replace the turbine(s) in the FLORIS model with the created one
+fi.reinitialize(
+    layout_x=[0],
+    layout_y=[0],
+    wind_speeds=wind_speeds,
+    turbine_type=[turbine_dict]
+)
+fi.calculate_wake()
+
+powers = fi.get_farm_power()
+
+fig, ax = plt.subplots(1,1)
+
+ax.scatter(wind_speeds, powers[0,:]/1000, color="C0", s=5, label="Test points")
+ax.scatter(turbine_data_dict["wind_speed"], turbine_data_dict["power_absolute"],
+    color="red", s=20, label="Specified points")
+
+ax.grid()
+ax.set_xlabel("Wind speed [m/s]")
+ax.set_ylabel("Power [kW]")
+ax.legend()
+
+plt.show()

floris/turbine_library/__init__.py

@@ -1 +1,2 @@
 from floris.turbine_library.turbine_previewer import TurbineInterface, TurbineLibrary
+from floris.turbine_library.turbine_utilities import build_turbine_dict

floris/turbine_library/turbine_previewer.py

@@ -12,6 +12,8 @@
 
 # See https://floris.readthedocs.io for documentation
 
+from __future__ import annotations
+
 from pathlib import Path
 
 import attrs

floris/turbine_library/turbine_utilities.py

@@ -0,0 +1,166 @@
+import os.path
+
+import numpy as np
+import yaml
+
+
+def build_turbine_dict(
+    turbine_data_dict,
+    turbine_name,
+    file_path=None,
+    generator_efficiency=1.0,
+    hub_height=90.0,
+    pP=1.88,
+    pT=1.88,
+    rotor_diameter=126.0,
+    TSR=8.0,
+    air_density=1.225,
+    ref_tilt_cp_ct=5.0
+):
+    """
+    Tool for formatting a full turbine dict from data formatted as a
+    dictionary.
+
+    Default value for turbine physical parameters are from the NREL 5MW reference
+    wind turbine.
+
+    Returns a turbine dictionary object as expected by FLORIS. Optionally,
+    prints the dictionary to a yaml to be included in a FLORIS wake model yaml.
+
+    turbine_data is a dictionary that contains keys specifying the
+    turbine power and thrust as a function of wind speed. The following keys
+    are possible:
+    - wind_speed [m/s]
+    - power_absolute [kW]
+    - power_coefficient [-]
+    - thrust_absolute [kN]
+    - thrust_coefficient [-]
+    Of these, wind_speed is required. One of power_absolute and power_coefficient
+    must be specified; and one of thrust_absolute and thrust_coefficient must be
+    specified. If both _absolute and _coefficient versions are specified, the
+    _coefficient entry will be used and the _absolute entry ignored.
+
+    Args:
+        turbine_data_dict (dict): Dictionary containing performance of the wind
+            turbine as a function of wind speed. Described in more detail above.
+        turbine_name (string): Name of the turbine, which will be used for the
+            turbine_type field as well as the filename.
+        file_path (): Path for placement of the produced yaml. Defaults to None,
+            in which case no yaml is written.
+        generator_efficiency (float): Generator efficiency [-]. Defaults to 1.0.
+        hub_height (float): Hub height [m]. Defaults to 90.0.
+        pP (float): Cosine exponent for power loss to yaw [-]. Defaults to 1.88.
+        pT (float): Cosine exponent for thrust loss to yaw [-]. Defaults to 1.88.
+        rotor_diameter (float). Rotor diameter [m]. Defaults to 126.0.
+        TSR (float). Turbine optimal tip-speed ratio [-]. Defaults to 8.0.
+        air_density (float). Air density used to specify power and thrust
+            curves [kg/m^3]. Defaults to 1.225.
+        ref_tilt_cp_ct (float). Rotor tilt (due to shaft tilt and/or platform
+            tilt) used when defining the power and thrust curves [deg]. Defaults
+            to 5.0.
+
+    Returns:
+        turbine_dict (dict): Formatted turbine dictionary as expected by FLORIS.
+    """
+
+    # Check that necessary columns are specified
+    if "wind_speed" not in turbine_data_dict:
+        raise KeyError("wind_speed column must be specified.")
+    u = np.array(turbine_data_dict["wind_speed"])
+    A = np.pi * rotor_diameter**2/4
+
+    # Construct the Cp curve
+    if "power_coefficient" in turbine_data_dict:
+        if "power_absolute" in turbine_data_dict:
+            print(
+                "Found both power_absolute and power_coefficient."
+                "Ignoring power_absolute."
+            )
+        Cp = np.array(turbine_data_dict["power_coefficient"])
+
+    elif "power_absolute" in turbine_data_dict:
+        P = np.array(turbine_data_dict["power_absolute"])
+        if _find_nearest_value_for_wind_speed(P, u, 10) > 20000 or \
+           _find_nearest_value_for_wind_speed(P, u, 10) < 1000:
+           print(
+               "Unusual power value detected. Please check that power_absolute",
+               "is specified in kW."
+           )
+
+        validity_mask = (P != 0) | (u != 0)
+        Cp = np.zeros_like(P, dtype=float)
+
+        Cp[validity_mask] = (P[validity_mask]*1000) / \
+                            (0.5*air_density*A*u[validity_mask]**3)
+
+    else:
+        raise KeyError(
+            "Either power_absolute or power_coefficient must be specified."
+        )
+
+    # Construct Ct curve
+    if "thrust_coefficient" in turbine_data_dict:
+        if "thrust_absolute" in turbine_data_dict:
+            print(
+                "Found both thrust_absolute and thrust_coefficient."
+                "Ignoring thrust_absolute."
+            )
+        Ct = np.array(turbine_data_dict["thrust_coefficient"])
+
+    elif "thrust_absolute" in turbine_data_dict:
+        T = np.array(turbine_data_dict["thrust_absolute"])
+        if _find_nearest_value_for_wind_speed(T, u, 10) > 3000 or \
+           _find_nearest_value_for_wind_speed(T, u, 10) < 100:
+           print(
+               "Unusual thrust value detected. Please check that thrust_absolute",
+               "is specified in kN."
+           )
+
+        validity_mask = (T != 0) | (u != 0)
+        Ct = np.zeros_like(T)
+
+        Ct[validity_mask] = (T[validity_mask]*1000)/\
+                            (0.5*air_density*A*u[validity_mask]**2)
+
+    else:
+        raise KeyError(
+            "Either thrust_absolute or thrust_coefficient must be specified."
+        )
+
+    # Build the turbine dict
+    power_thrust_dict = {
+        "wind_speed": u.tolist(),
+        "power": Cp.tolist(),
+        "thrust": Ct.tolist()
+    }
+
+    turbine_dict = {
+        "turbine_type": turbine_name,
+        "generator_efficiency": generator_efficiency,
+        "hub_height": hub_height,
+        "pP": pP,
+        "pT": pT,
+        "rotor_diameter": rotor_diameter,
+        "TSR": TSR,
+        "ref_density_cp_ct": air_density,
+        "ref_tilt_cp_ct": ref_tilt_cp_ct,
+        "power_thrust_table": power_thrust_dict
+    }
+
+    # Create yaml file
+    if file_path is not None:
+        full_name = os.path.join(file_path, turbine_name+".yaml")
+        yaml.dump(
+            turbine_dict,
+            open(full_name, "w"),
+            sort_keys=False
+        )
+
+        print(full_name, "created.")
+
+    return turbine_dict
+
+def _find_nearest_value_for_wind_speed(test_vals, ws_vals, ws):
+    errs = np.absolute(ws_vals-ws)
+    idx = errs.argmin()
+    return test_vals[idx]

tests/turbine_unit_test.py

@@ -13,9 +13,13 @@
 # See https://floris.readthedocs.io for documentation
 
 
+import os
+from pathlib import Path
+
 import attr
 import numpy as np
 import pytest
+import yaml
 from scipy.interpolate import interp1d
 
 from floris.simulation import (
@@ -31,6 +35,7 @@
     compute_tilt_angles_for_floating_turbines,
     PowerThrustTable,
 )
+from floris.turbine_library import build_turbine_dict
 from tests.conftest import SampleInputs, WIND_SPEEDS
 
 
@@ -563,3 +568,88 @@ def test_asdict(sample_inputs_fixture: SampleInputs):
     dict2 = new_turb.as_dict()
 
     assert dict1 == dict2
+
+def test_build_turbine_dict():
+
+    orig_file_path = Path(__file__).resolve().parent / "data" / "nrel_5MW_custom.yaml"
+    test_turb_name = "test_turbine_export"
+    test_file_path = "."
+
+    in_dict = yaml.safe_load( open(orig_file_path, "r") )
+
+    # Mocked up turbine data
+    turbine_data_dict = {
+        "wind_speed":in_dict["power_thrust_table"]["wind_speed"],
+        "power_coefficient":in_dict["power_thrust_table"]["power"],
+        "thrust_coefficient":in_dict["power_thrust_table"]["thrust"]
+    }
+
+    build_turbine_dict(
+        turbine_data_dict,
+        test_turb_name,
+        file_path=test_file_path,
+        generator_efficiency=in_dict["generator_efficiency"],
+        hub_height=in_dict["hub_height"],
+        pP=in_dict["pP"],
+        pT=in_dict["pT"],
+        rotor_diameter=in_dict["rotor_diameter"],
+        TSR=in_dict["TSR"],
+        air_density=in_dict["ref_density_cp_ct"],
+        ref_tilt_cp_ct=in_dict["ref_tilt_cp_ct"]
+    )
+
+    test_dict = yaml.safe_load(
+        open(os.path.join(test_file_path, test_turb_name+".yaml"), "r")
+    )
+
+    # Correct intended difference for test; assert equal
+    test_dict["turbine_type"] = in_dict["turbine_type"]
+    assert list(in_dict.keys()) == list(test_dict.keys())
+    assert in_dict == test_dict
+
+    # Now, in absolute values
+    Cp = np.array(in_dict["power_thrust_table"]["power"])
+    Ct = np.array(in_dict["power_thrust_table"]["thrust"])
+    ws = np.array(in_dict["power_thrust_table"]["wind_speed"])
+
+    P = 0.5 * in_dict["ref_density_cp_ct"] * (np.pi * in_dict["rotor_diameter"]**2/4) \
+        * Cp * ws**3
+    T = 0.5 * in_dict["ref_density_cp_ct"] * (np.pi * in_dict["rotor_diameter"]**2/4) \
+        * Ct * ws**2
+
+    turbine_data_dict = {
+        "wind_speed":in_dict["power_thrust_table"]["wind_speed"],
+        "power_absolute": P/1000,
+        "thrust_absolute": T/1000
+    }
+
+    build_turbine_dict(
+        turbine_data_dict,
+        test_turb_name,
+        file_path=test_file_path,
+        generator_efficiency=in_dict["generator_efficiency"],
+        hub_height=in_dict["hub_height"],
+        pP=in_dict["pP"],
+        pT=in_dict["pT"],
+        rotor_diameter=in_dict["rotor_diameter"],
+        TSR=in_dict["TSR"],
+        air_density=in_dict["ref_density_cp_ct"],
+        ref_tilt_cp_ct=in_dict["ref_tilt_cp_ct"]
+    )
+
+    test_dict = yaml.safe_load(
+        open(os.path.join(test_file_path, test_turb_name+".yaml"), "r")
+    )
+
+    test_dict["turbine_type"] = in_dict["turbine_type"]
+    assert list(in_dict.keys()) == list(test_dict.keys())
+    for k in in_dict.keys():
+        if type(in_dict[k]) is dict:
+            for k2 in in_dict[k].keys():
+                assert np.allclose(in_dict[k][k2], test_dict[k][k2])
+        elif type(in_dict[k]) is str:
+            assert in_dict[k] == test_dict[k]
+        else:
+            assert np.allclose(in_dict[k], test_dict[k])
+
+    os.remove( os.path.join(test_file_path, test_turb_name+".yaml") )