Skip to content

Commit

Permalink
enhancement: fix_farm_orientation for horizontal_plane
Browse files Browse the repository at this point in the history
Results plot with original layout (not rotate to west)
  • Loading branch information
dhcho347 committed Nov 14, 2022
1 parent 0c2adf3 commit ed93d54
Show file tree
Hide file tree
Showing 4 changed files with 108 additions and 6 deletions.
71 changes: 71 additions & 0 deletions examples/vis_03_horizontal_plane.txt
Original file line number Diff line number Diff line change
@@ -0,0 +1,71 @@
# 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

from floris.tools import FlorisInterface
from floris.tools.visualization import visualize_cut_plane
from floris.tools.visualization import plot_rotor_values

#dh
import numpy as np

fi = FlorisInterface("floris/examples/inputs/gch.yaml") # 3.2.1.2.1.1

# # Define 4 turbines
layout_x = np.array([3000.0, 0.0, 1500.0, 3000.0])
layout_y = np.array([800.0, 800.0, 800.0, 0.0])
turbine_type = ['nrel_5MW', 'nrel_5MW', 'iea_10MW', 'iea_15MW']
fi.reinitialize(layout_x=layout_x, layout_y=layout_y, turbine_type=turbine_type)

# sweep_wind_directions
# just with yawangle
if 1 : # just with yaw angle text
# select wind directions and wind speed for horizontal plot
wd = [[i] for i in np.arange(0,360,90)]; #change : wind directions
ws = [8.0]
# yaw angles: Change the yaw angles and configure the plot differently
n_wd=fi.floris.flow_field.n_wind_directions; n_ws=fi.floris.flow_field.n_wind_speeds
n_wtg=fi.floris.farm.n_turbines
yaw_angles = np.zeros((1, 1, n_wtg));
yaw_angles[:,:,:] = (0, 0, 15, -15)

# ready for plot
n_col=2 #change : graph's column
fig, ax_list = plt.subplots( round(len(wd)/n_col+0.5), n_col, figsize=(16, 8))
ax_list = ax_list.flatten()

horizontal_plane =[]; res=200;
# get DFs (x,y,z, u,v,w) for horizontal plane
for i in range(len(wd)):
horizontal_plane.append(fi.calculate_horizontal_plane(wd=wd[i], ws=ws, height=90.0, yaw_angles=yaw_angles, x_resolution=res, y_resolution=res), )

# plot DFs
for i in range(len(wd)):
ax=ax_list[i];
visualize_cut_plane(horizontal_plane[i], ax=ax, title="Wind direction "+str(wd[i])+"deg", color_bar=True);

# text on WTGs
turbine_yaw = yaw_angles.flatten()
turbine_type= fi.floris.farm.turbine_type

mytext = [f"yaw: {i:.1f}" for i in turbine_yaw]
if 1: mytext = [f"T{i:0d}: {turbine_type[i]} \n {mytext[i]}" for i in range(n_wtg)]
for j in range(fi.floris.farm.n_turbines):
ax.text(fi.layout_x[j], fi.layout_y[j], mytext[j], color='springgreen')

# text on Farm
ax.text(min(horizontal_plane[i].df.x1), min(horizontal_plane[i].df.x2), f' WD: {wd[i]}, WS: {ws[0]} \n',color='white')

