update power_xylk, aep_xylk and aep_xy in flow map to support multidimensional...

update power_xylk, aep_xylk and aep_xy in flow map to support multidimensional power curves. fixes #66

py_wake/flow_map.py

@@ -71,26 +71,22 @@ class FlowMap(FlowBox):
     def XY(self):
         return self.X, self.Y
 
-    def power_xylk(self, wt_type=0, with_wake_loss=True):
+    def power_xylk(self, with_wake_loss=True, **wt_kwargs):
         if with_wake_loss:
             ws = self.WS_eff_xylk
 
         else:
             ws = self.WS_xylk
 
-        type = {'type': wt_type} if wt_type != 0 else {}
-
-        power_xylk = self.windFarmModel.windTurbines.power(ws, **type)
+        power_xylk = self.windFarmModel.windTurbines.power(ws, **wt_kwargs)
         return xr.DataArray(power_xylk[:, :, na], self.coords, dims=['y', 'x', 'h', 'wd', 'ws'])
 
-    def aep_xylk(self, wt_type=0, normalize_probabilities=False, with_wake_loss=True):
+    def aep_xylk(self, normalize_probabilities=False, with_wake_loss=True, **wt_kwargs):
         """Anual Energy Production of a potential wind turbine at all grid positions (x,y)
         for all wind directions (l) and wind speeds (k) in GWh.
 
         Parameters
         ----------
-        wt_type : Optional int, defaults to 0
-            Type of potential wind turbine
         normalize_propabilities : Optional bool, defaults to False
             In case only a subset of all wind speeds and/or wind directions is simulated,
             this parameter determines whether the returned AEP represents the energy produced in the fraction
@@ -103,20 +99,22 @@ class FlowMap(FlowBox):
         with_wake_loss : Optional bool, defaults to True
             If True, wake loss is included, i.e. power is calculated using local effective wind speed\n
             If False, wake loss is neglected, i.e. power is calculated using local free flow wind speed
+        wt_type : Optional arguments
+            Additional required/optional arguments needed by the WindTurbines to computer power, e.g. type, Air_density
          """
-        power_xylk = self.power_xylk(wt_type, with_wake_loss)
+        power_xylk = self.power_xylk(with_wake_loss, **wt_kwargs)
         P_xylk = self.P_xylk  # .isel.ilk((1,) + power_xylk.shape[2:])
         if normalize_probabilities:
             P_xylk = P_xylk / P_xylk.sum(['wd', 'ws'])
         return power_xylk * P_xylk * 24 * 365 * 1e-9
 
-    def aep_xy(self, wt_type=0, normalize_probabilities=False, with_wake_loss=True):
+    def aep_xy(self, normalize_probabilities=False, with_wake_loss=True, **wt_kwargs):
         """Anual Energy Production of a potential wind turbine at all grid positions (x,y)
         (sum of all wind directions and wind speeds)  in GWh.
 
         see aep_xylk
         """
-        return self.aep_xylk(wt_type, normalize_probabilities, with_wake_loss).sum(['wd', 'ws'])
+        return self.aep_xylk(normalize_probabilities, with_wake_loss, **wt_kwargs).sum(['wd', 'ws'])
 
     def plot(self, data, clabel, levels=100, cmap=None, plot_colorbar=True, plot_windturbines=True,
              normalize_with=1, ax=None):

py_wake/tests/test_flow_map.py

@@ -8,6 +8,8 @@ from py_wake.examples.data.iea37 import IEA37Site, IEA37_WindTurbines
 from py_wake import IEA37SimpleBastankhahGaussian
 import pytest
 from py_wake.deflection_models.jimenez import JimenezWakeDeflection
+from py_wake.wind_turbines._wind_turbines import WindTurbines
+from py_wake.examples.data import wtg_path
 
 
 @pytest.fixture(autouse=True)
@@ -31,6 +33,35 @@ def test_power_xylk():
     npt.assert_array_almost_equal(fm.power_xylk(with_wake_loss=False)[:, :, 0, 0] * 1e-6, 3.35)
 
 
+def test_power_xylk_wt_args():
+    site = IEA37Site(16)
+    x, y = site.initial_position.T
+    windTurbines = WindTurbines.from_WAsP_wtg(wtg_path + "Vestas V112-3.0 MW.wtg", default_mode=None)
+
+    # NOJ wake model
+    wind_farm_model = IEA37SimpleBastankhahGaussian(site, windTurbines)
+    simulation_result = wind_farm_model(x, y, wd=[0, 270], ws=[6, 8, 10], mode=0)
+    fm = simulation_result.flow_map(XYGrid(resolution=3))
+    npt.assert_array_almost_equal(fm.power_xylk(mode=1).sum(['wd', 'ws']).isel(h=0),
+                                  [[7030000., 6378864., 7029974.],
+                                   [7030000., 6144918., 4902029.],
+                                   [7030000., 7030000., 7029974.]], 0)
+    npt.assert_array_almost_equal(fm.power_xylk(mode=8).sum(['wd', 'ws']).isel(h=0),
+                                  [[8330000., 7577910., 8329970.],
+                                   [8330000., 7304188., 5837139.],
+                                   [8330000., 8330000., 8329970.]], 0)
+    # print(np.round(fm.power_xylk(mode=8).sum(['wd', 'ws']).squeeze()))
+
+    npt.assert_array_almost_equal(fm.aep_xylk(mode=1).sum(['x', 'y']).isel(h=0),
+                                  [[10., 24., 47.],
+                                   [75., 191., 375.]], 0)
+
+    npt.assert_array_almost_equal(fm.aep_xy(mode=1).isel(h=0),
+                                  [[88., 86., 88.],
+                                   [88., 68., 40.],
+                                   [88., 88., 88.]], 0)
+
+
 def test_YZGrid_perpendicular():
 
     site = IEA37Site(16)