Convert turbulence intensity from single value to n_findex length arr…

…ay (NREL#782) * Convert turbulence_intensity to tubulence_intensities throughout the code and refactor all code to expect turbulence_intensities to be an array and not a float * Add additional tests of turbulence intensity to confirm correct behavior of new features * Complete WindRose and TimeSeries handling of turbulence intensities * Add helper functions to WindRose and TimeSeries which allow turbulence intensities to be generated, rather than provided, as a function of wind directions and wind speeds * Add additional examples of usage --------- Co-authored-by: misi9170 <[email protected]> Co-authored-by: Rafael M Mudafort <[email protected]> Co-authored-by: Eric Simley <[email protected]>
rafmudaf · Feb 3, 2024 · 61e1f13 · 61e1f13
1 parent 1eb76b1
commit 61e1f13
Show file tree

Hide file tree

Showing 36 changed files with 584 additions and 105 deletions.
diff --git a/examples/12_optimize_yaw_in_parallel.py b/examples/12_optimize_yaw_in_parallel.py
@@ -63,7 +63,7 @@ def load_windrose():
     fi_aep.reinitialize(
         wind_directions=wind_directions,
         wind_speeds=wind_speeds,
-        turbulence_intensity=0.08  # Assume 8% turbulence intensity
+        turbulence_intensities=[0.08],  # Assume 8% turbulence intensity
     )
 
     # Pour this into a parallel computing interface
@@ -105,7 +105,7 @@ def load_windrose():
     fi_opt.reinitialize(
         wind_directions=wind_directions,
         wind_speeds=wind_speeds,
-        turbulence_intensity=0.08  # Assume 8% turbulence intensity
+        turbulence_intensities=[0.08],  # Assume 8% turbulence intensity
     )
 
     # Pour this into a parallel computing interface

diff --git a/examples/19_streamlit_demo.py b/examples/19_streamlit_demo.py
@@ -124,7 +124,7 @@
         layout_y=Y,
         wind_speeds=[wind_speed],
         wind_directions=[wind_direction],
-        turbulence_intensity=turbulence_intensity
+        turbulence_intensities=[turbulence_intensity],
     )
 
     fi.calculate_wake(yaw_angles=yaw_angles_base)
@@ -168,7 +168,7 @@
         layout_y=Y,
         wind_speeds=[wind_speed],
         wind_directions=[wind_direction],
-        turbulence_intensity=turbulence_intensity
+        turbulence_intensities=[turbulence_intensity],
     )
 
     fi.calculate_wake(yaw_angles=yaw_angles_yaw)

diff --git a/examples/34_wind_data.py b/examples/34_wind_data.py
@@ -50,7 +50,7 @@
 
 
 # Build the time series
-time_series = TimeSeries(wd_array, ws_array)  # , turbulence_intensity=ti_array)
+time_series = TimeSeries(wd_array, ws_array, turbulence_intensities=ti_array)
 
 # Now build the wind rose
 wind_rose = time_series.to_wind_rose()

