From 039fb94e66b3730ae3f8b69590edcfe0c23f52f0 Mon Sep 17 00:00:00 2001 From: Riccardo Riva <ricriv@dtu.dk> Date: Wed, 26 Jan 2022 06:10:25 +0000 Subject: [PATCH] Fix #67 --- py_wake/deficit_models/hybridinduction.py | 5 +++-- py_wake/deficit_models/rankinehalfbody.py | 7 ++++-- py_wake/deficit_models/rathmann.py | 9 +++++--- py_wake/deficit_models/vortexcylinder.py | 7 ++++-- py_wake/deficit_models/vortexdipole.py | 7 ++++-- py_wake/deflection_models/gcl_hill_vortex.py | 3 ++- py_wake/deflection_models/jimenez.py | 3 ++- .../data/iea34_130rwt/_iea34_130rwt.py | 3 ++- py_wake/tests/test_utils/test_gradients.py | 13 ++++++++++- py_wake/turbulence_models/gcl_turb.py | 3 +-- py_wake/turbulence_models/stf.py | 7 +++--- py_wake/utils/area_overlapping_factor.py | 2 +- py_wake/utils/gradients.py | 22 +++++++++++++++++++ py_wake/validation/validation_lib.py | 6 +++-- .../wind_farm_models/engineering_models.py | 7 +++--- 15 files changed, 78 insertions(+), 26 deletions(-) diff --git a/py_wake/deficit_models/hybridinduction.py b/py_wake/deficit_models/hybridinduction.py index fcc82b100..02f8234dc 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 56b560464..a1bc71263 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 78139a5da..a4180174d 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 b02e3f4f5..640816137 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 e2fac2a48..b16f687a5 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 72c9a696c..b6ed3b75a 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 1c3d00dea..01cf73d63 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 4ba89ed3d..54c3c5cf9 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 526f2899b..de2a27bdb 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 fa05b05b0..c909f7d10 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 ef31fe09a..26c047969 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 f78cf4c74..e4d977afd 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 8ddae93b5..9fbeb0a77 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 2fba67d78..fb5f6f3bd 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 abf3a6768..68d9e59a9 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: -- GitLab