Skip to content
Snippets Groups Projects
Commit 4cbab6af authored by Mads M. Pedersen's avatar Mads M. Pedersen
Browse files

apply rotor average model to turbulence too

parent 6c569cd8
No related branches found
No related tags found
No related merge requests found
from py_wake.deficit_models.deficit_model import DeficitModel
import numpy as np
from abc import abstractmethod
from numpy import newaxis as na
class RotorAvgModel(DeficitModel):
class RotorAvgModel():
"""Wrap a DeficitModel.
The RotorAvgModel
- add an extra dimension (one or more points covering the downstream rotors)
......@@ -13,22 +11,19 @@ class RotorAvgModel(DeficitModel):
"""
args4rotor_avg_deficit = ['hcw_ijlk', 'dh_ijl', 'D_dst_ijl']
def calc_deficit(self, deficitModel, D_dst_ijl, **kwargs):
self.deficitModel = deficitModel
if np.all(D_dst_ijl == 0):
return self.deficitModel.calc_deficit(D_dst_ijl=D_dst_ijl, **kwargs)
else:
return self._calc_rotor_avg_deficit(D_dst_ijl=D_dst_ijl, **kwargs)
def calc_deficit_convection(self, deficitModel, D_dst_ijl, **kwargs):
self.deficitModel = deficitModel
return self.deficitModel.calc_deficit_convection(D_dst_ijl=D_dst_ijl, **kwargs)
@abstractmethod
def _calc_rotor_avg_deficit(self):
"""Similar to calc_deficit, but with an extra point dimension to calculate the
rotor average wind speed as a weighted average of a number of points instead of
only the rotor center"""
def __call__(self, func, D_dst_ijl, **kwargs):
# add extra dimension, p, with 40 points distributed over the destination rotors
kwargs = self._update_kwargs(D_dst_ijl=D_dst_ijl, **kwargs)
values_ijlkp = func(**kwargs)
# Calculate weighted sum of deficit over the destination rotors
if self.nodes_weight is None:
return np.mean(values_ijlkp, -1)
return np.sum(self.nodes_weight[na, na, na, na, :] * values_ijlkp, -1)
class RotorCenter(RotorAvgModel):
......@@ -37,8 +32,8 @@ class RotorCenter(RotorAvgModel):
nodes_y = [0]
nodes_weight = [1]
def _calc_rotor_avg_deficit(self, **kwargs):
return self.deficitModel.calc_deficit(**kwargs)
def __call__(self, func, **kwargs):
return func(**kwargs)
def _calc_layout_terms(self, deficitModel, **kwargs):
deficitModel._calc_layout_terms(**kwargs)
......@@ -59,26 +54,13 @@ class GridRotorAvg(RotorAvgModel):
hcw_ijlkp = hcw_ijlk[..., na] + R_dst_ijl[:, :, :, na, na] * self.nodes_x[na, na, na, na, :]
dh_ijlp = dh_ijl[..., na] + R_dst_ijl[:, :, :, na] * self.nodes_y[na, na, na, :]
new_kwargs = {'dh_ijl': dh_ijlp, 'hcw_ijlk': hcw_ijlkp, 'D_dst_ijl': D_dst_ijl[..., na]}
if 'cw_ijlk' in self.deficitModel.args4deficit:
cw_ijlkp = np.sqrt(hcw_ijlkp**2 + dh_ijlp[:, :, :, na]**2)
new_kwargs['cw_ijlk'] = cw_ijlkp
if 'D_dst_ijl' in self.deficitModel.args4deficit:
new_kwargs['D_dst_ijl'] = D_dst_ijl
new_kwargs['cw_ijlk'] = np.sqrt(hcw_ijlkp**2 + dh_ijlp[:, :, :, na]**2)
new_kwargs['D_dst_ijl'] = D_dst_ijl
new_kwargs.update({k: v[..., na] for k, v in kwargs.items() if k not in new_kwargs})
return new_kwargs
def _calc_rotor_avg_deficit(self, **kwargs):
# add extra dimension, p, with 40 points distributed over the destination rotors
kwargs = self._update_kwargs(**kwargs)
deficit_ijlkp = self.deficitModel.calc_deficit(**kwargs)
# Calculate weighted sum of deficit over the destination rotors
if self.nodes_weight is None:
return np.mean(deficit_ijlkp, -1)
return np.sum(self.nodes_weight[na, na, na, na, :] * deficit_ijlkp, -1)
def _calc_layout_terms(self, deficitModel, **kwargs):
self.deficitModel = deficitModel
self.deficitModel._calc_layout_terms(**self._update_kwargs(**kwargs))
......
......@@ -6,14 +6,11 @@ import matplotlib.pyplot as plt
from py_wake.flow_map import YZGrid
import numpy as np
from py_wake.tests import npt
from py_wake.wind_turbines import WindTurbines, OneTypeWindTurbines
from py_wake.wind_turbines import WindTurbines
from py_wake.superposition_models import LinearSum, SquaredSum
from py_wake.wind_farm_models.engineering_models import PropagateDownwind, All2AllIterative
import pytest
from py_wake.deficit_models.gaussian import BastankhahGaussianDeficit, IEA37SimpleBastankhahGaussianDeficit,\
ZongGaussianDeficit, NiayifarGaussianDeficit
from py_wake.deficit_models.fuga import FugaDeficit
from py_wake.deficit_models.gcl import GCLDeficit
from py_wake.deficit_models.gaussian import ZongGaussianDeficit
from py_wake.turbulence_models.stf import STF2017TurbulenceModel
from py_wake.ground_models.ground_models import MirrorSquaredSum
......@@ -42,7 +39,6 @@ def test_Mirror_NOJ():
plt.legend()
plt.show()
plt.close()
print(res)
npt.assert_array_equal(res, ref)
......
......@@ -13,7 +13,7 @@ import pytest
from py_wake.turbulence_models.stf import STF2017TurbulenceModel
def test_RotorGridAvg():
def test_RotorGridAvg_deficit():
site = IEA37Site(16)
x, y = site.initial_position.T
windTurbines = IEA37_WindTurbines()
......@@ -53,6 +53,49 @@ def test_RotorGridAvg():
plt.close()
def test_RotorGridAvg_ti():
site = IEA37Site(16)
x, y = site.initial_position.T
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 [
('RotorCenter', RotorCenter(), 0.22292190804089568),
('RotorGrid2', EqGridRotorAvg(2), 0.21135016790319944),
('RotorGrid3', EqGridRotorAvg(3), 0.20618819397725846),
('RotorGrid4', EqGridRotorAvg(4), 0.20324660406762962),
('RotorGrid100', EqGridRotorAvg(100), 0.1989725533174574),
('RotorGQGrid_4,3', GQGridRotorAvg(4, 3), 0.19874837617113356)]:
# test with PropagateDownwind
wfm = IEA37SimpleBastankhahGaussian(site,
windTurbines,
rotorAvgModel=rotorAvgModel,
turbulenceModel=STF2017TurbulenceModel())
sim_res = wfm([0, 500], [0, 0], wd=270, ws=10)
npt.assert_almost_equal(sim_res.TI_eff_ilk[1, 0, 0], ref1)
# test with All2AllIterative
wfm = All2AllIterative(site, windTurbines,
IEA37SimpleBastankhahGaussianDeficit(),
rotorAvgModel=rotorAvgModel,
superpositionModel=SquaredSum(),
turbulenceModel=STF2017TurbulenceModel())
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.TI_eff_ilk[1, 0, 0]] * 2, label=name)
if 0:
plt.legend()
plt.show()
plt.close()
def test_gauss_quadrature():
x, _, w = gauss_quadrature(4, 1)
......
......@@ -86,6 +86,7 @@ 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)
self.args4all = self.args4all | set(self.turbulenceModel.args4addturb)
if self.deflectionModel:
self.args4all = self.args4all | set(self.deflectionModel.args4deflection)
......@@ -134,9 +135,8 @@ class EngineeringWindFarmModel(WindFarmModel):
def _calc_deficit(self, dw_ijlk, **kwargs):
"""Calculate wake (and blockage) deficit"""
def f(**kwargs):
return self.rotorAvgModel.calc_deficit(self.wake_deficitModel, **kwargs)
deficit = self.groundModel(f, dw_ijlk=dw_ijlk, **kwargs)
deficit = self.groundModel(lambda **kwargs: self.rotorAvgModel(self.wake_deficitModel.calc_deficit, **kwargs),
dw_ijlk=dw_ijlk, **kwargs)
return self._add_blockage(deficit, dw_ijlk, **kwargs)
def _calc_deficit_convection(self, dw_ijlk, **kwargs):
......@@ -486,13 +486,16 @@ class PropagateDownwind(EngineeringWindFarmModel):
if self.turbulenceModel:
if 'wake_radius_ijlk' in self.turbulenceModel.args4addturb:
if 'wake_radius_ijlk' in self.args4addturb:
wake_radius_ijlk = self.wake_deficitModel.wake_radius(dw_ijlk=dw_ijlk, **args)
arg_funcs['wake_radius_ijlk'] = lambda: wake_radius_ijlk
turb_args = {k: arg_funcs[k]() for k in self.turbulenceModel.args4addturb if k != "dw_ijlk"}
turb_args = {k: arg_funcs[k]() for k in self.args4addturb
if k != "dw_ijlk"}
# Calculate added turbulence
add_turb_nk[n_dw] = self.turbulenceModel.calc_added_turbulence(dw_ijlk=dw_ijlk, **turb_args)
add_turb_nk[n_dw] = self.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)
......@@ -593,7 +596,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.turbulenceModel.calc_added_turbulence(**args)
add_turb_ijlk = self.rotorAvgModel(self.turbulenceModel.calc_added_turbulence, **args)
TI_eff_ilk = self.turbulenceModel.calc_effective_TI(lw.TI_ilk, add_turb_ijlk)
# Check if converged
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment