From 607e524e10ff1844c0d0a353acf8432027280cbd Mon Sep 17 00:00:00 2001 From: Julian Quick Date: Tue, 16 Dec 2025 03:48:17 -0800 Subject: [PATCH 1/4] support timeseries with turbine-specific speeds --- .../FLOW_toy_study_energy_resource.yaml | 3 + .../Stochastic_atHubHeight.nc | Bin 0 -> 9608 bytes .../FLOW_toy_study_energy_site.yaml | 6 ++ .../DTU_10MW_turbine.yaml | 10 ++ .../FLOW_toy_study_wind_farm.yaml | 6 ++ .../wind_energy_system/analysis.yaml | 67 ++++++++++++ .../wind_energy_system/system.yaml | 30 ++++++ tests/test_pywake.py | 102 ++++++++++++++++++ wifa/pywake_api.py | 75 +++++++++++-- 9 files changed, 288 insertions(+), 11 deletions(-) create mode 100644 examples/cases/turbine_specific_speeds_timeseries/plant_energy_resource/FLOW_toy_study_energy_resource.yaml create mode 100644 examples/cases/turbine_specific_speeds_timeseries/plant_energy_resource/Stochastic_atHubHeight.nc create mode 100644 examples/cases/turbine_specific_speeds_timeseries/plant_energy_site/FLOW_toy_study_energy_site.yaml create mode 100644 examples/cases/turbine_specific_speeds_timeseries/plant_energy_turbine/DTU_10MW_turbine.yaml create mode 100644 examples/cases/turbine_specific_speeds_timeseries/plant_wind_farm/FLOW_toy_study_wind_farm.yaml create mode 100644 examples/cases/turbine_specific_speeds_timeseries/wind_energy_system/analysis.yaml create mode 100644 examples/cases/turbine_specific_speeds_timeseries/wind_energy_system/system.yaml diff --git a/examples/cases/turbine_specific_speeds_timeseries/plant_energy_resource/FLOW_toy_study_energy_resource.yaml b/examples/cases/turbine_specific_speeds_timeseries/plant_energy_resource/FLOW_toy_study_energy_resource.yaml new file mode 100644 index 0000000..3866df6 --- /dev/null +++ b/examples/cases/turbine_specific_speeds_timeseries/plant_energy_resource/FLOW_toy_study_energy_resource.yaml @@ -0,0 +1,3 @@ +name: FLOW UQ vnv study on toy problem, 4 WT Wind Farm +wind_resource: + !include "Stochastic_atHubHeight.nc" diff --git a/examples/cases/turbine_specific_speeds_timeseries/plant_energy_resource/Stochastic_atHubHeight.nc b/examples/cases/turbine_specific_speeds_timeseries/plant_energy_resource/Stochastic_atHubHeight.nc new file mode 100644 index 0000000000000000000000000000000000000000..849c248c18d3f0779b03a1ea985d53234fc64d08 GIT binary patch literal 9608 zcmeHLO>7%g5T5ngTg%y`38V_Yii=E;phO8Zm8xQD%XUbF>y#K(6$G4+-NH(#6WN=h zAS9YrYPRnYpc|$ETb1%$D!pOK2TJNW$(4T_LIeZ@j;+Zhryp&x9#pj>qQIp!z6lZeKrwSHFXSUip1>Dc zO%&=u<4`Lbm-9}0&StalDH?W~Of{hk_d>fF^XBJI%+4>lEioj-c@sT-wapH`kg21UVARb5E}xnaqyAXy-|aoG3!0T+AM?E-boL-^17< zR498%@fo*T^A;+pXAV?E5T=sDhldXyc*{BUY-)6PG&y|ev-)8Y#<6k~f71Ae*9u;r zJ51TgH@OCCg|W-zzE4zy9~9*MxVlm?rOVZ1z*(yfp7t=G?MZ@o71K z`JX=du!UzL5~ec{4ZFG2{y;6g@2i{7w>r3Nwdj`%%U}HeTsayXX_8%#yp2U%^U+3{ zw_fLb`BZJMA#VhSn_bM8OXM*ee@FtbhyOu(HFa`Sc;iuG8ey(UV*=pXqU$>BKoOUK zC!9@O<=uzttKA|`auk!T^W;wCT6s$vw_d-uxkB09oGWB@+qu$pasLk;$rZWhja=dR zM?+l)Wu#v2=v#eS<<*Q@cvs)*)S1{WWbW+lv<{iZJq{O2ndWez*HGKxa%t_xBIN{q zU?T8qLW0W;sJ7I}m**$vm;aPXTK>L}oh%4JmGJ8OrruWo+_6*N&w$VGPxS>1_$;qd zJm9muS|;GLyn1H`e3oaUXILRKG}Q=b1T+E~0gZr0KqH_L&i#vr(xy#Dt25m#d;8rf*f3^7ATC>I9 dsUB2Z36RCsu2?LaD;wQ6p0mX~&v><4{{iOwUM>It literal 0 HcmV?d00001 diff --git a/examples/cases/turbine_specific_speeds_timeseries/plant_energy_site/FLOW_toy_study_energy_site.yaml b/examples/cases/turbine_specific_speeds_timeseries/plant_energy_site/FLOW_toy_study_energy_site.yaml new file mode 100644 index 0000000..ee4c593 --- /dev/null +++ b/examples/cases/turbine_specific_speeds_timeseries/plant_energy_site/FLOW_toy_study_energy_site.yaml @@ -0,0 +1,6 @@ +name: FLOW UQ vnv study on toy problem, 4 WT Wind Farm +boundaries: + polygons: + - x: [-90, 3834,3, 3834,3, -90] + y: [90, 90, -90, -90] +energy_resource: !include ../plant_energy_resource/FLOW_toy_study_energy_resource.yaml diff --git a/examples/cases/turbine_specific_speeds_timeseries/plant_energy_turbine/DTU_10MW_turbine.yaml b/examples/cases/turbine_specific_speeds_timeseries/plant_energy_turbine/DTU_10MW_turbine.yaml new file mode 100644 index 0000000..0320701 --- /dev/null +++ b/examples/cases/turbine_specific_speeds_timeseries/plant_energy_turbine/DTU_10MW_turbine.yaml @@ -0,0 +1,10 @@ +name: DTU 10MW Offshore Reference Turbine +performance: + power_curve: + power_values: [263388., 751154., 1440738., 2355734., 3506858., 4993092., 6849310., 9116402., 10000754., 10009590., 10000942., 10042678., 10003480., 10001600., 10001506., 10013632., 10007428., 10005360., 10002728., 10001130., 10004984., 9997558.] + power_wind_speeds: [4.,5.,6.,7.,8.,9.,10.,11.,12.,13.,14.,15.,16.,17.,18.,19.,20.,21.,22.,23.,24.,25.] + Ct_curve: + Ct_values: [0.923,0.919,0.904,0.858,0.814,0.814,0.814,0.814,0.577,0.419,0.323,0.259,0.211,0.175,0.148,0.126,0.109,0.095,0.084,0.074,0.066,0.059] + Ct_wind_speeds: [4.,5.,6.,7.,8.,9.,10.,11.,12.,13.,14.,15.,16.,17.,18.,19.,20.,21.,22.,23.,24.,25.] +hub_height: 119.0 +rotor_diameter: 178.3 diff --git a/examples/cases/turbine_specific_speeds_timeseries/plant_wind_farm/FLOW_toy_study_wind_farm.yaml b/examples/cases/turbine_specific_speeds_timeseries/plant_wind_farm/FLOW_toy_study_wind_farm.yaml new file mode 100644 index 0000000..25cf13e --- /dev/null +++ b/examples/cases/turbine_specific_speeds_timeseries/plant_wind_farm/FLOW_toy_study_wind_farm.yaml @@ -0,0 +1,6 @@ +name: FLOW UQ vnv study on toy problem, 4 WT Wind Farm +layouts: + - coordinates: + x: [0, 1248.1, 2496.2, 3744.3] + y: [0, 0, 0, 0] +turbines: !include ../plant_energy_turbine/DTU_10MW_turbine.yaml diff --git a/examples/cases/turbine_specific_speeds_timeseries/wind_energy_system/analysis.yaml b/examples/cases/turbine_specific_speeds_timeseries/wind_energy_system/analysis.yaml new file mode 100644 index 0000000..1278689 --- /dev/null +++ b/examples/cases/turbine_specific_speeds_timeseries/wind_energy_system/analysis.yaml @@ -0,0 +1,67 @@ +#pywake and foxes +wind_deficit_model: + name: Bastankhah2014 + wake_expansion_coefficient: # k = ka*ti + kb + k_a: 0.0 + k_b: 0.04 + free_stream_ti: false + ceps: 0.2 + use_effective_ws: true +axial_induction_model: Madsen +deflection_model: + name: None +turbulence_model: + name: CrespoHernandez +superposition_model: + ws_superposition: Linear + ti_superposition: Squared +rotor_averaging: + grid: grid + n_x_grid_points: 4 + n_y_grid_points: 4 + background_averaging: center + wake_averaging: center + wind_speed_exponent_for_power: 3 + wind_speed_exponent_for_ct: 2 +blockage_model: + name: None + +#wayve +layers_description: + farm_layer_height: 238. + number_of_fa_layers: 1 +APM_additional_terms: + momentum_entrainment: # Momentum flux parametrization + mfp_type: "constant_flux" # Options: None, "constant_flux". Default: None + apm_mfp_settings: # Only used when mfp_type is "constant_flux" + a_mfp: 0.120 # Default: 0.120 + d_mfp: 27.80 # Default: 27.80 +apm_grid: + Lx: 1.e6 + Ly: 1.e6 + dx: 500 + L_filter: 1.e3 +#wake_tool: "wayve" # Options: "wayve", "foxes" (pywake to do). Default: "wayve" +wm_coupling: + method: "PB" # Options: VM, PB, US (see https://dx.doi.org/10.1088/1742-6596/2767/9/092079) + +#code_saturne +run_type: "simulate" #"postprocess" +# +HPC_config: + run_node_number: 1 + run_ntasks_per_node: 48 + run_wall_time_hours: 6 + run_partition: "" + # + mesh_node_number: 1 + mesh_ntasks_per_node: 48 + mesh_wall_time_hours: 1 + #mesh_partition: "" + # + wckey: "" + + +mesh: + remesh: True + mesh_file_name: MESH/toy_mesh.med diff --git a/examples/cases/turbine_specific_speeds_timeseries/wind_energy_system/system.yaml b/examples/cases/turbine_specific_speeds_timeseries/wind_energy_system/system.yaml new file mode 100644 index 0000000..69a21c1 --- /dev/null +++ b/examples/cases/turbine_specific_speeds_timeseries/wind_energy_system/system.yaml @@ -0,0 +1,30 @@ +name: FLOW UQ vnv study on toy problem, 4 WT Wind Farm +site: !include ../plant_energy_site/FLOW_toy_study_energy_site.yaml +wind_farm: !include ../plant_wind_farm/FLOW_toy_study_wind_farm.yaml +attributes: + flow_model: + name: foxes + analysis: !include analysis.yaml + + model_outputs_specification: + output_folder: "results" + # + run_configuration: + times_run: + all_occurences: True + # + turbine_outputs: + turbine_nc_filename: 'turbine_data.nc' # dimension = states, turbine + output_variables: ['power', 'rotor_effective_velocity'] #'frequency' + # + flow_field: + report: false + flow_nc_filename: flow_field.nc + output_variables: ['wind_speed', 'wind_direction'] + z_planes: + z_sampling: "hub_heights" + xy_sampling: "grid" + x_bounds: [-1000, 5000] + y_bounds: [-1000, 1000] + dx: 150 + dy: 150 diff --git a/tests/test_pywake.py b/tests/test_pywake.py index ef3a809..3c01eb2 100644 --- a/tests/test_pywake.py +++ b/tests/test_pywake.py @@ -11,6 +11,7 @@ from py_wake.site import XRSite from py_wake.superposition_models import LinearSum from py_wake.tests import npt +from py_wake.turbulence_models import CrespoHernandez from py_wake.wind_turbines import WindTurbine from py_wake.wind_turbines.power_ct_functions import PowerCtFunctionList, PowerCtTabular from scipy.special import gamma @@ -339,6 +340,107 @@ def test_heterogeneous_wind_rose_arbitrary_points(): assert wifa_res == res_aep +def test_turbine_specific_speeds_timeseries(): + """ + Test case for time-series simulation where inflow wind speed/direction + is specific to each turbine (dimensions: time x turbine). + Validates that the API correctly reduces 2D inputs to 1D reference arrays + for the simulation call while preserving local site data. + """ + from windIO import load_yaml # Import helper to read YAML + + case_name = "turbine_specific_speeds_timeseries" + system_yaml = ( + test_path / f"../examples/cases/{case_name}/wind_energy_system/system.yaml" + ) + resource_nc = ( + test_path + / f"../examples/cases/{case_name}/plant_energy_resource/Stochastic_atHubHeight.nc" + ) + + # 1. Run via API + wifa_res = run_pywake(system_yaml) + + # 2. Run manually to verify logic + + # Load coordinates from YAML + sys_dat = load_yaml(system_yaml) + farm_layout = sys_dat["wind_farm"]["layouts"] + if isinstance(farm_layout, list): + coords = farm_layout[0]["coordinates"] + else: + coords = farm_layout["coordinates"] + x = coords["x"] + y = coords["y"] + + # Load and prepare Site Data + ds = xr.open_dataset(resource_nc) + ds = ds.rename( + { + "wind_direction": "WD", + "wind_speed": "WS", + "wind_turbine": "i", + "turbulence_intensity": "TI", + } + ) + + # FIX 1: Re-index 'i' to be 0-based to match PyWake's internal numbering + # The file has [1, 2, 3, 4], PyWake expects [0, 1, 2, 3] + ds = ds.assign_coords(i=np.arange(len(ds.i))) + + # FIX 2: Transpose to (i, time) for XRSite linear interpolator + ds = ds.transpose("i", "time") + + # FIX 3: Add uniform probability 'P' + n_time = len(ds.time) + ds["P"] = (("time"), np.ones(n_time) / n_time) + + # Initialize Site + site = XRSite(ds, interp_method="linear") + + # Define Turbine + turbine = WindTurbine( + name="test", + diameter=178.3, + hub_height=119.0, + powerCtFunction=POWER_CT_TABLE, + ) + + # Define Wake Model + wfm = BastankhahGaussian( + site, + turbine, + k=0.04, + ceps=0.2, + superpositionModel=LinearSum(), + use_effective_ws=True, + turbulenceModel=CrespoHernandez(), + ) + + # Calculate reference arrays (API Logic) + if "i" in ds.WS.dims: + ws_ref = ds.WS.mean(dim="i").values + else: + ws_ref = ds.WS.values + + if "i" in ds.WD.dims: + rads = np.deg2rad(ds.WD) + mean_sin = np.sin(rads).mean(dim="i") + mean_cos = np.cos(rads).mean(dim="i") + wd_ref = np.rad2deg(np.arctan2(mean_sin, mean_cos)) % 360 + wd_ref = wd_ref.values + else: + wd_ref = ds.WD.values + + # Manual Simulation + res_manual = wfm(x, y, time=ds.time, ws=ws_ref, wd=wd_ref) + + manual_aep = res_manual.aep(normalize_probabilities=False).sum() + + # 3. Assert match + npt.assert_allclose(wifa_res, manual_aep, rtol=1e-6) + + # if __name__ == "__main__": # test_heterogeneous_wind_rose() # simple_yaml_to_pywake('../examples/cases/windio_4turbines_multipleTurbines/plant_energy_turbine/IEA_10MW_turbine.yaml') diff --git a/wifa/pywake_api.py b/wifa/pywake_api.py index e9a2f98..bb6a99f 100644 --- a/wifa/pywake_api.py +++ b/wifa/pywake_api.py @@ -136,19 +136,36 @@ def run_pywake(yamlFile, output_dir="output"): def dict_to_site(resource_dict): resource_ds = dict_to_netcdf(resource_dict) - to_rename = { + rename_map = { "height": "h", - "wind_direction": "wd", - "wind_speed": "ws", "weibull_a": "Weibull_A", "weibull_k": "Weibull_k", "sector_probability": "Sector_frequency", "turbulence_intensity": "TI", "wind_turbine": "i", } - for name in to_rename: + + # Smart rename for wind_direction and wind_speed + for key, coord_name, var_name in [ + ("wind_direction", "wd", "WD"), + ("wind_speed", "ws", "WS"), + ]: + if key in resource_ds: + # If it's a coordinate (dimension), use lowercase (wd, ws) + # If it's a data variable (time series/map), use uppercase (WD, WS) + rename_map[key] = coord_name if key in resource_ds.coords else var_name + + for name in rename_map: if name in resource_ds: - resource_ds = resource_ds.rename({name: to_rename[name]}) + resource_ds = resource_ds.rename({name: rename_map[name]}) + + if "P" not in resource_ds and "time" in resource_ds.dims: + n_time = len(resource_ds.time) + # Create uniform probability array (1/N) + resource_ds["P"] = (("time",), np.ones(n_time) / n_time) + if "i" in resource_ds.dims: + other_dims = [d for d in resource_ds.dims if d != "i"] + resource_ds = resource_ds.transpose("i", *other_dims) print("making site with ", resource_ds) return XRSite(resource_ds) @@ -369,10 +386,40 @@ def dict_to_site(resource_dict): heights = resource_dat["wind_resource"]["height"] else: heights = None - ws = np.array(resource_dat["wind_resource"]["wind_speed"]["data"])[cases_idx] - wd = np.array(resource_dat["wind_resource"]["wind_direction"]["data"])[ - cases_idx - ] + + # Helper to get data and dimensions safely + def get_resource_data(var_name): + data_obj = resource_dat["wind_resource"][var_name] + # Handle windIO schema structure (data + dims) + vals = np.array(data_obj["data"]) + # Default to ["time"] if dims missing but data exists + dims = data_obj.get("dims", ["time"]) + return vals, dims + + # Extract Raw Data + ws_vals, ws_dims = get_resource_data("wind_speed") + wd_vals, wd_dims = get_resource_data("wind_direction") + + # Apply subsetting (cases_idx) + # Assuming 'time' is always the first dimension + ws_vals = ws_vals[cases_idx] + wd_vals = wd_vals[cases_idx] + + # PREPARE REFERENCE ARRAYS FOR SIMULATION CALL + # If data is turbine specific (2D: time x turbine), flatten to 1D (time) + if "wind_turbine" in ws_dims: + ws = np.mean(ws_vals, axis=1) + else: + ws = ws_vals + + if "wind_turbine" in wd_dims: + # Vector mean for direction to handle 360/0 boundary + rads = np.deg2rad(wd_vals) + mean_sin = np.mean(np.sin(rads), axis=1) + mean_cos = np.mean(np.cos(rads), axis=1) + wd = np.mod(np.rad2deg(np.arctan2(mean_sin, mean_cos)), 360) + else: + wd = wd_vals if "operating" in resource_dat["wind_resource"].keys(): operating = np.array(resource_dat["wind_resource"]["operating"]["data"])[ @@ -440,7 +487,13 @@ def dict_to_site(resource_dict): )(hh) assert len(np.array(times)[cases_idx]) == len(ws) assert len(wd) == len(ws) - site = Hornsrev1Site() + if "wind_turbine" in ws_dims or "wind_turbine" in wd_dims: + # Turbine-specific data -> Use XRSite via dict_to_site + site = dict_to_site(resource_dat["wind_resource"]) + else: + # Uniform data -> Use Hornsrev1Site (Old behavior) + # This preserves behavior for KUL and 4wts tests which use uniform inflow + site = Hornsrev1Site() if "turbulence_intensity" not in resource_dat["wind_resource"]: TI = 0.02 else: @@ -798,7 +851,7 @@ def dict_to_site(resource_dict): time=timeseries, ws=ws, wd=wd, - TI=TI, + # TI=TI, # TI is defined in XRSite yaw=0, tilt=0, operating=operating, From ee975bc7ee90e6c9a6ccc6d95a86793c79c8a159 Mon Sep 17 00:00:00 2001 From: Julian Quick Date: Tue, 16 Dec 2025 05:54:12 -0800 Subject: [PATCH 2/4] better comments for transpose Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- wifa/pywake_api.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/wifa/pywake_api.py b/wifa/pywake_api.py index bb6a99f..1fff792 100644 --- a/wifa/pywake_api.py +++ b/wifa/pywake_api.py @@ -165,6 +165,9 @@ def dict_to_site(resource_dict): resource_ds["P"] = (("time",), np.ones(n_time) / n_time) if "i" in resource_ds.dims: other_dims = [d for d in resource_ds.dims if d != "i"] + # The transpose operation ensures that 'i' (turbine index) is the first dimension. + # This is required for XRSite's linear interpolation, which expects the turbine index + # as the leading dimension. See test at line 392 for details. resource_ds = resource_ds.transpose("i", *other_dims) print("making site with ", resource_ds) return XRSite(resource_ds) From c7eb56ff32da3d7801b16ce366531df2efb8fa86 Mon Sep 17 00:00:00 2001 From: Julian Quick Date: Tue, 16 Dec 2025 05:55:18 -0800 Subject: [PATCH 3/4] correct commas/decimals in examples/cases/turbine_specific_speeds_timeseries/plant_energy_site/FLOW_toy_study_energy_site.yaml Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- .../plant_energy_site/FLOW_toy_study_energy_site.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/cases/turbine_specific_speeds_timeseries/plant_energy_site/FLOW_toy_study_energy_site.yaml b/examples/cases/turbine_specific_speeds_timeseries/plant_energy_site/FLOW_toy_study_energy_site.yaml index ee4c593..e76c696 100644 --- a/examples/cases/turbine_specific_speeds_timeseries/plant_energy_site/FLOW_toy_study_energy_site.yaml +++ b/examples/cases/turbine_specific_speeds_timeseries/plant_energy_site/FLOW_toy_study_energy_site.yaml @@ -1,6 +1,6 @@ name: FLOW UQ vnv study on toy problem, 4 WT Wind Farm boundaries: polygons: - - x: [-90, 3834,3, 3834,3, -90] + - x: [-90, 3834.3, 3834.3, -90] y: [90, 90, -90, -90] energy_resource: !include ../plant_energy_resource/FLOW_toy_study_energy_resource.yaml From 4d0a2256df6617bc1c349231bf705d920e544af2 Mon Sep 17 00:00:00 2001 From: Julian Quick Date: Tue, 16 Dec 2025 13:39:43 -0800 Subject: [PATCH 4/4] ensure TI timeseries are properly handled --- wifa/pywake_api.py | 30 ++++++++++++++++++------------ 1 file changed, 18 insertions(+), 12 deletions(-) diff --git a/wifa/pywake_api.py b/wifa/pywake_api.py index 1fff792..1c42edb 100644 --- a/wifa/pywake_api.py +++ b/wifa/pywake_api.py @@ -847,18 +847,24 @@ def get_resource_data(var_name): # noj = NOJ(site, turbine, turbulenceModel=None) # sim_res = noj(x, y) # sim_res = windFarmModel(x, y, type=turbine_types, time=timeseries, ws=ws, wd=wd, TI=0, yaw=0, tilt=0) - sim_res = windFarmModel( - x, - y, - type=turbine_types, - time=timeseries, - ws=ws, - wd=wd, - # TI=TI, # TI is defined in XRSite - yaw=0, - tilt=0, - operating=operating, - ) + sim_kwargs = { + "x": x, + "y": y, + "type": turbine_types, + "time": timeseries, + "ws": ws, + "wd": wd, + "yaw": 0, + "tilt": 0, + "operating": operating, + } + + # Check if TI is missing from the site's data variables + # If it's missing (common in Hornsrev1Site), pass the local TI variable + if "TI" not in site.ds.data_vars: + sim_kwargs["TI"] = TI + + sim_res = windFarmModel(**sim_kwargs) aep = sim_res.aep(normalize_probabilities=not timeseries).sum() print("aep is ", aep, "GWh") # print('aep is ', sim_res.aep().sum(), 'GWh')