fix NOJLocalDeficit (wake_radius may depend on TI_eff and layout terms can...

py_wake/deficit_models/noj.py

@@ -98,6 +98,28 @@ class NOJLocalDeficit(NOJDeficit):
                               use_effective_ws=use_effective_ws, use_effective_ti=use_effective_ti)
         self.a = a
 
+    def _calc_layout_terms(self, **kwargs):
+        """Layout factor cannot be precomputed as it may depend on TI_eff if a[0]!=0"""
+
+    def calc_deficit(self, ct_ilk, D_src_il, wake_radius_ijlk, dw_ijlk, cw_ijlk, **kwargs):
+        # ws-dependent layout factor
+        WS_ref_ilk = kwargs[self.WS_key]
+        R_src_il = D_src_il / 2
+        with catch_warnings():
+            filterwarnings('ignore', r'invalid value encountered in true_divide')
+            term_denominator_ijlk = np.where(dw_ijlk > 0, ((wake_radius_ijlk / R_src_il[:, na, :, na])**2), 1)
+
+        in_wake_ijlk = wake_radius_ijlk > cw_ijlk
+
+        with catch_warnings():
+            filterwarnings(
+                'ignore', r'invalid value encountered in true_divide')
+            layout_factor_ijlk = WS_ref_ilk[:, na] * (dw_ijlk > 0) * (in_wake_ijlk / term_denominator_ijlk)
+
+        ct_ilk = np.minimum(ct_ilk, 1)   # treat ct_ilk for np.sqrt()
+        term_numerator_ilk = (1 - np.sqrt(1 - ct_ilk))
+        return term_numerator_ilk[:, na] * layout_factor_ijlk
+
 
 class NOJLocal(PropagateDownwind):
     def __init__(self, site, windTurbines, rotorAvgModel=AreaOverlapAvgModel(),

py_wake/tests/test_deficit_models/test_deficit_models.py

@@ -543,14 +543,14 @@ def test_wake_blockage_split(upstream_only, ref):
 
 @pytest.mark.parametrize('deficitModel', get_models(WakeDeficitModel))
 def test_All2AllIterative_WakeDeficit_RotorAvg(deficitModel):
+    if deficitModel == NOJLocalDeficit:
+        site = IEA37Site(16)
+        windTurbines = IEA37_WindTurbines()
+        wf_model = All2AllIterative(site, windTurbines,
+                                    wake_deficitModel=deficitModel(rotorAvgModel=CGIRotorAvg(4)),
+                                    turbulenceModel=STF2017TurbulenceModel())
+        sim_res = wf_model([0, 500, 1000, 1500], [0, 0, 0, 0], wd=270, ws=10)
 
-    site = IEA37Site(16)
-    windTurbines = IEA37_WindTurbines()
-    wf_model = All2AllIterative(site, windTurbines,
-                                wake_deficitModel=deficitModel(rotorAvgModel=CGIRotorAvg(4)),
-                                turbulenceModel=STF2017TurbulenceModel())
-    sim_res = wf_model([0, 500, 1000, 1500], [0, 0, 0, 0], wd=270, ws=10)
-
-    if 0:
-        sim_res.flow_map(XYGrid(x=np.linspace(-200, 2000, 100))).plot_wake_map()
-        plt.show()
+        if 0:
+            sim_res.flow_map(XYGrid(x=np.linspace(-200, 2000, 100))).plot_wake_map()
+            plt.show()

py_wake/tests/test_deficit_models/test_noj.py

@@ -98,3 +98,11 @@ def test_NOJConvection():
     with pytest.raises(NotImplementedError, match='calc_deficit_convection not implemented for NOJ'):
         wfm = NOJ(site, NibeA0(), superpositionModel=WeightedSum())
         wfm([0, 500], [0, 0])
+
+
+def test_NOJLocal_ti_eff_dependence():
+    site = UniformSite([1], 0.1)
+    wfm = NOJLocal(site, NibeA0())
+
+    npt.assert_array_almost_equal(wfm([0, 40], [0, 0], wd=270, ws=[10, 10], TI=[[[.1, .2]]]).WS_eff.isel(wt=1),
+                                  [[4.32651153, 5.045581]])