Sg/kl model dev

examples/example_05_basic_windfarm.py

@@ -1,4 +1,5 @@
 from pathlib import Path
+import numpy as np
 
 import matplotlib.pyplot as plt
 from MITRotor import IEA15MW
@@ -12,10 +13,10 @@
 if __name__ == "__main__":
     windfarm = Windfarm(rotor_model=BEM(IEA15MW()))
 
-    setpoints = [
-        (0, 7, 0),
-        (0, 7, 0),
-        (0, 7, 0),
+    setpoints = [ # pitch, tsr, yaw, tilt
+        (0, 7, 0, 0),
+        (0, 7, 0, 0),
+        (0, 7, 0, 0),
     ]
 
     layout = Layout([0, 12 / 2, 24 / 2], [0, 0, 0])

examples/example_07_BEM_yaw_optimisation.py

@@ -6,43 +6,39 @@
 from mitwindfarm import Plotting
 from mitwindfarm.windfarm import Windfarm, WindfarmSolution, Layout
 from mitwindfarm.Rotor import BEM
+from MITRotor.Momentum import UnifiedMomentum
 from MITRotor.ReferenceTurbines import IEA15MW
 from rich import print
 import numpy as np
 
 FIGDIR = Path(__file__).parent.parent / "fig"
 FIGDIR.mkdir(exist_ok=True, parents=True)
 
-windfarm = Windfarm(rotor_model=BEM(IEA15MW()), TIamb=0.1)
-layout = Layout([0, 12, 24], [0.0, 0.5, 1.0])
-
+momentum_model= UnifiedMomentum(averaging = "rotor")
+windfarm = Windfarm(
+    rotor_model=BEM(IEA15MW(), momentum_model = momentum_model),
+    TIamb=0.1,
+)
+layout = Layout([0, 12, 24], [0.0, 0.5, 1.0], [0.0, 0.0, 0.0])
 
 def solve_for_setpoints(x) -> WindfarmSolution:
-    setpoints = [(x[0], x[1], x[2]), (x[3], x[4], x[5]), (x[6], x[7], x[8])]
+    setpoints = [(x[0], x[1], x[2], 0.0), (x[3], x[4], x[5], 0.0), (x[6], x[7], x[8], 0.0)]
     return windfarm(layout, setpoints)
 
-
 def objective_func(x):
     windfarm_sol = solve_for_setpoints(x)
     return -windfarm_sol.Cp
 
-
 if __name__ == "__main__":
-    x0 = [0, 7, 0, 0, 7, 0, 0, 7, 0]
+    nturbs = 3
+    angle_bounds = (-np.deg2rad(15), np.deg2rad(15))
+    tsr_bounds = (3, 10)
+    x0 = [0, 7, 0] * nturbs
+
     sol = minimize(
         objective_func,
         x0,
-        bounds=[
-            (-np.deg2rad(15), np.deg2rad(15)),
-            (3, 10),
-            (-np.deg2rad(15), np.deg2rad(15)),
-            (-np.deg2rad(15), np.deg2rad(15)),
-            (3, 10),
-            (-np.deg2rad(15), np.deg2rad(15)),
-            (-np.deg2rad(15), np.deg2rad(15)),
-            (3, 10),
-            (-np.deg2rad(15), np.deg2rad(15)),
-        ],
+        bounds=[angle_bounds,tsr_bounds,angle_bounds] * nturbs,
     )
     print(sol)
 
@@ -53,12 +49,17 @@ def objective_func(x):
     fig, axes = plt.subplots(2, 1)
     Plotting.plot_windfarm(windfarm_sol_ref, axes[0])
     Plotting.plot_windfarm(windfarm_sol_opt, axes[1])
+    fig.suptitle("Original and Optimized Wind Farm Layout for Yaw ($z = 0$)")
+    axes[0].set_xlabel("$x/D$")
+    axes[1].set_xlabel("$x/D$")
+    axes[0].set_ylabel("$y/D$")
+    axes[1].set_ylabel("$y/D$")
 
     axes[0].set_title(f"$C_P$: {windfarm_sol_ref.Cp:2.3f}")
     axes[1].set_title(
         f"$C_P$: {windfarm_sol_opt.Cp:2.3f} ({100*(windfarm_sol_opt.Cp/windfarm_sol_ref.Cp - 1):+2.1f}%)"
     )
-
+    plt.tight_layout()
     plt.savefig(
         FIGDIR / "example_07_BEM_yaw_optimisation.png", dpi=300, bbox_inches="tight"
     )

examples/example_08_curled_windfarm.py

@@ -20,7 +20,7 @@ def plot_example(powerlaw=False):
         zhub = 0
         base_windfield = Uniform(TIamb=0.05)  # 5% ambient TI
 
-    wf = CurledWindfarm(
+    wf_curled = CurledWindfarm(
         rotor_model=UnifiedAD_TI(),
         base_windfield=base_windfield,
         solver_kwargs=dict(
@@ -33,45 +33,67 @@ def plot_example(powerlaw=False):
     )
     wf_gauss = Windfarm(TIamb=0.05)  # 5% ambient TI, default model is Gaussian wake
     layout = Layout([0, 5, 10], [0, 0.4, 0.8], [zhub, zhub, zhub])  # non-dim by diameter D
-    setpoints = [  # for UnifiedAD_TI() rotor, set points are (Ctprime, yaw) tuple pairs
-        (2, np.radians(30)),
-        (2, np.radians(15)),
-        (2, 0),
-    ]  #  Example setpoints for two turbines
+    setpoints = [  # for UnifiedAD_TI() rotor, set points are (Ctprime, yaw, tilt) tuple pairs
+        (2, np.radians(30), np.radians(20)),
+        (2, np.radians(15), np.radians(10)),
+        (2, 0, 0),
+    ]  #  Example setpoints for three turbines
 
     # compute windfarm solutions (Cp)
     wf_solutions = []
-    for name, _wf in zip(["Curl", "Gauss"], [wf, wf_gauss]):
+    for name, _wf in zip(["Curl", "Gauss"], [wf_curled, wf_gauss]):
         time_st = time.time()
         sol = _wf(layout, setpoints)
         print(f"Windfarm solution for {name} in {time.time() - time_st:.2f} seconds")
         wf_solutions.append((name, sol))
 
-    # plot the comparison: wind speed and power
-    fig, axarr = plt.subplots(
-        figsize=(4 * len(wf_solutions), 4),
-        nrows=2,
+    # plot the comparison: wind speed (streamwise slice) and power
+    fig1, axarr1 = plt.subplots(
+        figsize=(4 * len(wf_solutions), 6),
+        nrows=3,
         ncols=len(wf_solutions),
         sharex=True,
         sharey="row",
-        height_ratios=(1, 2),
+        height_ratios=(1, 1, 2),
     )
-    for axs, (name, sol) in zip(axarr.T, wf_solutions):
-        plot_windfarm(sol, ax=axs[0], pad=2, axis=True)
+    for axs, (name, sol) in zip(axarr1.T, wf_solutions):
+        plot_windfarm(sol, ax=axs[0], z = zhub, pad=2, axis=True)
+        y_plot_val = 0
+        plot_windfarm(sol, ax=axs[1], y = y_plot_val, pad=2, axis=True)
         axs[0].set_xlabel("$x/D$")
-        axs[0].set_title(name)
+        axs[1].set_xlabel("$x/D$")
+        axs[0].set_title(name + f" at z = {zhub}")
+        axs[1].set_title(name + f" at y = {y_plot_val}")
         # plot power per turbine
-        axs[1].bar(layout.x, [r.Cp for r in sol.rotors], width=2)
-        if np.all(axs == (axarr.T)[0]):  # only label the first columns
+        axs[2].bar(layout.x, [r.Cp for r in sol.rotors], width=2)
+        if np.all(axs == (axarr1.T)[0]):  # only label the first columns
             axs[0].set_ylabel("$y/D$")
-            axs[1].set_ylabel("$C_P$")
-        axs[1].set_xticks(layout.x)
-        axs[1].set_xticklabels(np.arange(len(layout.x)) + 1)
-        axs[1].set_xlabel("Turbine row")
+            axs[1].set_ylabel("$z/D$")
+            axs[2].set_ylabel("$C_P$")
+        axs[2].set_xticks(layout.x)
+        axs[2].set_xticklabels(np.arange(len(layout.x)) + 1)
+        axs[2].set_xlabel("Turbine row")
+
+    plt.savefig(FIGDIR / f"{Path(__file__).stem}1.png", bbox_inches="tight")
 
-    plt.savefig(FIGDIR / f"{Path(__file__).stem}.png", bbox_inches="tight")
+    fig2, axarr2 = plt.subplots(
+        figsize=(4 * len(wf_solutions), 4),
+        nrows=1,
+        ncols=len(wf_solutions),
+        sharex=True,
+        sharey="row",
+    )
+    for axs, (name, sol) in zip(axarr2.T, wf_solutions):
+        # at first turbine's hub
+        x_plot_val = 4
+        plot_windfarm(sol, ax=axs, x = x_plot_val, pad=2, axis=True)
+        axs.set_title(name + f" at x = {x_plot_val}")
+        axs.set_xlabel("$y/D$")
+        axs.set_ylabel("$z/D$")
 
+
+    plt.savefig(FIGDIR / f"{Path(__file__).stem}2.png", bbox_inches="tight")
 
 if __name__ == "__main__":
-    plot_example(powerlaw=True)
+    plot_example(powerlaw=False)
     plt.close()

examples/example_09_BEM_tilt_optimization.py

@@ -0,0 +1,68 @@
+from pathlib import Path
+
+import matplotlib.pyplot as plt
+from scipy.optimize import minimize
+
+from mitwindfarm import Plotting
+from mitwindfarm.windfarm import Windfarm, WindfarmSolution, Layout
+from mitwindfarm.Rotor import BEM
+from MITRotor.ReferenceTurbines import IEA15MW
+from MITRotor.Momentum import UnifiedMomentum
+from rich import print
+import numpy as np
+
+FIGDIR = Path(__file__).parent.parent / "fig"
+FIGDIR.mkdir(exist_ok=True, parents=True)
+
+# Note that this setup is meant to mimic the example 7 yaw optimization setup.
+# Thus we are optimizing pitch, tsr, and tilt.
+# It becomes visually apparent that optimizing for yaw vs tilt return the same thing!
+momentum_model= UnifiedMomentum(averaging = "rotor")
+windfarm = Windfarm(
+    rotor_model=BEM(IEA15MW(), momentum_model = momentum_model),
+    TIamb=0.1,
+)
+layout = Layout([0, 12, 24], [0.0, 0.0, 0.0], [2.0, 1.5, 1.0])
+
+def solve_for_setpoints(x) -> WindfarmSolution:
+    setpoints = [(x[0], x[1], 0.0, x[2]), (x[3], x[4], 0.0, x[5]), (x[6], x[7], 0.0, x[8])]
+    return windfarm(layout, setpoints)
+
+def objective_func(x):
+    windfarm_sol = solve_for_setpoints(x)
+    return -windfarm_sol.Cp
+
+if __name__ == "__main__":
+    nturbs = 3
+    # note that we limit to positive tilt so the blades don't hit to shaft
+    half_angle_bounds = (0, np.deg2rad(15))
+    tsr_bounds = (3, 10)
+    x0 = [0, 7, 0] * nturbs
+
+    sol = minimize(
+        objective_func,
+        x0,
+        bounds=[half_angle_bounds,tsr_bounds,half_angle_bounds] * nturbs,
+    )
+    print(sol)
+
+    windfarm_sol_ref = solve_for_setpoints(x0)
+    windfarm_sol_opt = solve_for_setpoints(sol.x)
+    print(windfarm_sol_opt)
+
+    fig, axes = plt.subplots(2, 1)
+    Plotting.plot_windfarm(windfarm_sol_ref, axes[0], y = 0)
+    Plotting.plot_windfarm(windfarm_sol_opt, axes[1], y = 0)
+    fig.suptitle("Original and Optimized Wind Farm Layout for Tilt ($y = 0$)")
+    axes[0].set_xlabel("$x/D$")
+    axes[1].set_xlabel("$x/D$")
+    axes[0].set_ylabel("$y/D$")
+    axes[1].set_ylabel("$y/D$")
+    axes[0].set_title(f"$C_P$: {windfarm_sol_ref.Cp:2.3f}")
+    axes[1].set_title(
+        f"$C_P$: {windfarm_sol_opt.Cp:2.3f} ({100*(windfarm_sol_opt.Cp/windfarm_sol_ref.Cp - 1):+2.1f}%)"
+    )
+
+    plt.savefig(
+        FIGDIR / "example_09_BEM_tilt_optimisation.png", dpi=300, bbox_inches="tight"
+    )

examples/example_10_BEM_yaw_tilt.py

@@ -0,0 +1,73 @@
+from pathlib import Path
+
+import matplotlib.pyplot as plt
+from scipy.optimize import minimize
+
+from mitwindfarm import Plotting
+from mitwindfarm.windfarm import Windfarm, WindfarmSolution, Layout
+from mitwindfarm.Rotor import BEM
+from MITRotor.Momentum import UnifiedMomentum
+from MITRotor.ReferenceTurbines import IEA15MW
+from UnifiedMomentumModel.Utilities.Geometry import calc_eff_yaw
+from rich import print
+import numpy as np
+
+FIGDIR = Path(__file__).parent.parent / "fig"
+FIGDIR.mkdir(exist_ok=True, parents=True)
+
+momentum_model= UnifiedMomentum(averaging = "rotor")
+windfarm = Windfarm(
+    rotor_model=BEM(IEA15MW(), momentum_model = momentum_model),
+    TIamb=0.1,
+)
+layout = Layout([0, 12], [0.0, 0.0], [0.0, 0.0])
+
+def solve_for_setpoints(x) -> WindfarmSolution:
+    setpoints = [(x[0], x[1], x[2], 0.0), (x[3], x[4], x[5], 0.0), (x[6], x[7], x[8], 0.0)]
+    return windfarm(layout, setpoints)
+
+def objective_func(x):
+    windfarm_sol = solve_for_setpoints(x)
+    return -windfarm_sol.Cp
+
+if __name__ == "__main__":
+    pitch_val = 0
+    tsr_val = 7
+
+    yaw_tilt_val = np.deg2rad(15)
+    eff_yaw = calc_eff_yaw(yaw_tilt_val, yaw_tilt_val)
+    yaw_val = eff_yaw
+    tilt_val = eff_yaw
+
+    aligned_setpoints = (pitch_val, tsr_val, 0, 0)
+    setpoints_yaw = [(pitch_val, tsr_val, eff_yaw, 0), aligned_setpoints]
+    setpoints_tilt = [(pitch_val, tsr_val, 0, eff_yaw), aligned_setpoints]
+    setpoints_yaw_and_tilt = [(pitch_val, tsr_val, yaw_tilt_val, yaw_tilt_val), aligned_setpoints]   
+
+    yaw_sol = windfarm(layout, setpoints_yaw)
+    tilt_sol = windfarm(layout, setpoints_tilt) 
+    yaw_and_tilt_sol = windfarm(layout, setpoints_yaw_and_tilt) 
+
+    fig, axes = plt.subplots(2, 3)
+    Plotting.plot_windfarm(yaw_sol, axes[0, 0], z = 0)
+    Plotting.plot_windfarm(tilt_sol, axes[0, 1], z = 0)
+    Plotting.plot_windfarm(yaw_and_tilt_sol, axes[0, 2], z = 0)
+    axes[0, 0].set_ylabel("$z/D$")
+
+    Plotting.plot_windfarm(yaw_sol, axes[1, 0], y = 0)
+    Plotting.plot_windfarm(tilt_sol, axes[1, 1], y = 0)
+    Plotting.plot_windfarm(yaw_and_tilt_sol, axes[1, 2], y = 0)
+    axes[1, 0].set_ylabel("$y/D$")
+    axes[1, 0].set_xlabel("$x/D$")
+    axes[1, 1].set_xlabel("$x/D$")
+    axes[1, 2].set_xlabel("$x/D$")
+
+    fig.suptitle("$C_P$ for Equivalent Turbine Setups", size = 16)
+    deg_yaw_tilt, deg_eff_yaw = np.rad2deg(yaw_tilt_val), np.rad2deg(eff_yaw)
+    axes[0, 0].set_title(f"Yaw {np.round(deg_eff_yaw, decimals=1)}$^\circ$\n$C_P$: {yaw_sol.Cp:2.3f}")
+    axes[0, 1].set_title(f"Tilt {np.round(deg_eff_yaw, decimals=1)}$^\circ$\n$C_P$: {tilt_sol.Cp:2.3f}")
+    axes[0, 2].set_title(f"Yaw {np.round(deg_yaw_tilt, decimals=1)}$^\circ$ & Tilt {np.round(deg_yaw_tilt, decimals=1)}$^\circ$\n$C_P$: {yaw_and_tilt_sol.Cp:2.3f}")
+    plt.tight_layout()
+    plt.savefig(
+        FIGDIR / "example_10_yaw_tilt_comparison.png", dpi=300, bbox_inches="tight"
+    )