diff --git a/docs/examples.md b/docs/examples.md index a27517bea..e87627cd4 100644 --- a/docs/examples.md +++ b/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 diff --git a/examples/01_opening_floris_computing_power.py b/examples/01_opening_floris_computing_power.py index 049718f38..b006dfe4d 100644 --- a/examples/01_opening_floris_computing_power.py +++ b/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.]) diff --git a/examples/18_check_turbine.py b/examples/18_check_turbine.py index ec9d9b20e..e71f321ff 100644 --- a/examples/18_check_turbine.py +++ b/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)) diff --git a/examples/30_multi_dimensional_cp_ct.py b/examples/30_multi_dimensional_cp_ct.py new file mode 100644 index 000000000..2d2303018 --- /dev/null +++ b/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) diff --git a/examples/31_multi_dimensional_cp_ct_2Hs.py b/examples/31_multi_dimensional_cp_ct_2Hs.py new file mode 100644 index 000000000..6bbc31d6d --- /dev/null +++ b/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() diff --git a/examples/inputs/gch.yaml b/examples/inputs/gch.yaml index c30c35c3b..220fafeac 100644 --- a/examples/inputs/gch.yaml +++ b/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: diff --git a/examples/inputs/gch_multi_dim_cp_ct.yaml b/examples/inputs/gch_multi_dim_cp_ct.yaml new file mode 100644 index 000000000..8709fbcc7 --- /dev/null +++ b/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 diff --git a/floris/simulation/__init__.py b/floris/simulation/__init__.py index 12d30aab8..6da5c5ac5 100644 --- a/floris/simulation/__init__.py +++ b/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 diff --git a/floris/simulation/farm.py b/floris/simulation/farm.py index 211032ce1..dc1435a88 100644 --- a/floris/simulation/farm.py +++ b/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,22 @@ 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() \ + and self.turbine_definitions[0]['multi_dimensional_cp_ct'] is True: + 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 +272,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 +292,38 @@ def expand_farm_properties( sorted_coord_indices, axis=2 ) + if 'multi_dimensional_cp_ct' in self.turbine_definitions[0].keys() \ + and self.turbine_definitions[0]['multi_dimensional_cp_ct'] is True: + 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 +400,18 @@ 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() \ + and self.turbine_definitions[0]['multi_dimensional_cp_ct'] is True: + 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], diff --git a/floris/simulation/floris.py b/floris/simulation/floris.py index a5f96c2a4..09722a2d7 100644 --- a/floris/simulation/floris.py +++ b/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, diff --git a/floris/simulation/flow_field.py b/floris/simulation/flow_field.py index a8dbf7393..2781d8c12 100644 --- a/floris/simulation/flow_field.py +++ b/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) diff --git a/floris/simulation/solver.py b/floris/simulation/solver.py index 713f8aeb1..6e53a718a 100644 --- a/floris/simulation/solver.py +++ b/floris/simulation/solver.py @@ -27,6 +27,11 @@ TurbineGrid, ) from floris.simulation.turbine import average_velocity +from floris.simulation.turbine_multi_dim import ( + axial_induction_multidim, + Ct_multidim, + multidim_Ct_down_select, +) from floris.simulation.wake import WakeModelManager from floris.simulation.wake_deflection.empirical_gauss import yaw_added_wake_mixing from floris.simulation.wake_deflection.gauss import ( @@ -1504,3 +1509,208 @@ def full_flow_empirical_gauss_solver( flow_field.u_sorted = flow_field.u_initial_sorted - wake_field flow_field.v_sorted += v_wake flow_field.w_sorted += w_wake + + +def sequential_multidim_solver( + farm: Farm, + flow_field: FlowField, + grid: TurbineGrid, + model_manager: WakeModelManager +) -> None: + # Algorithm + # For each turbine, calculate its effect on every downstream turbine. + # For the current turbine, we are calculating the deficit that it adds to downstream turbines. + # Integrate this into the main data structure. + # Move on to the next turbine. + + # <> + deflection_model_args = model_manager.deflection_model.prepare_function(grid, flow_field) + deficit_model_args = model_manager.velocity_model.prepare_function(grid, flow_field) + downselect_turbine_fCts = multidim_Ct_down_select( + farm.turbine_fCts_sorted, + flow_field.multidim_conditions, + ) + + # This is u_wake + wake_field = np.zeros_like(flow_field.u_initial_sorted) + v_wake = np.zeros_like(flow_field.v_initial_sorted) + w_wake = np.zeros_like(flow_field.w_initial_sorted) + + turbine_turbulence_intensity = ( + flow_field.turbulence_intensity + * np.ones((flow_field.n_wind_directions, flow_field.n_wind_speeds, farm.n_turbines, 1, 1)) + ) + ambient_turbulence_intensity = flow_field.turbulence_intensity + + # Calculate the velocity deficit sequentially from upstream to downstream turbines + for i in range(grid.n_turbines): + + # Get the current turbine quantities + x_i = np.mean(grid.x_sorted[:, :, i:i+1], axis=(3, 4)) + x_i = x_i[:, :, :, None, None] + y_i = np.mean(grid.y_sorted[:, :, i:i+1], axis=(3, 4)) + y_i = y_i[:, :, :, None, None] + z_i = np.mean(grid.z_sorted[:, :, i:i+1], axis=(3, 4)) + z_i = z_i[:, :, :, None, None] + + u_i = flow_field.u_sorted[:, :, i:i+1] + v_i = flow_field.v_sorted[:, :, i:i+1] + + ct_i = Ct_multidim( + velocities=flow_field.u_sorted, + yaw_angle=farm.yaw_angles_sorted, + tilt_angle=farm.tilt_angles_sorted, + ref_tilt_cp_ct=farm.ref_tilt_cp_cts_sorted, + fCt=downselect_turbine_fCts, + tilt_interp=farm.turbine_fTilts, + correct_cp_ct_for_tilt=farm.correct_cp_ct_for_tilt_sorted, + turbine_type_map=farm.turbine_type_map_sorted, + ix_filter=[i], + average_method=grid.average_method, + cubature_weights=grid.cubature_weights + ) + # Since we are filtering for the i'th turbine in the Ct function, + # get the first index here (0:1) + ct_i = ct_i[:, :, 0:1, None, None] + axial_induction_i = axial_induction_multidim( + velocities=flow_field.u_sorted, + yaw_angle=farm.yaw_angles_sorted, + tilt_angle=farm.tilt_angles_sorted, + ref_tilt_cp_ct=farm.ref_tilt_cp_cts_sorted, + fCt=downselect_turbine_fCts, + tilt_interp=farm.turbine_fTilts, + correct_cp_ct_for_tilt=farm.correct_cp_ct_for_tilt_sorted, + turbine_type_map=farm.turbine_type_map_sorted, + ix_filter=[i], + average_method=grid.average_method, + cubature_weights=grid.cubature_weights + ) + # Since we are filtering for the i'th turbine in the axial induction function, + # get the first index here (0:1) + axial_induction_i = axial_induction_i[:, :, 0:1, None, None] + turbulence_intensity_i = turbine_turbulence_intensity[:, :, i:i+1] + yaw_angle_i = farm.yaw_angles_sorted[:, :, i:i+1, None, None] + hub_height_i = farm.hub_heights_sorted[:, :, i:i+1, None, None] + rotor_diameter_i = farm.rotor_diameters_sorted[:, :, i:i+1, None, None] + TSR_i = farm.TSRs_sorted[:, :, i:i+1, None, None] + + effective_yaw_i = np.zeros_like(yaw_angle_i) + effective_yaw_i += yaw_angle_i + + if model_manager.enable_secondary_steering: + added_yaw = wake_added_yaw( + u_i, + v_i, + flow_field.u_initial_sorted, + grid.y_sorted[:, :, i:i+1] - y_i, + grid.z_sorted[:, :, i:i+1], + rotor_diameter_i, + hub_height_i, + ct_i, + TSR_i, + axial_induction_i, + flow_field.wind_shear, + ) + effective_yaw_i += added_yaw + + # Model calculations + # NOTE: exponential + deflection_field = model_manager.deflection_model.function( + x_i, + y_i, + effective_yaw_i, + turbulence_intensity_i, + ct_i, + rotor_diameter_i, + **deflection_model_args, + ) + + if model_manager.enable_transverse_velocities: + v_wake, w_wake = calculate_transverse_velocity( + u_i, + flow_field.u_initial_sorted, + flow_field.dudz_initial_sorted, + grid.x_sorted - x_i, + grid.y_sorted - y_i, + grid.z_sorted, + rotor_diameter_i, + hub_height_i, + yaw_angle_i, + ct_i, + TSR_i, + axial_induction_i, + flow_field.wind_shear, + ) + + if model_manager.enable_yaw_added_recovery: + I_mixing = yaw_added_turbulence_mixing( + u_i, + turbulence_intensity_i, + v_i, + flow_field.w_sorted[:, :, i:i+1], + v_wake[:, :, i:i+1], + w_wake[:, :, i:i+1], + ) + gch_gain = 2 + turbine_turbulence_intensity[:, :, i:i+1] = turbulence_intensity_i + gch_gain * I_mixing + + # NOTE: exponential + velocity_deficit = model_manager.velocity_model.function( + x_i, + y_i, + z_i, + axial_induction_i, + deflection_field, + yaw_angle_i, + turbulence_intensity_i, + ct_i, + hub_height_i, + rotor_diameter_i, + **deficit_model_args, + ) + + wake_field = model_manager.combination_model.function( + wake_field, + velocity_deficit * flow_field.u_initial_sorted + ) + + wake_added_turbulence_intensity = model_manager.turbulence_model.function( + ambient_turbulence_intensity, + grid.x_sorted, + x_i, + rotor_diameter_i, + axial_induction_i, + ) + + # Calculate wake overlap for wake-added turbulence (WAT) + area_overlap = ( + np.sum(velocity_deficit * flow_field.u_initial_sorted > 0.05, axis=(3, 4)) + / (grid.grid_resolution * grid.grid_resolution) + ) + area_overlap = area_overlap[:, :, :, None, None] + + # Modify wake added turbulence by wake area overlap + downstream_influence_length = 15 * rotor_diameter_i + ti_added = ( + area_overlap + * np.nan_to_num(wake_added_turbulence_intensity, posinf=0.0) + * (grid.x_sorted > x_i) + * (np.abs(y_i - grid.y_sorted) < 2 * rotor_diameter_i) + * (grid.x_sorted <= downstream_influence_length + x_i) + ) + + # Combine turbine TIs with WAT + turbine_turbulence_intensity = np.maximum( + np.sqrt( ti_added ** 2 + ambient_turbulence_intensity ** 2 ), + turbine_turbulence_intensity + ) + + flow_field.u_sorted = flow_field.u_initial_sorted - wake_field + flow_field.v_sorted += v_wake + flow_field.w_sorted += w_wake + + flow_field.turbulence_intensity_field_sorted = turbine_turbulence_intensity + flow_field.turbulence_intensity_field_sorted_avg = np.mean( + turbine_turbulence_intensity, + axis=(3,4) + )[:, :, :, None, None] diff --git a/floris/simulation/turbine.py b/floris/simulation/turbine.py index b284b10cd..0c056322d 100644 --- a/floris/simulation/turbine.py +++ b/floris/simulation/turbine.py @@ -521,9 +521,9 @@ class PowerThrustTable(FromDictMixin): ValueError: Raised if the power, thrust, and wind_speed are not all 1-d array-like shapes. ValueError: Raised if power, thrust, and wind_speed don't have the same number of values. """ - power: NDArrayFloat = field(converter=floris_array_converter) - thrust: NDArrayFloat = field(converter=floris_array_converter) - wind_speed: NDArrayFloat = field(converter=floris_array_converter) + power: NDArrayFloat = field(default=[], converter=floris_array_converter) + thrust: NDArrayFloat = field(default=[], converter=floris_array_converter) + wind_speed: NDArrayFloat = field(default=[], converter=floris_array_converter) def __attrs_post_init__(self) -> None: # Validate the power, thrust, and wind speed inputs. @@ -624,9 +624,11 @@ class Turbine(BaseClass): generator_efficiency: float = field() ref_density_cp_ct: float = field() ref_tilt_cp_ct: float = field() - power_thrust_table: PowerThrustTable = field(converter=PowerThrustTable.from_dict) - floating_tilt_table = field(default=None) - floating_correct_cp_ct_for_tilt = field(default=None) + power_thrust_table: PowerThrustTable = field(default=None) + floating_tilt_table: TiltTable = field(default=None) + floating_correct_cp_ct_for_tilt: bool = field(default=None) + power_thrust_data_file: str = field(default=None) + multi_dimensional_cp_ct: bool = field(default=False) # rloc: float = float_attrib() # TODO: goes here or on the Grid? # use_points_on_perimeter: bool = bool_attrib() @@ -650,6 +652,7 @@ class Turbine(BaseClass): def __attrs_post_init__(self) -> None: # Post-init initialization for the power curve interpolation functions + self.power_thrust_table = PowerThrustTable.from_dict(self.power_thrust_table) wind_speeds = self.power_thrust_table.wind_speed self.fCp_interp = interp1d( wind_speeds, diff --git a/floris/simulation/turbine_multi_dim.py b/floris/simulation/turbine_multi_dim.py new file mode 100644 index 000000000..3e2cc7b8d --- /dev/null +++ b/floris/simulation/turbine_multi_dim.py @@ -0,0 +1,518 @@ +# 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 + +from __future__ import annotations + +import copy +from collections.abc import Iterable + +import numpy as np +import pandas as pd +from attrs import define, field +from flatten_dict import flatten +from scipy.interpolate import interp1d + +# import floris.simulation.turbine as turbine +from floris.simulation import ( + average_velocity, + compute_tilt_angles_for_floating_turbines, + TiltTable, + Turbine, +) +from floris.simulation.turbine import _filter_convert +from floris.type_dec import ( + NDArrayBool, + NDArrayFilter, + NDArrayFloat, + NDArrayInt, + NDArrayObject, +) +from floris.utilities import cosd + + +def power_multidim( + ref_density_cp_ct: float, + rotor_effective_velocities: NDArrayFloat, + power_interp: NDArrayObject, + ix_filter: NDArrayInt | Iterable[int] | None = None, +) -> NDArrayFloat: + """Power produced by a turbine defined with multi-dimensional + Cp/Ct values, adjusted for yaw and tilt. Value given in Watts. + + Args: + ref_density_cp_cts (NDArrayFloat[wd, ws, turbines]): The reference density for each turbine + rotor_effective_velocities (NDArrayFloat[wd, ws, turbines, grid1, grid2]): The rotor + effective velocities at a turbine. + power_interp (NDArrayObject[wd, ws, turbines]): The power interpolation function + for each turbine. + ix_filter (NDArrayInt, optional): The boolean array, or + integer indices to filter out before calculation. Defaults to None. + + Returns: + NDArrayFloat: The power, in Watts, for each turbine after adjusting for yaw and tilt. + """ + # TODO: Change the order of input arguments to be consistent with the other + # utility functions - velocities first... + # Update to power calculation which replaces the fixed pP exponent with + # an exponent pW, that changes the effective wind speed input to the power + # calculation, rather than scaling the power. This better handles power + # loss to yaw in above rated conditions + # + # based on the paper "Optimising yaw control at wind farm level" by + # Ervin Bossanyi + + # TODO: check this - where is it? + # P = 1/2 rho A V^3 Cp + + # Down-select inputs if ix_filter is given + if ix_filter is not None: + ix_filter = _filter_convert(ix_filter, rotor_effective_velocities) + power_interp = power_interp[:, :, ix_filter] + rotor_effective_velocities = rotor_effective_velocities[:, :, ix_filter] + # Loop over each turbine to get power for all turbines + p = np.zeros(np.shape(rotor_effective_velocities)) + for i, wd in enumerate(power_interp): + for j, ws in enumerate(wd): + for k, turb in enumerate(ws): + p[i, j, k] = power_interp[i, j, k](rotor_effective_velocities[i, j, k]) + + return p * ref_density_cp_ct + + +def Ct_multidim( + velocities: NDArrayFloat, + yaw_angle: NDArrayFloat, + tilt_angle: NDArrayFloat, + ref_tilt_cp_ct: NDArrayFloat, + fCt: list, + tilt_interp: NDArrayObject, + correct_cp_ct_for_tilt: NDArrayBool, + turbine_type_map: NDArrayObject, + ix_filter: NDArrayFilter | Iterable[int] | None = None, + average_method: str = "cubic-mean", + cubature_weights: NDArrayFloat | None = None +) -> NDArrayFloat: + + """Thrust coefficient of a turbine defined with multi-dimensional + Cp/Ct values, incorporating the yaw angle. The value is interpolated + from the coefficient of thrust vs wind speed table using the rotor + swept area average velocity. + + Args: + velocities (NDArrayFloat[wd, ws, turbines, grid1, grid2]): The velocity field at + a turbine. + yaw_angle (NDArrayFloat[wd, ws, turbines]): The yaw angle for each turbine. + tilt_angle (NDArrayFloat[wd, ws, turbines]): The tilt angle for each turbine. + ref_tilt_cp_ct (NDArrayFloat[wd, ws, turbines]): The reference tilt angle for each turbine + that the Cp/Ct tables are defined at. + fCt (list): The thrust coefficient interpolation functions for each turbine. + tilt_interp (Iterable[tuple]): The tilt interpolation functions for each + turbine. + correct_cp_ct_for_tilt (NDArrayBool[wd, ws, turbines]): Boolean for determining if the + turbines Cp and Ct should be corrected for tilt. + turbine_type_map: (NDArrayObject[wd, ws, turbines]): The Turbine type definition + for each turbine. + ix_filter (NDArrayFilter | Iterable[int] | None, optional): The boolean array, or + integer indices as an iterable of array to filter out before calculation. + Defaults to None. + + Returns: + NDArrayFloat: Coefficient of thrust for each requested turbine. + """ + + if isinstance(yaw_angle, list): + yaw_angle = np.array(yaw_angle) + + if isinstance(tilt_angle, list): + tilt_angle = np.array(tilt_angle) + + # Down-select inputs if ix_filter is given + if ix_filter is not None: + ix_filter = _filter_convert(ix_filter, yaw_angle) + velocities = velocities[:, :, ix_filter] + yaw_angle = yaw_angle[:, :, ix_filter] + tilt_angle = tilt_angle[:, :, ix_filter] + ref_tilt_cp_ct = ref_tilt_cp_ct[:, :, ix_filter] + fCt = fCt[:, :, ix_filter] + turbine_type_map = turbine_type_map[:, :, ix_filter] + correct_cp_ct_for_tilt = correct_cp_ct_for_tilt[:, :, ix_filter] + + average_velocities = average_velocity( + velocities, + method=average_method, + cubature_weights=cubature_weights + ) + + # Compute the tilt, if using floating turbines + old_tilt_angle = copy.deepcopy(tilt_angle) + tilt_angle = compute_tilt_angles_for_floating_turbines( + turbine_type_map, + tilt_angle, + tilt_interp, + average_velocities, + ) + # Only update tilt angle if requested (if the tilt isn't accounted for in the Ct curve) + tilt_angle = np.where(correct_cp_ct_for_tilt, tilt_angle, old_tilt_angle) + + # Loop over each turbine to get thrust coefficient for all turbines + thrust_coefficient = np.zeros(np.shape(average_velocities)) + for i, wd in enumerate(fCt): + for j, ws in enumerate(wd): + for k, turb in enumerate(ws): + thrust_coefficient[i, j, k] = fCt[i, j, k](average_velocities[i, j, k]) + thrust_coefficient = np.clip(thrust_coefficient, 0.0001, 0.9999) + effective_thrust = thrust_coefficient * cosd(yaw_angle) * cosd(tilt_angle - ref_tilt_cp_ct) + return effective_thrust + + +def axial_induction_multidim( + velocities: NDArrayFloat, # (wind directions, wind speeds, turbines, grid, grid) + yaw_angle: NDArrayFloat, # (wind directions, wind speeds, turbines) + tilt_angle: NDArrayFloat, # (wind directions, wind speeds, turbines) + ref_tilt_cp_ct: NDArrayFloat, + fCt: list, # (turbines) + tilt_interp: NDArrayObject, # (turbines) + correct_cp_ct_for_tilt: NDArrayBool, # (wind directions, wind speeds, turbines) + turbine_type_map: NDArrayObject, # (wind directions, 1, turbines) + ix_filter: NDArrayFilter | Iterable[int] | None = None, + average_method: str = "cubic-mean", + cubature_weights: NDArrayFloat | None = None +) -> NDArrayFloat: + """Axial induction factor of the turbines defined with multi-dimensional + Cp/Ct values, incorporating the thrust coefficient and yaw angle. + + Args: + velocities (NDArrayFloat): The velocity field at each turbine; should be shape: + (number of turbines, ngrid, ngrid), or (ngrid, ngrid) for a single turbine. + yaw_angle (NDArrayFloat[wd, ws, turbines]): The yaw angle for each turbine. + tilt_angle (NDArrayFloat[wd, ws, turbines]): The tilt angle for each turbine. + ref_tilt_cp_ct (NDArrayFloat[wd, ws, turbines]): The reference tilt angle for each turbine + that the Cp/Ct tables are defined at. + fCt (list): The thrust coefficient interpolation functions for each turbine. + tilt_interp (Iterable[tuple]): The tilt interpolation functions for each + turbine. + correct_cp_ct_for_tilt (NDArrayBool[wd, ws, turbines]): Boolean for determining if the + turbines Cp and Ct should be corrected for tilt. + turbine_type_map: (NDArrayObject[wd, ws, turbines]): The Turbine type definition + for each turbine. + ix_filter (NDArrayFilter | Iterable[int] | None, optional): The boolean array, or + integer indices (as an aray or iterable) to filter out before calculation. + Defaults to None. + + Returns: + Union[float, NDArrayFloat]: [description] + """ + + if isinstance(yaw_angle, list): + yaw_angle = np.array(yaw_angle) + + # TODO: Should the tilt_angle used for the return calculation be modified the same as the + # tilt_angle in Ct, if the user has supplied a tilt/wind_speed table? + if isinstance(tilt_angle, list): + tilt_angle = np.array(tilt_angle) + + # Get Ct first before modifying any data + thrust_coefficient = Ct_multidim( + velocities, + yaw_angle, + tilt_angle, + ref_tilt_cp_ct, + fCt, + tilt_interp, + correct_cp_ct_for_tilt, + turbine_type_map, + ix_filter, + average_method, + cubature_weights + ) + + # Then, process the input arguments as needed for this function + ix_filter = _filter_convert(ix_filter, yaw_angle) + if ix_filter is not None: + yaw_angle = yaw_angle[:, :, ix_filter] + tilt_angle = tilt_angle[:, :, ix_filter] + ref_tilt_cp_ct = ref_tilt_cp_ct[:, :, ix_filter] + + return ( + 0.5 + / (cosd(yaw_angle) + * cosd(tilt_angle - ref_tilt_cp_ct)) + * ( + 1 - np.sqrt( + 1 - thrust_coefficient * cosd(yaw_angle) * cosd(tilt_angle - ref_tilt_cp_ct) + ) + ) + ) + + +def multidim_Ct_down_select( + turbine_fCts, + conditions, +) -> list: + """ + Ct interpolants are down selected from the multi-dimensional Ct data + provided for the turbine based on the specified conditions. + + Args: + turbine_fCts (NDArray[wd, ws, turbines]): The Ct interpolants generated from the + multi-dimensional Ct turbine data for all specified conditions. + conditions (dict): The conditions at which to determine which Ct interpolant to use. + + Returns: + NDArray: The downselected Ct interpolants for the selected conditions. + """ + downselect_turbine_fCts = np.empty_like(turbine_fCts) + # Loop over the wind directions, wind speeds, and turbines, finding the Ct interpolant + # that is closest to the specified multi-dimensional condition. + for i, wd in enumerate(turbine_fCts): + for j, ws in enumerate(wd): + for k, turb in enumerate(ws): + # Get the interpolant keys in float type for comparison + keys_float = np.array([[float(v) for v in val] for val in turb.keys()]) + + # Find the nearest key to the specified conditions. + key_vals = [] + for ii, cond in enumerate(conditions.values()): + key_vals.append( + keys_float[:, ii][np.absolute(keys_float[:, ii] - cond).argmin()] + ) + + downselect_turbine_fCts[i, j, k] = turb[tuple(key_vals)] + + return downselect_turbine_fCts + + +def multidim_power_down_select( + power_interps, + conditions, +) -> list: + """ + Cp interpolants are down selected from the multi-dimensional Cp data + provided for the turbine based on the specified conditions. + + Args: + power_interps (NDArray[wd, ws, turbines]): The power interpolants generated from the + multi-dimensional Cp turbine data for all specified conditions. + conditions (dict): The conditions at which to determine which Ct interpolant to use. + + Returns: + NDArray: The downselected power interpolants for the selected conditions. + """ + downselect_power_interps = np.empty_like(power_interps) + # Loop over the wind directions, wind speeds, and turbines, finding the power interpolant + # that is closest to the specified multi-dimensional condition. + for i, wd in enumerate(power_interps): + for j, ws in enumerate(wd): + for k, turb in enumerate(ws): + # Get the interpolant keys in float type for comparison + keys_float = np.array([[float(v) for v in val] for val in turb.keys()]) + + # Find the nearest key to the specified conditions. + key_vals = [] + for ii, cond in enumerate(conditions.values()): + key_vals.append( + keys_float[:, ii][np.absolute(keys_float[:, ii] - cond).argmin()] + ) + + # Use the constructed key to choose the correct interpolant + downselect_power_interps[i, j, k] = turb[tuple(key_vals)] + + return downselect_power_interps + + +@define +class MultiDimensionalPowerThrustTable(): + """Helper class to convert the multi-dimensional inputs to a dictionary of objects. + """ + + @classmethod + def from_dataframe(self, df) -> None: + # Validate the dataframe + if not all(ele in df.columns.values.tolist() for ele in ["ws", "Cp", "Ct"]): + print(df.columns.values.tolist()) + raise ValueError("Multidimensional data missing required ws/Cp/Ct data.") + if df.columns.values[-3:].tolist() != ["ws", "Cp", "Ct"]: + print(df.columns.values[-3:].tolist()) + raise ValueError( + "Multidimensional data not in correct form. ws, Cp, and Ct must be " + "defined as the last 3 columns, in that order." + ) + + # Extract the supplied dimensions, minus the required ws, Cp, and Ct columns. + keys = df.columns.values[:-3].tolist() + values = [df[df.columns.values[i]].unique().tolist() for i in range(len(keys))] + values = [[str(val) for val in value] for value in values] + + # Functions for recursively building a nested dictionary from + # an arbitrary number of paired-inputs. + def add_level(obj, k, v): + tmp = {} + for val in v: + tmp.update({val: []}) + obj.update({k: tmp}) + return obj + + def add_sub_level(obj, k): + tmp = {} + for key in k: + tmp.update({key: obj}) + return tmp + + obj = {} + # Reverse the lists to start from the lowest level of the dictionary + keys.reverse() + values.reverse() + # Recursively build a nested dictionary from the user-supplied dimensions + for i, key in enumerate(keys): + if i == 0: + obj = add_level(obj, key, values[i]) + else: + obj = add_sub_level(obj, values[i]) + obj = {key: obj} + + return flatten(obj) + + +@define +class TurbineMultiDimensional(Turbine): + """ + Turbine is a class containing objects pertaining to the individual + turbines. + + Turbine is a model class representing a particular wind turbine. It + is largely a container of data and parameters, but also contains + methods to probe properties for output. + + Parameters: + rotor_diameter (:py:obj: float): The rotor diameter (m). + hub_height (:py:obj: float): The hub height (m). + pP (:py:obj: float): The cosine exponent relating the yaw + misalignment angle to power. + pT (:py:obj: float): The cosine exponent relating the rotor + tilt angle to power. + generator_efficiency (:py:obj: float): The generator + efficiency factor used to scale the power production. + ref_density_cp_ct (:py:obj: float): The density at which the provided + cp and ct is defined + power_thrust_table (PowerThrustTable): A dictionary containing the + following key-value pairs: + + power (:py:obj: List[float]): The coefficient of power at + different wind speeds. + thrust (:py:obj: List[float]): The coefficient of thrust + at different wind speeds. + wind_speed (:py:obj: List[float]): The wind speeds for + which the power and thrust values are provided (m/s). + ngrid (*int*, optional): The square root of the number + of points to use on the turbine grid. This number will be + squared so that the points can be evenly distributed. + Defaults to 5. + rloc (:py:obj: float, optional): A value, from 0 to 1, that determines + the width/height of the grid of points on the rotor as a ratio of + the rotor radius. + Defaults to 0.5. + """ + + power_thrust_data_file: str = field(default=None) + multi_dimensional_cp_ct: bool = field(default=False) + + # rloc: float = float_attrib() # TODO: goes here or on the Grid? + # use_points_on_perimeter: bool = bool_attrib() + + # Initialized in the post_init function + # rotor_radius: float = field(init=False) + # rotor_area: float = field(init=False) + # fCp_interp: interp1d = field(init=False) + # fCt_interp: interp1d = field(init=False) + # power_interp: interp1d = field(init=False) + # tilt_interp: interp1d = field(init=False) + + + # For the following parameters, use default values if not user-specified + # self.rloc = float(input_dictionary["rloc"]) if "rloc" in input_dictionary else 0.5 + # if "use_points_on_perimeter" in input_dictionary: + # self.use_points_on_perimeter = bool(input_dictionary["use_points_on_perimeter"]) + # else: + # self.use_points_on_perimeter = False + + def __attrs_post_init__(self) -> None: + + # Read in the multi-dimensional data supplied by the user. + df = pd.read_csv(self.power_thrust_data_file) + + # Build the multi-dimensional power/thrust table + self.power_thrust_data = MultiDimensionalPowerThrustTable.from_dataframe(df) + + # Create placeholders for the interpolation functions + self.fCt_interp = {} + self.power_interp = {} + + # Down-select the DataFrame to have just the ws, Cp, and Ct values + index_col = df.columns.values[:-3] + df2 = df.set_index(index_col.tolist()) + + # Loop over the multi-dimensional keys to get the correct ws/Cp/Ct data to make + # the Ct and power interpolants. + for key in df2.index.unique(): + # Select the correct ws/Cp/Ct data + data = df2.loc[key] + + # Build the interpolants + wind_speeds = data['ws'].values + self.fCp_interp = interp1d( + wind_speeds, + data['Cp'].values, + fill_value=(0.0, 1.0), + bounds_error=False, + ) + inner_power = ( + 0.5 * self.rotor_area + * self.fCp_interp(wind_speeds) + * self.generator_efficiency + * wind_speeds ** 3 + ) + self.power_interp.update({ + key: interp1d( + wind_speeds, + inner_power, + bounds_error=False, + fill_value=0 + ) + }) + self.fCt_interp.update({ + key: interp1d( + wind_speeds, + data['Ct'].values, + fill_value=(0.0001, 0.9999), + bounds_error=False, + ) + }) + + # If defined, create a tilt interpolation function for floating turbines. + # fill_value currently set to apply the min or max tilt angles if outside + # of the interpolation range. + if self.floating_tilt_table is not None: + self.floating_tilt_table = TiltTable.from_dict(self.floating_tilt_table) + self.fTilt_interp = interp1d( + self.floating_tilt_table.wind_speeds, + self.floating_tilt_table.tilt, + fill_value=(0.0, self.floating_tilt_table.tilt[-1]), + bounds_error=False, + ) + self.tilt_interp = self.fTilt_interp + self.correct_cp_ct_for_tilt = self.floating_correct_cp_ct_for_tilt + else: + self.fTilt_interp = None + self.tilt_interp = None + self.correct_cp_ct_for_tilt = False diff --git a/floris/simulation/wake.py b/floris/simulation/wake.py index a141c94be..558f6ecbe 100644 --- a/floris/simulation/wake.py +++ b/floris/simulation/wake.py @@ -65,7 +65,8 @@ "gauss": GaussVelocityDeficit, "jensen": JensenVelocityDeficit, "turbopark": TurbOParkVelocityDeficit, - "empirical_gauss": EmpiricalGaussVelocityDeficit + "empirical_gauss": EmpiricalGaussVelocityDeficit, + "multidim_cp_ct": GaussVelocityDeficit }, } diff --git a/floris/tools/floris_interface.py b/floris/tools/floris_interface.py index 76b2dcfee..c5404d0b2 100644 --- a/floris/tools/floris_interface.py +++ b/floris/tools/floris_interface.py @@ -29,6 +29,7 @@ power, rotor_effective_velocity, ) +from floris.simulation.turbine_multi_dim import multidim_power_down_select, power_multidim from floris.tools.cut_plane import CutPlane from floris.type_dec import NDArrayFloat @@ -586,7 +587,7 @@ def get_turbine_powers(self) -> NDArrayFloat: """Calculates the power at each turbine in the windfarm. Returns: - NDArrayFloat: [description] + NDArrayFloat: Powers at each turbine. """ # Confirm calculate wake has been run @@ -608,6 +609,37 @@ def get_turbine_powers(self) -> NDArrayFloat: ) return turbine_powers + def get_turbine_powers_multidim(self) -> NDArrayFloat: + """Calculates the power at each turbine in the windfarm + when using multi-dimensional Cp/Ct turbine definitions. + + Returns: + NDArrayFloat: Powers at each turbine. + """ + + # Confirm calculate wake has been run + if self.floris.state is not State.USED: + raise RuntimeError( + "Can't run function `FlorisInterface.get_turbine_powers_multidim` without " + "first running `FlorisInterface.calculate_wake`." + ) + # Check for negative velocities, which could indicate bad model + # parameters or turbines very closely spaced. + if (self.turbine_effective_velocities < 0.).any(): + self.logger.warning("Some rotor effective velocities are negative.") + + turbine_power_interps = multidim_power_down_select( + self.floris.farm.turbine_power_interps, + self.floris.flow_field.multidim_conditions + ) + + turbine_powers = power_multidim( + ref_density_cp_ct=self.floris.farm.ref_density_cp_cts, + rotor_effective_velocities=self.turbine_effective_velocities, + power_interp=turbine_power_interps, + ) + return turbine_powers + def get_turbine_Cts(self) -> NDArrayFloat: turbine_Cts = Ct( velocities=self.floris.flow_field.u, diff --git a/floris/turbine_library/iea_15MW_floating_multi_dim_cp_ct.yaml b/floris/turbine_library/iea_15MW_floating_multi_dim_cp_ct.yaml new file mode 100644 index 000000000..ea8623eee --- /dev/null +++ b/floris/turbine_library/iea_15MW_floating_multi_dim_cp_ct.yaml @@ -0,0 +1,29 @@ +turbine_type: 'iea_15MW_floating' +generator_efficiency: 1.0 +hub_height: 150.0 +pP: 1.88 +pT: 1.88 +rotor_diameter: 242.24 +TSR: 8.0 +ref_density_cp_ct: 1.225 +ref_tilt_cp_ct: 6.0 +multi_dimensional_cp_ct: True +power_thrust_data_file: '../floris/turbine_library/iea_15MW_multi_dim_Tp_Hs.csv' +floating_tilt_table: + tilt: + - 5.747296314800103 + - 7.2342400188651068 + - 9.0468701999352397 + - 9.762182013267733 + - 8.795649572299896 + - 8.089078308325314 + - 7.7229584934943614 + wind_speeds: + - 4.0 + - 6.0 + - 8.0 + - 10.0 + - 12.0 + - 14.0 + - 16.0 +floating_correct_cp_ct_for_tilt: True diff --git a/floris/turbine_library/iea_15MW_multi_dim_Tp_Hs.csv b/floris/turbine_library/iea_15MW_multi_dim_Tp_Hs.csv new file mode 100644 index 000000000..b30eac5a3 --- /dev/null +++ b/floris/turbine_library/iea_15MW_multi_dim_Tp_Hs.csv @@ -0,0 +1,213 @@ +Tp,Hs,ws,Cp,Ct +2,1,0,0,0 +2,1,3,0.049361236,0.817533319 +2,1,3.54953237,0.224324252,0.792115292 +2,1,4.067900771,0.312216418,0.786401899 +2,1,4.553906848,0.36009987,0.788898744 +2,1,5.006427063,0.38761204,0.790774576 +2,1,5.424415288,0.404010164,0.79208669 +2,1,5.806905228,0.413979324,0.79185809 +2,1,6.153012649,0.420083692,0.7903853 +2,1,6.461937428,0.423787764,0.788253035 +2,1,6.732965398,0.425977895,0.785845184 +2,1,6.965470002,0.427193272,0.783367164 +2,1,7.158913742,0.427183505,0.77853469 +2,1,7.312849418,0.426860928,0.77853469 +2,1,7.426921164,0.426617959,0.77853469 +2,1,7.500865272,0.426458783,0.77853469 +2,1,7.534510799,0.426385957,0.77853469 +2,1,7.541241633,0.426371389,0.77853469 +2,1,7.58833327,0.426268826,0.77853469 +2,1,7.675676842,0.426077456,0.77853469 +2,1,7.803070431,0.425795302,0.77853469 +2,1,7.970219531,0.425420049,0.77853469 +2,1,8.176737731,0.424948854,0.77853469 +2,1,8.422147605,0.424379028,0.77853469 +2,1,8.70588182,0.423707714,0.77853469 +2,1,9.027284445,0.422932811,0.77853469 +2,1,9.385612468,0.422052556,0.77853469 +2,1,9.780037514,0.421065815,0.77853469 +2,1,10.20964776,0.419972455,0.77853469 +2,1,10.67345004,0.419400676,0.781531069 +2,1,10.86770694,0.418981957,0.758935311 +2,1,11.17037214,0.385839135,0.614478855 +2,1,11.6992653,0.335840083,0.498687801 +2,1,12.25890683,0.29191329,0.416354609 +2,1,12.84800295,0.253572514,0.351944846 +2,1,13.46519181,0.220278082,0.299832337 +2,1,14.10904661,0.191477908,0.256956606 +2,1,14.77807889,0.166631343,0.221322169 +2,1,15.470742,0.145236797,0.19150758 +2,1,16.18543466,0.126834289,0.166435523 +2,1,16.92050464,0.111011925,0.145263684 +2,1,17.67425264,0.097406118,0.127319849 +2,1,18.44493615,0.085699408,0.11206048 +2,1,19.23077353,0.075616912,0.099042189 +2,1,20.02994808,0.066922115,0.087901155 +2,1,20.8406123,0.059412477,0.078337446 +2,1,21.66089211,0.052915227,0.07010295 +2,1,22.4888912,0.04728299,0.062991402 +2,1,23.32269542,0.042390922,0.056831647 +2,1,24.1603772,0.038132739,0.05148062 +2,1,25,0.03441828,0.046818787 +2,1,25.02,0,0 +2,1,50,0,0 +2,5,0,0,0 +2,5,3,0.024680618,0.40876666 +2,5,3.54953237,0.112162126,0.396057646 +2,5,4.067900771,0.156108209,0.39320095 +2,5,4.553906848,0.180049935,0.394449372 +2,5,5.006427063,0.19380602,0.395387288 +2,5,5.424415288,0.202005082,0.396043345 +2,5,5.806905228,0.206989662,0.395929045 +2,5,6.153012649,0.210041846,0.39519265 +2,5,6.461937428,0.211893882,0.394126518 +2,5,6.732965398,0.212988948,0.392922592 +2,5,6.965470002,0.213596636,0.391683582 +2,5,7.158913742,0.213591753,0.389267345 +2,5,7.312849418,0.213430464,0.389267345 +2,5,7.426921164,0.21330898,0.389267345 +2,5,7.500865272,0.213229392,0.389267345 +2,5,7.534510799,0.213192979,0.389267345 +2,5,7.541241633,0.213185695,0.389267345 +2,5,7.58833327,0.213134413,0.389267345 +2,5,7.675676842,0.213038728,0.389267345 +2,5,7.803070431,0.212897651,0.389267345 +2,5,7.970219531,0.212710025,0.389267345 +2,5,8.176737731,0.212474427,0.389267345 +2,5,8.422147605,0.212189514,0.389267345 +2,5,8.70588182,0.211853857,0.389267345 +2,5,9.027284445,0.211466406,0.389267345 +2,5,9.385612468,0.211026278,0.389267345 +2,5,9.780037514,0.210532908,0.389267345 +2,5,10.20964776,0.209986228,0.389267345 +2,5,10.67345004,0.209700338,0.390765535 +2,5,10.86770694,0.209490979,0.379467656 +2,5,11.17037214,0.192919568,0.307239428 +2,5,11.6992653,0.167920042,0.249343901 +2,5,12.25890683,0.145956645,0.208177305 +2,5,12.84800295,0.126786257,0.175972423 +2,5,13.46519181,0.110139041,0.149916169 +2,5,14.10904661,0.095738954,0.128478303 +2,5,14.77807889,0.083315672,0.110661085 +2,5,15.470742,0.072618399,0.09575379 +2,5,16.18543466,0.063417145,0.083217762 +2,5,16.92050464,0.055505963,0.072631842 +2,5,17.67425264,0.048703059,0.063659925 +2,5,18.44493615,0.042849704,0.05603024 +2,5,19.23077353,0.037808456,0.049521095 +2,5,20.02994808,0.033461058,0.043950578 +2,5,20.8406123,0.029706239,0.039168723 +2,5,21.66089211,0.026457614,0.035051475 +2,5,22.4888912,0.023641495,0.031495701 +2,5,23.32269542,0.021195461,0.028415824 +2,5,24.1603772,0.01906637,0.02574031 +2,5,25,0.01720914,0.023409394 +2,5,25.02,0,0 +2,5,50,0,0 +4,1,0,0,0 +4,1,3,0.012340309,0.20438333 +4,1,3.54953237,0.056081063,0.198028823 +4,1,4.067900771,0.078054105,0.196600475 +4,1,4.553906848,0.090024968,0.197224686 +4,1,5.006427063,0.09690301,0.197693644 +4,1,5.424415288,0.101002541,0.198021673 +4,1,5.806905228,0.103494831,0.197964523 +4,1,6.153012649,0.105020923,0.197596325 +4,1,6.461937428,0.105946941,0.197063259 +4,1,6.732965398,0.106494474,0.196461296 +4,1,6.965470002,0.106798318,0.195841791 +4,1,7.158913742,0.106795876,0.194633673 +4,1,7.312849418,0.106715232,0.194633673 +4,1,7.426921164,0.10665449,0.194633673 +4,1,7.500865272,0.106614696,0.194633673 +4,1,7.534510799,0.106596489,0.194633673 +4,1,7.541241633,0.106592847,0.194633673 +4,1,7.58833327,0.106567207,0.194633673 +4,1,7.675676842,0.106519364,0.194633673 +4,1,7.803070431,0.106448826,0.194633673 +4,1,7.970219531,0.106355012,0.194633673 +4,1,8.176737731,0.106237214,0.194633673 +4,1,8.422147605,0.106094757,0.194633673 +4,1,8.70588182,0.105926929,0.194633673 +4,1,9.027284445,0.105733203,0.194633673 +4,1,9.385612468,0.105513139,0.194633673 +4,1,9.780037514,0.105266454,0.194633673 +4,1,10.20964776,0.104993114,0.194633673 +4,1,10.67345004,0.104850169,0.195382767 +4,1,10.86770694,0.104745489,0.189733828 +4,1,11.17037214,0.096459784,0.153619714 +4,1,11.6992653,0.083960021,0.12467195 +4,1,12.25890683,0.072978323,0.104088652 +4,1,12.84800295,0.063393129,0.087986212 +4,1,13.46519181,0.055069521,0.074958084 +4,1,14.10904661,0.047869477,0.064239152 +4,1,14.77807889,0.041657836,0.055330542 +4,1,15.470742,0.036309199,0.047876895 +4,1,16.18543466,0.031708572,0.041608881 +4,1,16.92050464,0.027752981,0.036315921 +4,1,17.67425264,0.02435153,0.031829962 +4,1,18.44493615,0.021424852,0.02801512 +4,1,19.23077353,0.018904228,0.024760547 +4,1,20.02994808,0.016730529,0.021975289 +4,1,20.8406123,0.014853119,0.019584362 +4,1,21.66089211,0.013228807,0.017525738 +4,1,22.4888912,0.011820748,0.015747851 +4,1,23.32269542,0.010597731,0.014207912 +4,1,24.1603772,0.009533185,0.012870155 +4,1,25,0.00860457,0.011704697 +4,1,25.02,0,0 +4,1,50,0,0 +4,5,0,0,0 +4,5,3,0.006170155,0.102191665 +4,5,3.54953237,0.028040532,0.099014412 +4,5,4.067900771,0.039027052,0.098300238 +4,5,4.553906848,0.045012484,0.098612343 +4,5,5.006427063,0.048451505,0.098846822 +4,5,5.424415288,0.050501271,0.099010836 +4,5,5.806905228,0.051747416,0.098982261 +4,5,6.153012649,0.052510462,0.098798163 +4,5,6.461937428,0.052973471,0.09853163 +4,5,6.732965398,0.053247237,0.098230648 +4,5,6.965470002,0.053399159,0.097920896 +4,5,7.158913742,0.053397938,0.097316836 +4,5,7.312849418,0.053357616,0.097316836 +4,5,7.426921164,0.053327245,0.097316836 +4,5,7.500865272,0.053307348,0.097316836 +4,5,7.534510799,0.053298245,0.097316836 +4,5,7.541241633,0.053296424,0.097316836 +4,5,7.58833327,0.053283603,0.097316836 +4,5,7.675676842,0.053259682,0.097316836 +4,5,7.803070431,0.053224413,0.097316836 +4,5,7.970219531,0.053177506,0.097316836 +4,5,8.176737731,0.053118607,0.097316836 +4,5,8.422147605,0.053047379,0.097316836 +4,5,8.70588182,0.052963464,0.097316836 +4,5,9.027284445,0.052866602,0.097316836 +4,5,9.385612468,0.05275657,0.097316836 +4,5,9.780037514,0.052633227,0.097316836 +4,5,10.20964776,0.052496557,0.097316836 +4,5,10.67345004,0.052425085,0.097691384 +4,5,10.86770694,0.052372745,0.094866914 +4,5,11.17037214,0.048229892,0.076809857 +4,5,11.6992653,0.041980011,0.062335975 +4,5,12.25890683,0.036489161,0.052044326 +4,5,12.84800295,0.031696564,0.043993106 +4,5,13.46519181,0.02753476,0.037479042 +4,5,14.10904661,0.023934739,0.032119576 +4,5,14.77807889,0.020828918,0.027665271 +4,5,15.470742,0.0181546,0.023938448 +4,5,16.18543466,0.015854286,0.020804441 +4,5,16.92050464,0.013876491,0.018157961 +4,5,17.67425264,0.012175765,0.015914981 +4,5,18.44493615,0.010712426,0.01400756 +4,5,19.23077353,0.009452114,0.012380274 +4,5,20.02994808,0.008365265,0.010987645 +4,5,20.8406123,0.00742656,0.009792181 +4,5,21.66089211,0.006614404,0.008762869 +4,5,22.4888912,0.005910374,0.007873925 +4,5,23.32269542,0.005298865,0.007103956 +4,5,24.1603772,0.004766593,0.006435078 +4,5,25,0.004302285,0.005852349 +4,5,25.02,0,0 +4,5,50,0,0 diff --git a/floris/turbine_library/iea_15MW_multi_dim_cp_ct.yaml b/floris/turbine_library/iea_15MW_multi_dim_cp_ct.yaml new file mode 100644 index 000000000..51e0a83f6 --- /dev/null +++ b/floris/turbine_library/iea_15MW_multi_dim_cp_ct.yaml @@ -0,0 +1,11 @@ +turbine_type: 'iea_15MW_multi_dim_cp_ct' +generator_efficiency: 1.0 +hub_height: 150.0 +pP: 1.88 +pT: 1.88 +rotor_diameter: 242.24 +TSR: 8.0 +ref_density_cp_ct: 1.225 +ref_tilt_cp_ct: 6.0 +multi_dimensional_cp_ct: True +power_thrust_data_file: '../floris/turbine_library/iea_15MW_multi_dim_Tp_Hs.csv' diff --git a/floris/turbine_library/nrel_5MW.yaml b/floris/turbine_library/nrel_5MW.yaml index 70c34a9f4..653ef14c7 100644 --- a/floris/turbine_library/nrel_5MW.yaml +++ b/floris/turbine_library/nrel_5MW.yaml @@ -197,3 +197,16 @@ power_thrust_table: - 25.01 - 25.02 - 50.0 + +### +# A boolean flag used when the user wants FLORIS to use the user-supplied multi-dimensional +# Cp/Ct information. +multi_dimensional_cp_ct: False + +### +# The path to the .csv file that contains the multi-dimensional Cp/Ct data. The format of this +# file is such that any external conditions, such as wave height or wave period, that the +# Cp/Ct data is dependent on come first, in column format. The last three columns of the .csv +# file must be ``ws``, ``Cp``, and ``Ct``, in that order. An example of fictional data is given +# in ``floris/turbine_library/iea_15MW_multi_dim_Tp_Hs.csv``. +power_thrust_data_file: '../floris/turbine_library/iea_15MW_multi_dim_Tp_Hs.csv' diff --git a/setup.py b/setup.py index d3af38d21..0bab76eb1 100644 --- a/setup.py +++ b/setup.py @@ -42,6 +42,7 @@ # utilities "coloredlogs>=10.0", + "flatten_dict", ] # What packages are optional? diff --git a/tests/conftest.py b/tests/conftest.py index efa4fd13c..1c96948f8 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -361,6 +361,11 @@ def __init__(self): } self.turbine_floating["floating_correct_cp_ct_for_tilt"] = True + self.turbine_multi_dim = copy.deepcopy(self.turbine) + del self.turbine_multi_dim['power_thrust_table'] + self.turbine_multi_dim["multi_dimensional_cp_ct"] = True + self.turbine_multi_dim["power_thrust_data_file"] = "" + self.farm = { "layout_x": X_COORDS, "layout_y": Y_COORDS, diff --git a/tests/turbine_multi_dim_unit_test.py b/tests/turbine_multi_dim_unit_test.py new file mode 100644 index 000000000..fd6cdacce --- /dev/null +++ b/tests/turbine_multi_dim_unit_test.py @@ -0,0 +1,318 @@ +# 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 + + +from pathlib import Path + +import numpy as np +import pandas as pd +import pytest +from scipy.interpolate import interp1d + +from floris.simulation import ( + Turbine, + TurbineMultiDimensional, +) +from floris.simulation.turbine_multi_dim import ( + axial_induction_multidim, + Ct_multidim, + multidim_Ct_down_select, + multidim_power_down_select, + MultiDimensionalPowerThrustTable, + power_multidim, +) +from tests.conftest import SampleInputs, WIND_SPEEDS + + +TEST_DATA = Path(__file__).resolve().parent.parent / "floris" / "turbine_library" +CSV_INPUT = TEST_DATA / "iea_15MW_multi_dim_Tp_Hs.csv" + + +# size 3 x 4 x 1 x 1 x 1 +WIND_CONDITION_BROADCAST = np.stack( + ( + np.reshape(np.array(WIND_SPEEDS), (-1, 1, 1, 1)), # Wind direction 0 + np.reshape(np.array(WIND_SPEEDS), (-1, 1, 1, 1)), # Wind direction 1 + np.reshape(np.array(WIND_SPEEDS), (-1, 1, 1, 1)), # Wind direction 2 + ), + axis=0, +) +INDEX_FILTER = [0, 2] + + +def test_multidim_Ct_down_select(): + CONDITIONS = {'Tp': 2, 'Hs': 1} + + turbine_data = SampleInputs().turbine_multi_dim + turbine_data["power_thrust_data_file"] = CSV_INPUT + turbine = TurbineMultiDimensional.from_dict(turbine_data) + turbine_type_map = np.array([turbine.turbine_type]) + turbine_type_map = turbine_type_map[None, None, :] + + downselect_turbine_fCts = multidim_Ct_down_select([[[turbine.fCt_interp]]], CONDITIONS) + + assert downselect_turbine_fCts == turbine.fCt_interp[(2, 1)] + + +def test_multidim_power_down_select(): + CONDITIONS = {'Tp': 2, 'Hs': 1} + + turbine_data = SampleInputs().turbine_multi_dim + turbine_data["power_thrust_data_file"] = CSV_INPUT + turbine = TurbineMultiDimensional.from_dict(turbine_data) + turbine_type_map = np.array([turbine.turbine_type]) + turbine_type_map = turbine_type_map[None, None, :] + + downselect_power_interps = multidim_power_down_select([[[turbine.power_interp]]], CONDITIONS) + + assert downselect_power_interps == turbine.power_interp[(2, 1)] + + +def test_multi_dimensional_power_thrust_table(): + turbine_data = SampleInputs().turbine_multi_dim + turbine_data["power_thrust_data_file"] = CSV_INPUT + df_data = pd.read_csv(turbine_data["power_thrust_data_file"]) + flattened_dict = MultiDimensionalPowerThrustTable.from_dataframe(df_data) + flattened_dict_base = { + ('Tp', '2', 'Hs', '1'): [], + ('Tp', '2', 'Hs', '5'): [], + ('Tp', '4', 'Hs', '1'): [], + ('Tp', '4', 'Hs', '5'): [], + } + assert flattened_dict == flattened_dict_base + + # Test for initialization errors + for el in ("ws", "Cp", "Ct"): + df_data = pd.read_csv(turbine_data["power_thrust_data_file"]) + df = df_data.drop(el, axis=1) + with pytest.raises(ValueError): + MultiDimensionalPowerThrustTable.from_dataframe(df) + + +def test_turbine_init(): + turbine_data = SampleInputs().turbine_multi_dim + turbine_data["power_thrust_data_file"] = CSV_INPUT + turbine = TurbineMultiDimensional.from_dict(turbine_data) + assert turbine.rotor_diameter == turbine_data["rotor_diameter"] + assert turbine.hub_height == turbine_data["hub_height"] + assert turbine.pP == turbine_data["pP"] + assert turbine.pT == turbine_data["pT"] + assert turbine.generator_efficiency == turbine_data["generator_efficiency"] + + assert isinstance(turbine.power_thrust_data, dict) + assert isinstance(turbine.fCp_interp, interp1d) + assert isinstance(turbine.fCt_interp, dict) + assert isinstance(turbine.power_interp, dict) + assert turbine.rotor_radius == turbine_data["rotor_diameter"] / 2.0 + + +def test_ct(): + N_TURBINES = 4 + + turbine_data = SampleInputs().turbine_multi_dim + turbine_data["power_thrust_data_file"] = CSV_INPUT + turbine = TurbineMultiDimensional.from_dict(turbine_data) + turbine_type_map = np.array(N_TURBINES * [turbine.turbine_type]) + turbine_type_map = turbine_type_map[None, None, :] + + # Single turbine + # yaw angle / fCt are (n wind direction, n wind speed, n turbine) + wind_speed = 10.0 + thrust = Ct_multidim( + velocities=wind_speed * np.ones((1, 1, 1, 3, 3)), + yaw_angle=np.zeros((1, 1, 1)), + tilt_angle=np.ones((1, 1, 1)) * 5.0, + ref_tilt_cp_ct=np.ones((1, 1, 1)) * 5.0, + fCt=np.array([[[turbine.fCt_interp[(2, 1)]]]]), + tilt_interp=np.array([(turbine.turbine_type, None)]), + correct_cp_ct_for_tilt=np.array([[[False]]]), + turbine_type_map=turbine_type_map[:,:,0] + ) + + print(thrust) + np.testing.assert_allclose(thrust, np.array([[[0.77853469]]])) + + # Multiple turbines with index filter + # 4 turbines with 3 x 3 grid arrays + thrusts = Ct_multidim( + velocities=np.ones((N_TURBINES, 3, 3)) * WIND_CONDITION_BROADCAST, # 3 x 4 x 4 x 3 x 3 + yaw_angle=np.zeros((1, 1, N_TURBINES)), + tilt_angle=np.ones((1, 1, N_TURBINES)) * 5.0, + ref_tilt_cp_ct=np.ones((1, 1, N_TURBINES)) * 5.0, + fCt=np.tile( + [turbine.fCt_interp[(2, 1)]], + ( + np.shape(WIND_CONDITION_BROADCAST)[0], + np.shape(WIND_CONDITION_BROADCAST)[1], + N_TURBINES, + ) + ), + tilt_interp=np.array([(turbine.turbine_type, None)]), + correct_cp_ct_for_tilt=np.array([[[False] * N_TURBINES]]), + turbine_type_map=turbine_type_map, + ix_filter=INDEX_FILTER, + ) + assert len(thrusts[0, 0]) == len(INDEX_FILTER) + + print(thrusts) + + thrusts_truth = [ + [ + [0.77853469, 0.77853469], + [0.77853469, 0.77853469], + [0.77853469, 0.77853469], + [0.6957943, 0.6957943 ], + ], + [ + [0.77853469, 0.77853469], + [0.77853469, 0.77853469], + [0.77853469, 0.77853469], + [0.6957943, 0.6957943 ], + ], + [ + [0.77853469, 0.77853469], + [0.77853469, 0.77853469], + [0.77853469, 0.77853469], + [0.6957943, 0.6957943 ], + ], + ] + + np.testing.assert_allclose(thrusts, thrusts_truth) + + +def test_power(): + N_TURBINES = 4 + AIR_DENSITY = 1.225 + + turbine_data = SampleInputs().turbine_multi_dim + turbine_data["power_thrust_data_file"] = CSV_INPUT + turbine = TurbineMultiDimensional.from_dict(turbine_data) + turbine_type_map = np.array(N_TURBINES * [turbine.turbine_type]) + turbine_type_map = turbine_type_map[None, None, :] + + # Single turbine + wind_speed = 10.0 + p = power_multidim( + ref_density_cp_ct=AIR_DENSITY, + rotor_effective_velocities=wind_speed * np.ones((1, 1, 1, 3, 3)), + power_interp=np.array([[[turbine.power_interp[(2, 1)]]]]), + ) + + power_truth = [ + [ + [ + [ + [3215682.686486, 3215682.686486, 3215682.686486], + [3215682.686486, 3215682.686486, 3215682.686486], + [3215682.686486, 3215682.686486, 3215682.686486], + ] + ] + ] + ] + + np.testing.assert_allclose(p, power_truth ) + + # Multiple turbines with ix filter + rotor_effective_velocities = np.ones((N_TURBINES, 3, 3)) * WIND_CONDITION_BROADCAST + p = power_multidim( + ref_density_cp_ct=AIR_DENSITY, + rotor_effective_velocities=rotor_effective_velocities, + power_interp=np.tile( + [turbine.power_interp[(2, 1)]], + ( + np.shape(WIND_CONDITION_BROADCAST)[0], + np.shape(WIND_CONDITION_BROADCAST)[1], + N_TURBINES, + ) + ), + ix_filter=INDEX_FILTER, + ) + assert len(p[0, 0]) == len(INDEX_FILTER) + + unique_power = turbine.power_interp[(2, 1)]( + np.unique(rotor_effective_velocities) + ) * AIR_DENSITY + + power_truth = np.zeros_like(rotor_effective_velocities) + for i in range(3): + for j in range(4): + for k in range(4): + for m in range(3): + for n in range(3): + power_truth[i, j, k, m, n] = unique_power[j] + + np.testing.assert_allclose(p, power_truth[:, :, INDEX_FILTER[0]:INDEX_FILTER[1], :, :]) + + +def test_axial_induction(): + + N_TURBINES = 4 + + turbine_data = SampleInputs().turbine_multi_dim + turbine_data["power_thrust_data_file"] = CSV_INPUT + turbine = TurbineMultiDimensional.from_dict(turbine_data) + turbine_type_map = np.array(N_TURBINES * [turbine.turbine_type]) + turbine_type_map = turbine_type_map[None, None, :] + + baseline_ai = 0.2646995 + + # Single turbine + wind_speed = 10.0 + ai = axial_induction_multidim( + velocities=wind_speed * np.ones((1, 1, 1, 3, 3)), + yaw_angle=np.zeros((1, 1, 1)), + tilt_angle=np.ones((1, 1, 1)) * 5.0, + ref_tilt_cp_ct=np.ones((1, 1, 1)) * 5.0, + fCt=np.array([[[turbine.fCt_interp[(2, 1)]]]]), + tilt_interp=np.array([(turbine.turbine_type, None)]), + correct_cp_ct_for_tilt=np.array([[[False]]]), + turbine_type_map=turbine_type_map[0,0,0], + ) + np.testing.assert_allclose(ai, baseline_ai) + + # Multiple turbines with ix filter + ai = axial_induction_multidim( + velocities=np.ones((N_TURBINES, 3, 3)) * WIND_CONDITION_BROADCAST, # 3 x 4 x 4 x 3 x 3 + yaw_angle=np.zeros((1, 1, N_TURBINES)), + tilt_angle=np.ones((1, 1, N_TURBINES)) * 5.0, + ref_tilt_cp_ct=np.ones((1, 1, N_TURBINES)) * 5.0, + fCt=np.tile( + [turbine.fCt_interp[(2, 1)]], + ( + np.shape(WIND_CONDITION_BROADCAST)[0], + np.shape(WIND_CONDITION_BROADCAST)[1], + N_TURBINES, + ) + ), + tilt_interp=np.array([(turbine.turbine_type, None)] * N_TURBINES), + correct_cp_ct_for_tilt=np.array([[[False] * N_TURBINES]]), + turbine_type_map=turbine_type_map, + ix_filter=INDEX_FILTER, + ) + + assert len(ai[0, 0]) == len(INDEX_FILTER) + + # Test the 10 m/s wind speed to use the same baseline as above + np.testing.assert_allclose(ai[0,2], baseline_ai) + + +def test_asdict(sample_inputs_fixture: SampleInputs): + + turbine = Turbine.from_dict(sample_inputs_fixture.turbine) + dict1 = turbine.as_dict() + + new_turb = Turbine.from_dict(dict1) + dict2 = new_turb.as_dict() + + assert dict1 == dict2