Update uncertainty interface to 4d, new API (#821)

examples/20_calculate_farm_power_with_uncertainty.py

@@ -1,4 +1,3 @@
-
 import matplotlib.pyplot as plt
 import numpy as np
 
@@ -16,41 +15,106 @@
 """
 
 # Instantiate FLORIS using either the GCH or CC model
-fi = FlorisInterface("inputs/gch.yaml") # GCH model
-fi_unc = UncertaintyInterface("inputs/gch.yaml") # Add uncertainty with default settings
+fi = FlorisInterface("inputs/gch.yaml")  # GCH model
+fi_unc_3 = UncertaintyInterface(
+    "inputs/gch.yaml", verbose=True, wd_std=3
+)
+fi_unc_5 = UncertaintyInterface(
+    "inputs/gch.yaml", verbose=True, wd_std=5
+)
 
 # Define a two turbine farm
 D = 126.0
-layout_x = np.array([0, D*6, D*12])
-layout_y = [0, 0, 0]
-wd_array = np.arange(0.0, 360.0, 1.0)
-fi.reinitialize(layout_x=layout_x, layout_y=layout_y, wind_directions=wd_array)
-fi_unc.reinitialize(layout_x=layout_x, layout_y=layout_y, wind_directions=wd_array)
-
-# Define a matrix of yaw angles to be all 0
-# Note that yaw angles is now specified as a matrix whose dimesions are
-# wd/ws/turbine
-num_wd = len(wd_array) # Number of wind directions
-num_ws = 1 # Number of wind speeds
-num_turbine = len(layout_x) #  Number of turbines
-yaw_angles = np.zeros((num_wd, num_ws, num_turbine))
-
-# Calculate the nominal wake solution
-fi.calculate_wake(yaw_angles=yaw_angles)
-
-# Calculate the nominal wind farm power production
-farm_powers_nom = fi.get_farm_power() / 1e3
+layout_x = np.array([0, D * 6])
+layout_y = [0, 0]
+wd_array = np.arange(240.0, 300.0, 1.0)
+wind_speeds = 8.0 * np.ones_like(wd_array)
+fi.set(layout_x=layout_x, layout_y=layout_y, wind_directions=wd_array, wind_speeds=wind_speeds)
+fi_unc_3.set(
+    layout_x=layout_x, layout_y=layout_y, wind_directions=wd_array, wind_speeds=wind_speeds
+)
+fi_unc_5.set(
+    layout_x=layout_x, layout_y=layout_y, wind_directions=wd_array, wind_speeds=wind_speeds
+)
+
+
+# Run both models
+fi.run()
+fi_unc_3.run()
+fi_unc_5.run()
 
-# Calculate the wind farm power with uncertainty on the wind direction
-fi_unc.calculate_wake(yaw_angles=yaw_angles)
-farm_powers_unc = fi_unc.get_farm_power() / 1e3
+# Collect the nominal and uncertain farm power
+turbine_powers_nom = fi.get_turbine_powers() / 1e3
+turbine_powers_unc_3 = fi_unc_3.get_turbine_powers() / 1e3
+turbine_powers_unc_5 = fi_unc_5.get_turbine_powers() / 1e3
+farm_powers_nom = fi.get_farm_power() / 1e3
+farm_powers_unc_3 = fi_unc_3.get_farm_power() / 1e3
+farm_powers_unc_5 = fi_unc_5.get_farm_power() / 1e3
 
 # Plot results
-fig, ax = plt.subplots()
-ax.plot(wd_array, farm_powers_nom.flatten(), color='k',label='Nominal farm power')
-ax.plot(wd_array, farm_powers_unc.flatten(), color='r',label='Farm power with uncertainty')
+fig, axarr = plt.subplots(1, 3, figsize=(15, 5))
+ax = axarr[0]
+ax.plot(wd_array, turbine_powers_nom[:, 0].flatten(), color="k", label="Nominal power")
+ax.plot(
+    wd_array,
+    turbine_powers_unc_3[:, 0].flatten(),
+    color="r",
+    label="Power with uncertainty = 3 deg",
+)
+ax.plot(
+    wd_array, turbine_powers_unc_5[:, 0].flatten(), color="m", label="Power with uncertainty = 5deg"
+)
+ax.grid(True)
+ax.legend()
+ax.set_xlabel("Wind Direction (deg)")
+ax.set_ylabel("Power (kW)")
+ax.set_title("Upstream Turbine")
+
+ax = axarr[1]
+ax.plot(wd_array, turbine_powers_nom[:, 1].flatten(), color="k", label="Nominal power")
+ax.plot(
+    wd_array,
+    turbine_powers_unc_3[:, 1].flatten(),
+    color="r",
+    label="Power with uncertainty = 3 deg",
+)
+ax.plot(
+    wd_array,
+    turbine_powers_unc_5[:, 1].flatten(),
+    color="m",
+    label="Power with uncertainty = 5 deg",
+)
+ax.set_title("Downstream Turbine")
 ax.grid(True)
 ax.legend()
-ax.set_xlabel('Wind Direction (deg)')
-ax.set_ylabel('Power (kW)')
+ax.set_xlabel("Wind Direction (deg)")
+ax.set_ylabel("Power (kW)")
+
+ax = axarr[2]
+ax.plot(wd_array, farm_powers_nom.flatten(), color="k", label="Nominal farm power")
+ax.plot(
+    wd_array, farm_powers_unc_3.flatten(), color="r", label="Farm power with uncertainty = 3 deg"
+)
+ax.plot(
+    wd_array, farm_powers_unc_5.flatten(), color="m", label="Farm power with uncertainty = 5 deg"
+)
+ax.set_title("Farm Power")
+ax.grid(True)
+ax.legend()
+ax.set_xlabel("Wind Direction (deg)")
+ax.set_ylabel("Power (kW)")
+
+# Compare the AEP calculation
+freq = np.ones_like(wd_array)
+freq = freq / freq.sum()
+
+aep_nom = fi.get_farm_AEP(freq=freq)
+aep_unc_3 = fi_unc_3.get_farm_AEP(freq=freq)
+aep_unc_5 = fi_unc_5.get_farm_AEP(freq=freq)
+
+print(f"AEP without uncertainty {aep_nom}")
+print(f"AEP without uncertainty (3 deg) {aep_unc_3} ({100*aep_unc_3/aep_nom:.2f}%)")
+print(f"AEP without uncertainty (5 deg) {aep_unc_5} ({100*aep_unc_5/aep_nom:.2f}%)")
+
+
 plt.show()