power and ct funcs are now a list of functions (one element pr type) similar...

power and ct funcs are now a list of functions (one element pr type) similar to name, diameter and hub_height

py_wake/aep_calculator.py

@@ -103,4 +103,4 @@ class AEPCalculator():
         # return X_j, Y_j, aep_jlk.sum((1, 2)).reshape(X_j.shape)
 
         # same as above but requires less memory
-        return X_j, Y_j, ((self.windTurbines.power_func(WS_eff_jlk, type_j) * P_lk[na, :, :]).sum((1, 2)) * 24 * 365 * 1e-9).reshape(X_j.shape)
+        return X_j, Y_j, ((self.windTurbines.power(WS_eff_jlk, type_j) * P_lk[na, :, :]).sum((1, 2)) * 24 * 365 * 1e-9).reshape(X_j.shape)

py_wake/examples/data/hornsrev1.py

@@ -106,7 +106,8 @@ def main():
         print('Hub height', wt.hub_height())
         ws = np.arange(3, 25)
         import matplotlib.pyplot as plt
-        plt.plot(ws, wt.power_func(ws), '.-')
+        plt.plot(ws, wt.power(ws), '.-')
+        plt.show()
 
 
 main()

py_wake/tests/test_wake_models/test_noj.py

@@ -12,11 +12,11 @@ from py_wake.wind_turbines._wind_turbines import WindTurbines
 from py_wake.aep_calculator import AEPCalculator
 
 
-# Two turbines, 1: Nibe-A, 2:Ct=0
+# Two turbines, 0: Nibe-A, 1:Ct=0
 NibeA0 = WindTurbines(names=['Nibe-A'] * 2, diameters=[40] * 2,
                       hub_heights=[50] * 2,
-                      ct_func=lambda _, types: 8 / 9 * (types == 0),
-                      power_func=lambda *_: 0, power_unit='w')
+                      ct_funcs=[lambda _: 8 / 9, lambda _: 0],
+                      power_funcs=[lambda _: 0] * 2, power_unit='w')
 
 
 def test_NOJ_Nibe_result():

py_wake/wind_turbines/_wind_turbines.py

@@ -3,7 +3,7 @@ import numpy as np
 
 class WindTurbines():
 
-    def __init__(self, names, diameters, hub_heights, ct_func, power_func, power_unit):
+    def __init__(self, names, diameters, hub_heights, ct_funcs, power_funcs, power_unit):
         """Set of wind turbines
 
         Parameters
@@ -24,13 +24,12 @@ class WindTurbines():
         self._names = names
         self._diameters = np.array(diameters)
         self._hub_heights = np.array(hub_heights)
-        self.ct_func = ct_func
-
-        power_scale = {'w': 1, 'kw': 1e3, 'mw': 1e6, 'gw': 1e9}[power_unit.lower()]
-        if power_scale != 1:
-            self.power_func = lambda *args, **kwargs: power_func(*args, **kwargs) * power_scale
+        self.ct_funcs = ct_funcs
+        self.power_scale = {'w': 1, 'kw': 1e3, 'mw': 1e6, 'gw': 1e9}[power_unit.lower()]
+        if self.power_scale != 1:
+            self.power_funcs = [lambda ws: f(ws) * self.power_scale for f in power_funcs]
         else:
-            self.power_func = power_func
+            self.power_funcs = power_funcs
 
     def info(self, var, types):
         return var[np.asarray(types, np.int)]
@@ -44,8 +43,23 @@ class WindTurbines():
     def name(self, types=0):
         return self.info(self._names, types)
 
-    def ct_power(self, ws_i__, type_i=0):
-        return self.ct_func(ws_i__, type_i, ), self.power_func(ws_i__, type_i)
+    def power(self, ws_i, type_i):
+        return self.ct_power(ws_i, type_i)[1]
+
+    def ct(self, ws_i, type_i):
+        return self.ct_power(ws_i, type_i)[0]
+
+    def ct_power(self, ws_i, type_i=0):
+        if np.any(type_i != 0):
+            CT = np.zeros_like(ws_i)
+            P = np.zeros_like(ws_i)
+            for t in np.unique(type_i):
+                m = type_i == t
+                CT[m] = self.ct_funcs[t](ws_i[m])
+                P[m] = self.power_funcs[t](ws_i[m])
+            return CT, P
+        else:
+            return self.ct_funcs[0](ws_i), self.power_funcs[0](ws_i)
 
     def plot(self, x, y, types=None, ax=None):
         import matplotlib.pyplot as plt
@@ -72,26 +86,38 @@ class OneTypeWindTurbines(WindTurbines):
 
     def __init__(self, name, diameter, hub_height, ct_func, power_func, power_unit):
         WindTurbines.__init__(self, [name], [diameter], [hub_height],
-                              lambda ws, *_: ct_func(ws),
-                              lambda ws, *_: power_func(ws),
+                              [lambda ws: ct_func(ws)],
+                              [lambda ws: power_func(ws)],
                               power_unit)
 
+    def power(self, ws_i, type_i=0):
+        return self.ct_power(ws_i, type_i)[1]
+
+    def ct(self, ws_i, type_i=0):
+        return self.ct_power(ws_i, type_i)[0]
+
+
+def cube_power(ws_cut_in, ws_cut_out, ws_rated, power_rated):
+    def power_func(ws):
+        ws = np.asarray(ws)
+        power = np.zeros_like(ws, dtype=np.float)
+        m = (ws > ws_cut_in) & (ws < ws_rated)
+        power[m] = power_rated * ((ws[m] - ws_cut_in) / (ws_rated - ws_cut_in))**3
+        power[(ws >= ws_rated) & (ws <= ws_cut_out)] = power_rated
+        return power
+    return power_func
+
 
 def main():
     if __name__ == '__main__':
-        def power(ws, types=0):
-            """Calculate power in kW"""
-            rated = 2000 + (1000 * types)
-            return np.minimum(np.maximum((ws - 4)**3, 0), rated)
-
-        def ct(ws, types):
-            return 8 / 9
 
         wts = WindTurbines(names=['tb1', 'tb2'],
                            diameters=[80, 120],
                            hub_heights=[70, 110],
-                           ct_func=ct,
-                           power_func=power,
+                           ct_funcs=[lambda ws: 8 / 9,
+                                     lambda ws: 8 / 9],
+                           power_funcs=[cube_power(ws_cut_in=3, ws_cut_out=25, ws_rated=12, power_rated=2000),
+                                        cube_power(ws_cut_in=3, ws_cut_out=25, ws_rated=12, power_rated=3000)],
                            power_unit='kW')
 
         ws = np.arange(25)