apply rotor average model to turbulence too

py_wake/rotor_avg_models/rotor_avg_model.py

-from py_wake.deficit_models.deficit_model import DeficitModel
 import numpy as np
-from abc import abstractmethod
 from numpy import newaxis as na
 
 
-class RotorAvgModel(DeficitModel):
+class RotorAvgModel():
     """Wrap a DeficitModel.
     The RotorAvgModel
     - add an extra dimension (one or more points covering the downstream rotors)
@@ -13,22 +11,19 @@ class RotorAvgModel(DeficitModel):
     """
     args4rotor_avg_deficit = ['hcw_ijlk', 'dh_ijl', 'D_dst_ijl']
 
-    def calc_deficit(self, deficitModel, D_dst_ijl, **kwargs):
-        self.deficitModel = deficitModel
-        if np.all(D_dst_ijl == 0):
-            return self.deficitModel.calc_deficit(D_dst_ijl=D_dst_ijl, **kwargs)
-        else:
-            return self._calc_rotor_avg_deficit(D_dst_ijl=D_dst_ijl, **kwargs)
-
     def calc_deficit_convection(self, deficitModel, D_dst_ijl, **kwargs):
         self.deficitModel = deficitModel
         return self.deficitModel.calc_deficit_convection(D_dst_ijl=D_dst_ijl, **kwargs)
 
-    @abstractmethod
-    def _calc_rotor_avg_deficit(self):
-        """Similar to calc_deficit, but with an extra point dimension to calculate the
-        rotor average wind speed as a weighted average of a number of points instead of
-        only the rotor center"""
+    def __call__(self, func, D_dst_ijl, **kwargs):
+        # add extra dimension, p, with 40 points distributed over the destination rotors
+        kwargs = self._update_kwargs(D_dst_ijl=D_dst_ijl, **kwargs)
+
+        values_ijlkp = func(**kwargs)
+        # Calculate weighted sum of deficit over the destination rotors
+        if self.nodes_weight is None:
+            return np.mean(values_ijlkp, -1)
+        return np.sum(self.nodes_weight[na, na, na, na, :] * values_ijlkp, -1)
 
 
 class RotorCenter(RotorAvgModel):
@@ -37,8 +32,8 @@ class RotorCenter(RotorAvgModel):
     nodes_y = [0]
     nodes_weight = [1]
 
-    def _calc_rotor_avg_deficit(self, **kwargs):
-        return self.deficitModel.calc_deficit(**kwargs)
+    def __call__(self, func, **kwargs):
+        return func(**kwargs)
 
     def _calc_layout_terms(self, deficitModel, **kwargs):
         deficitModel._calc_layout_terms(**kwargs)
@@ -59,26 +54,13 @@ class GridRotorAvg(RotorAvgModel):
         hcw_ijlkp = hcw_ijlk[..., na] + R_dst_ijl[:, :, :, na, na] * self.nodes_x[na, na, na, na, :]
         dh_ijlp = dh_ijl[..., na] + R_dst_ijl[:, :, :, na] * self.nodes_y[na, na, na, :]
         new_kwargs = {'dh_ijl': dh_ijlp, 'hcw_ijlk': hcw_ijlkp, 'D_dst_ijl': D_dst_ijl[..., na]}
-        if 'cw_ijlk' in self.deficitModel.args4deficit:
-            cw_ijlkp = np.sqrt(hcw_ijlkp**2 + dh_ijlp[:, :, :, na]**2)
-            new_kwargs['cw_ijlk'] = cw_ijlkp
-        if 'D_dst_ijl' in self.deficitModel.args4deficit:
-            new_kwargs['D_dst_ijl'] = D_dst_ijl
+
+        new_kwargs['cw_ijlk'] = np.sqrt(hcw_ijlkp**2 + dh_ijlp[:, :, :, na]**2)
+        new_kwargs['D_dst_ijl'] = D_dst_ijl
 
         new_kwargs.update({k: v[..., na] for k, v in kwargs.items() if k not in new_kwargs})
         return new_kwargs
 
