From fc2898e1506b4df523b5880f0b2c913825982ee9 Mon Sep 17 00:00:00 2001
From: mmpe <mmpe@dtu.dk>
Date: Mon, 13 Feb 2023 14:23:23 +0100
Subject: [PATCH] rename PolarGridRotorAvg to PolarRotorAvg and add
PolarGridRotorAvg
---
.../notebooks/EngineeringWindFarmModels.ipynb | 2 +-
py_wake/rotor_avg_models/__init__.py | 3 +-
py_wake/rotor_avg_models/rotor_avg_model.py | 18 +++++++-
.../test_rotor_avg_models.py | 44 ++++++++++++++++++-
.../wind_farm_models/engineering_models.py | 8 ++--
py_wake/wind_farm_models/wind_farm_model.py | 16 ++++++-
6 files changed, 80 insertions(+), 11 deletions(-)
diff --git a/docs/notebooks/EngineeringWindFarmModels.ipynb b/docs/notebooks/EngineeringWindFarmModels.ipynb
index 36d236db9..e563e8b04 100644
--- a/docs/notebooks/EngineeringWindFarmModels.ipynb
+++ b/docs/notebooks/EngineeringWindFarmModels.ipynb
@@ -1651,7 +1651,7 @@
"source": [
"theta = np.linspace(-np.pi,np.pi,6, endpoint=False)\n",
"r = 2/3\n",
- "plot_rotor_avg_model(PolarGridRotorAvg(nodes_r =r, nodes_theta=theta, nodes_weight=None), 'PolarGrid_6')"
+ "plot_rotor_avg_model(PolarGridRotorAvg(r=r, theta=theta, r_weight=None, theta_weight=None), 'PolarGrid_6')"
]
},
{
diff --git a/py_wake/rotor_avg_models/__init__.py b/py_wake/rotor_avg_models/__init__.py
index e9bb91589..97077a75e 100644
--- a/py_wake/rotor_avg_models/__init__.py
+++ b/py_wake/rotor_avg_models/__init__.py
@@ -1,4 +1,5 @@
-from .rotor_avg_model import RotorAvgModel, RotorCenter, EqGridRotorAvg, GQGridRotorAvg, GridRotorAvg, PolarGridRotorAvg, \
+from .rotor_avg_model import RotorAvgModel, RotorCenter, EqGridRotorAvg, GQGridRotorAvg, GridRotorAvg, \
+ PolarGridRotorAvg, PolarRotorAvg, \
CGIRotorAvg, WSPowerRotorAvg, polar_gauss_quadrature, gauss_quadrature
from .area_overlap_model import AreaOverlapAvgModel
from .gaussian_overlap_model import GaussianOverlapAvgModel
diff --git a/py_wake/rotor_avg_models/rotor_avg_model.py b/py_wake/rotor_avg_models/rotor_avg_model.py
index 2e7aed4e7..00354237a 100644
--- a/py_wake/rotor_avg_models/rotor_avg_model.py
+++ b/py_wake/rotor_avg_models/rotor_avg_model.py
@@ -36,6 +36,8 @@ class NodeRotorAvgModel(RotorAvgModel):
"""
def __call__(self, func, D_dst_ijl, **kwargs):
+ if D_dst_ijl.shape == (1, 1, 1) and D_dst_ijl[0, 0, 0] == 0:
+ return func(**kwargs)
# add extra dimension, p, with 40 points distributed over the destination rotors
kwargs = self._update_kwargs(D_dst_ijl=D_dst_ijl, **kwargs)
@@ -97,13 +99,27 @@ class GQGridRotorAvg(GridRotorAvg):
GridRotorAvg.__init__(self, nodes_x=x[m], nodes_y=y[m], nodes_weight=w)
-class PolarGridRotorAvg(GridRotorAvg):
+class PolarRotorAvg(GridRotorAvg):
def __init__(self, nodes_r=2 / 3, nodes_theta=np.linspace(-np.pi, np.pi, 6, endpoint=False), nodes_weight=None):
self.nodes_x = nodes_r * np.cos(-nodes_theta - np.pi / 2)
self.nodes_y = nodes_r * np.sin(-nodes_theta - np.pi / 2)
self.nodes_weight = nodes_weight
+class PolarGridRotorAvg(PolarRotorAvg):
+ def __init__(self, r=[1 / 3, 2 / 3], theta=np.linspace(-np.pi, np.pi, 6, endpoint=False),
+ r_weight=[.5**2, 1 - .5**2], theta_weight=1 / 6):
+ assert (r_weight is None) == (theta_weight is None)
+ nodes_r, nodes_theta = np.meshgrid(r, theta)
+ nodes_weight = None
+ if r_weight is not None:
+ r_weight, theta_weight = np.zeros(len(r)) + r_weight, np.zeros(len(theta)) + theta_weight
+ nodes_weight = np.prod(np.meshgrid(r_weight, theta_weight), 0).flatten()
+
+ PolarRotorAvg.__init__(self, nodes_r=nodes_r.flatten(), nodes_theta=nodes_theta.flatten(),
+ nodes_weight=nodes_weight)
+
+
class CGIRotorAvg(GridRotorAvg):
"""Circular Gauss Integration"""
diff --git a/py_wake/tests/test_rotor_avg_models/test_rotor_avg_models.py b/py_wake/tests/test_rotor_avg_models/test_rotor_avg_models.py
index 6410dfa38..762ee3872 100644
--- a/py_wake/tests/test_rotor_avg_models/test_rotor_avg_models.py
+++ b/py_wake/tests/test_rotor_avg_models/test_rotor_avg_models.py
@@ -13,7 +13,7 @@ from py_wake.flow_map import HorizontalGrid, XYGrid
from py_wake.rotor_avg_models import gauss_quadrature, PolarGridRotorAvg, \
polar_gauss_quadrature, EqGridRotorAvg, GQGridRotorAvg, CGIRotorAvg, GridRotorAvg, WSPowerRotorAvg
from py_wake.rotor_avg_models.gaussian_overlap_model import GaussianOverlapAvgModel
-from py_wake.rotor_avg_models.rotor_avg_model import RotorAvgModel, RotorCenter
+from py_wake.rotor_avg_models.rotor_avg_model import RotorAvgModel, RotorCenter, PolarRotorAvg
from py_wake.site._site import UniformSite
from py_wake.superposition_models import SquaredSum, LinearSum, WeightedSum
from py_wake.tests import npt
@@ -21,6 +21,8 @@ from py_wake.turbulence_models.stf import STF2017TurbulenceModel
from py_wake.utils.model_utils import get_models
from py_wake.wind_farm_models.engineering_models import All2AllIterative, PropagateDownwind, EngineeringWindFarmModel
from py_wake.turbulence_models.turbulence_model import TurbulenceModel
+from py_wake.rotor_avg_models.area_overlap_model import AreaOverlapAvgModel
+from py_wake.deficit_models.noj import NOJ
EngineeringWindFarmModel.verbose = False
@@ -277,7 +279,7 @@ def test_RotorGaussQuadratureGridAvgModel():
def test_polar_gauss_quadrature():
- m = PolarGridRotorAvg(*polar_gauss_quadrature(4, 3))
+ m = PolarRotorAvg(*polar_gauss_quadrature(4, 3))
if 0:
for v in [m.nodes_x, m.nodes_y, m.nodes_weight]:
@@ -296,6 +298,26 @@ def test_polar_gauss_quadrature():
0.08, 0.05, 0.09, 0.09, 0.05], 2)
+def test_polar_grid():
+ m = PolarGridRotorAvg(r_weight=[.5**2, 1 - .5**2], theta_weight=1 / 6)
+
+ if 0:
+ for v in [m.nodes_x, m.nodes_y, m.nodes_weight]:
+ print(np.round(v, 2).tolist())
+ plt.scatter(m.nodes_x, m.nodes_y, c=m.nodes_weight)
+ plt.axis('equal')
+ plt.gca().add_artist(plt.Circle((0, 0), 1, fill=False))
+ plt.ylim([-1, 1])
+ plt.show()
+
+ npt.assert_array_almost_equal(m.nodes_x,
+ [0.0, 0.0, 0.29, 0.58, 0.29, 0.58, 0.0, 0.0, -0.29, -0.58, -0.29, -0.58], 2)
+ npt.assert_array_almost_equal(m.nodes_y,
+ [0.33, 0.67, 0.17, 0.33, -0.17, -0.33, -0.33, -0.67, -0.17, -0.33, 0.17, 0.33], 2)
+ npt.assert_array_almost_equal(m.nodes_weight,
+ [0.04, 0.12, 0.04, 0.12, 0.04, 0.12, 0.04, 0.12, 0.04, 0.12, 0.04, 0.12], 2)
+
+
@pytest.mark.parametrize('n,x,y,w', [
(4, [-0.5, -0.5, 0.5, 0.5], [-0.5, 0.5, -0.5, 0.5], [0.25, 0.25, 0.25, 0.25]),
(7, [0.0, -0.82, 0.82, -0.41, -0.41, 0.41, 0.41],
@@ -411,3 +433,21 @@ def test_WSPowerRotorAvgModel():
rotorAvgModel = WSPowerRotorAvg(GridRotorAvg(nodes_x=[-1, 0, 1], nodes_y=[0, 0, 0]), alpha=2)
wfm = BastankhahGaussian(UniformSite(), V80(), rotorAvgModel=rotorAvgModel)
npt.assert_almost_equal(wfm(x, y, wd=270).WS_eff.sel(wt=1).squeeze(), ws_eff_ref)
+
+
+def test_flow_map():
+ wfm = NOJ(UniformSite(), V80(), rotorAvgModel=CGIRotorAvg(7))
+ y = np.linspace(-110, -50, 10)
+ fm = wfm([0], [0], wd=270, ws=12).flow_map(XYGrid(x=400, y=y))
+ fm_avg = wfm([0], [0], wd=270, ws=12).flow_map(XYGrid(x=400, y=y), D_dst=40)
+ if 0:
+ print(list(np.round(fm.WS_eff.squeeze().values, 2)))
+ print(list(np.round(fm_avg.WS_eff.squeeze().values, 2)))
+ plt.plot(y, fm.WS_eff.squeeze())
+ plt.plot(y, fm_avg.WS_eff.squeeze())
+ plt.show()
+
+ npt.assert_array_almost_equal(fm.WS_eff.squeeze(),
+ [12.0, 12.0, 12.0, 12.0, 12.0, 10.62, 10.62, 10.62, 10.62, 10.62], 2)
+ npt.assert_array_almost_equal(fm_avg.WS_eff.squeeze(),
+ [12.0, 12.0, 12.0, 11.83, 11.48, 11.14, 10.79, 10.62, 10.62, 10.62], 2)
diff --git a/py_wake/wind_farm_models/engineering_models.py b/py_wake/wind_farm_models/engineering_models.py
index e1cbac619..59dce52eb 100644
--- a/py_wake/wind_farm_models/engineering_models.py
+++ b/py_wake/wind_farm_models/engineering_models.py
@@ -249,7 +249,7 @@ class EngineeringWindFarmModel(WindFarmModel):
def _calc_wt_interaction(self, **kwargs):
"""calculate WT interaction"""
- def get_map_args(self, x_j, y_j, h_j, sim_res_data):
+ def get_map_args(self, x_j, y_j, h_j, sim_res_data, D_dst=0):
wt_d_i = self.windTurbines.diameter(sim_res_data.type)
wd, ws = [np.atleast_1d(sim_res_data[k].values) for k in ['wd', 'ws']]
time = sim_res_data.get('time', False)
@@ -270,7 +270,7 @@ class EngineeringWindFarmModel(WindFarmModel):
for k in sim_res_data if k not in ['wd_bin_size', 'ws_l', 'ws_u']}
map_arg_funcs.update({
'D_src_il': lambda l: wt_d_i[:, na],
- 'D_dst_ijl': lambda l: np.zeros((I, J, 1)),
+ 'D_dst_ijl': lambda l: np.zeros((1, 1, 1)) + D_dst,
'IJLK': lambda l=slice(None), I=I, J=J, L=L, K=K: (I, J, len(np.arange(L)[l]), K)})
return map_arg_funcs, lw_j, wd, WD_il
@@ -341,9 +341,9 @@ class EngineeringWindFarmModel(WindFarmModel):
aep_j = (power_jlk * lw_j.P_ilk).sum((1, 2))
return aep_j * 365 * 24 * 1e-9
- def _flow_map(self, x_j, y_j, h_j, sim_res_data):
+ def _flow_map(self, x_j, y_j, h_j, sim_res_data, D_dst=0):
"""call this function via SimulationResult.flow_map"""
- arg_funcs, lw_j, wd, WD_il = self.get_map_args(x_j, y_j, h_j, sim_res_data)
+ arg_funcs, lw_j, wd, WD_il = self.get_map_args(x_j, y_j, h_j, sim_res_data, D_dst=D_dst)
I, J, L, K = arg_funcs['IJLK']()
if I == 0:
return (lw_j, np.broadcast_to(lw_j.WS_ilk, (len(x_j), L, K)).astype(float),
diff --git a/py_wake/wind_farm_models/wind_farm_model.py b/py_wake/wind_farm_models/wind_farm_model.py
index 44bab78b8..d666f95f1 100644
--- a/py_wake/wind_farm_models/wind_farm_model.py
+++ b/py_wake/wind_farm_models/wind_farm_model.py
@@ -727,7 +727,7 @@ class SimulationResult(xr.Dataset):
else: # pragma: no cover
raise NotImplementedError()
- def flow_map(self, grid=None, wd=None, ws=None, time=None):
+ def flow_map(self, grid=None, wd=None, ws=None, time=None, D_dst=0):
"""Return a FlowMap object with WS_eff and TI_eff of all grid points
Parameters
@@ -736,6 +736,18 @@ class SimulationResult(xr.Dataset):
Grid, e.g. HorizontalGrid or\n
tuple(X, Y, x, y, h) where X, Y is the meshgrid for visualizing data\n
and x, y, h are the flattened grid points
+ wd : int, float, array_like or None
+ Wind directions to include in the flow map (if more than one, an weighted average will be computed)
+ The simulation result must include the requested wind directions.
+ If None, an weighted average of all wind directions from the simulation results will be computed.
+ Note, computing a flow map with multiple wind directions may be slow
+ ws : int, float, array_like or None
+ Same as "wd", but for wind speed
+ ws : int, array_like or None
+ Same as "wd", but for time
+ D_dst : int, float or None
+ In combination with a rotor average model, D_dst defines the downstream rotor diameter
+ at which the deficits will be averaged
See Also
--------
@@ -750,7 +762,7 @@ class SimulationResult(xr.Dataset):
for k in self.__slots__:
setattr(sim_res, k, getattr(self, k))
- lw_j, WS_eff_jlk, TI_eff_jlk = self.windFarmModel._flow_map(x_j, y_j, h_j, sim_res)
+ lw_j, WS_eff_jlk, TI_eff_jlk = self.windFarmModel._flow_map(x_j, y_j, h_j, sim_res, D_dst=D_dst)
return FlowMap(sim_res, X, Y, lw_j, WS_eff_jlk, TI_eff_jlk, plane=plane)
def _wd_ws(self, wd, ws):
--
GitLab