Sg/floris interface

.gitignore

@@ -15,6 +15,11 @@ fig*
 *.tar
 *.whl
 
+# Allow these CSVs needed for the floris interface
+!MITRotor/FlorisInterface/IEA_15mw_rotor.csv
+
+Validation/
+
 ### Python ###
 # Byte-compiled / optimized / DLL files
 __pycache__/
@@ -188,4 +193,6 @@ poetry.toml
 # LSP config files
 pyrightconfig.json
 
-# End of https://www.toptal.com/developers/gitignore/api/python
+# End of https://www.toptal.com/developers/gitignore/api/python
+
+.DS_Store

MITRotor/Aerodynamics.py

@@ -7,7 +7,7 @@
 from numpy.typing import ArrayLike
 
 from .RotorDefinition import RotorDefinition
-from .Geometry import BEMGeometry
+from .Geometry import BEMGeometry, expand_to_Nr_Ntheta
 from UnifiedMomentumModel.Utilities.Geometry import calc_eff_yaw
 
 __all__ = [
@@ -269,13 +269,16 @@ def __call__(
             AerodynamicProperties: Calculated aerodynamic properties stored in AerodynamicProperties object.
 
         """
+        tsr = expand_to_Nr_Ntheta(tsr)
+        pitch = expand_to_Nr_Ntheta(pitch)
         # calculate values in "yaw-only" frame
         local_yaw = -self.eff_yaw
         Vax = U * ((1 - an) * np.cos(local_yaw))
+        theta_eff = geom.theta_mesh + self.delta_theta
         Vtan = (
             (1 + aprime) * tsr * geom.mu_mesh
             - U * (1 - an)
-            * np.cos(self.eff_theta_mesh)
+            * np.cos(theta_eff)
             * np.sin(local_yaw)
         )
 
@@ -291,8 +294,8 @@ def __call__(
             an = an,
             aprime = aprime,
             solidity = solidity,
-            U = U * np.ones(geom.shape),
-            wdir = wdir * np.ones(geom.shape),
+            U = U * np.ones_like(geom.mu_mesh),
+            wdir = wdir * np.ones_like(geom.mu_mesh),
             Vax = Vax,
             Vtan = Vtan,
             aoa = aoa,

MITRotor/BEMSolver.py

@@ -7,12 +7,13 @@
 
 from . import Momentum, TipLoss
 from .Aerodynamics import AerodynamicModel, AerodynamicProperties, DefaultAerodynamics
-from .Geometry import BEMGeometry
+from .Geometry import BEMGeometry, expand_to_Np, expand_to_Nr_Ntheta
 from .RotorDefinition import RotorDefinition
 from .TangentialInduction import DefaultTangentialInduction, TangentialInductionModel
 from UnifiedMomentumModel.Utilities.Geometry import calc_eff_yaw
 
 
+
 def average(geometry: BEMGeometry, value: ArrayLike, grid: Literal["sector", "annulus", "rotor"] = "rotor"):
     # Assuming the function returns a 2D grid of values
 
@@ -58,13 +59,13 @@ def solidity(self, grid: Literal["sector ", "annulus", "rotor"] = "rotor"):
     def U(self, grid: Literal["sector", "annulus", "rotor"] = "rotor"):
         return average(self.geom, self.aero_props.U, grid)
 
-    def wdir(self, grid: Literal["sector ", "annulus", "rotor"] = "rotor"):
+    def wdir(self, grid: Literal["sector", "annulus", "rotor"] = "rotor"):
         return average(self.geom, self.aero_props.wdir, grid)
 
     def Vax(self, grid: Literal["sector", "annulus", "rotor"] = "rotor"):
         return average(self.geom, self.aero_props.Vax, grid)
 
-    def Vtan(self, grid: Literal["sector ", "annulus", "rotor"] = "rotor"):
+    def Vtan(self, grid: Literal["sector", "annulus", "rotor"] = "rotor"):
         return average(self.geom, self.aero_props.Vtan, grid)
 
     def W(self, grid: Literal["sector", "annulus", "rotor"] = "rotor"):
@@ -91,29 +92,27 @@ def Ctan(self, grid: Literal["sector", "annulus", "rotor"] = "rotor"):
     def Cx(self, grid: Literal["sector", "annulus", "rotor"] = "rotor"):
         return average(self.geom, self.aero_props.C_x_corr, grid)
 
-    def Ctau(self, grid: Literal["sector ", "annulus", "rotor"] = "rotor"):
+    def Ctau(self, grid: Literal["sector", "annulus", "rotor"] = "rotor"):
         return average(self.geom, self.aero_props.C_tau_corr, grid)
 
-    def Ctau_uncorr(self, grid: Literal["sector ", "annulus", "rotor"] = "rotor"):
+    def Ctau_uncorr(self, grid: Literal["sector", "annulus", "rotor"] = "rotor"):
         return average(self.geom, self.aero_props.C_tau, grid)
 
     def F(self, grid: Literal["sector", "annulus", "rotor"] = "rotor"):
         return average(self.geom, self.aero_props.F, grid)
 
     def Cp(self, grid: Literal["sector", "annulus", "rotor"] = "rotor"):
-        dCp = (
-            self.tsr
-            * self.geom.mu_mesh
-            * self.Ctau_uncorr(grid="sector")
-        )
+        tsr = np.asarray(self.tsr)
+        if tsr.ndim == 1:
+            tsr = expand_to_Nr_Ntheta(tsr)
+        dCp = (tsr * self.geom.mu_mesh * self.Ctau_uncorr(grid="sector"))
         return average(self.geom, dCp, grid=grid)
 
     def Cp_corr(self, grid: Literal["sector", "annulus", "rotor"] = "rotor"):
-        dCp = (
-            self.tsr
-            * self.geom.mu_mesh
-            * self.Ctau(grid="sector")
-        )
+        tsr = np.asarray(self.tsr)
+        if tsr.ndim == 1:
+            tsr = expand_to_Nr_Ntheta(tsr)
+        dCp = (tsr * self.geom.mu_mesh * self.Ctau(grid="sector"))
         return average(self.geom, dCp, grid=grid)
 
     def Ct(self, grid: Literal["sector", "annulus", "rotor"] = "rotor"):
@@ -169,21 +168,37 @@ def sample_points(self, yaw: float = 0.0, tilt: float = 0.0) -> tuple[ArrayLike,
         return X, Y, Z
 
     def pre_process(self, pitch, tsr, yaw = 0, tilt = 0, **kwargs):
+        pitch, tsr = np.asarray(pitch), np.asarray(tsr)
+        yaw, tilt = np.asarray(yaw), np.asarray(tilt)
+        self.scalar_inputs = (pitch.ndim == 0) & (tsr.ndim == 0) & (yaw.ndim == 0) & (tilt.ndim == 0)
+        if not self.scalar_inputs:
+            assert len(pitch) == len(tsr) == len(yaw) == len(tilt), "Setpoint arrays should be the same lenght"
         # switch reference frame to a "yaw-only" frame where y' is aligned with the lateral wake
-        self.aerodynamic_model.eff_yaw = calc_eff_yaw(yaw, tilt)
-        if tilt == 0:
-            dtheta = 0
-        elif yaw == 0:
-            dtheta = np.pi / 2
-        else: # non-zero yaw and tilt
-            sin_eff = np.sin(self.aerodynamic_model.eff_yaw)
-            dtheta = np.arccos(np.sin(yaw) / sin_eff)
-        self.aerodynamic_model.eff_theta_mesh = self.geometry.theta_mesh + dtheta
+        yaw, tilt = np.broadcast_arrays(yaw, tilt)
+        eff_yaw = calc_eff_yaw(yaw, tilt)
+        # initialize dtheta
+        dtheta = np.zeros_like(eff_yaw, dtype=float)
+        # masks
+        yaw_zero = (yaw == 0)
+        not_tilt_zero = np.logical_not(tilt == 0)
+        # case1: yaw == 0 and tilt != 0
+        case1 = yaw_zero & not_tilt_zero
+        dtheta[case1] = np.pi / 2
+        # case2: yaw != 0 and tilt != 0
+        case2 = np.logical_not(yaw_zero) & not_tilt_zero
+        dtheta[case2] = np.arccos(
+            np.sin(yaw[case2]) / np.sin(eff_yaw[case2])
+        )
+        # expand everything to the correct dimensions to ensure broadcasting
+        self.aerodynamic_model.eff_yaw = expand_to_Nr_Ntheta(eff_yaw)
+        self.aerodynamic_model.delta_theta = expand_to_Nr_Ntheta(dtheta)
+        self.geometry.mu_mesh = expand_to_Np(self.geometry.mu_mesh)
+        self.geometry.theta_mesh = expand_to_Np(self.geometry.theta_mesh)
         return
 
     def initial_guess(self, *args, **kwargs) -> Tuple[ArrayLike, ...]:
-        a = (1 / 3) * np.ones(self.geometry.shape)
-        aprime = np.zeros(self.geometry.shape)
+        a = (1 / 3) * np.ones_like(self.geometry.mu_mesh)
+        aprime = np.zeros_like(self.geometry.mu_mesh)
 
         return a, aprime
 
@@ -198,8 +213,8 @@ def residual(
         tilt: ArrayLike = 0.0,
     ) -> Tuple[ArrayLike, ...]:
         an, aprime = x
-        U = np.ones(self.geometry.shape) if U is None else U
-        wdir = np.zeros(self.geometry.shape) if wdir is None else wdir
+        U = np.ones_like(self.geometry.mu_mesh) if U is None else U
+        wdir = np.zeros_like(self.geometry.mu_mesh) if wdir is None else wdir
 
         aero_props = self.aerodynamic_model(
             an = an, 
@@ -221,12 +236,25 @@ def residual(
         return e_an, e_aprime
 
     def post_process(self, result: FixedPointIterationResult, pitch, tsr, yaw = 0, U=None, wdir=None, tilt = 0.0) -> BEMSolution:
-        U = np.ones(self.geometry.shape) if U is None else U
-        wdir = np.zeros(self.geometry.shape) if wdir is None else wdir
+        U = np.ones_like(self.geometry.mu_mesh) if U is None else U
+        wdir = np.zeros_like(self.geometry.mu_mesh) if wdir is None else wdir
         an, aprime = result.x
         aero_props = self.aerodynamic_model(an, aprime, pitch, tsr, yaw, self.rotor, self.geometry, U, wdir, tilt = tilt)
         aero_props.F = self.tiploss_model(aero_props, pitch, tsr, yaw, self.rotor, self.geometry, tilt = tilt)
         avg_Ct = average(self.geometry, aero_props.C_x)
         u4,v4,w4 = self.momentum_model.compute_initial_wake_velocities(avg_Ct, yaw, tilt = tilt)
 
+        if self.scalar_inputs: # if all setpoints were scalars
+            # return single values as scalars
+            pitch, tsr, yaw, tilt = [np.asarray(x).item() for x in (pitch, tsr, yaw, tilt)]
+            u4, v4, w4 = [np.asarray(x).item() for x in (u4, v4, w4)]
+            # remove unneeded extra axis from Ntheta x Nr arrays
+            an = np.squeeze(an)
+            aprime = np.squeeze(aprime)
+            self.geometry.theta_mesh = np.squeeze(self.geometry.theta_mesh)
+            self.geometry.mu_mesh = np.squeeze(self.geometry.mu_mesh)
+            for key, value in vars(aero_props).items():
+                if isinstance(value, np.ndarray):
+                    setattr(aero_props, key, np.squeeze(value))
+
         return BEMSolution(pitch, tsr, yaw, aero_props, self.geometry, result.converged, result.niter, u4, v4, tilt = tilt, w4 = w4)

MITRotor/FlorisInterface/FlorisInterface.py

@@ -0,0 +1,136 @@
+import os
+import numpy as np
+import polars as pl
+from attrs import define, field
+from typing import Optional
+from scipy.interpolate import interp1d
+# FLORIS Imports
+from floris.type_dec import floris_float_type, NDArrayFloat
+from floris.core.turbine.operation_models import BaseOperationModel
+from floris.core.rotor_velocity import average_velocity
+# MITRotor / UMM Imports
+from MITRotor.ReferenceTurbines import IEA15MW
+from MITRotor.Momentum import UnifiedMomentum
+from MITRotor.Geometry import BEMGeometry
+from MITRotor.TipLoss import NoTipLoss
+from MITRotor.BEMSolver import BEM
+
+# default rotor if none provided by user (IEA 15MW)
+def default_bem_factory():
+    return BEM(
+        rotor=IEA15MW(),
+        momentum_model=UnifiedMomentum(averaging="rotor"),
+        geometry=BEMGeometry(Nr=10, Ntheta=20),
+        tiploss_model=NoTipLoss()
+    )
+# pitch vs windspeed interpolater if none provided by user
+# for IEA 15MW from figure 2 (https://docs.nrel.gov/docs/fy22osti/82134.pdf)
+def default_pitch_interp():
+    module_dir = os.path.dirname(__file__)
+    pitch_file = os.path.join(module_dir, "IEA_15mw_rotor.csv")
+    df = pl.read_csv(pitch_file)
+    wind_table = df["Wind [m/s]"].to_numpy()
+    pitch_table = df["Pitch [deg]"].to_numpy()
+    # TODO: should fill_value be extrapolate?
+    return interp1d(wind_table, pitch_table, kind="linear", fill_value="extrapolate", bounds_error=False)
+
+# tsr vs windspeed interpolater if none provided by user
+# for IEA 15MW from figure 2 (https://docs.nrel.gov/docs/fy22osti/82134.pdf)
+def default_tsr_interp():
+    module_dir = os.path.dirname(__file__)
+    tsr_file = os.path.join(module_dir, "IEA_15mw_rotor.csv")
+    df = pl.read_csv(tsr_file)
+    wind_table = df["Wind [m/s]"].to_numpy()
+    tip_speed_table = df["Tip Speed [m/s]"].to_numpy()
+    tsr_table = tip_speed_table / wind_table
+    # TODO: should fill_value be extrapolate?
+    return interp1d(wind_table, tsr_table, kind="linear", fill_value="extrapolate", bounds_error=False)
+
+@define
+class MITRotorTurbine(BaseOperationModel):
+    """
+    Turbine operation model as described by Liew et al. (2024).
+
+    Args:
+        bem_model (BEM): optional BEM model as defined in MITRotor, defaults to IEA15MW with UMM momentum model
+        pitch_csv (str): optional path to pitch trajectory based on wind speed, defaults to IEA15MW Figure 2 (https://docs.nrel.gov/docs/fy22osti/82134.pdf)
+        tsr_csv (str)): optional path to tsr trajectory based on wind speed, defaults to IEA15MW Figure 2 (https://docs.nrel.gov/docs/fy22osti/82134.pdf)
+
+    Methods:
+        power
+        thrust_coefficient
+        axial_induction
+    """
+    # user can define a BEM model if they want a different rotor, momentum model, or geometry
+    bem_model = field(init = True, factory = default_bem_factory, type = BEM)
+
+    # create interp objects based on pitch and tsr csvs
+    pitch_interp = field(init=True, factory=default_pitch_interp, type = interp1d, repr = False)
+    tsr_interp = field(init=True, factory=default_tsr_interp, type = interp1d, repr = False)
+
+    # save most recent solution by unique floris arguments
+    _last_key = field(init=False, default=None, type = bytes)
+    _a = field(init=False, default=None, type = NDArrayFloat)
+    _Ct = field(init=False, default=None, type = NDArrayFloat)
+    _power = field(init=False, default=None, type = NDArrayFloat)
+
+    def _get_state_key(self, velocities: np.ndarray, yaw_angles: np.ndarray, tilt_angles: np.ndarray) -> tuple:
+        # saves key to uniquely identify farm state -> avoids re-solving for calls to power, thrust, and induction for same state
+        return velocities.tobytes(), yaw_angles.tobytes(), tilt_angles.tobytes()
+
+    def _update_solution(self,
+        velocities: NDArrayFloat,
+        air_density: float,
+        yaw_angles: NDArrayFloat,
+        tilt_angles: NDArrayFloat,
+        average_method: str = "cubic-mean",
+        cubature_weights: Optional[NDArrayFloat] = None,
+        **_,
+    ):
+        # create cache key for current inputs
+        key = self._get_state_key(velocities, yaw_angles, tilt_angles)
+        # update solution if conditions are different
+        if key != self._last_key:
+            n_findex, n_turbines = yaw_angles.shape
+
+            # save new key and clear fields
+            self._last_key = key
+            self._a = np.empty((n_findex, n_turbines), dtype=floris_float_type)
+            self._Ct = np.empty((n_findex, n_turbines), dtype=floris_float_type)
+            self._power = np.empty((n_findex, n_turbines), dtype=floris_float_type)
+
+            # compute the power-effective wind speed across the rotor
+            rotor_average_velocities = average_velocity(
+                velocities=velocities,
+                method=average_method,
+                cubature_weights=cubature_weights,
+            )
+
+            # calculate rotor area
+            rotor_area = np.pi * self.bem_model.rotor.R**2 
+
+                        # get setpoints
+            yaw, tilt = np.deg2rad(yaw_angles), np.deg2rad(tilt_angles)
+            pitch = np.deg2rad(self.pitch_interp(rotor_average_velocities))
+            tsr = self.tsr_interp(rotor_average_velocities)
+            for tindex in range(n_turbines):
+                # solve BEM
+                bem_sol = self.bem_model(pitch[:, tindex], tsr[:, tindex], yaw = yaw[:, tindex], tilt = tilt[:, tindex])
+                # get induction and thrust coeff
+                self._a[:, tindex] = bem_sol.a()
+                self._Ct[:, tindex] = bem_sol.Ct()
+                # compute power
+                self._power[:, tindex] = 0.5 * bem_sol.Cp() * air_density * rotor_area * (rotor_average_velocities[:, tindex])**3
+        return
+
+    def power(self, **kwargs) -> NDArrayFloat:
+        self._update_solution(**kwargs)
+        return self._power
+
+    def thrust_coefficient(self, **kwargs) -> NDArrayFloat:
+        self._update_solution(**kwargs)
+        return self._Ct
+
+    def axial_induction(self, **kwargs) -> NDArrayFloat:
+        self._update_solution(**kwargs)
+        return self._a