diff --git a/examples/35_sweep_ti.py b/examples/35_sweep_ti.py
@@ -0,0 +1,62 @@
+# Copyright 2024 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
+
+from floris.tools import (
+    FlorisInterface,
+    TimeSeries,
+    WindRose,
+)
+from floris.utilities import wrap_360
+
+
+"""
+Demonstrate the new behavior in V4 where TI is an array rather than a float.
+Set up an array of two turbines and sweep TI while holding wd/ws constant.
+Use the TimeSeries object to drive the FLORIS calculations.
+"""
+
+
+# Generate a random time series of wind speeds, wind directions and turbulence intensities
+N = 50
+wd_array = 270.0 * np.ones(N)
+ws_array = 8.0 * np.ones(N)
+ti_array = np.linspace(0.03, 0.2, N)
+
+
+# Build the time series
+time_series = TimeSeries(wd_array, ws_array, turbulence_intensities=ti_array)
+
+
+# Now set up a FLORIS model and initialize it using the time
+fi = FlorisInterface("inputs/gch.yaml")
+fi.reinitialize(layout_x=[0, 500.0], layout_y=[0.0, 0.0], wind_data=time_series)
+fi.calculate_wake()
+turbine_power = fi.get_turbine_powers()
+
+fig, axarr = plt.subplots(2, 1, sharex=True, figsize=(6, 6))
+ax = axarr[0]
+ax.plot(ti_array*100, turbine_power[:, 0]/1000, color="k")
+ax.set_ylabel("Front turbine power [kW]")
+ax = axarr[1]
+ax.plot(ti_array*100, turbine_power[:, 1]/1000, color="k")
+ax.set_ylabel("Rear turbine power [kW]")
+ax.set_xlabel("Turbulence intensity [%]")
+
+for ax in axarr:
+    ax.grid(True)
+
+plt.show()
diff --git a/examples/36_generate_ti.py b/examples/36_generate_ti.py
@@ -0,0 +1,82 @@
+# Copyright 2024 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
+
+from floris.tools import (
+    FlorisInterface,
+    TimeSeries,
+    WindRose,
+)
+from floris.utilities import wrap_360
+
+
+"""
+Demonstrate usage of TI generating and plotting functionality in the WindRose
+and TimeSeries classes
+"""
+
+
+# Generate a random time series of wind speeds, wind directions and turbulence intensities
+wind_directions = np.array([250, 260, 270])
+wind_speeds = np.array([5, 6, 7, 8, 9, 10])
+
+# Declare a WindRose object
+wind_rose = WindRose(wind_directions=wind_directions, wind_speeds=wind_speeds)
+
+
+# Define a custom function where TI = 1 / wind_speed
+def custom_ti_func(wind_directions, wind_speeds):
+    return 1 / wind_speeds
+
+
+wind_rose.assign_ti_using_wd_ws_function(custom_ti_func)
+
+fig, ax = plt.subplots()
+wind_rose.plot_ti_over_ws(ax)
+ax.set_title("Turbulence Intensity defined by custom function")
+
+# Now use the normal turbulence model approach from the IEC 61400-1 standard,
+# wherein TI is defined as a function of wind speed:
+# Iref is defined as the TI value at 15 m/s. Note that Iref = 0.07 is lower
+# than the values of Iref used in the IEC standard, but produces TI values more
+# in line with those typically used in FLORIS (TI=8.6% at 8 m/s).
+Iref = 0.07
+wind_rose.assign_ti_using_IEC_method(Iref)
+fig, ax = plt.subplots()
+wind_rose.plot_ti_over_ws(ax)
+ax.set_title(f"Turbulence Intensity defined by Iref = {Iref:0.2}")
+
+
+# Demonstrate equivalent usage in time series
+N = 100
+wind_directions = 270 * np.ones(N)
+wind_speeds = np.linspace(5, 15, N)
+time_series = TimeSeries(wind_directions=wind_directions, wind_speeds=wind_speeds)
+time_series.assign_ti_using_IEC_method(Iref=Iref)
+
+fig, axarr = plt.subplots(2, 1, sharex=True, figsize=(7, 8))
+ax = axarr[0]
+ax.plot(wind_speeds)
+ax.set_ylabel("Wind Speeds (m/s)")
+ax.grid(True)
+ax = axarr[1]
+ax.plot(time_series.turbulence_intensities)
+ax.set_ylabel("Turbulence Intensity (-)")
+ax.grid(True)
+fig.suptitle("Generating TI in TimeSeries")
+
+
+plt.show()
diff --git a/examples/inputs/cc.yaml b/examples/inputs/cc.yaml
@@ -30,7 +30,8 @@ farm:
 flow_field:
   air_density: 1.225
   reference_wind_height: -1 # -1 is code for use the hub height
-  turbulence_intensity: 0.06
+  turbulence_intensities:
+  - 0.06
   wind_directions:
   - 270.0
   wind_shear: 0.12

diff --git a/examples/inputs/emgauss.yaml b/examples/inputs/emgauss.yaml
@@ -30,7 +30,8 @@ farm:
 flow_field:
   air_density: 1.225
   reference_wind_height: -1 # -1 is code for use the hub height
-  turbulence_intensity: 0.06
+  turbulence_intensities:
+  - 0.06
   wind_directions:
   - 270.0
   wind_shear: 0.12

