diff --git a/examples/02a_visualizations_small.py b/examples/02a_visualizations_small.py new file mode 100644 index 000000000..82c36131e --- /dev/null +++ b/examples/02a_visualizations_small.py @@ -0,0 +1,78 @@ +# Copyright 2021 NREL + +# Licensed under the Apache License, Version 2.0 (the "License"); you may not +# use this file except in compliance with the License. You may obtain a copy of +# the License at http://www.apache.org/licenses/LICENSE-2.0 + +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, WITHOUT +# WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the +# License for the specific language governing permissions and limitations under +# the License. + +# See https://floris.readthedocs.io for documentation + + +import matplotlib.pyplot as plt +import numpy as np + +import floris.tools.visualization as wakeviz +from floris.tools import FlorisInterface +from floris.tools.visualization import ( + calculate_horizontal_plane_with_turbines, + visualize_cut_plane, +) + + +""" +This example initializes the FLORIS software, and then uses internal +functions to run a simulation and plot the results. In this case, +we are plotting three slices of the resulting flow field: +1. Horizontal slice parallel to the ground and located at the hub height +2. Vertical slice of parallel with the direction of the wind +3. Veritical slice parallel to to the turbine disc plane + +Additionally, an alternative method of plotting a horizontal slice +is shown. Rather than calculating points in the domain behind a turbine, +this method adds an additional turbine to the farm and moves it to +locations throughout the farm while calculating the velocity at it's +rotor. +""" + +# Initialize FLORIS with the given input file via FlorisInterface. +# For basic usage, FlorisInterface provides a simplified and expressive +# entry point to the simulation routines. +fi = FlorisInterface("inputs/gch.yaml") + +# The rotor plots show what is happening at each turbine, but we do not +# see what is happening between each turbine. For this, we use a +# grid that has points regularly distributed throughout the fluid domain. +# The FlorisInterface contains functions for configuring the new grid, +# running the simulation, and generating plots of 2D slices of the +# flow field. + +# Note this visualization grid created within the calculate_horizontal_plane function will be reset +# to what existed previously at the end of the function + +fi.reinitialize(wind_directions=[270]) + +# fi.calculate_wake() + +# Using the FlorisInterface functions, get 2D slices. +horizontal_plane = fi.calculate_horizontal_plane( + x_resolution=200, + y_resolution=100, + height=90.0, + yaw_angles=np.array([[[0.,0.,0.]]]), +) + +# Create the plots +fig, ax = plt.subplots(1, 1, figsize=(10, 8)) + + + +wakeviz.visualize_cut_plane(horizontal_plane, ax=ax, title="Horizontal") + +print(fi.get_turbine_powers()) + +plt.show() diff --git a/examples/02b_test_added_points.py b/examples/02b_test_added_points.py new file mode 100644 index 000000000..33d0369df --- /dev/null +++ b/examples/02b_test_added_points.py @@ -0,0 +1,104 @@ +# Copyright 2021 NREL + +# Licensed under the Apache License, Version 2.0 (the "License"); you may not +# use this file except in compliance with the License. You may obtain a copy of +# the License at http://www.apache.org/licenses/LICENSE-2.0 + +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, WITHOUT +# WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the +# License for the specific language governing permissions and limitations under +# the License. + +# See https://floris.readthedocs.io for documentation + + +import matplotlib.pyplot as plt +import numpy as np + +import floris.tools.visualization as wakeviz +from floris.tools import FlorisInterface +from floris.tools.visualization import ( + calculate_horizontal_plane_with_turbines, + visualize_cut_plane, +) + + +""" +Test the use of add points +""" + +# Initialize FLORIS with the given input file via FlorisInterface. +# For basic usage, FlorisInterface provides a simplified and expressive +# entry point to the simulation routines. +fi = FlorisInterface("inputs/gch.yaml") + + +# Set up a two-turbine farm +D = 126 +fi.reinitialize(layout_x=[0, 7*D], layout_y=[0, 0]) + +# Grab the power before we change anything +fi.calculate_wake() +turbine_powers_before_added_points = np.array(fi.get_turbine_powers())/1000. + +# Set up test points to run through the turbines +points_x = np.arange(D*-1, D*30, D/4) +points_y = np.zeros_like(points_x) +points_z = 90 * np.ones_like(points_x) + +# Collect the points +u_at_points = fi.sample_flow_at_points(points_x = points_x, + points_y = points_y, + points_z = points_z) + +# Re-collect the turbine powers +turbine_powers_after_added_points = np.array(fi.get_turbine_powers())/1000. + +# Reinitialize and collect the turbine powers a third time +fi.calculate_wake() +turbine_powers_after_calc_again_points = np.array(fi.get_turbine_powers())/1000. + +print('Turbine power before points: ', turbine_powers_before_added_points) +print('Turbine power after points: ', turbine_powers_after_added_points) +print('Turbine power after points and calc again: ', turbine_powers_after_calc_again_points) + +# Get the same points using the turbine-based method +cut_plane_from_turbine_sweep = calculate_horizontal_plane_with_turbines(fi, + x_resolution=100, + y_resolution=1, + x_bounds=[-D, 30*D], + y_bounds=[0, 0],) + +# Get the same points using horizontal plane +horizontal_plane = fi.calculate_horizontal_plane( + x_resolution=200, + y_resolution=1, + x_bounds=[-D, 30*D], + y_bounds=[0, 0], + height=90.0, + yaw_angles=np.array([[[0.,0.,0.]]]), +) + + +# Plot the points +fig, ax = plt.subplots() +ax.plot(cut_plane_from_turbine_sweep.df.x1/D, + cut_plane_from_turbine_sweep.df.u, + label='Velocity using turbine method', + color='k') +ax.plot(horizontal_plane.df.x1/D, + horizontal_plane.df.u, + label='Velocity using horizontal plane', + ls= '--') +ax.plot(points_x/D, + u_at_points, + label='Velocity at points', + ls= '--', + color='r') +ax.set_xlabel('X/D') +ax.set_ylabel('U (m/s)') +ax.grid() +ax.legend() + +plt.show() diff --git a/examples/02c_test_added_points_wd.py b/examples/02c_test_added_points_wd.py new file mode 100644 index 000000000..79bd27a62 --- /dev/null +++ b/examples/02c_test_added_points_wd.py @@ -0,0 +1,90 @@ +# Copyright 2021 NREL + +# Licensed under the Apache License, Version 2.0 (the "License"); you may not +# use this file except in compliance with the License. You may obtain a copy of +# the License at http://www.apache.org/licenses/LICENSE-2.0 + +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, WITHOUT +# WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the +# License for the specific language governing permissions and limitations under +# the License. + +# See https://floris.readthedocs.io for documentation + + +import matplotlib.pyplot as plt +import numpy as np + +import floris.tools.visualization as wakeviz +from floris.tools import FlorisInterface +from floris.tools.visualization import ( + calculate_horizontal_plane_with_turbines, + visualize_cut_plane, +) + + +""" +Test the use of add points +""" + +# Initialize FLORIS with the given input file via FlorisInterface. +# For basic usage, FlorisInterface provides a simplified and expressive +# entry point to the simulation routines. +fi = FlorisInterface("inputs/gch.yaml") + + +# Set up a two-turbine farm +D = 126 +fi.reinitialize(layout_x=[0, 3*D], layout_y=[0, 3*D]) + +fig, ax = plt.subplots(1,2) +fig.set_size_inches(10,4) +ax[0].scatter(fi.layout_x, fi.layout_y, color="black", label="turbine") + +# Set the wind direction to run 360 degrees +wd_array = np.arange(0, 360, 1) +fi.reinitialize(wind_directions=wd_array) + +# Grab the power before we change anything +fi.calculate_wake() + +# Simulate a met mast in between the turbines +met_mast_option = 0 # or try 1, 2, 3 (requires python >= 3.10) +match met_mast_option: + case 0: + points_x = [3*D]*4 + points_y = [0]*4 + case 1: + points_x = [200.]*4 + points_y = [200.]*4 + case 2: + points_x = [20.]*4 + points_y = [20.]*4 + case 3: + points_x = [305.]*4 + points_y = [158.]*4 + +points_z = [30, 90, 150, 250] + +ax[0].scatter(points_x, points_y, color="red", marker="x", label="test point") +ax[0].grid() +ax[0].set_xlabel("x [m]") +ax[0].set_ylabel("y [m]") +ax[0].legend() + +# Collect the points +u_at_points = fi.sample_flow_at_points(points_x = points_x, + points_y = points_y, + points_z = points_z) + + +# Plot the velocities +for z_idx, z in enumerate(points_z): + ax[1].plot(wd_array, u_at_points[:, z_idx].flatten(), label=f'Speed at {z} m') +ax[1].grid() +ax[1].legend() +ax[1].set_xlabel('Wind Direction (deg)') +ax[1].set_ylabel('Wind Speed (m/s)') + +plt.show() diff --git a/floris/simulation/__init__.py b/floris/simulation/__init__.py index 2ae2a870d..02fa29b5f 100644 --- a/floris/simulation/__init__.py +++ b/floris/simulation/__init__.py @@ -39,7 +39,7 @@ from .base import BaseClass, BaseModel, State from .turbine import average_velocity, axial_induction, Ct, power, rotor_effective_velocity, Turbine from .farm import Farm -from .grid import FlowFieldGrid, FlowFieldPlanarGrid, Grid, TurbineGrid +from .grid import FlowFieldGrid, FlowFieldPlanarGrid, Grid, PointsGrid, TurbineGrid from .flow_field import FlowField from .wake import WakeModelManager from .solver import ( diff --git a/floris/simulation/floris.py b/floris/simulation/floris.py index c64fed99d..5d4721609 100644 --- a/floris/simulation/floris.py +++ b/floris/simulation/floris.py @@ -14,12 +14,13 @@ from __future__ import annotations +import copy from pathlib import Path import yaml from attrs import define, field -import floris.logging_manager as logging_manager +from floris import logging_manager from floris.simulation import ( BaseClass, cc_solver, @@ -31,6 +32,7 @@ full_flow_sequential_solver, full_flow_turbopark_solver, Grid, + PointsGrid, sequential_solver, State, TurbineGrid, @@ -103,7 +105,15 @@ def __attrs_post_init__(self) -> None: time_series=self.flow_field.time_series, ) elif self.solver["type"] == "flow_field_planar_grid": - self.grid = FlowFieldPlanarGrid( + self.grid = TurbineGrid( + turbine_coordinates=self.farm.coordinates, + reference_turbine_diameter=self.farm.rotor_diameters, + wind_directions=self.flow_field.wind_directions, + wind_speeds=self.flow_field.wind_speeds, + grid_resolution=3, + time_series=self.flow_field.time_series, + ) + self.field_grid = FlowFieldPlanarGrid( turbine_coordinates=self.farm.coordinates, reference_turbine_diameter=self.farm.rotor_diameters, wind_directions=self.flow_field.wind_directions, @@ -111,9 +121,30 @@ def __attrs_post_init__(self) -> None: normal_vector=self.solver["normal_vector"], planar_coordinate=self.solver["planar_coordinate"], grid_resolution=self.solver["flow_field_grid_points"], + time_series=self.flow_field.time_series, x1_bounds=self.solver["flow_field_bounds"][0], x2_bounds=self.solver["flow_field_bounds"][1], + ) + elif self.solver["type"] == "points_grid": + self.grid = TurbineGrid( + turbine_coordinates=self.farm.coordinates, + reference_turbine_diameter=self.farm.rotor_diameters, + wind_directions=self.flow_field.wind_directions, + wind_speeds=self.flow_field.wind_speeds, + grid_resolution=3, + time_series=self.flow_field.time_series, + ) + self.field_grid = PointsGrid( + turbine_coordinates=self.farm.coordinates, + reference_turbine_diameter=self.farm.rotor_diameters, + wind_directions=self.flow_field.wind_directions, + wind_speeds=self.flow_field.wind_speeds, + grid_resolution=1, time_series=self.flow_field.time_series, + points_x=self.solver["points_x"], + points_y=self.solver["points_y"], + points_z=self.solver["points_z"], + center_of_rotation=self.grid.center_of_rotation ) else: raise ValueError( @@ -212,7 +243,118 @@ def solve_for_viz(self): elif vel_model=="turbopark": full_flow_turbopark_solver(self.farm, self.flow_field, self.grid, self.wake) else: - full_flow_sequential_solver(self.farm, self.flow_field, self.grid, self.wake) + + turbine_grid_farm = copy.deepcopy(self.farm) + turbine_grid_flow_field = copy.deepcopy(self.flow_field) + + turbine_grid_farm.construct_turbine_map() + turbine_grid_farm.construct_turbine_fCts() + turbine_grid_farm.construct_turbine_power_interps() + turbine_grid_farm.construct_hub_heights() + turbine_grid_farm.construct_rotor_diameters() + turbine_grid_farm.construct_turbine_TSRs() + turbine_grid_farm.construct_turbine_pPs() + turbine_grid_farm.construct_turbine_pTs() + turbine_grid_farm.construct_turbine_ref_density_cp_cts() + turbine_grid_farm.construct_turbine_ref_tilt_cp_cts() + turbine_grid_farm.construct_turbine_fTilts() + turbine_grid_farm.construct_turbine_correct_cp_ct_for_tilt() + turbine_grid_farm.construct_coordinates() + turbine_grid_farm.set_tilt_to_ref_tilt(self.flow_field.n_wind_directions, + self.flow_field.n_wind_speeds) + + turbine_grid_farm.expand_farm_properties( + turbine_grid_flow_field.n_wind_directions, + turbine_grid_flow_field.n_wind_speeds, + self.grid.sorted_coord_indices + ) + + turbine_grid_flow_field.initialize_velocity_field(self.grid) + turbine_grid_farm.initialize(self.grid.sorted_indices) + sequential_solver(turbine_grid_farm, self.flow_field, self.grid, self.wake) + + turbine_grid_flow_field.turbulence_intensity_field_sorted_avg = \ + self.flow_field.turbulence_intensity_field_sorted_avg + + full_flow_sequential_solver(turbine_grid_farm, + self.flow_field, + self.field_grid, + self.wake, + self.grid, + turbine_grid_flow_field + ) + + def solve_for_points(self, points_x, points_y, points_z): + # Do the calculation with the TurbineGrid for a single wind speed + # and wind direction and 1 point on the grid. Then, use the result + # to construct the full flow field grid. + # This function call should be for a single wind direction and wind speed + # since the memory consumption is very large. + + # self.flow_field.initialize_velocity_field(self.grid) + + 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) + elif vel_model=="turbopark": + full_flow_turbopark_solver(self.farm, self.flow_field, self.grid, self.wake) + else: + + # Instantiate the flow_grid + field_grid = PointsGrid( + turbine_coordinates=self.farm.coordinates, + reference_turbine_diameter=self.farm.rotor_diameters, + wind_directions=self.flow_field.wind_directions, + wind_speeds=self.flow_field.wind_speeds, + grid_resolution=1, + time_series=self.flow_field.time_series, + points_x=points_x, + points_y=points_y, + points_z=points_z, + center_of_rotation=self.grid.center_of_rotation + ) + + # turbine_grid_farm = copy.deepcopy(self.farm) + # turbine_grid_flow_field = copy.deepcopy(self.flow_field) + + # turbine_grid_farm.construct_turbine_map() + # turbine_grid_farm.construct_turbine_fCts() + # turbine_grid_farm.construct_turbine_power_interps() + # turbine_grid_farm.construct_hub_heights() + # turbine_grid_farm.construct_rotor_diameters() + # turbine_grid_farm.construct_turbine_TSRs() + # turbine_grid_farm.construct_turbine_pPs() + # turbine_grid_farm.construct_turbine_pTs() + # turbine_grid_farm.construct_turbine_ref_density_cp_cts() + # turbine_grid_farm.construct_turbine_ref_tilt_cp_cts() + # turbine_grid_farm.construct_turbine_fTilts() + # turbine_grid_farm.construct_turbine_correct_cp_ct_for_tilt() + # turbine_grid_farm.construct_coordinates() + # turbine_grid_farm.set_tilt_to_ref_tilt(self.flow_field.n_wind_directions, + # self.flow_field.n_wind_speeds) + + # turbine_grid_farm.expand_farm_properties( + # turbine_grid_flow_field.n_wind_directions, + # turbine_grid_flow_field.n_wind_speeds, + # self.grid.sorted_coord_indices + # ) + + # DON"T RERUN sequential_solver + #sequential_solver(turbine_grid_farm, self.flow_field, self.grid, self.wake) + + turbine_grid_flow_field = copy.deepcopy(self.flow_field) + self.flow_field.initialize_velocity_field(field_grid) + + full_flow_sequential_solver(self.farm, + self.flow_field, + field_grid, + self.wake, + self.grid, + turbine_grid_flow_field) + + return self.flow_field.u_sorted.squeeze() #TODO: This returns an + # indeterminate number of dimensions, maybe bad def finalize(self): # Once the wake calculation is finished, unsort the values to match diff --git a/floris/simulation/grid.py b/floris/simulation/grid.py index 7bc6cf94a..70d9bf367 100644 --- a/floris/simulation/grid.py +++ b/floris/simulation/grid.py @@ -121,6 +121,8 @@ def grid_resolution_validator(self, instance: attrs.Attribute, value: int | Iter assert type(value[0]) is int assert type(value[1]) is int assert type(value[2]) is int + elif type(self) is PointsGrid: + return else: raise TypeError("`grid_resolution` must be of type int or Iterable(int,)") @@ -147,6 +149,7 @@ class TurbineGrid(Grid): sorted_indices: NDArrayInt = field(init=False) sorted_coord_indices: NDArrayInt = field(init=False) unsorted_indices: NDArrayInt = field(init=False) + center_of_rotation: NDArrayFloat = field(init=False) def __attrs_post_init__(self) -> None: super().__attrs_post_init__() @@ -206,7 +209,12 @@ def set_grid(self) -> None: # the foot of the turbine where the tower meets the ground. # These are the rotated coordinates of the wind turbines based on the wind direction - x, y, z = rotate_coordinates_rel_west(self.wind_directions, self.turbine_coordinates_array) + x, y, z, cor = rotate_coordinates_rel_west( + self.wind_directions, + self.turbine_coordinates_array, + return_center_of_rotation=True + ) + self.center_of_rotation = cor # - **rloc** (*float, optional): A value, from 0 to 1, that determines # the width/height of the grid of points on the rotor as a ratio of @@ -473,3 +481,36 @@ def set_grid(self) -> None: # xoffset * sind(angle) + yoffset * cosd(angle) + center_of_rotation[1] # ) # return rotated_x, rotated_y, self.z +@define +class PointsGrid(Grid): + """ + Args: + grid_resolution (`Vec3`): The number of grid points to be created in each direction. + turbine_coordinates (`list[Vec3]`): The collection of turbine coordinate (`Vec3`) objects. + reference_turbine_diameter (:py:obj:`float`): The reference turbine's rotor diameter. + grid_resolution (:py:obj:`int`): The number of points on each turbine + """ + points_x: NDArrayFloat = field() + points_y: NDArrayFloat = field() + points_z: NDArrayFloat = field() + center_of_rotation: NDArrayFloat | None = field(default=None) + + def __attrs_post_init__(self) -> None: + super().__attrs_post_init__() + self.set_grid() + + def set_grid(self) -> None: + """ + Set points for calculation based on a series of user-supplied coordinates. + """ + point_coordinates = np.array(list(zip(self.points_x, self.points_y, self.points_z))) + + # These are the rotated coordinates of the wind turbines based on the wind direction + x, y, z = rotate_coordinates_rel_west( + self.wind_directions, + point_coordinates, + self.center_of_rotation + ) + self.x_sorted = x[:,:,:,None,None] + self.y_sorted = y[:,:,:,None,None] + self.z_sorted = z[:,:,:,None,None] diff --git a/floris/simulation/solver.py b/floris/simulation/solver.py index f1c8e636e..fd2bf4256 100644 --- a/floris/simulation/solver.py +++ b/floris/simulation/solver.py @@ -10,9 +10,9 @@ # License for the specific language governing permissions and limitations under # the License. +from __future__ import annotations + import copy -import sys -import time import numpy as np @@ -22,6 +22,8 @@ Farm, FlowField, FlowFieldGrid, + FlowFieldPlanarGrid, + PointsGrid, TurbineGrid, ) from floris.simulation.turbine import average_velocity @@ -244,50 +246,18 @@ def sequential_solver( axis=(3,4) )[:, :, :, None, None] - def full_flow_sequential_solver( farm: Farm, flow_field: FlowField, - flow_field_grid: FlowFieldGrid, - model_manager: WakeModelManager + flow_field_grid: FlowFieldGrid | FlowFieldPlanarGrid | PointsGrid, + model_manager: WakeModelManager, + turbine_grid: TurbineGrid, + turbine_grid_flow_field: FlowField ) -> None: # Get the flow quantities and turbine performance turbine_grid_farm = copy.deepcopy(farm) - turbine_grid_flow_field = copy.deepcopy(flow_field) - - turbine_grid_farm.construct_turbine_map() - turbine_grid_farm.construct_turbine_fCts() - turbine_grid_farm.construct_turbine_power_interps() - turbine_grid_farm.construct_hub_heights() - turbine_grid_farm.construct_rotor_diameters() - turbine_grid_farm.construct_turbine_TSRs() - turbine_grid_farm.construct_turbine_pPs() - turbine_grid_farm.construct_turbine_pTs() - turbine_grid_farm.construct_turbine_ref_density_cp_cts() - turbine_grid_farm.construct_turbine_ref_tilt_cp_cts() - turbine_grid_farm.construct_turbine_fTilts() - turbine_grid_farm.construct_turbine_correct_cp_ct_for_tilt() - turbine_grid_farm.construct_coordinates() - turbine_grid_farm.set_tilt_to_ref_tilt(flow_field.n_wind_directions, flow_field.n_wind_speeds) - - - turbine_grid = TurbineGrid( - turbine_coordinates=turbine_grid_farm.coordinates, - reference_turbine_diameter=turbine_grid_farm.rotor_diameters, - wind_directions=turbine_grid_flow_field.wind_directions, - wind_speeds=turbine_grid_flow_field.wind_speeds, - grid_resolution=3, - time_series=turbine_grid_flow_field.time_series, - ) - turbine_grid_farm.expand_farm_properties( - turbine_grid_flow_field.n_wind_directions, - turbine_grid_flow_field.n_wind_speeds, - turbine_grid.sorted_coord_indices - ) - turbine_grid_flow_field.initialize_velocity_field(turbine_grid) - turbine_grid_farm.initialize(turbine_grid.sorted_indices) - sequential_solver(turbine_grid_farm, turbine_grid_flow_field, turbine_grid, model_manager) + # turbine_grid_flow_field = copy.deepcopy(flow_field) ### Referring to the quantities from above, calculate the wake in the full grid @@ -387,6 +357,7 @@ def full_flow_sequential_solver( if model_manager.enable_transverse_velocities: v_wake, w_wake = calculate_transverse_velocity( u_i, + # turbine_grid_flow_field.u_initial_sorted, flow_field.u_initial_sorted, flow_field.dudz_initial_sorted, flow_field_grid.x_sorted - x_i, diff --git a/floris/tools/floris_interface.py b/floris/tools/floris_interface.py index d093b4165..39495ad1f 100644 --- a/floris/tools/floris_interface.py +++ b/floris/tools/floris_interface.py @@ -347,6 +347,7 @@ def get_plane_of_points( return df + # New version using points def calculate_horizontal_plane( self, height, @@ -378,59 +379,183 @@ def calculate_horizontal_plane( :py:class:`~.tools.cut_plane.CutPlane`: containing values of x, y, u, v, w """ - # TODO update docstring + if wd is None: wd = self.floris.flow_field.wind_directions if ws is None: ws = self.floris.flow_field.wind_speeds self.check_wind_condition_for_viz(wd=wd, ws=ws) + self.reinitialize(wind_directions=wd, wind_speeds=ws) - # Store the current state for reinitialization - floris_dict = self.floris.as_dict() - current_yaw_angles = self.floris.farm.yaw_angles + # Use the provided yaw angles to calculate wake + #TODO: Tilt? + self.calculate_wake(yaw_angles=yaw_angles) - # Set the solver to a flow field planar grid - solver_settings = { - "type": "flow_field_planar_grid", - "normal_vector": "z", - "planar_coordinate": height, - "flow_field_grid_points": [x_resolution, y_resolution], - "flow_field_bounds": [x_bounds, y_bounds], - } - self.reinitialize(wind_directions=wd, wind_speeds=ws, solver_settings=solver_settings) - # TODO this has to be done here as it seems to be lost with reinitialize - if yaw_angles is not None: - self.floris.farm.yaw_angles = yaw_angles + # Get a grid of points test test + D = self.floris.farm.rotor_diameters_sorted[0, 0, 0] + if x_bounds is None: + x_bounds = (np.min(self.layout_x) - 2 * D, np.max(self.layout_x) + 10 * D) - # Calculate wake - self.floris.solve_for_viz() + if y_bounds is None: + y_bounds = (np.min(self.layout_y) - 2 * D, np.max(self.layout_y) + 2 * D) - # Get the points of data in a dataframe - # TODO this just seems to be flattening and storing the data in a df; is this necessary? - # It seems the biggest depenedcy is on CutPlane and the subsequent visualization tools. - df = self.get_plane_of_points( - normal_vector="z", - planar_coordinate=height, - ) + # Now generate a list of points + x = np.linspace(x_bounds[0], x_bounds[1], x_resolution) + y = np.linspace(y_bounds[0], y_bounds[1], y_resolution) - # Compute the cutplane - horizontal_plane = CutPlane( - df, - self.floris.grid.grid_resolution[0], - self.floris.grid.grid_resolution[1], - "z" - ) + points_x, points_y = np.meshgrid(x, y)#, indexing='ij') + points_x = points_x.flatten() + points_y = points_y.flatten() + points_z = height * np.ones_like(points_x) - # Reset the fi object back to the turbine grid configuration - self.floris = Floris.from_dict(floris_dict) - self.floris.flow_field.het_map = self.het_map + u = self.floris.solve_for_points(points_x, points_y, points_z) - # Run the simulation again for futher postprocessing (i.e. now we can get farm power) - self.calculate_wake(yaw_angles=current_yaw_angles) + # Make a dataframe + df = pd.DataFrame({ + 'x1':points_x, + 'x2':points_y, + 'x3':points_z, + 'u':u, + 'v':0, + 'w':0, + }) + + # Convert to a cut_plane + horizontal_plane = CutPlane(df, x_resolution, y_resolution, "z") return horizontal_plane + + + # # OLD VERSION + # def calculate_horizontal_plane( + # self, + # height, + # x_resolution=200, + # y_resolution=200, + # x_bounds=None, + # y_bounds=None, + # wd=None, + # ws=None, + # yaw_angles=None, + # ): + # """ + # Shortcut method to instantiate a :py:class:`~.tools.cut_plane.CutPlane` + # object containing the velocity field in a horizontal plane cut through + # the simulation domain at a specific height. + + # Args: + # height (float): Height of cut plane. Defaults to Hub-height. + # x_resolution (float, optional): Output array resolution. + # Defaults to 200 points. + # y_resolution (float, optional): Output array resolution. + # Defaults to 200 points. + # x_bounds (tuple, optional): Limits of output array (in m). + # Defaults to None. + # y_bounds (tuple, optional): Limits of output array (in m). + # Defaults to None. + + # Returns: + # :py:class:`~.tools.cut_plane.CutPlane`: containing values + # of x, y, u, v, w + # """ + # # TODO update docstring + # if wd is None: + # wd = self.floris.flow_field.wind_directions + # if ws is None: + # ws = self.floris.flow_field.wind_speeds + # self.check_wind_condition_for_viz(wd=wd, ws=ws) + + # # Store the current state for reinitialization + # floris_dict = self.floris.as_dict() + # current_yaw_angles = self.floris.farm.yaw_angles + + # # Set the solver to a flow field planar grid + # solver_settings = { + # "type": "flow_field_planar_grid", + # "normal_vector": "z", + # "planar_coordinate": height, + # "flow_field_grid_points": [x_resolution, y_resolution], + # "flow_field_bounds": [x_bounds, y_bounds], + # } + # self.reinitialize(wind_directions=wd, wind_speeds=ws, solver_settings=solver_settings) + + # # TODO this has to be done here as it seems to be lost with reinitialize + # if yaw_angles is not None: + # self.floris.farm.yaw_angles = yaw_angles + + # # Calculate wake + # self.floris.solve_for_viz() + + + # print(self.floris.flow_field.u_sorted.shape) + + # self.floris.grid = self.floris.field_grid + + # # Get the points of data in a dataframe + # # TODO this just seems to be flattening and storing the data in a df; is this necessary? + # # It seems the biggest depenedcy is on CutPlane and the subsequent visualization tools. + # df = self.get_plane_of_points( + # normal_vector="z", + # planar_coordinate=height, + # ) + + # # Compute the cutplane + # horizontal_plane = CutPlane( + # df, + # self.floris.grid.grid_resolution[0], + # self.floris.grid.grid_resolution[1], + # "z" + # ) + + # # Reset the fi object back to the turbine grid configuration + # self.floris = Floris.from_dict(floris_dict) + # self.floris.flow_field.het_map = self.het_map + + # # Run the simulation again for futher postprocessing (i.e. now we can get farm power) + # # TODO THIS SHOULD INCLUDE TILT? + # self.calculate_wake(yaw_angles=current_yaw_angles) + + # return horizontal_plane + + + def sample_flow_at_points( + self, + points_x: NDArrayFloat, + points_y: NDArrayFloat, + points_z: NDArrayFloat, + ): + """ + To be written + + Args: + probe_points (list[tuple], optional): List of probe points to add. + + Returns: + Nothing? + """ + + # Check that point_x, point_y, point_z are all the same length + if not ( + len(points_x) == len(points_y) + and len(points_x) == len(points_z) + ): + raise ValueError( + "The number of points in each coordinate direction must be the same." + ) + + # Confirm calculate wake has been run + if self.floris.state is not State.USED: + raise RuntimeError( + "Can't run function `FlorisInterface.sample_flow_at_points` without " + "first running `FlorisInterface.calculate_wake`." + ) + + # DON"T create a new FLORIS object. Use the one we have! + + return self.floris.solve_for_points(points_x, points_y, points_z) + def calculate_cross_plane( self, downstream_dist, diff --git a/floris/utilities.py b/floris/utilities.py index 733104e94..1c55d7224 100644 --- a/floris/utilities.py +++ b/floris/utilities.py @@ -220,7 +220,9 @@ def wind_delta(wind_directions: NDArrayFloat | float): def rotate_coordinates_rel_west( wind_directions: NDArrayFloat, - coordinates: NDArrayFloat + coordinates: NDArrayFloat, + center_of_rotation: NDArrayFloat | None = None, + return_center_of_rotation: bool = False ): """ This function rotates the given coordinates so that they are aligned with West (270) rather @@ -243,8 +245,16 @@ def rotate_coordinates_rel_west( x_coordinates, y_coordinates, z_coordinates = coordinates.T # Find center of rotation - this is the center of box bounding all of the turbines - x_center_of_rotation = (np.min(x_coordinates) + np.max(x_coordinates)) / 2 - y_center_of_rotation = (np.min(y_coordinates) + np.max(y_coordinates)) / 2 + if center_of_rotation is None: + x_center_of_rotation = (np.min(x_coordinates) + np.max(x_coordinates)) / 2 + y_center_of_rotation = (np.min(y_coordinates) + np.max(y_coordinates)) / 2 + else: + x_center_of_rotation = center_of_rotation[0] + y_center_of_rotation = center_of_rotation[1] + + #x_center_of_rotation = 378.0 + #y_center_of_rotation = 0 + #import ipdb; ipdb.set_trace() # Rotate turbine coordinates about the center x_coord_offset = x_coordinates - x_center_of_rotation @@ -260,7 +270,12 @@ def rotate_coordinates_rel_west( + y_center_of_rotation ) z_coord_rotated = np.ones_like(wind_deviation_from_west) * z_coordinates - return x_coord_rotated, y_coord_rotated, z_coord_rotated + + if return_center_of_rotation: + cor = np.array([x_center_of_rotation, y_center_of_rotation]) + return x_coord_rotated, y_coord_rotated, z_coord_rotated, cor + else: + return x_coord_rotated, y_coord_rotated, z_coord_rotated class Loader(yaml.SafeLoader): diff --git a/tests/conftest.py b/tests/conftest.py index 19dd3c369..852136afd 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -22,6 +22,7 @@ Floris, FlowField, FlowFieldGrid, + PointsGrid, TurbineGrid, ) from floris.utilities import Vec3 @@ -146,6 +147,25 @@ def flow_field_grid_fixture(sample_inputs_fixture) -> FlowFieldGrid: grid_resolution=[3,2,2] ) +@pytest.fixture +def points_grid_fixture(sample_inputs_fixture) -> PointsGrid: + turbine_coordinates = [Vec3(c) for c in list(zip(X_COORDS, Y_COORDS, Z_COORDS))] + rotor_diameters = ROTOR_DIAMETER * np.ones( (N_WIND_DIRECTIONS, N_WIND_SPEEDS, N_TURBINES) ) + points_x = [0.0, 10.0] + points_y = [0.0, 0.0] + points_z = [1.0, 2.0] + return PointsGrid( + turbine_coordinates=turbine_coordinates, + reference_turbine_diameter=rotor_diameters, + wind_directions=np.array(WIND_DIRECTIONS), + wind_speeds=np.array(WIND_SPEEDS), + grid_resolution=None, + time_series=False, + points_x=points_x, + points_y=points_y, + points_z=points_z, + ) + @pytest.fixture def floris_fixture(): sample_inputs = SampleInputs()