py-wake - npm Package Compare versions

Comparing version 2.6.16 to 2.6.17

METADATA

 Metadata-Version: 2.4
 Name: py_wake
-Version: 2.6.16
+Version: 2.6.17
 Summary: Open source static wake modeling framework from DTU
@@ -5,0 +5,0 @@ Author: DTU Wind Energy

py_wake/deficit_models/deficit_model.py

@@ -173,3 +173,3 @@ from abc import ABC, abstractmethod
         deficit_centre_ijlk : array_like
-            Wind speed deficit caused by the i'th turbine at j'th downstream location, without accounting for crosswind distance (ie cw = 0)
+            Wind speed centerline deficit caused by the i'th turbine at j'th downstream location, without accounting for crosswind distance (ie cw = 0)
         uc_ijlk : array_like
@@ -187,3 +187,3 @@ Convection velocity of the i'th turbine at locations j
                  method='linear', bounds='limit',
-                 rotorAvgModel=None, groundModel=None, use_effective_ws=True, use_effective_ti=False):
+                 superpositionModel=None, rotorAvgModel=None, groundModel=None, use_effective_ws=True, use_effective_ti=False):
         """
@@ -215,4 +215,6 @@ Parameters
         XRLUTModel.__init__(self, da, get_input, get_output, method=method, bounds=bounds)
-        BlockageDeficitModel.__init__(self, upstream_only=True, rotorAvgModel=rotorAvgModel, groundModel=groundModel)
-        WakeDeficitModel.__init__(self, rotorAvgModel=rotorAvgModel, groundModel=groundModel,
+        BlockageDeficitModel.__init__(self, upstream_only=True, superpositionModel=superpositionModel,
+                                      rotorAvgModel=rotorAvgModel, groundModel=groundModel)
+        WakeDeficitModel.__init__(self, superpositionModel=superpositionModel, rotorAvgModel=rotorAvgModel,
+                                  groundModel=groundModel,
                                   use_effective_ws=use_effective_ws, use_effective_ti=use_effective_ti)
@@ -219,0 +221,0 @@

py_wake/deficit_models/rans_lut.py

 import xarray as xr
-from py_wake.deficit_models.deficit_model import XRLUTDeficitModel
+from py_wake.deficit_models.deficit_model import XRLUTDeficitModel, ConvectionDeficitModel
+from py_wake.utils.model_utils import XRLUTModel
 from py_wake import np
-from py_wake.superposition_models import LinearSum
+from py_wake.superposition_models import LinearSum, WeightedSum
 from py_wake.wind_farm_models.engineering_models import All2AllIterative
@@ -10,6 +11,7 @@ from py_wake.turbulence_models.rans_lut_turb import RANSLUTTurbulence
 from pathlib import Path
+from numpy import newaxis as na
 
 
 class RANSLUT(All2AllIterative):
-    def __init__(self, lut, site, windTurbines, rotorAvgModel=None):
+    def __init__(self, lut, site, windTurbines, convlut=None, rotorAvgModel=None):
         """
@@ -39,10 +41,24 @@ Parameters
         lut_lst = [load_lut(lut) for lut in lut_lst]
-        deficit = RANSLUTDeficit(lut=[lut.deficits for lut in lut_lst], rotorAvgModel=rotorAvgModel)
+
+        if convlut is not None:
+            # Use convection wake deficit model for WeigthedSum superposition
+            superpositionModel = WeightedSum(weight_limit=1.0)
+            if not isinstance(convlut, (list, tuple)):
+                convlut = [convlut]
+            convlut_lst = convlut
+            convlut_lst = [load_lut(convlut) for convlut in convlut_lst]
+            wake_deficit = RANSLUTConvDeficit(lut=[lut.deficits for lut in lut_lst], convlut=convlut_lst, rotorAvgModel=rotorAvgModel, superpositionModel=WeightedSum(weight_limit=1.0))
+        else:
+            # Use linear wake superposition
+            superpositionModel = LinearSum()
+            wake_deficit = RANSLUTDeficit(lut=[lut.deficits for lut in lut_lst], rotorAvgModel=rotorAvgModel, superpositionModel=LinearSum())
+
+        blockage_deficit = RANSLUTDeficit(lut=[lut.deficits for lut in lut_lst], rotorAvgModel=rotorAvgModel, superpositionModel=LinearSum())
         turb = RANSLUTTurbulence(lut=[lut.added_ti for lut in lut_lst], rotorAvgModel=rotorAvgModel)
 
         All2AllIterative.__init__(self, site, windTurbines,
-                                  wake_deficitModel=deficit,
-                                  blockage_deficitModel=deficit,
+                                  wake_deficitModel=wake_deficit,
+                                  blockage_deficitModel=blockage_deficit,
                                   turbulenceModel=turb,
-                                  superpositionModel=LinearSum())
+                                  superpositionModel=superpositionModel)
 
@@ -53,7 +69,7 @@
 
