Feature: Add multidimensional Cp/Ct turbine definition support (#711)

* Add multi-dimensional Cp/Ct turbine support

* Update Turbine to take None/[ ] values for Cp/Ct variables

* Add necessary farm expansions for multidim Cp/Ct

* Add multidim specific solver and wrap gauss wake model

* Add multidim conditions to FlowField

* Add get_turbine_powers_multidim to FlorisInterface

* Add flatten_dict to required pacakges

* Updating to support the full wd and ws dimensions

* Updated scope of multidim functions

* Adding example for multi-dimensional Cp/Ct turbines

* Adding test for multi-dimensional class and functions

* Add an example that uses two different wave height

* Updated examples and documentation

* Remove multi-dim turbines from example 18

* Updating input reference guide for multi-dimensional Cp/Ct data

---------

Co-authored-by: Paul <paul.fleming@nrel.gov>

docs/examples.md

@@ -158,6 +158,20 @@ For further examples on floating wind turbines, see also examples
 24 (effects of tilt on turbine power and thrust coefficients) and 25
 (vertical wake deflection by a forced tilt angle).
 
+### 30_multi_dimensional_cp_ct.py
+
+This example showcases the capability of using multi-dimensional Cp/Ct data in turbine defintions
+dependent on external conditions. Specifically, fictional data for varying Cp/Ct values based on
+wave period, Ts, and wave height, Hs, is used, showing the user how to setup the turbine
+definition and input file. Also demonstrated is the different method for getting turbine
+powers when using multi-dimensional Cp/Ct data.
+
+### 31_multi_dimensional_cp_ct_2Hs.py
+
+This example builds on example 30. Specifically, fictional data for varying Cp/Ct values based on
+wave period, Ts, and wave height, Hs, is used to show the difference in power performance for
+different wave heights.
+
 ## Optimization
 
 These examples demonstrate use of the optimization routines

examples/01_opening_floris_computing_power.py

@@ -51,7 +51,7 @@
 print(turbine_powers)
 print("Shape: ",turbine_powers.shape)
 
-# Single wind speed and wind direction
+# Single wind speed and multiple wind directions
 print('\n========================= Single Wind Direction and Multiple Wind Speeds ===============')
 
 
@@ -64,7 +64,7 @@
 print(turbine_powers)
 print("Shape: ",turbine_powers.shape)
 
-# Single wind speed and wind direction
+# Multiple wind speeds and multiple wind directions
 print('\n========================= Multiple Wind Directions and Multiple Wind Speeds ============')
 
 wind_directions = np.array([260., 270., 280.])

examples/18_check_turbine.py

@@ -42,6 +42,8 @@
 turbines = os.listdir('../floris/turbine_library')
 turbines = [t for t in turbines if 'yaml' in t]
 turbines = [t.strip('.yaml') for t in turbines]
+# Remove multi-dimensional Cp/Ct turbine definitions as they require different handling
+turbines = [i for i in turbines if ('multi_dim' not in i)]
 
 # Declare a set of figures for comparing cp and ct across models
 fig_cp_ct, axarr_cp_ct = plt.subplots(2,1,sharex=True,figsize=(10,10))

examples/30_multi_dimensional_cp_ct.py

@@ -0,0 +1,104 @@
+# Copyright 2023 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 numpy as np
+
+from floris.tools import FlorisInterface
+
+
+"""
+This example follows the same setup as example 01 to createa a FLORIS instance and:
+1) Makes a two-turbine layout
+2) Demonstrates single ws/wd simulations
+3) Demonstrates mulitple ws/wd simulations
+
+with the modification of using a turbine definition that has a multi-dimensional Cp/Ct table.
+
+In the input file `gch_multi_dim_cp_ct.yaml`, the turbine_type points to a turbine definition,
+iea_15MW_floating_multi_dim_cp_ct.yaml located in the turbine_library,
+that supplies a multi-dimensional Cp/Ct data file in the form of a .csv file. This .csv file
+contains two additional conditions to define Cp/Ct values for: Tp for wave period, and Hs for wave
+height. For every combination of Tp and Hs defined, a Cp/Ct/Wind speed table of values is also
+defined. It is required for this .csv file to have the last 3 columns be ws, Cp, and Ct. In order
+for this table to be used, the flag 'multi_dimensional_cp_ct' must be present and set to true in
+the turbine definition. Also of note is the 'velocity_model' must be set to 'multidim_cp_ct' in
+the main input file. With both of these values provided, the solver will downselect to use the
+interpolant defined at the closest conditions. The user must supply these conditions in the
+main input file under the 'flow_field' section, e.g.:
+
+NOTE: The multi-dimensional Cp/Ct data used in this example is fictional for the purposes of
+facilitating this example. The Cp/Ct values for the different wave conditions are scaled
+values of the original Cp/Ct data for the IEA 15MW turbine.
+
+flow_field:
+  multidim_conditions:
+    Tp: 2.5
+    Hs: 3.01
+
+The solver will then use the nearest-neighbor interpolant. These conditions are currently global
+and used to select the interpolant at each turbine.
+
+Also note in the example below that there is a specific method for computing powers when
+using turbines with multi-dimensional Cp/Ct data under FlorisInterface, called
+'get_turbine_powers_multidim'. The normal 'get_turbine_powers' method will not work.
+"""
+
+# Initialize FLORIS with the given input file via FlorisInterface.
+fi = FlorisInterface("inputs/gch_multi_dim_cp_ct.yaml")
+
+# Convert to a simple two turbine layout
+fi.reinitialize(layout_x=[0., 500.], layout_y=[0., 0.])
+
+# Single wind speed and wind direction
+print('\n========================= Single Wind Direction and Wind Speed =========================')
+
+# Get the turbine powers assuming 1 wind speed and 1 wind direction
+fi.reinitialize(wind_directions=[270.], wind_speeds=[8.0])
+
+# Set the yaw angles to 0
+yaw_angles = np.zeros([1,1,2]) # 1 wind direction, 1 wind speed, 2 turbines
+fi.calculate_wake(yaw_angles=yaw_angles)
+
+# Get the turbine powers
+turbine_powers = fi.get_turbine_powers_multidim()/1000.
+print('The turbine power matrix should be of dimensions 1 WD X 1 WS X 2 Turbines')
+print(turbine_powers)
+print("Shape: ",turbine_powers.shape)
+
+# Single wind speed and multiple wind directions
+print('\n========================= Single Wind Direction and Multiple Wind Speeds ===============')
+
+
+wind_speeds = np.array([8.0, 9.0, 10.0])
+fi.reinitialize(wind_speeds=wind_speeds)
+yaw_angles = np.zeros([1,3,2]) # 1 wind direction, 3 wind speeds, 2 turbines
+fi.calculate_wake(yaw_angles=yaw_angles)
+turbine_powers = fi.get_turbine_powers_multidim()/1000.
+print('The turbine power matrix should be of dimensions 1 WD X 3 WS X 2 Turbines')
+print(turbine_powers)
+print("Shape: ",turbine_powers.shape)
+
+# Multiple wind speeds and multiple wind directions
+print('\n========================= Multiple Wind Directions and Multiple Wind Speeds ============')
+
+wind_directions = np.array([260., 270., 280.])
+wind_speeds = np.array([8.0, 9.0, 10.0])
+fi.reinitialize(wind_directions=wind_directions, wind_speeds=wind_speeds)
+yaw_angles = np.zeros([1,3,2]) # 1 wind direction, 3 wind speeds, 2 turbines
+fi.calculate_wake(yaw_angles=yaw_angles)
+turbine_powers = fi.get_turbine_powers_multidim()/1000.
+print('The turbine power matrix should be of dimensions 3 WD X 3 WS X 2 Turbines')
+print(turbine_powers)
+print("Shape: ",turbine_powers.shape)

examples/31_multi_dimensional_cp_ct_2Hs.py

@@ -0,0 +1,76 @@
+# Copyright 2023 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
+
+
+"""
+This example follows after example 30 but shows the effect of changing the Hs setting.
+
+NOTE: The multi-dimensional Cp/Ct data used in this example is fictional for the purposes of
+facilitating this example. The Cp/Ct values for the different wave conditions are scaled
+values of the original Cp/Ct data for the IEA 15MW turbine.
+"""
+
+# Initialize FLORIS with the given input file via FlorisInterface.
+fi = FlorisInterface("inputs/gch_multi_dim_cp_ct.yaml")
+
+# Make a second FLORIS interface with a different setting for Hs.
+# Note the multi-cp-ct file (iea_15MW_multi_dim_Tp_Hs.csv)
+# for the turbine model iea_15MW_floating_multi_dim_cp_ct.yaml
+# Defines Hs at 1 and 5.
+# The value in gch_multi_dim_cp_ct.yaml is 3.01 which will map
+# to 5 as the nearer value, so we set the other case to 1
+# for contrast.
+fi_dict_mod = fi.floris.as_dict()
+fi_dict_mod['flow_field']['multidim_conditions']['Hs'] = 1.0
+fi_hs_1 = FlorisInterface(fi_dict_mod)
+
+# Set both cases to 3 turbine layout
+fi.reinitialize(layout_x=[0., 500., 1000.], layout_y=[0., 0., 0.])
+fi_hs_1.reinitialize(layout_x=[0., 500., 1000.], layout_y=[0., 0., 0.])
+
+# Use a sweep of wind speeds
+wind_speeds = np.arange(5,20,1.)
+fi.reinitialize(wind_directions=[270.], wind_speeds=wind_speeds)
+fi_hs_1.reinitialize(wind_directions=[270.], wind_speeds=wind_speeds)
+
+# Calculate wakes with baseline yaw
+fi.calculate_wake()
+fi_hs_1.calculate_wake()
+
+# Collect the turbine powers in kW
+turbine_powers = fi.get_turbine_powers_multidim()/1000.
+turbine_powers_hs_1 = fi_hs_1.get_turbine_powers_multidim()/1000.
+
+# Plot the power in each case and the difference in power
+fig, axarr = plt.subplots(1,3,sharex=True,figsize=(12,4))
+
+for t_idx in range(3):
+    ax = axarr[t_idx]
+    ax.plot(wind_speeds, turbine_powers[0,:,t_idx], color='k', label='Hs=3.1 (5)')
+    ax.plot(wind_speeds, turbine_powers_hs_1[0,:,t_idx], color='r', label='Hs=1.0')
+    ax.grid(True)
+    ax.set_xlabel('Wind Speed (m/s)')
+    ax.set_title(f'Turbine {t_idx}')
+
+axarr[0].set_ylabel('Power (kW)')
+axarr[0].legend()
+fig.suptitle('Power of each turbine')
+
+plt.show()

examples/inputs/gch.yaml

@@ -134,6 +134,15 @@ flow_field:
   # The wind veer as a constant value for all points in the grid.
   wind_veer: 0.0
 
+  ###
+  # The conditions that are specified for use with the multi-dimensional Cp/Ct capbility.
+  # These conditions are external to FLORIS and specified by the user. They are used internally
+  # through a nearest-neighbor selection process to choose the correct Cp/Ct interpolants
+  # to use. These conditions are only used with the ``multidim_cp_ct`` velocity deficit model.
+  multidim_conditions:
+    Tp: 2.5
+    Hs: 3.01
+
 ###
 # Configure the wake model.
 wake:

examples/inputs/gch_multi_dim_cp_ct.yaml

@@ -0,0 +1,92 @@
+
+name: GCH multi dimensional Cp/Ct
+description: Three turbines using GCH model
+floris_version: v3.0.0
+
+logging:
+  console:
+    enable: true
+    level: WARNING
+  file:
+    enable: false
+    level: WARNING
+
+solver:
+  type: turbine_grid
+  turbine_grid_points: 3
+
+farm:
+  layout_x:
+  - 0.0
+  - 630.0
+  - 1260.0
+  layout_y:
+  - 0.0
+  - 0.0
+  - 0.0
+  turbine_type:
+  - iea_15MW_floating_multi_dim_cp_ct
+
+flow_field:
+  multidim_conditions:
+    Tp: 2.5
+    Hs: 3.01
+  air_density: 1.225
+  reference_wind_height: -1 # -1 is code for use the hub height
+  turbulence_intensity: 0.06
+  wind_directions:
+  - 270.0
+  wind_shear: 0.12
+  wind_speeds:
+  - 8.0
+  wind_veer: 0.0
+
+wake:
+  model_strings:
+    combination_model: sosfs
+    deflection_model: gauss
+    turbulence_model: crespo_hernandez
+    velocity_model: multidim_cp_ct
+
+  enable_secondary_steering: true
+  enable_yaw_added_recovery: true
+  enable_transverse_velocities: true
+
+  wake_deflection_parameters:
+    gauss:
+      ad: 0.0
+      alpha: 0.58
+      bd: 0.0
+      beta: 0.077
+      dm: 1.0
+      ka: 0.38
+      kb: 0.004
+    jimenez:
+      ad: 0.0
+      bd: 0.0
+      kd: 0.05
+
+  wake_velocity_parameters:
+    cc:
+      a_s: 0.179367259
+      b_s: 0.0118889215
+      c_s1: 0.0563691592
+      c_s2: 0.13290157
+      a_f: 3.11
+      b_f: -0.68
+      c_f: 2.41
+      alpha_mod: 1.0
+    multidim_cp_ct:
+      alpha: 0.58
+      beta: 0.077
+      ka: 0.38
+      kb: 0.004
+    jensen:
+      we: 0.05
+
+  wake_turbulence_parameters:
+    crespo_hernandez:
+      initial: 0.01
+      constant: 0.9
+      ai: 0.83
+      downstream: -0.25

floris/simulation/__init__.py

@@ -37,7 +37,21 @@
 import floris.logging_manager
 
 from .base import BaseClass, BaseModel, State
-from .turbine import average_velocity, axial_induction, Ct, power, rotor_effective_velocity, Turbine
+from .turbine import (
+    average_velocity,
+    axial_induction,
+    compute_tilt_angles_for_floating_turbines,
+    Ct,
+    power,
+    rotor_effective_velocity,
+    TiltTable,
+    Turbine
+)
+from .turbine_multi_dim import (
+    axial_induction_multidim,
+    Ct_multidim,
+    TurbineMultiDimensional
+)
 from .farm import Farm
 from .grid import (
     FlowFieldGrid,
@@ -57,6 +71,7 @@
     full_flow_sequential_solver,
     full_flow_turbopark_solver,
     sequential_solver,
+    sequential_multidim_solver,
     turbopark_solver,
 )
 from .floris import Floris