import numpy as np from py_wake.site.shear import PowerShear from py_wake.site.xrsite import XRSite import xarray as xr from py_wake.tests import npt import pytest import matplotlib.pyplot as plt from numpy import newaxis as na from py_wake.tests.test_files import tfp from py_wake.examples.data.hornsrev1 import Hornsrev1Site, V80, wt9_x, wt9_y from py_wake.examples.data.iea37._iea37 import IEA37_WindTurbines from py_wake.wind_turbines import WindTurbines from py_wake.deficit_models.gaussian import BastankhahGaussian f = [0.035972, 0.039487, 0.051674, 0.070002, 0.083645, 0.064348, 0.086432, 0.117705, 0.151576, 0.147379, 0.10012, 0.05166] A = [9.176929, 9.782334, 9.531809, 9.909545, 10.04269, 9.593921, 9.584007, 10.51499, 11.39895, 11.68746, 11.63732, 10.08803] k = [2.392578, 2.447266, 2.412109, 2.591797, 2.755859, 2.595703, 2.583984, 2.548828, 2.470703, 2.607422, 2.626953, 2.326172] ti = .1 @pytest.fixture def uniform_site(): ti = 0.1 ds = xr.Dataset( data_vars={'WS': 10, 'P': ('wd', f), 'TI': ti}, coords={'wd': np.linspace(0, 360, len(f), endpoint=False)}) return XRSite(ds, shear=PowerShear(h_ref=100, alpha=.2)) @pytest.fixture def uniform_weibull_site(): ti = 0.1 ds = xr.Dataset( data_vars={'Sector_frequency': ('wd', f), 'Weibull_A': ('wd', A), 'Weibull_k': ('wd', k), 'TI': ti}, coords={'wd': np.linspace(0, 360, len(f), endpoint=False)}) return XRSite(ds, shear=PowerShear(h_ref=100, alpha=.2)) @pytest.fixture def complex_fixed_pos_site(): ti = 0.1 ds = xr.Dataset( data_vars={'Speedup': ('i', np.arange(.8, 1.3, .1)), 'Turning': ('i', np.arange(-2, 3)), 'Sector_frequency': ('wd', f), 'Weibull_A': ('wd', A), 'Weibull_k': ('wd', k), 'TI': ti}, coords={'i': np.arange(5), 'wd': np.linspace(0, 360, len(f), endpoint=False)}) x_i = np.arange(5) return XRSite(ds, initial_position=np.array([x_i, x_i + 1]).T, shear=PowerShear(h_ref=100, alpha=.2)) @pytest.fixture def complex_grid_site(): ti = 0.1 ds = xr.Dataset( data_vars={'Speedup': (['x', 'y'], np.arange(.8, 1.4, .1).reshape((3, 2))), 'Turning': (['x', 'y'], np.arange(-2, 4, 1).reshape((3, 2))), 'Sector_frequency': ('wd', f), 'Weibull_A': ('wd', A), 'Weibull_k': ('wd', k), 'TI': ti}, coords={'x': [0, 5, 10], 'y': [0, 5], 'wd': np.linspace(0, 360, len(f), endpoint=False)}) return XRSite(ds, shear=PowerShear(h_ref=100, alpha=.2), interp_method='linear') def test_uniform_local_wind(uniform_site): site = uniform_site x_i = y_i = np.arange(5) wdir_lst = np.arange(0, 360, 90) wsp_lst = np.arange(3, 6) lw = site.local_wind(x_i=x_i, y_i=y_i, h_i=100, wd=wdir_lst, ws=wsp_lst) npt.assert_array_equal(lw.WS, 10) npt.assert_array_equal(lw.WD, wdir_lst) npt.assert_array_equal(lw.TI, 0.1) npt.assert_array_equal(lw.P, np.array([0.035972, 0.070002, 0.086432, 0.147379]) / 30) npt.assert_array_equal(site.elevation(x_i, y_i), 0) lw = site.local_wind(x_i=x_i, y_i=y_i, h_i=100) npt.assert_array_equal(lw.WD_ilk.shape, (1, 360, 1)) z = np.arange(1, 200) zero = [0] * len(z) ws100_2 = site.local_wind(x_i=zero, y_i=zero, h_i=z, wd=[0], ws=[10]).WS_ilk[:, 0, 0] site.shear = PowerShear(70, alpha=.3) ws70_3 = site.local_wind(x_i=zero, y_i=zero, h_i=z, wd=[0], ws=[10]).WS_ilk[:, 0, 0] if 0: plt.plot(ws100_2, z) plt.plot(ws70_3, z) plt.show() npt.assert_array_equal(10 * (z / 100)**.2, ws100_2) npt.assert_array_equal(10 * (z / 70)**.3, ws70_3) def test_uniform_weibull_local_wind(uniform_weibull_site): site = uniform_weibull_site x_i = y_i = np.arange(5) wdir_lst = np.arange(0, 360, 90) wsp_lst = np.arange(3, 6) lw = site.local_wind(x_i=x_i, y_i=y_i, h_i=100, wd=wdir_lst, ws=wsp_lst) npt.assert_array_equal(lw.WS, [3, 4, 5]) npt.assert_array_equal(lw.WD, wdir_lst) npt.assert_array_equal(lw.TI, 0.1) # ref_site = UniformWeibullSite(p_wd=f, a=A, k=k, ti=0.1, shear=site.shear) # lw_ref = ref_site.local_wind(x_i=x_i, y_i=y_i, h_i=100, wd=wdir_lst, ws=wsp_lst) # print(lw_ref.P) npt.assert_array_almost_equal(lw.P, [[6.14690773e-05, 8.55914207e-05, 1.05950962e-04], [8.72311046e-05, 1.30862281e-04, 1.73054727e-04], [1.17587820e-04, 1.75168689e-04, 2.29638490e-04], [1.19999621e-04, 1.84091959e-04, 2.50831917e-04]]) def test_complex_fixed_pos_local_wind(complex_fixed_pos_site): site = complex_fixed_pos_site x_i, y_i = site.initial_position.T npt.assert_array_equal(x_i, np.arange(5)) npt.assert_array_equal(y_i, np.arange(5) + 1) wdir_lst = np.arange(0, 360, 90) wsp_lst = np.arange(3, 6) lw = site.local_wind(x_i=x_i, y_i=y_i, h_i=100, wd=wdir_lst, ws=wsp_lst) npt.assert_array_equal(lw.WS, [3, 4, 5] * np.arange(.8, 1.3, .1)[:, na]) npt.assert_array_equal(lw.WD, (wdir_lst + np.arange(-2, 3)[:, na]) % 360) npt.assert_array_equal(lw.TI, 0.1) npt.assert_array_almost_equal(lw.P, [[6.14690773e-05, 8.55914207e-05, 1.05950962e-04], [8.72311046e-05, 1.30862281e-04, 1.73054727e-04], [1.17587820e-04, 1.75168689e-04, 2.29638490e-04], [1.19999621e-04, 1.84091959e-04, 2.50831917e-04]]) def test_complex_grid_local_wind(complex_grid_site): site = complex_grid_site y = np.arange(5) x = y * 2 X, Y = np.meshgrid(x, y) x_i, y_i = X.flatten(), Y.flatten() wdir_lst = np.arange(0, 360, 90) wsp_lst = np.arange(3, 6) lw = site.local_wind(x_i=x_i, y_i=y_i, h_i=100, wd=wdir_lst, ws=wsp_lst) if 0: c = plt.contourf(X, Y, lw.WS.sel(ws=5).data.reshape(X.shape)) plt.colorbar(c) plt.figure() c = plt.contourf(X, Y, lw.WD.sel(wd=90).data.reshape(X.shape)) plt.colorbar(c) plt.show() lw = site.local_wind(x_i=[2.5, 7.5], y_i=[2.5, 2.5], h_i=100, wd=wdir_lst, ws=wsp_lst) npt.assert_array_almost_equal(lw.WS, [3, 4, 5] * np.array([0.95, 1.15])[:, na]) npt.assert_array_equal(lw.WD, (wdir_lst + np.array([-.5, 1.5])[:, na]) % 360) npt.assert_array_equal(lw.TI, 0.1) npt.assert_array_almost_equal(lw.P, [[6.14690773e-05, 8.55914207e-05, 1.05950962e-04], [8.72311046e-05, 1.30862281e-04, 1.73054727e-04], [1.17587820e-04, 1.75168689e-04, 2.29638490e-04], [1.19999621e-04, 1.84091959e-04, 2.50831917e-04]]) def test_wd_independent_site(): ti = 0.1 ds = xr.Dataset( data_vars={ 'WS': 10, 'Sector_frequency': 1, 'Weibull_A': 4, 'Weibull_k': 2, 'TI': ti}, coords={}) site = XRSite(ds, shear=None) npt.assert_equal(site.ds.sector_width, 360) def test_load_save(complex_grid_site): complex_grid_site.save(tfp + "tmp.nc") site = XRSite.load(tfp + "tmp.nc", interp_method='linear') test_complex_grid_local_wind(site) def test_elevation(): ti = 0.1 ds = xr.Dataset( data_vars={'Elevation': (['x', 'y'], np.arange(.8, 1.4, .1).reshape((3, 2))), 'P': 1, 'TI': ti}, coords={'x': [0, 5, 10], 'y': [0, 5]}) site = XRSite(ds) npt.assert_array_almost_equal(site.elevation([2.5, 7.5], [2.5, 2.5]), [0.95, 1.15]) def test_cyclic_interpolation(uniform_site): site = uniform_site lw = site.local_wind([0], [0], 100, wd=np.arange(360), ws=10) if 0: plt.plot(np.linspace(0, 360, len(f) + 1), np.r_[f, f[0]], '.') plt.plot(lw.wd, lw.P) plt.show() npt.assert_array_almost_equal(lw.P[::30], np.array(f) / 30) def test_from_flow_box_2wt(): site = Hornsrev1Site() windTurbines = V80() # simulate current and neighbour wt wfm = BastankhahGaussian(site, windTurbines) wd = np.arange(30) sim_res = wfm([0, 0], [0, 500], wd=wd) ref_aep = sim_res.aep().sel(wt=0) wt_x, wt_y = [0], [0] neighbour_x, neighbour_y = [0], [500] # make site with effects of neighbour wt sim_res = wfm(neighbour_x, neighbour_y, wd=wd) e = 100 box = sim_res.flow_box(x=np.linspace(min(wt_x) - e, max(wt_x) + e, 21), y=np.linspace(min(wt_y) - e, max(wt_y) + e, 21), h=windTurbines.hub_height(windTurbines.types())) site = XRSite.from_flow_box(box) # Simujlate current wt and compare aep wfm = BastankhahGaussian(site, windTurbines) sim_res = wfm(wt_x, wt_y, wd=wd) aep = sim_res.aep() if 0: site.ds.WS.sel(ws=10, wd=3).plot() windTurbines.plot(wt_x, wt_y) windTurbines.plot(neighbour_x, neighbour_y) plt.show() npt.assert_array_almost_equal(ref_aep, aep.sel(wt=0))