-    def __init__(self, lut, rotorAvgModel=None, groundModel=None, use_effective_ws=True,
+    def __init__(self, lut, superpositionModel=None, rotorAvgModel=None, groundModel=None, use_effective_ws=True,
                  use_effective_ti=False):
         assert use_effective_ws, "RANSLUTDeficit only makes sense when scaling with effective wind speed"
         XRLUTDeficitModel.__init__(self, self.get_lut(lut, 'deficits'),
-                                   bounds='limit', rotorAvgModel=rotorAvgModel, groundModel=groundModel,
+                                   bounds='limit', superpositionModel=superpositionModel, rotorAvgModel=rotorAvgModel, groundModel=groundModel,
                                    use_effective_ws=True, use_effective_ti=use_effective_ti)
@@ -67,4 +83,57 @@
 
+class RANSLUTConvDeficit(RANSLUTModel, XRLUTDeficitModel, ConvectionDeficitModel):
+    """Expects LUT velocity deficit and convection variables xarray"""
+
+    def __init__(self, lut, convlut, superpositionModel=None, rotorAvgModel=None, groundModel=None, use_effective_ws=True,
+                 use_effective_ti=False):
+        assert use_effective_ws, "RANSLUTDeficit only makes sense when scaling with effective wind speed"
+        XRLUTDeficitModel.__init__(self, self.get_lut(lut, 'deficits'),
+                                   bounds='limit', superpositionModel=superpositionModel, rotorAvgModel=rotorAvgModel, groundModel=groundModel,
+                                   use_effective_ws=True, use_effective_ti=use_effective_ti)
+        self.sigma_lut = XRLUTModel(self.get_lut([lut.sigma for lut in convlut], 'sigma'),
+                                    self.get_conv_input, bounds='limit')
+        self.center_def_lut = XRLUTModel(self.get_lut([lut.centerline_deficit for lut in convlut], 'centerline_deficit'),
+                                         self.get_conv_input, bounds='limit')
+
+    def get_conv_input(self, dw_ijlk, D_src_il, TI_eff_ilk, ct_ilk, type_il, **kwargs):
+        kwargs = dict(x=dw_ijlk / D_src_il[:, na, :, na],
+                      ti=TI_eff_ilk[:, na],
+                      ct=ct_ilk[:, na],
+                      type_i=type_il[:, na, :, na])
+
+        return [kwargs[k] for k in self.sigma_lut.da_lst[0].dims]
+
+    def sigma_ijlk(self, D_src_il, dw_ijlk, ct_ilk, WS_ref_ijlk, **kwargs):
+        # dimensional wake expansion
+        sigma_ijlk = self.sigma_lut.__call__(D_src_il=D_src_il, dw_ijlk=dw_ijlk, ct_ilk=ct_ilk, **kwargs)
+        # Replace undefined sigma values (upstream of WT) with zeros
+        sigma_ijlk[:, :, :][dw_ijlk[:, :, :, 0] < 0] = 0.0
+        return sigma_ijlk * D_src_il[:, na, :, na]
+
+    def _calc_deficit(self, D_src_il, dw_ijlk, ct_ilk, WS_ref_ijlk, **kwargs):
+        # dimensional wake expansion
+        sigma_sqr_ijlk = (self.sigma_ijlk(D_src_il=D_src_il, dw_ijlk=dw_ijlk, ct_ilk=ct_ilk, WS_ref_ijlk=WS_ref_ijlk, **kwargs))**2
+        deficit_centre_ijlk = WS_ref_ijlk * np.minimum(1.0, self.center_def_lut.__call__(D_src_il=D_src_il, dw_ijlk=dw_ijlk, ct_ilk=ct_ilk, **kwargs))
+        # Replace undefined center deficit values (upstream of WT) with zeros
+        deficit_centre_ijlk[:, :, :][dw_ijlk[:, :, :, 0] < 0] = 0.0
+        return WS_ref_ijlk, sigma_sqr_ijlk, deficit_centre_ijlk
+
+    def calc_deficit_convection(self, D_src_il, dw_ijlk, cw_ijlk, ct_ilk, **kwargs):
+        if self.groundModel:  # pragma: no cover
+            raise NotImplementedError(
+                "calc_deficit_convection (WeightedSum) cannot be used in combination with GroundModels")
+        WS_ref_ijlk, sigma_sqr_ijlk, deficit_centre_ijlk = self._calc_deficit(
+            D_src_il, dw_ijlk, ct_ilk, **kwargs, **self.get_WS_ref_kwargs(kwargs))
+        # Convection velocity
+        uc_ijlk = WS_ref_ijlk - 0.5 * deficit_centre_ijlk
+        sigma_sqr_ijlk = np.broadcast_to(sigma_sqr_ijlk, deficit_centre_ijlk.shape)
+        return deficit_centre_ijlk, uc_ijlk, sigma_sqr_ijlk
+
+    @property
+    def args4deficit(self):
+        return XRLUTModel.args4model.fget(self) | ConvectionDeficitModel.args4deficit.fget(self)
+
+
 class RANSLUTDemoDeficit(RANSLUTDeficit):
-    def __init__(self, rotorAvgModel=None, groundModel=None, use_effective_ws=True, use_effective_ti=False):
+    def __init__(self, superpositionModel=None, rotorAvgModel=None, groundModel=None, use_effective_ws=True, use_effective_ti=False):
         # Load default RANS LUT file based on a demo V80 LUT
@@ -74,3 +143,3 @@ demo_lut = ptf('ranslut/V80_ranslut_demo.nc',
 
-        RANSLUTDeficit.__init__(self, demo_lut, rotorAvgModel=rotorAvgModel, groundModel=groundModel,
+        RANSLUTDeficit.__init__(self, demo_lut, superpositionModel=superpositionModel, rotorAvgModel=rotorAvgModel, groundModel=groundModel,
                                 use_effective_ws=use_effective_ws, use_effective_ti=use_effective_ti)

py_wake/site/shear.py

@@ -53,3 +53,4 @@ from py_wake import np
         L_inv = self.h_zeta / self.h_ref  # 1 / Obukhov length
-        return (np.log(h[:, na, na] / z0) - psi(h[:, na, na] * L_inv, Cm1=self.Cm1, Cm2=self.Cm2)) / \
+        h = np.reshape(h, (h.shape + (1, 1))[:3])
+        return (np.log(h / z0) - psi(h * L_inv, Cm1=self.Cm1, Cm2=self.Cm2)) / \
             (np.log(self.h_ref / z0) - psi(self.h_zeta, Cm1=self.Cm1, Cm2=self.Cm2)) * WS_ilk
@@ -56,0 +57,0 @@

py_wake/superposition_models.py

@@ -82,3 +82,3 @@ from py_wake import np
 
-    def __init__(self, delta=0.01, max_err=1e-3, max_iter=5):
+    def __init__(self, delta=0.01, max_err=1e-3, max_iter=5, weight_limit=10 ** 10):
         # minimum deficit (as fraction of free-stream) to invoke weighted summation
@@ -90,10 +90,12 @@ self.delta = delta
         self.max_iter = max_iter
+        # Limit weights to this value. A limit of 1.0 means not exceeding linear superposition
+        self.weight_limit = weight_limit
 
-    def __call__(self, centerline_deficit_jxxx, WS_xxx,
+    def __call__(self, deficit_jxxx, WS_xxx,
                  convection_velocity_jxxx,
-                 sigma_sqr_jxxx, cw_jxxx, hcw_jxxx, dh_jxxx):
+                 sigma_sqr_jxxx, cw_jxxx, hcw_jxxx, dh_jxxx, deficit_centre_jxxx):
 
-        Ws = WS_xxx + np.zeros(centerline_deficit_jxxx.shape[1:])
+        Ws = WS_xxx + np.zeros(deficit_centre_jxxx.shape[1:])
 
-        usc = centerline_deficit_jxxx
+        usc = deficit_centre_jxxx
         uc = convection_velocity_jxxx
@@ -107,3 +109,4 @@ sigma_sqr = sigma_sqr_jxxx
         # Local deficit
-        us = usc * np.exp(-1 / (2 * sigma_sqr) * cw**2)
+        us = usc * np.exp(-1 / (2 * (sigma_sqr + 1e-30)) * cw**2)
+        us[usc == 0] = 0.0
 
@@ -115,2 +118,3 @@ # Set lower deficit limit below which deficits are linearly added
 
+        ucn = np.ones_like(uc)
         # Only start weighting computation if at least two deficits need to be combined
@@ -153,3 +157,4 @@ # Computations are performed where velocities exceed the specified limit (us_lim).
                 Inz = Uc != 0
-                ucn[:, Inz] = uc[:, Inz] / Uc[Inz]
+                # Limit weights to self.weight_limit
+                ucn[:, Inz] = np.minimum(uc[:, Inz] / Uc[Inz], self.weight_limit)
 
@@ -200,3 +205,7 @@ # Combined local deficit
                 count += 1
-        return Us + np.sum(np.where(~Il, us, 0), axis=0)
+        # Replace nan with 1.0 for cases where convection variable are undefined (i.e. upstream of wt)
+        ucn = np.nan_to_num(ucn, nan=1.0)
+        # Reset weights to linear sum where WeightedSum is not applied
+        ucn[~Il] = 1.0
+        return np.sum(deficit_jxxx * ucn, axis=0)
 
@@ -203,0 +212,0 @@

py_wake/tests/test_deficit_models/test_rans_lut.py

@@ -14,3 +14,3 @@ import matplotlib.pyplot as plt
 from py_wake.utils.profiling import timeit
-from py_wake.utils.rans_lut_utils import ADControl, get_Ellipsys_equivalent_output
+from py_wake.utils.rans_lut_utils import ADControl, get_Ellipsys_equivalent_output, fit_gauss
 from py_wake.wind_farm_models.engineering_models import All2AllIterative
@@ -20,2 +20,3 @@ from py_wake.wind_turbines import WindTurbine
 from py_wake.wind_turbines.power_ct_functions import PowerCtTabular
+from py_wake.rotor_avg_models.rotor_avg_model import GQGridRotorAvg
 
@@ -108,3 +109,3 @@
 
-    npt.assert_array_almost_equal(aDControl.U_CT_CP_AD, aDControl_saved.U_CT_CP_AD, 10)
+    npt.assert_array_almost_equal(aDControl.U_CT_CP_AD_ws, aDControl_saved.U_CT_CP_AD_ws, 10)
 
@@ -225,3 +226,3 @@ wfm = RANSLUT(dataset, site, wts)
 
-    npt.assert_array_almost_equal(aDControl.U_CT_CP_AD, aDControl_saved.U_CT_CP_AD, 10)
+    npt.assert_array_almost_equal(aDControl.U_CT_CP_AD_ws, aDControl_saved.U_CT_CP_AD_ws, 10)
 
@@ -310,1 +311,75 @@ wfm = RANSLUT(dataset, site, wts)
         ellipsys_power, WS_eff_star, ct_star = get_Ellipsys_equivalent_output(sim_res, aDControl)
+
+
+def test_rans_conv_lut():
+    # move turbine 1 600 300
+    wt_x = [-250, 600, -500, 0, 500, -250, 250]
+    wt_y = [433, 300, 0, 0, 0, -433, -433]
+    wts = HornsrevV80()
+    site = UniformSite([1, 0, 0, 0], ti=0.075 * 0.8)
+    Dref = wts.diameter()
+    demo_lut = ptf('ranslut/V80_ranslut_demo.nc',
+                   '846213eb655255f6e2201a47c2406f9e77f243f369398cb389bf7320b457dea8')
+    # Fit Gaussian deficits to hub height profiles for WeightedSum superposition
+    demo_convlut = 'V80_gaussian.nc'
+    fit_gauss(demo_lut, 'V80_gaussian.nc', ymin=-10, ymax=10)
+    wfm = RANSLUT(demo_lut, site, wts, convlut=demo_convlut, rotorAvgModel=GQGridRotorAvg(4, 3))
+
+    simres, _ = timeit(wfm.__call__, verbose=0, line_profile=0,
+                       profile_funcs=[GridInterpolator.__call__])(x=wt_x, y=wt_y, wd=[30], ws=[10])
+    npt.assert_array_almost_equal(simres.WS_eff_ilk.flatten(), [9.97150681, 9.93711854, 7.97991956, 9.99272109, 9.20975831,
+                                                                7.98196409, 7.46986422])
+    npt.assert_array_almost_equal(simres.ct_ilk.flatten(), [0.7933989, 0.79388034, 0.80597992, 0.7931019, 0.80406338,
+                                                            0.80598196, 0.80546986])
+    ds = xr.load_dataset(demo_lut)
+    aDControl = ADControl.from_lut(ds.deficits, wts, ws_cutin=4, ws_cutout=25, dws=1.0, cal_TI=0.06)
+    _, WS_eff_star, ct_star = get_Ellipsys_equivalent_output(simres, aDControl, rotorAvgModel=None)
+    npt.assert_array_almost_equal(WS_eff_star.flatten(), [7.69759095, 7.66977332, 6.1306985, 7.71475178, 7.08139117, 6.13226504, 5.7398877])
+    npt.assert_array_almost_equal(ct_star.flatten(), [1.33143388, 1.33273949, 1.36552774, 1.33062844, 1.36035502, 1.36553333, 1.36413403])
+
+
+def test_rans_conv_lut_multiturbine():
+    # Test with provided control file
+    demo_lut = ptf('ranslut/V80_ranslut_demo.nc',
+                   '846213eb655255f6e2201a47c2406f9e77f243f369398cb389bf7320b457dea8')
+    ds = xr.load_dataset(demo_lut)
+
+    # Test multi lut (using the same file for now with 50% wake deficit)
+    lut_V120 = ds.copy(deep=True)
+    lut_V120.deficits[:] *= .5
+
+    # Fit Gaussian deficits to hub height profiles for WeightedSum superposition
+    demo_convlut = 'V80_gaussian.nc'
+    fit_gauss(demo_lut, 'V80_gaussian.nc', ymin=-10, ymax=10, xfits=np.linspace(0.0, 100.0, 11))
+    lut_V120.to_netcdf('V120.nc')
+    fit_gauss('V120.nc', 'V120_gaussian.nc', ymin=-10, ymax=10, xfits=np.linspace(0.0, 100.0, 11))
+
+    v80 = HornsrevV80()
+    v120 = WindTurbine('V120', 120, 70, powerCtFunction=PowerCtTabular(
+        hornsrev1.power_curve[:, 0], hornsrev1.power_curve[:, 1], 'w', hornsrev1.ct_curve[:, 1]))
+    wts = WindTurbine.from_WindTurbines([v80, v120])
+
+    wfm = RANSLUT([ds, lut_V120], UniformSite(ti=0.075 * 0.8), wts, convlut=['V80_gaussian.nc', 'V120_gaussian.nc'], rotorAvgModel=GQGridRotorAvg(4, 3))
+    type_i = np.array([0, 0, 1, 1])
+    sim_res = wfm([0, 500, 1000, 1500], [0, 0, 0, 0], type=type_i, wd=[90, 270], ws=10.0)
+    if 0:
+        axes = plt.subplots(2, 1)[1]
+        sim_res.flow_map(wd=90).plot_wake_map(ax=axes[0], levels=np.linspace(3, 10.5))
+        sim_res.flow_map(wd=270).plot_wake_map(ax=axes[1], levels=np.linspace(3, 10.5))
+        plt.show()
+
+    WS_eff_ref = [[7.3973466, 7.90998823, 8.56042535, 9.97638533], [9.97756182, 7.94966889, 7.78868653, 7.50030898]]
+    ct_ref = [[0.80539735, 0.80590999, 0.80656043, 0.79333061], [0.79331413, 0.80594967, 0.80578869, 0.80550031]]
+    power_ref = [[0.5537738, 0.67475722, 0.8641276, 1.33285294], [1.33325883, 0.68412186, 0.64613002, 0.57807292]]
+    TI_eff_ref = [[0.26789829, 0.21534878, 0.14177617, 0.06003113], [0.0600239, 0.13786305, 0.18521198, 0.23600163]]
+    npt.assert_array_almost_equal(sim_res.WS_eff_ilk.squeeze().T, WS_eff_ref, 6)
+    npt.assert_array_almost_equal(sim_res.ct_ilk.squeeze().T, ct_ref, 6)
+    npt.assert_array_almost_equal(sim_res.power_ilk.squeeze().T * 1e-6, power_ref, 6)
+    npt.assert_array_almost_equal(sim_res.TI_eff.values.squeeze().T, TI_eff_ref, 6)
+
+    aDControl = ADControl.from_lut([ds, lut_V120], wts, ws_cutin=[3, 4], ws_cutout=[26, 25], dws=1.0, cal_TI=0.06)
+    _, WS_eff_star, ct_star = get_Ellipsys_equivalent_output(sim_res, aDControl, rotorAvgModel=None)
+    WS_eff_ref = [[5.68432379, 6.07711624, 7.50664894, 8.76770315], [7.70248902, 6.10752006, 6.83079781, 6.57818957]]
+    ct_star_ref = [[1.36393587, 1.36533666, 1.0488932, 1.0271518], [1.33120399, 1.36544508, 1.04761978, 1.04714416]]
+    npt.assert_array_almost_equal(WS_eff_star.squeeze().T, WS_eff_ref, 6)
+    npt.assert_array_almost_equal(ct_star.squeeze().T, ct_star_ref, 6)

py_wake/utils/model_utils.py

@@ -181,2 +181,4 @@ import inspect
     from py_wake.superposition_models import LinearSum
+    from py_wake.deficit_models.deficit_model import XRLUTDeficitModel
+    from py_wake.deficit_models.rans_lut import RANSLUTDeficit, RANSLUTConvDeficit
     from py_wake.turbulence_models.rans_lut_turb import RANSLUTTurbulence
@@ -194,6 +196,6 @@ from py_wake.turbulence_models.turbulence_model import XRLUTTurbulenceModel
                          [RotorAvgModel], NOJDeficit),
-        "WakeDeficitModel": ([ConvectionDeficitModel, XRLUTDeficitModel, RANSLUTDeficit], [RotorAvgModel], NOJDeficit),
+        "WakeDeficitModel": ([ConvectionDeficitModel, XRLUTDeficitModel, RANSLUTDeficit, RANSLUTConvDeficit], [RotorAvgModel], NOJDeficit),
         "RotorAvgModel": ([NodeRotorAvgModel], [], None),
         "SuperpositionModel": ([], [], LinearSum),
-        "BlockageDeficitModel": ([XRLUTDeficitModel, RANSLUTDeficit], [], None),
+        "BlockageDeficitModel": ([XRLUTDeficitModel, RANSLUTDeficit], [RANSLUTConvDeficit], None),
         "DeflectionModel": ([DeflectionIntegrator], [], None),
@@ -200,0 +202,0 @@ "TurbulenceModel": ([XRLUTTurbulenceModel, RANSLUTTurbulence], [], None),

py_wake/utils/rans_lut_utils.py

@@ -9,2 +9,5 @@ import xarray as xr
 from py_wake.utils.most import phi, psi, phi_eps
+import sys
+import matplotlib.pyplot as plt
+from scipy.optimize import curve_fit
 
@@ -33,11 +36,13 @@
 class ADControl():
-    def __init__(self, UAD_CTstar_CPstar_lst, generate_args=None):
-        self.U_CT_CP_AD = UAD_CTstar_CPstar_lst
+    def __init__(self, UAD_CTstar_CPstar_ws_lst, generate_args=None):
+        self.U_CT_CP_AD_ws = UAD_CTstar_CPstar_ws_lst
         self.generate_args = generate_args
         self.fint_CTstar = []
         self.fint_CPstar = []
-        for UAD, CTstar, CPstar in UAD_CTstar_CPstar_lst:
-            assert len(UAD) == len(CTstar) == len(CPstar)
+        self.fint_UAD = []
+        for UAD, CTstar, CPstar, ws in UAD_CTstar_CPstar_ws_lst:
+            assert len(UAD) == len(CTstar) == len(CPstar) == len(ws)
             self.fint_CPstar.append(interp1d(UAD, CPstar, fill_value=0.0, bounds_error=False))
             self.fint_CTstar.append(interp1d(UAD, CTstar, fill_value=0.0, bounds_error=False))
+            self.fint_UAD.append(interp1d(ws, UAD, fill_value=0.0, bounds_error=False))
 
@@ -48,5 +53,5 @@ @staticmethod
         # the disk averaged velocity UAD provided by the user
-        UAD_CTstar_CPstar_lst = [np.genfromtxt(ADcontrolfile, skip_header=True).T[:3]
-                                 for ADcontrolfile in ADcontrolfiles]
-        return ADControl(UAD_CTstar_CPstar_lst)
+        UAD_CTstar_CPstar_ws_lst = [np.genfromtxt(ADcontrolfile, skip_header=True).T[:4]
+                                    for ADcontrolfile in ADcontrolfiles]
+        return ADControl(UAD_CTstar_CPstar_ws_lst)
 
@@ -57,3 +62,3 @@ def save(self, folder='.'):
 
-        for type_, U_CT_CP_AD in enumerate(self.U_CT_CP_AD):
+        for type_, U_CT_CP_AD_ws in enumerate(self.U_CT_CP_AD_ws):
             Dref = windTurbines.diameter(type_)
@@ -67,3 +72,3 @@ cal_TI = cal_TI
                 #        (0.95 * U_CT_CP_AD[0, 0], U_CT_CP_AD[0, 1], U_CT_CP_AD[0, 2], 0.95 * ws[0]))
-                for (U, CT, CP), ws_ in zip(U_CT_CP_AD.T, ws):
+                for U, CT, CP, ws_ in U_CT_CP_AD_ws.T:
                     f.write('%16.14f %16.14f %16.14f %16.14f\n' % (U, CT, CP, ws_))
@@ -147,7 +152,8 @@ # f.write('%16.14f %16.14f %16.14f %16.14f\n' % (ws[-1], 0.0, 0.0, ws[-1]))
             ratio = ws / UAD
-            UAD_CTstar_CPstar = np.array([UAD, cts * ratio ** 2, cps * ratio ** 3])
-            UAD_CTstar_CPstar[0, 0] = 0.95 * UAD[1]
-            UAD_CTstar_CPstar[:, -1] = ws[-1], 0.0, 0.0
+            UAD_CTstar_CPstar_ws = np.array([UAD, cts * ratio ** 2, cps * ratio ** 3, ws])
+            UAD_CTstar_CPstar_ws[0, 0] = 0.95 * UAD[1]
+            UAD_CTstar_CPstar_ws[3, 0] = 0.95 * ws[1]
+            UAD_CTstar_CPstar_ws[:, -1] = ws[-1], 0.0, 0.0, ws[-1]
 
-            return UAD_CTstar_CPstar
+            return UAD_CTstar_CPstar_ws
         return ADControl([get_UAD_CTstar_CPstar(t) for t in types], (windTurbines, cal_TI, ws_lst))
@@ -158,2 +164,13 @@
                                    rotorAvgModel=EllipSysPolygonRotorAvg(n_r=9, n_theta=32)):
+    if rotorAvgModel is not None:
+        # Use flow map to integrate rotor-averaged wind speed
+        return get_Ellipsys_equivalent_output_flowmap(sim_res, aDControl, flowmap_maxpoints=flowmap_maxpoints,
+                                                      rotorAvgModel=rotorAvgModel)
+    else:
+        # Interpolate rotor-averaged wind speed, applicable when a rotorAvgModel is used inside RANSLUT model
+        return get_Ellipsys_equivalent_output_simple(sim_res, aDControl)
+
+
+def get_Ellipsys_equivalent_output_flowmap(sim_res, aDControl, flowmap_maxpoints=None,
+                                           rotorAvgModel=EllipSysPolygonRotorAvg(n_r=9, n_theta=32)):
     """
@@ -249,3 +266,3 @@ Recalculate wt power and thrust by the integrating rotor averaged wind speed
 
-    U_tab_lst = [U_CT_CP_AD[0, :] for U_CT_CP_AD in aDControl.U_CT_CP_AD]
+    U_tab_lst = [U_CT_CP_AD[0, :] for U_CT_CP_AD in aDControl.U_CT_CP_AD_ws[:3]]
 
@@ -278,1 +295,86 @@ Prated = [wts.power(U_tab,
     return power_ilk, WS_eff_star_ilk, ct_star_ilk
+
+
+def get_Ellipsys_equivalent_output_simple(sim_res, aDControl):
+    """
+    Calculate CTstar, the thrust coefficient based rotor averaged wind speed, U_AD, without
+    flow_map integration, by interpolating a U_AD, CTstar = f(WS_eff) relation from RANS,
+    as done in EllipSys.
+    """
+    type_i = sim_res.type.values
+    I = sim_res.WS_eff_ilk.shape[0]
+
+    def get_UAD_ct(i, t):
+        UAD = aDControl.fint_UAD[t](sim_res.WS_eff_ilk[i, ::])
+        ct_star = aDControl.fint_CTstar[t](UAD)
+        return UAD, ct_star
+
+    if np.all(type_i == 0):
+        WS_eff_star_ilk, ct_star_ilk = map(np.array, get_UAD_ct(slice(None), 0))
+    else:
+        WS_eff_star_ilk, ct_star_ilk = map(np.array, zip(*[get_UAD_ct(i, type_i[i]) for i in range(I)]))
+
+    # Set ct to zero for wts that were originally shut down
+    operating = sim_res.power_ilk > 0
+    ct_star_ilk *= operating
+    return sim_res.power_ilk, WS_eff_star_ilk, ct_star_ilk
+
+
+def gauss(x, a, x0, sigma):
+    return a * np.exp(-(x - x0) ** 2 / (2 * sigma ** 2))
+
+
+def fit_gauss(path_LUT, database_name, ymin=-10, ymax=10, xfits=None):
+    # Fit a Gaussian for each single wake shape and store sigma and magnitude
+    ds = xr.open_dataset(path_LUT).interp(z=0.0)
+    ds = ds.where(ds.y >= ymin, drop=True).where(ds.y <= ymax, drop=True)
+    if xfits is None:
+        xmin = 2.0
+        xmax = 100.0
+        xfits = np.insert(np.append(ds.x.where(ds.x > xmin, drop=True).where(ds.x < xmax, drop=True), xmax), 0, 0.0)
+    centerline_deficit = np.zeros((len(ds.ct), len(ds.ti), len(xfits) + 1))
+    sigma = np.zeros((len(ds.ct), len(ds.ti), len(xfits) + 1))
+    for k in range(len(ds.ct)):
+        for j in range(len(ds.ti)):
+            deficit = ds.isel(ct=k, ti=j)
+            for i in range(len(xfits)):
+                u = deficit.interp(x=xfits[i])
+                try:
+                    # Substract outer values to make fit possible when velocity deficit at hub height is completely negative
+                    # in the very far wake due to vertical wake center displacement
+                    popt, pcov = curve_fit(gauss, u.y, u.deficits - u.deficits.isel(y=0), p0=[1.0, 0.0, 1.0], maxfev=800)
+                except Exception:  # pragma: no cover
+                    plt.plot(u.y, u.deficits, 'b+:', label='data')
+                    plt.plot(u.y, gauss(u.y, *popt), 'ro:', label='fit')
+                    plt.savefig('fittest.pdf')
+                    sys.exit()
+                centerline_deficit[k, j, i] = popt[0]
+                sigma[k, j, i] = popt[2]
+        # Limit centerline deficit in near wake to 1D momentum value, 2*a_x, since fitting a gaussian to a non-gaussian near wake
+        # is prone to overpredicting the deficit
+        centerline_deficit[k, :, :] = np.minimum(centerline_deficit[k, :, :], (1.0 - np.sqrt(1.0 - np.minimum(ds.ct.isel(ct=k).values, 1.0)))[na, na])
+        # Likewise, limit deficit at x=0 with 1D momentum value, a_x
+        centerline_deficit[k, :, 0] = np.minimum(centerline_deficit[k, :, 0], (0.5 * (1.0 - np.sqrt(1.0 - np.minimum(ds.ct.isel(ct=k).values, 1.0))))[na])
+
+    centerline_deficit[:, :, -1] = centerline_deficit[:, :, -2]
+    sigma[:, :, -1] = sigma[:, :, -2]
+    xfits = np.append(xfits, 10**10)
+    dsout = xr.Dataset({'centerline_deficit': xr.DataArray(centerline_deficit,
+                        [('ct', ds.ct.values),
+                         ('ti', ds.ti.values),
+                         ('x', xfits)],
+                        attrs={'long_name': 'Normalized streamwise centerline deficit from fitted Gaussian', 'units': '[-]'}),
+                       'sigma': xr.DataArray(sigma,
+                        [('ct', ds.ct.values),
+                         ('ti', ds.ti.values),
+                         ('x', xfits)],
+                        attrs={'long_name': 'Standard deviation normalized by rotor diameter from fitted Gaussian', 'units': '[-]'})})
+    encoding = {}
+    vars = ['centerline_deficit', 'sigma']
+    for var in vars:
+        encoding[var] = {'dtype': 'float32'}
+    vars = ['x', 'ct', 'ti']
+    for var in vars:
+        encoding[var] = {'dtype': 'float64'}
+    # Write netcdf file
+    dsout.to_netcdf(database_name, encoding=encoding)

py_wake/version.py

@@ -31,5 +31,5 @@ # file generated by setuptools-scm
 
-__version__ = version = '2.6.16'
-__version_tuple__ = version_tuple = (2, 6, 16)
+__version__ = version = '2.6.17'
+__version_tuple__ = version_tuple = (2, 6, 17)
 
 __commit_id__ = commit_id = None

py_wake/wind_farm_models/external_wind_farm_models.py

@@ -8,3 +8,2 @@ from abc import ABC, abstractmethod
 from scipy.interpolate import RegularGridInterpolator as RGI
-from tqdm import tqdm
 
@@ -153,3 +152,3 @@ from py_wake.flow_map import Points
         from py_wake.site.streamline_distance import StreamlineDistance
-        return cls(name, ds, windFarmModel.windTurbines, wt_x, wt_y,
+        return cls(name, ds.sortby('wd'), windFarmModel.windTurbines, wt_x, wt_y,
                    relative_distance=not isinstance(windFarmModel.site.distance, StreamlineDistance))
@@ -221,3 +220,3 @@
         from py_wake.site.streamline_distance import StreamlineDistance
-        return cls(name, ds, windFarmModel.windTurbines, wt_x, wt_y,
+        return cls(name, ds.sortby('wd'), windFarmModel.windTurbines, wt_x, wt_y,
                    relative_distance=not isinstance(windFarmModel.site.distance, StreamlineDistance))
@@ -224,0 +223,0 @@

RECORD

@@ -8,6 +8,6 @@ docs/__init__.py,sha256=47DEQpj8HBSa-_TImW-5JCeuQeRkm5NMpJWZG3hSuFU,0
 py_wake/flow_map.py,sha256=pZCV-mbPDtVLuvdGejKJhUnR52VorJ7fobP10Oj7ZL4,22441
-py_wake/superposition_models.py,sha256=mlZSy-T6YiAhdSgAu-reP1YECRsPpgHh8GDCfPAsZ6w,10124
-py_wake/version.py,sha256=w6Wi-Vv2osckU1Br_lqmghSD8v0m4QyQ_pMHdGIEQi4,706
+py_wake/superposition_models.py,sha256=wJ-JJElwdWiNOwv1v0NXjbTV1ZbbBSmC-GJziO1YPJs,10669
+py_wake/version.py,sha256=NauhVIERykrLkuAo4rI2RyJQXOSgmqEnQesJmE0DWyY,706
 py_wake/deficit_models/__init__.py,sha256=w1S4F3FmAgtwqGO90fASYoqxOb6SsgjbasiD3V8_1uk,709
-py_wake/deficit_models/deficit_model.py,sha256=iLb9wUO8wDeYAEBAoAlS2W4Q4B16ZSO0L4DhVq3NxYQ,10280
+py_wake/deficit_models/deficit_model.py,sha256=ly8A6pZQ6kt39fzsAz9OVw_5ohr7b26I5iHRYni2QrM,10466
 py_wake/deficit_models/fuga.py,sha256=pZ9wUI5VBPcJM-2FWGUVYwQIQ6oWKFlOw0qSyu6V4hI,13907
@@ -20,3 +20,3 @@ py_wake/deficit_models/gaussian.py,sha256=_ntryBrucuo9EMSLnyJH3dt9mM9SnB0DxP9wzH_w7Zo,36119
 py_wake/deficit_models/rankinehalfbody.py,sha256=idu2_roqZQWDzkcvRMhUdv4-KWuVUlCgSsxb2N2wbm8,5408
-py_wake/deficit_models/rans_lut.py,sha256=0zSDrQSgtjjhjzo0-Oiz5UFRN2vecJrT1hVfkq6c_jQ,3451
+py_wake/deficit_models/rans_lut.py,sha256=-sFLM53WbH0N8NgF9bdVa9mDxVZC2STLkdVlHK0Yyg8,7918
 py_wake/deficit_models/rathmann.py,sha256=CO2zQXCrjiyoh9Q3Pds8ddHKL_G5M1DEOLnGz9SeOgI,15368
@@ -289,3 +289,3 @@ py_wake/deficit_models/selfsimilarity.py,sha256=l4hFh9l4WkWlJxgEhkXMnGSmiSAXCQlg_7_VaYEnmoc,13534
 py_wake/site/jit_streamline_distance.py,sha256=GJo2mn_bTjuGX3-gSGkv1dkRfo7hGDVKDmefRWYk5o8,159
-py_wake/site/shear.py,sha256=HVrUmwzUiTAEKGP6eVJQKOODv1FCbI30CJaPSMj4H0s,4248
+py_wake/site/shear.py,sha256=7Nzg3wBTAHIknX_yLhq5VNZGNIy0oR5rJfk2BWb5wMY,4276
 py_wake/site/streamline_distance.py,sha256=01NdwM__ITY0w4kWavUx1arKxCvtc42cEteugQhN7mo,4518
@@ -310,3 +310,3 @@ py_wake/site/wasp_grid_site.py,sha256=d-vonrkXMuVVgbq9OG0ri3XHYekdDfvfjdndkigCOrQ,10227
 py_wake/tests/test_deficit_models/test_noj.py,sha256=qJVv3q3Lp_QILFn-6xdslWevvR4Dmnx0AWj39XCIdXw,5495
-py_wake/tests/test_deficit_models/test_rans_lut.py,sha256=gpskfpoK6HBpW-8dFB5dHgD5NKtrZle_2q8ftx8fmlM,17135
+py_wake/tests/test_deficit_models/test_rans_lut.py,sha256=9u94ZVzzGIRZEdJ0CftcN8LX1htqoXrw_U8y35Hi2eM,21881
 py_wake/tests/test_deficit_models/test_selfsimilarity.py,sha256=mI5DUstBbaC36JQ-eCOMrE88jubvZjb3tp_DUQQzU2Y,4446
@@ -407,3 +407,3 @@ py_wake/tests/test_deflection_models/__init__.py,sha256=47DEQpj8HBSa-_TImW-5JCeuQeRkm5NMpJWZG3hSuFU,0
 py_wake/utils/maps.py,sha256=cvz1j14Is_EdsXX4wkautCXm_OcDaSwDk4a6JJdhGs4,593
-py_wake/utils/model_utils.py,sha256=UP3H4XlAQsE45zl3QFje1zroF7OaqCUDEm6TZ81OrEc,15521
+py_wake/utils/model_utils.py,sha256=BdM_RoEFXfxP2Dha2TUPvqJalhKznBPdLouePR1VkTU,15713
 py_wake/utils/most.py,sha256=9DEZ5C06qdttfGdfUCyaFCGt-YT8xPnQ3URIuHpbcdE,1208
@@ -415,3 +415,3 @@ py_wake/utils/notebook_extensions.py,sha256=E0UOiQ1AAMaIbnnzsyj-gGv8oRWD8wAvHsj2BN10LXU,704
 py_wake/utils/profiling.py,sha256=g6aDNtZg_7lcxUg4WdAawGMysaoZW87iPyTWwiPvWdA,4668
-py_wake/utils/rans_lut_utils.py,sha256=-zs25tCdg2heaoCjyEDdk0vSDTImdhjtTpbDH_YEClI,12423
+py_wake/utils/rans_lut_utils.py,sha256=BjEurt2nM7ZO2VLQN1bDdaCtig1N3ksqafYcxlOTQxQ,17596
 py_wake/utils/streamline.py,sha256=QSfp_kLVxihTlWq3P0AHvpf_OS-QUIj17xoNhhI8QXk,2835
@@ -427,4 +427,4 @@ py_wake/utils/tensorflow_surrogate_utils.py,sha256=QdIPqxoYiuXFvieHQyP-TE6gS8syXbxoUZ1liCBFil0,7198
 py_wake/wind_farm_models/__init__.py,sha256=cl0nC9vqhXgpRuaCFkQpDET2Yyzbwnd5Kv2wgZHUI6E,137
-py_wake/wind_farm_models/engineering_models.py,sha256=Mq3kh5aFfO4YSrBHzg3ukPN-RBZb8XKBMT9vJESptp4,63776
-py_wake/wind_farm_models/external_wind_farm_models.py,sha256=QRSwdq45J4hwmYANBCv87HFGAho04_gkTWY9rUUUazM,12587
+py_wake/wind_farm_models/engineering_models.py,sha256=co-te7NzyLkz3PchizBj2_-YGU6wFlWl6ZbnTgKSLdg,64395
+py_wake/wind_farm_models/external_wind_farm_models.py,sha256=bFAmIJ7bRHEZEIo_wyAabefw4OYNKZ1-p5OoSAca7Ok,12591
 py_wake/wind_farm_models/minimalistic_wind_farm_model.py,sha256=9BepPjAo_WDMcG4feUQo7YDY2WSFJOpNMmzgT0yNBn4,11160
@@ -438,6 +438,6 @@ py_wake/wind_farm_models/wind_farm_model.py,sha256=8rQ0tbfZDMs8_9Ag_c1UQz28O7xF7tEBS43wey_BWfs,42573
 py_wake/wind_turbines/wind_turbines_deprecated.py,sha256=HpNmBR8CJL4-8JBaygDI0t086qfw5bR2DOQI8Ox4AZ4,6250
-py_wake-2.6.16.dist-info/licenses/LICENSE,sha256=XE2CGPqQgzSXqIajXpAVYJ5SRNmaWOIeMePK6MocsuY,1084
-py_wake-2.6.16.dist-info/METADATA,sha256=gVr6wRfNyRfmVfe0Mbzgi9H8-1EDANsVPMZO5O-SI5Q,3760
-py_wake-2.6.16.dist-info/WHEEL,sha256=_zCd3N1l69ArxyTb8rzEoP9TpbYXkqRFSNOD5OuxnTs,91
-py_wake-2.6.16.dist-info/top_level.txt,sha256=GsaXU4YwyMkZZ6dkb4h0FMc5RaLIT2Qns_YoScKoXdk,20
-py_wake-2.6.16.dist-info/RECORD,,
+py_wake-2.6.17.dist-info/licenses/LICENSE,sha256=XE2CGPqQgzSXqIajXpAVYJ5SRNmaWOIeMePK6MocsuY,1084
+py_wake-2.6.17.dist-info/METADATA,sha256=bA0FK31ceJP7fXVWs2kbF5heDUstSVqr3ZgQmmbYy_E,3760
+py_wake-2.6.17.dist-info/WHEEL,sha256=_zCd3N1l69ArxyTb8rzEoP9TpbYXkqRFSNOD5OuxnTs,91
+py_wake-2.6.17.dist-info/top_level.txt,sha256=GsaXU4YwyMkZZ6dkb4h0FMc5RaLIT2Qns_YoScKoXdk,20
+py_wake-2.6.17.dist-info/RECORD,,

Alert delta unavailable

Currently unable to show alert delta for PyPI packages.

No metric changes