Feature: farm orientation Fix(north up) and matching plot turbine with fi

examples/vis_03_horizontal_plane.py

@@ -0,0 +1,72 @@
+# 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("inputs/gch.yaml")
+
+# # 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])
+if 0 : turbine_type = ['nrel_5MW', 'nrel_5MW', 'nrel_5MW', 'nrel_5MW'] # same WTGs
+if 1 : turbine_type = ['nrel_5MW', 'nrel_5MW', 'iea_10MW', 'iea_15MW'] # mix WTGs
+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_names_sorted
+
+        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(); 

floris/simulation/grid.py

@@ -354,17 +354,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(
@@ -378,6 +391,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, :, :, :]

floris/tools/floris_interface.py

@@ -31,6 +31,8 @@
 from floris.tools.cut_plane import CutPlane
 from floris.type_dec import NDArrayFloat
 
+#dh
+from floris.utilities import rotate_coordinates_rel_west
 
 class FlorisInterface(LoggerBase):
     """
@@ -297,6 +299,14 @@ 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)
+        fix_orientation = True #dh. TODO add variable in class
+        if fix_orientation :
+            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(
@@ -342,7 +352,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
 

floris/tools/visualization.py

@@ -56,6 +56,9 @@ def plot_turbines(
     if color is None:
         color = "k"
 
+    fix_orientation=1 #dh. how to 
+    if fix_orientation : wind_direction = np.ones_like(wind_direction)*270 #dh
+
     coordinates_array = np.array([[x, y, 0.0] for x, y in list(zip(layout_x, layout_y))])
     layout_x, layout_y, _ = rotate_coordinates_rel_west(
         np.array([wind_direction]),
@@ -86,12 +89,20 @@ def plot_turbines_with_fi(fi: FlorisInterface, ax=None, color=None, yaw_angles=N
         fig, ax = plt.subplots()
     if yaw_angles is None:
         yaw_angles = fi.floris.farm.yaw_angles
+
+    #dh. after cal_wake. rotor_diameters shape has changed.   
+    try:
+        np.shape(fi.floris.farm.rotor_diameters)[1]
+        rd=fi.floris.farm.rotor_diameters[0,0] # after calculate_wake with wd or ws is over 2
+    except:
+        rd=fi.floris.farm.rotor_diameters # after FI, reinitialize, etc...
+
     plot_turbines(
         ax,
         fi.layout_x,
         fi.layout_y,
         yaw_angles[0, 0],
-        fi.floris.farm.rotor_diameters[0, 0],
+        rotor_diameters=rd, #dh
         color=color,
         wind_direction=fi.floris.flow_field.wind_directions[0],
     )

floris/utilities.py

@@ -195,13 +195,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