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/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..718cf84 --- /dev/null +++ b/tests/test_pywake_submodels.py @@ -0,0 +1,590 @@ +"""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 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 +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), + ("Jensen_1983", NOJDeficit), + ("jensen_1983", NOJDeficit), + ("JENSEN_1983", NOJDeficit), + ("NOJDeficit", NOJDeficit), + ("nojdeficit", NOJDeficit), + ("NOJDEFICIT", NOJDeficit), + ], +) +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), + ("Jensen_1983", NOJDeficit), + ("NOJDeficit", NOJDeficit), + ], +) +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_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( + "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, + ), + ( + "Jensen_1983", + {"wake_expansion_coefficient": {"k": 0.04}}, + NOJDeficit, + ), + ], +) +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..a44729d 100644 --- a/wifa/pywake_api.py +++ b/wifa/pywake_api.py @@ -1,19 +1,22 @@ 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 # 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 +46,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 +121,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 +237,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 +272,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 +329,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 +441,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 +587,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 +599,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 +619,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 +630,111 @@ 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 NOJDeficit, NOJLocalDeficit, TurboNOJDeficit - wake_deficit_key = None model_name = wind_deficit_data["name"] + normalized = _normalize_name(model_name) + deficit_args = {"use_effective_ws": True} + + 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 model_name == "Jensen": + 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["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"] - if "ceps" in analysis.get("wind_deficit_model", {}): - deficit_args["ceps"] = analysis["wind_deficit_model"]["ceps"] - elif model_name == "SuperGaussian": + 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 +758,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 +792,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 + + normalized = _normalize_name(name) + if normalized == "none": return None - elif name == "SelfSimilarityDeficit2020": - return SelfSimilarityDeficit2020(ss_alpha=blockage_data["ss_alpha"]) - elif name.upper() == "FUGA": + + # 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 +936,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 +972,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 +995,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 +1004,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 +1022,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 +1037,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)