-    def _calc_rotor_avg_deficit(self, **kwargs):
-
-        # add extra dimension, p, with 40 points distributed over the destination rotors
-        kwargs = self._update_kwargs(**kwargs)
-
-        deficit_ijlkp = self.deficitModel.calc_deficit(**kwargs)
-        # Calculate weighted sum of deficit over the destination rotors
-        if self.nodes_weight is None:
-            return np.mean(deficit_ijlkp, -1)
-        return np.sum(self.nodes_weight[na, na, na, na, :] * deficit_ijlkp, -1)
-
     def _calc_layout_terms(self, deficitModel, **kwargs):
         self.deficitModel = deficitModel
         self.deficitModel._calc_layout_terms(**self._update_kwargs(**kwargs))

py_wake/tests/test_ground_models/test_mirror.py

@@ -6,14 +6,11 @@ import matplotlib.pyplot as plt
 from py_wake.flow_map import YZGrid
 import numpy as np
 from py_wake.tests import npt
-from py_wake.wind_turbines import WindTurbines, OneTypeWindTurbines
+from py_wake.wind_turbines import WindTurbines
 from py_wake.superposition_models import LinearSum, SquaredSum
 from py_wake.wind_farm_models.engineering_models import PropagateDownwind, All2AllIterative
 import pytest
-from py_wake.deficit_models.gaussian import BastankhahGaussianDeficit, IEA37SimpleBastankhahGaussianDeficit,\
-    ZongGaussianDeficit, NiayifarGaussianDeficit
-from py_wake.deficit_models.fuga import FugaDeficit
-from py_wake.deficit_models.gcl import GCLDeficit
+from py_wake.deficit_models.gaussian import ZongGaussianDeficit
 from py_wake.turbulence_models.stf import STF2017TurbulenceModel
 from py_wake.ground_models.ground_models import MirrorSquaredSum
 
@@ -42,7 +39,6 @@ def test_Mirror_NOJ():
         plt.legend()
         plt.show()
     plt.close()
-    print(res)
     npt.assert_array_equal(res, ref)
 
 

py_wake/tests/test_deficit_models/test_rotor_avg_models.py → py_wake/tests/test_rotor_avg_models.py

@@ -13,7 +13,7 @@ import pytest
 from py_wake.turbulence_models.stf import STF2017TurbulenceModel
 
 
-def test_RotorGridAvg():
+def test_RotorGridAvg_deficit():
     site = IEA37Site(16)
     x, y = site.initial_position.T
     windTurbines = IEA37_WindTurbines()
@@ -53,6 +53,49 @@ def test_RotorGridAvg():
     plt.close()
 
 
+def test_RotorGridAvg_ti():
+    site = IEA37Site(16)
+    x, y = site.initial_position.T
+    windTurbines = IEA37_WindTurbines()
+    wfm = IEA37SimpleBastankhahGaussian(site,
+                                        windTurbines,
+                                        turbulenceModel=STF2017TurbulenceModel())
+    flow_map = wfm([0, 500], [0, 0], wd=270, ws=10).flow_map(HorizontalGrid(x=[500], y=np.arange(-100, 100)))
+    plt.plot(flow_map.Y[:, 0], flow_map.TI_eff_xylk[:, 0, 0, 0])
+    R = windTurbines.diameter() / 2
+
+    for name, rotorAvgModel, ref1 in [
+            ('RotorCenter', RotorCenter(), 0.22292190804089568),
+            ('RotorGrid2', EqGridRotorAvg(2), 0.21135016790319944),
+            ('RotorGrid3', EqGridRotorAvg(3), 0.20618819397725846),
+            ('RotorGrid4', EqGridRotorAvg(4), 0.20324660406762962),
+            ('RotorGrid100', EqGridRotorAvg(100), 0.1989725533174574),
+            ('RotorGQGrid_4,3', GQGridRotorAvg(4, 3), 0.19874837617113356)]:
+
+        # test with PropagateDownwind
+        wfm = IEA37SimpleBastankhahGaussian(site,
+                                            windTurbines,
+                                            rotorAvgModel=rotorAvgModel,
+                                            turbulenceModel=STF2017TurbulenceModel())
+        sim_res = wfm([0, 500], [0, 0], wd=270, ws=10)
+        npt.assert_almost_equal(sim_res.TI_eff_ilk[1, 0, 0], ref1)
+
+        # test with All2AllIterative
+        wfm = All2AllIterative(site, windTurbines,
+                               IEA37SimpleBastankhahGaussianDeficit(),
+                               rotorAvgModel=rotorAvgModel,
+                               superpositionModel=SquaredSum(),
+                               turbulenceModel=STF2017TurbulenceModel())
+        sim_res = wfm([0, 500], [0, 0], wd=270, ws=10)
+        npt.assert_almost_equal(sim_res.TI_eff_ilk[1, 0, 0], ref1)
+
+        plt.plot([-R, R], [sim_res.TI_eff_ilk[1, 0, 0]] * 2, label=name)
+    if 0:
+        plt.legend()
+        plt.show()
+    plt.close()
+
+
 def test_gauss_quadrature():
 
     x, _, w = gauss_quadrature(4, 1)