diff --git a/examples/inputs/gch.yaml b/examples/inputs/gch.yaml
@@ -112,7 +112,8 @@ flow_field:
 
   ###
   # The level of turbulence intensity level in the wind.
-  turbulence_intensity: 0.06
+  turbulence_intensities:
+  - 0.06
 
   ###
   # The wind directions to include in the simulation.

diff --git a/examples/inputs/gch_heterogeneous_inflow.yaml b/examples/inputs/gch_heterogeneous_inflow.yaml
@@ -44,7 +44,8 @@ flow_field:
       - -300.
       - 300.
   reference_wind_height: -1
-  turbulence_intensity: 0.06
+  turbulence_intensities:
+  - 0.06
   wind_directions:
   - 270.0
   wind_shear: 0.12

diff --git a/examples/inputs/gch_multi_dim_cp_ct.yaml b/examples/inputs/gch_multi_dim_cp_ct.yaml
@@ -33,7 +33,8 @@ flow_field:
     Hs: 3.01
   air_density: 1.225
   reference_wind_height: -1 # -1 is code for use the hub height
-  turbulence_intensity: 0.06
+  turbulence_intensities:
+  - 0.06
   wind_directions:
   - 270.0
   wind_shear: 0.12

diff --git a/examples/inputs/gch_multiple_turbine_types.yaml b/examples/inputs/gch_multiple_turbine_types.yaml
@@ -29,7 +29,8 @@ farm:
 flow_field:
   air_density: 1.225
   reference_wind_height: 90.0 # Since multiple defined turbines, must specify explicitly the reference wind height
-  turbulence_intensity: 0.06
+  turbulence_intensities:
+  - 0.06
   wind_directions:
   - 270.0
   wind_shear: 0.12

diff --git a/examples/inputs/jensen.yaml b/examples/inputs/jensen.yaml
@@ -30,7 +30,8 @@ farm:
 flow_field:
   air_density: 1.225
   reference_wind_height: -1 # -1 is code for use the hub height
-  turbulence_intensity: 0.06
+  turbulence_intensities:
+  - 0.06
   wind_directions:
   - 270.0
   wind_shear: 0.12

diff --git a/examples/inputs/turbopark.yaml b/examples/inputs/turbopark.yaml
@@ -30,7 +30,8 @@ farm:
 flow_field:
   air_density: 1.225
   reference_wind_height: 90.0
-  turbulence_intensity: 0.06
+  turbulence_intensities:
+  - 0.06
   wind_directions:
   - 270.0
   wind_shear: 0.12

diff --git a/examples/inputs_floating/emgauss_fixed.yaml b/examples/inputs_floating/emgauss_fixed.yaml
@@ -30,7 +30,8 @@ farm:
 flow_field:
   air_density: 1.225
   reference_wind_height: -1 # -1 is code for use the hub height
-  turbulence_intensity: 0.06
+  turbulence_intensities:
+  - 0.06
   wind_directions:
   - 270.0
   wind_shear: 0.12

diff --git a/examples/inputs_floating/emgauss_floating.yaml b/examples/inputs_floating/emgauss_floating.yaml
@@ -30,7 +30,8 @@ farm:
 flow_field:
   air_density: 1.225
   reference_wind_height: -1 # -1 is code for use the hub height
-  turbulence_intensity: 0.06
+  turbulence_intensities:
+  - 0.06
   wind_directions:
   - 270.0
   wind_shear: 0.12

diff --git a/examples/inputs_floating/emgauss_floating_fixedtilt15.yaml b/examples/inputs_floating/emgauss_floating_fixedtilt15.yaml
@@ -26,7 +26,8 @@ farm:
 flow_field:
   air_density: 1.225
   reference_wind_height: -1 # -1 is code for use the hub height
-  turbulence_intensity: 0.06
+  turbulence_intensities:
+  - 0.06
   wind_directions:
   - 270.0
   wind_shear: 0.12

diff --git a/examples/inputs_floating/emgauss_floating_fixedtilt5.yaml b/examples/inputs_floating/emgauss_floating_fixedtilt5.yaml
@@ -26,7 +26,8 @@ farm:
 flow_field:
   air_density: 1.225
   reference_wind_height: -1 # -1 is code for use the hub height
