From 33c848b58f4798e731e2e1a68182cd6112f3fa4d Mon Sep 17 00:00:00 2001 From: Rob Hammond <13874373+RHammond2@users.noreply.github.com> Date: Mon, 19 Sep 2022 11:50:33 -0600 Subject: [PATCH 01/20] generalize the x, y, z mean calc and start sequential solver class --- floris/simulation/solver.py | 76 ++++++++++++++++++++----------------- 1 file changed, 42 insertions(+), 34 deletions(-) diff --git a/floris/simulation/solver.py b/floris/simulation/solver.py index 255e49fcf..622760238 100644 --- a/floris/simulation/solver.py +++ b/floris/simulation/solver.py @@ -11,11 +11,13 @@ # the License. import copy +from typing import Type + +import attrs import numpy as np -import time -import sys +from attrs import define, field -from floris.simulation import Farm +from floris.simulation import Farm, flow_field from floris.simulation import TurbineGrid, FlowFieldGrid from floris.simulation import Ct, axial_induction from floris.simulation import FlowField @@ -26,6 +28,24 @@ wake_added_yaw, yaw_added_turbulence_mixing ) +from floris.type_dec import NDArrayFloat + + +def _get_mean_grid(grid: TurbineGrid | FlowFieldGrid, i: int) -> tuple[NDArrayFloat, NDArrayFloat, NDArrayFloat]: + """Calculates the mean of the `grid`'s `x_sorted`, `y_sorted`, and `z_sorted` attributes and + returns them as a new 5-dimensional object to align with expected internal structures. + + Args: + grid (TurbineGrid | FlowFieldGrid): The solver's grid object. + i (int): The turbine index that the mean x, y, and z will be computed over. + + Returns: + tuple[NDArrayFloat, NDArrayFloat, NDArrayFloat]: The mean over axis 3 and 4 at turbine i. + """ + x_i = np.mean(grid.x_sorted[:, :, i:i+1], axis=(3, 4))[:, :, :, None, None] + y_i = np.mean(grid.y_sorted[:, :, i:i+1], axis=(3, 4))[:, :, :, None, None] + z_i = np.mean(grid.z_sorted[:, :, i:i+1], axis=(3, 4))[:, :, :, None, None] + return x_i, y_i, z_i def calculate_area_overlap(wake_velocities, freestream_velocities, y_ngrid, z_ngrid): @@ -43,7 +63,20 @@ def calculate_area_overlap(wake_velocities, freestream_velocities, y_ngrid, z_ng return np.sum(freestream_velocities - wake_velocities > 0.05, axis=(3, 4)) / (y_ngrid * z_ngrid) -# @profile +@define(auto_attribs=True) +class SequentialSolver: + farm: Farm = field(validator=attrs.validators.instance_of(Farm)) + flow_field: FlowField = field(validator=attrs.validators.instance_of(Farm)) + grid: TurbineGrid | FlowFieldGrid = field(validator=attrs.validators.instance_of((FlowFieldGrid, TurbineGrid))) + model_manager: WakeModelManager = field(validator=attrs.validators.instance_of(WakeModelManager)) + full_flow: bool = field(validator=attrs.validators.instance_of(bool)) + + def _full_flow_init(self): + if not isinstance(grid, FlowFieldGrid): + raise TypeError("When `full_flow` is True, `grid` must be a `FlowFieldGrid`, not `TurbineGrid`.") + + + def sequential_solver(farm: Farm, flow_field: FlowField, grid: TurbineGrid, model_manager: WakeModelManager) -> None: # Algorithm # For each turbine, calculate its effect on every downstream turbine. @@ -67,12 +100,7 @@ def sequential_solver(farm: Farm, flow_field: FlowField, grid: TurbineGrid, mode 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] + x_i, y_i, z_i = _get_mean_grid(grid, i) u_i = flow_field.u_sorted[:, :, i:i+1] v_i = flow_field.v_sorted[:, :, i:i+1] @@ -257,12 +285,7 @@ def full_flow_sequential_solver(farm: Farm, flow_field: FlowField, flow_field_gr for i in range(flow_field_grid.n_turbines): # Get the current turbine quantities - x_i = np.mean(turbine_grid.x_sorted[:, :, i:i+1], axis=(3, 4)) - x_i = x_i[:, :, :, None, None] - y_i = np.mean(turbine_grid.y_sorted[:, :, i:i+1], axis=(3, 4)) - y_i = y_i[:, :, :, None, None] - z_i = np.mean(turbine_grid.z_sorted[:, :, i:i+1], axis=(3, 4)) - z_i = z_i[:, :, :, None, None] + x_i, y_i, z_i = _get_mean_grid(turbine_grid, i) u_i = turbine_grid_flow_field.u_sorted[:, :, i:i+1] v_i = turbine_grid_flow_field.v_sorted[:, :, i:i+1] @@ -386,12 +409,7 @@ def cc_solver(farm: Farm, flow_field: FlowField, grid: TurbineGrid, model_manage 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] + x_i, y_i, z_i = _get_mean_grid(grid, i) mask2 = np.array(grid.x_sorted < x_i + 0.01) * np.array(grid.x_sorted > x_i - 0.01) * np.array(grid.y_sorted < y_i + 0.51*126.0) * np.array(grid.y_sorted > y_i - 0.51*126.0) # mask2 = np.logical_and(np.logical_and(np.logical_and(grid.x_sorted < x_i + 0.01, grid.x_sorted > x_i - 0.01), grid.y_sorted < y_i + 0.51*126.0), grid.y_sorted > y_i - 0.51*126.0) @@ -590,12 +608,7 @@ def full_flow_cc_solver(farm: Farm, flow_field: FlowField, flow_field_grid: Flow for i in range(flow_field_grid.n_turbines): # Get the current turbine quantities - x_i = np.mean(turbine_grid.x_sorted[:, :, i:i+1], axis=(3, 4)) - x_i = x_i[:, :, :, None, None] - y_i = np.mean(turbine_grid.y_sorted[:, :, i:i+1], axis=(3, 4)) - y_i = y_i[:, :, :, None, None] - z_i = np.mean(turbine_grid.z_sorted[:, :, i:i+1], axis=(3, 4)) - z_i = z_i[:, :, :, None, None] + x_i, y_i, z_i = _get_mean_grid(turbine_grid, i) u_i = turbine_grid_flow_field.u_sorted[:, :, i:i+1] v_i = turbine_grid_flow_field.v_sorted[:, :, i:i+1] @@ -718,12 +731,7 @@ def turbopark_solver(farm: Farm, flow_field: FlowField, grid: TurbineGrid, model # 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] + x_i, y_i, z_i = _get_mean_grid(grid, i) u_i = flow_field.u_sorted[:, :, i:i+1] v_i = flow_field.v_sorted[:, :, i:i+1] From 082722d4ebd7116680f565260e9d6ed4624c4501 Mon Sep 17 00:00:00 2001 From: Rob Hammond <13874373+RHammond2@users.noreply.github.com> Date: Mon, 19 Sep 2022 17:03:21 -0600 Subject: [PATCH 02/20] work towards base solver class, complete sequential solver, and start cc solver --- floris/simulation/solver.py | 427 +++++++++++++++++++++++++++++++++++- 1 file changed, 417 insertions(+), 10 deletions(-) diff --git a/floris/simulation/solver.py b/floris/simulation/solver.py index 622760238..1008bd49f 100644 --- a/floris/simulation/solver.py +++ b/floris/simulation/solver.py @@ -10,6 +10,7 @@ # License for the specific language governing permissions and limitations under # the License. +from abc import abstractmethod import copy from typing import Type @@ -64,16 +65,422 @@ def calculate_area_overlap(wake_velocities, freestream_velocities, y_ngrid, z_ng @define(auto_attribs=True) -class SequentialSolver: - farm: Farm = field(validator=attrs.validators.instance_of(Farm)) - flow_field: FlowField = field(validator=attrs.validators.instance_of(Farm)) - grid: TurbineGrid | FlowFieldGrid = field(validator=attrs.validators.instance_of((FlowFieldGrid, TurbineGrid))) - model_manager: WakeModelManager = field(validator=attrs.validators.instance_of(WakeModelManager)) - full_flow: bool = field(validator=attrs.validators.instance_of(bool)) - - def _full_flow_init(self): - if not isinstance(grid, FlowFieldGrid): - raise TypeError("When `full_flow` is True, `grid` must be a `FlowFieldGrid`, not `TurbineGrid`.") +class Solver: + farm: Farm = field(converter=copy.deepcopy, validator=attrs.validators.instance_of(Farm)) + flow_field: FlowField = field(converter=copy.deepcopy, validator=attrs.validators.instance_of(Farm)) + grid: TurbineGrid | FlowFieldGrid = field(converter=copy.deepcopy, validator=attrs.validators.instance_of((FlowFieldGrid, TurbineGrid))) + model_manager: WakeModelManager = field(converter=copy.deepcopy, validator=attrs.validators.instance_of(WakeModelManager)) + + @abstractmethod + def solve(self, *, full_flow: bool = False, farm: Farm, flow_field: FlowFieldGrid = None, grid: TurbineGrid | FlowFieldGrid = None) -> None: + pass + + def full_flow_solve(self): + """Initializes all the additional attributes to compute the full flow field, then runs the + sequential solver with a turbine grid, then again with the `self.grid` object. + + Raises: + TypeError: Raised if initialized `grid` value is not a `FlowFieldGrid` object. + """ + if not isinstance(self.grid, FlowFieldGrid): + raise TypeError("Cannot run `full_flow_solve` with a `TurbineGrid` object, reinitialize with the input to `grid` being a `FlowFieldGrid` object.") + + + # TODO: Why are we making copies of the originals, and then never using the original, + # can we just use the originally provided values and get on with the calculating? + farm = copy.deepcopy(self.farm) + flow_field = copy.deepcopy(self.flow_field) + + farm.construct_turbine_map() + farm.construct_turbine_fCts() + farm.construct_turbine_fCps() + farm.construct_turbine_power_interps() + farm.construct_hub_heights() + farm.construct_rotor_diameters() + farm.construct_turbine_TSRs() + farm.construc_turbine_pPs() + farm.construc_turbine_ref_density_cp_cts() + farm.construct_coordinates() + + turbine_grid = TurbineGrid( + turbine_coordinates=farm.coordinates, + reference_turbine_diameter=farm.rotor_diameters, + wind_directions=flow_field.wind_directions, + wind_speeds=flow_field.wind_speeds, + grid_resolution=3, + time_series=flow_field.time_series, + ) + farm.expand_farm_properties( + flow_field.n_wind_directions, flow_field.n_wind_speeds, turbine_grid.sorted_coord_indices + ) + flow_field.initialize_velocity_field(turbine_grid) + farm.initialize(turbine_grid.sorted_indices) + self.solve(full_flow=False, farm=farm, flow_field=flow_field, grid=turbine_grid) + self.solve(full_flow=True, farm=farm, flow_field=flow_field, grid=turbine_grid) + + +@define(auto_attribs=True) +class SequentialSolver(Solver): + farm: Farm = field(converter=copy.deepcopy, validator=attrs.validators.instance_of(Farm)) + flow_field: FlowField = field(converter=copy.deepcopy, validator=attrs.validators.instance_of(Farm)) + grid: TurbineGrid | FlowFieldGrid = field(converter=copy.deepcopy, validator=attrs.validators.instance_of((FlowFieldGrid, TurbineGrid))) + model_manager: WakeModelManager = field(converter=copy.deepcopy, validator=attrs.validators.instance_of(WakeModelManager)) + + def solve(self, *, full_flow: bool = False, grid: TurbineGrid | FlowFieldGrid = None) -> None: + """Runs the sequential sover methodology, or full flow sequential solver methodology for a + wind farm. + + Args: + full_flow (bool, optional): Runs the full flow solver when True, and the standard + sequential solver, when False. Defaults to False. + grid (TurbineGrid | FlowFieldGrid, optional): Allows for a non-initialized `grid` object + to be used. It should be noted that this functionality is intended for use with + `full_flow_solve`, which computes over a `TurbineGrid`, then a `FlowFieldGrid`. Defaults to None. + """ + if farm is None: + farm = self.farm + if flow_field is None: + flow_field = self.flow_field + if grid is None: + grid = self.grid + + gch_gain = 2 + + deflection_model_args = self.model_manager.deflection_model.prepare_function(grid, flow_field) + deficit_model_args = self.model_manager.velocity_model.prepare_function(grid, flow_field) + + 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) + + if not full_flow: + 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, y_i, z_i = _get_mean_grid(grid, i) + u_i = flow_field.u_sorted[:, :, i:i+1] + v_i = flow_field.v_sorted[:, :, i:i+1] + + # Since we are filtering for the ith turbine in the Ct function, get the first index here (0:1) + ct_i = Ct( + velocities=flow_field.u_sorted, + yaw_angle=farm.yaw_angles_sorted, + fCt=farm.turbine_fCts, + turbine_type_map=farm.turbine_type_map_sorted, + ix_filter=[i], + )[:, :, 0:1, None, None] + + axial_induction_i = axial_induction( + velocities=flow_field.u_sorted, + yaw_angle=farm.yaw_angles_sorted, + fCt=farm.turbine_fCts, + turbine_type_map=farm.turbine_type_map_sorted, + ix_filter=[i], + )[:, :, 0:1, None, None] # Since we are filtering for the i'th turbine in the axial induction function, get the first index here (0:1) + turbulence_intensity_i = flow_field.turbulence_intensity_field[:, :, 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 self.model_manager.enable_secondary_steering: + effective_yaw_i += 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 + ) + + deflection_field = self.model_manager.deflection_model.function( + x_i, + y_i, + effective_yaw_i, + turbulence_intensity_i, + ct_i, + rotor_diameter_i, + **deflection_model_args + ) + + if self.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 + ) + if not full_flow: + if self.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], + ) + turbine_turbulence_intensity[:, :, i:i+1] = turbulence_intensity_i + gch_gain * I_mixing + + # NOTE: exponential + velocity_deficit = self.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 = self.model_manager.combination_model.function( + wake_field, + velocity_deficit * flow_field.u_initial_sorted + ) + + if not full_flow: + wake_added_turbulence_intensity = self.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) + * np.array(grid.x_sorted > x_i) + * np.array(np.abs(y_i - grid.y_sorted) < 2 * rotor_diameter_i) + * np.array(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 + + if not full_flow: + flow_field.turbulence_intensity_field = np.mean(turbine_turbulence_intensity, axis=(3, 4)) + flow_field.turbulence_intensity_field = flow_field.turbulence_intensity_field[:, :, :, None, None] + + +@define(auto_attribs=True) +class CCSolver(Solver): + farm: Farm = field(converter=copy.deepcopy, validator=attrs.validators.instance_of(Farm)) + flow_field: FlowField = field(converter=copy.deepcopy, validator=attrs.validators.instance_of(Farm)) + grid: TurbineGrid | FlowFieldGrid = field(converter=copy.deepcopy, validator=attrs.validators.instance_of((FlowFieldGrid, TurbineGrid))) + model_manager: WakeModelManager = field(converter=copy.deepcopy, validator=attrs.validators.instance_of(WakeModelManager)) + + def solve(self, *, full_flow: bool = False, grid: TurbineGrid | FlowFieldGrid = None) -> None: + + #****** + # TODO: Have only added self. to the model_manager, and not compared to the full flow verion + #****** + + # <> + deflection_model_args = self.model_manager.deflection_model.prepare_function(grid, flow_field) + deficit_model_args = self.model_manager.velocity_model.prepare_function(grid, flow_field) + + # This is u_wake + v_wake = np.zeros_like(flow_field.v_initial_sorted) + w_wake = np.zeros_like(flow_field.w_initial_sorted) + turb_u_wake = np.zeros_like(flow_field.u_initial_sorted) + turb_inflow_field = copy.deepcopy(flow_field.u_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 + + shape = (farm.n_turbines,) + np.shape(flow_field.u_initial_sorted) + Ctmp = np.zeros((shape)) + # Ctmp = np.zeros((len(x_coord), len(wd), len(ws), len(x_coord), y_ngrid, z_ngrid)) + + sigma_i = np.zeros((shape)) + # sigma_i = np.zeros((len(x_coord), len(wd), len(ws), len(x_coord), y_ngrid, z_ngrid)) + + # Calculate the velocity deficit sequentially from upstream to downstream turbines + for i in range(grid.n_turbines): + + # Get the current turbine quantities + x_i, y_i, z_i = _get_mean_grid(grid, i) + + mask2 = np.array(grid.x_sorted < x_i + 0.01) * np.array(grid.x_sorted > x_i - 0.01) * np.array(grid.y_sorted < y_i + 0.51*126.0) * np.array(grid.y_sorted > y_i - 0.51*126.0) + # mask2 = np.logical_and(np.logical_and(np.logical_and(grid.x_sorted < x_i + 0.01, grid.x_sorted > x_i - 0.01), grid.y_sorted < y_i + 0.51*126.0), grid.y_sorted > y_i - 0.51*126.0) + turb_inflow_field = turb_inflow_field * ~mask2 + (flow_field.u_initial_sorted - turb_u_wake) * mask2 + + turb_avg_vels = average_velocity(turb_inflow_field) + turb_Cts = Ct( + turb_avg_vels, + farm.yaw_angles_sorted, + farm.turbine_fCts, + turbine_type_map=farm.turbine_type_map_sorted, + ) + turb_Cts = turb_Cts[:, :, :, None, None] + turb_aIs = axial_induction( + turb_avg_vels, + farm.yaw_angles_sorted, + farm.turbine_fCts, + turbine_type_map=farm.turbine_type_map_sorted, + ix_filter=[i], + ) + turb_aIs = turb_aIs[:, :, :, None, None] + + u_i = turb_inflow_field[:, :, i:i+1] + v_i = flow_field.v_sorted[:, :, i:i+1] + + axial_induction_i = axial_induction( + velocities=flow_field.u_sorted, + yaw_angle=farm.yaw_angles_sorted, + fCt=farm.turbine_fCts, + turbine_type_map=farm.turbine_type_map_sorted, + ix_filter=[i], + ) + + axial_induction_i = axial_induction_i[:, :, :, 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, + turb_Cts[:, :, i:i+1], + TSR_i, + axial_induction_i, + scale=2.0, + ) + effective_yaw_i += added_yaw + + # Model calculations + # NOTE: exponential + deflection_field = self.model_manager.deflection_model.function( + x_i, + y_i, + effective_yaw_i, + turbulence_intensity_i, + turb_Cts[:, :, i:i+1], + rotor_diameter_i, + **deflection_model_args + ) + + if self.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, + turb_Cts[:, :, i:i+1], + TSR_i, + axial_induction_i, + scale=2.0 + ) + + if self.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 = 1.0 + turbine_turbulence_intensity[:, :, i:i+1] = turbulence_intensity_i + gch_gain * I_mixing + + turb_u_wake, Ctmp = self.model_manager.velocity_model.function( + i, + x_i, + y_i, + z_i, + u_i, + deflection_field, + yaw_angle_i, + turbine_turbulence_intensity, + turb_Cts, + farm.rotor_diameters_sorted[:, :, :, None, None], + turb_u_wake, + Ctmp, + **deficit_model_args + ) + + wake_added_turbulence_intensity = self.model_manager.turbulence_model.function( + ambient_turbulence_intensity, + grid.x_sorted, + x_i, + rotor_diameter_i, + turb_aIs + ) + + # Calculate wake overlap for wake-added turbulence (WAT) + area_overlap = 1 - np.sum(turb_u_wake <= 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) + * np.array(grid.x_sorted > x_i) + * np.array(np.abs(y_i - grid.y_sorted) < 2 * rotor_diameter_i) + * np.array(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.v_sorted += v_wake + flow_field.w_sorted += w_wake + flow_field.u_sorted = turb_inflow_field + + flow_field.turbulence_intensity_field = np.mean(turbine_turbulence_intensity, axis=(3,4)) + flow_field.turbulence_intensity_field = flow_field.turbulence_intensity_field[:,:,:,None,None] From 1f61279234ac99b9358e4a345b6f401f7ece155d Mon Sep 17 00:00:00 2001 From: Rob Hammond <13874373+RHammond2@users.noreply.github.com> Date: Tue, 20 Sep 2022 11:20:22 -0600 Subject: [PATCH 03/20] refactor cc solver --- floris/simulation/solver.py | 170 +++++++++++++++++++++--------------- 1 file changed, 99 insertions(+), 71 deletions(-) diff --git a/floris/simulation/solver.py b/floris/simulation/solver.py index 1008bd49f..646370f36 100644 --- a/floris/simulation/solver.py +++ b/floris/simulation/solver.py @@ -12,7 +12,6 @@ from abc import abstractmethod import copy -from typing import Type import attrs import numpy as np @@ -72,7 +71,14 @@ class Solver: model_manager: WakeModelManager = field(converter=copy.deepcopy, validator=attrs.validators.instance_of(WakeModelManager)) @abstractmethod - def solve(self, *, full_flow: bool = False, farm: Farm, flow_field: FlowFieldGrid = None, grid: TurbineGrid | FlowFieldGrid = None) -> None: + def solve( + self, + *, + full_flow: bool = False, + farm: Farm, + flow_field: FlowFieldGrid = None, + grid: TurbineGrid | FlowFieldGrid = None + ) -> None: pass def full_flow_solve(self): @@ -126,7 +132,14 @@ class SequentialSolver(Solver): grid: TurbineGrid | FlowFieldGrid = field(converter=copy.deepcopy, validator=attrs.validators.instance_of((FlowFieldGrid, TurbineGrid))) model_manager: WakeModelManager = field(converter=copy.deepcopy, validator=attrs.validators.instance_of(WakeModelManager)) - def solve(self, *, full_flow: bool = False, grid: TurbineGrid | FlowFieldGrid = None) -> None: + def solve( + self, + *, + full_flow: bool = False, + farm: Farm, + flow_field: FlowFieldGrid = None, + grid: TurbineGrid | FlowFieldGrid = None + ) -> None: """Runs the sequential sover methodology, or full flow sequential solver methodology for a wind farm. @@ -302,11 +315,21 @@ class CCSolver(Solver): grid: TurbineGrid | FlowFieldGrid = field(converter=copy.deepcopy, validator=attrs.validators.instance_of((FlowFieldGrid, TurbineGrid))) model_manager: WakeModelManager = field(converter=copy.deepcopy, validator=attrs.validators.instance_of(WakeModelManager)) - def solve(self, *, full_flow: bool = False, grid: TurbineGrid | FlowFieldGrid = None) -> None: + def solve( + self, + *, + full_flow: bool = False, + farm: Farm, + flow_field: FlowFieldGrid = None, + grid: TurbineGrid | FlowFieldGrid = None + ) -> None: - #****** - # TODO: Have only added self. to the model_manager, and not compared to the full flow verion - #****** + if farm is None: + farm = self.farm + if flow_field is None: + flow_field = self.flow_field + if grid is None: + grid = self.grid # <> deflection_model_args = self.model_manager.deflection_model.prepare_function(grid, flow_field) @@ -321,12 +344,7 @@ def solve(self, *, full_flow: bool = False, grid: TurbineGrid | FlowFieldGrid = 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 - shape = (farm.n_turbines,) + np.shape(flow_field.u_initial_sorted) - Ctmp = np.zeros((shape)) - # Ctmp = np.zeros((len(x_coord), len(wd), len(ws), len(x_coord), y_ngrid, z_ngrid)) - - sigma_i = np.zeros((shape)) - # sigma_i = np.zeros((len(x_coord), len(wd), len(ws), len(x_coord), y_ngrid, z_ngrid)) + Ctmp = np.zeros((farm.n_turbines,) + np.shape(flow_field.u_initial_sorted)) # Calculate the velocity deficit sequentially from upstream to downstream turbines for i in range(grid.n_turbines): @@ -334,9 +352,18 @@ def solve(self, *, full_flow: bool = False, grid: TurbineGrid | FlowFieldGrid = # Get the current turbine quantities x_i, y_i, z_i = _get_mean_grid(grid, i) - mask2 = np.array(grid.x_sorted < x_i + 0.01) * np.array(grid.x_sorted > x_i - 0.01) * np.array(grid.y_sorted < y_i + 0.51*126.0) * np.array(grid.y_sorted > y_i - 0.51*126.0) - # mask2 = np.logical_and(np.logical_and(np.logical_and(grid.x_sorted < x_i + 0.01, grid.x_sorted > x_i - 0.01), grid.y_sorted < y_i + 0.51*126.0), grid.y_sorted > y_i - 0.51*126.0) - turb_inflow_field = turb_inflow_field * ~mask2 + (flow_field.u_initial_sorted - turb_u_wake) * mask2 + if not full_flow: + mask = ( + (grid.x_sorted < x_i + 0.01) + * (grid.x_sorted > x_i - 0.01) + * (grid.y_sorted < y_i + 0.51 * 126.0) + * (grid.y_sorted > y_i - 0.51 * 126.0) + ) + turb_inflow_field = turb_inflow_field * ~mask + (flow_field.u_initial_sorted - turb_u_wake) * mask + + turbine_inflow_field = flow_field.u_sorted if full_flow else turb_inflow_field + u_i = turbine_inflow_field[:, :, i:i+1] + v_i = flow_field.v_sorted[:, :, i:i+1] turb_avg_vels = average_velocity(turb_inflow_field) turb_Cts = Ct( @@ -344,19 +371,14 @@ def solve(self, *, full_flow: bool = False, grid: TurbineGrid | FlowFieldGrid = farm.yaw_angles_sorted, farm.turbine_fCts, turbine_type_map=farm.turbine_type_map_sorted, - ) - turb_Cts = turb_Cts[:, :, :, None, None] + )[:, :, :, None, None] turb_aIs = axial_induction( turb_avg_vels, farm.yaw_angles_sorted, farm.turbine_fCts, turbine_type_map=farm.turbine_type_map_sorted, ix_filter=[i], - ) - turb_aIs = turb_aIs[:, :, :, None, None] - - u_i = turb_inflow_field[:, :, i:i+1] - v_i = flow_field.v_sorted[:, :, i:i+1] + )[:, :, :, None, None] axial_induction_i = axial_induction( velocities=flow_field.u_sorted, @@ -364,21 +386,17 @@ def solve(self, *, full_flow: bool = False, grid: TurbineGrid | FlowFieldGrid = fCt=farm.turbine_fCts, turbine_type_map=farm.turbine_type_map_sorted, ix_filter=[i], - ) - - axial_induction_i = axial_induction_i[:, :, :, None, None] + )[:, :, :, 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 + 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 = yaw_angle_i.copy() - if model_manager.enable_secondary_steering: - added_yaw = wake_added_yaw( + if self.model_manager.enable_secondary_steering: + effective_yaw_i += wake_added_yaw( u_i, v_i, flow_field.u_initial_sorted, @@ -391,7 +409,6 @@ def solve(self, *, full_flow: bool = False, grid: TurbineGrid | FlowFieldGrid = axial_induction_i, scale=2.0, ) - effective_yaw_i += added_yaw # Model calculations # NOTE: exponential @@ -422,17 +439,20 @@ def solve(self, *, full_flow: bool = False, grid: TurbineGrid | FlowFieldGrid = scale=2.0 ) - if self.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 = 1.0 - turbine_turbulence_intensity[:, :, i:i+1] = turbulence_intensity_i + gch_gain * I_mixing + if full_flow: + turbine_turbulence_intensity = flow_field.turbulence_intensity_field + else: + if self.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 = 1.0 + turbine_turbulence_intensity[:, :, i:i+1] = turbulence_intensity_i + gch_gain * I_mixing turb_u_wake, Ctmp = self.model_manager.velocity_model.function( i, @@ -450,37 +470,45 @@ def solve(self, *, full_flow: bool = False, grid: TurbineGrid | FlowFieldGrid = **deficit_model_args ) - wake_added_turbulence_intensity = self.model_manager.turbulence_model.function( - ambient_turbulence_intensity, - grid.x_sorted, - x_i, - rotor_diameter_i, - turb_aIs - ) + if not full_flow: + wake_added_turbulence_intensity = self.model_manager.turbulence_model.function( + ambient_turbulence_intensity, + grid.x_sorted, + x_i, + rotor_diameter_i, + turb_aIs + ) - # Calculate wake overlap for wake-added turbulence (WAT) - area_overlap = 1 - np.sum(turb_u_wake <= 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) - * np.array(grid.x_sorted > x_i) - * np.array(np.abs(y_i - grid.y_sorted) < 2 * rotor_diameter_i) - * np.array(grid.x_sorted <= downstream_influence_length + x_i) - ) + # Calculate wake overlap for wake-added turbulence (WAT) + area_overlap = ( + 1 + - np.sum(turb_u_wake <= 0.05, axis=(3, 4)) + / (grid.grid_resolution * grid.grid_resolution) + )[:, :, :, None, None] - # Combine turbine TIs with WAT - turbine_turbulence_intensity = np.maximum( np.sqrt( ti_added ** 2 + ambient_turbulence_intensity ** 2 ) , turbine_turbulence_intensity ) + # 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.v_sorted += v_wake flow_field.w_sorted += w_wake - flow_field.u_sorted = turb_inflow_field + + flow_field.u_sorted = flow_field.u_initial_sorted - turb_u_wake if full_flow else turb_inflow_field - flow_field.turbulence_intensity_field = np.mean(turbine_turbulence_intensity, axis=(3,4)) - flow_field.turbulence_intensity_field = flow_field.turbulence_intensity_field[:,:,:,None,None] + if not full_flow: + flow_field.turbulence_intensity_field = np.mean(turbine_turbulence_intensity, axis=(3, 4))[:, :, :, None, None] From 35259098fedfd97326b1d45098fec941e9fdc599 Mon Sep 17 00:00:00 2001 From: Rob Hammond <13874373+RHammond2@users.noreply.github.com> Date: Tue, 20 Sep 2022 11:47:18 -0600 Subject: [PATCH 04/20] generalize long computations that naturally break onto multiple lines --- floris/simulation/solver.py | 82 +++++++++++++++++++++++-------------- 1 file changed, 51 insertions(+), 31 deletions(-) diff --git a/floris/simulation/solver.py b/floris/simulation/solver.py index 646370f36..da5f8134c 100644 --- a/floris/simulation/solver.py +++ b/floris/simulation/solver.py @@ -31,21 +31,31 @@ from floris.type_dec import NDArrayFloat -def _get_mean_grid(grid: TurbineGrid | FlowFieldGrid, i: int) -> tuple[NDArrayFloat, NDArrayFloat, NDArrayFloat]: - """Calculates the mean of the `grid`'s `x_sorted`, `y_sorted`, and `z_sorted` attributes and - returns them as a new 5-dimensional object to align with expected internal structures. +def _expansion_mean(x: NDArrayFloat) -> NDArrayFloat: + """Calculates the mean of the `x` over axis (3, 4) and returns the result as a new 5-dimensional + object to align with expected internal structures. Args: - grid (TurbineGrid | FlowFieldGrid): The solver's grid object. - i (int): The turbine index that the mean x, y, and z will be computed over. + x (NDArrayFloat): An array to calculate the mean. Returns: - tuple[NDArrayFloat, NDArrayFloat, NDArrayFloat]: The mean over axis 3 and 4 at turbine i. + NDArrayFloat: The mean over axis 3 and 4 at turbine i and expandd to 5 dimensions. """ - x_i = np.mean(grid.x_sorted[:, :, i:i+1], axis=(3, 4))[:, :, :, None, None] - y_i = np.mean(grid.y_sorted[:, :, i:i+1], axis=(3, 4))[:, :, :, None, None] - z_i = np.mean(grid.z_sorted[:, :, i:i+1], axis=(3, 4))[:, :, :, None, None] - return x_i, y_i, z_i + return np.mean(x, axis=(3, 4))[:, :, :, None, None] + + +def _expansion_mean_i(x: NDArrayFloat, i: int) -> NDArrayFloat: + """Calculates the mean of the `x` over axis (3, 4) at turbine index `i` and returns the result + as a new 5-dimensional object to align with expected internal structures. + + Args: + x (NDArrayFloat): An array to calculate the mean. + i (int): The turbine index for where to compute the mean. + + Returns: + NDArrayFloat: The mean over axis 3 and 4 at turbine i and expandd to 5 dimensions. + """ + return _expansion_mean(x[:, :, i:i+1]) def calculate_area_overlap(wake_velocities, freestream_velocities, y_ngrid, z_ngrid): @@ -173,7 +183,9 @@ def solve( # Calculate the velocity deficit sequentially from upstream to downstream turbines for i in range(grid.n_turbines): # Get the current turbine quantities - x_i, y_i, z_i = _get_mean_grid(grid, i) + x_i = _expansion_mean_i(grid.x_sorted, i) + y_i = _expansion_mean_i(grid.y_sorted, i) + z_i = _expansion_mean_i(grid.z_sorted, i) u_i = flow_field.u_sorted[:, :, i:i+1] v_i = flow_field.v_sorted[:, :, i:i+1] @@ -304,8 +316,7 @@ def solve( flow_field.w_sorted += w_wake if not full_flow: - flow_field.turbulence_intensity_field = np.mean(turbine_turbulence_intensity, axis=(3, 4)) - flow_field.turbulence_intensity_field = flow_field.turbulence_intensity_field[:, :, :, None, None] + flow_field.turbulence_intensity_field = _expansion_mean(turbine_turbulence_intensity) @define(auto_attribs=True) @@ -320,10 +331,10 @@ def solve( *, full_flow: bool = False, farm: Farm, - flow_field: FlowFieldGrid = None, + flow_field: FlowField = None, grid: TurbineGrid | FlowFieldGrid = None ) -> None: - + if farm is None: farm = self.farm if flow_field is None: @@ -331,6 +342,9 @@ def solve( if grid is None: grid = self.grid + gch_gain = 1.0 + scale_factor = 2.0 + # <> deflection_model_args = self.model_manager.deflection_model.prepare_function(grid, flow_field) deficit_model_args = self.model_manager.velocity_model.prepare_function(grid, flow_field) @@ -339,9 +353,13 @@ def solve( v_wake = np.zeros_like(flow_field.v_initial_sorted) w_wake = np.zeros_like(flow_field.w_initial_sorted) turb_u_wake = np.zeros_like(flow_field.u_initial_sorted) - turb_inflow_field = copy.deepcopy(flow_field.u_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)) + # Not needed for full flow solve, but isn't necessary to have in an if statement + turb_inflow_field = copy.deepcopy(flow_field.u_initial_sorted) + turbine_turbulence_intensity = np.full( + (flow_field.n_wind_directions, flow_field.n_wind_speeds, farm.n_turbines, 1, 1), + flow_field.turbulence_intensity + ) ambient_turbulence_intensity = flow_field.turbulence_intensity Ctmp = np.zeros((farm.n_turbines,) + np.shape(flow_field.u_initial_sorted)) @@ -351,6 +369,8 @@ def solve( # Get the current turbine quantities x_i, y_i, z_i = _get_mean_grid(grid, i) + u_i = turbine_inflow_field[:, :, i:i+1] + v_i = flow_field.v_sorted[:, :, i:i+1] if not full_flow: mask = ( @@ -359,11 +379,9 @@ def solve( * (grid.y_sorted < y_i + 0.51 * 126.0) * (grid.y_sorted > y_i - 0.51 * 126.0) ) - turb_inflow_field = turb_inflow_field * ~mask + (flow_field.u_initial_sorted - turb_u_wake) * mask + turb_inflow_field *= ~mask + (flow_field.u_initial_sorted - turb_u_wake) * mask turbine_inflow_field = flow_field.u_sorted if full_flow else turb_inflow_field - u_i = turbine_inflow_field[:, :, i:i+1] - v_i = flow_field.v_sorted[:, :, i:i+1] turb_avg_vels = average_velocity(turb_inflow_field) turb_Cts = Ct( @@ -372,13 +390,15 @@ def solve( farm.turbine_fCts, turbine_type_map=farm.turbine_type_map_sorted, )[:, :, :, None, None] - turb_aIs = axial_induction( - turb_avg_vels, - farm.yaw_angles_sorted, - farm.turbine_fCts, - turbine_type_map=farm.turbine_type_map_sorted, - ix_filter=[i], - )[:, :, :, None, None] + + if not full_flow: + turb_aIs = axial_induction( + turb_avg_vels, + farm.yaw_angles_sorted, + farm.turbine_fCts, + turbine_type_map=farm.turbine_type_map_sorted, + ix_filter=[i], + )[:, :, :, None, None] axial_induction_i = axial_induction( velocities=flow_field.u_sorted, @@ -407,7 +427,7 @@ def solve( turb_Cts[:, :, i:i+1], TSR_i, axial_induction_i, - scale=2.0, + scale=scale_factor, ) # Model calculations @@ -436,7 +456,7 @@ def solve( turb_Cts[:, :, i:i+1], TSR_i, axial_induction_i, - scale=2.0 + scale=scale_factor, ) if full_flow: @@ -451,9 +471,9 @@ def solve( v_wake[:, :, i:i+1], w_wake[:, :, i:i+1], ) - gch_gain = 1.0 turbine_turbulence_intensity[:, :, i:i+1] = turbulence_intensity_i + gch_gain * I_mixing + # NOTE: exponential turb_u_wake, Ctmp = self.model_manager.velocity_model.function( i, x_i, @@ -508,7 +528,7 @@ def solve( flow_field.u_sorted = flow_field.u_initial_sorted - turb_u_wake if full_flow else turb_inflow_field if not full_flow: - flow_field.turbulence_intensity_field = np.mean(turbine_turbulence_intensity, axis=(3, 4))[:, :, :, None, None] + flow_field.turbulence_intensity_field = _expansion_mean(turbine_turbulence_intensity) From bb28fdc340bbf4a4c39d1cf7c5fec14248eea84f Mon Sep 17 00:00:00 2001 From: Rob Hammond <13874373+RHammond2@users.noreply.github.com> Date: Tue, 20 Sep 2022 14:45:42 -0600 Subject: [PATCH 05/20] add turbopark refactor --- floris/simulation/solver.py | 312 +++++++++++++++++++++++++++++++++--- 1 file changed, 289 insertions(+), 23 deletions(-) diff --git a/floris/simulation/solver.py b/floris/simulation/solver.py index da5f8134c..6dfe43b8d 100644 --- a/floris/simulation/solver.py +++ b/floris/simulation/solver.py @@ -85,10 +85,11 @@ def solve( self, *, full_flow: bool = False, - farm: Farm, + farm: Farm = None, flow_field: FlowFieldGrid = None, grid: TurbineGrid | FlowFieldGrid = None ) -> None: + # TODO: Update with the new logger functionality that is on its way pass def full_flow_solve(self): @@ -100,8 +101,7 @@ def full_flow_solve(self): """ if not isinstance(self.grid, FlowFieldGrid): raise TypeError("Cannot run `full_flow_solve` with a `TurbineGrid` object, reinitialize with the input to `grid` being a `FlowFieldGrid` object.") - - + # TODO: Why are we making copies of the originals, and then never using the original, # can we just use the originally provided values and get on with the calculating? farm = copy.deepcopy(self.farm) @@ -156,9 +156,16 @@ def solve( Args: full_flow (bool, optional): Runs the full flow solver when True, and the standard sequential solver, when False. Defaults to False. + farm (Farm, optional): Allows for a non-initialized `farm` object to be used. It should + be noted that this functionality is intended for use with `full_flow_solve`. + Defaults to None. + flow_field (FlowField, optional): Allows for a non-initialized `flow_field` object to be + used. It should be noted that this functionality is intended for use with + `full_flow_solve`.Defaults to None. grid (TurbineGrid | FlowFieldGrid, optional): Allows for a non-initialized `grid` object - to be used. It should be noted that this functionality is intended for use with - `full_flow_solve`, which computes over a `TurbineGrid`, then a `FlowFieldGrid`. Defaults to None. + to be used. It should be noted that this functionality is intended for use with + `full_flow_solve`, which computes over a `TurbineGrid`, then a `FlowFieldGrid`. + Defaults to None. """ if farm is None: farm = self.farm @@ -330,11 +337,26 @@ def solve( self, *, full_flow: bool = False, - farm: Farm, + farm: Farm = None, flow_field: FlowField = None, grid: TurbineGrid | FlowFieldGrid = None ) -> None: + """Runs the CC sover methodology, or full flow CC solver methodology for a wind farm. + Args: + full_flow (bool, optional): Runs the full flow solver when True, and the standard CC + solver, when False. Defaults to False. + farm (Farm, optional): Allows for a non-initialized `farm` object to be used. It should + be noted that this functionality is intended for use with `full_flow_solve`. + Defaults to None. + flow_field (FlowField, optional): Allows for a non-initialized `flow_field` object to be + used. It should be noted that this functionality is intended for use with + `full_flow_solve`.Defaults to None. + grid (TurbineGrid | FlowFieldGrid, optional): Allows for a non-initialized `grid` object + to be used. It should be noted that this functionality is intended for use with + `full_flow_solve`, which computes over a `TurbineGrid`, then a `FlowFieldGrid`. + Defaults to None. + """ if farm is None: farm = self.farm if flow_field is None: @@ -368,7 +390,9 @@ def solve( for i in range(grid.n_turbines): # Get the current turbine quantities - x_i, y_i, z_i = _get_mean_grid(grid, i) + x_i = _expansion_mean_i(grid.x_sorted, i) + y_i = _expansion_mean_i(grid.y_sorted, i) + z_i = _expansion_mean_i(grid.z_sorted, i) u_i = turbine_inflow_field[:, :, i:i+1] v_i = flow_field.v_sorted[:, :, i:i+1] @@ -531,6 +555,248 @@ def solve( flow_field.turbulence_intensity_field = _expansion_mean(turbine_turbulence_intensity) +@define(auto_attribs=True) +class TurbOParkSolver(Solver): + farm: Farm = field(converter=copy.deepcopy, validator=attrs.validators.instance_of(Farm)) + flow_field: FlowField = field(converter=copy.deepcopy, validator=attrs.validators.instance_of(Farm)) + grid: TurbineGrid | FlowFieldGrid = field(converter=copy.deepcopy, validator=attrs.validators.instance_of((FlowFieldGrid, TurbineGrid))) + model_manager: WakeModelManager = field(converter=copy.deepcopy, validator=attrs.validators.instance_of(WakeModelManager)) + + def full_flow_solve(self): + # TODO: Update with the new logger functionality that is on its way + raise NotImplementedError("TurbOParkSolver has no full flow field solving capabilities yet.") + + def solve( + self, + *, + full_flow: bool = False, + farm: Farm = None, + flow_field: FlowField = None, + grid: TurbineGrid | FlowFieldGrid = None + ) -> None: + """Runs the TurbOPark sover methodology, or full flow TurbOPark solver methodology for a + wind farm. + + Args: + full_flow (bool, optional): Runs the full flow solver when True, and the standard + TurbOPark solver, when False. Defaults to False. + farm (Farm, optional): Allows for a non-initialized `farm` object to be used. It should + be noted that this functionality is intended for use with `full_flow_solve`. + Defaults to None. + flow_field (FlowField, optional): Allows for a non-initialized `flow_field` object to be + used. It should be noted that this functionality is intended for use with + `full_flow_solve`.Defaults to None. + grid (TurbineGrid | FlowFieldGrid, optional): Allows for a non-initialized `grid` object + to be used. It should be noted that this functionality is intended for use with + `full_flow_solve`, which computes over a `TurbineGrid`, then a `FlowFieldGrid`. + Defaults to 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. + + if farm is None: + farm = self.farm + if flow_field is None: + flow_field = self.flow_field + if grid is None: + grid = self.grid + + gch_gain = 2 + + # <> + deflection_model_args = self.model_manager.deflection_model.prepare_function(grid, flow_field) + deficit_model_args = self.model_manager.velocity_model.prepare_function(grid, flow_field) + + # 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) + shape = (farm.n_turbines,) + np.shape(flow_field.u_initial_sorted) + velocity_deficit = np.zeros(shape) + deflection_field = np.zeros_like(flow_field.u_initial_sorted) + + turbine_turbulence_intensity = np.full( + (flow_field.n_wind_directions, flow_field.n_wind_speeds, farm.n_turbines, 1, 1), + flow_field.turbulence_intensity + ) + 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 = _expansion_mean_i(grid.x_sorted, i) + y_i = _expansion_mean_i(grid.y_sorted, i) + z_i = _expansion_mean_i(grid.z_sorted, i) + u_i = flow_field.u_sorted[:, :, i:i+1] + v_i = flow_field.v_sorted[:, :, i:i+1] + + Cts = Ct( + velocities=flow_field.u_sorted, + yaw_angle=farm.yaw_angles_sorted, + fCt=farm.turbine_fCts, + turbine_type_map=farm.turbine_type_map_sorted, + ) + + # Since we are filtering for the ith turbine in the Ct function, get the first index here (0:1) + ct_i = Ct( + velocities=flow_field.u_sorted, + yaw_angle=farm.yaw_angles_sorted, + fCt=farm.turbine_fCts, + turbine_type_map=farm.turbine_type_map_sorted, + ix_filter=[i], + )[:, :, 0:1, None, None] + + # Since we are filtering for the ith turbine in the axial induction function, get the first index here (0:1) + axial_induction_i = axial_induction( + velocities=flow_field.u_sorted, + yaw_angle=farm.yaw_angles_sorted, + fCt=farm.turbine_fCts, + turbine_type_map=farm.turbine_type_map_sorted, + ix_filter=[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 self.model_manager.enable_secondary_steering: + effective_yaw_i += 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 + ) + + # Model calculations + # NOTE: exponential + if not np.all(farm.yaw_angles_sorted): + self.model_manager.deflection_model.logger.warning("WARNING: Deflection with the TurbOPark model has not been fully validated. This is an initial implementation, and we advise you use at your own risk and perform a thorough examination of the results.") + for ii in range(i): + x_ii = _expansion_mean_i(grid.x_sorted, ii) + y_ii = _expansion_mean_i(grid.y_sorted, ii) + + yaw_ii = farm.yaw_angles_sorted[:, :, ii:ii+1, None, None] + turbulence_intensity_ii = turbine_turbulence_intensity[:, :, ii:ii+1] + ct_ii = Ct( + velocities=flow_field.u_sorted, + yaw_angle=farm.yaw_angles_sorted, + fCt=farm.turbine_fCts, + turbine_type_map=farm.turbine_type_map_sorted, + ix_filter=[ii] + )[:, :, 0:1, None, None] + rotor_diameter_ii = farm.rotor_diameters_sorted[:, :, ii:ii+1, None, None] + + deflection_field_ii = self.model_manager.deflection_model.function( + x_ii, + y_ii, + yaw_ii, + turbulence_intensity_ii, + ct_ii, + rotor_diameter_ii, + **deflection_model_args + ) + + deflection_field[:, :, ii:ii+1, :, :] = deflection_field_ii[:, :, i:i+1, :, :] + + if self.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 + ) + + if self.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], + ) + turbine_turbulence_intensity[:, :, i:i+1] = turbulence_intensity_i + gch_gain * I_mixing + + # NOTE: exponential + velocity_deficit = self.model_manager.velocity_model.function( + x_i, + y_i, + z_i, + turbine_turbulence_intensity, + Cts[:, :, :, None, None], + rotor_diameter_i, + farm.rotor_diameters_sorted[:, :, :, None, None], + i, + deflection_field, + **deficit_model_args + ) + + wake_field = self.model_manager.combination_model.function( + wake_field, + velocity_deficit * flow_field.u_initial_sorted + ) + + wake_added_turbulence_intensity = self.model_manager.turbulence_model.function( + ambient_turbulence_intensity, + grid.x_sorted, + x_i, + rotor_diameter_i, + axial_induction_i + ) + + # TODO: leaving this in for GCH quantities; will need to find another way to compute area_overlap + # as the current wake deficit is solved for only upstream turbines; could use WAT_upstream + # 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 ** 2) + )[:, :, :, 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 = _expansion_mean(turbine_turbulence_intensity) + def sequential_solver(farm: Farm, flow_field: FlowField, grid: TurbineGrid, model_manager: WakeModelManager) -> None: # Algorithm @@ -578,9 +844,9 @@ def sequential_solver(farm: Farm, flow_field: FlowField, grid: TurbineGrid, mode axial_induction_i = axial_induction_i[:, :, 0:1, None, None] # Since we are filtering for the i'th turbine in the axial induction function, get the first index here (0:1) 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] + 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 @@ -763,9 +1029,9 @@ def full_flow_sequential_solver(farm: Farm, flow_field: FlowField, flow_field_gr axial_induction_i = axial_induction_i[:, :, 0:1, None, None] # Since we are filtering for the i'th turbine in the axial induction function, get the first index here (0:1) turbulence_intensity_i = turbine_grid_flow_field.turbulence_intensity_field[:, :, i:i+1] yaw_angle_i = turbine_grid_farm.yaw_angles_sorted[:, :, i:i+1, None, None] - hub_height_i = turbine_grid_farm.hub_heights_sorted[: ,:, i:i+1, None, None] - rotor_diameter_i = turbine_grid_farm.rotor_diameters_sorted[: ,:, i:i+1, None, None] - TSR_i = turbine_grid_farm.TSRs_sorted[: ,:, i:i+1, None, None] + hub_height_i = turbine_grid_farm.hub_heights_sorted[:, :, i:i+1, None, None] + rotor_diameter_i = turbine_grid_farm.rotor_diameters_sorted[:, :, i:i+1, None, None] + TSR_i = turbine_grid_farm.TSRs_sorted[:, :, i:i+1, None, None] effective_yaw_i = np.zeros_like(yaw_angle_i) effective_yaw_i += yaw_angle_i @@ -902,9 +1168,9 @@ def cc_solver(farm: Farm, flow_field: FlowField, grid: TurbineGrid, model_manage 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] + 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 @@ -1088,9 +1354,9 @@ def full_flow_cc_solver(farm: Farm, flow_field: FlowField, flow_field_grid: Flow turbulence_intensity_i = turbine_grid_flow_field.turbulence_intensity_field[:, :, i:i+1] yaw_angle_i = turbine_grid_farm.yaw_angles_sorted[:, :, i:i+1, None, None] - hub_height_i = turbine_grid_farm.hub_heights_sorted[: ,:, i:i+1, None, None] - rotor_diameter_i = turbine_grid_farm.rotor_diameters_sorted[: ,:, i:i+1, None, None] - TSR_i = turbine_grid_farm.TSRs_sorted[: ,:, i:i+1, None, None] + hub_height_i = turbine_grid_farm.hub_heights_sorted[:, :, i:i+1, None, None] + rotor_diameter_i = turbine_grid_farm.rotor_diameters_sorted[:, :, i:i+1, None, None] + TSR_i = turbine_grid_farm.TSRs_sorted[:, :, i:i+1, None, None] effective_yaw_i = np.zeros_like(yaw_angle_i) effective_yaw_i += yaw_angle_i @@ -1216,9 +1482,9 @@ def turbopark_solver(farm: Farm, flow_field: FlowField, grid: TurbineGrid, model axial_induction_i = axial_induction_i[:, :, 0:1, None, None] # Since we are filtering for the i'th turbine in the axial induction function, get the first index here (0:1) 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] + 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 @@ -1258,7 +1524,7 @@ def turbopark_solver(farm: Farm, flow_field: FlowField, grid: TurbineGrid, model ix_filter=[ii] ) ct_ii = ct_ii[:, :, 0:1, None, None] - rotor_diameter_ii = farm.rotor_diameters_sorted[: ,:, ii:ii+1, None, None] + rotor_diameter_ii = farm.rotor_diameters_sorted[:, :, ii:ii+1, None, None] deflection_field_ii = model_manager.deflection_model.function( x_ii, From e03f3be49985d8a4ec9d90cd9f80db6b67b8c649 Mon Sep 17 00:00:00 2001 From: Rob Hammond <13874373+RHammond2@users.noreply.github.com> Date: Tue, 20 Sep 2022 14:59:20 -0600 Subject: [PATCH 06/20] improve spacing and run linter --- floris/simulation/solver.py | 137 ++++++++++++++++++++++-------------- 1 file changed, 84 insertions(+), 53 deletions(-) diff --git a/floris/simulation/solver.py b/floris/simulation/solver.py index 6dfe43b8d..fcc1b1dfd 100644 --- a/floris/simulation/solver.py +++ b/floris/simulation/solver.py @@ -10,25 +10,29 @@ # License for the specific language governing permissions and limitations under # the License. -from abc import abstractmethod import copy +from abc import abstractmethod import attrs import numpy as np -from attrs import define, field +from attrs import field, define -from floris.simulation import Farm, flow_field -from floris.simulation import TurbineGrid, FlowFieldGrid -from floris.simulation import Ct, axial_induction -from floris.simulation import FlowField -from floris.simulation.turbine import average_velocity +from floris.type_dec import NDArrayFloat +from floris.simulation import ( + Ct, + Farm, + FlowField, + TurbineGrid, + FlowFieldGrid, + axial_induction, +) from floris.simulation.wake import WakeModelManager +from floris.simulation.turbine import average_velocity from floris.simulation.wake_deflection.gauss import ( - calculate_transverse_velocity, wake_added_yaw, - yaw_added_turbulence_mixing + yaw_added_turbulence_mixing, + calculate_transverse_velocity, ) -from floris.type_dec import NDArrayFloat def _expansion_mean(x: NDArrayFloat) -> NDArrayFloat: @@ -77,8 +81,13 @@ def calculate_area_overlap(wake_velocities, freestream_velocities, y_ngrid, z_ng class Solver: farm: Farm = field(converter=copy.deepcopy, validator=attrs.validators.instance_of(Farm)) flow_field: FlowField = field(converter=copy.deepcopy, validator=attrs.validators.instance_of(Farm)) - grid: TurbineGrid | FlowFieldGrid = field(converter=copy.deepcopy, validator=attrs.validators.instance_of((FlowFieldGrid, TurbineGrid))) - model_manager: WakeModelManager = field(converter=copy.deepcopy, validator=attrs.validators.instance_of(WakeModelManager)) + grid: TurbineGrid | FlowFieldGrid = field( + converter=copy.deepcopy, + validator=attrs.validators.instance_of((FlowFieldGrid, TurbineGrid)) + ) + model_manager: WakeModelManager = field( + converter=copy.deepcopy, validator=attrs.validators.instance_of(WakeModelManager) + ) @abstractmethod def solve( @@ -100,7 +109,7 @@ def full_flow_solve(self): TypeError: Raised if initialized `grid` value is not a `FlowFieldGrid` object. """ if not isinstance(self.grid, FlowFieldGrid): - raise TypeError("Cannot run `full_flow_solve` with a `TurbineGrid` object, reinitialize with the input to `grid` being a `FlowFieldGrid` object.") + raise TypeError("Cannot run `full_flow_solve` with a `TurbineGrid` object, reinitialize with the input to `grid` being a `FlowFieldGrid` object.") # noqa: E501 # TODO: Why are we making copies of the originals, and then never using the original, # can we just use the originally provided values and get on with the calculating? @@ -137,10 +146,6 @@ def full_flow_solve(self): @define(auto_attribs=True) class SequentialSolver(Solver): - farm: Farm = field(converter=copy.deepcopy, validator=attrs.validators.instance_of(Farm)) - flow_field: FlowField = field(converter=copy.deepcopy, validator=attrs.validators.instance_of(Farm)) - grid: TurbineGrid | FlowFieldGrid = field(converter=copy.deepcopy, validator=attrs.validators.instance_of((FlowFieldGrid, TurbineGrid))) - model_manager: WakeModelManager = field(converter=copy.deepcopy, validator=attrs.validators.instance_of(WakeModelManager)) def solve( self, @@ -184,7 +189,10 @@ def solve( w_wake = np.zeros_like(flow_field.w_initial_sorted) if not full_flow: - turbine_turbulence_intensity = flow_field.turbulence_intensity * np.ones((flow_field.n_wind_directions, flow_field.n_wind_speeds, farm.n_turbines, 1, 1)) + turbine_turbulence_intensity = np.full( + (flow_field.n_wind_directions, flow_field.n_wind_speeds, farm.n_turbines, 1, 1), + flow_field.turbulence_intensity + ) ambient_turbulence_intensity = flow_field.turbulence_intensity # Calculate the velocity deficit sequentially from upstream to downstream turbines @@ -205,13 +213,14 @@ def solve( ix_filter=[i], )[:, :, 0:1, None, None] + # Since we are filtering for the ith turbine in the axial induction function, get the first index here (0:1) axial_induction_i = axial_induction( velocities=flow_field.u_sorted, yaw_angle=farm.yaw_angles_sorted, fCt=farm.turbine_fCts, turbine_type_map=farm.turbine_type_map_sorted, ix_filter=[i], - )[:, :, 0:1, None, None] # Since we are filtering for the i'th turbine in the axial induction function, get the first index here (0:1) + )[:, :, 0:1, None, None] turbulence_intensity_i = flow_field.turbulence_intensity_field[:, :, 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] @@ -302,8 +311,10 @@ def solve( ) # 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] + area_overlap = ( + np.sum(velocity_deficit * flow_field.u_initial_sorted > 0.05, axis=(3, 4)) + / (grid.grid_resolution * grid.grid_resolution) + )[:, :, :, None, None] # Modify wake added turbulence by wake area overlap downstream_influence_length = 15 * rotor_diameter_i @@ -316,7 +327,10 @@ def solve( ) # Combine turbine TIs with WAT - turbine_turbulence_intensity = np.maximum( np.sqrt( ti_added ** 2 + ambient_turbulence_intensity ** 2 ) , turbine_turbulence_intensity ) + 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 @@ -328,10 +342,6 @@ def solve( @define(auto_attribs=True) class CCSolver(Solver): - farm: Farm = field(converter=copy.deepcopy, validator=attrs.validators.instance_of(Farm)) - flow_field: FlowField = field(converter=copy.deepcopy, validator=attrs.validators.instance_of(Farm)) - grid: TurbineGrid | FlowFieldGrid = field(converter=copy.deepcopy, validator=attrs.validators.instance_of((FlowFieldGrid, TurbineGrid))) - model_manager: WakeModelManager = field(converter=copy.deepcopy, validator=attrs.validators.instance_of(WakeModelManager)) def solve( self, @@ -377,7 +387,7 @@ def solve( turb_u_wake = np.zeros_like(flow_field.u_initial_sorted) # Not needed for full flow solve, but isn't necessary to have in an if statement - turb_inflow_field = copy.deepcopy(flow_field.u_initial_sorted) + turbine_inflow_field = copy.deepcopy(flow_field.u_initial_sorted) turbine_turbulence_intensity = np.full( (flow_field.n_wind_directions, flow_field.n_wind_speeds, farm.n_turbines, 1, 1), flow_field.turbulence_intensity @@ -403,18 +413,18 @@ def solve( * (grid.y_sorted < y_i + 0.51 * 126.0) * (grid.y_sorted > y_i - 0.51 * 126.0) ) - turb_inflow_field *= ~mask + (flow_field.u_initial_sorted - turb_u_wake) * mask + turbine_inflow_field *= ~mask + (flow_field.u_initial_sorted - turb_u_wake) * mask - turbine_inflow_field = flow_field.u_sorted if full_flow else turb_inflow_field + turbine_inflow_field = flow_field.u_sorted if full_flow else turbine_inflow_field - turb_avg_vels = average_velocity(turb_inflow_field) + turb_avg_vels = average_velocity(turbine_inflow_field) turb_Cts = Ct( turb_avg_vels, farm.yaw_angles_sorted, farm.turbine_fCts, turbine_type_map=farm.turbine_type_map_sorted, )[:, :, :, None, None] - + if not full_flow: turb_aIs = axial_induction( turb_avg_vels, @@ -548,8 +558,8 @@ def solve( flow_field.v_sorted += v_wake flow_field.w_sorted += w_wake - - flow_field.u_sorted = flow_field.u_initial_sorted - turb_u_wake if full_flow else turb_inflow_field + + flow_field.u_sorted = flow_field.u_initial_sorted - turb_u_wake if full_flow else turbine_inflow_field if not full_flow: flow_field.turbulence_intensity_field = _expansion_mean(turbine_turbulence_intensity) @@ -557,10 +567,6 @@ def solve( @define(auto_attribs=True) class TurbOParkSolver(Solver): - farm: Farm = field(converter=copy.deepcopy, validator=attrs.validators.instance_of(Farm)) - flow_field: FlowField = field(converter=copy.deepcopy, validator=attrs.validators.instance_of(Farm)) - grid: TurbineGrid | FlowFieldGrid = field(converter=copy.deepcopy, validator=attrs.validators.instance_of((FlowFieldGrid, TurbineGrid))) - model_manager: WakeModelManager = field(converter=copy.deepcopy, validator=attrs.validators.instance_of(WakeModelManager)) def full_flow_solve(self): # TODO: Update with the new logger functionality that is on its way @@ -596,7 +602,7 @@ def solve( # 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. - + if farm is None: farm = self.farm if flow_field is None: @@ -626,7 +632,7 @@ def solve( # Calculate the velocity deficit sequentially from upstream to downstream turbines for i in range(grid.n_turbines): - + # Get the current turbine quantities x_i = _expansion_mean_i(grid.x_sorted, i) y_i = _expansion_mean_i(grid.y_sorted, i) @@ -649,7 +655,7 @@ def solve( turbine_type_map=farm.turbine_type_map_sorted, ix_filter=[i], )[:, :, 0:1, None, None] - + # Since we are filtering for the ith turbine in the axial induction function, get the first index here (0:1) axial_induction_i = axial_induction( velocities=flow_field.u_sorted, @@ -685,11 +691,11 @@ def solve( # Model calculations # NOTE: exponential if not np.all(farm.yaw_angles_sorted): - self.model_manager.deflection_model.logger.warning("WARNING: Deflection with the TurbOPark model has not been fully validated. This is an initial implementation, and we advise you use at your own risk and perform a thorough examination of the results.") + self.model_manager.deflection_model.logger.warning("WARNING: Deflection with the TurbOPark model has not been fully validated. This is an initial implementation, and we advise you use at your own risk and perform a thorough examination of the results.") # noqa: #501 for ii in range(i): x_ii = _expansion_mean_i(grid.x_sorted, ii) y_ii = _expansion_mean_i(grid.y_sorted, ii) - + yaw_ii = farm.yaw_angles_sorted[:, :, ii:ii+1, None, None] turbulence_intensity_ii = turbine_turbulence_intensity[:, :, ii:ii+1] ct_ii = Ct( @@ -797,7 +803,7 @@ def solve( flow_field.turbulence_intensity_field = _expansion_mean(turbine_turbulence_intensity) - +# flake8: noqa def sequential_solver(farm: Farm, flow_field: FlowField, grid: TurbineGrid, model_manager: WakeModelManager) -> None: # Algorithm # For each turbine, calculate its effect on every downstream turbine. @@ -821,7 +827,12 @@ def sequential_solver(farm: Farm, flow_field: FlowField, grid: TurbineGrid, mode for i in range(grid.n_turbines): # Get the current turbine quantities - x_i, y_i, z_i = _get_mean_grid(grid, i) + 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] @@ -1006,7 +1017,12 @@ def full_flow_sequential_solver(farm: Farm, flow_field: FlowField, flow_field_gr for i in range(flow_field_grid.n_turbines): # Get the current turbine quantities - x_i, y_i, z_i = _get_mean_grid(turbine_grid, i) + 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 = turbine_grid_flow_field.u_sorted[:, :, i:i+1] v_i = turbine_grid_flow_field.v_sorted[:, :, i:i+1] @@ -1122,7 +1138,7 @@ def cc_solver(farm: Farm, flow_field: FlowField, grid: TurbineGrid, model_manage shape = (farm.n_turbines,) + np.shape(flow_field.u_initial_sorted) Ctmp = np.zeros((shape)) # Ctmp = np.zeros((len(x_coord), len(wd), len(ws), len(x_coord), y_ngrid, z_ngrid)) - + sigma_i = np.zeros((shape)) # sigma_i = np.zeros((len(x_coord), len(wd), len(ws), len(x_coord), y_ngrid, z_ngrid)) @@ -1130,7 +1146,12 @@ def cc_solver(farm: Farm, flow_field: FlowField, grid: TurbineGrid, model_manage for i in range(grid.n_turbines): # Get the current turbine quantities - x_i, y_i, z_i = _get_mean_grid(grid, i) + 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] mask2 = np.array(grid.x_sorted < x_i + 0.01) * np.array(grid.x_sorted > x_i - 0.01) * np.array(grid.y_sorted < y_i + 0.51*126.0) * np.array(grid.y_sorted > y_i - 0.51*126.0) # mask2 = np.logical_and(np.logical_and(np.logical_and(grid.x_sorted < x_i + 0.01, grid.x_sorted > x_i - 0.01), grid.y_sorted < y_i + 0.51*126.0), grid.y_sorted > y_i - 0.51*126.0) @@ -1143,7 +1164,7 @@ def cc_solver(farm: Farm, flow_field: FlowField, grid: TurbineGrid, model_manage farm.turbine_fCts, turbine_type_map=farm.turbine_type_map_sorted, ) - turb_Cts = turb_Cts[:, :, :, None, None] + turb_Cts = turb_Cts[:, :, :, None, None] turb_aIs = axial_induction( turb_avg_vels, farm.yaw_angles_sorted, @@ -1329,7 +1350,12 @@ def full_flow_cc_solver(farm: Farm, flow_field: FlowField, flow_field_grid: Flow for i in range(flow_field_grid.n_turbines): # Get the current turbine quantities - x_i, y_i, z_i = _get_mean_grid(turbine_grid, i) + 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 = turbine_grid_flow_field.u_sorted[:, :, i:i+1] v_i = turbine_grid_flow_field.v_sorted[:, :, i:i+1] @@ -1452,7 +1478,12 @@ def turbopark_solver(farm: Farm, flow_field: FlowField, grid: TurbineGrid, model # Calculate the velocity deficit sequentially from upstream to downstream turbines for i in range(grid.n_turbines): # Get the current turbine quantities - x_i, y_i, z_i = _get_mean_grid(grid, i) + 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] @@ -1511,7 +1542,7 @@ def turbopark_solver(farm: Farm, flow_field: FlowField, grid: TurbineGrid, model for ii in range(i): x_ii = np.mean(grid.x_sorted[:, :, ii:ii+1], axis=(3, 4)) x_ii = x_ii[:, :, :, None, None] - y_ii = np.mean(grid.y_sorted[:, :, ii:ii+1], axis=(3, 4)) + y_ii = np.mean(grid.y_sorted[:, :, ii:ii+1], axis=(3, 4)) y_ii = y_ii[:, :, :, None, None] yaw_ii = farm.yaw_angles_sorted[:, :, ii:ii+1, None, None] @@ -1625,7 +1656,7 @@ def full_flow_turbopark_solver(farm: Farm, flow_field: FlowField, flow_field_gri # TODO: Below is a first attempt at plotting, and uses just the values on the rotor. The current TurbOPark model requires that # points to be calculated are only at turbine locations. Modification will be required to allow for full flow field calculations. - + # # Get the flow quantities and turbine performance # turbine_grid_farm = copy.deepcopy(farm) # turbine_grid_flow_field = copy.deepcopy(flow_field) @@ -1655,7 +1686,7 @@ def full_flow_turbopark_solver(farm: Farm, flow_field: FlowField, flow_field_gri # turbine_grid_farm.initialize(turbine_grid.sorted_indices) # turbopark_solver(turbine_grid_farm, turbine_grid_flow_field, turbine_grid, model_manager) - + # flow_field.u = copy.deepcopy(turbine_grid_flow_field.u) # flow_field.v = copy.deepcopy(turbine_grid_flow_field.v) From 5eb59db6b9dd2a4547b8ed02f2f39aad41950405 Mon Sep 17 00:00:00 2001 From: Rob Hammond <13874373+RHammond2@users.noreply.github.com> Date: Tue, 20 Sep 2022 14:59:59 -0600 Subject: [PATCH 07/20] add comment --- floris/simulation/solver.py | 1 + 1 file changed, 1 insertion(+) diff --git a/floris/simulation/solver.py b/floris/simulation/solver.py index fcc1b1dfd..19becd50f 100644 --- a/floris/simulation/solver.py +++ b/floris/simulation/solver.py @@ -803,6 +803,7 @@ def solve( flow_field.turbulence_intensity_field = _expansion_mean(turbine_turbulence_intensity) +# Turn off flake8 for the original code # flake8: noqa def sequential_solver(farm: Farm, flow_field: FlowField, grid: TurbineGrid, model_manager: WakeModelManager) -> None: # Algorithm From a68d2cdbc88c313a4f9ed8b2349c15ac15883450 Mon Sep 17 00:00:00 2001 From: Rob Hammond <13874373+RHammond2@users.noreply.github.com> Date: Tue, 7 Feb 2023 09:05:41 -0700 Subject: [PATCH 08/20] enable future typing patterns --- floris/simulation/solver.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/floris/simulation/solver.py b/floris/simulation/solver.py index a64e0a539..eaa3387f0 100644 --- a/floris/simulation/solver.py +++ b/floris/simulation/solver.py @@ -10,6 +10,8 @@ # License for the specific language governing permissions and limitations under # the License. +from __future__ import annotations + import copy from abc import abstractmethod From 9ed1b17fe43279261cf8ae3e31e71dcd7595fe7d Mon Sep 17 00:00:00 2001 From: Rob Hammond <13874373+RHammond2@users.noreply.github.com> Date: Tue, 7 Feb 2023 09:59:55 -0700 Subject: [PATCH 09/20] fix bad grid reference --- floris/simulation/solver.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/floris/simulation/solver.py b/floris/simulation/solver.py index eaa3387f0..62f314894 100644 --- a/floris/simulation/solver.py +++ b/floris/simulation/solver.py @@ -1047,11 +1047,11 @@ def full_flow_sequential_solver( for i in range(flow_field_grid.n_turbines): # Get the current turbine quantities - x_i = np.mean(grid.x_sorted[:, :, i:i+1], axis=(3, 4)) + x_i = np.mean(turbine_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 = np.mean(turbine_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 = np.mean(turbine_grid.z_sorted[:, :, i:i+1], axis=(3, 4)) z_i = z_i[:, :, :, None, None] u_i = turbine_grid_flow_field.u_sorted[:, :, i:i+1] From e3d462dba423f02df1fb6233d2bf88a0d41f8f00 Mon Sep 17 00:00:00 2001 From: RHammond2 <13874373+RHammond2@users.noreply.github.com> Date: Thu, 13 Jul 2023 08:35:59 -0700 Subject: [PATCH 10/20] bring the refactored solvers up to date with all past changes --- floris/simulation/solver.py | 128 ++++++++++++++++++++++++++++-------- 1 file changed, 101 insertions(+), 27 deletions(-) diff --git a/floris/simulation/solver.py b/floris/simulation/solver.py index fb8790787..570f8d407 100644 --- a/floris/simulation/solver.py +++ b/floris/simulation/solver.py @@ -102,7 +102,7 @@ def solve( *, full_flow: bool = False, farm: Farm = None, - flow_field: FlowFieldGrid = None, + flow_field: FlowFieldGrid | FlowFieldPlanarGrid | PointsGrid = None, grid: TurbineGrid | FlowFieldGrid = None ) -> None: # TODO: Update with the new logger functionality that is on its way @@ -118,21 +118,23 @@ def full_flow_solve(self): if not isinstance(self.grid, FlowFieldGrid): raise TypeError("Cannot run `full_flow_solve` with a `TurbineGrid` object, reinitialize with the input to `grid` being a `FlowFieldGrid` object.") # noqa: E501 - # TODO: Why are we making copies of the originals, and then never using the original, - # can we just use the originally provided values and get on with the calculating? farm = copy.deepcopy(self.farm) flow_field = copy.deepcopy(self.flow_field) farm.construct_turbine_map() farm.construct_turbine_fCts() - farm.construct_turbine_fCps() farm.construct_turbine_power_interps() farm.construct_hub_heights() farm.construct_rotor_diameters() farm.construct_turbine_TSRs() - farm.construc_turbine_pPs() - farm.construc_turbine_ref_density_cp_cts() + farm.construct_turbine_pPs() + farm.construct_turbine_pTs() + farm.construct_turbine_ref_density_cp_cts() + farm.construct_turbine_ref_tilt_cp_cts() + farm.construct_turbine_fTilts() + farm.construct_turbine_correct_cp_ct_for_tilt() farm.construct_coordinates() + farm.set_tilt_to_ref_tilt(flow_field.n_wind_directions, flow_field.n_wind_speeds) turbine_grid = TurbineGrid( turbine_coordinates=farm.coordinates, @@ -158,7 +160,7 @@ def solve( self, *, full_flow: bool = False, - farm: Farm, + farm: Farm | None = None, flow_field: FlowFieldGrid = None, grid: TurbineGrid | FlowFieldGrid = None ) -> None: @@ -215,20 +217,37 @@ def solve( ct_i = Ct( 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=farm.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, )[:, :, 0:1, None, None] # Since we are filtering for the ith turbine in the axial induction function, get the first index here (0:1) axial_induction_i = axial_induction( 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=farm.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, )[:, :, 0:1, None, None] - turbulence_intensity_i = flow_field.turbulence_intensity_field[:, :, i:i+1] + + # TODO: Should the solve and full flow solver actually use two different intensity fields + if full_flow: + turbulence_intensity_i = flow_field.turbulence_intensity_field[:, :, i:i+1] + else: + turbulence_intensity_i = flow_field.turbulence_intensity_field_sorted_avg[:, :, 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] @@ -248,7 +267,8 @@ def solve( hub_height_i, ct_i, TSR_i, - axial_induction_i + axial_induction_i, + flow_field.wind_shear, ) deflection_field = self.model_manager.deflection_model.function( @@ -274,7 +294,8 @@ def solve( yaw_angle_i, ct_i, TSR_i, - axial_induction_i + axial_induction_i, + flow_field.wind_shear, ) if not full_flow: if self.model_manager.enable_yaw_added_recovery: @@ -328,9 +349,9 @@ def solve( ti_added = ( area_overlap * np.nan_to_num(wake_added_turbulence_intensity, posinf=0.0) - * np.array(grid.x_sorted > x_i) - * np.array(np.abs(y_i - grid.y_sorted) < 2 * rotor_diameter_i) - * np.array(grid.x_sorted <= downstream_influence_length + x_i) + * (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 @@ -344,7 +365,8 @@ def solve( flow_field.w_sorted += w_wake if not full_flow: - flow_field.turbulence_intensity_field = _expansion_mean(turbine_turbulence_intensity) + flow_field.turbulence_intensity_field_sorted = turbine_turbulence_intensity + flow_field.flow_field.turbulence_intensity_field_sorted_avg = _expansion_mean(turbine_turbulence_intensity) @define(auto_attribs=True) @@ -414,11 +436,12 @@ def solve( v_i = flow_field.v_sorted[:, :, i:i+1] if not full_flow: + rotor_diameter_i = farm.rotor_diameters_sorted[: ,:, i:i+1, None, None] mask = ( (grid.x_sorted < x_i + 0.01) * (grid.x_sorted > x_i - 0.01) - * (grid.y_sorted < y_i + 0.51 * 126.0) - * (grid.y_sorted > y_i - 0.51 * 126.0) + * (grid.y_sorted < y_i + 0.51 * rotor_diameter_i) + * (grid.y_sorted > y_i - 0.51 * rotor_diameter_i) ) turbine_inflow_field *= ~mask + (flow_field.u_initial_sorted - turb_u_wake) * mask @@ -426,30 +449,52 @@ def solve( turb_avg_vels = average_velocity(turbine_inflow_field) turb_Cts = Ct( - turb_avg_vels, - farm.yaw_angles_sorted, - farm.turbine_fCts, + velocities=turb_avg_vels, + yaw_angle=farm.yaw_angles_sorted, + tilt_angle=farm.tilt_angles_sorted, + ref_tilt_cp_ct=farm.ref_tilt_cp_cts_sorted, + fCt=farm.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, + average_method=grid.average_method, + cubature_weights=grid.cubature_weights, )[:, :, :, None, None] if not full_flow: turb_aIs = axial_induction( - turb_avg_vels, - farm.yaw_angles_sorted, - farm.turbine_fCts, + velocities=turb_avg_vels, + yaw_angle=farm.yaw_angles_sorted, + tilt_angle=farm.tilt_angles_sorted, + ref_tilt_cp_ct=farm.ref_tilt_cp_cts_sorted, + fCt=farm.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, )[:, :, :, None, None] axial_induction_i = axial_induction( 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=farm.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, )[:, :, :, None, None] - turbulence_intensity_i = turbine_turbulence_intensity[:, :, i:i+1] + # TODO: Should the solve and full flow solver actually use two different intensity fields + if full_flow: + turbulence_intensity_i = flow_field.turbulence_intensity_field[:, :, i:i+1] + else: + turbulence_intensity_i = flow_field.turbulence_intensity_field_sorted_avg[:, :, 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] @@ -468,6 +513,7 @@ def solve( turb_Cts[:, :, i:i+1], TSR_i, axial_induction_i, + flow_field.wind_shear, scale=scale_factor, ) @@ -497,6 +543,7 @@ def solve( turb_Cts[:, :, i:i+1], TSR_i, axial_induction_i, + flow_field.wind_shear, scale=scale_factor, ) @@ -650,26 +697,44 @@ def solve( Cts = Ct( 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=farm.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, + average_method=grid.average_method, + cubature_weights=grid.cubature_weights, ) # Since we are filtering for the ith turbine in the Ct function, get the first index here (0:1) ct_i = Ct( 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=farm.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, )[:, :, 0:1, None, None] # Since we are filtering for the ith turbine in the axial induction function, get the first index here (0:1) axial_induction_i = axial_induction( 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=farm.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, )[:, :, 0:1, None, None] turbulence_intensity_i = turbine_turbulence_intensity[:, :, i:i+1] @@ -692,7 +757,8 @@ def solve( hub_height_i, ct_i, TSR_i, - axial_induction_i + axial_induction_i, + flow_field.wind_shear, ) # Model calculations @@ -708,9 +774,15 @@ def solve( ct_ii = Ct( 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=farm.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=[ii] + ix_filter=[ii], + average_method=grid.average_method, + cubature_weights=grid.cubature_weights, )[:, :, 0:1, None, None] rotor_diameter_ii = farm.rotor_diameters_sorted[:, :, ii:ii+1, None, None] @@ -739,7 +811,8 @@ def solve( yaw_angle_i, ct_i, TSR_i, - axial_induction_i + axial_induction_i, + flow_field.wind_shear, ) if self.model_manager.enable_yaw_added_recovery: @@ -808,7 +881,8 @@ def solve( flow_field.v_sorted += v_wake flow_field.w_sorted += w_wake - flow_field.turbulence_intensity_field = _expansion_mean(turbine_turbulence_intensity) + flow_field.turbulence_intensity_field_sorted = turbine_turbulence_intensity + flow_field.flow_field.turbulence_intensity_field_sorted_avg = _expansion_mean(turbine_turbulence_intensity) # Turn off flake8 for the original code # flake8: noqa From e78158f04ee622e5f8b8328814d06c231435c886 Mon Sep 17 00:00:00 2001 From: RHammond2 <13874373+RHammond2@users.noreply.github.com> Date: Thu, 13 Jul 2023 08:58:29 -0700 Subject: [PATCH 11/20] add the trailing commas back --- floris/simulation/solver.py | 43 ++++++++++++++++++------------------- 1 file changed, 21 insertions(+), 22 deletions(-) diff --git a/floris/simulation/solver.py b/floris/simulation/solver.py index 570f8d407..452fcf49c 100644 --- a/floris/simulation/solver.py +++ b/floris/simulation/solver.py @@ -38,7 +38,6 @@ yaw_added_turbulence_mixing, ) from floris.type_dec import NDArrayFloat -from floris.utilities import cosd def _expansion_mean(x: NDArrayFloat) -> NDArrayFloat: @@ -90,7 +89,7 @@ class Solver: flow_field: FlowField = field(converter=copy.deepcopy, validator=attrs.validators.instance_of(Farm)) grid: TurbineGrid | FlowFieldGrid = field( converter=copy.deepcopy, - validator=attrs.validators.instance_of((FlowFieldGrid, TurbineGrid)) + validator=attrs.validators.instance_of((FlowFieldGrid, TurbineGrid)), ) model_manager: WakeModelManager = field( converter=copy.deepcopy, validator=attrs.validators.instance_of(WakeModelManager) @@ -101,9 +100,9 @@ def solve( self, *, full_flow: bool = False, - farm: Farm = None, - flow_field: FlowFieldGrid | FlowFieldPlanarGrid | PointsGrid = None, - grid: TurbineGrid | FlowFieldGrid = None + farm: Farm | None = None, + flow_field: FlowFieldGrid | None = None, + grid: TurbineGrid | FlowFieldGrid | None = None, ) -> None: # TODO: Update with the new logger functionality that is on its way pass @@ -162,7 +161,7 @@ def solve( full_flow: bool = False, farm: Farm | None = None, flow_field: FlowFieldGrid = None, - grid: TurbineGrid | FlowFieldGrid = None + grid: TurbineGrid | FlowFieldGrid = None, ) -> None: """Runs the sequential sover methodology, or full flow sequential solver methodology for a wind farm. @@ -278,7 +277,7 @@ def solve( turbulence_intensity_i, ct_i, rotor_diameter_i, - **deflection_model_args + **deflection_model_args, ) if self.model_manager.enable_transverse_velocities: @@ -321,12 +320,12 @@ def solve( ct_i, hub_height_i, rotor_diameter_i, - **deficit_model_args + **deficit_model_args, ) wake_field = self.model_manager.combination_model.function( wake_field, - velocity_deficit * flow_field.u_initial_sorted + velocity_deficit * flow_field.u_initial_sorted, ) if not full_flow: @@ -335,7 +334,7 @@ def solve( grid.x_sorted, x_i, rotor_diameter_i, - axial_induction_i + axial_induction_i, ) # Calculate wake overlap for wake-added turbulence (WAT) @@ -357,7 +356,7 @@ def solve( # Combine turbine TIs with WAT turbine_turbulence_intensity = np.maximum( np.sqrt(ti_added ** 2 + ambient_turbulence_intensity ** 2), - turbine_turbulence_intensity + turbine_turbulence_intensity, ) flow_field.u_sorted = flow_field.u_initial_sorted - wake_field @@ -419,7 +418,7 @@ def solve( turbine_inflow_field = copy.deepcopy(flow_field.u_initial_sorted) turbine_turbulence_intensity = np.full( (flow_field.n_wind_directions, flow_field.n_wind_speeds, farm.n_turbines, 1, 1), - flow_field.turbulence_intensity + flow_field.turbulence_intensity, ) ambient_turbulence_intensity = flow_field.turbulence_intensity @@ -526,7 +525,7 @@ def solve( turbulence_intensity_i, turb_Cts[:, :, i:i+1], rotor_diameter_i, - **deflection_model_args + **deflection_model_args, ) if self.model_manager.enable_transverse_velocities: @@ -575,7 +574,7 @@ def solve( farm.rotor_diameters_sorted[:, :, :, None, None], turb_u_wake, Ctmp, - **deficit_model_args + **deficit_model_args, ) if not full_flow: @@ -584,7 +583,7 @@ def solve( grid.x_sorted, x_i, rotor_diameter_i, - turb_aIs + turb_aIs, ) # Calculate wake overlap for wake-added turbulence (WAT) @@ -607,7 +606,7 @@ def solve( # Combine turbine TIs with WAT turbine_turbulence_intensity = np.maximum( np.sqrt(ti_added ** 2 + ambient_turbulence_intensity ** 2), - turbine_turbulence_intensity + turbine_turbulence_intensity, ) flow_field.v_sorted += v_wake @@ -632,7 +631,7 @@ def solve( full_flow: bool = False, farm: Farm = None, flow_field: FlowField = None, - grid: TurbineGrid | FlowFieldGrid = None + grid: TurbineGrid | FlowFieldGrid = None, ) -> None: """Runs the TurbOPark sover methodology, or full flow TurbOPark solver methodology for a wind farm. @@ -680,7 +679,7 @@ def solve( turbine_turbulence_intensity = np.full( (flow_field.n_wind_directions, flow_field.n_wind_speeds, farm.n_turbines, 1, 1), - flow_field.turbulence_intensity + flow_field.turbulence_intensity, ) ambient_turbulence_intensity = flow_field.turbulence_intensity @@ -793,7 +792,7 @@ def solve( turbulence_intensity_ii, ct_ii, rotor_diameter_ii, - **deflection_model_args + **deflection_model_args, ) deflection_field[:, :, ii:ii+1, :, :] = deflection_field_ii[:, :, i:i+1, :, :] @@ -837,7 +836,7 @@ def solve( farm.rotor_diameters_sorted[:, :, :, None, None], i, deflection_field, - **deficit_model_args + **deficit_model_args, ) wake_field = self.model_manager.combination_model.function( @@ -850,7 +849,7 @@ def solve( grid.x_sorted, x_i, rotor_diameter_i, - axial_induction_i + axial_induction_i, ) # TODO: leaving this in for GCH quantities; will need to find another way to compute area_overlap @@ -874,7 +873,7 @@ def solve( # Combine turbine TIs with WAT turbine_turbulence_intensity = np.maximum( np.sqrt(ti_added ** 2 + ambient_turbulence_intensity ** 2), - turbine_turbulence_intensity + turbine_turbulence_intensity, ) flow_field.u_sorted = flow_field.u_initial_sorted - wake_field From 80c07ec96d6ff143fb871c7f8cbe545fe8303991 Mon Sep 17 00:00:00 2001 From: RHammond2 <13874373+RHammond2@users.noreply.github.com> Date: Thu, 13 Jul 2023 10:23:04 -0700 Subject: [PATCH 12/20] udpate the solver method typing --- floris/simulation/solver.py | 21 ++++++++++++--------- 1 file changed, 12 insertions(+), 9 deletions(-) diff --git a/floris/simulation/solver.py b/floris/simulation/solver.py index c880f8e53..9b84c7228 100644 --- a/floris/simulation/solver.py +++ b/floris/simulation/solver.py @@ -160,8 +160,8 @@ def solve( *, full_flow: bool = False, farm: Farm | None = None, - flow_field: FlowFieldGrid = None, - grid: TurbineGrid | FlowFieldGrid = None, + flow_field: FlowField = None, + grid: TurbineGrid | FlowFieldGrid | FlowFieldPlanarGrid | PointsGrid = None, ) -> None: """Runs the sequential sover methodology, or full flow sequential solver methodology for a wind farm. @@ -174,11 +174,12 @@ def solve( Defaults to None. flow_field (FlowField, optional): Allows for a non-initialized `flow_field` object to be used. It should be noted that this functionality is intended for use with - `full_flow_solve`.Defaults to None. + `full_flow_solve`. Defaults to None. grid (TurbineGrid | FlowFieldGrid, optional): Allows for a non-initialized `grid` object to be used. It should be noted that this functionality is intended for use with - `full_flow_solve`, which computes over a `TurbineGrid`, then a `FlowFieldGrid`. - Defaults to None. + `full_flow_solve`, which computes over a `TurbineGrid`, then a `FlowFieldGrid`. If + `full_flow=True`, then `grid` should be one of `FlowFieldGrid`, `FlowFieldPlanarGrid`, + or `PointsGrid`. Defaults to None. """ if farm is None: farm = self.farm @@ -392,8 +393,9 @@ def solve( `full_flow_solve`.Defaults to None. grid (TurbineGrid | FlowFieldGrid, optional): Allows for a non-initialized `grid` object to be used. It should be noted that this functionality is intended for use with - `full_flow_solve`, which computes over a `TurbineGrid`, then a `FlowFieldGrid`. - Defaults to None. + `full_flow_solve`, which computes over a `TurbineGrid`, then a `FlowFieldGrid`. If + `full_flow=False`, this should be a `TurbineGrid`, and if `full_flow=True`, this + should be a `FlowFieldGrid`. Defaults to None. """ if farm is None: farm = self.farm @@ -647,8 +649,9 @@ def solve( `full_flow_solve`.Defaults to None. grid (TurbineGrid | FlowFieldGrid, optional): Allows for a non-initialized `grid` object to be used. It should be noted that this functionality is intended for use with - `full_flow_solve`, which computes over a `TurbineGrid`, then a `FlowFieldGrid`. - Defaults to None. + `full_flow_solve`, which computes over a `TurbineGrid`, then a `FlowFieldGrid`. If + `full_flow=False`, this should be a `TurbineGrid`, and if `full_flow=True`, this + should be a `FlowFieldGrid`. Defaults to None. """ # Algorithm # For each turbine, calculate its effect on every downstream turbine. From ea9d239407ffa01529370e856e938d17aaa3f697 Mon Sep 17 00:00:00 2001 From: RHammond2 <13874373+RHammond2@users.noreply.github.com> Date: Thu, 13 Jul 2023 11:57:11 -0700 Subject: [PATCH 13/20] integrate the sequential solver and fix a logic error --- floris/simulation/__init__.py | 3 +++ floris/simulation/floris.py | 23 +++++++++++++++-------- floris/simulation/solver.py | 10 +++++----- 3 files changed, 23 insertions(+), 13 deletions(-) diff --git a/floris/simulation/__init__.py b/floris/simulation/__init__.py index 12d30aab8..6af12a4da 100644 --- a/floris/simulation/__init__.py +++ b/floris/simulation/__init__.py @@ -58,6 +58,9 @@ full_flow_turbopark_solver, sequential_solver, turbopark_solver, + SequentialSolver, + CCSolver, + TurbOParkSolver ) from .floris import Floris diff --git a/floris/simulation/floris.py b/floris/simulation/floris.py index a5f96c2a4..75120648b 100644 --- a/floris/simulation/floris.py +++ b/floris/simulation/floris.py @@ -35,6 +35,7 @@ Grid, PointsGrid, sequential_solver, + SequentialSolver, State, TurbineCubatureGrid, TurbineGrid, @@ -244,12 +245,14 @@ def steady_state_atmospheric_condition(self): self.wake ) else: - sequential_solver( - self.farm, - self.flow_field, - self.grid, - self.wake - ) + # sequential_solver( + # self.farm, + # self.flow_field, + # self.grid, + # self.wake + # ) + solver = SequentialSolver(self.farm, self.flow_field, self.grid, self.wake) + solver.solve() # end = time.time() # elapsed_time = end - start @@ -274,7 +277,9 @@ def solve_for_viz(self): elif vel_model=="empirical_gauss": full_flow_empirical_gauss_solver(self.farm, self.flow_field, self.grid, self.wake) else: - full_flow_sequential_solver(self.farm, self.flow_field, self.grid, self.wake) + # full_flow_sequential_solver(self.farm, self.flow_field, self.grid, self.wake) + solver = SequentialSolver(self.farm, self.flow_field, self.grid, self.wake) + solver.solve(full_flow=True) def solve_for_points(self, x, y, z): # Do the calculation with the TurbineGrid for a single wind speed @@ -310,7 +315,9 @@ def solve_for_points(self, x, y, z): elif vel_model == "empirical_gauss": full_flow_empirical_gauss_solver(self.farm, self.flow_field, field_grid, self.wake) else: - full_flow_sequential_solver(self.farm, self.flow_field, field_grid, self.wake) + # full_flow_sequential_solver(self.farm, self.flow_field, field_grid, self.wake) + solver = SequentialSolver(self.farm, self.flow_field, self.grid, self.wake) + solver.solve(full_flow=True) return self.flow_field.u_sorted[:,:,:,0,0] # Remove turbine grid dimensions diff --git a/floris/simulation/solver.py b/floris/simulation/solver.py index 9b84c7228..deb2d7617 100644 --- a/floris/simulation/solver.py +++ b/floris/simulation/solver.py @@ -86,7 +86,7 @@ def calculate_area_overlap(wake_velocities, freestream_velocities, y_ngrid, z_ng @define(auto_attribs=True) class Solver: farm: Farm = field(converter=copy.deepcopy, validator=attrs.validators.instance_of(Farm)) - flow_field: FlowField = field(converter=copy.deepcopy, validator=attrs.validators.instance_of(Farm)) + flow_field: FlowField = field(converter=copy.deepcopy, validator=attrs.validators.instance_of(FlowField)) grid: TurbineGrid | FlowFieldGrid = field( converter=copy.deepcopy, validator=attrs.validators.instance_of((FlowFieldGrid, TurbineGrid)), @@ -245,9 +245,9 @@ def solve( # TODO: Should the solve and full flow solver actually use two different intensity fields if full_flow: - turbulence_intensity_i = flow_field.turbulence_intensity_field[:, :, i:i+1] - else: turbulence_intensity_i = flow_field.turbulence_intensity_field_sorted_avg[:, :, i:i+1] + else: + 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] @@ -366,7 +366,7 @@ def solve( if not full_flow: flow_field.turbulence_intensity_field_sorted = turbine_turbulence_intensity - flow_field.flow_field.turbulence_intensity_field_sorted_avg = _expansion_mean(turbine_turbulence_intensity) + flow_field.turbulence_intensity_field_sorted_avg = _expansion_mean(turbine_turbulence_intensity) @define(auto_attribs=True) @@ -884,7 +884,7 @@ def solve( flow_field.w_sorted += w_wake flow_field.turbulence_intensity_field_sorted = turbine_turbulence_intensity - flow_field.flow_field.turbulence_intensity_field_sorted_avg = _expansion_mean(turbine_turbulence_intensity) + flow_field.turbulence_intensity_field_sorted_avg = _expansion_mean(turbine_turbulence_intensity) # Turn off flake8 for the original code # flake8: noqa From 0c4067d5813b23c2934f31cd46c3262c3cd5b490 Mon Sep 17 00:00:00 2001 From: RHammond2 <13874373+RHammond2@users.noreply.github.com> Date: Thu, 13 Jul 2023 14:52:18 -0700 Subject: [PATCH 14/20] reduce small piece to one line and fix logic bug in if/else --- floris/simulation/solver.py | 8 +++----- 1 file changed, 3 insertions(+), 5 deletions(-) diff --git a/floris/simulation/solver.py b/floris/simulation/solver.py index deb2d7617..1b243a36a 100644 --- a/floris/simulation/solver.py +++ b/floris/simulation/solver.py @@ -253,9 +253,7 @@ def solve( 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 - + effective_yaw_i = yaw_angle_i if self.model_manager.enable_secondary_steering: effective_yaw_i += wake_added_yaw( u_i, @@ -493,9 +491,9 @@ def solve( # TODO: Should the solve and full flow solver actually use two different intensity fields if full_flow: - turbulence_intensity_i = flow_field.turbulence_intensity_field[:, :, i:i+1] - else: turbulence_intensity_i = flow_field.turbulence_intensity_field_sorted_avg[:, :, i:i+1] + else: + turbulence_intensity_i = flow_field.turbulence_intensity_field[:, :, 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] From 61c15a8ecb6b218a1120558a2de27aa5ac15f39a Mon Sep 17 00:00:00 2001 From: RHammond2 <13874373+RHammond2@users.noreply.github.com> Date: Thu, 13 Jul 2023 16:06:11 -0700 Subject: [PATCH 15/20] fix remaining inconsistencies in the sequential solver and bring up to spec --- floris/simulation/__init__.py | 5 ++- floris/simulation/floris.py | 81 +++++++++++++++++++---------------- floris/simulation/solver.py | 17 +++++++- 3 files changed, 64 insertions(+), 39 deletions(-) diff --git a/floris/simulation/__init__.py b/floris/simulation/__init__.py index 6af12a4da..403945dc8 100644 --- a/floris/simulation/__init__.py +++ b/floris/simulation/__init__.py @@ -51,7 +51,9 @@ from .wake import WakeModelManager from .solver import ( cc_solver, + CCSolver, empirical_gauss_solver, + EmpiricalGaussSolver, full_flow_cc_solver, full_flow_empirical_gauss_solver, full_flow_sequential_solver, @@ -59,8 +61,7 @@ sequential_solver, turbopark_solver, SequentialSolver, - CCSolver, - TurbOParkSolver + TurbOParkSolver, ) from .floris import Floris diff --git a/floris/simulation/floris.py b/floris/simulation/floris.py index 75120648b..a41292179 100644 --- a/floris/simulation/floris.py +++ b/floris/simulation/floris.py @@ -23,7 +23,9 @@ from floris.simulation import ( BaseClass, cc_solver, + CCSolver, empirical_gauss_solver, + EmpiricalGaussSolver, Farm, FlowField, FlowFieldGrid, @@ -40,6 +42,7 @@ TurbineCubatureGrid, TurbineGrid, turbopark_solver, + TurbOParkSolver, WakeModelManager, ) from floris.utilities import load_yaml @@ -66,6 +69,7 @@ class Floris(BaseClass): floris_version: str = field(converter=str) grid: Grid = field(init=False) + solve: SequentialSolver | CCSolver | TurbOParkSolver | EmpiricalGaussSolver = field(init=False) def __attrs_post_init__(self) -> None: @@ -224,35 +228,16 @@ def steady_state_atmospheric_condition(self): ) if vel_model=="cc": - cc_solver( - self.farm, - self.flow_field, - self.grid, - self.wake - ) + self.solve = CCSolver(self.farm, self.flow_field, self.grid, self.wake) elif vel_model=="turbopark": - turbopark_solver( - self.farm, - self.flow_field, - self.grid, - self.wake - ) + self.solve = TurbOParkSolver(self.farm, self.flow_field, self.grid, self.wake) elif vel_model=="empirical_gauss": - empirical_gauss_solver( - self.farm, - self.flow_field, - self.grid, - self.wake - ) + self.solve = EmpiricalGaussSolver(self.farm, self.flow_field, self.grid, self.wake) else: - # sequential_solver( - # self.farm, - # self.flow_field, - # self.grid, - # self.wake - # ) - solver = SequentialSolver(self.farm, self.flow_field, self.grid, self.wake) - solver.solve() + self.solve = SequentialSolver(self.farm, self.flow_field, self.grid, self.wake) + + self.solve.solve() + self.post_solve_update_flow_field() # end = time.time() # elapsed_time = end - start @@ -271,15 +256,16 @@ def solve_for_viz(self): vel_model = self.wake.model_strings["velocity_model"] if vel_model=="cc": - full_flow_cc_solver(self.farm, self.flow_field, self.grid, self.wake) + self.solve = CCSolver(self.farm, self.flow_field, self.grid, self.wake) elif vel_model=="turbopark": - full_flow_turbopark_solver(self.farm, self.flow_field, self.grid, self.wake) + self.solve = TurbOParkSolver(self.farm, self.flow_field, self.grid, self.wake) elif vel_model=="empirical_gauss": - full_flow_empirical_gauss_solver(self.farm, self.flow_field, self.grid, self.wake) + self.solve = EmpiricalGaussSolver(self.farm, self.flow_field, self.grid, self.wake) else: - # full_flow_sequential_solver(self.farm, self.flow_field, self.grid, self.wake) - solver = SequentialSolver(self.farm, self.flow_field, self.grid, self.wake) - solver.solve(full_flow=True) + self.solve = SequentialSolver(self.farm, self.flow_field, self.grid, self.wake) + + self.solve.solve(full_flow=True) + self.post_solve_update_flow_field() def solve_for_points(self, x, y, z): # Do the calculation with the TurbineGrid for a single wind speed @@ -312,15 +298,38 @@ def solve_for_points(self, x, y, z): "solve_for_points is currently only available with the "+\ "gauss, jensen, and empirical_guass models." ) - elif vel_model == "empirical_gauss": - full_flow_empirical_gauss_solver(self.farm, self.flow_field, field_grid, self.wake) + + if vel_model == "empirical_gauss": + self.solve = EmpiricalGaussSolver(self.farm, self.flow_field, field_grid, self.wake) else: # full_flow_sequential_solver(self.farm, self.flow_field, field_grid, self.wake) - solver = SequentialSolver(self.farm, self.flow_field, self.grid, self.wake) - solver.solve(full_flow=True) + self.solve = SequentialSolver(self.farm, self.flow_field, self.grid, self.wake) + + self.solve.solve(full_flow=True) + self.post_solve_update_flow_field() return self.flow_field.u_sorted[:,:,:,0,0] # Remove turbine grid dimensions + def post_solve_update_flow_field(self): + """Updates the `flow_field` values with those that were solved during the steady state + calculation. + """ + self.flow_field.u = self.solve.flow_field.u + self.flow_field.v = self.solve.flow_field.v + self.flow_field.w = self.solve.flow_field.w + self.flow_field.u_sorted = self.solve.flow_field.u_sorted + self.flow_field.v_sorted = self.solve.flow_field.v_sorted + self.flow_field.w_sorted = self.solve.flow_field.w_sorted + self.flow_field.turbulence_intensity_field = ( + self.solve.flow_field.turbulence_intensity_field + ) + self.flow_field.turbulence_intensity_field_sorted = ( + self.solve.flow_field.turbulence_intensity_field_sorted + ) + self.flow_field.turbulence_intensity_field_sorted_avg = ( + self.solve.flow_field.turbulence_intensity_field_sorted_avg + ) + def finalize(self): # Once the wake calculation is finished, unsort the values to match # the user-supplied order of things. diff --git a/floris/simulation/solver.py b/floris/simulation/solver.py index 1b243a36a..98016f345 100644 --- a/floris/simulation/solver.py +++ b/floris/simulation/solver.py @@ -253,7 +253,7 @@ def solve( 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 = yaw_angle_i + effective_yaw_i = yaw_angle_i.copy() if self.model_manager.enable_secondary_steering: effective_yaw_i += wake_added_yaw( u_i, @@ -884,6 +884,21 @@ def solve( flow_field.turbulence_intensity_field_sorted = turbine_turbulence_intensity flow_field.turbulence_intensity_field_sorted_avg = _expansion_mean(turbine_turbulence_intensity) + +@define(auto_attribs=True) +class EmpiricalGaussSolver(Solver): + + def solve( + self, + *, + full_flow: bool = False, + farm: Farm = None, + flow_field: FlowField = None, + grid: TurbineGrid | FlowFieldGrid = None, + ) -> None: + ... + + # Turn off flake8 for the original code # flake8: noqa def sequential_solver(farm: Farm, flow_field: FlowField, grid: TurbineGrid, model_manager: WakeModelManager) -> None: From 00990fdd3fc4a586e0fff31ca245947c100af5e8 Mon Sep 17 00:00:00 2001 From: RHammond2 <13874373+RHammond2@users.noreply.github.com> Date: Fri, 14 Jul 2023 08:54:05 -0700 Subject: [PATCH 16/20] fix small logic bugs in cc solver --- floris/simulation/solver.py | 19 +++++++++++-------- 1 file changed, 11 insertions(+), 8 deletions(-) diff --git a/floris/simulation/solver.py b/floris/simulation/solver.py index 98016f345..b768ff8a2 100644 --- a/floris/simulation/solver.py +++ b/floris/simulation/solver.py @@ -434,7 +434,9 @@ def solve( u_i = turbine_inflow_field[:, :, i:i+1] v_i = flow_field.v_sorted[:, :, i:i+1] - if not full_flow: + if full_flow: + turbine_inflow_field = flow_field.u_sorted + else: rotor_diameter_i = farm.rotor_diameters_sorted[: ,:, i:i+1, None, None] mask = ( (grid.x_sorted < x_i + 0.01) @@ -444,8 +446,6 @@ def solve( ) turbine_inflow_field *= ~mask + (flow_field.u_initial_sorted - turb_u_wake) * mask - turbine_inflow_field = flow_field.u_sorted if full_flow else turbine_inflow_field - turb_avg_vels = average_velocity(turbine_inflow_field) turb_Cts = Ct( velocities=turb_avg_vels, @@ -493,7 +493,7 @@ def solve( if full_flow: turbulence_intensity_i = flow_field.turbulence_intensity_field_sorted_avg[:, :, i:i+1] else: - turbulence_intensity_i = flow_field.turbulence_intensity_field[:, :, i:i+1] + 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] @@ -546,9 +546,7 @@ def solve( scale=scale_factor, ) - if full_flow: - turbine_turbulence_intensity = flow_field.turbulence_intensity_field - else: + if not full_flow: if self.model_manager.enable_yaw_added_recovery: I_mixing = yaw_added_turbulence_mixing( u_i, @@ -561,6 +559,10 @@ def solve( turbine_turbulence_intensity[:, :, i:i+1] = turbulence_intensity_i + gch_gain * I_mixing # NOTE: exponential + if full_flow: + ti = turbine_grid_flow_field.turbulence_intensity_field_sorted_avg + else: + ti = turbine_turbulence_intensity turb_u_wake, Ctmp = self.model_manager.velocity_model.function( i, x_i, @@ -569,7 +571,7 @@ def solve( u_i, deflection_field, yaw_angle_i, - turbine_turbulence_intensity, + ti, turb_Cts, farm.rotor_diameters_sorted[:, :, :, None, None], turb_u_wake, @@ -615,6 +617,7 @@ def solve( flow_field.u_sorted = flow_field.u_initial_sorted - turb_u_wake if full_flow else turbine_inflow_field if not full_flow: + flow_field.turbulence_intensity_field_sorted = turbine_turbulence_intensity flow_field.turbulence_intensity_field = _expansion_mean(turbine_turbulence_intensity) From 9ed672b9d9a028814ded95f068dbe2b39f6f1ade Mon Sep 17 00:00:00 2001 From: RHammond2 <13874373+RHammond2@users.noreply.github.com> Date: Fri, 14 Jul 2023 11:24:55 -0700 Subject: [PATCH 17/20] convert empirical gauss solver to EmpiricalGauss.solve() --- floris/simulation/solver.py | 191 +++++++++++++++++++++++++++++++++++- 1 file changed, 187 insertions(+), 4 deletions(-) diff --git a/floris/simulation/solver.py b/floris/simulation/solver.py index b768ff8a2..e8e0e4260 100644 --- a/floris/simulation/solver.py +++ b/floris/simulation/solver.py @@ -746,8 +746,7 @@ def solve( 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 + effective_yaw_i = yaw_angle_i.copy() if self.model_manager.enable_secondary_steering: effective_yaw_i += wake_added_yaw( @@ -897,9 +896,193 @@ def solve( full_flow: bool = False, farm: Farm = None, flow_field: FlowField = None, - grid: TurbineGrid | FlowFieldGrid = None, + grid: TurbineGrid | FlowFieldGrid = None ) -> None: - ... + """Runs the Empirical Gauss sover methodology, or full flow Empirical Gauss solver + methodology for a wind farm. + + Args: + full_flow (bool, optional): Runs the full flow solver when True, and the standard + Empirical Gauss solver, when False. Defaults to False. + farm (Farm, optional): Allows for a non-initialized `farm` object to be used. It should + be noted that this functionality is intended for use with `full_flow_solve`. + Defaults to None. + flow_field (FlowField, optional): Allows for a non-initialized `flow_field` object to be + used. It should be noted that this functionality is intended for use with + `full_flow_solve`.Defaults to None. + grid (TurbineGrid | FlowFieldGrid, optional): Allows for a non-initialized `grid` object + to be used. It should be noted that this functionality is intended for use with + `full_flow_solve`, which computes over a `TurbineGrid`, then a `FlowFieldGrid`. If + `full_flow=False`, this should be a `TurbineGrid`, and if `full_flow=True`, this + should be a `FlowFieldGrid`. Defaults to None. + """ + if farm is None: + farm = self.farm + if flow_field is None: + flow_field = self.flow_field + if grid is None: + grid = self.grid + + gch_gain = 1.0 + scale_factor = 2.0 + + # <> + deflection_model_args = self.model_manager.deflection_model.prepare_function(grid, flow_field) + deficit_model_args = self.model_manager.velocity_model.prepare_function(grid, flow_field) + + # 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) + + x_locs = np.mean(grid.x_sorted, axis=(3, 4))[:,:,:,None] + downstream_distance_D = np.maximum( + (x_locs - np.transpose(x_locs, axes=(0,1,3,2))) + / np.repeat(farm.rotor_diameters_sorted[:, :, :, None], grid.n_turbines, axis=-1), + 0.1 + ) # Max for ease + mixing_factor = np.zeros_like(downstream_distance_D) + mixing_factor[:,:,:,:] = ( + self.model_manager.turbulence_model.atmospheric_ti_gain + * flow_field.turbulence_intensity + * np.eye(grid.n_turbines) + ) + + # Calculate the velocity deficit sequentially from upstream to downstream turbines + for i in range(grid.n_turbines): + + # Get the current turbine quantities + x_i = _expansion_mean_i(grid.x_sorted, i) + y_i = _expansion_mean_i(grid.y_sorted, i) + z_i = _expansion_mean_i(grid.z_sorted, i) + u_i = flow_field.u_sorted[:, :, i:i+1] + v_i = flow_field.v_sorted[:, :, i:i+1] + + # Since we are filtering for the ith turbine in the Ct function, get the first index here (0:1) + ct_i = Ct( + 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=farm.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, + )[:, :, 0:1, None, None] + + # Since we are filtering for the ith turbine in the axial induction function, get the first index here (0:1) + axial_induction_i = axial_induction( + 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=farm.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, + )[:, :, 0:1, None, None] + + 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] + + effective_yaw_i = yaw_angle_i.copy() + + average_velocities = average_velocity( + flow_field.u_sorted, + method=grid.average_method, + cubature_weights=grid.cubature_weights + ) + tilt_angle_i = farm.calculate_tilt_for_eff_velocities(average_velocities)[:, :, i:i+1, None, None] + + if self.model_manager.enable_secondary_steering: + raise NotImplementedError( + "Secondary steering not available for this model.") + + if self.model_manager.enable_transverse_velocities: + raise NotImplementedError( + "Transverse velocities not used in this model.") + + if self.model_manager.enable_yaw_added_recovery: + # Influence of yawing on turbine's own wake + mixing_factor[:, :, i:i+1, i] += \ + yaw_added_wake_mixing( + axial_induction_i, + yaw_angle_i, + 1, + self.model_manager.deflection_model.yaw_added_mixing_gain + ) + + # Extract total wake induced mixing for turbine i + mixing_i = np.linalg.norm( + mixing_factor[:, :, i:i+1, :, None], + ord=2, axis=3, keepdims=True + ) + + # Model calculations + # NOTE: exponential + deflection_field_y, deflection_field_z = self.model_manager.deflection_model.function( + x_i, + y_i, + effective_yaw_i, + tilt_angle_i, + mixing_i, + ct_i, + rotor_diameter_i, + **deflection_model_args + ) + + # NOTE: exponential + velocity_deficit = self.model_manager.velocity_model.function( + x_i, + y_i, + z_i, + axial_induction_i, + deflection_field_y, + deflection_field_z, + yaw_angle_i, + tilt_angle_i, + mixing_i, + ct_i, + hub_height_i, + rotor_diameter_i, + **deficit_model_args + ) + + wake_field = self.model_manager.combination_model.function( + wake_field, + velocity_deficit * flow_field.u_initial_sorted + ) + + # 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) + + # Compute wake induced mixing factor + mixing_factor[:, :, :, i] += \ + area_overlap * self.model_manager.turbulence_model.function( + axial_induction_i, downstream_distance_D[:,:,:,i] + ) + if self.model_manager.enable_yaw_added_recovery: + mixing_factor[:,:,:,i] += \ + area_overlap * yaw_added_wake_mixing( + axial_induction_i, + yaw_angle_i, + downstream_distance_D[:,:,:,i], + self.model_manager.deflection_model.yaw_added_mixing_gain + ) + + flow_field.u_sorted = flow_field.u_initial_sorted - wake_field + flow_field.v_sorted += v_wake + flow_field.w_sorted += w_wake + + return mixing_factor # Turn off flake8 for the original code From 211f4bdfc8cc7bce6a8c08d4500c01b370aed2fb Mon Sep 17 00:00:00 2001 From: RHammond2 <13874373+RHammond2@users.noreply.github.com> Date: Fri, 14 Jul 2023 14:26:48 -0700 Subject: [PATCH 18/20] initialize solver with a value and better check types --- floris/simulation/floris.py | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/floris/simulation/floris.py b/floris/simulation/floris.py index a41292179..f380c24b8 100644 --- a/floris/simulation/floris.py +++ b/floris/simulation/floris.py @@ -16,6 +16,7 @@ from pathlib import Path +import attrs import yaml from attrs import define, field @@ -69,7 +70,13 @@ class Floris(BaseClass): floris_version: str = field(converter=str) grid: Grid = field(init=False) - solve: SequentialSolver | CCSolver | TurbOParkSolver | EmpiricalGaussSolver = field(init=False) + solve: SequentialSolver | CCSolver | TurbOParkSolver | EmpiricalGaussSolver = field( + default=None, + init=False, + validator=attrs.validators.instance_of( + (SequentialSolver, CCSolver, TurbOParkSolver, EmpiricalGaussSolver, type(None)) + ) + ) def __attrs_post_init__(self) -> None: From 6ec0b5fb87e27466abf8ee586669002b92fa82ae Mon Sep 17 00:00:00 2001 From: RHammond2 <13874373+RHammond2@users.noreply.github.com> Date: Fri, 14 Jul 2023 14:54:02 -0700 Subject: [PATCH 19/20] copy.deepcopy is not a valid converter, Chris --- floris/simulation/solver.py | 10 +++------- 1 file changed, 3 insertions(+), 7 deletions(-) diff --git a/floris/simulation/solver.py b/floris/simulation/solver.py index e8e0e4260..1c99edf97 100644 --- a/floris/simulation/solver.py +++ b/floris/simulation/solver.py @@ -85,14 +85,13 @@ def calculate_area_overlap(wake_velocities, freestream_velocities, y_ngrid, z_ng @define(auto_attribs=True) class Solver: - farm: Farm = field(converter=copy.deepcopy, validator=attrs.validators.instance_of(Farm)) - flow_field: FlowField = field(converter=copy.deepcopy, validator=attrs.validators.instance_of(FlowField)) + farm: Farm = field(validator=attrs.validators.instance_of(Farm)) + flow_field: FlowField = field(validator=attrs.validators.instance_of(FlowField)) grid: TurbineGrid | FlowFieldGrid = field( - converter=copy.deepcopy, validator=attrs.validators.instance_of((FlowFieldGrid, TurbineGrid)), ) model_manager: WakeModelManager = field( - converter=copy.deepcopy, validator=attrs.validators.instance_of(WakeModelManager) + validator=attrs.validators.instance_of(WakeModelManager) ) @abstractmethod @@ -923,9 +922,6 @@ def solve( if grid is None: grid = self.grid - gch_gain = 1.0 - scale_factor = 2.0 - # <> deflection_model_args = self.model_manager.deflection_model.prepare_function(grid, flow_field) deficit_model_args = self.model_manager.velocity_model.prepare_function(grid, flow_field) From 2797f50624fc4e03a10a84fd747cbbdf243e98c7 Mon Sep 17 00:00:00 2001 From: RHammond2 <13874373+RHammond2@users.noreply.github.com> Date: Fri, 14 Jul 2023 15:20:02 -0700 Subject: [PATCH 20/20] add FlowFieldPlanarGrid as valid solver --- floris/simulation/solver.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/floris/simulation/solver.py b/floris/simulation/solver.py index 1c99edf97..bf103348c 100644 --- a/floris/simulation/solver.py +++ b/floris/simulation/solver.py @@ -87,8 +87,8 @@ def calculate_area_overlap(wake_velocities, freestream_velocities, y_ngrid, z_ng class Solver: farm: Farm = field(validator=attrs.validators.instance_of(Farm)) flow_field: FlowField = field(validator=attrs.validators.instance_of(FlowField)) - grid: TurbineGrid | FlowFieldGrid = field( - validator=attrs.validators.instance_of((FlowFieldGrid, TurbineGrid)), + grid: TurbineGrid | FlowFieldGrid | FlowFieldPlanarGrid = field( + validator=attrs.validators.instance_of((FlowFieldGrid, TurbineGrid, FlowFieldPlanarGrid)), ) model_manager: WakeModelManager = field( validator=attrs.validators.instance_of(WakeModelManager)