Add and rearange initial example

.gitignore

@@ -24,6 +24,8 @@ pip-wheel-metadata
 examples/cp_ct_cq_lut.p
 _site/
 .jekyll-cache/
+examples/hist.hist
+examples/SLSQP.out
 
 # Log files
-*.log
+*.log

examples/01_opening_floris_computing_power.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
+
+from floris.tools import FlorisInterface
+
+"""
+This example creates a FLORIS instance
+1) Makes a two-turbine layout
+2) Demonstrates single ws/wd simulations
+3) Demonstrates mulitple ws/wd simulations
+
+Main concept is introduce FLORIS and illustrate essential structure of most-used FLORIS calls
+"""
+
+# 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")
+
+# Convert to a simple two turbine layout
+fi.reinitialize( layout=( [0, 500.], [0., 0.] ) )
+
+# Single wind speed and wind direction
+print('\n============================= Single Wind Direction and Wind Speed =============================')
+
+# Get the turbine powers assuming 1 wind speed and 1 wind direction
+fi.reinitialize(wind_directions=[270.], wind_speeds=[8.0])
+
+# Set the yaw angles to 0
+yaw_angles = np.zeros([1,1,2]) # 1 wind direction, 1 wind speed, 2 turbines
+fi.calculate_wake(yaw_angles=yaw_angles)
+
+# Get the turbine powers
+turbine_powers = fi.get_turbine_powers()/1000.
+print('The turbine power matrix should be of dimensions 1 WD X 1 WS X 2 Turbines')
+print(turbine_powers)
+print("Shape: ",turbine_powers.shape)
+
+# Single wind speed and wind direction
+print('\n============================= Single Wind Direction and Multiple Wind Speeds =============================')
+
+
+wind_speeds = np.array([8.0, 9.0, 10.0])
+fi.reinitialize(wind_speeds=wind_speeds)
+yaw_angles = np.zeros([1,3,2]) # 1 wind direction, 3 wind speeds, 2 turbines
+fi.calculate_wake(yaw_angles=yaw_angles)
+turbine_powers = fi.get_turbine_powers()/1000.
+print('The turbine power matrix should be of dimensions 1 WD X 3 WS X 2 Turbines')
+print(turbine_powers)
+print("Shape: ",turbine_powers.shape)
+
+# Single wind speed and wind direction
+print('\n============================= Multiple Wind Directions and Multiple Wind Speeds =============================')
+
+wind_directions = np.array([260., 270., 280.])
+wind_speeds = np.array([8.0, 9.0, 10.0])
+fi.reinitialize(wind_directions=wind_directions, wind_speeds=wind_speeds)
+yaw_angles = np.zeros([1,3,2]) # 1 wind direction, 3 wind speeds, 2 turbines
+fi.calculate_wake(yaw_angles=yaw_angles)
+turbine_powers = fi.get_turbine_powers()/1000.
+print('The turbine power matrix should be of dimensions 3 WD X 3 WS X 2 Turbines')
+print(turbine_powers)
+print("Shape: ",turbine_powers.shape)

examples/01_plot_wakes.py → examples/02_visualizations.py

@@ -33,47 +33,60 @@
 # entry point to the simulation routines.
 fi = FlorisInterface("inputs/gch.yaml")
 
-# FLORIS supports multiple types of grids for capturing wind speed
-# information. The current input file is configured with a square grid
-# placed on each rotor plane with 9 points in a 3x3 layout. This can
-# be plotted to show the wind conditions that each turbine is experiencing.
-# However, 9 points is very coarse for visualization so let's first
-# increase the rotor grid to 10x10 points.
-
-# Create a solver settings dictionary with the new number of points
-solver_settings = {
-  "type": "turbine_grid",
-  "turbine_grid_points": 10
-}
-
-# Since we already have a FlorisInterface (fi), simply reinitialize it
-# with the new configuration
-fi.reinitialize(solver_settings=solver_settings)
-
-# Run the wake calculation to get the turbine-turbine interfactions
-# on the turbine grids
-fi.calculate_wake()
-
-# Plot the values at each rotor
-plot_rotor_values(fi.floris.flow_field.u, wd_index=0, ws_index=0, n_rows=1, n_cols=3)
-
 # 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
+
 # Using the FlorisInterface functions, get 2D slices.
 horizontal_plane = fi.calculate_horizontal_plane(x_resolution=200, y_resolution=100, height=90.0)
 y_plane = fi.calculate_y_plane(x_resolution=200, z_resolution=100, crossstream_dist=630.0)
 cross_plane = fi.calculate_cross_plane(y_resolution=100, z_resolution=100, downstream_dist=630.0)
 
+
 # Create the plots
 fig, ax_list = plt.subplots(3, 1, figsize=(10, 8))
 ax_list = ax_list.flatten()
 visualize_cut_plane(horizontal_plane, ax=ax_list[0], title="Horizontal")
 visualize_cut_plane(y_plane, ax=ax_list[1], title="Streamwise profile")
 visualize_cut_plane(cross_plane, ax=ax_list[2], title="Spanwise profile")
 
