Support timeseries with turbine-specific speeds

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"

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

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

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

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

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

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')

wifa/pywake_api.py

@@ -136,19 +136,39 @@ 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"]
+            # 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)
 
@@ -369,10 +389,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 +490,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:
@@ -791,18 +847,24 @@ def dict_to_site(resource_dict):
     # 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,
-        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')