Feature: Add multidimensional Cp/Ct turbine definition support

examples/30_multi_dimensional_cp_ct.py

@@ -0,0 +1,101 @@
+# 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 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.:
+
+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.
+# turbine_powers = fi.get_turbine_powers()/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 wind direction
+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)
+
+# Single wind speed and wind direction
+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/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,22 @@
 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 (
+    _filter_convert,
+    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 +72,7 @@
     full_flow_sequential_solver,
     full_flow_turbopark_solver,
     sequential_solver,
+    sequential_multidim_solver,
     turbopark_solver,
 )
 from .floris import Floris

floris/simulation/farm.py

@@ -24,6 +24,7 @@
     BaseClass,
     State,
     Turbine,
+    TurbineMultiDimensional,
 )
 from floris.simulation.turbine import compute_tilt_angles_for_floating_turbines
 from floris.type_dec import (
@@ -247,13 +248,21 @@ def construct_turbine_correct_cp_ct_for_tilt(self):
         )
 
     def construct_turbine_map(self):
-        self.turbine_map = [Turbine.from_dict(turb) for turb in self.turbine_definitions]
+        if 'multi_dimensional_cp_ct' in self.turbine_definitions[0].keys():
+            self.turbine_map = [
+                TurbineMultiDimensional.from_dict(turb) for turb in self.turbine_definitions
+            ]
+        else:
+            self.turbine_map = [Turbine.from_dict(turb) for turb in self.turbine_definitions]
 
     def construct_turbine_fCts(self):
         self.turbine_fCts = {
             turb.turbine_type: turb.fCt_interp for turb in self.turbine_map
         }
 
+    def construct_multidim_turbine_fCts(self):
+        self.turbine_fCts = [turb.fCt_interp for turb in self.turbine_map]
+
     def construct_turbine_fTilts(self):
         self.turbine_fTilts = [(turb.turbine_type, turb.fTilt_interp) for turb in self.turbine_map]
 
@@ -262,6 +271,9 @@ def construct_turbine_power_interps(self):
             turb.turbine_type: turb.power_interp for turb in self.turbine_map
         }
 
+    def construct_multidim_turbine_power_interps(self):
+        self.turbine_power_interps = [turb.power_interp for turb in self.turbine_map]
+
     def construct_coordinates(self):
         self.coordinates = np.array([
             Vec3([x, y, z]) for x, y, z in zip(self.layout_x, self.layout_y, self.hub_heights)
@@ -279,6 +291,37 @@ def expand_farm_properties(
             sorted_coord_indices,
             axis=2
         )
+        if 'multi_dimensional_cp_ct' in self.turbine_definitions[0].keys():
+            wd_dim = np.shape(template_shape)[0]
+            ws_dim = np.shape(template_shape)[1]
+            if wd_dim != 1 | ws_dim != 0:
+                self.turbine_fCts_sorted = np.take_along_axis(
+                    np.reshape(
+                        np.repeat(self.turbine_fCts, wd_dim * ws_dim),
+                        np.shape(template_shape)
+                    ),
+                    sorted_coord_indices,
+                    axis=2
+                )
+                self.turbine_power_interps_sorted = np.take_along_axis(
+                    np.reshape(
+                        np.repeat(self.turbine_power_interps, wd_dim * ws_dim),
+                        np.shape(template_shape)
+                    ),
+                    sorted_coord_indices,
+                    axis=2
+                )
+            else:
+                self.turbine_fCts_sorted = np.take_along_axis(
+                    np.reshape(self.turbine_fCts, np.shape(template_shape)),
+                    sorted_coord_indices,
+                    axis=2
+                )
+                self.turbine_power_interps_sorted = np.take_along_axis(
+                    np.reshape(self.turbine_power_interps, np.shape(template_shape)),
+                    sorted_coord_indices,
+                    axis=2
+                )
         self.rotor_diameters_sorted = np.take_along_axis(
             self.rotor_diameters * template_shape,
             sorted_coord_indices,
@@ -355,6 +398,17 @@ def calculate_tilt_for_eff_velocities(self, rotor_effective_velocities):
         return tilt_angles
 
     def finalize(self, unsorted_indices):
+        if 'multi_dimensional_cp_ct' in self.turbine_definitions[0].keys():
+            self.turbine_fCts = np.take_along_axis(
+                self.turbine_fCts_sorted,
+                unsorted_indices[:,:,:,0,0],
+                axis=2
+            )
+            self.turbine_power_interps = np.take_along_axis(
+                self.turbine_power_interps_sorted,
+                unsorted_indices[:,:,:,0,0],
+                axis=2
+            )
         self.yaw_angles = np.take_along_axis(
             self.yaw_angles_sorted,
             unsorted_indices[:,:,:,0,0],

floris/simulation/floris.py

@@ -34,6 +34,7 @@
     full_flow_turbopark_solver,
     Grid,
     PointsGrid,
+    sequential_multidim_solver,
     sequential_solver,
     State,
     TurbineCubatureGrid,
@@ -82,8 +83,12 @@ def __attrs_post_init__(self) -> None:
 
         # Initialize farm quanitities that depend on other objects
         self.farm.construct_turbine_map()
-        self.farm.construct_turbine_fCts()
-        self.farm.construct_turbine_power_interps()
+        if self.wake.model_strings['velocity_model'] == 'multidim_cp_ct':
+            self.farm.construct_multidim_turbine_fCts()
+            self.farm.construct_multidim_turbine_power_interps()
+        else:
+            self.farm.construct_turbine_fCts()
+            self.farm.construct_turbine_power_interps()
         self.farm.construct_hub_heights()
         self.farm.construct_rotor_diameters()
         self.farm.construct_turbine_TSRs()
@@ -243,6 +248,13 @@ def steady_state_atmospheric_condition(self):
                 self.grid,
                 self.wake
             )
+        elif vel_model=="multidim_cp_ct":
+            sequential_multidim_solver(
+                self.farm,
+                self.flow_field,
+                self.grid,
+                self.wake
+            )
         else:
             sequential_solver(
                 self.farm,

floris/simulation/flow_field.py

@@ -41,8 +41,9 @@ class FlowField(BaseClass):
     air_density: float = field(converter=float)
     turbulence_intensity: float = field(converter=float)
     reference_wind_height: float = field(converter=float)
-    time_series : bool = field(default=False)
+    time_series: bool = field(default=False)
     heterogenous_inflow_config: dict = field(default=None)
+    multidim_conditions: dict = field(default=None)
 
     n_wind_speeds: int = field(init=False)
     n_wind_directions: int = field(init=False)