+# FLORIS further includes visualization methods for visualing the rotor plane of each
+# Turbine in the simulation
+
+# Run the wake calculation to get the turbine-turbine interfactions
+# on the turbine grids
+fi.calculate_wake()
+
+# Plot the values at each rotor
+fig, axes, _ , _ = plot_rotor_values(fi.floris.flow_field.u, wd_index=0, ws_index=0, n_rows=1, n_cols=3, return_fig_objects=True)
+fig.suptitle("Rotor Plane Visualization, Original Resolution")
+
+# FLORIS supports multiple types of grids for capturing wind speed
+# information. The current input file is configured with a square grid
+# placed on each rotor plane with 9 points in a 3x3 layout. For visualization,
+# this resolution can be increased.  Note this operation, unlike the 
+# calc_x_plane above operations does not automatically reset the grid to
+# the initial status as definied by the input file
+
+# Increase the resolution of points on each turbien plane
+solver_settings = {
+  "type": "turbine_grid",
+  "turbine_grid_points": 10
+}
+fi.reinitialize(solver_settings=solver_settings)
+
+# Run the wake calculation to get the turbine-turbine interfactions
+# on the turbine grids
+fi.calculate_wake()
+
+# Plot the values at each rotor
+fig, axes, _ , _ = plot_rotor_values(fi.floris.flow_field.u, wd_index=0, ws_index=0, n_rows=1, n_cols=3, return_fig_objects=True)
+fig.suptitle("Rotor Plane Visualization, 10x10 Resolution")
+
 plt.show()

examples/02_basic_adjustments.py → examples/03_making_adjustments.py

@@ -48,7 +48,7 @@
 
 # Change the wind shear, reset the wind speed, and plot a vertical slice
 fi.reinitialize( wind_shear=0.2, wind_speeds=[8.0] )
-y_plane = fi.calculate_y_plane()
+y_plane = fi.calculate_y_plane(crossstream_dist=0.0)
 visualize_cut_plane(y_plane, ax=axarr[2], title="Wind shear at 0.2", minSpeed=MIN_WS, maxSpeed=MAX_WS)
 
 
@@ -61,7 +61,6 @@
 fi.reinitialize( layout=( X.flatten(), Y.flatten() ) )
 horizontal_plane = fi.calculate_horizontal_plane(height=90.0)
 visualize_cut_plane(horizontal_plane, ax=axarr[3], title="3x3 Farm", minSpeed=MIN_WS, maxSpeed=MAX_WS)
-# plot_turbines_with_fi(axarr[7], fi)
 
 
 # Change the yaw angles and configure the plot differently
@@ -88,10 +87,4 @@
 axarr[5].invert_xaxis()
 
 
-plt.show()
-
-
-# # Change the veer
-# fi.reinitialize(wind_veer=5.0)
-# fi.floris.solve_for_viz()
-# plot_slice_shortcut(fi, axarr[5], "Veer=5")
+plt.show()

examples/04_sweep_wind_directions.py

@@ -72,4 +72,4 @@
 ax.legend()
 ax.set_xlabel('Wind Direction (deg)')
 ax.set_ylabel('Power (kW)')
-plt.show()
+plt.show()

examples/05_sweep_wind_speeds.py

@@ -72,4 +72,4 @@
 ax.legend()
 ax.set_xlabel('Wind Speed (m/s)')
 ax.set_ylabel('Power (kW)')
-plt.show()
+plt.show()

examples/06_sweep_wind_conditions.py

@@ -44,7 +44,6 @@
 layout_y = [0, 0, 0, 0, 0]
 fi.reinitialize(layout = [layout_x, layout_y])
 
-
 # Define a ws and wd to sweep
 # Note that all combinations will be computed
 ws_array = np.arange(6, 9, 1.)
@@ -65,7 +64,6 @@
 # Collect the turbine powers
 turbine_powers = fi.get_turbine_powers() / 1E3 # In kW
 
-
 # Show results by ws and wd
 fig, axarr = plt.subplots(num_ws, 1, sharex=True,sharey=True,figsize=(6,10))
 for ws_idx, ws in enumerate(ws_array):
@@ -78,12 +76,10 @@
     ax.set_ylabel('Power (kW)')
 ax.set_xlabel('Wind Direction (deg)')
 
-
 # Sum across wind speeds and directions to show energy produced by turbine as bar plot
 energy_by_turbine = np.sum(turbine_powers, axis=(0,1)) # Sum over wind direction (0-axis) and wind speed (1-axis)
 fig, ax = plt.subplots()
 ax.bar(['T%d' % t for t in range(num_turbine)],energy_by_turbine)
 ax.set_title('Energy Produced by Turbine')
 
 plt.show()