py_wake/wind_farm_models/engineering_models.py

@@ -86,6 +86,7 @@ class EngineeringWindFarmModel(WindFarmModel):
             self.args4deficit = set(self.args4deficit) | set(self.groundModel.args4deficit)
         self.args4all = set(self.args4deficit)
         if self.turbulenceModel:
+            self.args4addturb = set(self.turbulenceModel.args4addturb) | set(self.rotorAvgModel.args4rotor_avg_deficit)
             self.args4all = self.args4all | set(self.turbulenceModel.args4addturb)
         if self.deflectionModel:
             self.args4all = self.args4all | set(self.deflectionModel.args4deflection)
@@ -134,9 +135,8 @@ class EngineeringWindFarmModel(WindFarmModel):
 
     def _calc_deficit(self, dw_ijlk, **kwargs):
         """Calculate wake (and blockage) deficit"""
-        def f(**kwargs):
-            return self.rotorAvgModel.calc_deficit(self.wake_deficitModel, **kwargs)
-        deficit = self.groundModel(f, dw_ijlk=dw_ijlk, **kwargs)
+        deficit = self.groundModel(lambda **kwargs: self.rotorAvgModel(self.wake_deficitModel.calc_deficit, **kwargs),
+                                   dw_ijlk=dw_ijlk, **kwargs)
         return self._add_blockage(deficit, dw_ijlk, **kwargs)
 
     def _calc_deficit_convection(self, dw_ijlk, **kwargs):
@@ -486,13 +486,16 @@ class PropagateDownwind(EngineeringWindFarmModel):
 
                 if self.turbulenceModel:
 
-                    if 'wake_radius_ijlk' in self.turbulenceModel.args4addturb:
+                    if 'wake_radius_ijlk' in self.args4addturb:
                         wake_radius_ijlk = self.wake_deficitModel.wake_radius(dw_ijlk=dw_ijlk, **args)
                         arg_funcs['wake_radius_ijlk'] = lambda: wake_radius_ijlk
 
-                    turb_args = {k: arg_funcs[k]() for k in self.turbulenceModel.args4addturb if k != "dw_ijlk"}
+                    turb_args = {k: arg_funcs[k]() for k in self.args4addturb
+                                 if k != "dw_ijlk"}
+
                     # Calculate added turbulence
-                    add_turb_nk[n_dw] = self.turbulenceModel.calc_added_turbulence(dw_ijlk=dw_ijlk, **turb_args)
+                    add_turb_nk[n_dw] = self.rotorAvgModel(self.turbulenceModel.calc_added_turbulence,
+                                                           dw_ijlk=dw_ijlk, **turb_args)
 
         WS_eff_jlk, power_jlk, ct_jlk = np.array(WS_eff_mk), np.array(power_jlk), np.array(ct_jlk)
 
@@ -593,7 +596,7 @@ class All2AllIterative(EngineeringWindFarmModel):
             # Calculate effective wind speed
             WS_eff_ilk = self.superpositionModel.calc_effective_WS(lw.WS_ilk, deficit_iilk)
             if self.turbulenceModel:
-                add_turb_ijlk = self.turbulenceModel.calc_added_turbulence(**args)
+                add_turb_ijlk = self.rotorAvgModel(self.turbulenceModel.calc_added_turbulence, **args)
                 TI_eff_ilk = self.turbulenceModel.calc_effective_TI(lw.TI_ilk, add_turb_ijlk)
 
             # Check if converged