From 758a1c85840087af24567e608a60bc0a8891d00d Mon Sep 17 00:00:00 2001
From: mmpe <mmpe@dtu.dk>
Date: Tue, 8 Dec 2020 11:43:48 +0100
Subject: [PATCH] Allow different rotor average model for turbulence
---
py_wake/deficit_models/no_wake.py | 4 ++
.../test_deficit_models.py | 5 +-
py_wake/tests/test_rotor_avg_models.py | 6 +-
.../test_turbulence_models.py | 67 ++++++++++++++++++-
py_wake/turbulence_models/crespo.py | 4 +-
py_wake/turbulence_models/gcl.py | 4 +-
py_wake/turbulence_models/stf.py | 5 +-
py_wake/turbulence_models/turbulence_model.py | 3 +-
.../wind_farm_models/engineering_models.py | 11 +--
9 files changed, 92 insertions(+), 17 deletions(-)
diff --git a/py_wake/deficit_models/no_wake.py b/py_wake/deficit_models/no_wake.py
index 2d09f86a8..77025dc4b 100644
--- a/py_wake/deficit_models/no_wake.py
+++ b/py_wake/deficit_models/no_wake.py
@@ -1,7 +1,11 @@
from py_wake.deficit_models import DeficitModel
from numpy import newaxis as na
+import numpy as np
class NoWakeDeficit(DeficitModel):
def calc_deficit(self, WS_ilk, dw_ijlk, **_):
return (WS_ilk)[:, na] * (dw_ijlk > 0) * 0
+
+ def wake_radius(self, dw_ijlk, **_):
+ return np.zeros_like(dw_ijlk)
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 ba6a02a33..d7a1ab6a1 100644
--- a/py_wake/tests/test_deficit_models/test_deficit_models.py
+++ b/py_wake/tests/test_deficit_models/test_deficit_models.py
@@ -28,6 +28,7 @@ from py_wake.turbulence_models.gcl import GCLTurbulence
from py_wake.wind_farm_models.engineering_models import PropagateDownwind, All2AllIterative
from py_wake.turbulence_models.stf import STF2017TurbulenceModel
from py_wake.examples.data.hornsrev1 import Hornsrev1Site
+from py_wake.deficit_models.selfsimilarity import SelfSimilarityDeficit
class GCLLocalDeficit(GCLDeficit):
@@ -229,9 +230,9 @@ def test_wake_radius_not_implemented():
site = IEA37Site(16)
x, y = site.initial_position.T
windTurbines = IEA37_WindTurbines()
- wfm = PropagateDownwind(site, windTurbines, wake_deficitModel=NoWakeDeficit(),
+ wfm = PropagateDownwind(site, windTurbines, wake_deficitModel=SelfSimilarityDeficit(),
turbulenceModel=GCLTurbulence())
- with pytest.raises(NotImplementedError, match="wake_radius not implemented for NoWakeDeficit"):
+ with pytest.raises(NotImplementedError, match="wake_radius not implemented for SelfSimilarityDeficit"):
wfm(x, y)
diff --git a/py_wake/tests/test_rotor_avg_models.py b/py_wake/tests/test_rotor_avg_models.py
index 272820464..4ac04a986 100644
--- a/py_wake/tests/test_rotor_avg_models.py
+++ b/py_wake/tests/test_rotor_avg_models.py
@@ -29,7 +29,11 @@ def test_RotorGridAvg_deficit():
('RotorGrid3', EqGridRotorAvg(3), 7.633415167369133),
('RotorGrid4', EqGridRotorAvg(4), 7.710215921858325),
('RotorGrid100', EqGridRotorAvg(100), 7.820762402628349),
- ('RotorGQGrid_4,3', GQGridRotorAvg(4, 3), 7.826105012683896)]:
+ ('RotorGQGrid_4,3', GQGridRotorAvg(4, 3), 7.826105012683896),
+ ('RotorCGI4', CGIRotorAvg(4), 7.848406907726826),
+ ('RotorCGI4', CGIRotorAvg(7), 7.819900693605533),
+ ('RotorCGI4', CGIRotorAvg(9), 7.82149363932618),
+ ('RotorCGI4', CGIRotorAvg(21), 7.821558905416136)]:
# test with PropagateDownwind
wfm = IEA37SimpleBastankhahGaussian(site,
diff --git a/py_wake/tests/test_turbulence_models/test_turbulence_models.py b/py_wake/tests/test_turbulence_models/test_turbulence_models.py
index 7a03dc21d..4fc6908ce 100644
--- a/py_wake/tests/test_turbulence_models/test_turbulence_models.py
+++ b/py_wake/tests/test_turbulence_models/test_turbulence_models.py
@@ -18,12 +18,15 @@ from py_wake.wind_farm_models.engineering_models import PropagateDownwind, All2A
from py_wake.turbulence_models.gcl import GCLTurbulence
import matplotlib.pyplot as plt
from py_wake.turbulence_models.crespo import CrespoHernandez
-from py_wake.deficit_models.gaussian import BastankhahGaussian, IEA37SimpleBastankhahGaussianDeficit
+from py_wake.deficit_models.gaussian import BastankhahGaussian, IEA37SimpleBastankhahGaussianDeficit,\
+ IEA37SimpleBastankhahGaussian
import xarray as xr
import os
import pkgutil
import inspect
from py_wake.examples.data.hornsrev1 import Hornsrev1Site
+from py_wake.rotor_avg_models.rotor_avg_model import RotorCenter, EqGridRotorAvg, GQGridRotorAvg, CGIRotorAvg
+from py_wake.deficit_models.gcl import GCLDeficit
def get_all_turbulence_models():
@@ -35,7 +38,7 @@ def get_all_turbulence_models():
v = _module.__dict__[n]
if inspect.isclass(v) and issubclass(v, TurbulenceModel) and \
v not in [TurbulenceModel]:
- all_turbulence_modules.append(v())
+ all_turbulence_modules.append(v)
return all_turbulence_modules
@@ -146,6 +149,64 @@ def test_own_turbulence_is_zero(turbulenceModel):
site = Hornsrev1Site()
windTurbines = IEA37_WindTurbines()
wf_model = All2AllIterative(site, windTurbines, wake_deficitModel=IEA37SimpleBastankhahGaussianDeficit(),
- turbulenceModel=turbulenceModel)
+ turbulenceModel=turbulenceModel())
sim_res = wf_model([0], [0])
npt.assert_array_equal(sim_res.TI_eff, sim_res.TI.broadcast_like(sim_res.TI_eff))
+
+
+def test_RotorAvg_deficit():
+ site = IEA37Site(16)
+ 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 [
+ ('None', None, 0.22292190804089568),
+ ('RotorCenter', RotorCenter(), 0.22292190804089568),
+ ('RotorGrid100', EqGridRotorAvg(100), 0.1989725533174574),
+ ('RotorGQGrid_4,3', GQGridRotorAvg(4, 3), 0.19874837617113356),
+ ('RotorCGI4', CGIRotorAvg(4), 0.19822024411411204),
+ ('RotorCGI4', CGIRotorAvg(21), 0.1989414764606653)]:
+
+ # test with PropagateDownwind
+ wfm = IEA37SimpleBastankhahGaussian(site,
+ windTurbines,
+ turbulenceModel=STF2017TurbulenceModel(rotorAvgModel=rotorAvgModel))
+ sim_res = wfm([0, 500], [0, 0], wd=270, ws=10)
+ npt.assert_almost_equal(sim_res.TI_eff_ilk[1, 0, 0], ref1, err_msg=name)
+
+ # test with All2AllIterative
+ wfm = All2AllIterative(site, windTurbines,
+ IEA37SimpleBastankhahGaussianDeficit(),
+ turbulenceModel=STF2017TurbulenceModel(rotorAvgModel=rotorAvgModel),
+ superpositionModel=SquaredSum())
+ 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.WS_eff_ilk[1, 0, 0]] * 2, label=name)
+ if 0:
+ plt.legend()
+ plt.show()
+ plt.close()
+
+
+@pytest.mark.parametrize('WFM', [All2AllIterative, PropagateDownwind])
+@pytest.mark.parametrize('turbulenceModel', get_all_turbulence_models())
+def test_with_all_turbulence_models(WFM, turbulenceModel):
+ site = IEA37Site(16)
+ windTurbines = IEA37_WindTurbines()
+
+ wfm = WFM(site, windTurbines, wake_deficitModel=NoWakeDeficit(),
+ rotorAvgModel=CGIRotorAvg(4),
+ superpositionModel=LinearSum(),
+ turbulenceModel=turbulenceModel())
+
+ wfm2 = WFM(site, windTurbines, wake_deficitModel=NoWakeDeficit(),
+ superpositionModel=LinearSum(),
+ turbulenceModel=turbulenceModel(rotorAvgModel=CGIRotorAvg(4)))
+ kwargs = {'x': [0, 0, 500, 500], 'y': [0, 500, 0, 500], 'wd': [0], 'ws': [8]}
+ npt.assert_array_equal(wfm(**kwargs).TI_eff, wfm2(**kwargs).TI_eff, turbulenceModel.__name__)
diff --git a/py_wake/turbulence_models/crespo.py b/py_wake/turbulence_models/crespo.py
index 5ffc01519..02e50b3b5 100644
--- a/py_wake/turbulence_models/crespo.py
+++ b/py_wake/turbulence_models/crespo.py
@@ -14,8 +14,8 @@ class CrespoHernandez(TurbulenceModel, AreaOverlappingFactor):
"""
args4addturb = ['dw_ijlk', 'cw_ijlk', 'D_src_il', 'ct_ilk', 'TI_ilk', 'D_dst_ijl', 'wake_radius_ijlk']
- def __init__(self, addedTurbulenceSuperpositionModel=SqrMaxSum()):
- TurbulenceModel.__init__(self, addedTurbulenceSuperpositionModel)
+ def __init__(self, addedTurbulenceSuperpositionModel=SqrMaxSum(), **kwargs):
+ TurbulenceModel.__init__(self, addedTurbulenceSuperpositionModel, **kwargs)
def calc_added_turbulence(self, dw_ijlk, cw_ijlk, D_src_il, ct_ilk, TI_ilk, D_dst_ijl, wake_radius_ijlk, **_):
""" Calculate the added turbulence intensity at locations specified by
diff --git a/py_wake/turbulence_models/gcl.py b/py_wake/turbulence_models/gcl.py
index 56c8724fe..b7fc558bb 100644
--- a/py_wake/turbulence_models/gcl.py
+++ b/py_wake/turbulence_models/gcl.py
@@ -15,8 +15,8 @@ class GCLTurbulence(TurbulenceModel, AreaOverlappingFactor):
"""
args4addturb = ['D_src_il', 'dw_ijlk', 'ct_ilk', 'D_dst_ijl', 'cw_ijlk', 'wake_radius_ijlk']
- def __init__(self, addedTurbulenceSuperpositionModel=SqrMaxSum()):
- TurbulenceModel.__init__(self, addedTurbulenceSuperpositionModel)
+ def __init__(self, addedTurbulenceSuperpositionModel=SqrMaxSum(), **kwargs):
+ TurbulenceModel.__init__(self, addedTurbulenceSuperpositionModel, **kwargs)
def calc_added_turbulence(self, dw_ijlk, D_src_il, ct_ilk, wake_radius_ijlk,
D_dst_ijl, cw_ijlk, **_):
diff --git a/py_wake/turbulence_models/stf.py b/py_wake/turbulence_models/stf.py
index d01f9fc49..23f5329d4 100644
--- a/py_wake/turbulence_models/stf.py
+++ b/py_wake/turbulence_models/stf.py
@@ -8,8 +8,9 @@ class STF2017TurbulenceModel(TurbulenceModel):
args4addturb = ['dw_ijlk', 'cw_ijlk', 'D_src_il', 'ct_ilk', 'TI_ilk']
- def __init__(self, addedTurbulenceSuperpositionModel=LinearSum()):
- TurbulenceModel.__init__(self, addedTurbulenceSuperpositionModel)
+ def __init__(self, addedTurbulenceSuperpositionModel=LinearSum(),
+ **kwargs):
+ TurbulenceModel.__init__(self, addedTurbulenceSuperpositionModel, **kwargs)
def weight(self, dw_ijlk, cw_ijlk, D_src_il):
# The weight is given by the exponential term in Eq 3.16 and accounts
diff --git a/py_wake/turbulence_models/turbulence_model.py b/py_wake/turbulence_models/turbulence_model.py
index d6a035b21..7c39f8f4f 100644
--- a/py_wake/turbulence_models/turbulence_model.py
+++ b/py_wake/turbulence_models/turbulence_model.py
@@ -4,9 +4,10 @@ from abc import abstractmethod
class TurbulenceModel():
- def __init__(self, addedTurbulenceSuperpositionModel):
+ def __init__(self, addedTurbulenceSuperpositionModel, rotorAvgModel=None):
assert isinstance(addedTurbulenceSuperpositionModel, AddedTurbulenceSuperpositionModel)
self.addedTurbulenceSuperpositionModel = addedTurbulenceSuperpositionModel
+ self.rotorAvgModel = rotorAvgModel
@abstractmethod
def calc_added_turbulence(self):
diff --git a/py_wake/wind_farm_models/engineering_models.py b/py_wake/wind_farm_models/engineering_models.py
index 4f5baffed..082c362d1 100644
--- a/py_wake/wind_farm_models/engineering_models.py
+++ b/py_wake/wind_farm_models/engineering_models.py
@@ -86,7 +86,10 @@ 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)
+ if self.turbulenceModel.rotorAvgModel is None:
+ self.turbulenceModel.rotorAvgModel = rotorAvgModel
+ self.args4addturb = set(self.turbulenceModel.args4addturb) | set(
+ self.turbulenceModel.rotorAvgModel.args4rotor_avg_deficit)
self.args4all = self.args4all | set(self.turbulenceModel.args4addturb)
if self.deflectionModel:
self.args4all = self.args4all | set(self.deflectionModel.args4deflection)
@@ -494,8 +497,8 @@ class PropagateDownwind(EngineeringWindFarmModel):
if k != "dw_ijlk"}
# Calculate added turbulence
- add_turb_nk[n_dw] = self.rotorAvgModel(self.turbulenceModel.calc_added_turbulence,
- dw_ijlk=dw_ijlk, **turb_args)
+ add_turb_nk[n_dw] = self.turbulenceModel.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)
@@ -596,7 +599,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.rotorAvgModel(self.turbulenceModel.calc_added_turbulence, **args)
+ add_turb_ijlk = self.turbulenceModel.rotorAvgModel(self.turbulenceModel.calc_added_turbulence, **args)
TI_eff_ilk = self.turbulenceModel.calc_effective_TI(lw.TI_ilk, add_turb_ijlk)
# Check if converged
--
GitLab