Add MIT yaw correction model (3rd pass)

examples/examples_turbine/004_compare_yaw_loss.py

@@ -0,0 +1,87 @@
+"""
+Example: Change operation model and compare power loss in yaw.
+
+This example illustrates how to define different operational models and compares
+the power loss resulting from yaw misalignment across these various models.
+"""
+
+import itertools
+
+import matplotlib.pyplot as plt
+import numpy as np
+
+from floris import FlorisModel, TimeSeries
+
+
+# Parameters
+N = 101  # How many steps to cover yaw range in
+yaw_max = 30  # Maximum yaw angle to test
+
+# Set up the yaw angle sweep
+yaw_angles = np.zeros((N, 1))
+yaw_angles[:, 0] = np.linspace(-yaw_max, yaw_max, N)
+print(yaw_angles.shape)
+
+
+def evaluate_yawed_power(wsp: float, op_model: str) -> float:
+    print(f"Evaluating model: {op_model}   wind speed: {wsp} m/s")
+
+    # Grab model of FLORIS
+    fmodel = FlorisModel("../inputs/gch.yaml")
+
+    # Run N cases by setting up a TimeSeries (which is just several independent simulations)
+    wind_directions = np.ones(N) * 270.0
+    fmodel.set(
+        wind_data=TimeSeries(
+            wind_speeds=wsp,
+            wind_directions=wind_directions,
+            turbulence_intensities=0.06,
+        )
+    )
+
+    yaw_angles = np.array(
+        [(yaw, 0.0, 0.0) for yaw in np.linspace(-yaw_max, yaw_max, N)]
+    )
+    fmodel.set_operation_model(op_model)
+    fmodel.set(yaw_angles=yaw_angles)
+    fmodel.run()
+
+    # Save the power output results in kW
+    return fmodel.get_turbine_powers()[:, 0] / 1000
+
+
+# Loop over the operational models and wind speeds to compare
+op_models = ["simple", "cosine-loss", "mit-loss"]
+wind_speeds = [11.0, 11.5, 15.0]
+results = {}
+for op_model, wsp in itertools.product(op_models, wind_speeds):
+
+    # Save the power output results in kW
+    results[(op_model, wsp)] = evaluate_yawed_power(wsp, op_model)
+# Plot the results
+fig, axes = plt.subplots(1, len(wind_speeds), sharey=True)
+
+colors = ["C0", "k", "r"]
+linestyles = ["solid", "dashed", "dotted"]
+for wsp, ax in zip(wind_speeds, axes):
+    ax.set_title(f"wsp: {wsp} m/s")
+    ax.set_xlabel("Yaw angle [deg]")
+    ax.grid(True)
+    for op_model, c, ls in zip(op_models, colors, linestyles):
+
+        upstream_yaw_angle = yaw_angles[:, 0]
+        central_power = results[(op_model, wsp)][upstream_yaw_angle == 0]
+        ax.plot(
+            upstream_yaw_angle,
+            results[(op_model, wsp)] / central_power,
+            label=op_model,
+            color=c,
+            linestyle=ls,
+        )
+
+ax.grid(True)
+ax.legend()
+axes[0].set_xlabel("Yaw angle [deg]")
+axes[0].set_ylabel("Normalized turbine power [-]")
+
+plt.show()

floris/core/rotor_velocity.py

@@ -15,6 +15,7 @@
     NDArrayObject,
 )
 from floris.utilities import cosd