-  turbulence_intensity: 0.06
+  turbulence_intensities:
+  - 0.06
   wind_directions:
   - 270.0
   wind_shear: 0.12

diff --git a/examples/inputs_floating/gch_fixed.yaml b/examples/inputs_floating/gch_fixed.yaml
@@ -26,7 +26,8 @@ farm:
 flow_field:
   air_density: 1.225
   reference_wind_height: -1
-  turbulence_intensity: 0.06
+  turbulence_intensities:
+  - 0.06
   wind_directions:
   - 270.0
   wind_shear: 0.12

diff --git a/examples/inputs_floating/gch_floating.yaml b/examples/inputs_floating/gch_floating.yaml
@@ -27,7 +27,8 @@ farm:
 flow_field:
   air_density: 1.225
   reference_wind_height: -1
-  turbulence_intensity: 0.06
+  turbulence_intensities:
+  - 0.06
   wind_directions:
   - 270.0
   wind_shear: 0.12

diff --git a/examples/inputs_floating/gch_floating_defined_floating.yaml b/examples/inputs_floating/gch_floating_defined_floating.yaml
@@ -26,7 +26,8 @@ farm:
 flow_field:
   air_density: 1.225
   reference_wind_height: -1
-  turbulence_intensity: 0.06
+  turbulence_intensities:
+  - 0.06
   wind_directions:
   - 270.0
   wind_shear: 0.12

diff --git a/floris/simulation/flow_field.py b/floris/simulation/flow_field.py
@@ -39,7 +39,7 @@ class FlowField(BaseClass):
     wind_veer: float = field(converter=float)
     wind_shear: float = field(converter=float)
     air_density: float = field(converter=float)
-    turbulence_intensity: float = field(converter=float)
+    turbulence_intensities: NDArrayFloat = field(converter=floris_array_converter)
     reference_wind_height: float = field(converter=float)
     time_series: bool = field(default=False)
     heterogenous_inflow_config: dict = field(default=None)
@@ -66,6 +66,17 @@ class FlowField(BaseClass):
         init=False, factory=lambda: np.array([])
     )
 
+    @turbulence_intensities.validator
+    def turbulence_intensities_validator(
+        self, instance: attrs.Attribute, value: NDArrayFloat
+    ) -> None:
+
+        # Check the turbulence intensity is either length 1 or n_findex
+        if len(value) != 1 and len(value) != self.n_findex:
+            raise ValueError("turbulence_intensities should either be length 1 or n_findex")
+
+
+
     @wind_directions.validator
     def wind_directions_validator(self, instance: attrs.Attribute, value: NDArrayFloat) -> None:
         """Using the validator method to keep the `n_findex` attribute up to date."""
@@ -108,6 +119,10 @@ def __attrs_post_init__(self) -> None:
         if self.heterogenous_inflow_config is not None:
             self.generate_heterogeneous_wind_map()
 
+        # If turbulence_intensity is length 1, then convert it to a uniform array of
+        # length n_findex
+        if len(self.turbulence_intensities) == 1:
+            self.turbulence_intensities = self.turbulence_intensities[0] * np.ones(self.n_findex)
 
     def initialize_velocity_field(self, grid: Grid) -> None:
 
@@ -197,14 +212,13 @@ def initialize_velocity_field(self, grid: Grid) -> None:
         self.v_sorted = self.v_initial_sorted.copy()
         self.w_sorted = self.w_initial_sorted.copy()
 
-        self.turbulence_intensity_field = self.turbulence_intensity * np.ones(
-            (
-                self.n_findex,
-                grid.n_turbines,
-                1,
-                1,
-            )
+        self.turbulence_intensity_field = self.turbulence_intensities[:, None, None, None]
+        self.turbulence_intensity_field = np.repeat(
+            self.turbulence_intensity_field,
+            grid.n_turbines,
+            axis=1
         )
+
         self.turbulence_intensity_field_sorted = self.turbulence_intensity_field.copy()
 
     def finalize(self, unsorted_indices):