-

examples/07_calc_aep_from_rose.py

@@ -12,6 +12,7 @@
 
 # See https://floris.readthedocs.io for documentation
 
+
 import numpy as np
 import pandas as pd
 from scipy.interpolate import NearestNDInterpolator
@@ -29,11 +30,19 @@
 # Read the windrose information file & normalize wind rose frequencies
 fn = "inputs/wind_rose.csv"
 df_wr = pd.read_csv(fn)
-df_wr["freq_val"] = df_wr["freq_val"] / df_wr["freq_val"].sum()
+
+# Normalize the frequencies
+df_wr["freq_val"] = df_wr["freq_val"].copy() / df_wr["freq_val"].sum()
+
+# Split the wind rose into wind speeds above (df_wr_op)
+# and below cut-in (df_wr_below_cut_in) as below cut-in winds are 0 power and
+# 0 ambient wind speed generates warnings and nans in some of the wake models
+df_wr_below_cut_in = df_wr[df_wr.ws < 3.0].copy()
+df_wr_op = df_wr[df_wr.ws >= 3.0].copy()
 
 # Derive the wind directions and speeds we need to evaluate in FLORIS
-wd_array = np.array(df_wr["wd"].unique(), dtype=float)
-ws_array = np.array(df_wr["ws"].unique(), dtype=float)
+wd_array = np.array(df_wr_op["wd"].unique(), dtype=float)
+ws_array = np.array(df_wr_op["ws"].unique(), dtype=float)
 
 # Load the default example FLORIS object
 fi = FlorisInterface("inputs/gch.yaml") # GCH model matched to the default "legacy_gauss" of V2
@@ -54,18 +63,23 @@
 farm_power_array = fi.get_farm_power()
 
 # Now map the FLORIS solutions to the wind rose dataframe
-X, Y = np.meshgrid(wd_array, ws_array, indexing='ij')
+wd_grid, ws_grid = np.meshgrid(wd_array, ws_array, indexing='ij')
 interpolant = NearestNDInterpolator(
-    np.vstack([X.flatten(), Y.flatten()]).T, 
+    np.vstack([wd_grid.flatten(), ws_grid.flatten()]).T, 
     farm_power_array.flatten()
 )
-df_wr["farm_power"] = interpolant(df_wr[["wd", "ws"]])
-df_wr["farm_power"] = df_wr["farm_power"].fillna(0.0)
 
-print("Farm solutions:")
-print(df_wr)
+# Use an interpolant to map the results back to the wind rose dataframe
+# Technically this could be done directly but an interpolant is safer
+# in the event the wind rose is irregular and/or ordered differently
+# than floris
+df_wr_op["farm_power"] = interpolant(df_wr_op[["wd", "ws"]])
+
+# Recombine with the below cut-in data
+df_wr_below_cut_in["farm_power"] = 0.
+df_wr = df_wr_below_cut_in.append(df_wr_op).sort_values(["wd","ws"])
 
 # Finally, calculate AEP in GWh
 aep = np.dot(df_wr["freq_val"], df_wr["farm_power"]) * 365 * 24
 
-print("Farm AEP: {:.2f} GWh".format(aep / 1.0e9))
+print("Farm AEP: {:.3f} GWh".format(aep / 1.0e9))

examples/08_opt_yaw_single_ws.py

@@ -36,6 +36,7 @@
     wind_directions=np.arange(0.0, 360.0, 3.0), 
     wind_speeds=[8.0],
 )
+print(fi.floris.farm.rotor_diameters)
 
 # Initialize optimizer object and run optimization using the Serial-Refine method
 yaw_opt = YawOptimizationSR(fi)#, exploit_layout_symmetry=False)

examples/09_optimize_yaw.py → examples/10_optimize_yaw.py

@@ -42,8 +42,8 @@ def load_floris():
     # Specify wind farm layout and update in the floris object
     N = 5  # number of turbines per row and per column
     X, Y = np.meshgrid(
-        5.0 * fi.floris.grid.reference_turbine_diameter * np.arange(0, N, 1),
-        5.0 * fi.floris.grid.reference_turbine_diameter * np.arange(0, N, 1),
+        5.0 * fi.floris.farm.rotor_diameters[0][0][0] * np.arange(0, N, 1),
+        5.0 * fi.floris.farm.rotor_diameters[0][0][0] * np.arange(0, N, 1),
     )
     fi.reinitialize(layout=(X.flatten(), Y.flatten()))
 
@@ -134,7 +134,7 @@ def calculate_aep(fi, df_windrose, column_name="farm_power"):
         exploit_layout_symmetry=True,
     )
 
-    df_opt = yaw_opt._optimize()
+    df_opt = yaw_opt.optimize()
     end_time = timerpc()
     t_tot = end_time - start_time
     t_fi = yaw_opt.time_spent_in_floris