+from floris.core.turbine.mit_models.momentum import Heck
 
 
 def rotor_velocity_yaw_cosine_correction(
@@ -239,3 +240,49 @@ def rotor_velocity_air_density_correction(
     # Produce equivalent velocities at the reference air density
 
     return (air_density/ref_air_density)**(1/3) * velocities
+
+
+def mit_rotor_axial_induction(Cts: NDArrayFloat, yaw_angles: NDArrayFloat)-> NDArrayFloat:
+    """
+    Computes the axial induction of a yawed rotor given the yaw-aligned thrust
+    coefficient and yaw angles using the yawed actuator disk model developed at
+    MIT as described in Heck et al. 2023. Assumes the modified thrust
+    coefficient, C_T', is invariant to yaw misalignment angle.
+
+    Args
+        Cts (NDArrayFloat): Yaw-aligned thrust coefficient(s).
+        yaw_angles (NDArrayFloat): Rotor yaw angle(s) in degrees.
+
+    Returns: NDArrayFloat: Axial induction factor(s) of the yawed rotor.
+    """
+    Ctprime = - 4 * (Cts + 2 * np.sqrt(1 - Cts) - 2) / Cts
+    sol = Heck()(Ctprime, np.deg2rad(yaw_angles))
+
+    return sol.an
+
+
+def mit_rotor_velocity_yaw_correction(
+        Cts: NDArrayFloat,
+        yaw_angles: NDArrayFloat,
+        axial_inductions: NDArrayFloat,
+        rotor_effective_velocities: NDArrayFloat,
+) -> NDArrayFloat:
+    """
+    Computes adjusted rotor wind speeds given the yaw-aligned thrust
+    coefficient, yaw angles, and axial induction values using the yawed actuator
+    disk model developed at MIT as described in Heck et al. 2023. Assumes the
+    modified thrust coefficient, C_T', is invariant to yaw misalignment angle.
+
+    Args
+        Cts (NDArrayFloat): Yaw-aligned thrust coefficient(s).
+        yaw_angles (NDArrayFloat): Rotor yaw angle(s) in degrees.
+        axial_induction (NDArrayFloat): Rotor axial induction(s).
+        rotor_effective_velocities (NDArrayFloat) rotor effective wind speed(s) at the rotor. .
+
+    Returns: NDArrayFloat: corrected rotor effective wind speed(s) of the yawed rotor.
+    """
+    Ctprime = - 4 * (Cts + 2 * np.sqrt(1 - Cts) - 2) / Cts
+
+    ratio = (1 + 0.25 * Ctprime) * (1 - axial_inductions) * np.cos(np.deg2rad(yaw_angles))
+
+    return ratio * rotor_effective_velocities

floris/core/turbine/__init__.py

@@ -6,4 +6,5 @@
     PeakShavingTurbine,
     SimpleDeratingTurbine,
     SimpleTurbine,
+    MITTurbine,
 )

floris/core/turbine/mit_models/FixedPointIteration.py

@@ -0,0 +1,165 @@
+from dataclasses import dataclass
+from typing import Any, Callable, List, Protocol, Tuple
+
+import numpy as np
+from numpy.typing import ArrayLike
+
+
+class FixedPointIterationCompatible(Protocol):
+    def residual(self, *args, **kwargs) -> Tuple[ArrayLike]: ...
+
+    def initial_guess(self, *args, **kwargs) -> Tuple[ArrayLike]: ...
+
+
+@dataclass
+class FixedPointIterationResult:
+    converged: bool
+    niter: int
+    relax: float
+    max_resid: float
+    x: ArrayLike
+
+
+def _fixedpointiteration(
+    f: Callable[[ArrayLike, Any], np.ndarray],
+    x0: np.ndarray,
+    args=(),
+    kwargs={},
+    eps=0.00001,
+    maxiter=100,
+    relax=0,
+    callback=None,
+) -> FixedPointIterationResult:
+    """
+    Performs fixed-point iteration on function f until residuals converge or max
+    iterations is reached.
+
+    Args:
+        f (Callable): residual function of form f(x, *args, **kwargs) -> np.ndarray
+        x0 (np.ndarray): Initial guess
+        args (tuple): arguments to pass to residual function. Defaults to ().
+        kwargs (dict): keyword arguments to pass to residual function. Defaults to {}.
+        eps (float): Convergence tolerance. Defaults to 0.000001.
+        maxiter (int): Maximum number of iterations. Defaults to 100.
+        callback (Callable): optional callback function at each iteration of the form f(x0) -> None
+
+    Returns:
+        FixedPointIterationResult: Solution to residual function.
+    """
+
+    for c in range(maxiter):
+        residuals = f(x0, *args, **kwargs)
+
+        x0 = [_x0 + (1 - relax) * _r for _x0, _r in zip(x0, residuals)]
+        max_resid = [np.nanmax(np.abs(_r)) for _r in residuals]
+
+        if callback:
+            callback(x0)
+
+        if all(_r < eps for _r in max_resid):
+            converged = True
+            break
+    else:
+        converged = False
+
+    if maxiter == 0:
+        return FixedPointIterationResult(False, 0, np.nan, np.nan, x0)
+    return FixedPointIterationResult(converged, c, relax, max_resid, x0)
+
+
+def fixedpointiteration(
+    max_iter: int = 100, tolerance: float = 1e-6, relaxation: float = 0.0
+) -> FixedPointIterationCompatible:
+    """
+    Class decorator which adds a __call__ method to the class which performs
+    fixed-point iteration.
+
+    Args:
+        max_iter (int): Maximum number of iterations (default: 100)
+        tolerance (float): Convergence criteria (default: 1e-6)
+        relaxation (float): Relaxation factor between 0 and 1 (default: 0.0)
+
+    The class must contain 2 mandatory methods and 3
+    optional method:
+
+    mandatory:
+    initial_guess(self, *args, **kwargs)
+    residual(self, x, *args, **kwargs)
+
+    optional:
+    pre_process(self, *args, **kwargs) # Optional
+    post_process(self, result:FixedPointIterationResult) # Optional
+    callback(self, x) # Optional
+
+    """
+
+    def decorator(cls: FixedPointIterationCompatible) -> Callable:
+        def call(self, *args, **kwargs):
+            if hasattr(self, "pre_process"):
+                self.pre_process(*args, **kwargs)
+
+            callback = self.callback if hasattr(self, "callback") else None
+
+            x0 = self.initial_guess(*args, **kwargs)
+            result = _fixedpointiteration(
+                self.residual,
+                x0,
+                args=args,
+                kwargs=kwargs,
+                eps=tolerance,
+                maxiter=max_iter,
+                relax=relaxation,
+                callback=callback,
+            )
+
+            if hasattr(self, "post_process"):
+                return self.post_process(result, *args, **kwargs)
+            else:
+                return result
+
+        setattr(cls, "__call__", call)
+        return cls
+
+    return decorator
+
+
+def adaptivefixedpointiteration(
+    max_iter: int = 100, tolerance: float = 1e-6, relaxations: List[float] = [0.0]
+) -> Callable:
+    """
+    Class decorator which adds a __call__ method to the class which performs
+    fixed-point iteration. Same as `fixedpointiteration`, but takes a list of
+    relaxation factors, and iterates over all of them in order until convergence
+    is reached.
+    """
+
+    def decorator(cls: FixedPointIterationCompatible) -> Callable:
+        def call(self, *args, **kwargs):
+            if hasattr(self, "pre_process"):
+                self.pre_process(*args, **kwargs)
+            callback = self.callback if hasattr(self, "callback") else None
+
+            for relaxation in relaxations:
+                x0 = self.initial_guess(*args, **kwargs)
+                result = _fixedpointiteration(
+                    self.residual,
+                    x0,
+                    args=args,
+                    kwargs=kwargs,
+                    eps=tolerance,
+                    maxiter=max_iter,
+                    relax=relaxation,
+                    callback=callback,
+                )
+                if result.converged:
+                    break
+
+            if hasattr(self, "post_process"):
+                return self.post_process(result, *args, **kwargs)
+            else:
+                return result
+
+        setattr(cls, "__call__", call)
+        return cls
+
+    return decorator