From b9a1b22feaa50c95e5a492f6369995beeeeaeec4 Mon Sep 17 00:00:00 2001 From: Bjarke Tobias Olsen Date: Fri, 13 Mar 2026 12:36:09 +0100 Subject: [PATCH 01/10] Expand PyWake submodel support with case-insensitive lookup, optional deps, and NOJLocalDeficit alias - Make all wrapped tools (floris, pywake, wayve, foxes) optional dependencies - Add case-insensitive submodel name lookup across deficit, deflection, turbulence, superposition, rotor averaging, and blockage models - Add NOJLocalDeficit as a deficit model name alias alongside Jensen - Fix CI: skip floris on Python <3.10, fix numpy 2.x compat in wayve - Add comprehensive parametrized tests for all submodel configuration functions Co-Authored-By: Claude Opus 4.6 --- pyproject.toml | 1 + tests/test_floris.py | 1 + tests/test_pywake_submodels.py | 566 ++++++++++++++++++++++++++++++++ wifa/pywake_api.py | 577 +++++++++++++++++---------------- wifa/wayve_api.py | 2 +- 5 files changed, 875 insertions(+), 272 deletions(-) create mode 100644 tests/test_pywake_submodels.py diff --git a/pyproject.toml b/pyproject.toml index 1205899..06b5f4c 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -70,6 +70,7 @@ dev = [ "ncplot", "nctoolkit", "cartopy", + "pre-commit", ] docs = [ "sphinx>=7.0", diff --git a/tests/test_floris.py b/tests/test_floris.py index f95b706..f2477fc 100644 --- a/tests/test_floris.py +++ b/tests/test_floris.py @@ -9,6 +9,7 @@ import os import shutil +import sys from pathlib import Path import numpy as np diff --git a/tests/test_pywake_submodels.py b/tests/test_pywake_submodels.py new file mode 100644 index 0000000..8c61c27 --- /dev/null +++ b/tests/test_pywake_submodels.py @@ -0,0 +1,566 @@ +"""Parametrized unit tests for PyWake submodel configuration functions. + +Tests each _configure_*() function in wifa/pywake_api.py to verify: +- Correct PyWake class is returned for each model name +- Parameters are passed through correctly +- NotImplementedError raised for unsupported names +- Case-insensitive matching works +""" + +import pytest +from py_wake.deficit_models import ( + HybridInduction, + RankineHalfBody, + SelfSimilarityDeficit, + SelfSimilarityDeficit2020, + VortexCylinder, + VortexDipole, +) +from py_wake.deficit_models.gaussian import ( + BastankhahGaussianDeficit, + BlondelSuperGaussianDeficit2020, + BlondelSuperGaussianDeficit2023, + CarbajofuertesGaussianDeficit, + NiayifarGaussianDeficit, + TurboGaussianDeficit, + ZongGaussianDeficit, +) +from py_wake.deficit_models.gcl import GCLDeficit +from py_wake.deficit_models.noj import NOJLocalDeficit, TurboNOJDeficit +from py_wake.deficit_models.rathmann import Rathmann +from py_wake.deflection_models import JimenezWakeDeflection +from py_wake.deflection_models.gcl_hill_vortex import GCLHillDeflection +from py_wake.rotor_avg_models import ( + CGIRotorAvg, + EqGridRotorAvg, + GQGridRotorAvg, + GridRotorAvg, + PolarGridRotorAvg, + RotorCenter, +) +from py_wake.superposition_models import ( + CumulativeWakeSum, + LinearSum, + MaxSum, + SquaredSum, + WeightedSum, +) +from py_wake.turbulence_models import ( + CrespoHernandez, + STF2005TurbulenceModel, + STF2017TurbulenceModel, +) +from py_wake.turbulence_models.gcl_turb import GCLTurbulence + +from wifa.pywake_api import ( + DEFAULTS, + _configure_blockage_model, + _configure_deficit_model, + _configure_deflection_model, + _configure_rotor_averaging, + _configure_superposition_model, + _configure_turbulence_model, + configure_wake_model, + get_with_default, +) + +# Default rotor diameter and hub height for deficit model tests +_RD = 126.0 +_HH = 90.0 + + +def _call_deficit(name, analysis_extra=None): + """Helper to call _configure_deficit_model with minimal boilerplate.""" + wind_deficit_model = {"name": name, **(analysis_extra or {})} + analysis = {"wind_deficit_model": wind_deficit_model} + return _configure_deficit_model({"name": name}, analysis, _RD, _HH) + + +# --------------------------------------------------------------------------- +# Deficit model tests +# --------------------------------------------------------------------------- + + +@pytest.mark.parametrize( + "name,expected_class", + [ + ("Jensen", NOJLocalDeficit), + ("jensen", NOJLocalDeficit), + ("JENSEN", NOJLocalDeficit), + ("Bastankhah2014", BastankhahGaussianDeficit), + ("bastankhah2014", BastankhahGaussianDeficit), + ("BASTANKHAH2014", BastankhahGaussianDeficit), + ("SuperGaussian", BlondelSuperGaussianDeficit2020), + ("supergaussian", BlondelSuperGaussianDeficit2020), + ("SuperGaussian2023", BlondelSuperGaussianDeficit2023), + ("TurboPark", TurboGaussianDeficit), + ("TurbOPark", TurboGaussianDeficit), + ("turbopark", TurboGaussianDeficit), + ("Niayifar2016", NiayifarGaussianDeficit), + ("niayifar2016", NiayifarGaussianDeficit), + ("Zong2020", ZongGaussianDeficit), + ("zong2020", ZongGaussianDeficit), + ("Carbajofuertes2018", CarbajofuertesGaussianDeficit), + ("carbajofuertes2018", CarbajofuertesGaussianDeficit), + ("TurboNOJ", TurboNOJDeficit), + ("turbonoj", TurboNOJDeficit), + ("GCL", GCLDeficit), + ("gcl", GCLDeficit), + ("NOJLocalDeficit", NOJLocalDeficit), + ("nojlocaldeficit", NOJLocalDeficit), + ("NOJLOCALDEFICIT", NOJLocalDeficit), + ], +) +def test_configure_deficit_model(name, expected_class): + cls, args = _call_deficit(name) + assert cls is expected_class + + +@pytest.mark.parametrize( + "name,expected_class", + [ + ("Jensen", NOJLocalDeficit), + ("Bastankhah2014", BastankhahGaussianDeficit), + ("SuperGaussian", BlondelSuperGaussianDeficit2020), + ("SuperGaussian2023", BlondelSuperGaussianDeficit2023), + ("TurboPark", TurboGaussianDeficit), + ("Niayifar2016", NiayifarGaussianDeficit), + ("Zong2020", ZongGaussianDeficit), + ("Carbajofuertes2018", CarbajofuertesGaussianDeficit), + ("TurboNOJ", TurboNOJDeficit), + ("GCL", GCLDeficit), + ("NOJLocalDeficit", NOJLocalDeficit), + ], +) +def test_configure_deficit_model_instantiation(name, expected_class): + """Verify returned kwargs can actually instantiate the model without TypeError.""" + cls, args = _call_deficit(name) + instance = cls(**args) + assert isinstance(instance, expected_class) + + +def test_configure_deficit_model_bastankhah2014_params(): + """Verify wake expansion and ceps params are passed for Bastankhah2014.""" + cls, args = _call_deficit( + "Bastankhah2014", + {"wake_expansion_coefficient": {"k": 0.04}, "ceps": 0.2}, + ) + assert cls is BastankhahGaussianDeficit + assert args["k"] == 0.04 + assert args["ceps"] == 0.2 + + +def test_configure_deficit_model_bastankhah2014_k_b(): + """Verify k_b wake expansion is used when present.""" + _, args = _call_deficit( + "Bastankhah2014", + {"wake_expansion_coefficient": {"k_a": 0.38, "k_b": 0.004}}, + ) + assert args["k"] == 0.004 + + +def test_configure_deficit_model_jensen_k_b(): + """Verify Jensen k_a/k_b expansion params.""" + _, args = _call_deficit( + "Jensen", + {"wake_expansion_coefficient": {"k_a": 0.38, "k_b": 0.004}}, + ) + assert args["a"] == [0.38, 0.004] + + +def test_configure_deficit_model_nojlocaldeficit_k_b(): + """Verify NOJLocalDeficit k_a/k_b expansion params.""" + _, args = _call_deficit( + "NOJLocalDeficit", + {"wake_expansion_coefficient": {"k_a": 0.38, "k_b": 0.004}}, + ) + assert args["a"] == [0.38, 0.004] + + +def test_configure_deficit_model_gaussian_params_niayifar(): + """Verify Gaussian params pass through for Niayifar2016.""" + cls, args = _call_deficit( + "Niayifar2016", + { + "wake_expansion_coefficient": {"k_a": 0.38, "k_b": 0.004}, + "ceps": 0.3, + "use_effective_ti": True, + }, + ) + assert cls is NiayifarGaussianDeficit + assert args["a"] == [0.38, 0.004] + assert args["ceps"] == 0.3 + assert args["use_effective_ti"] is True + + +def test_configure_deficit_model_zong_no_ceps(): + """Verify Zong2020 does not pass ceps (unsupported).""" + _, args = _call_deficit( + "Zong2020", + {"ceps": 0.3, "use_effective_ti": False}, + ) + assert "ceps" not in args + assert args["use_effective_ti"] is False + + +def test_configure_deficit_model_bastankhah2014_no_effective_ti(): + """Verify Bastankhah2014 does not pass use_effective_ti (unsupported).""" + _, args = _call_deficit( + "Bastankhah2014", + {"wake_expansion_coefficient": {"k": 0.04}, "use_effective_ti": True}, + ) + assert "use_effective_ti" not in args + assert args["k"] == 0.04 + + +def test_configure_deficit_model_a_param_warns_on_scalar_k(): + """Verify warning when scalar k is provided for a=[k_a, k_b] models.""" + with pytest.warns(UserWarning, match="uses a="): + _, args = _call_deficit( + "Niayifar2016", {"wake_expansion_coefficient": {"k": 0.05}} + ) + assert "k" not in args + assert "a" not in args + + +def test_configure_deficit_model_a_param_warns_on_missing_k_a(): + """Verify warning when k_b is provided without k_a.""" + with pytest.warns(UserWarning, match="k_a not specified"): + _, args = _call_deficit( + "Zong2020", {"wake_expansion_coefficient": {"k_b": 0.004}} + ) + assert args["a"] == [0, 0.004] + + +@pytest.mark.parametrize( + "name,extra,expected_class", + [ + ( + "Bastankhah2014", + {"wake_expansion_coefficient": {"k": 0.04}, "ceps": 0.2}, + BastankhahGaussianDeficit, + ), + ( + "Niayifar2016", + { + "wake_expansion_coefficient": {"k_a": 0.38, "k_b": 0.004}, + "ceps": 0.3, + "use_effective_ti": True, + }, + NiayifarGaussianDeficit, + ), + ("Zong2020", {"use_effective_ti": False}, ZongGaussianDeficit), + ( + "Jensen", + {"wake_expansion_coefficient": {"k_a": 0.38, "k_b": 0.004}}, + NOJLocalDeficit, + ), + ( + "NOJLocalDeficit", + {"wake_expansion_coefficient": {"k_a": 0.38, "k_b": 0.004}}, + NOJLocalDeficit, + ), + ], +) +def test_configure_deficit_model_instantiation_with_params(name, extra, expected_class): + """Verify models with user-specified params can be instantiated.""" + cls, args = _call_deficit(name, extra) + assert isinstance(cls(**args), expected_class) + + +def test_configure_deficit_model_turbonoj_A_param(): + """Verify TurboNOJ passes through the A parameter and can be instantiated.""" + cls, args = _call_deficit("TurboNOJ", {"A": 0.6}) + assert cls is TurboNOJDeficit + assert args["A"] == 0.6 + instance = cls(**args) + assert isinstance(instance, TurboNOJDeficit) + + +@pytest.mark.parametrize("name", ["Bastankhah2016", "bastankhah2016"]) +def test_configure_deficit_model_bastankhah2016_not_implemented(name): + with pytest.raises(NotImplementedError, match="Bastankhah2016"): + _call_deficit(name) + + +def test_configure_deficit_model_unknown(): + with pytest.raises(NotImplementedError, match="NonexistentModel"): + _call_deficit("NonexistentModel") + + +# --------------------------------------------------------------------------- +# Deflection model tests +# --------------------------------------------------------------------------- + + +@pytest.mark.parametrize( + "name,expected_class", + [ + ("Jimenez", JimenezWakeDeflection), + ("jimenez", JimenezWakeDeflection), + ("JIMENEZ", JimenezWakeDeflection), + ("GCLHill", GCLHillDeflection), + ("gclhill", GCLHillDeflection), + ("GCLhill", GCLHillDeflection), + ], +) +def test_configure_deflection_model(name, expected_class): + model = _configure_deflection_model({"name": name, "beta": 0.1}) + assert isinstance(model, expected_class) + + +@pytest.mark.parametrize("name", [None, "None", "none", "NONE"]) +def test_configure_deflection_model_none(name): + assert _configure_deflection_model({"name": name, "beta": 0.1}) is None + + +def test_configure_deflection_model_bastankhah2016(): + with pytest.raises(NotImplementedError, match="Bastankhah2016"): + _configure_deflection_model({"name": "Bastankhah2016", "beta": 0.1}) + + +def test_configure_deflection_model_unknown(): + with pytest.raises(NotImplementedError, match="UnknownDeflection"): + _configure_deflection_model({"name": "UnknownDeflection", "beta": 0.1}) + + +# --------------------------------------------------------------------------- +# Turbulence model tests +# --------------------------------------------------------------------------- + + +@pytest.mark.parametrize( + "name,expected_class", + [ + ("STF2005", STF2005TurbulenceModel), + ("stf2005", STF2005TurbulenceModel), + ("STF2017", STF2017TurbulenceModel), + ("stf2017", STF2017TurbulenceModel), + ("CrespoHernandez", CrespoHernandez), + ("crespohernandez", CrespoHernandez), + ("CRESPOHERNANDEZ", CrespoHernandez), + ("IEC-TI-2019", STF2017TurbulenceModel), + ("iec-ti-2019", STF2017TurbulenceModel), + ("GCL", GCLTurbulence), + ("gcl", GCLTurbulence), + ], +) +def test_configure_turbulence_model(name, expected_class): + data = {"name": name, "c1": 1.0, "c2": 1.0} + assert isinstance(_configure_turbulence_model(data), expected_class) + + +@pytest.mark.parametrize("name", [None, "None", "none", "NONE"]) +def test_configure_turbulence_model_none(name): + assert _configure_turbulence_model({"name": name, "c1": 1.0, "c2": 1.0}) is None + + +def test_configure_turbulence_model_unknown(): + with pytest.raises(NotImplementedError, match="UnknownTurb"): + _configure_turbulence_model({"name": "UnknownTurb", "c1": 1.0, "c2": 1.0}) + + +# --------------------------------------------------------------------------- +# Superposition model tests +# --------------------------------------------------------------------------- + + +@pytest.mark.parametrize( + "name,expected_class", + [ + ("Linear", LinearSum), + ("linear", LinearSum), + ("LINEAR", LinearSum), + ("Squared", SquaredSum), + ("squared", SquaredSum), + ("Max", MaxSum), + ("max", MaxSum), + ("Weighted", WeightedSum), + ("weighted", WeightedSum), + ("Cumulative", CumulativeWakeSum), + ("cumulative", CumulativeWakeSum), + ], +) +def test_configure_superposition_model(name, expected_class): + assert isinstance( + _configure_superposition_model({"ws_superposition": name}), expected_class + ) + + +def test_configure_superposition_model_product_not_implemented(): + with pytest.raises(NotImplementedError, match="Product"): + _configure_superposition_model({"ws_superposition": "Product"}) + + +def test_configure_superposition_model_unknown(): + with pytest.raises(NotImplementedError, match="UnknownSuper"): + _configure_superposition_model({"ws_superposition": "UnknownSuper"}) + + +# --------------------------------------------------------------------------- +# Rotor averaging tests +# --------------------------------------------------------------------------- + + +@pytest.mark.parametrize( + "name,expected_class", + [ + ("Center", RotorCenter), + ("center", RotorCenter), + ("CENTER", RotorCenter), + ("avg_deficit", GridRotorAvg), + ("Avg_Deficit", GridRotorAvg), + ("EqGrid", EqGridRotorAvg), + ("eqgrid", EqGridRotorAvg), + ("GQGrid", GQGridRotorAvg), + ("gqgrid", GQGridRotorAvg), + ("PolarGrid", PolarGridRotorAvg), + ("polargrid", PolarGridRotorAvg), + ("CGI", CGIRotorAvg), + ("cgi", CGIRotorAvg), + ], +) +def test_configure_rotor_averaging(name, expected_class): + assert isinstance(_configure_rotor_averaging({"name": name}), expected_class) + + +def test_configure_rotor_averaging_eqgrid_n(): + assert isinstance( + _configure_rotor_averaging({"name": "EqGrid", "n": 9}), EqGridRotorAvg + ) + + +def test_configure_rotor_averaging_gqgrid_params(): + data = {"name": "GQGrid", "n_x_grid_points": 3, "n_y_grid_points": 5} + model = _configure_rotor_averaging(data) + assert isinstance(model, GQGridRotorAvg) + # Verify custom params produced different nodes than defaults + default = GQGridRotorAvg() + assert len(model.nodes_x) != len(default.nodes_x) + + +def test_configure_rotor_averaging_cgi_n(): + assert isinstance(_configure_rotor_averaging({"name": "CGI", "n": 7}), CGIRotorAvg) + + +def test_configure_rotor_averaging_unknown(): + with pytest.raises(NotImplementedError, match="UnknownRotor"): + _configure_rotor_averaging({"name": "UnknownRotor"}) + + +# --------------------------------------------------------------------------- +# Blockage model tests +# --------------------------------------------------------------------------- + + +@pytest.mark.parametrize( + "name,expected_class", + [ + ("SelfSimilarityDeficit2020", SelfSimilarityDeficit2020), + ("selfsimilaritydeficit2020", SelfSimilarityDeficit2020), + ("SelfSimilarityDeficit", SelfSimilarityDeficit), + ("selfsimilaritydeficit", SelfSimilarityDeficit), + ("RankineHalfBody", RankineHalfBody), + ("rankinehalfbody", RankineHalfBody), + ("Rathmann", Rathmann), + ("rathmann", Rathmann), + ("VortexCylinder", VortexCylinder), + ("vortexcylinder", VortexCylinder), + ("VortexDipole", VortexDipole), + ("vortexdipole", VortexDipole), + ("HybridInduction", HybridInduction), + ("hybridinduction", HybridInduction), + ], +) +def test_configure_blockage_model(name, expected_class): + model = _configure_blockage_model({"name": name, "ss_alpha": 0.888}, {}) + assert isinstance(model, expected_class) + + +@pytest.mark.parametrize("name", [None, "None", "none", "NONE"]) +def test_configure_blockage_model_none(name): + assert _configure_blockage_model({"name": name}, {}) is None + + +def test_configure_blockage_model_unknown(): + with pytest.raises(NotImplementedError, match="UnknownBlockage"): + _configure_blockage_model({"name": "UnknownBlockage"}, {}) + + +# --------------------------------------------------------------------------- +# get_with_default preserves extra user keys +# --------------------------------------------------------------------------- + + +def test_get_with_default_preserves_extra_keys(): + """Verify that get_with_default merges defaults without dropping user keys.""" + analysis = { + "rotor_averaging": { + "name": "GQGrid", + "n_x_grid_points": 3, + "n_y_grid_points": 5, + }, + } + result = get_with_default(analysis, "rotor_averaging", DEFAULTS) + assert result["name"] == "GQGrid" + assert result["n_x_grid_points"] == 3 + assert result["n_y_grid_points"] == 5 + + +def test_get_with_default_rotor_avg_eqgrid_n(): + """Verify EqGrid 'n' param survives through get_with_default.""" + analysis = {"rotor_averaging": {"name": "EqGrid", "n": 9}} + result = get_with_default(analysis, "rotor_averaging", DEFAULTS) + model = _configure_rotor_averaging(result) + assert isinstance(model, EqGridRotorAvg) + assert result["n"] == 9 + + +def test_get_with_default_fills_missing_keys(): + """Verify that missing keys are filled from defaults.""" + # deflection_model defaults have beta=0.1; user only provides name + analysis = {"deflection_model": {"name": "Jimenez"}} + result = get_with_default(analysis, "deflection_model", DEFAULTS) + assert result["name"] == "Jimenez" + assert result["beta"] == 0.1 + + +def test_get_with_default_recursive_nested_dicts(): + """Verify recursive merge fills deep missing keys while preserving user extras.""" + nested_defaults = { + "model": { + "params": {"a": 1, "b": 2}, + "name": "default", + } + } + # User provides partial nested dict (missing key "b") plus an extra key "c" + data = { + "model": { + "params": {"a": 10, "c": 99}, + "name": "custom", + } + } + result = get_with_default(data, "model", nested_defaults) + assert result["name"] == "custom" + assert result["params"]["a"] == 10 # user value preserved + assert result["params"]["b"] == 2 # missing key filled from defaults + assert result["params"]["c"] == 99 # extra user key preserved + + +# --------------------------------------------------------------------------- +# configure_wake_model return contract +# --------------------------------------------------------------------------- + + +def test_configure_wake_model_returns_wake_deficit_key(): + """Verify configure_wake_model returns wake_deficit_key for API compat.""" + system_dat = { + "attributes": { + "analysis": { + "wind_deficit_model": {"name": "Jensen"}, + } + } + } + config = configure_wake_model(system_dat, rotor_diameter=126.0, hub_height=90.0) + assert "wake_deficit_key" in config + assert config["wake_deficit_key"] is None diff --git a/wifa/pywake_api.py b/wifa/pywake_api.py index 670ab1b..5e473b0 100644 --- a/wifa/pywake_api.py +++ b/wifa/pywake_api.py @@ -1,19 +1,24 @@ import argparse -import os import warnings from pathlib import Path import numpy as np import xarray as xr -import yaml from scipy.interpolate import interp1d from scipy.special import gamma -from windIO import dict_to_netcdf, load_yaml -from windIO import validate as validate_yaml + +from wifa._optional import require from wifa._optional import require # Define default values for wind_deficit_model parameters + + +def _normalize_name(name): + """Normalize model name for case-insensitive matching.""" + return name.strip().lower().replace("-", "").replace("_", "") + + DEFAULTS = { "wind_deficit_model": { "name": "Jensen", @@ -43,16 +48,25 @@ def get_with_default(data, key, defaults): If the value is a dictionary, apply the same process recursively. """ if key not in data: - print("WARNING: Using default value for ", key) + warnings.warn(f"Using default value for {key}") return defaults[key] - elif isinstance(data[key], dict): - # For nested dictionaries, ensure all subkeys are checked for defaults - return { - sub_key: get_with_default(data[key], sub_key, defaults[key]) - for sub_key in defaults[key] - } - else: - return data[key] + + if isinstance(data[key], dict): + # Merge defaults into user dict: fill missing keys from defaults, + # but preserve all extra user-provided keys (e.g. n, n_x_grid_points). + # Recurse when both user and default values are dicts. + merged = dict(data[key]) + for sub_key in defaults[key]: + if sub_key not in merged: + warnings.warn(f"Using default value for {sub_key}") + merged[sub_key] = defaults[key][sub_key] + elif isinstance(merged[sub_key], dict) and isinstance( + defaults[key][sub_key], dict + ): + merged[sub_key] = get_with_default(data[key], sub_key, defaults[key]) + return merged + + return data[key] def load_and_validate_config(yaml_input, default_output_dir="output"): @@ -109,9 +123,7 @@ def create_turbines(farm_dat): turbine_dats = [farm_dat["turbines"]] type_names = "0" else: - turbine_dats = [ - farm_dat["turbine_types"][key] for key in farm_dat["turbine_types"] - ] + turbine_dats = list(farm_dat["turbine_types"].values()) type_names = list(farm_dat["turbine_types"].keys()) turbines = [] @@ -227,7 +239,6 @@ def dict_to_site(resource_dict): # This is required for XRSite's linear interpolation, which expects the turbine index # as the leading dimension. resource_ds = resource_ds.transpose("i", *other_dims) - print("making site with ", resource_ds) return XRSite(resource_ds) @@ -263,68 +274,49 @@ def construct_site(system_dat, resource_dat, hub_heights, x_positions): dict with keys: site, ws, wd, TI, timeseries, operating, additional_heights, cases_idx, flow_bounds """ - from py_wake.examples.data.hornsrev1 import Hornsrev1Site - from py_wake.site import XRSite - from windIO import dict_to_netcdf - - # Get flow field bounds from config or site boundaries - boundaries = system_dat["site"]["boundaries"]["polygons"][0] - WFXLB = np.min(boundaries["x"]) - WFXUB = np.max(boundaries["x"]) - WFYLB = np.min(boundaries["y"]) - WFYUB = np.max(boundaries["y"]) - - # Override with explicit flow field bounds if specified - WFXLB = get_flow_field_param(system_dat, "xlb", WFXLB) - WFXUB = get_flow_field_param(system_dat, "xub", WFXUB) - WFYLB = get_flow_field_param(system_dat, "ylb", WFYLB) - WFYUB = get_flow_field_param(system_dat, "yub", WFYUB) - WFDX = get_flow_field_param(system_dat, "dx", (WFXUB - WFXLB) / 100) - WFDY = get_flow_field_param(system_dat, "dy", (WFYUB - WFYLB) / 100) - - flow_bounds = { - "xlb": WFXLB, - "xub": WFXUB, - "ylb": WFYLB, - "yub": WFYUB, - "dx": WFDX, - "dy": WFDY, - } + # Compute flow field bounds from site boundaries, with optional overrides + flow_bounds = _compute_flow_bounds(system_dat) # Determine site type and construct accordingly - if "time" in resource_dat["wind_resource"]: - # Timeseries site + wind_resource = resource_dat["wind_resource"] + if "time" in wind_resource: result = _construct_timeseries_site( system_dat, resource_dat, hub_heights, x_positions ) - result["flow_bounds"] = flow_bounds - return result - - elif "weibull_k" in resource_dat["wind_resource"]: - # Weibull distribution site + elif "weibull_k" in wind_resource: result = _construct_weibull_site(resource_dat, hub_heights, x_positions) - result["flow_bounds"] = flow_bounds - return result - else: - # Simple probability-based site - ws = resource_dat["wind_resource"]["wind_speed"] - wd = resource_dat["wind_resource"]["wind_direction"] - site = dict_to_site(resource_dat["wind_resource"]) - TI = resource_dat["wind_resource"]["turbulence_intensity"]["data"] - - return { - "site": site, - "ws": ws, - "wd": wd, - "TI": TI, + result = { + "site": dict_to_site(wind_resource), + "ws": wind_resource["wind_speed"], + "wd": wind_resource["wind_direction"], + "TI": wind_resource["turbulence_intensity"]["data"], "timeseries": False, "operating": np.ones((len(x_positions), 1)), "additional_heights": [], "cases_idx": np.ones(1).astype(bool), - "flow_bounds": flow_bounds, } + result["flow_bounds"] = flow_bounds + return result + + +def _compute_flow_bounds(system_dat): + """Compute flow field bounds from site boundaries with optional overrides.""" + boundaries = system_dat["site"]["boundaries"]["polygons"][0] + xlb = get_flow_field_param(system_dat, "xlb", np.min(boundaries["x"])) + xub = get_flow_field_param(system_dat, "xub", np.max(boundaries["x"])) + ylb = get_flow_field_param(system_dat, "ylb", np.min(boundaries["y"])) + yub = get_flow_field_param(system_dat, "yub", np.max(boundaries["y"])) + return { + "xlb": xlb, + "xub": xub, + "ylb": ylb, + "yub": yub, + "dx": get_flow_field_param(system_dat, "dx", (xub - xlb) / 100), + "dy": get_flow_field_param(system_dat, "dy", (yub - ylb) / 100), + } + def _construct_timeseries_site(system_dat, resource_dat, hub_heights, x_positions): """Construct site from timeseries data. @@ -339,14 +331,14 @@ def _construct_timeseries_site(system_dat, resource_dat, hub_heights, x_position cases_idx = np.ones(len(times)).astype(bool) # Check for subset configuration - output_spec = system_dat["attributes"].get("model_outputs_specification", {}) - if "run_configuration" in output_spec: - run_config = output_spec["run_configuration"] - if "times_run" in run_config and not run_config["times_run"].get( - "all_occurences", True - ): - if "subset" in run_config["times_run"]: - cases_idx = run_config["times_run"]["subset"] + times_run = ( + system_dat["attributes"] + .get("model_outputs_specification", {}) + .get("run_configuration", {}) + .get("times_run", {}) + ) + if not times_run.get("all_occurences", True) and "subset" in times_run: + cases_idx = times_run["subset"] heights = wind_resource.get("height") @@ -451,7 +443,6 @@ def get_resource_data(var_name): ) else: # Single turbine type - print(np.array(ws).shape, np.array(heights).shape) if heights: ws, wd = _interpolate_wind_data(heights, ws, wd, hh) @@ -598,27 +589,11 @@ def configure_wake_model(system_dat, rotor_diameter, hub_height): turbulence_model, superposition_model, rotor_averaging, blockage_model, solver_class, solver_args """ - from py_wake.deficit_models import SelfSimilarityDeficit2020 - from py_wake.deficit_models.fuga import FugaDeficit - from py_wake.deficit_models.gaussian import ( - BastankhahGaussianDeficit, - BlondelSuperGaussianDeficit2020, - TurboGaussianDeficit, - ) - from py_wake.deficit_models.noj import NOJLocalDeficit - from py_wake.deflection_models import JimenezWakeDeflection - from py_wake.rotor_avg_models import GridRotorAvg, RotorCenter - from py_wake.superposition_models import LinearSum, SquaredSum - from py_wake.turbulence_models import ( - CrespoHernandez, - STF2005TurbulenceModel, - STF2017TurbulenceModel, - ) from py_wake.wind_farm_models import All2AllIterative, PropagateDownwind analysis = system_dat["attributes"]["analysis"] - # Get model configurations with defaults + # Resolve each submodel config, filling missing keys from DEFAULTS wind_deficit_data = get_with_default(analysis, "wind_deficit_model", DEFAULTS) deflection_data = get_with_default(analysis, "deflection_model", DEFAULTS) turbulence_data = get_with_default(analysis, "turbulence_model", DEFAULTS) @@ -626,35 +601,16 @@ def configure_wake_model(system_dat, rotor_diameter, hub_height): rotor_avg_data = get_with_default(analysis, "rotor_averaging", DEFAULTS) blockage_data = get_with_default(analysis, "blockage_model", DEFAULTS) - # Configure wind deficit model - deficit_args = {"use_effective_ws": True} - wake_deficit_key = None - - print("Running deficit ", wind_deficit_data) - - wake_model_class, deficit_args, wake_deficit_key = _configure_deficit_model( - wind_deficit_data, analysis, rotor_diameter, hub_height, deficit_args + wake_model_class, deficit_args = _configure_deficit_model( + wind_deficit_data, analysis, rotor_diameter, hub_height ) - - print("deficit args ", deficit_args) - - # Configure deflection model deflection_model = _configure_deflection_model(deflection_data) - - # Configure turbulence model turbulence_model = _configure_turbulence_model(turbulence_data) - - # Configure superposition model superposition_model = _configure_superposition_model(superposition_data) - print("using superposition ", superposition_data) - - # Configure rotor averaging rotor_averaging = _configure_rotor_averaging(rotor_avg_data) - - # Configure blockage model blockage_model = _configure_blockage_model(blockage_data, deficit_args) - # Determine solver based on blockage + # Blockage requires All2AllIterative solver solver_args = {} if blockage_model is not None: solver_class = All2AllIterative @@ -665,7 +621,7 @@ def configure_wake_model(system_dat, rotor_diameter, hub_height): return { "wake_model_class": wake_model_class, "deficit_args": deficit_args, - "wake_deficit_key": wake_deficit_key, + "wake_deficit_key": None, # Deprecated: kept for API compatibility "deflection_model": deflection_model, "turbulence_model": turbulence_model, "superposition_model": superposition_model, @@ -676,58 +632,105 @@ def configure_wake_model(system_dat, rotor_diameter, hub_height): } -def _configure_deficit_model( - wind_deficit_data, analysis, rotor_diameter, hub_height, deficit_args -): +def _configure_deficit_model(wind_deficit_data, analysis, rotor_diameter, hub_height): """Configure the wind deficit model. Returns: - tuple: (wake_model_class, deficit_args, wake_deficit_key) + tuple: (wake_model_class, deficit_args) """ from py_wake.deficit_models.fuga import FugaDeficit from py_wake.deficit_models.gaussian import ( BastankhahGaussianDeficit, BlondelSuperGaussianDeficit2020, + BlondelSuperGaussianDeficit2023, + CarbajofuertesGaussianDeficit, + NiayifarGaussianDeficit, TurboGaussianDeficit, + ZongGaussianDeficit, ) - from py_wake.deficit_models.noj import NOJLocalDeficit + from py_wake.deficit_models.gcl import GCLDeficit + from py_wake.deficit_models.noj import NOJLocalDeficit, TurboNOJDeficit - wake_deficit_key = None model_name = wind_deficit_data["name"] + normalized = _normalize_name(model_name) + deficit_args = {"use_effective_ws": True} - if model_name == "Jensen": + wind_deficit_cfg = analysis.get("wind_deficit_model", {}) + wake_expansion = wind_deficit_cfg.get("wake_expansion_coefficient", {}) + + GAUSSIAN_MODELS = { + "bastankhah2014": BastankhahGaussianDeficit, + "niayifar2016": NiayifarGaussianDeficit, + "zong2020": ZongGaussianDeficit, + "carbajofuertes2018": CarbajofuertesGaussianDeficit, + } + # Models that accept a=[k_a, k_b] instead of k (scalar) + A_PARAM_MODELS = {"niayifar2016", "zong2020", "carbajofuertes2018"} + + if normalized in ("jensen", "nojlocaldeficit"): wake_model_class = NOJLocalDeficit - wake_expansion = analysis.get("wind_deficit_model", {}).get( - "wake_expansion_coefficient", {} - ) - if "k_b" in wake_expansion: - k_a = wake_expansion.get("k_a", 0) - k_b = wake_expansion["k_b"] - deficit_args["a"] = [k_a, k_b] - - elif model_name.lower() == "bastankhah2014": - wake_model_class = BastankhahGaussianDeficit - wake_expansion = analysis.get("wind_deficit_model", {}).get( - "wake_expansion_coefficient", {} - ) if "k_b" in wake_expansion: - deficit_args["k"] = wake_expansion["k_b"] - elif "k" in wake_expansion: - deficit_args["k"] = wake_expansion["k"] - if "ceps" in analysis.get("wind_deficit_model", {}): - deficit_args["ceps"] = analysis["wind_deficit_model"]["ceps"] - - elif model_name == "SuperGaussian": + deficit_args["a"] = [wake_expansion.get("k_a", 0), wake_expansion["k_b"]] + + elif normalized in GAUSSIAN_MODELS: + wake_model_class = GAUSSIAN_MODELS[normalized] + if normalized in A_PARAM_MODELS: + # Niayifar, Zong, Carbajofuertes use a=[k_a, k_b] + if "k" in wake_expansion: + warnings.warn( + f"{model_name} uses a=[k_a, k_b] for wake expansion, not scalar k. " + f"Scalar 'k' is ignored; specify k_a/k_b instead." + ) + if "k_b" in wake_expansion: + if "k_a" not in wake_expansion: + warnings.warn( + f"k_a not specified for {model_name}, defaulting to 0" + ) + deficit_args["a"] = [ + wake_expansion.get("k_a", 0), + wake_expansion["k_b"], + ] + else: + # Bastankhah2014 uses k (scalar) + if "k_b" in wake_expansion: + deficit_args["k"] = wake_expansion["k_b"] + elif "k" in wake_expansion: + deficit_args["k"] = wake_expansion["k"] + # ceps: valid for Bastankhah, Niayifar, Carbajofuertes (not Zong) + if normalized != "zong2020" and "ceps" in wind_deficit_cfg: + deficit_args["ceps"] = wind_deficit_cfg["ceps"] + # use_effective_ti: valid for Niayifar, Zong, Carbajofuertes (not Bastankhah) + if normalized != "bastankhah2014" and "use_effective_ti" in wind_deficit_cfg: + deficit_args["use_effective_ti"] = wind_deficit_cfg["use_effective_ti"] + + elif normalized == "supergaussian": wake_model_class = BlondelSuperGaussianDeficit2020 - elif model_name == "TurboPark": + elif normalized == "supergaussian2023": + wake_model_class = BlondelSuperGaussianDeficit2023 + + elif normalized == "turbopark": wake_model_class = TurboGaussianDeficit - elif model_name.upper() == "FUGA": + elif normalized == "turbonoj": + wake_model_class = TurboNOJDeficit + if "A" in wind_deficit_cfg: + deficit_args["A"] = wind_deficit_cfg["A"] + + elif normalized == "gcl": + wake_model_class = GCLDeficit + + elif normalized == "bastankhah2016": + raise NotImplementedError( + "Bastankhah2016 is not available in PyWake. Use flow_model 'foxes', " + "or choose Bastankhah2014/Zong2020 for PyWake." + ) + + elif normalized == "fuga": wake_model_class = FugaDeficit from pyfuga import get_luts - lut = get_luts( + get_luts( folder="luts", zeta0=0, nkz0=8, @@ -751,28 +754,31 @@ def _configure_deficit_model( else: raise NotImplementedError(f"Wake model '{model_name}' is not supported") - # Handle k/k2 format conversion - if "k2" in deficit_args: - k = deficit_args.pop("k") - k2 = deficit_args.pop("k2") - deficit_args["a"] = [k2, k] - - return wake_model_class, deficit_args, wake_deficit_key + return wake_model_class, deficit_args def _configure_deflection_model(deflection_data): """Configure the wake deflection model.""" from py_wake.deflection_models import JimenezWakeDeflection + from py_wake.deflection_models.gcl_hill_vortex import GCLHillDeflection + + name = deflection_data.get("name") + if name is None: + return None - name = deflection_data["name"].lower() - if name == "none": + normalized = _normalize_name(name) + if normalized == "none": return None - elif name == "jimenez": + if normalized == "jimenez": return JimenezWakeDeflection(beta=deflection_data["beta"]) - else: + if normalized == "gclhill": + return GCLHillDeflection() + if normalized == "bastankhah2016": raise NotImplementedError( - f"Deflection model '{deflection_data['name']}' is not supported" + "Bastankhah2016 deflection is not available in PyWake. Use flow_model " + "'foxes', or choose Jimenez/GCLHill for PyWake." ) + raise NotImplementedError(f"Deflection model '{name}' is not supported") def _configure_turbulence_model(turbulence_data): @@ -782,67 +788,132 @@ def _configure_turbulence_model(turbulence_data): STF2005TurbulenceModel, STF2017TurbulenceModel, ) + from py_wake.turbulence_models.gcl_turb import GCLTurbulence + + name = turbulence_data.get("name") + if name is None: + return None - name = turbulence_data["name"].upper() - if turbulence_data["name"].lower() == "none": + normalized = _normalize_name(name) + if normalized == "none": return None - elif name == "STF2005": - return STF2005TurbulenceModel(c=[turbulence_data["c1"], turbulence_data["c2"]]) - elif name == "STF2017": - return STF2017TurbulenceModel(c=[turbulence_data["c1"], turbulence_data["c2"]]) - elif name == "CRESPOHERNANDEZ": + + STF_MODELS = { + "stf2005": STF2005TurbulenceModel, + "stf2017": STF2017TurbulenceModel, + "iecti2019": STF2017TurbulenceModel, + } + + if normalized in STF_MODELS: + c = [turbulence_data.get("c1", 1.0), turbulence_data.get("c2", 1.0)] + return STF_MODELS[normalized](c=c) + if normalized == "crespohernandez": return CrespoHernandez() - else: - raise NotImplementedError( - f"Turbulence model '{turbulence_data['name']}' is not supported" - ) + if normalized == "gcl": + return GCLTurbulence() + raise NotImplementedError(f"Turbulence model '{name}' is not supported") def _configure_superposition_model(superposition_data): """Configure the superposition model.""" - from py_wake.superposition_models import LinearSum, SquaredSum + from py_wake.superposition_models import ( + CumulativeWakeSum, + LinearSum, + MaxSum, + SquaredSum, + WeightedSum, + ) - name = superposition_data["ws_superposition"].lower() - if name == "linear": - return LinearSum() - elif name == "squared": - return SquaredSum() - else: - raise NotImplementedError( - f"Superposition model '{superposition_data['ws_superposition']}' is not supported" - ) + name = superposition_data["ws_superposition"] + normalized = _normalize_name(name) + + SUPERPOSITION_MODELS = { + "linear": LinearSum, + "squared": SquaredSum, + "max": MaxSum, + "weighted": WeightedSum, + "cumulative": CumulativeWakeSum, + } + + if normalized in SUPERPOSITION_MODELS: + return SUPERPOSITION_MODELS[normalized]() + if normalized == "product": + raise NotImplementedError("Product superposition is not available in PyWake.") + raise NotImplementedError(f"Superposition model '{name}' is not supported") def _configure_rotor_averaging(rotor_avg_data): """Configure the rotor averaging model.""" - from py_wake.rotor_avg_models import GridRotorAvg, RotorCenter + from py_wake.rotor_avg_models import ( + CGIRotorAvg, + EqGridRotorAvg, + GQGridRotorAvg, + GridRotorAvg, + PolarGridRotorAvg, + RotorCenter, + ) - name = rotor_avg_data["name"].lower() - if name == "center": - print("Using Center Average") + name = rotor_avg_data["name"] + normalized = _normalize_name(name) + + if normalized == "center": return RotorCenter() - elif name == "avg_deficit": + if normalized == "avgdeficit": return GridRotorAvg() - else: - raise NotImplementedError( - f"Rotor averaging model '{rotor_avg_data['name']}' is not supported" + if normalized == "eqgrid": + return EqGridRotorAvg(n=rotor_avg_data.get("n", 4)) + if normalized == "gqgrid": + return GQGridRotorAvg( + n_x=rotor_avg_data.get("n_x_grid_points", 4), + n_y=rotor_avg_data.get("n_y_grid_points", 4), ) + if normalized == "polargrid": + return PolarGridRotorAvg() + if normalized == "cgi": + return CGIRotorAvg(n=rotor_avg_data.get("n", 4)) + raise NotImplementedError(f"Rotor averaging model '{name}' is not supported") def _configure_blockage_model(blockage_data, deficit_args): """Configure the blockage model.""" - from py_wake.deficit_models import SelfSimilarityDeficit2020 + from py_wake.deficit_models import ( + HybridInduction, + RankineHalfBody, + SelfSimilarityDeficit, + SelfSimilarityDeficit2020, + VortexCylinder, + VortexDipole, + ) from py_wake.deficit_models.fuga import FugaDeficit + from py_wake.deficit_models.rathmann import Rathmann name = blockage_data["name"] - if name == "None" or name is None: + if name is None: return None - elif name == "SelfSimilarityDeficit2020": - return SelfSimilarityDeficit2020(ss_alpha=blockage_data["ss_alpha"]) - elif name.upper() == "FUGA": + + normalized = _normalize_name(name) + if normalized == "none": + return None + + # Models that take no constructor arguments + SIMPLE_BLOCKAGE_MODELS = { + "selfsimilaritydeficit": SelfSimilarityDeficit, + "rankinehalfbody": RankineHalfBody, + "rathmann": Rathmann, + "vortexcylinder": VortexCylinder, + "vortexdipole": VortexDipole, + "hybridinduction": HybridInduction, + } + + if normalized == "selfsimilaritydeficit2020": + return SelfSimilarityDeficit2020( + ss_alpha=blockage_data.get("ss_alpha", 0.8888888888888888) + ) + if normalized in SIMPLE_BLOCKAGE_MODELS: + return SIMPLE_BLOCKAGE_MODELS[normalized]() + if normalized == "fuga": return FugaDeficit(deficit_args["LUT_path"]) - else: - raise ValueError(f"Unknown blockage model: {name}") + raise NotImplementedError(f"Blockage model '{name}' is not supported") def run_simulation(site, turbine, wake_config, site_data, x, y, turbine_types): @@ -861,16 +932,12 @@ def run_simulation(site, turbine, wake_config, site_data, x, y, turbine_types): dict with keys: sim_res, aep, aep_per_turbine """ # Build deficit model - print("Running ", wake_config["wake_model_class"], wake_config["deficit_args"]) deficit_model = wake_config["wake_model_class"]( rotorAvgModel=wake_config["rotor_averaging"], groundModel=None, **wake_config["deficit_args"], ) - if wake_config["wake_deficit_key"]: - deficit_model.WS_key = wake_config["wake_deficit_key"] - # Build wind farm model wind_farm_model = wake_config["solver_class"]( site, @@ -901,20 +968,10 @@ def run_simulation(site, turbine, wake_config, site_data, x, y, turbine_types): # Run simulation sim_res = wind_farm_model(**sim_kwargs) - aep = sim_res.aep(normalize_probabilities=not site_data["timeseries"]).sum() - print("aep is ", aep, "GWh") - - # Calculate per-turbine AEP - if site_data["timeseries"]: - aep_per_turbine = ( - sim_res.aep(normalize_probabilities=True).sum(["time"]).to_numpy() - ) - else: - aep_per_turbine = ( - sim_res.aep(normalize_probabilities=True).sum(["ws", "wd"]).to_numpy() - ) - - print(sim_res) + is_timeseries = site_data["timeseries"] + aep = sim_res.aep(normalize_probabilities=not is_timeseries).sum() + sum_dims = ["time"] if is_timeseries else ["ws", "wd"] + aep_per_turbine = sim_res.aep(normalize_probabilities=True).sum(sum_dims).to_numpy() return {"sim_res": sim_res, "aep": aep, "aep_per_turbine": aep_per_turbine} @@ -934,9 +991,8 @@ def generate_outputs(sim_results, system_dat, site_data, hub_heights, output_dir """ sim_res = sim_results["sim_res"] flow_bounds = site_data["flow_bounds"] - - # Ensure output directory exists - os.makedirs(output_dir, exist_ok=True) + output_path = Path(output_dir) + output_path.mkdir(parents=True, exist_ok=True) # Write turbine outputs if requested output_spec = system_dat["attributes"].get("model_outputs_specification", {}) @@ -944,20 +1000,17 @@ def generate_outputs(sim_results, system_dat, site_data, hub_heights, output_dir sim_res_formatted = sim_res[["Power", "WS_eff"]].rename( {"Power": "power", "WS_eff": "effective_wind_speed", "wt": "turbine"} ) - turbine_nc_filename = str( - output_spec.get("turbine_outputs", {}).get( - "turbine_nc_filename", "PowerTable.nc" - ) + turbine_nc_filename = output_spec["turbine_outputs"].get( + "turbine_nc_filename", "PowerTable.nc" ) - turbine_nc_filepath = Path(output_dir) / turbine_nc_filename - sim_res_formatted.to_netcdf(turbine_nc_filepath) + sim_res_formatted.to_netcdf(output_path / turbine_nc_filename) # Flow field handling flow_map = _generate_flow_field( sim_res, system_dat, site_data, hub_heights, flow_bounds ) - if flow_map: + if flow_map is not None: flow_map = flow_map[["WS_eff", "TI_eff"]].rename( { "h": "z", @@ -965,7 +1018,7 @@ def generate_outputs(sim_results, system_dat, site_data, hub_heights, output_dir "TI_eff": "turbulence_intensity", } ) - flow_map.to_netcdf(Path(output_dir) / "FarmFlow.nc") + flow_map.to_netcdf(output_path / "FarmFlow.nc") # Write YAML output _write_yaml_output(output_dir) @@ -980,70 +1033,52 @@ def _generate_flow_field(sim_res, system_dat, site_data, hub_heights, flow_bound Flow map xarray or None """ output_spec = system_dat["attributes"].get("model_outputs_specification", {}) - timeseries = site_data["timeseries"] - - WFXLB, WFXUB = flow_bounds["xlb"], flow_bounds["xub"] - WFYLB, WFYUB = flow_bounds["ylb"], flow_bounds["yub"] - WFDX, WFDY = flow_bounds["dx"], flow_bounds["dy"] + if "flow_field" not in output_spec: + return None - flow_map = None + x_range = np.arange( + flow_bounds["xlb"], flow_bounds["xub"] + flow_bounds["dx"], flow_bounds["dx"] + ) + y_range = np.arange( + flow_bounds["ylb"], flow_bounds["yub"] + flow_bounds["dy"], flow_bounds["dy"] + ) - if "flow_field" in output_spec and not timeseries: + if not site_data["timeseries"]: flow_map = sim_res.flow_box( - x=np.arange(WFXLB, WFXUB + WFDX, WFDX), - y=np.arange(WFYLB, WFYUB + WFDY, WFDY), + x=x_range, + y=y_range, h=list(hub_heights.values()), ) - # Warn if user requests unsupported outputs requested_vars = output_spec["flow_field"].get("output_variables", []) - if any( - var not in ["velocity_u", "turbulence_intensity"] for var in requested_vars - ): + unsupported = {"velocity_u", "turbulence_intensity"} + if any(var not in unsupported for var in requested_vars): warnings.warn("PyWake can only output velocity_u and turbulence_intensity") + return flow_map - elif "flow_field" in output_spec and timeseries: - flow_field_spec = output_spec["flow_field"] - if flow_field_spec.get("report") is not False: - z_list = flow_field_spec.get("z_list", sorted(list(hub_heights.values()))) - flow_map = sim_res.flow_box( - x=np.arange(WFXLB, WFXUB + WFDX, WFDX), - y=np.arange(WFYLB, WFYUB + WFDY, WFDY), - h=z_list, - time=sim_res.time.values, - ) + # Timeseries flow field + flow_field_spec = output_spec["flow_field"] + if flow_field_spec.get("report") is False: + return None - return flow_map + z_list = flow_field_spec.get("z_list", sorted(hub_heights.values())) + return sim_res.flow_box( + x=x_range, + y=y_range, + h=z_list, + time=sim_res.time.values, + ) def _write_yaml_output(output_dir): """Write the output YAML file with include directives.""" - data = { - "wind_energy_system": "INCLUDE_YAML_PLACEHOLDER", - "power_table": "INCLUDE_POWER_TABLE_PLACEHOLDER", - "flow_field": "INCLUDE_FLOW_FIELD_PLACEHOLDER", - } - - output_yaml_name = Path(output_dir) / "output.yaml" - with open(output_yaml_name, "w") as file: - yaml.dump(data, file, default_flow_style=False, allow_unicode=True) - - # Replace placeholders with include directives - with open(output_yaml_name, "r") as file: - yaml_content = file.read() - - yaml_content = yaml_content.replace( - "INCLUDE_YAML_PLACEHOLDER", "!include recorded_inputs.yaml" + # Write directly with !include tags (avoids round-trip through yaml.dump) + content = ( + "flow_field: !include FarmFlow.nc\n" + "power_table: !include PowerTable.nc\n" + "wind_energy_system: !include recorded_inputs.yaml\n" ) - yaml_content = yaml_content.replace( - "INCLUDE_POWER_TABLE_PLACEHOLDER", "!include PowerTable.nc" - ) - yaml_content = yaml_content.replace( - "INCLUDE_FLOW_FIELD_PLACEHOLDER", "!include FarmFlow.nc" - ) - - with open(output_yaml_name, "w") as file: - file.write(yaml_content) + (Path(output_dir) / "output.yaml").write_text(content) def run_pywake(yaml_input, output_dir="output"): diff --git a/wifa/wayve_api.py b/wifa/wayve_api.py index f731d19..83920bb 100644 --- a/wifa/wayve_api.py +++ b/wifa/wayve_api.py @@ -177,7 +177,7 @@ def run_wayve(yamlFile, output_dir="output", debug_mode=False): for time_index, time in enumerate(times): if debug_mode: # Print timestep - print(f"time {time_index+1}/{len(times)}") + print(f"time {time_index + 1}/{len(times)}") try: # Set up ABL abl = flow_io_abl(resource_dat["wind_resource"], time_index, hh, h1) From f8b7533ef05591191d9216e12f90b78b598f16db Mon Sep 17 00:00:00 2001 From: Bjarke Tobias Olsen Date: Fri, 13 Mar 2026 13:15:13 +0100 Subject: [PATCH 02/10] WIP --- tests/test_pywake_submodels.py | 26 +++++++++++++++++++++++++- wifa/pywake_api.py | 8 +++++++- 2 files changed, 32 insertions(+), 2 deletions(-) diff --git a/tests/test_pywake_submodels.py b/tests/test_pywake_submodels.py index 8c61c27..718cf84 100644 --- a/tests/test_pywake_submodels.py +++ b/tests/test_pywake_submodels.py @@ -26,7 +26,7 @@ ZongGaussianDeficit, ) from py_wake.deficit_models.gcl import GCLDeficit -from py_wake.deficit_models.noj import NOJLocalDeficit, TurboNOJDeficit +from py_wake.deficit_models.noj import NOJDeficit, NOJLocalDeficit, TurboNOJDeficit from py_wake.deficit_models.rathmann import Rathmann from py_wake.deflection_models import JimenezWakeDeflection from py_wake.deflection_models.gcl_hill_vortex import GCLHillDeflection @@ -109,6 +109,12 @@ def _call_deficit(name, analysis_extra=None): ("NOJLocalDeficit", NOJLocalDeficit), ("nojlocaldeficit", NOJLocalDeficit), ("NOJLOCALDEFICIT", NOJLocalDeficit), + ("Jensen_1983", NOJDeficit), + ("jensen_1983", NOJDeficit), + ("JENSEN_1983", NOJDeficit), + ("NOJDeficit", NOJDeficit), + ("nojdeficit", NOJDeficit), + ("NOJDEFICIT", NOJDeficit), ], ) def test_configure_deficit_model(name, expected_class): @@ -130,6 +136,8 @@ def test_configure_deficit_model(name, expected_class): ("TurboNOJ", TurboNOJDeficit), ("GCL", GCLDeficit), ("NOJLocalDeficit", NOJLocalDeficit), + ("Jensen_1983", NOJDeficit), + ("NOJDeficit", NOJDeficit), ], ) def test_configure_deficit_model_instantiation(name, expected_class): @@ -177,6 +185,17 @@ def test_configure_deficit_model_nojlocaldeficit_k_b(): assert args["a"] == [0.38, 0.004] +def test_configure_deficit_model_jensen_1983_k(): + """Verify Jensen_1983 passes scalar k and does not pass use_effective_ws.""" + cls, args = _call_deficit( + "Jensen_1983", + {"wake_expansion_coefficient": {"k": 0.04}}, + ) + assert cls is NOJDeficit + assert args["k"] == 0.04 + assert "use_effective_ws" not in args + + def test_configure_deficit_model_gaussian_params_niayifar(): """Verify Gaussian params pass through for Niayifar2016.""" cls, args = _call_deficit( @@ -260,6 +279,11 @@ def test_configure_deficit_model_a_param_warns_on_missing_k_a(): {"wake_expansion_coefficient": {"k_a": 0.38, "k_b": 0.004}}, NOJLocalDeficit, ), + ( + "Jensen_1983", + {"wake_expansion_coefficient": {"k": 0.04}}, + NOJDeficit, + ), ], ) def test_configure_deficit_model_instantiation_with_params(name, extra, expected_class): diff --git a/wifa/pywake_api.py b/wifa/pywake_api.py index 5e473b0..2bd38e7 100644 --- a/wifa/pywake_api.py +++ b/wifa/pywake_api.py @@ -649,7 +649,7 @@ def _configure_deficit_model(wind_deficit_data, analysis, rotor_diameter, hub_he ZongGaussianDeficit, ) from py_wake.deficit_models.gcl import GCLDeficit - from py_wake.deficit_models.noj import NOJLocalDeficit, TurboNOJDeficit + from py_wake.deficit_models.noj import NOJDeficit, NOJLocalDeficit, TurboNOJDeficit model_name = wind_deficit_data["name"] normalized = _normalize_name(model_name) @@ -672,6 +672,12 @@ def _configure_deficit_model(wind_deficit_data, analysis, rotor_diameter, hub_he if "k_b" in wake_expansion: deficit_args["a"] = [wake_expansion.get("k_a", 0), wake_expansion["k_b"]] + elif normalized in ("jensen1983", "nojdeficit"): + wake_model_class = NOJDeficit + deficit_args.pop("use_effective_ws", None) + if "k" in wake_expansion: + deficit_args["k"] = wake_expansion["k"] + elif normalized in GAUSSIAN_MODELS: wake_model_class = GAUSSIAN_MODELS[normalized] if normalized in A_PARAM_MODELS: From 56cff6ee394a8f9a70e2e5388f376d6028b930eb Mon Sep 17 00:00:00 2001 From: Bjarke Tobias Olsen Date: Fri, 13 Mar 2026 13:33:25 +0100 Subject: [PATCH 03/10] WIP --- pyproject.toml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pyproject.toml b/pyproject.toml index 06b5f4c..dc5c211 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -38,8 +38,8 @@ classifiers = [ ] requires-python = ">=3.10,<3.13" dependencies = [ - "windIO @ git+https://github.com/EUFlow/windIO.git", - "xarray>=2022.0.0,<2025", + "windIO @ git+https://github.com/bjarketol/windIO.git@expand-pywake-submodel-schema", + "xarray", "scipy", "pyyaml", ] From 9865a5896b9685f6dddb145e9dedb120145d0c57 Mon Sep 17 00:00:00 2001 From: Bjarke Tobias Olsen Date: Fri, 13 Mar 2026 13:42:05 +0100 Subject: [PATCH 04/10] WIP --- wifa/main_api.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/wifa/main_api.py b/wifa/main_api.py index dd267a9..fa102cb 100644 --- a/wifa/main_api.py +++ b/wifa/main_api.py @@ -17,7 +17,7 @@ def run_api(yaml_input): # validate input - validate_yaml(yaml_input, windIO.__path__[0] + "/plant/wind_energy_system.yaml") + validate_yaml(yaml_input, windIO.__path__[0] + "/schemas/plant/wind_energy_system.yaml") # get number of turbines if isinstance(yaml_input, dict): From f053e5283fccc41dde962d6783f74a5b2a99860c Mon Sep 17 00:00:00 2001 From: Bjarke Tobias Olsen Date: Fri, 13 Mar 2026 13:58:24 +0100 Subject: [PATCH 05/10] WIP --- wifa/main_api.py | 2 +- wifa/pywake_api.py | 4 ++++ 2 files changed, 5 insertions(+), 1 deletion(-) diff --git a/wifa/main_api.py b/wifa/main_api.py index fa102cb..87bda99 100644 --- a/wifa/main_api.py +++ b/wifa/main_api.py @@ -17,7 +17,7 @@ def run_api(yaml_input): # validate input - validate_yaml(yaml_input, windIO.__path__[0] + "/schemas/plant/wind_energy_system.yaml") + validate_yaml(yaml_input, "plant/wind_energy_system") # get number of turbines if isinstance(yaml_input, dict): diff --git a/wifa/pywake_api.py b/wifa/pywake_api.py index 2bd38e7..91e9f4c 100644 --- a/wifa/pywake_api.py +++ b/wifa/pywake_api.py @@ -229,6 +229,10 @@ def dict_to_site(resource_dict): if name in resource_ds: resource_ds = resource_ds.rename({name: rename_map[name]}) + if "time" in resource_ds.dims: + # Convert time coordinate to integer indices for GridInterpolator compatibility + # (string or datetime time coords cannot be interpolated numerically) + resource_ds = resource_ds.assign_coords(time=np.arange(len(resource_ds.time))) if "P" not in resource_ds and "time" in resource_ds.dims: n_time = len(resource_ds.time) # Create uniform probability array (1/N) From 49b0cd8115c13ccd50fa21d438b8cf2a4bb95702 Mon Sep 17 00:00:00 2001 From: Bjarke Tobias Olsen Date: Fri, 13 Mar 2026 15:15:38 +0100 Subject: [PATCH 06/10] Eliminate redundant YAML+netCDF loads in run_api MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit run_api() was loading the full YAML+netCDF 4 times per simulation: validate_yaml (load 1, result discarded), load_yaml (load 2), then the downstream model runner called validate_yaml (load 3) and load_yaml (load 4) again internally. Now validate_yaml is called once (its return value is the loaded dict), and the dict is passed to model runners — which skip their own validate+load when they receive a dict. Reduces peak transient memory from ~800 MB to ~200 MB for time-series wind resources. Co-Authored-By: Claude Opus 4.6 --- wifa/main_api.py | 20 ++++++-------------- 1 file changed, 6 insertions(+), 14 deletions(-) diff --git a/wifa/main_api.py b/wifa/main_api.py index 87bda99..1f54552 100644 --- a/wifa/main_api.py +++ b/wifa/main_api.py @@ -3,7 +3,6 @@ import sys import windIO -from windIO import load_yaml from windIO import validate as validate_yaml from .cs_api.cs_modules.csLaunch.cs_run_function import run_code_saturne @@ -16,40 +15,33 @@ def run_api(yaml_input): - # validate input - validate_yaml(yaml_input, "plant/wind_energy_system") - - # get number of turbines if isinstance(yaml_input, dict): yaml_dat = yaml_input else: - yaml_dat = load_yaml(yaml_input) + yaml_dat = validate_yaml(yaml_input, "plant/wind_energy_system") model_name = yaml_dat["attributes"]["flow_model"]["name"] if model_name.lower() == "pywake": - pywake_aep = run_pywake(yaml_input) + run_pywake(yaml_dat) elif model_name.lower() == "foxes": - foxes_aep = run_foxes(yaml_input) + run_foxes(yaml_dat) elif model_name.lower() == "floris": - floris_aep = run_floris(yaml_input) + run_floris(yaml_dat) elif model_name.lower() == "wayve": - # Output directory - # yaml_input_no_ext = os.path.splitext(yaml_input)[0] # Remove the file extension - # output_dir_name = 'output_wayve' + yaml_input_no_ext.replace(os.sep, '_') # Replace directory separators output_dir_name = yaml_dat["attributes"]["model_outputs_specification"][ "output_folder" ] if not os.path.exists(output_dir_name): os.makedirs(output_dir_name) - run_wayve(yaml_input, output_dir_name) + run_wayve(yaml_dat, output_dir_name) elif model_name.lower() == "codesaturne": - run_code_saturne(yaml_input, test_mode=True) + run_code_saturne(yaml_dat, test_mode=True) else: print("Invalid Model") From 1a912a056928770b8b6aa549b1e2b6ee86ce30bc Mon Sep 17 00:00:00 2001 From: btol Date: Tue, 17 Mar 2026 11:02:53 +0100 Subject: [PATCH 07/10] Use DensityCompensation instead of DensityScale for PyWake air density correction Switch from the default DensityScale (scales power/CT after lookup) to DensityCompensation (corrects wind speed before power curve lookup) to match foxes' air density handling approach: ws *= (rho/rho_ref)^(1/3). Co-Authored-By: Claude Opus 4.6 (1M context) --- wifa/pywake_api.py | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) diff --git a/wifa/pywake_api.py b/wifa/pywake_api.py index 91e9f4c..fbd1b60 100644 --- a/wifa/pywake_api.py +++ b/wifa/pywake_api.py @@ -114,8 +114,10 @@ def create_turbines(farm_dat): """ from py_wake.wind_turbines import WindTurbine, WindTurbines from py_wake.wind_turbines.power_ct_functions import ( + DensityCompensation, PowerCtFunctionList, PowerCtTabular, + SimpleYawModel, ) # Handle single vs multiple turbine types @@ -162,11 +164,18 @@ def create_turbines(farm_dat): cutin = turbine_dat["performance"].get("cutin_wind_speed", 0) cutout = turbine_dat["performance"].get("cutout_wind_speed") + # Use DensityCompensation (wind speed correction before lookup) to + # match foxes' air density handling: ws *= (rho/rho_ref)^(1/3) + density_models = [SimpleYawModel(exp=2), DensityCompensation(1.225)] + this_turbine = WindTurbine( name=turbine_dat["name"], diameter=rd, hub_height=hh, - powerCtFunction=PowerCtTabular(speeds, powers, power_unit="W", ct=cts_int), + powerCtFunction=PowerCtTabular( + speeds, powers, power_unit="W", ct=cts_int, + additional_models=density_models, + ), ws_cutin=cutin, ws_cutout=cutout, ) From de6923c45ccbcadc9a928373c1bcef9604ba88ef Mon Sep 17 00:00:00 2001 From: btol Date: Wed, 25 Mar 2026 21:29:09 +0100 Subject: [PATCH 08/10] Fix Speedup axis ordering and add proper Weibull flow case sampling - Fix _construct_weibull_site() Speedup axis bug: use dim-name lookup instead of hardcoded axis=0, and set Speedup dims from input data's actual ordering. Previously, flow_model_chain's (wind_direction, wind_turbine) ordering caused Speedup to be silently ignored by PyWake, removing terrain-induced wind speed inhomogeneity and inflating wake losses from ~10% to ~39%. - Add automatic flow case computation when wind_speed is absent from the windIO dict: - WS range: 0 to Speedup-adjusted 99.9% Weibull CDF point, 0.5 m/s steps - WD sub-sectors: 5 per sector (matching pywasp), letting PyWake's probability partitioning handle sector probability distribution - Handle missing turbulence_intensity gracefully (default 0.06) - Add test_weibull_speedup_dim_ordering regression test verifying both (wind_direction, wind_turbine) and (wind_turbine, wind_direction) orderings produce identical AEP - Update test_heterogeneous_wind_rose_grid baseline to match new auto-computed ws range and wd sub-sectors Co-Authored-By: Claude Opus 4.6 (1M context) --- tests/test_pywake.py | 127 ++++++++++++++++++++++++++++++++++++++++++- wifa/pywake_api.py | 97 +++++++++++++++++++++++++++------ 2 files changed, 205 insertions(+), 19 deletions(-) diff --git a/tests/test_pywake.py b/tests/test_pywake.py index 2605cf5..ca04bdf 100644 --- a/tests/test_pywake.py +++ b/tests/test_pywake.py @@ -274,9 +274,27 @@ def test_heterogeneous_wind_rose_grid(): x = [0, 1248.1, 2496.2, 3744.3] y = [0, 0, 0, 0] + # Compute Speedup-adjusted ws range (same logic as WIFA) + A_vals = dat["Weibull_A"].values + k_vals = dat["Weibull_k"].values + ws_999 = A_vals * (-np.log(0.001)) ** (1.0 / k_vals) + min_su = np.min(speedup) + ws_max_ref = np.max(ws_999) / max(min_su, 0.1) + ws_range = np.arange(0, np.ceil(ws_max_ref) + 0.5, 0.5) + + # Compute sub-sector wd (same logic as WIFA) + wd_sectors = dat["wd"].values + if len(wd_sectors) > 1 and np.isclose(wd_sectors[-1], 360.0): + wd_sectors = wd_sectors[:-1] + n_sub = 5 + sw = 360.0 / len(wd_sectors) + ssw = sw / n_sub + offsets = np.linspace(-sw / 2 + ssw / 2, sw / 2 - ssw / 2, n_sub) + wd_fine = np.sort((wd_sectors[:, np.newaxis] + offsets[np.newaxis, :]).ravel() % 360) + # compute AEP with PyWake res_aep = ( - wfm(x, y, ws=np.arange(2, 30, 1), wd=dat["wd"]) + wfm(x, y, ws=ws_range, wd=wd_fine) .aep(normalize_probabilities=True) .sum() ) @@ -466,6 +484,113 @@ def test_pywake_dict_timeseries_per_turbine_with_density(tmp_path): assert aep_with != aep_without +def test_weibull_speedup_dim_ordering(tmp_path): + """Regression test: per-turbine Weibull Speedup with both dim orderings. + + flow_model_chain (via windkit) writes wind_resource.nc with dims + (wind_direction, wind_turbine), while WIFA's own test fixtures use + (wind_turbine, wind_direction). A bug in _construct_weibull_site() + previously hardcoded axis=0 for the Speedup normalisation, which only + worked for the turbine-first ordering. With direction-first data the + Speedup dims were silently swapped and PyWake ignored the variable, + removing all terrain-induced wind speed inhomogeneity from the wake + simulation and inflating wake losses from ~10 % to ~39 %. + + This test runs the same per-turbine Weibull case with BOTH dim + orderings and asserts identical AEP. + """ + from conftest import _ANALYSIS, _TURBINE + + n_wd = 4 + n_wt = 4 + wd_vals = [0.0, 90.0, 180.0, 270.0] + ws_vals = list(np.arange(4.0, 26.0, 1.0).tolist()) + + # Per-turbine, per-sector Weibull A — turbine 3 is windiest + # Shape: (wind_direction, wind_turbine) = (4, 4) + A_data = [ + [7.0, 8.0, 9.0, 10.0], # sector 0° + [6.5, 7.5, 8.5, 9.5], # sector 90° + [8.0, 9.0, 10.0, 11.0], # sector 180° + [6.0, 7.0, 8.0, 9.0], # sector 270° + ] + k_data = [[2.0] * n_wt] * n_wd + freq_data = [[1.0 / n_wd] * n_wt] * n_wd + ti_data = [[0.06] * n_wt] * n_wd + + common_site = { + "name": "Test site", + "boundaries": { + "polygons": [{"x": [-90, 5000, 5000, -90], "y": [90, 90, -90, -90]}] + }, + } + common_farm = { + "name": "Test farm", + "layouts": [ + {"coordinates": {"x": [0, 1248.1, 2496.2, 3744.3], "y": [0, 0, 0, 0]}} + ], + "turbines": _TURBINE, + } + common_attrs = { + "flow_model": {"name": "pywake"}, + "analysis": _ANALYSIS, + "model_outputs_specification": { + "turbine_outputs": { + "turbine_nc_filename": "PowerTable.nc", + "output_variables": ["power"], + }, + }, + } + + def _make_system(data, dims, name): + return { + "name": name, + "site": { + **common_site, + "energy_resource": { + "name": "Test resource", + "wind_resource": { + "wind_direction": wd_vals, + "wind_speed": ws_vals, + "wind_turbine": list(range(n_wt)), + "reference_height": 119.0, + "weibull_a": {"data": data["A"], "dims": dims}, + "weibull_k": {"data": data["k"], "dims": dims}, + "sector_probability": {"data": data["f"], "dims": dims}, + "turbulence_intensity": {"data": data["ti"], "dims": dims}, + }, + }, + }, + "wind_farm": common_farm, + "attributes": common_attrs, + } + + # --- 1. Direction-first ordering (flow_model_chain convention) -------- + wd_first = _make_system( + {"A": A_data, "k": k_data, "f": freq_data, "ti": ti_data}, + ["wind_direction", "wind_turbine"], + "Direction-first", + ) + aep_wd_first = run_pywake(wd_first, output_dir=str(tmp_path / "wd_first")) + assert np.isfinite(aep_wd_first) and aep_wd_first > 0 + + # --- 2. Turbine-first ordering (WIFA test-fixture convention) --------- + A_T = np.array(A_data).T.tolist() + k_T = np.array(k_data).T.tolist() + freq_T = np.array(freq_data).T.tolist() + ti_T = np.array(ti_data).T.tolist() + + wt_first = _make_system( + {"A": A_T, "k": k_T, "f": freq_T, "ti": ti_T}, + ["wind_turbine", "wind_direction"], + "Turbine-first", + ) + aep_wt_first = run_pywake(wt_first, output_dir=str(tmp_path / "wt_first")) + + # Both orderings must produce identical AEP + npt.assert_allclose(aep_wd_first, aep_wt_first, 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 fbd1b60..979d9b2 100644 --- a/wifa/pywake_api.py +++ b/wifa/pywake_api.py @@ -490,27 +490,40 @@ def get_resource_data(var_name): } -def _construct_weibull_site(resource_dat, hub_heights, x_positions): +def _construct_weibull_site(resource_dat, hub_heights, x_positions, n_subsector=5): """Construct site from Weibull distribution data. Internal helper for construct_site(). + + Parameters + ---------- + resource_dat : dict + Energy resource dictionary from windIO. + hub_heights : dict + Mapping of turbine type names to hub heights. + x_positions : list + Turbine x positions (for operating array sizing). + n_subsector : int + Number of sub-directions per wind direction sector. Higher values + smooth directional wake effects. Default 5 (matching pywasp). """ from windIO import dict_to_netcdf wind_resource = resource_dat["wind_resource"] A = wind_resource["weibull_a"] k = wind_resource["weibull_k"] - wd = wind_resource["wind_direction"] - ws = wind_resource.get("wind_speed", np.arange(2, 30, 1)) + wd_raw = wind_resource["wind_direction"] + # --- Speedup computation ------------------------------------------------ # Handle turbine-specific Weibull if "wind_turbine" in wind_resource["sector_probability"]["dims"]: mean_ws = np.array(A["data"]) * gamma(1 + 1.0 / np.array(k["data"])) - max_mean = np.max(mean_ws, axis=0) + wt_axis = list(A["dims"]).index("wind_turbine") + max_mean = np.max(mean_ws, axis=wt_axis, keepdims=True) Speedup = mean_ws / max_mean wind_resource["Speedup"] = { - "dims": ["wind_turbine", "wd"], - "data": Speedup, + "dims": list(A["dims"]), + "data": Speedup.tolist(), } # Handle spatial Weibull @@ -523,21 +536,69 @@ def _construct_weibull_site(resource_dat, hub_heights, x_positions): "data": Speedup, } + # --- Flow case computation ----------------------------------------------- + # When wind_speed is absent from the windIO dict, WIFA computes optimal + # flow cases: a Speedup-adjusted ws range and sub-sector wd values. + # When wind_speed IS present, the user has chosen explicit flow cases + # and both ws and wd are used as-is. + ws = wind_resource.get("wind_speed", None) + if ws is None: + # -- Wind speed range from Weibull + Speedup -------------------------- + A_arr = np.asarray(A["data"], dtype=float) + k_arr = np.asarray(k["data"], dtype=float) + # Weibull inverse CDF at 99.9 %: ws = A * (-ln(0.001))^(1/k) + ws_999 = A_arr * (-np.log(0.001)) ** (1.0 / k_arr) + ws_max_local = float(np.max(ws_999)) + # Extend for speed-downs so the reference WS grid covers every + # turbine's distribution after Speedup scaling + if "Speedup" in wind_resource: + min_speedup = float(np.min(wind_resource["Speedup"]["data"])) + ws_max_ref = ws_max_local / max(min_speedup, 0.1) + else: + ws_max_ref = ws_max_local + ws = np.arange(0, np.ceil(ws_max_ref) + 0.5, 0.5) + + # -- Wind direction sub-sectors --------------------------------------- + # Strip 360° wrap-around before computing sub-sectors + wd_sectors = np.asarray(wd_raw, dtype=float) + if len(wd_sectors) > 1 and np.isclose(wd_sectors[-1], 360.0): + wd_sectors = wd_sectors[:-1] + if n_subsector > 1 and len(wd_sectors) >= 4: + n_sectors = len(wd_sectors) + sector_width = 360.0 / n_sectors + subsector_width = sector_width / n_subsector + offsets = np.linspace( + -sector_width / 2 + subsector_width / 2, + sector_width / 2 - subsector_width / 2, + n_subsector, + ) + wd = np.sort( + (wd_sectors[:, np.newaxis] + offsets[np.newaxis, :]).ravel() % 360 + ) + else: + wd = wd_sectors + else: + # Explicit wind_speed provided: use original wd as-is + wd = wd_raw + + # --- Site and TI -------------------------------------------------------- site = dict_to_site(wind_resource) - # Handle TI - site_ds = dict_to_netcdf(wind_resource) - if "x" in site_ds.turbulence_intensity.dims: - interpolated_ti = site_ds.turbulence_intensity.interp( - x=x_positions, y=x_positions - ) - if "height" in interpolated_ti.dims: - interpolated_ti = interpolated_ti.interp(height=hub_heights["0"]) - TI = np.array( - [interpolated_ti.isel(x=i, y=i).values for i in range(len(x_positions))] - ) + if "turbulence_intensity" in wind_resource: + site_ds = dict_to_netcdf(wind_resource) + if "x" in site_ds.turbulence_intensity.dims: + interpolated_ti = site_ds.turbulence_intensity.interp( + x=x_positions, y=x_positions + ) + if "height" in interpolated_ti.dims: + interpolated_ti = interpolated_ti.interp(height=hub_heights["0"]) + TI = np.array( + [interpolated_ti.isel(x=i, y=i).values for i in range(len(x_positions))] + ) + else: + TI = wind_resource["turbulence_intensity"]["data"] else: - TI = wind_resource["turbulence_intensity"]["data"] + TI = 0.06 # default when TI is absent from wind resource return { "site": site, From 64a87c64f2f818247ea263570cb343e32c5bc314 Mon Sep 17 00:00:00 2001 From: Bjarke Tobias Olsen Date: Thu, 26 Mar 2026 15:35:41 +0100 Subject: [PATCH 09/10] Support PyWake >= 2.6: SimpleYawModel no longer accepts exp parameter MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit SimpleYawModel(exp=2) was removed in PyWake 2.6 — the exp=2 behavior is now the default. Use try/except to support both old and new versions. Co-Authored-By: Claude Opus 4.6 (1M context) --- wifa/pywake_api.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/wifa/pywake_api.py b/wifa/pywake_api.py index 979d9b2..724e135 100644 --- a/wifa/pywake_api.py +++ b/wifa/pywake_api.py @@ -166,7 +166,11 @@ def create_turbines(farm_dat): # Use DensityCompensation (wind speed correction before lookup) to # match foxes' air density handling: ws *= (rho/rho_ref)^(1/3) - density_models = [SimpleYawModel(exp=2), DensityCompensation(1.225)] + try: + yaw_model = SimpleYawModel(exp=2) # PyWake < 2.6 + except TypeError: + yaw_model = SimpleYawModel() # PyWake >= 2.6 + density_models = [yaw_model, DensityCompensation(1.225)] this_turbine = WindTurbine( name=turbine_dat["name"], From 5bb32bdd8426e6c2a3bd89c7033a21ce77d20734 Mon Sep 17 00:00:00 2001 From: btol Date: Fri, 27 Mar 2026 09:41:19 +0100 Subject: [PATCH 10/10] Fix pre-commit --- .../KUL_LES/wind_energy_system/analysis_US.yaml | 2 +- tests/test_pywake.py | 16 +++++++--------- wifa/pywake_api.py | 7 ++++--- 3 files changed, 12 insertions(+), 13 deletions(-) diff --git a/examples/cases/KUL_LES/wind_energy_system/analysis_US.yaml b/examples/cases/KUL_LES/wind_energy_system/analysis_US.yaml index dca8fcc..a6ca679 100644 --- a/examples/cases/KUL_LES/wind_energy_system/analysis_US.yaml +++ b/examples/cases/KUL_LES/wind_energy_system/analysis_US.yaml @@ -64,7 +64,7 @@ HPC_config: mesh_node_number: 2 mesh_ntasks_per_node: 48 mesh_wall_time_hours: 1 - run_partition: "" + mesh_partition: "" # wckey: "" diff --git a/tests/test_pywake.py b/tests/test_pywake.py index ca04bdf..21c0f88 100644 --- a/tests/test_pywake.py +++ b/tests/test_pywake.py @@ -290,14 +290,12 @@ def test_heterogeneous_wind_rose_grid(): sw = 360.0 / len(wd_sectors) ssw = sw / n_sub offsets = np.linspace(-sw / 2 + ssw / 2, sw / 2 - ssw / 2, n_sub) - wd_fine = np.sort((wd_sectors[:, np.newaxis] + offsets[np.newaxis, :]).ravel() % 360) + wd_fine = np.sort( + (wd_sectors[:, np.newaxis] + offsets[np.newaxis, :]).ravel() % 360 + ) # compute AEP with PyWake - res_aep = ( - wfm(x, y, ws=ws_range, wd=wd_fine) - .aep(normalize_probabilities=True) - .sum() - ) + res_aep = wfm(x, y, ws=ws_range, wd=wd_fine).aep(normalize_probabilities=True).sum() # compute AEP with API wifa_res = run_pywake( @@ -509,10 +507,10 @@ def test_weibull_speedup_dim_ordering(tmp_path): # Per-turbine, per-sector Weibull A — turbine 3 is windiest # Shape: (wind_direction, wind_turbine) = (4, 4) A_data = [ - [7.0, 8.0, 9.0, 10.0], # sector 0° - [6.5, 7.5, 8.5, 9.5], # sector 90° + [7.0, 8.0, 9.0, 10.0], # sector 0° + [6.5, 7.5, 8.5, 9.5], # sector 90° [8.0, 9.0, 10.0, 11.0], # sector 180° - [6.0, 7.0, 8.0, 9.0], # sector 270° + [6.0, 7.0, 8.0, 9.0], # sector 270° ] k_data = [[2.0] * n_wt] * n_wd freq_data = [[1.0 / n_wd] * n_wt] * n_wd diff --git a/wifa/pywake_api.py b/wifa/pywake_api.py index 724e135..5f9e56f 100644 --- a/wifa/pywake_api.py +++ b/wifa/pywake_api.py @@ -9,8 +9,6 @@ from wifa._optional import require -from wifa._optional import require - # Define default values for wind_deficit_model parameters @@ -177,7 +175,10 @@ def create_turbines(farm_dat): diameter=rd, hub_height=hh, powerCtFunction=PowerCtTabular( - speeds, powers, power_unit="W", ct=cts_int, + speeds, + powers, + power_unit="W", + ct=cts_int, additional_models=density_models, ), ws_cutin=cutin,