plt.tight_layout(); plt.savefig("fix_orient.png"); plt.show();
23 changes: 20 additions & 3 deletions floris/simulation/grid.py
Original file line number Diff line number Diff line change
Expand Up @@ -330,17 +330,30 @@ def set_grid(self) -> None:
First, sort the turbines so that we know the bounds in the correct orientation.
Then, create the grid based on this wind-from-left orientation
"""

fix_orientation = True #dh. TODO add variable in class
if fix_orientation : wd = np.ones_like(self.wind_directions)*270 #dh. do not rotate
else : wd = self.wind_directions #dh. do rotate

# 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 = rotate_coordinates_rel_west(wd, self.turbine_coordinates_array)

max_diameter = np.max(self.reference_turbine_diameter)

if self.normal_vector == "z": # Rules of thumb for horizontal plane
if self.x1_bounds is None:
self.x1_bounds = (np.min(x) - 2 * max_diameter, np.max(x) + 10 * max_diameter)
# dh. broaden the flowfiled_planar
if fix_orientation :
self.x1_bounds = (np.min(x) - 10 * max_diameter, np.max(x) + 10 * max_diameter)
else :
self.x1_bounds = (np.min(x) - 2 * max_diameter, np.max(x) + 10 * max_diameter)

if self.x2_bounds is None:
self.x2_bounds = (np.min(y) - 2 * max_diameter, np.max(y) + 2 * max_diameter)
# dh
if fix_orientation :
self.x2_bounds = (np.min(y) - 10 * max_diameter, np.max(y) + 10 * max_diameter)
else :
self.x2_bounds = (np.min(y) - 2 * max_diameter, np.max(y) + 2 * max_diameter)

# TODO figure out proper z spacing for GCH, currently set to +/- 10.0
x_points, y_points, z_points = np.meshgrid(
Expand All @@ -350,6 +363,10 @@ def set_grid(self) -> None:
indexing="ij"
)

#dh. rotating meshgrid to west
if fix_orientation :
x_points, y_points, z_points = rotate_coordinates_rel_west(self.wind_directions, (x_points, y_points, z_points), inv_rot=False )

self.x_sorted = x_points[None, None, :, :, :]
self.y_sorted = y_points[None, None, :, :, :]
self.z_sorted = z_points[None, None, :, :, :]
Expand Down
11 changes: 10 additions & 1 deletion floris/tools/floris_interface.py
Original file line number Diff line number Diff line change
Expand Up @@ -29,6 +29,8 @@
from floris.tools.cut_plane import CutPlane
from floris.simulation.turbine import Ct, power, axial_induction, average_velocity

#dh
from floris.utilities import rotate_coordinates_rel_west

class FlorisInterface(LoggerBase):
"""
Expand Down Expand Up @@ -267,6 +269,13 @@ def get_plane_of_points(
v_flat = self.floris.flow_field.v_sorted[0, 0].flatten()
w_flat = self.floris.flow_field.w_sorted[0, 0].flatten()

#dh. inverse rotation of cal. results (x,y,z)
if 1 :
x_flat2, y_flat2, z_flat2 = rotate_coordinates_rel_west(self.floris.flow_field.wind_directions, (x_flat, y_flat, z_flat), inv_rot=True )
x_flat=x_flat2[0,0].flatten();
y_flat=y_flat2[0,0].flatten();
z_flat=z_flat2[0,0].flatten();

# Create a df of these
if normal_vector == "z":
df = pd.DataFrame(
Expand Down Expand Up @@ -312,7 +321,7 @@ def get_plane_of_points(
df = df.drop_duplicates()

# Sort values of df to make sure plotting is acceptable
df = df.sort_values(["x2", "x1"]).reset_index(drop=True)
#df = df.sort_values(["x2", "x1"]).reset_index(drop=True) #dh. deactivate

return df

Expand Down
9 changes: 7 additions & 2 deletions floris/utilities.py
Original file line number Diff line number Diff line change
Expand Up @@ -193,13 +193,18 @@ def wind_delta(wind_directions):
return ((wind_directions - 270) % 360 + 360) % 360


def rotate_coordinates_rel_west(wind_directions, coordinates):
def rotate_coordinates_rel_west(wind_directions, coordinates, inv_rot=False): #dh. add inv_rot
# Calculate the difference in given wind direction from 270 / West
wind_deviation_from_west = wind_delta(wind_directions)
wind_deviation_from_west = np.reshape(wind_deviation_from_west, (len(wind_directions), 1, 1))

if inv_rot : wind_deviation_from_west = -wind_deviation_from_west #dh. for inverse rotation

# Construct the arrays storing the turbine locations
x_coordinates, y_coordinates, z_coordinates = coordinates.T
if isinstance(coordinates, (np.ndarray, np.generic) ) :
x_coordinates, y_coordinates, z_coordinates = coordinates.T
else :
x_coordinates, y_coordinates, z_coordinates = coordinates #dh. to handle additional input type (mesh grid)

# 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
Expand Down

0 comments on commit ed93d54

Please sign in to comment.