From 451a390d257538240bf83ad9847499c0e30b093d Mon Sep 17 00:00:00 2001
From: "Mads M. Pedersen" <mmpe@dtu.dk>
Date: Tue, 24 Aug 2021 13:09:18 +0000
Subject: [PATCH] fix All2AllIterative. Before layout terms was only
recalculated if deflection...
---
py_wake/deficit_models/deficit_model.py | 2 +-
py_wake/deflection_models/jimenez.py | 2 +-
.../test_deflection_models.py | 72 ++++++++++++++++++-
.../test_enginering_wind_farm_model.py | 29 ++++++--
.../wind_farm_models/engineering_models.py | 9 ++-
5 files changed, 102 insertions(+), 12 deletions(-)
diff --git a/py_wake/deficit_models/deficit_model.py b/py_wake/deficit_models/deficit_model.py
index 1e49a11ce..b022b161d 100644
--- a/py_wake/deficit_models/deficit_model.py
+++ b/py_wake/deficit_models/deficit_model.py
@@ -7,7 +7,7 @@ class DeficitModel(ABC):
deficit_initalized = False
def _calc_layout_terms(self, **_):
- pass
+ """Calculate layout dependent terms, which is not updated during simulation"""
@abstractmethod
def calc_deficit(self):
diff --git a/py_wake/deflection_models/jimenez.py b/py_wake/deflection_models/jimenez.py
index 1bb475f41..1c3d00dea 100644
--- a/py_wake/deflection_models/jimenez.py
+++ b/py_wake/deflection_models/jimenez.py
@@ -22,7 +22,7 @@ class JimenezWakeDeflection(DeflectionModel):
theta_ilk = np.sqrt(theta_yaw_ilk**2 + theta_tilt_ilk**2)
theta_deflection_ilk = np.arctan2(theta_tilt_ilk, theta_yaw_ilk)
denominator_ilk = np.cos(theta_ilk)**2 * np.sin(theta_ilk) * (ct_ilk / 2)
- nominator_ijxl = (1 + (self.beta / D_src_il)[:, na, na, :] * dw_ijxl)**2
+ nominator_ijxl = (1 + (self.beta / D_src_il)[:, na, na, :] * np.maximum(dw_ijxl, 0))**2
alpha = denominator_ilk[:, na, na] / nominator_ijxl[..., na]
deflection_ijlk = np.trapz(np.sin(alpha), dw_ijxl[..., na], axis=2)
self.hcw_ijlk = hcw_ijl[..., na] + deflection_ijlk * np.cos(theta_deflection_ilk[:, na])
diff --git a/py_wake/tests/test_deflection_models/test_deflection_models.py b/py_wake/tests/test_deflection_models/test_deflection_models.py
index 13794e280..66aa0136f 100644
--- a/py_wake/tests/test_deflection_models/test_deflection_models.py
+++ b/py_wake/tests/test_deflection_models/test_deflection_models.py
@@ -11,6 +11,8 @@ from py_wake.examples.data.hornsrev1 import V80
from py_wake.deflection_models.deflection_model import DeflectionModel
from py_wake.utils.model_utils import get_models
from py_wake.tests.test_files import tfp
+from py_wake.wind_farm_models.engineering_models import PropagateDownwind, All2AllIterative
+from py_wake.deficit_models.gaussian import IEA37SimpleBastankhahGaussianDeficit
@pytest.mark.parametrize('deflectionModel,dy10d', [
@@ -19,7 +21,8 @@ from py_wake.tests.test_files import tfp
((lambda: FugaDeflection(tfp + 'fuga/2MW/Z0=0.00408599Zi=00400Zeta0=0.00E+00/')), 0.37719329354768527),
((lambda: FugaDeflection(tfp + 'fuga/2MW/Z0=0.03000000Zi=00401Zeta0=0.00E+00/')), 0.32787746772608933),
])
-def test_deflection_model(deflectionModel, dy10d):
+def test_deflection_model_dy10d(deflectionModel, dy10d):
+ # center line deflection 10d downstream
site = IEA37Site(16)
x, y = [0], [0]
windTurbines = V80()
@@ -42,6 +45,73 @@ def test_deflection_model(deflectionModel, dy10d):
npt.assert_almost_equal(min_WS_line.interp(x=10 * D).item() / D, dy10d)
+@pytest.mark.parametrize('deflectionModel,dy', [
+ (JimenezWakeDeflection,
+ [2.0, 12.0, 20.0, 26.0, 2.0, -5.0, -11.0, -16.0, -0.0, 8.0, 15.0, 20.0]),
+ ((lambda: FugaDeflection(tfp + 'fuga/2MW/Z0=0.00001000Zi=00400Zeta0=0.00E+00/')),
+ [1.0, 6.0, 12.0, 18.0, 2.0, -0.0, -4.0, -7.0, -1.0, 2.0, 4.0, 7.0]),
+ # ((lambda: FugaDeflection(tfp + 'fuga/2MW/Z0=0.00408599Zi=00400Zeta0=0.00E+00/')), 0.37719329354768527),
+ ((lambda: FugaDeflection(tfp + 'fuga/2MW/Z0=0.03000000Zi=00401Zeta0=0.00E+00/')),
+ [1.0, 6.0, 11.0, 15.0, 2.0, -0.0, -3.0, -5.0, -1.0, 2.0, 4.0, 6.0]),
+])
+def test_deflection_model(deflectionModel, dy):
+ # center line deflection 10d downstream
+ site = IEA37Site(16)
+ x, y = [0, 400, 800], [0, 0, 0]
+ windTurbines = V80()
+ D = windTurbines.diameter()
+ wfm = PropagateDownwind(site, windTurbines, IEA37SimpleBastankhahGaussianDeficit(),
+ deflectionModel=deflectionModel())
+
+ yaw = [-30, 30, -30]
+
+ sim_res = wfm(x, y, yaw=yaw, wd=270, ws=10)
+ fm = sim_res.flow_map(XYGrid(x=np.arange(-D, 15 * D + 10, 10)))
+ min_WS_line = fm.min_WS_eff()
+ if 0:
+ plt.figure(figsize=(14, 3))
+ fm.plot_wake_map()
+ min_WS_line[::10].plot(ls='-', marker='.')
+ print(np.round(min_WS_line.values[::10][1:]).tolist())
+ plt.legend()
+ plt.show()
+
+ npt.assert_almost_equal(min_WS_line.values[::10][1:], dy, 0)
+
+
+@pytest.mark.parametrize('deflectionModel,dy', [
+ (JimenezWakeDeflection,
+ [2.0, 12.0, 20.0, 26.0, 2.0, -5.0, -11.0, -16.0, -0.0, 8.0, 15.0, 20.0]),
+ ((lambda: FugaDeflection(tfp + 'fuga/2MW/Z0=0.00001000Zi=00400Zeta0=0.00E+00/')),
+ [1.0, 6.0, 12.0, 18.0, 2.0, -0.0, -4.0, -7.0, -1.0, 2.0, 4.0, 7.0]),
+ ((lambda: FugaDeflection(tfp + 'fuga/2MW/Z0=0.03000000Zi=00401Zeta0=0.00E+00/')),
+ [1.0, 6.0, 11.0, 15.0, 2.0, -0.0, -3.0, -5.0, -1.0, 2.0, 4.0, 6.0]),
+])
+def test_deflection_model_All2AllIterative(deflectionModel, dy):
+ # center line deflection 10d downstream
+ site = IEA37Site(16)
+ x, y = [0, 400, 800], [0, 0, 0]
+ windTurbines = V80()
+ D = windTurbines.diameter()
+ wfm = All2AllIterative(site, windTurbines, IEA37SimpleBastankhahGaussianDeficit(),
+ deflectionModel=deflectionModel())
+
+ yaw = [-30, 30, -30]
+
+ sim_res = wfm(x, y, yaw=yaw, wd=270, ws=10)
+ fm = sim_res.flow_map(XYGrid(x=np.arange(-D, 15 * D + 10, 10)))
+ min_WS_line = fm.min_WS_eff()
+ if 0:
+ plt.figure(figsize=(14, 3))
+ fm.plot_wake_map()
+ min_WS_line[::10].plot(ls='-', marker='.')
+ print(np.round(min_WS_line.values[::10][1:]).tolist())
+ plt.legend()
+ plt.show()
+
+ npt.assert_almost_equal(min_WS_line.values[::10][1:], dy, 0)
+
+
@pytest.mark.parametrize('deflectionModel', [m for m in get_models(DeflectionModel) if m is not None])
def test_plot_deflection_grid(deflectionModel):
site = IEA37Site(16)
diff --git a/py_wake/tests/test_wind_farm_models/test_enginering_wind_farm_model.py b/py_wake/tests/test_wind_farm_models/test_enginering_wind_farm_model.py
index 4105ef5e5..c1326f15e 100644
--- a/py_wake/tests/test_wind_farm_models/test_enginering_wind_farm_model.py
+++ b/py_wake/tests/test_wind_farm_models/test_enginering_wind_farm_model.py
@@ -27,6 +27,8 @@ from py_wake.wind_turbines import WindTurbines
from py_wake.wind_turbines.wind_turbines_deprecated import DeprecatedOneTypeWindTurbines
import pandas as pd
import os
+from py_wake.deficit_models.deficit_model import BlockageDeficitModel
+from py_wake.rotor_avg_models.rotor_avg_model import CGIRotorAvg
WindFarmModel.verbose = False
@@ -165,16 +167,31 @@ def test_aep_mixed_type():
2 * wfm([0], [0], wd=270).aep(with_wake_loss=False).sum())
-def test_All2AllIterativeDeflection():
+@pytest.mark.parametrize('deflection_model,count',
+ [(None, 1),
+ (JimenezWakeDeflection(), 4)])
+def test_All2AllIterativeDeflection(deflection_model, count):
+
+ class FugaDeficitCount(FugaDeficit):
+ counter = 0
+
+ def _calc_layout_terms(self, dw_ijlk, hcw_ijlk, h_il, dh_ijlk, D_src_il, **_):
+ self.counter += 1
+ return FugaDeficit._calc_layout_terms(self, dw_ijlk, hcw_ijlk, h_il, dh_ijlk, D_src_il, **_)
+
site = IEA37Site(16)
windTurbines = IEA37_WindTurbines()
+ deficit_model = FugaDeficitCount()
wf_model = All2AllIterative(site, windTurbines,
- wake_deficitModel=NOJDeficit(),
- superpositionModel=SquaredSum(),
- deflectionModel=JimenezWakeDeflection())
- sim_res = wf_model([0, 500], [0, 0], wd=270, ws=10, yaw=30)
+ wake_deficitModel=deficit_model,
+ superpositionModel=LinearSum(),
+ blockage_deficitModel=SelfSimilarityDeficit(),
+ rotorAvgModel=CGIRotorAvg(4),
+ deflectionModel=deflection_model, convergence_tolerance=None)
+ sim_res = wf_model([0, 500, 1000, 1500], [0, 0, 0, 0], wd=270, ws=10, yaw=[30, -30, 30, -30])
+ assert wf_model.wake_deficitModel.counter == count
if 0:
- sim_res.flow_map(XYGrid(x=np.linspace(-200, 1000, 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/wind_farm_models/engineering_models.py b/py_wake/wind_farm_models/engineering_models.py
index 84b37dd1f..6a6be023a 100644
--- a/py_wake/wind_farm_models/engineering_models.py
+++ b/py_wake/wind_farm_models/engineering_models.py
@@ -644,6 +644,8 @@ class All2AllIterative(EngineeringWindFarmModel):
'tilt_ilk': tilt_ilk,
'D_src_il': D_src_il,
'D_dst_ijl': D_src_il[na],
+ 'dw_ijlk': dw_iil[..., na],
+ 'hcw_ijlk': hcw_iil[..., na],
'cw_ijlk': np.sqrt(hcw_iil**2 + dh_iil**2)[..., na],
'dh_ijlk': dh_iil[..., na],
'h_il': h_i[:, na]
@@ -652,7 +654,7 @@ class All2AllIterative(EngineeringWindFarmModel):
# Iterate until convergence
for j in tqdm(range(I), disable=I <= 1 or not self.verbose, desc="Calculate flow interaction", unit="wt"):
- ct_ilk = self.windTurbines.ct(WS_eff_ilk, **kwargs)
+ ct_ilk = self.windTurbines.ct(np.maximum(WS_eff_ilk, 0), **kwargs)
args['ct_ilk'] = ct_ilk
args['WS_eff_ilk'] = WS_eff_ilk
if self.deflectionModel:
@@ -660,9 +662,10 @@ class All2AllIterative(EngineeringWindFarmModel):
dw_ijl=dw_iil, hcw_ijl=hcw_iil, dh_ijl=dh_iil, **args)
args.update({'dw_ijlk': dw_ijlk, 'hcw_ijlk': hcw_ijlk, 'dh_ijlk': dh_ijlk,
'cw_ijlk': np.hypot(dh_iil[..., na], hcw_ijlk)})
- else:
- args.update({'dw_ijlk': dw_iil[..., na], 'hcw_ijlk': hcw_iil[..., na], 'dh_ijlk': dh_iil[..., na]})
+ self._reset_deficit()
+ elif j == 0:
self._init_deficit(**args)
+
if self.turbulenceModel:
args['TI_eff_ilk'] = TI_eff_ilk
if 'wake_radius_ijlk' in self.turbulenceModel.args4addturb:
--
GitLab