from py_wake import np
from abc import ABC, abstractmethod
from numpy import newaxis as na
from py_wake.utils.gradients import cabs
class SuperpositionModel(ABC):
def __init__(self):
pass
@abstractmethod
def __call__(self, value_jxxx):
"""Calculate the sum of jxxx
This method must be overridden by subclass
Parameters
----------
deficit_jxxx : array_like
deficit caused by source turbines(j) on xxx (xxx optionally includes
destination turbine/site, wind directions, wind speeds
Returns
-------
sum_xxx : array_like
sum for xxx (see above)
"""
class AddedTurbulenceSuperpositionModel():
def calc_effective_TI(self, TI_xxx, add_turb_jxxx):
"""Calculate effective turbulence intensity
Parameters
----------
TI_xxx : array_like
Local turbulence intensity. xxx optionally includes destination turbine/site, wind directions, wind speeds
add_turb_jxxx : array_like
added turbulence caused by source turbines(j) on xxx (see above)
Returns
-------
TI_eff_xxx : array_like
Effective turbulence intensity xxx (see TI_xxx)
"""
return TI_xxx + self(add_turb_jxxx)
class SquaredSum(SuperpositionModel, AddedTurbulenceSuperpositionModel):
def __call__(self, value_jxxx):
return np.sqrt(np.sum(value_jxxx**2, 0))
class LinearSum(SuperpositionModel, AddedTurbulenceSuperpositionModel):
def __call__(self, value_jxxx):
return np.sum(value_jxxx, 0)
class MaxSum(SuperpositionModel, AddedTurbulenceSuperpositionModel):
def __call__(self, value_jxxx):
return np.max(value_jxxx, 0)
class SqrMaxSum(AddedTurbulenceSuperpositionModel):
def calc_effective_TI(self, TI_xxx, add_turb_jxxx):
return np.sqrt(TI_xxx**2 + np.max(add_turb_jxxx, 0)**2)
class WeightedSum(SuperpositionModel):
"""
Implemented according to the paper by:
Haohua Zong and Fernando Porté-Agel
A momentum-conserving wake superposition method for wind farm power prediction
J. Fluid Mech. (2020), vol. 889, A8; doi:10.1017/jfm.2020.77
"""
def __init__(self, delta=0.01, max_err=1e-3, max_iter=5):
# minimum deficit (as fraction of free-stream) to invoke weighted summation
self.delta = delta
# convergence limit for computing global convection velocity
self.max_err = max_err
# maximum number of iterations used in computing weights
self.max_iter = max_iter
def __call__(self, WS_xxx, centerline_deficit_jxxx,
convection_velocity_jxxx,
sigma_sqr_jxxx, cw_jxxx, hcw_jxxx, dh_jxxx):
Ws = WS_xxx + np.zeros(centerline_deficit_jxxx.shape[1:])
usc = centerline_deficit_jxxx
uc = convection_velocity_jxxx
sigma_sqr = sigma_sqr_jxxx
cw = cw_jxxx
hcw = hcw_jxxx * np.ones_like(usc)
dh = dh_jxxx * np.ones_like(usc)
Us = np.zeros_like(Ws)
# Determine non-centreline deficit ratio
# Local deficit
us = usc * np.exp(-1 / (2 * sigma_sqr) * cw**2)
# Set lower deficit limit below which deficits are linearly added
us_lim = max(self.delta, 1e-30) * Ws[na]
# Get indices
Il = us >= us_lim
# Only start weighting computation if at least two deficits need to be combined
# Computations are performed where velocities exceed the specified limit (us_lim).
# The more complex indexing leads to more complex code.
if Il.any():
# Get indices where deficits need to be combined
Ilx = Il.any(axis=0)
# Total cross-wind integrated deficit
us_int = np.zeros_like(us)
us_int[Il] = usc[Il] * 2 * np.pi * sigma_sqr[Il]
# initialize combined quanatities
Uc = Ws.copy()
Uc_star = Ws.copy()
Us = np.zeros_like(Ws)
Us_int = np.zeros_like(Ws)
# Initialize
count = 0
Uc_star = 10 * Uc
if np.iscomplexobj(Ws) or np.iscomplexobj(uc):
Uc_star = Uc_star.astype(np.complex128)
tmp1, tmp2 = np.zeros_like(us), np.zeros_like(us)
tmpUS, tmpUSint = np.zeros_like(us), np.zeros_like(us)
sum1, sum2, ucn = np.zeros_like(Us), np.zeros_like(Us), np.ones_like(uc)
# sum1 precomputed part
sum1_pre = usc[Il]**2 * np.pi * sigma_sqr[Il]
n_wt = us.shape[0]
# Iterate until combined convection velocity converges
while (np.max(cabs((Uc[Ilx] - Uc_star[Ilx]) / Uc_star[Ilx])) > self.max_err) and (count < self.max_iter):
# Initialize combined convection velocity
if count == 0:
# Take maximum across all turbines to initilize global convection velocity
Uc = np.max(np.where(Il, uc, 0), 0)
else:
Uc = Uc_star.copy()
# Initialize and avoid division by zero
ucn[:] = 1.
Inz = Uc != 0
ucn[:, Inz] = uc[:, Inz] / Uc[Inz]
# Combined local deficit
# dummy matrix to keep original matrix shape
tmpUS[:] = 0
tmpUS[Il] = ucn[Il] * us[Il]
Us = np.sum(tmpUS, axis=0)
# Combined deficit integrated to infinity in cross-wind direction
tmpUSint[:] = 0
tmpUSint[Il] = ucn[Il] * us_int[Il]
Us_int = np.sum(tmpUSint, axis=0)
# First sum of momentum deficit
# sum_i^N (uc_i u_i)^2
sum1[:], tmp1[:] = .0, 0.
tmp1[Il] = ucn[Il]**2 * sum1_pre
sum1 = np.sum(tmp1, axis=0)
# Second sum which represents the cross terms
# 2 sum i>j (uc_i u_i) (uc_j u_j)
if n_wt > 1:
# Initialize
sum2[:] = .0
for j in range(n_wt - 1):
# Only cross with larger indices
k = np.arange(j + 1, n_wt, dtype=int)
# Find indices where deficits need to be combined
Ilxx = Il[j][na] & Il[k]
if Ilxx.any():
# To keep the shape, arrays are repeated and a dummy initilized
# Instead of a dummy one could use a loop, but this seemed faster
tmp2 = np.zeros(((len(k),) + sigma_sqr.shape[1:]), dtype=sigma_sqr.dtype)
s1, s2 = np.repeat(sigma_sqr[j][na], len(k), axis=0)[Ilxx], sigma_sqr[k][Ilxx]
w2w_hcw = cabs(hcw[j][na] - hcw[k])[Ilxx]
w2w_dh = cabs(dh[j][na] - dh[k])[Ilxx]
cross_sigma_jk = 2 * np.exp(-(w2w_hcw**2 + w2w_dh**2) /
(2 * (s1 + s2))) * np.pi * s1 * s2 / (s1 + s2)
tmp2[Ilxx] = 2 * np.repeat((ucn[j] * usc[j])[na], len(k), axis=0)[Ilxx] * \
(ucn[k][Ilxx] * usc[k][Ilxx]) * cross_sigma_jk
sum2 += np.sum(tmp2, axis=0)
# Avoid division by zero
Us_int[Us_int == 0] = 1
# Update combined convection velocity
Uc_star[Ilx] = Ws[Ilx] - (sum1 + sum2)[Ilx] / Us_int[Ilx]
count += 1
return Us + np.sum(np.where(~Il, us, 0), axis=0)