diff --git a/py_wake/deficit_models/hybridinduction.py b/py_wake/deficit_models/hybridinduction.py index fcc82b10057c5d08cd9285735a45d08598905098..02f8234dc13faedd02a5aab6cddb850697aeba22 100644 --- a/py_wake/deficit_models/hybridinduction.py +++ b/py_wake/deficit_models/hybridinduction.py @@ -5,6 +5,7 @@ from py_wake.deficit_models.vortexdipole import VortexDipole from py_wake.deficit_models import DeficitModel from py_wake.deficit_models import BlockageDeficitModel from py_wake.ground_models.ground_models import NoGround +from py_wake.utils.gradients import hypot class HybridInduction(BlockageDeficitModel): @@ -50,8 +51,8 @@ class HybridInduction(BlockageDeficitModel): # apply deficits in specified regions R_il = (D_src_il / 2) # radial distance from rotor centre - r_ijlk = np.hypot(dw_ijlk, cw_ijlk) - rcut_ijlk = (R_il * self.switch_radius)[:, na, :, na] * np.ones_like(dff_ijlk) + r_ijlk = hypot(dw_ijlk, cw_ijlk) + rcut_ijlk = np.broadcast_to((R_il * self.switch_radius)[:, na, :, na], dff_ijlk.shape) # region where to apply the far-field deficit iff = (r_ijlk > rcut_ijlk) | (dw_ijlk > 0) deficit_ijlk = dnr_ijlk.copy() diff --git a/py_wake/deficit_models/rankinehalfbody.py b/py_wake/deficit_models/rankinehalfbody.py index 56b56046493fdbb1f4c5a1694fe0348cd99e2f6a..a1bc71263dadd54c6eb632b29915cff3cb336435 100644 --- a/py_wake/deficit_models/rankinehalfbody.py +++ b/py_wake/deficit_models/rankinehalfbody.py @@ -3,6 +3,7 @@ from numpy import newaxis as na from py_wake.deficit_models import DeficitModel from py_wake.deficit_models import BlockageDeficitModel from py_wake.ground_models.ground_models import NoGround +from py_wake.utils.gradients import hypot class RankineHalfBody(BlockageDeficitModel): @@ -33,7 +34,9 @@ class RankineHalfBody(BlockageDeficitModel): """ BEM axial induction approximation by Madsen (1997). """ - a0_ilk = self.a0p[2] * ct_ilk**3 + self.a0p[1] * ct_ilk**2 + self.a0p[0] * ct_ilk + # Evaluate with Horner's rule. + # a0_ilk = self.a0p[2] * ct_ilk**3 + self.a0p[1] * ct_ilk**2 + self.a0p[0] * ct_ilk + a0_ilk = ct_ilk * (self.a0p[0] + ct_ilk * (self.a0p[1] + ct_ilk * self.a0p[2])) return a0_ilk def outside_body(self, WS_ilk, a0_ilk, R_il, dw_ijlk, cw_ijlk, r_ijlk): @@ -55,7 +58,7 @@ class RankineHalfBody(BlockageDeficitModel): R_il = D_src_il / 2. m_ilk = 2. * WS_ilk * a0_ilk * np.pi * R_il[:, :, na]**2 # radial distance - r_ijlk = np.hypot(dw_ijlk, cw_ijlk) + r_ijlk = hypot(dw_ijlk, cw_ijlk) # find points lying outside RHB, the only ones to be computed # remove singularities r_ijlk[2 * r_ijlk / D_src_il[:, na, :, na] < self.limiter] = np.inf diff --git a/py_wake/deficit_models/rathmann.py b/py_wake/deficit_models/rathmann.py index 78139a5dabcedbffe9d910d8d6a1cc07823a437b..a4180174d78df12aac99a8118aa9f25fd8c3d888 100644 --- a/py_wake/deficit_models/rathmann.py +++ b/py_wake/deficit_models/rathmann.py @@ -3,6 +3,7 @@ from numpy import newaxis as na from py_wake.deficit_models import DeficitModel from py_wake.deficit_models import BlockageDeficitModel from py_wake.ground_models.ground_models import NoGround +from py_wake.utils.gradients import hypot class Rathmann(BlockageDeficitModel): @@ -44,7 +45,7 @@ class Rathmann(BlockageDeficitModel): rho_ijlk = cw_ijlk / R_ijlk xi_ijlk = dw_ijlk / R_ijlk # mirror the bahaviour in the rotor-plane - xi_ijlk[xi_ijlk > 0] = -xi_ijlk[xi_ijlk > 0] + np.negative(xi_ijlk, out=xi_ijlk, where=xi_ijlk > 0) # centerline shape function dmu_ijlk = self.dmu(xi_ijlk) # radial shape function @@ -56,7 +57,9 @@ class Rathmann(BlockageDeficitModel): """ BEM axial induction approximation by Madsen (1997). """ - a0_ilk = self.a0p[2] * ct_ilk**3 + self.a0p[1] * ct_ilk**2 + self.a0p[0] * ct_ilk + # Evaluate with Horner's rule. + # a0_ilk = self.a0p[2] * ct_ilk**3 + self.a0p[1] * ct_ilk**2 + self.a0p[0] * ct_ilk + a0_ilk = ct_ilk * (self.a0p[0] + ct_ilk * (self.a0p[1] + ct_ilk * self.a0p[2])) return a0_ilk def dmu(self, xi_ijlk): @@ -144,7 +147,7 @@ class RathmannScaled(Rathmann): Scaling function defined in [1], Eq. 11-13 forcing the output closer to the CFD results. """ - r = np.hypot(cw_ijlk, dw_ijlk) / D_src_il[:, na, :, na] + r = hypot(cw_ijlk, dw_ijlk) / D_src_il[:, na, :, na] fac = 1. + self.sd[4] * (r - 5.) fac[fac > 1.] = 1. mval = 1. - 4. * self.sd[4] diff --git a/py_wake/deficit_models/vortexcylinder.py b/py_wake/deficit_models/vortexcylinder.py index b02e3f4f525117d2d43822e28776324c573e3550..64081613792e82f18e3d25895f52001434502eb2 100644 --- a/py_wake/deficit_models/vortexcylinder.py +++ b/py_wake/deficit_models/vortexcylinder.py @@ -5,6 +5,7 @@ from py_wake.ground_models.ground_models import NoGround from py_wake.utils.elliptic import ellipticPiCarlson from py_wake.deficit_models import DeficitModel from py_wake.deficit_models import BlockageDeficitModel +from py_wake.utils.gradients import hypot class VortexCylinder(BlockageDeficitModel): @@ -38,7 +39,9 @@ class VortexCylinder(BlockageDeficitModel): """ BEM axial induction approximation by Madsen (1997). """ - a0_ilk = self.a0p[2] * ct_ilk**3 + self.a0p[1] * ct_ilk**2 + self.a0p[0] * ct_ilk + # Evaluate with Horner's rule. + # a0_ilk = self.a0p[2] * ct_ilk**3 + self.a0p[1] * ct_ilk**2 + self.a0p[0] * ct_ilk + a0_ilk = ct_ilk * (self.a0p[0] + ct_ilk * (self.a0p[1] + ct_ilk * self.a0p[2])) return a0_ilk def calc_deficit(self, WS_ilk, D_src_il, dw_ijlk, cw_ijlk, ct_ilk, **_): @@ -52,7 +55,7 @@ class VortexCylinder(BlockageDeficitModel): R_il = D_src_il / 2 # radial distance from turbine centre - r_ijlk = np.hypot(dw_ijlk, cw_ijlk) + r_ijlk = hypot(dw_ijlk, cw_ijlk) # circulation/strength of vortex cylinder gammat_ilk = WS_ilk * 2. * self.a0(ct_ilk) # initialize deficit diff --git a/py_wake/deficit_models/vortexdipole.py b/py_wake/deficit_models/vortexdipole.py index e2fac2a48ef1cdf944fe688e1bf779508b00df92..b16f687a5111df1c0742ebd67d7746674987feff 100644 --- a/py_wake/deficit_models/vortexdipole.py +++ b/py_wake/deficit_models/vortexdipole.py @@ -3,6 +3,7 @@ from numpy import newaxis as na from py_wake.deficit_models import DeficitModel from py_wake.deficit_models import BlockageDeficitModel from py_wake.ground_models.ground_models import NoGround +from py_wake.utils.gradients import hypot class VortexDipole(BlockageDeficitModel): @@ -40,7 +41,9 @@ class VortexDipole(BlockageDeficitModel): """ BEM axial induction approximation by Madsen (1997). """ - a0_ilk = self.a0p[2] * ct_ilk**3 + self.a0p[1] * ct_ilk**2 + self.a0p[0] * ct_ilk + # Evaluate with Horner's rule. + # a0_ilk = self.a0p[2] * ct_ilk**3 + self.a0p[1] * ct_ilk**2 + self.a0p[0] * ct_ilk + a0_ilk = ct_ilk * (self.a0p[0] + ct_ilk * (self.a0p[1] + ct_ilk * self.a0p[2])) return a0_ilk def calc_deficit(self, WS_ilk, D_src_il, dw_ijlk, cw_ijlk, ct_ilk, **_): @@ -49,7 +52,7 @@ class VortexDipole(BlockageDeficitModel): """ R_il = (D_src_il / 2) # radial distance - r_ijlk = np.hypot(dw_ijlk, cw_ijlk) + r_ijlk = hypot(dw_ijlk, cw_ijlk) # circulation/strength of vortex dipole Eq. (1) in [1] gammat_ilk = WS_ilk * 2. * self.a0(ct_ilk * self.sct) # Eq. (2) in [1], induced velocities away from centreline, however diff --git a/py_wake/deflection_models/gcl_hill_vortex.py b/py_wake/deflection_models/gcl_hill_vortex.py index 72c9a696c4e7ad64b1bd51b990cbc28587e0d154..b6ed3b75adc7f135f85ead231537857f260dcdcc 100644 --- a/py_wake/deflection_models/gcl_hill_vortex.py +++ b/py_wake/deflection_models/gcl_hill_vortex.py @@ -1,6 +1,7 @@ from numpy import newaxis as na import numpy as np from py_wake.deflection_models import DeflectionModel +from py_wake.utils.gradients import hypot class GCLHillDeflection(DeflectionModel): @@ -47,7 +48,7 @@ class GCLHillDeflection(DeflectionModel): **{k: v[..., na] for k, v in kwargs.items()}) theta_yaw_ilk, theta_tilt_ilk = np.deg2rad(yaw_ilk), np.deg2rad(-tilt_ilk) - theta_ilk = np.sqrt(theta_yaw_ilk**2 + theta_tilt_ilk**2) + theta_ilk = hypot(theta_yaw_ilk, theta_tilt_ilk) theta_deflection_ilk = np.arctan2(theta_tilt_ilk, theta_yaw_ilk) U_d_ijlkx = -0.4 * U_w_ijlx * np.sin(theta_ilk)[:, na, :, :, na] diff --git a/py_wake/deflection_models/jimenez.py b/py_wake/deflection_models/jimenez.py index 1c3d00deadda82cbcaac774421c60a2ba9c80f49..01cf73d63b298b8fd537da2114a7e11b6d0ac6fd 100644 --- a/py_wake/deflection_models/jimenez.py +++ b/py_wake/deflection_models/jimenez.py @@ -1,6 +1,7 @@ from numpy import newaxis as na import numpy as np from py_wake.deflection_models import DeflectionModel +from py_wake.utils.gradients import hypot class JimenezWakeDeflection(DeflectionModel): @@ -19,7 +20,7 @@ class JimenezWakeDeflection(DeflectionModel): dw_lst = (np.logspace(0, 1.1, self.N) - 1) / (10**1.1 - 1) dw_ijxl = dw_ijl[:, :, na] * dw_lst[na, na, :, na] theta_yaw_ilk, theta_tilt_ilk = np.deg2rad(yaw_ilk), np.deg2rad(-tilt_ilk) - theta_ilk = np.sqrt(theta_yaw_ilk**2 + theta_tilt_ilk**2) + theta_ilk = hypot(theta_yaw_ilk, theta_tilt_ilk) 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, :] * np.maximum(dw_ijxl, 0))**2 diff --git a/py_wake/examples/data/iea34_130rwt/_iea34_130rwt.py b/py_wake/examples/data/iea34_130rwt/_iea34_130rwt.py index 4ba89ed3d0524336578aab33897c232d1fb7d7f9..54c3c5cf9b50512d096193464983a197ba268a5d 100644 --- a/py_wake/examples/data/iea34_130rwt/_iea34_130rwt.py +++ b/py_wake/examples/data/iea34_130rwt/_iea34_130rwt.py @@ -7,6 +7,7 @@ from py_wake.wind_turbines.power_ct_functions import PowerCtSurrogate from py_wake.wind_turbines.wind_turbine_functions import FunctionSurrogates from py_wake.examples.data import example_data_path from py_wake.utils.model_utils import fix_shape +from py_wake.utils.gradients import hypot class IEA34_130_PowerCtSurrogate(PowerCtSurrogate): @@ -109,7 +110,7 @@ class IEA34_130_2WT_Surrogate(IEA34_130_Base): def get_input(self, ws, dw_ijl, hcw_ijl, TI, Alpha=0): # ['ws','ti', 'shear', 'wdir', 'dist'] - dist = np.atleast_1d((np.hypot(dw_ijl, hcw_ijl) / 130)) + dist = np.atleast_1d((hypot(dw_ijl, hcw_ijl) / 130)) wd = np.atleast_1d(np.rad2deg(np.arctan2(hcw_ijl, dw_ijl))) unwaked = (dist == 0) | (dist > self.max_dist) | (np.abs(wd) > self.max_angle) dist[unwaked] = 20 diff --git a/py_wake/tests/test_utils/test_gradients.py b/py_wake/tests/test_utils/test_gradients.py index 526f2899b3f0f52fe35599ca979c959d622293fd..de2a27bdba3a33d55ae4c6e26c035001ee5d5f22 100644 --- a/py_wake/tests/test_utils/test_gradients.py +++ b/py_wake/tests/test_utils/test_gradients.py @@ -2,7 +2,7 @@ import matplotlib.pyplot as plt import numpy as np from autograd import numpy as anp from py_wake.examples.data.iea37._iea37 import IEA37_WindTurbines -from py_wake.utils.gradients import use_autograd_in, autograd, plot_gradients, fd, cs +from py_wake.utils.gradients import use_autograd_in, autograd, plot_gradients, fd, cs, hypot from py_wake.tests import npt from py_wake.wind_turbines import WindTurbines from py_wake.wind_turbines import _wind_turbines @@ -218,3 +218,14 @@ def test_plot_gradients(): if 0: plt.show() plt.close('all') + + +def test_hypot(): + # Test real. + a = np.array([3, 9]) + b = np.array([4, 40]) + npt.assert_equal(hypot(a, b), np.array([5, 41])) + # Test complex. + a = 3 + 4j + b = 1 - 2j + npt.assert_array_almost_equal(hypot(a, b), 2.486028939392892 + 4.022479320953552j) diff --git a/py_wake/turbulence_models/gcl_turb.py b/py_wake/turbulence_models/gcl_turb.py index fa05b05b0eeb460110cfc81aa92e51e8c15273d9..c909f7d10ec982af6ba67606b62525ce294b862b 100644 --- a/py_wake/turbulence_models/gcl_turb.py +++ b/py_wake/turbulence_models/gcl_turb.py @@ -41,8 +41,7 @@ class GCLTurbulence(TurbulenceModel, AreaOverlappingFactor): Added turbulence intensity [-] """ dw_ijlk_gt0 = np.maximum(dw_ijlk, 1e-10) - r = 0.29 * np.sqrt(1 - np.sqrt(1 - ct_ilk))[:, na] / \ - ((dw_ijlk_gt0 / D_src_il[:, na, :, na])**(1 / 3)) + r = 0.29 * np.sqrt(1 - np.sqrt(1 - ct_ilk))[:, na] / np.cbrt(dw_ijlk_gt0 / D_src_il[:, na, :, na]) area_overlap_ijlk = self.overlapping_area_factor(wake_radius_ijlk, dw_ijlk, cw_ijlk, D_src_il, D_dst_ijl) return area_overlap_ijlk * r * (dw_ijlk > 0) diff --git a/py_wake/turbulence_models/stf.py b/py_wake/turbulence_models/stf.py index ef31fe09a374dbd026100c6597a65112325a40da..26c047969d4d712653fafe028e278bb2b58a7133 100644 --- a/py_wake/turbulence_models/stf.py +++ b/py_wake/turbulence_models/stf.py @@ -2,6 +2,7 @@ from numpy import newaxis as na import numpy as np from py_wake.turbulence_models.turbulence_model import TurbulenceModel from py_wake.superposition_models import LinearSum +from py_wake.utils.gradients import hypot class FrandsenWeight(): @@ -29,7 +30,7 @@ class FrandsenWeight(): # I0 is added in the LinearSum TurbulenceSuperpositionModel # so we need to multiply the weight to I0 * alpha # 3.17: I0 * alpha = sqrt(Iadd^2 + I0^2) - I0 - return weights_ijlk * (np.sqrt(TI_add_ijlk**2 + TI_ilk[:, na]**2) - TI_ilk[:, na]) + return weights_ijlk * (hypot(TI_add_ijlk, TI_ilk[:, na]) - TI_ilk[:, na]) class IECWeight(): @@ -50,7 +51,7 @@ class IECWeight(): angleSpread = 21.6 / 2 # half angle r = np.tan(angleSpread * np.pi / 180.0) * dw_ijlk weights_ijlk = ((np.abs(cw_ijlk) < np.abs(r)) & (dw_ijlk > -1e-10) & - (np.sqrt(dw_ijlk**2 + cw_ijlk**2) < (self.dist_limit * D_src_il)[:, na, :, na])) + (hypot(dw_ijlk, cw_ijlk) < (self.dist_limit * D_src_il)[:, na, :, na])) return TI_add_ijlk * weights_ijlk @@ -66,7 +67,7 @@ class STF2017TurbulenceModel(TurbulenceModel): self.apply_weight = weight_function.apply_weight def max_centre_wake_turbulence(self, dw_ijlk, cw_ijlk, D_src_il, ct_ilk, **_): - dist_ijlk = np.sqrt(dw_ijlk**2 + cw_ijlk**2) / D_src_il[:, na, :, na] + dist_ijlk = hypot(dw_ijlk, cw_ijlk) / D_src_il[:, na, :, na] # In the standard (see page 103), the maximal added TI is calculated as # TI_add = 1/(1.5 + 0.8*d/sqrt(Ct)) return 1 / (self.c[0] + self.c[1] * dist_ijlk / np.sqrt(ct_ilk)[:, na]) diff --git a/py_wake/utils/area_overlapping_factor.py b/py_wake/utils/area_overlapping_factor.py index f78cf4c7449812baaad59bc8b5034178c46ac4c0..e4d977afd2359f484e4167e4bc8be10ca5325eea 100644 --- a/py_wake/utils/area_overlapping_factor.py +++ b/py_wake/utils/area_overlapping_factor.py @@ -89,7 +89,7 @@ class AreaOverlappingFactor(): # -1.0, cause problem to arccos(), resulting nan values, here fix this # issue. def arccos_lim(x): - return np.arccos(np.maximum(np.minimum(x, 1), -1)) + return np.arccos(np.clip(x, -1.0, +1.0)) alpha = arccos_lim((Rmax[mask]**2.0 + d[mask]**2 - Rmin[mask]**2) / (2.0 * Rmax[mask] * d[mask])) diff --git a/py_wake/utils/gradients.py b/py_wake/utils/gradients.py index 8ddae93b557b3aaa9dd7d388a8d00b0575d78fc2..9fbeb0a774813bbeb5d3b6dbba5985789537ce12 100644 --- a/py_wake/utils/gradients.py +++ b/py_wake/utils/gradients.py @@ -106,3 +106,25 @@ def plot_gradients(f, dfdx, x, label, step=1, ax=None): if label not in color_dict: color_dict[label] = c plt.legend() + + +def hypot(a, b): + """ + Given the “legs†of a right triangle, return its hypotenuse. + + Calls numpy.hypot(a, b) for real arguments and np.sqrt(a**2 + b**2) for complex arguments. + + Parameters + ---------- + a, b : real or complex array_like + Leg of the triangle(s). + + Returns + ------- + c : real or complex array_like + The hypotenuse of the triangle(s). + """ + if np.isrealobj(a) and np.isrealobj(b): + return np.hypot(a, b) + else: + return np.sqrt(a**2 + b**2) diff --git a/py_wake/validation/validation_lib.py b/py_wake/validation/validation_lib.py index 2fba67d78fabecfb73a284cf1d9ab962d17f60a1..fb5f6f3bd868959ce230e045489f20e0332998b6 100644 --- a/py_wake/validation/validation_lib.py +++ b/py_wake/validation/validation_lib.py @@ -11,6 +11,8 @@ from py_wake.deficit_models.gaussian import BastankhahGaussianDeficit from py_wake.superposition_models import SquaredSum from py_wake.rotor_avg_models import RotorCenter import xarray as xr +from py_wake.utils.gradients import hypot + # ----------------------------------------------------- # Default values @@ -66,8 +68,8 @@ def gauss(mu, sigma, x): def sigmaVarDist(x, y, xref, yref): - dist = np.hypot(x - xref, y - yref) - sigma = np.hypot(0.00035 * dist + 2.1, 2.5) + dist = hypot(x - xref, y - yref) + sigma = hypot(0.00035 * dist + 2.1, 2.5) return sigma diff --git a/py_wake/wind_farm_models/engineering_models.py b/py_wake/wind_farm_models/engineering_models.py index abf3a67684612afafd174a8243298714a8949e57..68d9e59a98439d555b082485d246470bc3f21ce7 100644 --- a/py_wake/wind_farm_models/engineering_models.py +++ b/py_wake/wind_farm_models/engineering_models.py @@ -12,6 +12,7 @@ from tqdm import tqdm from py_wake.wind_turbines._wind_turbines import WindTurbines from py_wake.utils.model_utils import check_model from py_wake.utils.functions import mean_deg +from py_wake.utils.gradients import hypot class EngineeringWindFarmModel(WindFarmModel): @@ -299,7 +300,7 @@ class EngineeringWindFarmModel(WindFarmModel): ** {k: arg_funcs[k]() for k in self.deflectionModel.args4deflection}) else: dw_ijlk, hcw_ijlk, dh_ijlk = dw_ijl[..., na], hcw_ijl[..., na], dh_ijl[..., na] - arg_funcs.update({'cw_ijlk': lambda: np.hypot(dh_ijlk, hcw_ijlk), + arg_funcs.update({'cw_ijlk': lambda: hypot(dh_ijlk, hcw_ijlk), 'dw_ijlk': lambda: dw_ijlk, 'hcw_ijlk': lambda: hcw_ijlk, 'dh_ijlk': lambda: dh_ijlk}) args = {k: arg_funcs[k]() for k in self.args4deficit if k != 'dw_ijlk'} @@ -347,7 +348,7 @@ class EngineeringWindFarmModel(WindFarmModel): l_ = [l, 0][lw_j.WS_ilk.shape[1] == 1] if isinstance(self.superpositionModel, WeightedSum): - cw_ijk = np.hypot(dh_ijl[..., na], hcw_ijlk)[:, :, 0] + cw_ijk = hypot(dh_ijl[..., na], hcw_ijlk)[:, :, 0] hcw_ijk, dh_ijk = hcw_ijlk[:, :, 0], dh_ijl[:, :, 0, na] WS_eff_jlk[:, l] = lw_j.WS_ilk[:, l_] - self.superpositionModel(lw_j.WS_ilk[:, l_], deficit_ijk, uc_ijk, sigma_sqr_ijk, cw_ijk, hcw_ijk, dh_ijk) @@ -689,7 +690,7 @@ class All2AllIterative(EngineeringWindFarmModel): dw_ijlk, hcw_ijlk, dh_ijlk = self.deflectionModel.calc_deflection( 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)}) + 'cw_ijlk': hypot(dh_iil[..., na], hcw_ijlk)}) self._reset_deficit() if self.turbulenceModel: