From f58e8bf7719b107e73bfe86284adccb4a71a0447 Mon Sep 17 00:00:00 2001
From: "Mads M. Pedersen" <mmpe@dtu.dk>
Date: Thu, 6 Mar 2025 12:03:09 +0000
Subject: [PATCH] fix fuga wake radius (before it was fixed to D) + switch
bounds from limit to check in Fuga + tests
---
py_wake/deficit_models/fuga.py | 51 +++++++++++++++----
.../test_blockage_models.py | 3 +-
.../test_deficit_models.py | 4 +-
.../tests/test_deficit_models/test_fuga.py | 44 +++++++++++++++-
4 files changed, 87 insertions(+), 15 deletions(-)
diff --git a/py_wake/deficit_models/fuga.py b/py_wake/deficit_models/fuga.py
index 952da611d..8b8af5ecc 100644
--- a/py_wake/deficit_models/fuga.py
+++ b/py_wake/deficit_models/fuga.py
@@ -19,7 +19,7 @@ import os
class Fuga(PropagateDownwind, DeprecatedModel):
def __init__(self, LUT_path, site, windTurbines,
rotorAvgModel=None, deflectionModel=None, turbulenceModel=None,
- smooth2zero_x=None, smooth2zero_y=None, remove_wriggles=False):
+ bounds='limit', smooth2zero_x=None, smooth2zero_y=None, remove_wriggles=False):
"""
Parameters
----------
@@ -39,7 +39,7 @@ class Fuga(PropagateDownwind, DeprecatedModel):
Model describing the amount of added turbulence in the wake
"""
PropagateDownwind.__init__(self, site, windTurbines,
- wake_deficitModel=FugaDeficit(LUT_path,
+ wake_deficitModel=FugaDeficit(LUT_path, bounds=bounds,
smooth2zero_x=smooth2zero_x,
smooth2zero_y=smooth2zero_y,
remove_wriggles=remove_wriggles),
@@ -50,7 +50,7 @@ class Fuga(PropagateDownwind, DeprecatedModel):
class FugaBlockage(All2AllIterative, DeprecatedModel):
def __init__(self, LUT_path, site, windTurbines, rotorAvgModel=None,
- deflectionModel=None, turbulenceModel=None, convergence_tolerance=1e-6, remove_wriggles=False):
+ deflectionModel=None, turbulenceModel=None, convergence_tolerance=1e-6, bounds='limit', remove_wriggles=False):
"""
Parameters
----------
@@ -69,7 +69,7 @@ class FugaBlockage(All2AllIterative, DeprecatedModel):
turbulenceModel : TurbulenceModel
Model describing the amount of added turbulence in the wake
"""
- fuga_deficit = FugaDeficit(LUT_path, remove_wriggles=remove_wriggles)
+ fuga_deficit = FugaDeficit(LUT_path, bounds=bounds, remove_wriggles=remove_wriggles)
All2AllIterative.__init__(self, site, windTurbines, wake_deficitModel=fuga_deficit,
rotorAvgModel=rotorAvgModel, superpositionModel=LinearSum(),
deflectionModel=deflectionModel, blockage_deficitModel=fuga_deficit,
@@ -135,9 +135,33 @@ class FugaMultiLUTDeficit(XRLUTDeficitModel):
dims = ['d_h', 'zeta0', 'zi', 'z0']
return dims
- def wake_radius(self, D_src_il, dw_ijlk, **_):
- # Set at twice the source radius for now
- return np.zeros_like(dw_ijlk) + D_src_il[:, na, :, na]
+ def wake_radius(self, dw_ijlk, D_src_il, h_ilk, TI_ilk, hcw_ijlk=None, z_ijlk=None, **kwargs):
+ z_ijlk = h_ilk[:, na]
+ lim = self.da.y.max().item()
+
+ def get_mdu(hcw_ijlk):
+ hcw_ijlk = np.nan_to_num(np.clip(hcw_ijlk, -lim, lim), nan=lim)
+ return self._calc_mdu(D_src_il=D_src_il, dw_ijlk=dw_ijlk,
+ hcw_ijlk=hcw_ijlk, h_ilk=h_ilk, z_ijlk=z_ijlk, TI_ilk=TI_ilk, **kwargs)
+ mdu_target = get_mdu(dw_ijlk * 0) * np.exp(-2) # corresponding to 2 sigma
+
+ def get_err(hcw):
+ return get_mdu(hcw) - mdu_target
+
+ def get_wake_radius(wake_radius_ijlk):
+ # Newton Raphson
+ for _ in range(4):
+ err = get_err(wake_radius_ijlk)
+ derr = get_err(wake_radius_ijlk + 1) - err
+ wake_radius_ijlk = wake_radius_ijlk - err / derr
+ return np.abs(wake_radius_ijlk)
+ wake_radius_ijlk = get_wake_radius(D_src_il[:, na, :, na]) # diameter as initial guess
+
+ if np.any(kwargs.get('yaw_ilk', [0])) and 'UT' in self.da.variable_names:
+ # mean of positive and negative side
+ wake_radius_ijlk = (wake_radius_ijlk + get_wake_radius(-D_src_il[:, na, :, na]))
+
+ return wake_radius_ijlk
def calc_deficit(self, WS_ilk, WS_eff_ilk, dw_ijlk, hcw_ijlk, z_ijlk, ct_ilk, D_src_il, **kwargs):
# bypass XRLUTDeficitModel.calc_deficit
@@ -145,19 +169,24 @@ class FugaMultiLUTDeficit(XRLUTDeficitModel):
self._calc_layout_terms(dw_ijlk=dw_ijlk, hcw_ijlk=hcw_ijlk, z_ijlk=z_ijlk, D_src_il=D_src_il, **kwargs)
return self.mdu_ijlk * (ct_ilk * WS_eff_ilk**2 / WS_ilk)[:, na]
- def _calc_layout_terms(self, **kwargs):
+ def _calc_mdu(self, **kwargs):
# 0 = UL
variables0 = np.reshape(0, np.ones_like(kwargs['dw_ijlk'].shape).astype(int))
- self.mdu_ijlk = XRLUTDeficitModel.calc_deficit(self, **kwargs, variables=variables0)
+ mdu_ijlk = XRLUTDeficitModel.calc_deficit(self, **kwargs, variables=variables0)
if 'yaw_ilk' in kwargs:
theta_yaw_ijlk = gradients.deg2rad(kwargs['yaw_ilk'])[:, na]
if list(self.da.variable_names) == ['UL', 'UT']:
# 1 = UT
mdUT_ijlk = XRLUTDeficitModel.calc_deficit(self, **kwargs, variables=variables0 + 1)
mdUT_ijlk = np.negative(mdUT_ijlk, out=mdUT_ijlk, where=kwargs['hcw_ijlk'] >= 0)
- self.mdu_ijlk = self.mdu_ijlk * np.cos(theta_yaw_ijlk) + np.sin(theta_yaw_ijlk) * mdUT_ijlk
+ mdu_ijlk = mdu_ijlk * np.cos(theta_yaw_ijlk) + np.sin(theta_yaw_ijlk) * mdUT_ijlk
else:
- self.mdu_ijlk *= np.cos(theta_yaw_ijlk)
+ mdu_ijlk *= np.cos(theta_yaw_ijlk)
+ return mdu_ijlk
+
+ def _calc_layout_terms(self, **kwargs):
+ self.mdu_ijlk = self._calc_mdu(**kwargs)
+ return self.mdu_ijlk
def get_input(self, D_src_il, TI_ilk, h_ilk, dw_ijlk, hcw_ijlk, z_ijlk, **kwargs):
user = {'zeta0': lambda: kwargs['zeta0_ilk'][:, na],
diff --git a/py_wake/tests/test_deficit_models/test_blockage_models.py b/py_wake/tests/test_deficit_models/test_blockage_models.py
index 50f767250..3db379c4c 100644
--- a/py_wake/tests/test_deficit_models/test_blockage_models.py
+++ b/py_wake/tests/test_deficit_models/test_blockage_models.py
@@ -189,8 +189,7 @@ def test_All2AllIterative_all_blockage_DeficitModels_with_RotorAvg(deficitModel)
sim_res = wf_model([0, 500, 1000, 1500], [0, 0, 0, 0], yaw=0, wd=270, ws=10)
if 0:
- sim_res.flow_map(
- XYGrid(x=np.linspace(-200, 2000, 100))).plot_wake_map()
+ sim_res.flow_map(XYGrid(x=np.linspace(-200, 2000, 100))).plot_wake_map()
plt.show()
diff --git a/py_wake/tests/test_deficit_models/test_deficit_models.py b/py_wake/tests/test_deficit_models/test_deficit_models.py
index d3839d3c9..cdb9a3b75 100644
--- a/py_wake/tests/test_deficit_models/test_deficit_models.py
+++ b/py_wake/tests/test_deficit_models/test_deficit_models.py
@@ -254,7 +254,8 @@ def test_deficitModel_wake_map(deficitModel, ref):
(TurboNOJDeficit(), [99.024477, 61.553917, 123.107833, 92.439673, 97.034049]),
(BastankhahGaussianDeficit(), [83.336286, 57.895893, 115.791786, 75.266662, 83.336286]),
(IEA37SimpleBastankhahGaussianDeficit(), [103.166178, 67.810839, 135.621678, 103.166178, 103.166178]),
- (FugaDeficit(LUT_path=tfp + 'fuga/2MW/Z0=0.00408599Zi=00400Zeta0=0.00E+00.nc'), [100, 50, 100, 100, 100]),
+ (FugaDeficit(LUT_path=tfp + 'fuga/2MW/Z0=0.00408599Zi=00400Zeta0=0.00E+00.nc'),
+ [88.044245, 88.044245, 110.523363, 88.044245, 88.044245]),
(GCLDeficit(), [156.949964, 97.763333, 195.526667, 113.225695, 111.340236]),
(GCLLocalDeficit(), [156.949964, 97.763333, 195.526667, 113.225695, 111.340236]),
(ZongGaussianDeficit(eps_coeff=0.35), [91.15734, 66.228381, 132.456762, 94.90156, 79.198215]),
@@ -272,6 +273,7 @@ def test_wake_radius(deficitModel, wake_radius_ref):
D_src_il=np.reshape([100, 50, 100, 100, 100], (5, 1)),
dw_ijlk=np.reshape([500, 500, 1000, 500, 500], (5, 1, 1, 1)),
WS_ilk=np.reshape([10, 10, 10, 10, 10], (5, 1, 1)),
+ h_ilk=np.reshape([70, 70, 70, 70, 70], (5, 1, 1)),
ct_ilk=np.reshape([.8, .8, .8, .4, .8], (5, 1, 1)),
TI_ilk=np.reshape([.1, .1, .1, .1, .05], (5, 1, 1)),
TI_eff_ilk=np.reshape([.1, .1, .1, .1, .05], (5, 1, 1)))[:, 0, 0, 0],
diff --git a/py_wake/tests/test_deficit_models/test_fuga.py b/py_wake/tests/test_deficit_models/test_fuga.py
index 6d9c074dc..b8bf4b8fa 100644
--- a/py_wake/tests/test_deficit_models/test_fuga.py
+++ b/py_wake/tests/test_deficit_models/test_fuga.py
@@ -18,7 +18,9 @@ from py_wake.wind_turbines._wind_turbines import WindTurbine, WindTurbines
from py_wake.utils.profiling import timeit
import warnings
from py_wake.utils import fuga_utils
-from py_wake.utils.fuga_utils import FugaXRLUT, FugaUtils
+from py_wake.utils.fuga_utils import FugaUtils
+from numpy import newaxis as na
+from scipy.optimize._minpack_py import curve_fit
def test_fuga():
@@ -469,6 +471,46 @@ def test_verify_FugaMultiLUT():
npt.assert_array_almost_equal(fm.WS_eff.sel(time=t).squeeze(), fm_ref.WS_eff.squeeze(), 5)
+@pytest.mark.parametrize('yaw', [0, 30, 60])
+def test_fuga_wake_radius(yaw):
+ wts = HornsrevV80()
+
+ path = tfp + 'fuga/2MW/Z0=0.00408599Zi=00400Zeta0=0.00E+00.nc'
+ site = UniformSite()
+ wfm = Fuga(path, site, wts)
+ sim_res = wfm([0], [0], wd=270, ws=10, yaw=yaw)
+
+ y = np.linspace(-300, 300, 500)
+ du = 10 - sim_res.flow_map(XYGrid(x=500, y=y)).WS_eff.squeeze()
+ wake_radius = wfm.wake_deficitModel.wake_radius(dw_ijlk=np.array([500])[na, :, na, na], D_src_il=wts.diameter([0])[:, na],
+ h_ilk=wts.hub_height([0])[:, na, na],
+ TI_ilk=np.array([0.1])[:, na, na]).squeeze()
+
+ def gauss(x, *p):
+ A, mu, sigma = p
+ return A * np.exp(-(x - mu)**2 / (2. * sigma**2))
+
+ # p0 is the initial guess for the fitting coefficients (A, mu and sigma above)
+ coeff, var_matrix = curve_fit(gauss, y, du, p0=[1, 0, 1])
+ A, mu, sigma = coeff
+ if 0:
+ x = np.linspace(-200, 1200)
+ sim_res.flow_map(XYGrid(x=x)).plot_wake_map()
+ wake_radius_map = wfm.wake_deficitModel.wake_radius(dw_ijlk=x[na, :, na, na], D_src_il=wts.diameter([0])[:, na],
+ h_ilk=wts.hub_height([0])[:, na, na],
+ TI_ilk=np.array([0.1])[:, na, na]).squeeze()
+ plt.plot(x, wake_radius_map)
+ plt.figure()
+ plt.plot(y, du, label='Fuga')
+ plt.axvline(wake_radius, label='wake radius')
+ plt.plot(y, gauss(y, *coeff), label=f'Gauss({coeff})')
+ plt.axvline(sigma * 2, ls='--', color='r', label='2*sigma')
+ plt.legend()
+ plt.show()
+ npt.assert_allclose(sigma * 2, wake_radius, rtol=0.02)
+
+ plt.show()
+
# def cmp_fuga_with_colonel():
# from py_wake.aep_calculator import AEPCalculator
#
--
GitLab