diff --git a/examples/cases/AWAKEN/AWAKEN/wind_energy_system/analysis.yaml b/examples/cases/AWAKEN/AWAKEN/wind_energy_system/analysis.yaml index 4e9086d..ef2bb49 100644 --- a/examples/cases/AWAKEN/AWAKEN/wind_energy_system/analysis.yaml +++ b/examples/cases/AWAKEN/AWAKEN/wind_energy_system/analysis.yaml @@ -9,9 +9,9 @@ wake_tool: "foxes" #pywake and foxes wind_deficit_model: name: Bastankhah2014 - wake_expansion_coefficient: # k = ka*ti + kb - k_a: 0.0 - k_b: 0.04 + wake_expansion_coefficient: # k = k_a + k_b*ti + k_a: 0.04 + k_b: 0.0 free_stream_ti: false ceps: 0.2 use_effective_ws: true diff --git a/examples/cases/AWAKEN/wind_energy_system/analysis.yaml b/examples/cases/AWAKEN/wind_energy_system/analysis.yaml index 801f19a..5b9db18 100644 --- a/examples/cases/AWAKEN/wind_energy_system/analysis.yaml +++ b/examples/cases/AWAKEN/wind_energy_system/analysis.yaml @@ -21,9 +21,9 @@ wm_coupling: wind_deficit_model: #name: SuperGaussian name: Bastankhah2014 - wake_expansion_coefficient: # k = ka*ti + kb - k_a: 0.0 - k_b: 0.04 + wake_expansion_coefficient: # k = k_a + k_b*ti + k_a: 0.04 + k_b: 0.0 free_stream_ti: false ceps: 0.2 use_effective_ws: true diff --git a/examples/cases/KUL_LES/wind_energy_system/analysis_PB.yaml b/examples/cases/KUL_LES/wind_energy_system/analysis_PB.yaml index 8ada863..215b20f 100644 --- a/examples/cases/KUL_LES/wind_energy_system/analysis_PB.yaml +++ b/examples/cases/KUL_LES/wind_energy_system/analysis_PB.yaml @@ -19,9 +19,9 @@ wm_coupling: #pywake and foxes wind_deficit_model: name: Bastankhah2014 - wake_expansion_coefficient: # k = ka*ti + kb - k_a: 0.0 - k_b: 0.04 + wake_expansion_coefficient: # k = k_a + k_b*ti + k_a: 0.04 + k_b: 0.0 free_stream_ti: false ceps: 0.2 use_effective_ws: true 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..d575d16 100644 --- a/examples/cases/KUL_LES/wind_energy_system/analysis_US.yaml +++ b/examples/cases/KUL_LES/wind_energy_system/analysis_US.yaml @@ -27,9 +27,9 @@ wm_coupling: #pywake and foxes wind_deficit_model: name: Bastankhah2014 - wake_expansion_coefficient: # k = ka*ti + kb - k_a: 0.0 - k_b: 0.04 + wake_expansion_coefficient: # k = k_a + k_b*ti + k_a: 0.04 + k_b: 0.0 free_stream_ti: false ceps: 0.2 use_effective_ws: true @@ -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/examples/cases/KUL_LES/wind_energy_system/analysis_VM.yaml b/examples/cases/KUL_LES/wind_energy_system/analysis_VM.yaml index 0e1cace..8e8599a 100644 --- a/examples/cases/KUL_LES/wind_energy_system/analysis_VM.yaml +++ b/examples/cases/KUL_LES/wind_energy_system/analysis_VM.yaml @@ -27,9 +27,9 @@ wm_coupling: #pywake and foxes wind_deficit_model: name: Bastankhah2014 - wake_expansion_coefficient: # k = ka*ti + kb - k_a: 0.0 - k_b: 0.04 + wake_expansion_coefficient: # k = k_a + k_b*ti + k_a: 0.04 + k_b: 0.0 free_stream_ti: false ceps: 0.2 use_effective_ws: true diff --git a/examples/cases/KUL_LES/wind_energy_system/analysis_pywake.yaml b/examples/cases/KUL_LES/wind_energy_system/analysis_pywake.yaml index 27736b9..745e1d8 100644 --- a/examples/cases/KUL_LES/wind_energy_system/analysis_pywake.yaml +++ b/examples/cases/KUL_LES/wind_energy_system/analysis_pywake.yaml @@ -19,9 +19,9 @@ wm_coupling: #pywake and foxes wind_deficit_model: name: Bastankhah2014 - wake_expansion_coefficient: # k = ka*ti + kb - k_a: 0.0 - k_b: 0.04 + wake_expansion_coefficient: # k = k_a + k_b*ti + k_a: 0.04 + k_b: 0.0 free_stream_ti: false ceps: 0.2 use_effective_ws: true diff --git a/examples/cases/heterogeneous_wind_rose_at_turbines/wind_energy_system/analysis.yaml b/examples/cases/heterogeneous_wind_rose_at_turbines/wind_energy_system/analysis.yaml index 7f66d39..dd44c88 100644 --- a/examples/cases/heterogeneous_wind_rose_at_turbines/wind_energy_system/analysis.yaml +++ b/examples/cases/heterogeneous_wind_rose_at_turbines/wind_energy_system/analysis.yaml @@ -1,9 +1,9 @@ #pywake and foxes wind_deficit_model: name: Bastankhah2014 - wake_expansion_coefficient: # k = ka*ti + kb - k_a: 0.0 - k_b: 0.04 + wake_expansion_coefficient: # k = k_a + k_b*ti + k_a: 0.04 + k_b: 0.0 free_stream_ti: false ceps: 0.2 use_effective_ws: true diff --git a/examples/cases/heterogeneous_wind_rose_map/wind_energy_system/analysis.yaml b/examples/cases/heterogeneous_wind_rose_map/wind_energy_system/analysis.yaml index 7f66d39..dd44c88 100644 --- a/examples/cases/heterogeneous_wind_rose_map/wind_energy_system/analysis.yaml +++ b/examples/cases/heterogeneous_wind_rose_map/wind_energy_system/analysis.yaml @@ -1,9 +1,9 @@ #pywake and foxes wind_deficit_model: name: Bastankhah2014 - wake_expansion_coefficient: # k = ka*ti + kb - k_a: 0.0 - k_b: 0.04 + wake_expansion_coefficient: # k = k_a + k_b*ti + k_a: 0.04 + k_b: 0.0 free_stream_ti: false ceps: 0.2 use_effective_ws: true diff --git a/examples/cases/open_source_scada/wind_energy_system/analysis.yaml b/examples/cases/open_source_scada/wind_energy_system/analysis.yaml index 264184f..79592bf 100644 --- a/examples/cases/open_source_scada/wind_energy_system/analysis.yaml +++ b/examples/cases/open_source_scada/wind_energy_system/analysis.yaml @@ -1,9 +1,9 @@ #pywake and foxes wind_deficit_model: name: Bastankhah2014 - wake_expansion_coefficient: # k = ka*ti + kb - k_a: 0.00 - k_b: 0.04 + wake_expansion_coefficient: # k = k_a + k_b*ti + k_a: 0.04 + k_b: 0.0 free_stream_ti: false ceps: 0.2 use_effective_ws: true diff --git a/examples/cases/simple_wind_rose/wind_energy_system/analysis.yaml b/examples/cases/simple_wind_rose/wind_energy_system/analysis.yaml index 829c278..907bcc1 100644 --- a/examples/cases/simple_wind_rose/wind_energy_system/analysis.yaml +++ b/examples/cases/simple_wind_rose/wind_energy_system/analysis.yaml @@ -1,9 +1,9 @@ #pywake and foxes wind_deficit_model: name: Bastankhah2014 - wake_expansion_coefficient: # k = ka*ti + kb - k_a: 0.0 - k_b: 0.04 + wake_expansion_coefficient: # k = k_a + k_b*ti + k_a: 0.04 + k_b: 0.0 free_stream_ti: false ceps: 0.2 use_effective_ws: true diff --git a/examples/cases/timeseries_with_operating_flag/wind_energy_system/analysis.yaml b/examples/cases/timeseries_with_operating_flag/wind_energy_system/analysis.yaml index 1278689..8545947 100644 --- a/examples/cases/timeseries_with_operating_flag/wind_energy_system/analysis.yaml +++ b/examples/cases/timeseries_with_operating_flag/wind_energy_system/analysis.yaml @@ -1,9 +1,9 @@ #pywake and foxes wind_deficit_model: name: Bastankhah2014 - wake_expansion_coefficient: # k = ka*ti + kb - k_a: 0.0 - k_b: 0.04 + wake_expansion_coefficient: # k = k_a + k_b*ti + k_a: 0.04 + k_b: 0.0 free_stream_ti: false ceps: 0.2 use_effective_ws: true diff --git a/examples/cases/turbine_specific_speeds_timeseries/wind_energy_system/analysis.yaml b/examples/cases/turbine_specific_speeds_timeseries/wind_energy_system/analysis.yaml index 1278689..8545947 100644 --- a/examples/cases/turbine_specific_speeds_timeseries/wind_energy_system/analysis.yaml +++ b/examples/cases/turbine_specific_speeds_timeseries/wind_energy_system/analysis.yaml @@ -1,9 +1,9 @@ #pywake and foxes wind_deficit_model: name: Bastankhah2014 - wake_expansion_coefficient: # k = ka*ti + kb - k_a: 0.0 - k_b: 0.04 + wake_expansion_coefficient: # k = k_a + k_b*ti + k_a: 0.04 + k_b: 0.0 free_stream_ti: false ceps: 0.2 use_effective_ws: true diff --git a/examples/cases/windio_4turbines/wind_energy_system/analysis.yaml b/examples/cases/windio_4turbines/wind_energy_system/analysis.yaml index 234eb2a..85977ee 100644 --- a/examples/cases/windio_4turbines/wind_energy_system/analysis.yaml +++ b/examples/cases/windio_4turbines/wind_energy_system/analysis.yaml @@ -1,9 +1,9 @@ #pywake and foxes wind_deficit_model: name: Bastankhah2014 - wake_expansion_coefficient: # k = ka*ti + kb - k_a: 0.0 - k_b: 0.04 + wake_expansion_coefficient: # k = k_a + k_b*ti + k_a: 0.04 + k_b: 0.0 free_stream_ti: false ceps: 0.2 use_effective_ws: true diff --git a/examples/cases/windio_4turbines_ABL/wind_energy_system/analysis.yaml b/examples/cases/windio_4turbines_ABL/wind_energy_system/analysis.yaml index c510bf0..06bff47 100644 --- a/examples/cases/windio_4turbines_ABL/wind_energy_system/analysis.yaml +++ b/examples/cases/windio_4turbines_ABL/wind_energy_system/analysis.yaml @@ -19,9 +19,9 @@ wm_coupling: #pywake and foxes wind_deficit_model: name: Bastankhah2014 - wake_expansion_coefficient: # k = ka*ti + kb - k_a: 0.0 - k_b: 0.04 + wake_expansion_coefficient: # k = k_a + k_b*ti + k_a: 0.04 + k_b: 0.0 free_stream_ti: false ceps: 0.2 use_effective_ws: true diff --git a/examples/cases/windio_4turbines_ABL_stable/wind_energy_system/analysis.yaml b/examples/cases/windio_4turbines_ABL_stable/wind_energy_system/analysis.yaml index 6344b3f..8a53136 100644 --- a/examples/cases/windio_4turbines_ABL_stable/wind_energy_system/analysis.yaml +++ b/examples/cases/windio_4turbines_ABL_stable/wind_energy_system/analysis.yaml @@ -17,9 +17,9 @@ apm_grid: #pywake and foxes wind_deficit_model: name: Bastankhah2014 - wake_expansion_coefficient: # k = ka*ti + kb - k_a: 0.0 - k_b: 0.04 + wake_expansion_coefficient: # k = k_a + k_b*ti + k_a: 0.04 + k_b: 0.0 free_stream_ti: false ceps: 0.2 use_effective_ws: true diff --git a/examples/cases/windio_4turbines_multipleTurbines/wind_energy_system/analysis.yaml b/examples/cases/windio_4turbines_multipleTurbines/wind_energy_system/analysis.yaml index 9278f22..b2d3cf1 100644 --- a/examples/cases/windio_4turbines_multipleTurbines/wind_energy_system/analysis.yaml +++ b/examples/cases/windio_4turbines_multipleTurbines/wind_energy_system/analysis.yaml @@ -1,9 +1,9 @@ #pywake and foxes wind_deficit_model: name: Bastankhah2014 - wake_expansion_coefficient: # k = ka*ti + kb - k_a: 0.0 - k_b: 0.04 + wake_expansion_coefficient: # k = k_a + k_b*ti + k_a: 0.04 + k_b: 0.0 free_stream_ti: false ceps: 0.2 use_effective_ws: true diff --git a/examples/cases/windio_4turbines_profiles_stable/wind_energy_system/analysis.yaml b/examples/cases/windio_4turbines_profiles_stable/wind_energy_system/analysis.yaml index 287a31c..23bc1c4 100644 --- a/examples/cases/windio_4turbines_profiles_stable/wind_energy_system/analysis.yaml +++ b/examples/cases/windio_4turbines_profiles_stable/wind_energy_system/analysis.yaml @@ -1,9 +1,9 @@ #pywake and foxes wind_deficit_model: name: Bastankhah2014 - wake_expansion_coefficient: # k = ka*ti + kb - k_a: 0.0 - k_b: 0.04 + wake_expansion_coefficient: # k = k_a + k_b*ti + k_a: 0.04 + k_b: 0.0 free_stream_ti: false ceps: 0.2 use_effective_ws: true diff --git a/examples/cases/windio_4turbines_profiles_stable/wind_energy_system/analysis_grid.yaml b/examples/cases/windio_4turbines_profiles_stable/wind_energy_system/analysis_grid.yaml index 205b38e..35a48cf 100644 --- a/examples/cases/windio_4turbines_profiles_stable/wind_energy_system/analysis_grid.yaml +++ b/examples/cases/windio_4turbines_profiles_stable/wind_energy_system/analysis_grid.yaml @@ -1,9 +1,9 @@ #pywake and foxes wind_deficit_model: name: Bastankhah2014 - wake_expansion_coefficient: # k = ka*ti + kb - k_a: 0.0 - k_b: 0.04 + wake_expansion_coefficient: # k = k_a + k_b*ti + k_a: 0.04 + k_b: 0.0 free_stream_ti: false ceps: 0.2 use_effective_ws: true diff --git a/pyproject.toml b/pyproject.toml index a707e99..249789a 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -38,7 +38,7 @@ classifiers = [ ] requires-python = ">=3.10,<3.13" dependencies = [ - "windIO @ git+https://github.com/EUFlow/windIO.git", + "windIO @ git+https://github.com/bjarketol/windIO.git@flow-model-chain", "xarray", "scipy", "pyyaml", @@ -70,6 +70,7 @@ dev = [ "ncplot", "nctoolkit", "cartopy", + "pre-commit", ] docs = [ "sphinx>=7.0", diff --git a/tests/conftest.py b/tests/conftest.py index 0534492..c5d40a8 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -136,7 +136,7 @@ def cleanup_test_outputs(): _ANALYSIS = { "wind_deficit_model": { "name": "Jensen", - "wake_expansion_coefficient": {"k_a": 0.0, "k_b": 0.04}, + "wake_expansion_coefficient": {"k_a": 0.04, "k_b": 0.0}, }, "deflection_model": {"name": "None"}, "turbulence_model": {"name": "STF2005", "c1": 1.0, "c2": 1.0}, @@ -149,6 +149,179 @@ def cleanup_test_outputs(): } +def make_mixed_type_timeseries_system_dict(flow_model_name): + """Build a system dict with TWO turbine types at TWO hub heights. + + 3x3 layout (9 turbines), per-turbine timeseries inflow with dims + ["time", "wind_turbine"] and NO vertical ``height`` profile (the wind + speeds are already at each turbine's own hub height, as a microscale + terrain-speedup model would produce). This exercises the mixed + hub-height + per-turbine path in the pyWake API. + + Type 0 is the DTU 10MW at 119 m; type 1 is a smaller turbine (half the + power, 120 m rotor, 90 m hub) so both the type assignment and the two + hub heights are genuinely distinct. Wind speeds are kept in 8-11 m/s, + away from cut-in/out, so the comparison is unaffected by curve edges. + """ + n_times = 5 + n_turbines = 9 + + # 3x3 grid, ~5D spacing of the larger rotor (5 * 178.3) + spacing = 5 * _TURBINE["rotor_diameter"] + axis = [0.0, spacing, 2 * spacing] + grid_x = [axis[i] for _ in range(3) for i in range(3)] + grid_y = [axis[j] for j in range(3) for _ in range(3)] + + # Alternating type assignment -> [0, 1, 0, 1, 0, 1, 0, 1, 0] + turbine_type_idx = [k % 2 for k in range(n_turbines)] + + perf = _TURBINE["performance"] + pc_ws = perf["power_curve"]["power_wind_speeds"] + pc_p = perf["power_curve"]["power_values"] + ct_ws = perf["Ct_curve"]["Ct_wind_speeds"] + ct_v = perf["Ct_curve"]["Ct_values"] + type1_power = [0.5 * p for p in pc_p] + + # Deterministic per-turbine timeseries (no RNG), wind from ~270 deg so the + # rows along x interact. + base_ws = [8.0, 9.0, 10.0, 8.5, 11.0] + base_wd = [270.0, 268.0, 272.0, 270.5, 269.0] + ws_data, wd_data, ti_data = [], [], [] + for t in range(n_times): + ws_data.append([base_ws[t] + 0.1 * i for i in range(n_turbines)]) + wd_data.append([base_wd[t] + 0.2 * (i - 4) for i in range(n_turbines)]) + ti_data.append([0.06 + 0.002 * i for i in range(n_turbines)]) + + return { + "name": "Dict test: mixed turbine types and hub heights", + "site": { + "name": "Test site", + "boundaries": { + "polygons": [ + { + "x": [-spacing, 3 * spacing, 3 * spacing, -spacing], + "y": [3 * spacing, 3 * spacing, -spacing, -spacing], + } + ] + }, + "energy_resource": { + "name": "Test resource", + "wind_resource": { + "time": list(range(n_times)), + "wind_turbine": list(range(n_turbines)), + "wind_speed": {"data": ws_data, "dims": ["time", "wind_turbine"]}, + "wind_direction": { + "data": wd_data, + "dims": ["time", "wind_turbine"], + }, + "turbulence_intensity": { + "data": ti_data, + "dims": ["time", "wind_turbine"], + }, + }, + }, + }, + "wind_farm": { + "name": "Test farm", + "layouts": [ + { + "coordinates": {"x": grid_x, "y": grid_y}, + "turbine_types": turbine_type_idx, + } + ], + "turbine_types": { + 0: { + "name": "type_0_DTU10MW", + "hub_height": 119.0, + "rotor_diameter": 178.3, + "performance": { + "power_curve": { + "power_wind_speeds": pc_ws, + "power_values": pc_p, + }, + "Ct_curve": {"Ct_wind_speeds": ct_ws, "Ct_values": ct_v}, + }, + }, + 1: { + "name": "type_1_small", + "hub_height": 90.0, + "rotor_diameter": 120.0, + "performance": { + "power_curve": { + "power_wind_speeds": pc_ws, + "power_values": type1_power, + }, + "Ct_curve": {"Ct_wind_speeds": ct_ws, "Ct_values": ct_v}, + }, + }, + }, + }, + "attributes": { + "flow_model": {"name": flow_model_name}, + "analysis": { + "wind_deficit_model": { + "name": "Bastankhah2014", + "wake_expansion_coefficient": { + "k_a": 0.04, + "k_b": 0.0, + "free_stream_ti": False, + }, + "ceps": 0.2, + "use_effective_ws": True, + }, + "axial_induction_model": "Madsen", + "deflection_model": {"name": "None"}, + "turbulence_model": {"name": "CrespoHernandez"}, + "superposition_model": { + "ws_superposition": "Linear", + "ti_superposition": "Squared", + }, + "rotor_averaging": {"name": "Center"}, + "blockage_model": {"name": "None"}, + }, + "model_outputs_specification": { + "turbine_outputs": { + "turbine_nc_filename": "turbine_data.nc", + "output_variables": ["power"], + }, + }, + }, + } + + +def make_mixed_type_profile_system_dict(flow_model_name): + """Mixed types/hub heights with a (time, height) vertical wind profile. + + Same farm as :func:`make_mixed_type_timeseries_system_dict` (two types at + 119 m and 90 m), but the inflow is given as a vertical profile over a + ``height`` dimension that brackets both hub heights, so WIFA must + interpolate the profile to each turbine's hub height. The wind direction + veers around 360 deg to exercise circular interpolation. + """ + system = make_mixed_type_timeseries_system_dict(flow_model_name) + + n_times = 5 + heights = [60.0, 90.0, 120.0, 150.0] + href, alpha = 119.0, 0.14 # power-law shear reference + exponent + ws_base = [8.0, 9.0, 10.0, 8.5, 11.0] + wd_base = [358.0, 359.0, 357.0, 0.5, 359.5] # straddles 360 deg + + ws_data, wd_data, ti_data = [], [], [] + for t in range(n_times): + ws_data.append([ws_base[t] * (h / href) ** alpha for h in heights]) + wd_data.append([(wd_base[t] + 0.05 * (h - href)) % 360 for h in heights]) + ti_data.append([max(0.08 - 0.0001 * (h - 60.0), 0.02) for h in heights]) + + system["site"]["energy_resource"]["wind_resource"] = { + "time": list(range(n_times)), + "height": heights, + "wind_speed": {"data": ws_data, "dims": ["time", "height"]}, + "wind_direction": {"data": wd_data, "dims": ["time", "height"]}, + "turbulence_intensity": {"data": ti_data, "dims": ["time", "height"]}, + } + return system + + def make_timeseries_per_turbine_system_dict(flow_model_name): """Build a complete system dict with per-turbine timeseries data including density. diff --git a/tests/mem_bench.py b/tests/mem_bench.py new file mode 100644 index 0000000..3b0e63b --- /dev/null +++ b/tests/mem_bench.py @@ -0,0 +1,266 @@ +"""Memory benchmark for WIFA site-data construction. + +Measures the peak Python allocation of turning the windIO ``wind_resource`` +dict into a pyWake site, isolating the copies this optimization targets: + * the dict-of-lists footprint (what windIO's _ds2yml currently produces) + * the extra peak allocated inside ``construct_site`` (np.array / [cases_idx] + / dict_to_netcdf copies) + * the Weibull path's duplicate ``dict_to_netcdf`` build + +Run: pixi run python tests/mem_bench.py +""" + +import gc +import tempfile +import tracemalloc +from pathlib import Path + +import numpy as np +import xarray as xr + +from wifa.pywake_api import construct_site, create_turbines + +_TURBINE = { + "name": "generic", + "hub_height": 100.0, + "rotor_diameter": 120.0, + "performance": { + "power_curve": { + "power_wind_speeds": [3.0, 8.0, 12.0, 25.0], + "power_values": [0.0, 1.0e6, 3.0e6, 3.0e6], + }, + "Ct_curve": { + "Ct_wind_speeds": [3.0, 8.0, 12.0, 25.0], + "Ct_values": [0.8, 0.8, 0.4, 0.2], + }, + }, +} + +_ANALYSIS = { + "wind_deficit_model": { + "name": "Jensen", + "wake_expansion_coefficient": {"k_a": 0.04, "k_b": 0.0}, + }, + "deflection_model": {"name": "None"}, + "turbulence_model": {"name": "STF2005", "c1": 1.0, "c2": 1.0}, + "superposition_model": { + "ws_superposition": "Linear", + "ti_superposition": "Squared", + }, + "rotor_averaging": {"name": "Center"}, + "blockage_model": {"name": "None"}, +} + + +def _attrs(): + return { + "flow_model": {"name": "pywake"}, + "analysis": _ANALYSIS, + "model_outputs_specification": { + "turbine_outputs": { + "turbine_nc_filename": "turbine_data.nc", + "output_variables": ["power"], + } + }, + } + + +def _layout(n_turbines, spacing=600.0): + return { + "coordinates": { + "x": [i * spacing for i in range(n_turbines)], + "y": [0.0] * n_turbines, + } + } + + +def _bounds(n_turbines, spacing=600.0): + hi = n_turbines * spacing + 100 + return {"polygons": [{"x": [-100, hi, hi, -100], "y": [100, 100, -100, -100]}]} + + +def timeseries_system(n_times, n_turbines, arrays=False): + """Per-turbine time-series system. + + ``arrays=False`` emits a dict-of-Python-lists (what windIO's _ds2yml + produces today). ``arrays=True`` emits a dict-of-ndarrays (what the + Phase 2 ``to_dict(data="array")`` loader would produce) to preview the win. + """ + rng = np.random.default_rng(0) + _c = (lambda a: a) if arrays else (lambda a: a.tolist()) + ws = _c(8.0 + 4.0 * rng.random((n_times, n_turbines))) + wd = _c(270.0 + 10.0 * rng.random((n_times, n_turbines))) + ti = _c(0.06 + 0.04 * rng.random((n_times, n_turbines))) + rho = _c(1.18 + 0.06 * rng.random((n_times, n_turbines))) + op = _c(np.ones((n_times, n_turbines), dtype=int)) + dims = ["time", "wind_turbine"] + return { + "name": "bench timeseries", + "site": { + "name": "s", + "boundaries": _bounds(n_turbines), + "energy_resource": { + "name": "r", + "wind_resource": { + "time": list(range(n_times)), + "wind_turbine": list(range(n_turbines)), + "wind_speed": {"data": ws, "dims": dims}, + "wind_direction": {"data": wd, "dims": dims}, + "turbulence_intensity": {"data": ti, "dims": dims}, + "density": {"data": rho, "dims": dims}, + "operating": {"data": op, "dims": dims}, + }, + }, + }, + "wind_farm": { + "name": "f", + "layouts": [_layout(n_turbines)], + "turbines": _TURBINE, + }, + "attributes": _attrs(), + } + + +def weibull_system(n_dirs, n_turbines): + """Per-turbine Weibull system as a dict-of-lists (exercises the Weibull path).""" + rng = np.random.default_rng(1) + a = (8.0 + 2.0 * rng.random((n_turbines, n_dirs))).tolist() + k = (2.0 + 0.3 * rng.random((n_turbines, n_dirs))).tolist() + p = rng.random((n_turbines, n_dirs)) + p = (p / p.sum(axis=1, keepdims=True)).tolist() + wd = np.linspace(0, 360, n_dirs, endpoint=False).tolist() + dims = ["wind_turbine", "wind_direction"] + return { + "name": "bench weibull", + "site": { + "name": "s", + "boundaries": _bounds(n_turbines), + "energy_resource": { + "name": "r", + "wind_resource": { + "wind_direction": wd, + "wind_turbine": list(range(n_turbines)), + "weibull_a": {"data": a, "dims": dims}, + "weibull_k": {"data": k, "dims": dims}, + "sector_probability": {"data": p, "dims": dims}, + "turbulence_intensity": { + "data": (0.07 * np.ones((n_turbines, n_dirs))).tolist(), + "dims": dims, + }, + }, + }, + }, + "wind_farm": { + "name": "f", + "layouts": [_layout(n_turbines)], + "turbines": _TURBINE, + }, + "attributes": _attrs(), + } + + +def _peak_of(fn): + tracemalloc.start() + tracemalloc.clear_traces() + obj = fn() + cur, peak = tracemalloc.get_traced_memory() + tracemalloc.stop() + return obj, cur, peak + + +def _construct_peak(system): + farm = system["wind_farm"] + _, _, hub_heights = create_turbines(farm) + x = farm["layouts"][0]["coordinates"]["x"] + resource = system["site"]["energy_resource"] + tracemalloc.start() + tracemalloc.reset_peak() + base, _ = tracemalloc.get_traced_memory() + construct_site(system, resource, hub_heights, x) + _, peak = tracemalloc.get_traced_memory() + tracemalloc.stop() + return peak - base + + +def nc_load_peaks(n_times, n_turbines): + """Peak memory of windIO.load_yaml on a real wind_resource.nc, list vs array. + + Requires a windIO build with the ``nc_data`` loader option; returns None if + unavailable so the benchmark still runs the in-process cases. + """ + import windIO + + d = Path(tempfile.mkdtemp()) + rng = np.random.default_rng(0) + ds = xr.Dataset( + { + v: (("time", "wind_turbine"), rng.random((n_times, n_turbines))) + for v in ("wind_speed", "wind_direction", "turbulence_intensity", "density") + }, + coords={"time": np.arange(n_times), "wind_turbine": np.arange(n_turbines)}, + ) + ds["operating"] = ( + ("time", "wind_turbine"), + np.ones((n_times, n_turbines), dtype=int), + ) + ds.to_netcdf(d / "wind_resource.nc") + sysf = d / "system.yaml" + sysf.write_text("wind_resource: !include wind_resource.nc\n") + + def peak(mode): + gc.collect() + tracemalloc.start() + tracemalloc.clear_traces() + obj = windIO.load_yaml(sysf, nc_data=mode) + _, pk = tracemalloc.get_traced_memory() + tracemalloc.stop() + del obj + gc.collect() + return pk + + try: + return peak("list"), peak("array") + except TypeError: + return None # windIO without the nc_data option + + +def main(): + MB = 1e6 + print(f"{'case':<28}{'input dict':>14}{'construct_site':>18}") + print("-" * 60) + + n_times, n_wt = 4000, 100 + nbytes = n_times * n_wt * 8 * 5 / MB # 5 float arrays + sysd, _, dpeak = _peak_of(lambda: timeseries_system(n_times, n_wt)) + cs = _construct_peak(sysd) + print( + f"{'timeseries 4000x100 (lists)':<28}{dpeak / MB:>11.1f}MB{cs / MB:>15.1f}MB" + f" (raw ndarray ~{nbytes:.1f}MB)" + ) + + sysa, _, apeak = _peak_of(lambda: timeseries_system(n_times, n_wt, arrays=True)) + csa = _construct_peak(sysa) + print( + f"{'timeseries 4000x100 (arrays)':<28}{apeak / MB:>11.1f}MB{csa / MB:>15.1f}MB" + f" <- Phase 2 preview" + ) + + n_dirs, n_wt_w = 12, 2000 + sysw, _, wpeak = _peak_of(lambda: weibull_system(n_dirs, n_wt_w)) + csw = _construct_peak(sysw) + print(f"{'weibull 2000wt x 12dir':<28}{wpeak / MB:>11.1f}MB{csw / MB:>15.1f}MB") + + print() + peaks = nc_load_peaks(n_times, n_wt) + if peaks is None: + print("load_yaml(.nc): windIO has no nc_data option (skipped)") + else: + pl, pa = peaks + print( + f"load_yaml(.nc) 4000x100 list={pl / MB:.1f}MB " + f"array={pa / MB:.1f}MB ({pl / pa:.1f}x lower peak)" + ) + + +if __name__ == "__main__": + main() 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_memory.py b/tests/test_memory.py new file mode 100644 index 0000000..e455f61 --- /dev/null +++ b/tests/test_memory.py @@ -0,0 +1,71 @@ +"""Regression guard for site-data load memory. + +Keeping an included ``wind_resource.nc`` as numpy arrays (windIO's +``nc_data="array"``) must stay dramatically cheaper than the dict-of-lists +default for a large resource. See tests/mem_bench.py for the full breakdown. +""" + +import gc +import tracemalloc + +import numpy as np +import pytest +import windIO +import xarray as xr + +pytest.importorskip("netCDF4") + + +def _write_resource(tmp_path, n_times, n_turbines): + rng = np.random.default_rng(0) + ds = xr.Dataset( + { + v: (("time", "wind_turbine"), rng.random((n_times, n_turbines))) + for v in ("wind_speed", "wind_direction", "turbulence_intensity", "density") + }, + coords={"time": np.arange(n_times), "wind_turbine": np.arange(n_turbines)}, + ) + ds["operating"] = ( + ("time", "wind_turbine"), + np.ones((n_times, n_turbines), dtype=int), + ) + ds.to_netcdf(tmp_path / "wind_resource.nc") + sysf = tmp_path / "system.yaml" + sysf.write_text("wind_resource: !include wind_resource.nc\n") + return sysf + + +def _peak_load(sysf, mode): + gc.collect() + tracemalloc.start() + tracemalloc.clear_traces() + obj = windIO.load_yaml(sysf, nc_data=mode) + _, peak = tracemalloc.get_traced_memory() + tracemalloc.stop() + del obj + gc.collect() + return peak + + +def test_array_loader_reduces_peak_memory(tmp_path): + """Array-backed load must peak well below the dict-of-lists load.""" + sysf = _write_resource(tmp_path, 2000, 60) + + list_peak = _peak_load(sysf, "list") + array_peak = _peak_load(sysf, "array") + + # Measured ~4x; require a comfortable >=2x to avoid flakiness. + assert array_peak < 0.5 * list_peak, ( + f"array load peak {array_peak / 1e6:.1f}MB not < 0.5 x list " + f"{list_peak / 1e6:.1f}MB" + ) + + +def test_array_loader_keeps_ndarrays(tmp_path): + """The array loader must actually keep numpy arrays (not lists).""" + sysf = _write_resource(tmp_path, 100, 4) + wr = windIO.load_yaml(sysf, nc_data="array")["wind_resource"] + assert isinstance(wr["wind_speed"]["data"], np.ndarray) + # default stays list-backed (backwards compatible) + wr_list = windIO.load_yaml(sysf)["wind_resource"] + assert isinstance(wr_list["wind_speed"]["data"], list) diff --git a/tests/test_pywake.py b/tests/test_pywake.py index 2605cf5..c3c71b2 100644 --- a/tests/test_pywake.py +++ b/tests/test_pywake.py @@ -18,7 +18,7 @@ from py_wake.superposition_models import LinearSum from py_wake.tests import npt from py_wake.turbulence_models import CrespoHernandez -from py_wake.wind_turbines import WindTurbine +from py_wake.wind_turbines import WindTurbine, WindTurbines from py_wake.wind_turbines.power_ct_functions import PowerCtFunctionList, PowerCtTabular from scipy.special import gamma from windIO import __path__ as wiop @@ -274,13 +274,32 @@ def test_heterogeneous_wind_rose_grid(): x = [0, 1248.1, 2496.2, 3744.3] y = [0, 0, 0, 0] - # compute AEP with PyWake - res_aep = ( - wfm(x, y, ws=np.arange(2, 30, 1), wd=dat["wd"]) - .aep(normalize_probabilities=True) - .sum() + # 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 grid starts at 0.5 (not 0): WIFA drops the degenerate ws=0 reference + # case, which breaks the WeightedSum superposition (see + # test_weibull_ws_grid_excludes_zero_for_weightedsum). + ws_range = np.arange(0.5, 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=ws_range, wd=wd_fine).aep(normalize_probabilities=True).sum() + # compute AEP with API wifa_res = run_pywake( test_path @@ -447,6 +466,225 @@ def test_turbine_specific_speeds_timeseries(): npt.assert_allclose(wifa_res, manual_aep, rtol=1e-6) +def test_pywake_mixed_turbine_types_hub_heights(tmp_path): + """Two turbine types at two hub heights with per-turbine timeseries inflow. + + WIFA must reproduce a hand-built pyWake simulation. This is the golden + equivalence check for the mixed-type / mixed-hub-height path: the API + builds the per-turbine XRSite (via dict_to_site) inside the + ``len(hub_heights) > 1`` branch and assigns turbine types, and the + resulting per-turbine power time series must match raw pyWake. + """ + from conftest import make_mixed_type_timeseries_system_dict + + system_dict = make_mixed_type_timeseries_system_dict("pywake") + output_dir = tmp_path / "output_pywake_mixed" + aep_wifa = run_pywake(system_dict, output_dir=str(output_dir)) + + # --- Build the reference pyWake simulation directly --- + wr = system_dict["site"]["energy_resource"]["wind_resource"] + farm = system_dict["wind_farm"] + coords = farm["layouts"][0]["coordinates"] + x, y = coords["x"], coords["y"] + type_list = farm["layouts"][0]["turbine_types"] + + ws = np.array(wr["wind_speed"]["data"]) # (time, turbine) + wd = np.array(wr["wind_direction"]["data"]) + ti = np.array(wr["turbulence_intensity"]["data"]) + n_time, n_wt = ws.shape + + # Per-turbine site (mirrors dict_to_site: wind_turbine -> i, i leading dim, + # uniform P, integer time). + ds = xr.Dataset( + { + "WS": (("time", "i"), ws), + "WD": (("time", "i"), wd), + "TI": (("time", "i"), ti), + }, + coords={"time": np.arange(n_time), "i": np.arange(n_wt)}, + ).transpose("i", "time") + ds["P"] = (("time",), np.ones(n_time) / n_time) + site = XRSite(ds, interp_method="linear") + + # Two turbine types matching the windIO definitions (WIFA interpolates the + # curves onto integer wind speeds; the nodes here are already integer). + tdefs = farm["turbine_types"] + turbines = [] + for k in sorted(tdefs): + td = tdefs[k] + pc = td["performance"]["power_curve"] + ctc = td["performance"]["Ct_curve"] + speeds = np.arange( + min(pc["power_wind_speeds"][0], ctc["Ct_wind_speeds"][0]), + max(pc["power_wind_speeds"][-1], ctc["Ct_wind_speeds"][-1]) + 1, + 1, + ) + powers = np.interp(speeds, pc["power_wind_speeds"], pc["power_values"]) + cts = np.interp(speeds, ctc["Ct_wind_speeds"], ctc["Ct_values"]) + turbines.append( + WindTurbine( + name=td["name"], + diameter=td["rotor_diameter"], + hub_height=td["hub_height"], + powerCtFunction=PowerCtTabular(speeds, powers, "W", ct=cts), + ) + ) + turbine = WindTurbines.from_WindTurbine_lst(turbines) + + wfm = BastankhahGaussian( + site, + turbine, + k=0.04, + ceps=0.2, + superpositionModel=LinearSum(), + use_effective_ws=True, + turbulenceModel=CrespoHernandez(), + ) + + # Reference inflow arrays, exactly as the API reduces them (mean over + # turbines, vector mean for direction). + ws_ref = ds.WS.mean(dim="i").values + rads = np.deg2rad(ds.WD) + wd_ref = ( + np.rad2deg(np.arctan2(np.sin(rads).mean("i"), np.cos(rads).mean("i"))) % 360 + ).values + + res = wfm(x, y, type=type_list, time=True, ws=ws_ref, wd=wd_ref) + ref_power = res.Power.transpose("wt", "time").values + ref_aep = res.aep(normalize_probabilities=False).sum() + + # --- Compare WIFA output (per-turbine power) to the reference --- + with xr.open_dataset(output_dir / "turbine_data.nc") as out: + wifa_power = out["power"].transpose("turbine", "time").values + + npt.assert_allclose(wifa_power, ref_power, rtol=1e-6, atol=1.0) + npt.assert_allclose(aep_wifa, ref_aep, rtol=1e-6) + + # The case really does exercise two distinct hub heights. + assert {td["hub_height"] for td in tdefs.values()} == {119.0, 90.0} + + +def test_pywake_mixed_types_vertical_profile(tmp_path): + """Two turbine types at two hub heights driven by a (time, height) vertical + profile. WIFA interpolates the profile to each turbine's hub height and must + match a hand-built per-turbine pyWake reference (per-turbine power).""" + from conftest import make_mixed_type_profile_system_dict + + system_dict = make_mixed_type_profile_system_dict("pywake") + output_dir = tmp_path / "output_pywake_profile" + aep_wifa = run_pywake(system_dict, output_dir=str(output_dir)) + + wr = system_dict["site"]["energy_resource"]["wind_resource"] + farm = system_dict["wind_farm"] + coords = farm["layouts"][0]["coordinates"] + x, y = coords["x"], coords["y"] + type_list = farm["layouts"][0]["turbine_types"] + + heights = np.array(wr["height"], dtype=float) + ws_prof = np.array(wr["wind_speed"]["data"]) # (time, height) + wd_prof = np.array(wr["wind_direction"]["data"]) + ti_prof = np.array(wr["turbulence_intensity"]["data"]) + n_time = ws_prof.shape[0] + + tdefs = farm["turbine_types"] + ordered_keys = sorted(tdefs) + ordered_hh = [tdefs[k]["hub_height"] for k in ordered_keys] + + # Interpolate the profile to each turbine's hub height: linear WS/TI, + # vector (sin/cos) WD — mirroring the WIFA helpers. + def _lin(prof, hub): + return np.array([np.interp(hub, heights, prof[t]) for t in range(n_time)]) + + def _dir(prof, hub): + rad = np.deg2rad(prof) + s = np.array([np.interp(hub, heights, np.sin(rad)[t]) for t in range(n_time)]) + c = np.array([np.interp(hub, heights, np.cos(rad)[t]) for t in range(n_time)]) + return np.mod(np.rad2deg(np.arctan2(s, c)), 360.0) + + n_wt = len(type_list) + hub_of = [ordered_hh[type_list[i]] for i in range(n_wt)] + ws_i = np.array([_lin(ws_prof, hub_of[i]) for i in range(n_wt)]) + wd_i = np.array([_dir(wd_prof, hub_of[i]) for i in range(n_wt)]) + ti_i = np.maximum(np.array([_lin(ti_prof, hub_of[i]) for i in range(n_wt)]), 0.02) + + ds = xr.Dataset( + { + "WS": (("i", "time"), ws_i), + "WD": (("i", "time"), wd_i), + "TI": (("i", "time"), ti_i), + }, + coords={"i": np.arange(n_wt), "time": np.arange(n_time)}, + ) + ds["P"] = (("time",), np.ones(n_time) / n_time) + site = XRSite(ds, interp_method="linear") + + turbines = [] + for k in ordered_keys: + td = tdefs[k] + pc = td["performance"]["power_curve"] + ctc = td["performance"]["Ct_curve"] + speeds = np.arange( + min(pc["power_wind_speeds"][0], ctc["Ct_wind_speeds"][0]), + max(pc["power_wind_speeds"][-1], ctc["Ct_wind_speeds"][-1]) + 1, + 1, + ) + powers = np.interp(speeds, pc["power_wind_speeds"], pc["power_values"]) + cts = np.interp(speeds, ctc["Ct_wind_speeds"], ctc["Ct_values"]) + turbines.append( + WindTurbine( + name=td["name"], + diameter=td["rotor_diameter"], + hub_height=td["hub_height"], + powerCtFunction=PowerCtTabular(speeds, powers, "W", ct=cts), + ) + ) + turbine = WindTurbines.from_WindTurbine_lst(turbines) + + wfm = BastankhahGaussian( + site, + turbine, + k=0.04, + ceps=0.2, + superpositionModel=LinearSum(), + use_effective_ws=True, + turbulenceModel=CrespoHernandez(), + ) + ws_ref = ds.WS.mean("i").values + rads = np.deg2rad(ds.WD) + wd_ref = ( + np.rad2deg(np.arctan2(np.sin(rads).mean("i"), np.cos(rads).mean("i"))) % 360 + ).values + + res = wfm(x, y, type=type_list, time=True, ws=ws_ref, wd=wd_ref) + ref_power = res.Power.transpose("wt", "time").values + ref_aep = res.aep(normalize_probabilities=False).sum() + + with xr.open_dataset(output_dir / "turbine_data.nc") as out: + wifa_power = out["power"].transpose("turbine", "time").values + + npt.assert_allclose(wifa_power, ref_power, rtol=1e-6, atol=1.0) + npt.assert_allclose(aep_wifa, ref_aep, rtol=1e-6) + + # Two distinct hub heights, at least one strictly between profile levels + # (so the interpolation is genuinely exercised, not just node selection). + assert len(set(ordered_hh)) == 2 + assert any(h not in set(heights.tolist()) for h in ordered_hh) + + +def test_pywake_mixed_height_and_per_turbine_raises(tmp_path): + """A resource carrying BOTH a height profile and a wind_turbine dimension is + ambiguous and must fail loudly rather than guess.""" + from conftest import make_mixed_type_timeseries_system_dict + + system_dict = make_mixed_type_timeseries_system_dict("pywake") + # The per-turbine resource already has wind_turbine-dimensioned ws/wd; + # add a height axis alongside it to create the ambiguous combination. + system_dict["site"]["energy_resource"]["wind_resource"]["height"] = [80.0, 120.0] + + with pytest.raises(NotImplementedError, match="both a 'height' profile"): + run_pywake(system_dict, output_dir=str(tmp_path / "output_pywake_ambig")) + + def test_pywake_dict_timeseries_per_turbine_with_density(tmp_path): from conftest import make_timeseries_per_turbine_system_dict @@ -466,6 +704,231 @@ 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) + + +def test_weibull_ws_grid_excludes_zero_for_weightedsum(tmp_path): + """Regression: the auto-generated Weibull ws grid must exclude ws=0. + + When no explicit ``wind_speed`` is given, ``_construct_weibull_site`` + builds a reference ws grid. It used to start at 0 m/s. A ws=0 flow case + carries zero energy for every model, but it is degenerate for the + ``WeightedSum`` superposition (Zong 2020), whose convection-velocity + iteration divides by the convection speed and is undefined at zero wind + speed. Including ws=0 silently corrupted the WeightedSum AEP — collapsing + the apparent wake loss (e.g. Zong fell to ~3.5 % while the near-identical + LinearSum Niayifar stayed at ~8 %). LinearSum and the other superpositions + were unaffected, so the bug only surfaced on the distributions path with + Weighted superposition. + + Guards both the grid (ws[0] > 0) and the behaviour (Weighted must not + diverge from Linear for an otherwise-identical Zong farm). + """ + from conftest import _TURBINE + from wifa.pywake_api import ( + construct_site, + create_turbines, + ) + + n_wd, n_wt = 4, 5 + wd_vals = [0.0, 90.0, 180.0, 270.0] + A = [[9.0] * n_wt for _ in range(n_wd)] + k = [[2.0] * n_wt for _ in range(n_wd)] + freq = [[1.0 / n_wd] * n_wt for _ in range(n_wd)] + ti = [[0.06] * n_wt for _ in range(n_wd)] + + def make_system(superposition, rotor): + return { + "name": "ws0-regression", + "site": { + "name": "Test site", + "boundaries": { + "polygons": [ + {"x": [-90, 5000, 5000, -90], "y": [90, 90, -90, -90]} + ] + }, + "energy_resource": { + "name": "Test resource", + "wind_resource": { + # NOTE: deliberately NO "wind_speed" -> auto ws grid + "wind_direction": wd_vals, + "wind_turbine": list(range(n_wt)), + "reference_height": 119.0, + "weibull_a": { + "data": A, + "dims": ["wind_direction", "wind_turbine"], + }, + "weibull_k": { + "data": k, + "dims": ["wind_direction", "wind_turbine"], + }, + "sector_probability": { + "data": freq, + "dims": ["wind_direction", "wind_turbine"], + }, + "turbulence_intensity": { + "data": ti, + "dims": ["wind_direction", "wind_turbine"], + }, + }, + }, + }, + "wind_farm": { + "name": "Test farm", + "layouts": [ + { + "coordinates": { + # 5-turbine row at 5D spacing -> strong aligned wakes + "x": [i * 5 * _TURBINE["rotor_diameter"] for i in range(5)], + "y": [0.0] * 5, + } + } + ], + "turbines": _TURBINE, + }, + "attributes": { + "flow_model": {"name": "pywake"}, + "analysis": { + "wind_deficit_model": { + "name": "Zong2020", + "wake_expansion_coefficient": {"k_a": 0.004, "k_b": 0.38}, + }, + "deflection_model": {"name": "None"}, + "turbulence_model": {"name": "CrespoHernandez"}, + "superposition_model": { + "ws_superposition": superposition, + "ti_superposition": "Squared", + }, + "rotor_averaging": {"name": rotor}, + "blockage_model": {"name": "None"}, + "axial_induction_model": "Madsen", + }, + }, + } + + # 1. The constructed ws grid must not contain a 0 m/s reference case. + weighted = make_system("Weighted", "none") + turbine, _types, hub_heights = create_turbines(weighted["wind_farm"]) + x = weighted["wind_farm"]["layouts"][0]["coordinates"]["x"] + site_data = construct_site( + weighted, weighted["site"]["energy_resource"], hub_heights, x + ) + ws_grid = np.asarray(site_data["ws"]) + assert ws_grid[0] > 0.0, f"ws grid must exclude 0; starts at {ws_grid[0]}" + + # 2. Behaviour: WeightedSum must not collapse relative to LinearSum on the + # same Zong farm (pre-fix the Weighted AEP was inflated by the ws=0 bug). + aep_weighted = run_pywake(weighted, output_dir=str(tmp_path / "weighted")) + aep_linear = run_pywake( + make_system("Linear", "Center"), output_dir=str(tmp_path / "linear") + ) + npt.assert_allclose(aep_weighted, aep_linear, rtol=0.03) + + # 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/tests/test_pywake_eddy_viscosity.py b/tests/test_pywake_eddy_viscosity.py new file mode 100644 index 0000000..9e5ab9a --- /dev/null +++ b/tests/test_pywake_eddy_viscosity.py @@ -0,0 +1,104 @@ +"""Tests for the Eddy Viscosity (EV) bundle dispatch in wifa/pywake_api.py. + +The EV model exists only on pyWake's unmerged ``cj_add_eddy_viscosity_model`` +branch (mirrored at eu-flow/wp5/PyWake@flow-model-chain), so this file lives +apart from test_pywake_submodels.py and every test is skipped on a released +pyWake. +""" + +import pytest + +eddy_viscosity = pytest.importorskip( + "py_wake.deficit_models.eddy_viscosity", + reason="EV requires pyWake's cj_add_eddy_viscosity_model branch", +) + +from py_wake.rotor_avg_models.simplified_gaussian_rotor_average_model import ( # noqa: E402 + SimplifiedGaussianRotorAverageModel, +) +from py_wake.turbulence_models.quarton_and_ainslie import ( # noqa: E402 + ModifiedQuartonAndAinslieTurbulenceModel, +) + +from wifa.pywake_api import ( # noqa: E402 + _configure_deficit_model, + _configure_rotor_averaging, + _configure_turbulence_model, +) + +_RD = 126.0 +_HH = 90.0 + + +def _call_deficit(name, analysis_extra=None): + wind_deficit_model = {"name": name, **(analysis_extra or {})} + analysis = {"wind_deficit_model": wind_deficit_model} + cls, args, _post = _configure_deficit_model({"name": name}, analysis, _RD, _HH) + return cls, args + + +@pytest.mark.parametrize( + "name", ["EddyViscosity", "eddyviscosity", "EDDYVISCOSITY", "eddy_viscosity"] +) +def test_configure_deficit_model_eddy_viscosity(name): + cls, args = _call_deficit(name) + assert cls is eddy_viscosity.EddyViscosityDeficitModel + # EV pairs with MaxSum: deficits are referenced to the free stream, so the + # usual use_effective_ws=True default flips to False for this model. + assert args["use_effective_ws"] is False + + +def test_eddy_viscosity_use_effective_ws_override(): + _cls, args = _call_deficit("EddyViscosity", {"use_effective_ws": True}) + assert args["use_effective_ws"] is True + + +def test_eddy_viscosity_free_stream_ti_mapped(): + # EV is TI-capable: windIO's free_stream_ti flag maps to use_effective_ti. + _cls, args = _call_deficit( + "EddyViscosity", {"wake_expansion_coefficient": {"free_stream_ti": True}} + ) + assert args["use_effective_ti"] is False + + +def test_eddy_viscosity_no_ct2a_injected(): + # EddyViscosityDeficitModel has no ct2a parameter; the axial-induction + # mapping must skip it rather than pass an unexpected kwarg. + wind_deficit_model = {"name": "EddyViscosity"} + analysis = { + "wind_deficit_model": wind_deficit_model, + "axial_induction_model": "1D", + } + _cls, args, _post = _configure_deficit_model( + {"name": "EddyViscosity"}, analysis, _RD, _HH + ) + assert "ct2a" not in args + + +@pytest.mark.parametrize( + "name", + [ + "QuartonAndAinslie", + "quarton_and_ainslie", + "ModifiedQuartonAndAinslie", + "modifiedquartonandainslie", + ], +) +def test_configure_turbulence_model_quarton_and_ainslie(name): + model = _configure_turbulence_model({"name": name}) + assert isinstance(model, ModifiedQuartonAndAinslieTurbulenceModel) + + +@pytest.mark.parametrize("name", ["simplified_gaussian", "SimplifiedGaussian"]) +def test_configure_rotor_averaging_simplified_gaussian(name): + model = _configure_rotor_averaging({"name": name}) + assert isinstance(model, SimplifiedGaussianRotorAverageModel) + + +def test_eddy_viscosity_deficit_constructs(): + # The dispatched class + args must actually build (triggers the ~45 ms + # LUT generation) with the EV bundle defaults. + cls, args = _call_deficit("EddyViscosity") + deficit = cls(**args) + assert deficit.use_effective_ws is False + assert deficit.use_effective_ti is True diff --git a/tests/test_pywake_histogram.py b/tests/test_pywake_histogram.py new file mode 100644 index 0000000..b625985 --- /dev/null +++ b/tests/test_pywake_histogram.py @@ -0,0 +1,318 @@ +"""Tests for the histogram (wind_direction x wind_speed probability) site path.""" + +import os +import shutil +from pathlib import Path + +import numpy as np +import pytest + +pytest.importorskip( + "py_wake", reason="py_wake not installed, install with: pip install wifa[pywake]" +) + +import xarray as xr +from py_wake.deficit_models.gaussian import BastankhahGaussian +from py_wake.examples.data.dtu10mw._dtu10mw import DTU10MW +from py_wake.tests import npt + +from wifa.pywake_api import ( + _subsector_wind_directions, + construct_site, + dict_to_site, + run_pywake, +) + +test_path = Path(os.path.dirname(__file__)) + +N_SECTORS = 12 +WS_CENTERS = np.arange(0.5, 26.0, 1.0) # 26 one-m/s bins, 0.5 .. 25.5 +WD_SECTORS = np.arange(0.0, 360.0, 360.0 / N_SECTORS) + + +def _conditional_ws_pdf(mean_ws): + """Discretized Rayleigh-like pdf on WS_CENTERS, normalized to sum to 1.""" + p = WS_CENTERS * np.exp(-((WS_CENTERS / mean_ws) ** 2)) + return p / p.sum() + + +def _histogram_resource(n_turbines=2, with_ti=True, uniform_wd=False): + """Build a windIO wind_resource dict with a per-turbine 2D histogram. + + ``probability`` is the per-sector conditional over wind speed (windkit + convention), ``sector_probability`` the direction marginal. + """ + rng = np.random.default_rng(42) + if uniform_wd: + sector_p = np.full((n_turbines, N_SECTORS), 1.0 / N_SECTORS) + else: + raw = 1.0 + rng.random((n_turbines, N_SECTORS)) + sector_p = raw / raw.sum(axis=1, keepdims=True) + + cond = np.empty((n_turbines, N_SECTORS, WS_CENTERS.size)) + for t in range(n_turbines): + for s in range(N_SECTORS): + mean_ws = 8.0 + 0.5 * t if uniform_wd else 7.0 + t + 0.2 * s + cond[t, s] = _conditional_ws_pdf(mean_ws) + + resource = { + "wind_direction": WD_SECTORS.tolist(), + "wind_speed": WS_CENTERS.tolist(), + "wind_turbine": list(range(n_turbines)), + "probability": { + "dims": ["wind_turbine", "wind_direction", "wind_speed"], + "data": cond.tolist(), + }, + "sector_probability": { + "dims": ["wind_turbine", "wind_direction"], + "data": sector_p.tolist(), + }, + } + if with_ti: + resource["turbulence_intensity"] = { + "dims": ["wind_turbine", "wind_direction"], + "data": np.full((n_turbines, N_SECTORS), 0.08).tolist(), + } + return resource + + +# --- dict_to_site: joint P construction -------------------------------------- + + +def test_dict_to_site_builds_joint_p(): + resource = _histogram_resource(n_turbines=2) + site = dict_to_site(resource) + + assert "P" in site.ds + assert "probability" not in site.ds + assert "Sector_frequency" not in site.ds + + # Joint = conditional x marginal, and sums to 1 per turbine + cond = np.array(resource["probability"]["data"]) + sector_p = np.array(resource["sector_probability"]["data"]) + # XRSite appends a wd=360 wrap slice for circular interpolation; compare + # on the original sectors. + expected = cond * sector_p[:, :, np.newaxis] + p = site.ds["P"].sel(wd=WD_SECTORS).transpose("i", "wd", "ws").values + npt.assert_array_almost_equal(p, expected, 12) + npt.assert_array_almost_equal(p.sum(axis=(1, 2)), np.ones(2), 12) + + # TI carried into the site dataset (run_simulation prefers it there) + assert "TI" in site.ds.data_vars + + +def test_dict_to_site_bare_probability_is_joint(): + resource = _histogram_resource(n_turbines=2) + cond = np.array(resource["probability"]["data"]) + sector_p = np.array(resource["sector_probability"]["data"]) + joint = cond * sector_p[:, :, np.newaxis] + resource["probability"]["data"] = joint.tolist() + del resource["sector_probability"] + + site = dict_to_site(resource) + assert "P" in site.ds + p = site.ds["P"].sel(wd=WD_SECTORS).transpose("i", "wd", "ws").values + npt.assert_array_almost_equal(p, joint, 12) + + +def test_dict_to_site_windkit_dim_order(): + # windkit exports dims as (wind_speed, wind_direction, wind_turbine); + # the joint-P construction must be dim-order independent. + resource = _histogram_resource(n_turbines=2) + cond = np.array(resource["probability"]["data"]) + resource["probability"] = { + "dims": ["wind_speed", "wind_direction", "wind_turbine"], + "data": cond.transpose(2, 1, 0).tolist(), + } + site = dict_to_site(resource) + + sector_p = np.array(resource["sector_probability"]["data"]) + expected = cond * sector_p[:, :, np.newaxis] + p = site.ds["P"].sel(wd=WD_SECTORS).transpose("i", "wd", "ws").values + npt.assert_array_almost_equal(p, expected, 12) + + +# --- construct_site dispatch -------------------------------------------------- + + +def _system_dat(): + return { + "attributes": {}, + "site": {"boundaries": {"polygons": [{"x": [0, 1000], "y": [0, 1000]}]}}, + } + + +def test_construct_site_dispatches_histogram(): + resource_dat = {"wind_resource": _histogram_resource(n_turbines=2)} + site_data = construct_site(_system_dat(), resource_dat, {"0": 119.0}, [0.0, 500.0]) + + assert site_data["timeseries"] is False + # ws verbatim: pyWake cannot interpolate a ws-dependent P + npt.assert_array_equal(site_data["ws"], WS_CENTERS) + # wd refined to sub-sectors + assert len(site_data["wd"]) == N_SECTORS * 5 + # operating must be (n_turbines, 1) to pass pyWake's arg2ilk validation + assert site_data["operating"].shape == (2, 1) + assert "P" in site_data["site"].ds + + +def test_histogram_ti_default_when_absent(): + resource_dat = {"wind_resource": _histogram_resource(n_turbines=2, with_ti=False)} + site_data = construct_site(_system_dat(), resource_dat, {"0": 119.0}, [0, 500]) + + assert site_data["TI"] == 0.06 + assert "TI" not in site_data["site"].ds.data_vars + + +# --- sub-sector helper -------------------------------------------------------- + + +def test_subsector_wind_directions(): + wd = _subsector_wind_directions(WD_SECTORS, 5) + assert len(wd) == N_SECTORS * 5 + assert np.all(np.diff(wd) > 0) + npt.assert_array_almost_equal(np.diff(wd), np.full(N_SECTORS * 5 - 1, 6.0), 10) + # Sector-0 sub-sectors wrap across north onto a 0-anchored 6° grid + assert wd[0] == pytest.approx(0.0) + assert wd[-1] == pytest.approx(354.0) + + # Trailing 360 wrap value is stripped before refinement + wd_wrapped = _subsector_wind_directions(np.append(WD_SECTORS, 360.0), 5) + npt.assert_array_almost_equal(wd_wrapped, wd, 10) + + # No refinement for n_subsector=1 or fewer than 4 sectors + npt.assert_array_equal(_subsector_wind_directions(WD_SECTORS, 1), WD_SECTORS) + npt.assert_array_equal( + _subsector_wind_directions([0.0, 120.0, 240.0], 5), [0.0, 120.0, 240.0] + ) + + +# --- end-to-end through pyWake ------------------------------------------------ + + +def _run_histogram_simulation(resource, x, y): + """Mirror run_simulation()'s kwargs for a histogram site.""" + resource_dat = {"wind_resource": resource} + site_data = construct_site(_system_dat(), resource_dat, {"0": 119.0}, x) + # No ``operating`` here: that kwarg only exists for WIFA-built turbines + # (PowerCtFunctionList wrapper); the run_pywake test covers it end-to-end. + wfm = BastankhahGaussian(site_data["site"], DTU10MW()) + sim_res = wfm( + x, + y, + type=0, + time=site_data["timeseries"], + ws=site_data["ws"], + wd=site_data["wd"], + yaw=0, + tilt=0, + ) + return sim_res + + +def test_histogram_single_turbine_aep_matches_analytic(): + # wd-uniform resource: the sub-sector interpolation of P is exact, so the + # no-wake AEP must equal sum_ws pdf(ws) * power(ws) * 8760 h / 1e9. + resource = _histogram_resource(n_turbines=1, uniform_wd=True) + sim_res = _run_histogram_simulation(resource, [0.0], [0.0]) + aep = float(sim_res.aep(normalize_probabilities=True).sum()) + + pdf = np.array(resource["probability"]["data"])[0, 0] # same for every sector + power_w = DTU10MW().power(WS_CENTERS) + expected = float((pdf * power_w).sum() * 24 * 365 / 1e9) + npt.assert_array_almost_equal(aep, expected, 6) + + +def test_histogram_two_turbine_wake_loss(): + # Two turbines aligned east-west, 3 diameters apart: with per-turbine P + # (i dim) the simulation runs end-to-end and produces a wake loss. + resource = _histogram_resource(n_turbines=2) + sim_res = _run_histogram_simulation(resource, [0.0, 3 * 178.3], [0.0, 0.0]) + aep_per_turbine = ( + sim_res.aep(normalize_probabilities=True).sum(["ws", "wd"]).to_numpy() + ) + + assert np.all(np.isfinite(aep_per_turbine)) + assert np.all(aep_per_turbine > 0) + total = float(sim_res.aep(normalize_probabilities=True).sum()) + npt.assert_array_almost_equal(total, aep_per_turbine.sum(), 8) + # Wake interaction must cost something relative to two isolated turbines + solo = [ + float( + _run_histogram_simulation( + { + **resource, + "wind_turbine": [0], + "probability": { + "dims": resource["probability"]["dims"], + "data": [resource["probability"]["data"][t]], + }, + "sector_probability": { + "dims": resource["sector_probability"]["dims"], + "data": [resource["sector_probability"]["data"][t]], + }, + "turbulence_intensity": { + "dims": resource["turbulence_intensity"]["dims"], + "data": [resource["turbulence_intensity"]["data"][t]], + }, + }, + [0.0], + [0.0], + ) + .aep(normalize_probabilities=True) + .sum() + ) + for t in range(2) + ] + assert total < sum(solo) + + +def test_run_pywake_histogram_case(tmp_path): + """Full run_pywake() on a copy of the heterogeneous case with a histogram + resource replacing the per-turbine Weibull.""" + src = test_path / "../examples/cases/heterogeneous_wind_rose_at_turbines" + case = tmp_path / "histogram_case" + shutil.copytree(src, case) + + # Swap the flow model to pywake + system_yaml = case / "wind_energy_system/system.yaml" + system_yaml.write_text( + system_yaml.read_text().replace("name: foxes", "name: PyWake") + ) + + # Replace the Weibull resource with a per-turbine histogram, keeping the + # original turbine positions/heights + nc_path = case / "plant_energy_resource/WTResource.nc" + with xr.open_dataset(nc_path) as old: + old_vars = {v: old[v].values.copy() for v in ["x", "y", "height"]} + n_turbines = old.sizes["wind_turbine"] + + resource = _histogram_resource(n_turbines=n_turbines) + ds = xr.Dataset( + { + "probability": ( + tuple(resource["probability"]["dims"]), + np.array(resource["probability"]["data"]), + ), + "sector_probability": ( + tuple(resource["sector_probability"]["dims"]), + np.array(resource["sector_probability"]["data"]), + ), + "turbulence_intensity": ( + tuple(resource["turbulence_intensity"]["dims"]), + np.array(resource["turbulence_intensity"]["data"]), + ), + **{name: (("wind_turbine",), vals) for name, vals in old_vars.items()}, + }, + coords={ + "wind_turbine": np.arange(n_turbines), + "wind_direction": WD_SECTORS, + "wind_speed": WS_CENTERS, + }, + ) + nc_path.unlink() + ds.to_netcdf(nc_path) + + aep = run_pywake(system_yaml, output_dir=str(tmp_path / "output")) + assert np.isfinite(float(aep)) + assert float(aep) > 0 diff --git a/tests/test_pywake_submodels.py b/tests/test_pywake_submodels.py new file mode 100644 index 0000000..657d3a3 --- /dev/null +++ b/tests/test_pywake_submodels.py @@ -0,0 +1,1198 @@ +"""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.deficit_models.utils import ct2a_madsen, ct2a_mom1d +from py_wake.deflection_models import JimenezWakeDeflection +from py_wake.deflection_models.gcl_hill_vortex import GCLHillDeflection +from py_wake.ground_models.ground_models import Mirror +from py_wake.rotor_avg_models import ( + AreaOverlapAvgModel, + CGIRotorAvg, + EqGridRotorAvg, + GaussianOverlapAvgModel, + GQGridRotorAvg, + GridRotorAvg, + PolarGridRotorAvg, + RotorCenter, +) +from py_wake.superposition_models import ( + CumulativeWakeSum, + LinearSum, + MaxSum, + SqrMaxSum, + 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, + _fuga_atmosphere, + _fuga_z0_sweep, + 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, analysis_top=None): + """Helper to call _configure_deficit_model with minimal boilerplate. + + Returns ``(wake_model_class, deficit_args)``; the third element + (``deficit_post_attrs``) is dropped for the common case. ``analysis_top`` + merges keys at the ``analysis`` level (e.g. ``axial_induction_model``). + """ + wind_deficit_model = {"name": name, **(analysis_extra or {})} + analysis = {"wind_deficit_model": wind_deficit_model, **(analysis_top or {})} + cls, args, _post = _configure_deficit_model({"name": name}, analysis, _RD, _HH) + return cls, args + + +def _call_deficit_full(name, analysis_extra=None, analysis_top=None): + """Like ``_call_deficit`` but returns the full 3-tuple including post-attrs.""" + wind_deficit_model = {"name": name, **(analysis_extra or {})} + analysis = {"wind_deficit_model": wind_deficit_model, **(analysis_top or {})} + 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_a(): + """The windIO constant k_a is Bastankhah2014's scalar k; a nonzero TI + coefficient k_b cannot be represented and warns.""" + with pytest.warns(UserWarning, match="k_b=0.38 is ignored"): + _, args = _call_deficit( + "Bastankhah2014", + {"wake_expansion_coefficient": {"k_a": 0.004, "k_b": 0.38}}, + ) + assert args["k"] == 0.004 + + +def test_configure_deficit_model_jensen_k_a_k_b(): + """Jensen: windIO k = k_a + k_b*TI maps to PyWake a = [k_b, k_a].""" + _, args = _call_deficit( + "Jensen", + {"wake_expansion_coefficient": {"k_a": 0.004, "k_b": 0.38}}, + ) + assert args["a"] == [0.38, 0.004] + + +def test_configure_deficit_model_nojlocaldeficit_k_a_k_b(): + """NOJLocalDeficit: windIO k = k_a + k_b*TI maps to PyWake a = [k_b, k_a].""" + _, args = _call_deficit( + "NOJLocalDeficit", + {"wake_expansion_coefficient": {"k_a": 0.004, "k_b": 0.38}}, + ) + 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_jensen_1983_k_a(): + """Jensen_1983 (NOJDeficit) takes its scalar k from windIO's constant k_a, + since the wake_expansion_coefficient schema has no scalar k field.""" + cls, args = _call_deficit( + "Jensen_1983", + {"wake_expansion_coefficient": {"k_a": 0.1, "k_b": 0.0}}, + ) + assert cls is NOJDeficit + assert args["k"] == 0.1 + + +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.004, + "k_b": 0.38, + "free_stream_ti": False, + }, + "ceps": 0.3, + }, + ) + 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, "wake_expansion_coefficient": {"free_stream_ti": True}}, + ) + 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, "free_stream_ti": False}}, + ) + 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 k_a/k_b models.""" + with pytest.warns(UserWarning, match="uses k_a/k_b"): + _, 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_b(): + """Verify warning when only the constant k_a is given to a TI-dependent + expansion model (k_b defaults to 0).""" + with pytest.warns(UserWarning, match="k_b not specified"): + _, args = _call_deficit( + "Zong2020", {"wake_expansion_coefficient": {"k_a": 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.004, + "k_b": 0.38, + "free_stream_ti": False, + }, + "ceps": 0.3, + }, + NiayifarGaussianDeficit, + ), + ( + "Zong2020", + { + "wake_expansion_coefficient": { + "k_a": 0.004, + "k_b": 0.38, + "free_stream_ti": True, + } + }, + ZongGaussianDeficit, + ), + ( + "Jensen", + {"wake_expansion_coefficient": {"k_a": 0.004, "k_b": 0.38}}, + NOJLocalDeficit, + ), + ( + "NOJLocalDeficit", + {"wake_expansion_coefficient": {"k_a": 0.004, "k_b": 0.38}}, + 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") + + +# --------------------------------------------------------------------------- +# TI reference flag (windIO free_stream_ti -> PyWake use_effective_ti) +# --------------------------------------------------------------------------- + + +@pytest.mark.parametrize( + "name", + [ + "Jensen", + "Niayifar2016", + "Carbajofuertes2018", + "Zong2020", + "TurbOPark", + "SuperGaussian", + "SuperGaussian2023", + ], +) +@pytest.mark.parametrize("free_stream_ti,expected", [(False, True), (True, False)]) +def test_free_stream_ti_inverts_to_use_effective_ti(name, free_stream_ti, expected): + """free_stream_ti maps to use_effective_ti with inverted polarity, for every + TI-capable deficit (including SuperGaussian and TurbOPark, which the previous + narrow handling missed).""" + _, args = _call_deficit( + name, {"wake_expansion_coefficient": {"free_stream_ti": free_stream_ti}} + ) + assert args["use_effective_ti"] is expected + # kwargs must actually instantiate the model + cls, _ = _call_deficit(name) + cls(**args) + + +@pytest.mark.parametrize("name", ["Bastankhah2014"]) +def test_free_stream_ti_ignored_for_non_ti_capable(name): + """Deficits without a use_effective_ti param must not receive it, even if + free_stream_ti is present (would raise TypeError on instantiation). + + (GCL was moved to TI_CAPABLE — GCLDeficit accepts use_effective_ti, as + GCLLocal demonstrates.)""" + _, args = _call_deficit( + name, {"wake_expansion_coefficient": {"k_b": 0.04, "free_stream_ti": True}} + ) + assert "use_effective_ti" not in args + + +def test_use_effective_ti_key_is_ignored(): + """The deprecated top-level use_effective_ti key is no longer read.""" + _, args = _call_deficit( + "Niayifar2016", + { + "wake_expansion_coefficient": {"k_a": 0.38, "k_b": 0.004}, + "use_effective_ti": False, + }, + ) + assert "use_effective_ti" not in args + + +# --------------------------------------------------------------------------- +# 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_vector_not_implemented(): + """Vector superposition is foxes-only; the pyWake path rejects it.""" + with pytest.raises(NotImplementedError, match="Vector"): + _configure_superposition_model({"ws_superposition": "Vector"}) + + +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), + ("grid", GridRotorAvg), + ("avg_deficit", GridRotorAvg), + ("Avg_Deficit", GridRotorAvg), + ("gaussian_overlap", GaussianOverlapAvgModel), + ("gaussianoverlap", GaussianOverlapAvgModel), + ("area_overlap", AreaOverlapAvgModel), + ("areaoverlap", AreaOverlapAvgModel), + ("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"}, {}) + + +@pytest.mark.parametrize( + "name,expected_class", + [ + ("SelfSimilarityDeficit2020", SelfSimilarityDeficit2020), + ("SelfSimilarityDeficit", SelfSimilarityDeficit), + ("RankineHalfBody", RankineHalfBody), + ("Rathmann", Rathmann), + ("VortexCylinder", VortexCylinder), + ("VortexDipole", VortexDipole), + ("HybridInduction", HybridInduction), + ], +) +def test_configure_blockage_model_ground_mirror(name, expected_class): + model = _configure_blockage_model({"name": name, "ground_mirror": True}, {}) + assert isinstance(model, expected_class) + assert isinstance(model.groundModel, Mirror) + + +def test_configure_blockage_model_no_ground_mirror_by_default(): + model = _configure_blockage_model({"name": "SelfSimilarityDeficit2020"}, {}) + assert model.groundModel is None + + model = _configure_blockage_model( + {"name": "SelfSimilarityDeficit2020", "ground_mirror": False}, {} + ) + assert model.groundModel is None + + +def test_configure_blockage_model_ground_mirror_fuga_warns(tmp_path): + # FUGA LUTs already include the ground; the flag must warn and be ignored. + # A real LUT is not needed to hit the warning path, but FugaDeficit + # requires a valid path, so only the warning is asserted before the + # constructor fails on the dummy path. + with pytest.warns(UserWarning, match="ground_mirror is ignored for FUGA"): + try: + _configure_blockage_model( + {"name": "FUGA", "ground_mirror": True}, + {"LUT_path": str(tmp_path / "missing.nc")}, + ) + except Exception: + pass + + +def test_configure_blockage_model_ground_mirror_increases_deficit(): + # Physics sanity check: the image rotor must increase the upstream + # blockage deficit relative to the unmirrored model. + import numpy as np + from py_wake.deficit_models.no_wake import NoWakeDeficit + from py_wake.examples.data.hornsrev1 import V80 + from py_wake.site._site import UniformSite + from py_wake.wind_farm_models.engineering_models import All2AllIterative + + site = UniformSite([1], ti=0.06) + wt = V80() + + def upstream_ws(cfg): + wfm = All2AllIterative( + site, + wt, + NoWakeDeficit(), + blockage_deficitModel=_configure_blockage_model(cfg, {}), + ) + # second turbine acts as a probe 3D upstream (wd=270 blows along +x) + sim_res = wfm([0, -240], [0, 0], wd=270, ws=10) + return sim_res.WS_eff.sel(wt=1).item() + + ws_plain = upstream_ws({"name": "SelfSimilarityDeficit2020"}) + ws_mirror = upstream_ws( + {"name": "SelfSimilarityDeficit2020", "ground_mirror": True} + ) + assert ws_mirror < ws_plain < 10.0 + assert np.isclose(ws_plain, ws_mirror, atol=0.1) # small correction, same order + + +def test_ground_mirror_works_with_squared_wake_superposition(): + # Regression: Mirror without an explicit superposition falls back to the + # wind farm model's superposition. With SquaredSum (e.g. the TurbOPark + # recipe) that asserts on the speed-up (negative deficit) regions every + # blockage model produces ("SquaredSum only works for deficit - not + # speedups"). The mirror must therefore sum real+image linearly. + from py_wake.deficit_models.utils import ct2a_madsen + from py_wake.examples.data.hornsrev1 import V80 + from py_wake.site._site import UniformSite + from py_wake.superposition_models import SquaredSum + from py_wake.turbulence_models.stf import STF2017TurbulenceModel + from py_wake.wind_farm_models.engineering_models import All2AllIterative + from py_wake.deficit_models.gaussian import TurboGaussianDeficit + + site = UniformSite([1], ti=0.06) + wt = V80() + wfm = All2AllIterative( + site, + wt, + TurboGaussianDeficit(ct2a=ct2a_madsen), + superpositionModel=SquaredSum(), + turbulenceModel=STF2017TurbulenceModel(), + blockage_deficitModel=_configure_blockage_model( + {"name": "VortexCylinder", "ground_mirror": True}, {} + ), + ) + sim_res = wfm([0, 400], [0, 0], wd=270, ws=10) + assert (sim_res.WS_eff < 10.0).all() + + +# --------------------------------------------------------------------------- +# 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 + + +def _weighted_system(rotor_name, ws_superposition="Weighted"): + return { + "attributes": { + "analysis": { + "wind_deficit_model": { + "name": "Zong2020", + "wake_expansion_coefficient": {"k_a": 0.38, "k_b": 0.004}, + }, + "superposition_model": {"ws_superposition": ws_superposition}, + "rotor_averaging": {"name": rotor_name}, + "deflection_model": {"name": "None"}, + "turbulence_model": {"name": "None"}, + "blockage_model": {"name": None}, + } + } + } + + +@pytest.mark.parametrize("rotor_name", ["gaussian_overlap", "area_overlap", "center"]) +@pytest.mark.parametrize("ws_superposition", ["Weighted", "Cumulative"]) +def test_weighted_superposition_requires_node_rotor_avg(rotor_name, ws_superposition): + """WeightedSum/CumulativeWakeSum with a non-node rotor-averaging model raises + a clear ValueError instead of PyWake's deep AssertionError.""" + with pytest.raises(ValueError, match="node"): + configure_wake_model( + _weighted_system(rotor_name, ws_superposition), + rotor_diameter=126.0, + hub_height=90.0, + ) + + +@pytest.mark.parametrize("rotor_name", ["grid", "eq_grid", "gq_grid", "cgi"]) +def test_weighted_superposition_allows_node_rotor_avg(rotor_name): + """A node rotor-averaging model is accepted with Weighted superposition.""" + config = configure_wake_model( + _weighted_system(rotor_name), rotor_diameter=126.0, hub_height=90.0 + ) + assert isinstance(config["superposition_model"], WeightedSum) + + +def _weighted_deficit_system(deficit_name): + """Weighted superposition + node rotor-avg, varying only the deficit.""" + return { + "attributes": { + "analysis": { + "wind_deficit_model": { + "name": deficit_name, + "wake_expansion_coefficient": {"free_stream_ti": True}, + }, + "superposition_model": {"ws_superposition": "Weighted"}, + "rotor_averaging": {"name": "grid"}, + "deflection_model": {"name": "None"}, + "turbulence_model": {"name": "None"}, + "blockage_model": {"name": None}, + } + } + } + + +@pytest.mark.parametrize("deficit_name", ["SuperGaussian", "SuperGaussian2023", "GCL"]) +def test_weighted_superposition_requires_convection_deficit(deficit_name): + """WeightedSum/CumulativeWakeSum with a non-ConvectionDeficitModel deficit + (super-Gaussian, GCL) raises a clear ValueError instead of PyWake's deep + AssertionError, even when the rotor-averaging model is a node model.""" + with pytest.raises(ValueError, match="ConvectionDeficitModel"): + configure_wake_model( + _weighted_deficit_system(deficit_name), + rotor_diameter=126.0, + hub_height=90.0, + ) + + +@pytest.mark.parametrize("deficit_name", ["Zong2020", "Niayifar2016", "Bastankhah2014"]) +def test_weighted_superposition_allows_convection_deficit(deficit_name): + """Convection-based deficits are accepted with Weighted superposition.""" + config = configure_wake_model( + _weighted_deficit_system(deficit_name), rotor_diameter=126.0, hub_height=90.0 + ) + assert isinstance(config["superposition_model"], WeightedSum) + + +# --------------------------------------------------------------------------- +# CrespoHernandez calibration coefficients (Phase 2 / Fix C) +# --------------------------------------------------------------------------- + + +def test_crespo_default_without_c(): + """No c -> the PyWake-default CrespoHernandez (ct2a_madsen).""" + tm = _configure_turbulence_model({"name": "CrespoHernandez"}) + assert isinstance(tm, CrespoHernandez) + assert tm.ct2a is ct2a_madsen + + +def test_crespo_with_c_uses_literature_recipe(): + """c -> CrespoHernandez with those coefficients, 1D induction and SqrMaxSum.""" + c = [0.73, 0.83, 0.03, -0.32] + tm = _configure_turbulence_model({"name": "CrespoHernandez", "c": c}) + assert isinstance(tm, CrespoHernandez) + assert list(tm.c) == c + assert tm.ct2a is ct2a_mom1d + assert isinstance(tm.addedTurbulenceSuperpositionModel, SqrMaxSum) + + +# --------------------------------------------------------------------------- +# 'none' rotor averaging + WeightedSum (Phase 2 / Fix D) +# --------------------------------------------------------------------------- + + +def test_rotor_averaging_none(): + assert _configure_rotor_averaging({"name": "none"}) is None + + +def test_weighted_superposition_allows_none_rotor(): + """WeightedSum accepts rotorAvgModel=None (rotor centre), as Zong (2020) uses.""" + sd = _weighted_deficit_system("Zong2020") + sd["attributes"]["analysis"]["rotor_averaging"] = {"name": "none"} + config = configure_wake_model(sd, rotor_diameter=126.0, hub_height=90.0) + assert config["rotor_averaging"] is None + assert isinstance(config["superposition_model"], WeightedSum) + + +def test_weighted_superposition_rejects_center_rotor(): + """A non-node, non-None rotor model is still rejected for WeightedSum.""" + sd = _weighted_deficit_system("Zong2020") + sd["attributes"]["analysis"]["rotor_averaging"] = {"name": "center"} + with pytest.raises(ValueError, match="node"): + configure_wake_model(sd, rotor_diameter=126.0, hub_height=90.0) + + +# --------------------------------------------------------------------------- +# Zong ceps -> eps_coeff (Phase 2) +# --------------------------------------------------------------------------- + + +def test_zong_ceps_maps_to_eps_coeff(): + """Zong's near-wake epsilon is named eps_coeff (not ceps) in PyWake.""" + cls, args = _call_deficit( + "Zong2020", + {"wake_expansion_coefficient": {"k_a": 0.38, "k_b": 0.004}, "ceps": 0.35}, + ) + assert cls is ZongGaussianDeficit + assert args["eps_coeff"] == 0.35 + assert "ceps" not in args + + +# --------------------------------------------------------------------------- +# Axial induction -> ct2a (honoring axial_induction_model) +# --------------------------------------------------------------------------- + + +@pytest.mark.parametrize( + "axial, expected", + [("1D", ct2a_mom1d), ("Madsen", ct2a_madsen), ("madsen", ct2a_madsen)], +) +def test_axial_induction_sets_ct2a(axial, expected): + """axial_induction_model maps to the deficit's ct2a on a ct2a-capable model.""" + _, args = _call_deficit( + "Bastankhah2014", analysis_top={"axial_induction_model": axial} + ) + assert args["ct2a"] is expected + + +def test_axial_induction_default_absent_keeps_deficit_default(): + """No axial_induction_model -> no ct2a override (deficit keeps its default).""" + _, args = _call_deficit("Bastankhah2014") + assert "ct2a" not in args + + +def test_axial_induction_skipped_when_deficit_has_no_ct2a(): + """Blondel2020 has no ct2a parameter, so it is left untouched.""" + _, args = _call_deficit( + "SuperGaussian", analysis_top={"axial_induction_model": "1D"} + ) + assert "ct2a" not in args + + +def test_axial_induction_ct2a_instantiates(): + """The injected ct2a is a valid constructor argument.""" + cls, args = _call_deficit( + "Niayifar2016", analysis_top={"axial_induction_model": "1D"} + ) + assert isinstance(cls(**args), NiayifarGaussianDeficit) + + +# --------------------------------------------------------------------------- +# TurbOPark canonical recipe (Nygaard 2022) +# --------------------------------------------------------------------------- + + +def test_turbopark_recipe_ground_and_ctlim(): + """TurbOPark gets a Mirror ground model and ctlim=0.96 as constructor args.""" + cls, args, post = _call_deficit_full("TurbOPark") + assert cls is TurboGaussianDeficit + assert isinstance(args["groundModel"], Mirror) + assert args["ctlim"] == 0.96 + + +def test_turbopark_recipe_ws_key_post_attr(): + """TurbOPark scales by the downstream ambient WS via WS_key='WS_jlk'.""" + _, _, post = _call_deficit_full("TurbOPark") + assert post == {"WS_key": "WS_jlk"} + + +def test_turbopark_recipe_instantiates_and_applies_ws_key(): + """The recipe builds a valid deficit and WS_key is set post-construction.""" + cls, args, post = _call_deficit_full("TurbOPark") + deficit = cls(**args) + for attr, value in post.items(): + setattr(deficit, attr, value) + assert deficit.WS_key == "WS_jlk" + + +def test_non_turbopark_has_no_post_attrs(): + """Only TurbOPark carries post-construction attributes.""" + _, _, post = _call_deficit_full("Bastankhah2014") + assert post == {} + + +# --------------------------------------------------------------------------- +# use_effective_ws / use_effective_ti for GCL (Phase 3) +# --------------------------------------------------------------------------- + + +def test_use_effective_ws_honored(): + """The windIO use_effective_ws flag is passed through (not hardcoded).""" + _, args = _call_deficit("Bastankhah2014", {"use_effective_ws": False}) + assert args["use_effective_ws"] is False + + +def test_use_effective_ws_defaults_true(): + _, args = _call_deficit("Bastankhah2014") + assert args["use_effective_ws"] is True + + +def test_gcl_honors_free_stream_ti(): + """GCLDeficit accepts use_effective_ti (GCLLocal); free_stream_ti is honored.""" + _, args = _call_deficit( + "GCL", {"wake_expansion_coefficient": {"free_stream_ti": False}} + ) + assert args["use_effective_ti"] is True + + +# --------------------------------------------------------------------------- +# Fuga LUT-atmosphere helpers (pure functions; no LUT generation) +# --------------------------------------------------------------------------- +import numpy as np # noqa: E402 + + +def _resource(ti=None, abl=None, z0=None): + wr = {} + if ti is not None: + wr["turbulence_intensity"] = {"data": np.asarray(ti, dtype=float)} + if abl is not None: + wr["ABL_height"] = {"data": np.asarray(abl, dtype=float)} + if z0 is not None: + wr["z0"] = {"data": np.asarray(z0, dtype=float)} + return {"wind_resource": wr} + + +def test_fuga_atmosphere_derives_z0_from_ti(): + # Neutral inversion: z0 = zhub * exp(-1/TI). Mean TI 0.10 at hub 78. + z0, zi, zeta0, ti = _fuga_atmosphere( + _resource(ti=np.full(50, 0.10), abl=np.full(50, 600.0)), {}, 78.0 + ) + assert zeta0 == 0.0 + assert zi == 600.0 # from ABL_height + assert ti == pytest.approx(0.10) + assert z0 == pytest.approx(78.0 * np.exp(-1 / 0.10), rel=1e-3) + + +def test_fuga_atmosphere_defaults_without_resource(): + z0, zi, zeta0, ti = _fuga_atmosphere(None, {}, 78.0) + assert z0 == 0.03 and zi == 500.0 and ti is None + + +def test_fuga_atmosphere_config_overrides_site(): + z0, zi, _, _ = _fuga_atmosphere( + _resource(ti=np.full(10, 0.10)), {"z0": 0.001, "zi": 800.0}, 78.0 + ) + assert z0 == 0.001 and zi == 800.0 + + +def test_fuga_z0_sweep_spans_distribution_and_is_sorted(): + rng = np.random.default_rng(0) + ti = np.clip(rng.normal(0.10, 0.03, 5000), 0.03, 0.25) + z0s = _fuga_z0_sweep(_resource(ti=ti), {"n_z0": 5}, 78.0, 0.0, 0.03) + assert len(z0s) >= 2 + assert z0s == sorted(z0s) + # all physical: clamp keeps z0 within ~[1e-5, 0.3] m + assert z0s[0] >= 1e-5 and z0s[-1] <= 0.31 + + +def test_fuga_z0_sweep_clamps_high_ti_to_physical_z0(): + # All-high TI must not produce absurd roughness (TI 0.30 -> z0 ~2.8 m). + z0s = _fuga_z0_sweep(_resource(ti=np.full(200, 0.30)), {"n_z0": 5}, 78.0, 0.0, 0.5) + assert max(z0s) <= 0.31 + + +def test_fuga_z0_sweep_single_when_n1_or_no_ti(): + assert _fuga_z0_sweep( + _resource(ti=np.full(10, 0.1)), {"n_z0": 1}, 78.0, 0.0, 0.04 + ) == [0.04] + assert _fuga_z0_sweep(None, {"n_z0": 5}, 78.0, 0.0, 0.04) == [0.04] + + +def test_fuga_z0_sweep_explicit_z0_overrides(): + z0s = _fuga_z0_sweep( + _resource(ti=np.full(10, 0.1)), {"z0": [0.01, 0.05]}, 78.0, 0.0, 0.04 + ) + assert z0s == [0.01, 0.05] + + +def test_fuga_atmosphere_ignores_list_z0(): + # An explicit z0 list is a sweep (handled by _fuga_z0_sweep); _fuga_atmosphere + # must not choke on it when computing the scalar fallback. + z0, zi, _, _ = _fuga_atmosphere( + _resource(ti=np.full(10, 0.10)), {"z0": [1e-4, 1e-2]}, 78.0 + ) + assert isinstance(z0, float) # derived from TI, not the list + + +# --------------------------------------------------------------------------- +# Fuga mixed turbine geometries (LUT generation intercepted; no pyfuga) +# +# Regression for mixed-rotor-diameter/hub-height layouts (neighbor farms): +# single-hub-height LUTs at different hub heights turn FugaMultiLUTDeficit's +# merged table all-NaN (xarray cannot interpolate a size-1 z axis), which +# surfaced as zero production at every waked turbine. Mixed layouts must get +# LUTs spanning all hub heights, one shared x/y grid, and a z_lst pinning the +# merged z axis to the hub heights. +# --------------------------------------------------------------------------- +import wifa.pywake_api as _pywake_api # noqa: E402 + + +def _call_fuga(monkeypatch, turbine_geometries=None): + """Call the FUGA deficit branch with _ensure_fuga_luts intercepted.""" + captured = {} + + def fake_ensure_luts(**kwargs): + captured.update(kwargs) + n = len(kwargs["geometries"]) * len(kwargs["z0_list"]) + return [f"/fake/lut_{i}.nc" for i in range(n)] + + monkeypatch.setattr(_pywake_api, "_ensure_fuga_luts", fake_ensure_luts) + analysis = { + "wind_deficit_model": {"name": "FUGA", "fuga": {"z0": 0.03, "n_z0": 1}} + } + _, args, _ = _configure_deficit_model( + {"name": "FUGA"}, analysis, _RD, _HH, turbine_geometries=turbine_geometries + ) + return args, captured + + +def test_fuga_mixed_geometries_span_hub_heights_and_share_grid(monkeypatch): + args, luts = _call_fuga(monkeypatch, [(80.0, 80.0), (100.0, 90.0)]) + # Every LUT spans all hub heights so the merged z axis has no NaN slices. + assert luts["zlow"] == 80.0 and luts["zhigh"] == 90.0 + # One shared x/y grid (finest natural resolution: min diameter / 4, / 16). + assert luts["dx"] == 20.0 and luts["dy"] == 5.0 + # The merged z axis is pinned to exactly the hub heights. + assert args["z_lst"] == [80.0, 90.0] + assert isinstance(args["LUT_path"], list) and len(args["LUT_path"]) == 2 + + +def test_fuga_mixed_diameters_same_hub_share_grid_only(monkeypatch): + args, luts = _call_fuga(monkeypatch, [(80.0, 80.0), (100.0, 80.0)]) + assert luts["dx"] == 20.0 and luts["dy"] == 5.0 + # Same hub height everywhere: single-level LUTs stay valid, no z axis. + assert "zlow" not in luts and "z_lst" not in args + + +def test_fuga_single_geometry_keeps_hub_level_lut(monkeypatch): + args, luts = _call_fuga(monkeypatch, None) + # Unchanged cheap path: per-geometry defaults, single plain LUT path. + assert "zlow" not in luts and "dx" not in luts + assert "z_lst" not in args + assert isinstance(args["LUT_path"], str) + + +def test_fuga_duplicate_geometries_deduped(monkeypatch): + # Two turbine types sharing one geometry must not produce duplicate d_h + # coordinates (they break FugaMultiLUTDeficit's merge) nor a second LUT. + args, luts = _call_fuga(monkeypatch, [(80.0, 80.0), (80.0, 80.0)]) + assert luts["geometries"] == [(80.0, 80.0)] + assert isinstance(args["LUT_path"], str) diff --git a/tests/test_wayve.py b/tests/test_wayve.py index 4957f04..43c46d5 100644 --- a/tests/test_wayve.py +++ b/tests/test_wayve.py @@ -27,5 +27,856 @@ def test_wayve_4wts(): run_wayve(yaml_input, output_dir=output_dir_name, debug_mode=True) +def test_wayve_k_convention_lanzilao(): + """windIO defines k = k_a + k_b*TI; wayve's Lanzilao computes + kwake = ka*TI + kb, so the reader must map swapped: ka <- k_b, kb <- k_a. + Regression for the unswapped pass-through that silently turned the + constant expansion into a TI multiplier (and vice versa).""" + from wifa.wayve_api import wake_model_setup + + wm = wake_model_setup( + { + "wind_deficit_model": { + "wake_expansion_coefficient": {"k_a": 0.004, "k_b": 0.38}, + "ceps": 0.2, + } + } + ) + assert wm.ka == 0.38 # Lanzilao TI multiplier holds windIO k_b + assert wm.kb == 0.004 # Lanzilao constant holds windIO k_a + + +def test_wayve_scalar_k_is_constant_expansion(): + """A scalar k is a constant expansion: Lanzilao kb, with no TI term.""" + from wifa.wayve_api import wake_model_setup + + wm = wake_model_setup( + { + "wind_deficit_model": { + "wake_expansion_coefficient": {"k": 0.05}, + "ceps": 0.2, + } + } + ) + assert wm.ka == 0.0 + assert wm.kb == 0.05 + + +_LANZILAO_ANALYSIS = { + "wind_deficit_model": { + "wake_expansion_coefficient": {"k_a": 0.004, "k_b": 0.38}, + "ceps": 0.2, + } +} + + +def test_wake_tool_read_from_wm_coupling(): + """The windIO schema nests `wake_tool` under `wm_coupling` (a top-level + `analysis.wake_tool` fails validation). Reading only the analysis level made + the foxes coupling unreachable from any schema-valid file.""" + from wifa.wayve_api import wake_model_setup + + analysis = {**_LANZILAO_ANALYSIS, "wm_coupling": {"wake_tool": "wayve"}} + wm = wake_model_setup(analysis) + assert type(wm).__name__ == "Lanzilao" + + analysis = {**_LANZILAO_ANALYSIS, "wm_coupling": {"wake_tool": "nonsense"}} + with pytest.raises(NotImplementedError, match="nonsense"): + wake_model_setup(analysis) + + +def test_wake_tool_analysis_level_fallback(): + """Legacy hand-written yamls put `wake_tool` at the analysis level.""" + from wifa.wayve_api import wake_model_setup + + with pytest.raises(NotImplementedError, match="nonsense"): + wake_model_setup({**_LANZILAO_ANALYSIS, "wake_tool": "nonsense"}) + + +def test_wake_tool_defaults_to_wayve(): + from wifa.wayve_api import wake_model_setup + + wm = wake_model_setup(dict(_LANZILAO_ANALYSIS)) + assert type(wm).__name__ == "Lanzilao" + + +def _scalar_resource(ti): + """Minimal single-point (no `height`) wind resource for flow_io_abl.""" + return { + "wind_speed": {"data": [9.0]}, + "wind_direction": {"data": [270.0]}, + "turbulence_intensity": {"data": [ti]}, + "z0": {"data": [0.03]}, + } + + +def test_scalar_resource_ti_is_a_fraction(): + """windIO carries TI as a fraction. The scalar branch used to divide it by + 100 (and read it only when `z0` happened to be present), collapsing the + TI-proportional wake expansion by two orders of magnitude.""" + from wifa.wayve_api import flow_io_abl + + abl, _ = flow_io_abl(_scalar_resource(0.08), 0, zh=78.0, h1=156.0) + assert abl.TI == pytest.approx(0.08) + + +def test_scalar_resource_ti_read_without_z0(): + """TI must be read whenever it is present, not only when `z0` is.""" + from wifa.wayve_api import flow_io_abl + + resource = _scalar_resource(0.08) + del resource["z0"] + abl, _ = flow_io_abl(resource, 0, zh=78.0, h1=156.0) + assert abl.TI == pytest.approx(0.08) + + +def test_scalar_resource_ti_defaults_when_absent(): + from wifa.wayve_api import flow_io_abl + + resource = _scalar_resource(0.08) + del resource["turbulence_intensity"] + abl, _ = flow_io_abl(resource, 0, zh=78.0, h1=156.0) + assert abl.TI == pytest.approx(0.04) + + +def test_scalar_resource_air_density_is_per_state(): + """A per-state `air_density` in the wind resource sets abl.rho, which + linearly scales the reported turbine power (the APM velocity solution + does not depend on it).""" + from wifa.wayve_api import flow_io_abl + + resource = _scalar_resource(0.08) + resource["air_density"] = {"data": [1.29, 1.11]} + resource["wind_speed"]["data"].append(9.0) + resource["wind_direction"]["data"].append(270.0) + resource["turbulence_intensity"]["data"].append(0.08) + resource["z0"]["data"].append(0.03) + assert flow_io_abl(resource, 0, zh=78.0, h1=156.0)[0].rho == pytest.approx(1.29) + assert flow_io_abl(resource, 1, zh=78.0, h1=156.0)[0].rho == pytest.approx(1.11) + + +def test_scalar_resource_air_density_defaults_when_absent(): + from wifa.wayve_api import flow_io_abl + + abl, _ = flow_io_abl(_scalar_resource(0.08), 0, zh=78.0, h1=156.0) + assert abl.rho == pytest.approx(1.225) + + +def _timeseries_system(wind_speeds, subset): + """Single-turbine system whose wind speed differs at every timestep.""" + n = len(wind_speeds) + return { + "site": { + "energy_resource": { + "wind_resource": { + "time": list(range(n)), + "wind_speed": {"data": list(wind_speeds)}, + "wind_direction": {"data": [270.0] * n}, + "turbulence_intensity": {"data": [0.08] * n}, + "z0": {"data": [0.03] * n}, + "fc": {"data": [1.0e-4] * n}, + } + } + }, + "wind_farm": { + "layouts": [{"coordinates": {"x": [0.0], "y": [0.0]}}], + "turbines": { + "performance": { + "power_curve": { + "power_values": [0, 1e6, 3e6, 3e6], + "power_wind_speeds": [3, 8, 12, 25], + }, + "Ct_curve": { + "Ct_values": [0.8, 0.8, 0.4, 0.2], + "Ct_wind_speeds": [3, 8, 12, 25], + }, + }, + "hub_height": 80.0, + "rotor_diameter": 80.0, + }, + }, + "attributes": { + "flow_model": {"name": "wayve"}, + "analysis": { + "wind_deficit_model": { + "wake_expansion_coefficient": {"k_a": 0.04, "k_b": 0.0}, + "ceps": 0.2, + }, + "superposition_model": {"ws_superposition": "Product"}, + "wm_coupling": {"method": "PB"}, + }, + "model_outputs_specification": { + "output_folder": "output", + "run_configuration": { + "times_run": {"all_occurences": False, "subset": subset} + }, + "turbine_outputs": { + "turbine_nc_filename": "turbine_data.nc", + "output_variables": ["rotor_effective_velocity"], + }, + }, + }, + } + + +def test_times_run_subset_selects_the_requested_rows(tmp_path): + """Regression: the subsetted timestamps were re-enumerated from 0, so the + ABL was built from rows 0..n-1 of the wind resource while being *labelled* + with the requested timestamps — silently simulating the wrong states.""" + import xarray as xr + + from wifa.wayve_api import run_wayve + + wind_speeds = [6.0, 9.0, 12.0] + subset = [2] + # debug_mode skips the APM solve: the wake model still sees the full ABL. + run_wayve( + _timeseries_system(wind_speeds, subset), + output_dir=tmp_path, + debug_mode=True, + ) + + ds = xr.load_dataset(tmp_path / "turbine_data.nc") + assert ds["states"].values.tolist() == subset + # A lone turbine is unwaked, so its rotor-effective velocity is the + # background speed of the *requested* state (12 m/s), not of row 0 (6 m/s). + rews = float(ds["rotor_effective_velocity"].isel(states=0, turbine=0)) + assert rews == pytest.approx(wind_speeds[subset[0]], rel=0.02) + + if __name__ == "__main__": test_wayve_4wts() + + +# --------------------------------------------------------------------------- +# Solver-frame rotation (westerly normalization) and profile-resolving +# free atmosphere +# --------------------------------------------------------------------------- + +_CT_CURVE = { + "Ct_values": [0.8, 0.8, 0.4, 0.2], + "Ct_wind_speeds": [3, 8, 12, 25], +} +_POWER_CURVE = { + "power_values": [0, 1e6, 3e6, 3e6], + "power_wind_speeds": [3, 8, 12, 25], +} + + +def _solver_system(x, y, wind_resource, flow_field=None, layers_description=None): + """Multi-turbine system for full (non-debug) APM solves on a small grid.""" + analysis = { + "wind_deficit_model": { + "wake_expansion_coefficient": {"k_a": 0.04, "k_b": 0.0}, + "ceps": 0.2, + }, + "superposition_model": {"ws_superposition": "Product"}, + "wm_coupling": {"method": "PB"}, + "apm_grid": {"Lx": 1.0e5, "Ly": 1.0e5, "dx": 1.0e3}, + } + if layers_description is not None: + analysis["layers_description"] = layers_description + outputs = { + "output_folder": "output", + "run_configuration": {"times_run": {"all_occurences": True}}, + "turbine_outputs": { + "turbine_nc_filename": "turbine_data.nc", + "output_variables": ["power", "rotor_effective_velocity"], + }, + } + if flow_field is not None: + outputs["flow_field"] = flow_field + return { + "site": {"energy_resource": {"wind_resource": wind_resource}}, + "wind_farm": { + "layouts": [{"coordinates": {"x": list(x), "y": list(y)}}], + "turbines": { + "performance": {"power_curve": _POWER_CURVE, "Ct_curve": _CT_CURVE}, + "hub_height": 80.0, + "rotor_diameter": 80.0, + }, + }, + "attributes": { + "flow_model": {"name": "wayve"}, + "analysis": analysis, + "model_outputs_specification": outputs, + }, + } + + +def _scalar_wind_resource(wd, n=1): + return { + "time": list(range(n)), + "wind_speed": {"data": [9.0] * n}, + "wind_direction": {"data": [float(wd)] * n}, + "turbulence_intensity": {"data": [0.08] * n}, + "z0": {"data": [0.03] * n}, + "fc": {"data": [1.0e-4] * n}, + } + + +def _flow_field_spec(x_bounds, y_bounds, n): + return { + "report": True, + "flow_nc_filename": "flow_field.nc", + "output_variables": ["wind_speed", "wind_direction"], + "z_planes": { + "xy_sampling": "grid", + "x_bounds": list(x_bounds), + "y_bounds": list(y_bounds), + "Nx": n[0], + "Ny": n[1], + "z_sampling": "hub_heights", + }, + } + + +def test_scalar_resource_is_solver_frame_aligned(): + """The ABL comes out with the hub-height wind along +x (wayve's westerly + convention) for any input direction, and the undone rotation is + reported.""" + import numpy as np + + from wifa.wayve_api import flow_io_abl + + resource = _scalar_resource(0.08) + resource["wind_direction"]["data"] = [225.0] + abl, rotation = flow_io_abl(resource, 0, zh=78.0, h1=156.0) + assert rotation == pytest.approx(np.deg2rad(45.0)) + u_hub = np.interp(78.0, abl.zs, abl.us) + v_hub = np.interp(78.0, abl.zs, abl.vs) + assert u_hub > 0.0 + assert v_hub == pytest.approx(0.0, abs=1.0e-8 * u_hub) + + +@pytest.fixture(scope="module") +def _rotation_pair(tmp_path_factory): + """Physically identical cases under rotation: a westerly-aligned turbine + row (a); the same row and grids rotated 90 degrees with a southerly wind, + run as two identical states (b); and the row rotated by a generic 37 + degrees (c). 90 degrees maps the square FFT grid onto itself, so only the + generic angle pins the solver-frame rotation itself (see the tests).""" + import numpy as np + import xarray as xr + + from wifa.wayve_api import run_wayve + + span = 400.0 + xa = [0.0, span, 2 * span] + ff_a = _flow_field_spec((-1000.0, 1800.0), (-1200.0, 1200.0), (15, 13)) + ff_b = _flow_field_spec((-1200.0, 1200.0), (-1000.0, 1800.0), (13, 15)) + system_a = _solver_system(xa, [0.0, 0.0, 0.0], _scalar_wind_resource(270.0), ff_a) + system_b = _solver_system( + [0.0, 0.0, 0.0], xa, _scalar_wind_resource(180.0, n=2), ff_b + ) + # Case c: rotate the westerly row into the frame of a 233 deg (270 - 37) + # wind: earth layout = R(+rotation) applied to the westerly layout. + rot_c = np.deg2rad(270.0 - 233.0) + c, s = np.cos(rot_c), np.sin(rot_c) + xc = [c * x for x in xa] + yc = [s * x for x in xa] + system_c = _solver_system(xc, yc, _scalar_wind_resource(233.0)) + out_a = tmp_path_factory.mktemp("wayve_rot_a") + out_b = tmp_path_factory.mktemp("wayve_rot_b") + out_c = tmp_path_factory.mktemp("wayve_rot_c") + run_wayve(system_a, output_dir=out_a) + run_wayve(system_b, output_dir=out_b) + run_wayve(system_c, output_dir=out_c) + return { + "turbines_a": xr.load_dataset(out_a / "turbine_data.nc"), + "turbines_b": xr.load_dataset(out_b / "turbine_data.nc"), + "turbines_c": xr.load_dataset(out_c / "turbine_data.nc"), + "flow_a": xr.load_dataset(out_a / "flow_field.nc"), + "flow_b": xr.load_dataset(out_b / "flow_field.nc"), + } + + +def test_rotational_invariance_of_turbine_outputs(_rotation_pair): + """A southerly case with the layout turned 90 degrees is the same physical + system as the westerly case, so per-turbine outputs must match. Before the + solver-frame rotation, the gravity-wave solve saw the farm footprint (and + any anisotropic grid) at the wrong angle for non-westerly states.""" + import numpy as np + + power_a = _rotation_pair["turbines_a"]["power"].values + power_b = _rotation_pair["turbines_b"]["power"].values + assert np.all(np.isfinite(power_a)) and np.all(power_a > 0.0) + np.testing.assert_allclose(power_b[0], power_a[0], rtol=1.0e-6) + rews_a = _rotation_pair["turbines_a"]["rotor_effective_velocity"].values + rews_b = _rotation_pair["turbines_b"]["rotor_effective_velocity"].values + np.testing.assert_allclose(rews_b[0], rews_a[0], rtol=1.0e-6) + # Downstream turbines are waked in both frames + assert power_a[0, 0] > power_a[0, -1] + + +def test_rotational_invariance_at_generic_angle(_rotation_pair): + """A generic 37-degree rotation is what actually pins the solver-frame + normalization: 90-degree rotations map the square FFT grid onto itself, + so the pre-rotation code already passed the 90-degree pair exactly, while + it fails this case at ~4e-4 relative. At a generic angle both cos and sin + terms of every rotation formula are exercised (90 deg blinds the cos + terms, 45 deg hides cos/sin swaps).""" + import numpy as np + + power_a = _rotation_pair["turbines_a"]["power"].values + power_c = _rotation_pair["turbines_c"]["power"].values + np.testing.assert_allclose(power_c, power_a, rtol=1.0e-6) + + +def test_repeated_states_give_identical_output(_rotation_pair): + """Two identical southerly states must give identical per-turbine output. + Guards the per-state wind-farm rebuild: a cumulative in-place rotation of + farm_dat would rotate state 1 twice yet still produce finite, positive, + plausibly-waked power.""" + import numpy as np + + power = _rotation_pair["turbines_b"]["power"] + np.testing.assert_allclose( + power.isel(states=1).values, power.isel(states=0).values, rtol=1.0e-12 + ) + + +def test_flow_field_is_reported_in_the_earth_frame(_rotation_pair): + """flow_field.nc stays on the requested east/north grid: the southerly + case reports ~180 deg background flow, with the wake north (downstream) + of the farm.""" + import numpy as np + + ds = _rotation_pair["flow_b"].isel(states=0, z=0) + # Upstream (south) edge, away from the farm axis: background direction + wd_upstream = float( + ds["wind_direction"].sel(x=-1200.0, y=-1000.0, method="nearest") + ) + assert wd_upstream % 360.0 == pytest.approx(180.0, abs=3.0) + # The deepest deficit sits downstream of the last turbine (north) + ws = ds["wind_speed"] + loc = ws.argmin(dim=["x", "y"]) + assert float(ws.y[int(loc["y"])]) >= 800.0 + assert abs(float(ws.x[int(loc["x"])])) < 400.0 + + +def test_flow_fields_map_onto_each_other_under_rotation(_rotation_pair): + """The southerly flow field equals the westerly one rotated by 90 deg: + ws_b(x, y) == ws_a(y, -x) on the matching grids.""" + import numpy as np + + ws_a = _rotation_pair["flow_a"]["wind_speed"].isel(states=0, z=0).values + ws_b = _rotation_pair["flow_b"]["wind_speed"].isel(states=0, z=0).values + # grids: a is (15 x, 13 y), b is (13 x, 15 y); earth point (x, y) in b + # maps to solver/earth point (y, -x) in a; x_b symmetric -> index flip + mapped = np.transpose(ws_a)[::-1, :] + np.testing.assert_allclose(ws_b, mapped, rtol=1.0e-6) + wd_a = _rotation_pair["flow_a"]["wind_direction"].isel(states=0, z=0).values + wd_b = _rotation_pair["flow_b"]["wind_direction"].isel(states=0, z=0).values + mapped_wd = (np.transpose(wd_a)[::-1, :] - 90.0) % 360.0 + np.testing.assert_allclose(wd_b % 360.0, mapped_wd, rtol=0, atol=1.0e-6) + + +def test_xy_points_matches_xy_plane(): + """_xy_plane_points on an axis-aligned meshgrid reproduces wayve's + xy_plane exactly (it mirrors the same code with the meshgrid hoisted + out) — the guard against drift between the fork helper and wayve.""" + import numpy as np + from wayve.apm import APM + from wayve.grid.grid import Stat2Dgrid + from wayve.momentum_flux_parametrizations import FrictionCoefficients + from wayve.solvers import FixedPointIteration + + from wifa.wayve_api import ( + _pressure_for_state, + _xy_plane_points, + flow_io_abl, + wf_setup, + ) + + system = _solver_system( + [0.0, 400.0, 800.0], [0.0, 0.0, 0.0], _scalar_wind_resource(270.0) + ) + farm_dat = system["wind_farm"] + analysis_dat = system["attributes"]["analysis"] + resource = system["site"]["energy_resource"]["wind_resource"] + abl, rotation = flow_io_abl(resource, 0, zh=80.0, h1=160.0) + assert rotation == pytest.approx(0.0) + wind_farm, forcing, _, _ = wf_setup( + farm_dat, analysis_dat, 1.0e3, rotation=rotation + ) + grid = Stat2Dgrid(1.0e5, 100, 1.0e5, 100) + model = APM(grid, forcing, abl, FrictionCoefficients(), _pressure_for_state(abl, 1)) + model.solve(method=FixedPointIteration(tol=5.0e-3, relax=0.7)) + coupling = wind_farm.coupling + wake_model = coupling.wake_model + u_bg_evaluator = coupling.set_up_u_bg_evaluator(abl) + apm_evaluator = coupling.apm_evaluator + xs = np.linspace(-2000.0, 4000.0, 25) + ys = np.linspace(-1500.0, 1500.0, 21) + ref = wake_model.xy_plane( + wind_farm, abl, u_bg_evaluator, apm_evaluator, xs, ys, 80.0 + ) + Xs, Ys = np.meshgrid(xs, ys, indexing="ij") + new = _xy_plane_points( + wake_model, wind_farm, abl, u_bg_evaluator, apm_evaluator, Xs, Ys, 80.0 + ) + for ref_field, new_field in zip(ref, new): + np.testing.assert_array_equal(new_field, ref_field) + + +def _profile_wind_resource(wd_hub=225.0, veer=10.0, n=1, spelling="flat"): + """Truncated (1.7 km) mesoscale-style profile resource: winds/theta + profiles plus ERA5-style surface scalars and an explicit CI block. + + ``spelling`` selects how the capping inversion is written: ``"flat"`` + (the windIO schema keys, which is what a netCDF ``!include`` can carry) + or ``"nested"`` (the older thermal_stratification block). The two must + produce identical ABL states — see + ``test_flat_and_nested_capping_inversion_agree``.""" + import numpy as np + + zs = np.linspace(10.0, 1700.0, 18) + ws = 8.0 + 3.0 * zs / 1700.0 + # Small linear veer across the profile around the hub-height direction + hub_frac = np.interp(80.0, zs, np.linspace(0.0, 1.0, len(zs))) + wd = wd_hub + veer * (np.linspace(0.0, 1.0, len(zs)) - hub_frac) + theta = np.where( + zs < 1000.0, + 290.0, + np.where( + zs <= 1200.0, + 290.0 + 3.0 * (zs - 1000.0) / 200.0, + 293.0 + 3.5e-3 * (zs - 1200.0), + ), + ) + per_state = lambda arr: {"data": [list(arr)] * n} # noqa: E731 + scalar = lambda val: {"data": [val] * n} # noqa: E731 + if spelling == "flat": + ci = { + "ABL_height": scalar(1100.0), + "capping_inversion_thickness": scalar(200.0), + "capping_inversion_strength": scalar(3.0), + "lapse_rate": scalar(3.5e-3), + } + elif spelling == "nested": + ci = { + "thermal_stratification": { + "capping_inversion": { + "ABL_height": scalar(1100.0), + "dH": scalar(200.0), + "dtheta": scalar(3.0), + "lapse_rate": scalar(3.5e-3), + } + } + } + else: + raise ValueError(f"unknown spelling {spelling!r}") + return { + "time": list(range(n)), + "height": list(zs), + "wind_speed": per_state(ws), + "wind_direction": per_state(wd), + "potential_temperature": per_state(theta), + "turbulence_intensity": scalar(0.08), + "air_density": scalar(1.2), + "z0": scalar(0.03), + "fc": scalar(1.0e-4), + "friction_velocity": scalar(0.35), + "boundary_layer_height": scalar(900.0), + "LMO": scalar(5.0e3), + **ci, + } + + +def test_flat_and_nested_capping_inversion_agree(): + """The two windIO spellings of the capping inversion must be the same input. + + The flat keys are the windIO schema spelling (and what code_saturne's + adapter already reads); the nested thermal_stratification block is the + older one, kept so existing files keep working. A netCDF ``!include`` + can only carry the flat one, so the whole point is that switching a file + over changes nothing about the atmosphere it describes. + """ + import numpy as np + + from wifa.wayve_api import flow_io_abl + + flat, _ = flow_io_abl(_profile_wind_resource(spelling="flat"), 0, zh=80.0, h1=160.0) + nested, _ = flow_io_abl( + _profile_wind_resource(spelling="nested"), 0, zh=80.0, h1=160.0 + ) + # The four scalars drive exactly these: the upper-layer depth, the + # reduced gravity, the free-atmosphere stratification and the inversion + # bounds. Exact equality — the values travel through unchanged. + for field in ("H", "H2", "gprime", "N", "inv_bottom", "inv_top", "U1", "U2", "U3"): + assert getattr(flat, field) == getattr(nested, field), field + for field in ("zs", "us", "vs", "tauxs", "tauys"): + np.testing.assert_array_equal( + getattr(flat, field), getattr(nested, field), err_msg=field + ) + + +def test_capping_inversion_rejects_both_spellings_at_once(): + """A file carrying both spellings is half-migrated; picking one silently + would pick a physics answer for the user.""" + from wifa.wayve_api import capping_inversion_spelling + + resource = _profile_wind_resource(spelling="flat") + resource |= { + k: v + for k, v in _profile_wind_resource(spelling="nested").items() + if k == "thermal_stratification" + } + with pytest.raises(ValueError, match="twice"): + capping_inversion_spelling(resource) + + +@pytest.mark.parametrize( + "spelling, dropped", + [ + ("flat", "capping_inversion_strength"), + ("flat", "ABL_height"), + ("nested", "dtheta"), + ], +) +def test_capping_inversion_rejects_a_partial_block(spelling, dropped): + """A partial block is not completed from defaults. + + dtheta alone sets the amplitude of the gravity-wave forcing, so pairing a + real ABL_height with an invented one yields a result that looks fitted and + is not. (WIFA's code_saturne adapter does fill per-key defaults; that + behaviour is not copied here deliberately.) + """ + from wifa.wayve_api import capping_inversion_spelling + + resource = _profile_wind_resource(spelling=spelling) + if spelling == "flat": + del resource[dropped] + else: + del resource["thermal_stratification"]["capping_inversion"][dropped] + with pytest.raises(ValueError, match="Incomplete"): + capping_inversion_spelling(resource) + + +def test_no_capping_inversion_falls_back_to_the_profile_fit(): + """Omitting all four scalars is a valid, distinct request: fit the + inversion from the potential-temperature profile.""" + from wifa.wayve_api import capping_inversion_spelling, flow_io_abl + + resource = _profile_wind_resource(spelling="flat") + for key in ( + "ABL_height", + "capping_inversion_thickness", + "capping_inversion_strength", + "lapse_rate", + ): + del resource[key] + assert capping_inversion_spelling(resource) is None + + fitted, _ = flow_io_abl(resource, 0, zh=80.0, h1=160.0) + explicit, _ = flow_io_abl( + _profile_wind_resource(spelling="flat"), 0, zh=80.0, h1=160.0 + ) + # The fit finds the 1000-1200 m jump in the test profile, not the + # explicit block's 1100 +/- 100 m; if these agreed the test would be moot. + assert fitted.H != explicit.H + + +def test_truncated_profile_free_atmosphere_capped_at_data_top(): + """The free atmosphere of a truncated profile is capped at the top of the + data: NonUniform's layers then span actual profile points instead of the + default 10 km tropopause (constant-fill splines above ~1.7 km).""" + import numpy as np + + from wifa.wayve_api import flow_io_abl + + abl, rotation = flow_io_abl(_profile_wind_resource(), 0, zh=80.0, h1=160.0) + assert abl.h_strat == pytest.approx(1700.0) + assert abl.Uinf == pytest.approx(abl.us[-1]) + assert abl.Vinf == pytest.approx(abl.vs[-1]) + assert abl.Ninf == pytest.approx(abl.N) + # Solver frame: hub-height wind along +x, veer preserved aloft. + # (Interpolating components vs. directions differs by ~1e-3 m/s here.) + u_hub = np.interp(80.0, abl.zs, abl.us) + v_hub = np.interp(80.0, abl.zs, abl.vs) + assert u_hub > 0.0 + assert v_hub == pytest.approx(0.0, abs=1.0e-3) + assert rotation == pytest.approx(np.deg2rad(270.0 - 225.0), abs=0.02) + assert np.abs(abl.vs).max() > 0.1 # veer survives the rotation + + +def test_pressure_for_state_selects_nonuniform_only_with_data(): + """NonUniform needs profile points between the inversion top and the + (capped) free-atmosphere top; scalar states fall back to Uniform.""" + from wifa.wayve_api import _pressure_for_state, flow_io_abl + + abl_profile, _ = flow_io_abl(_profile_wind_resource(), 0, zh=80.0, h1=160.0) + assert type(_pressure_for_state(abl_profile, 6)).__name__ == "NonUniform" + assert type(_pressure_for_state(abl_profile, 1)).__name__ == "Uniform" + # Scalar branch: the Nieuwstadt profile stops at the inversion -> no + # data-backed free atmosphere -> warn and fall back + abl_scalar, _ = flow_io_abl(_scalar_resource(0.08), 0, zh=78.0, h1=156.0) + with pytest.warns(UserWarning, match="falling back to the Uniform"): + assert type(_pressure_for_state(abl_scalar, 6)).__name__ == "Uniform" + + +def test_run_wayve_profile_nonuniform_end_to_end(tmp_path): + """Non-westerly mesoscale profiles with a profile-resolving free + atmosphere run end-to-end and produce finite, waked turbine output — + and NonUniform is genuinely active (no silent Uniform fallback), so the + number_of_fa_layers wiring through run_wayve is pinned, not just the + _pressure_for_state unit behavior.""" + import warnings as _warnings + + import numpy as np + import xarray as xr + + from wifa.wayve_api import run_wayve + + # Row aligned with the 225 deg flow (SW): downstream is to the NE + span = 400.0 / np.sqrt(2.0) + system = _solver_system( + [0.0, span, 2 * span], + [0.0, span, 2 * span], + _profile_wind_resource(wd_hub=225.0), + layers_description={"farm_layer_height": 160.0, "number_of_fa_layers": 6}, + ) + with _warnings.catch_warnings(record=True) as caught: + _warnings.simplefilter("always") + run_wayve(system, output_dir=tmp_path) + fallbacks = [w for w in caught if "falling back to the Uniform" in str(w.message)] + assert not fallbacks, "NonUniform silently fell back to Uniform" + ds = xr.load_dataset(tmp_path / "turbine_data.nc") + power = ds["power"].values + assert power.shape == (1, 3) + assert np.all(np.isfinite(power)) and np.all(power > 0.0) + # SW flow: the NE-most turbine is waked + assert power[0, 0] > power[0, -1] + + +# Full windIO-shaped analysis block for the foxes wake tool (mirrors what the +# flow_model_chain Experiment serializer emits; foxes' windIO reader needs the +# named model blocks, which the Lanzilao path simply ignores). +_FOXES_ANALYSIS = { + "wind_deficit_model": { + "name": "Bastankhah2014", + "wake_expansion_coefficient": {"k_a": 0.04, "k_b": 0.0}, + "use_effective_ws": True, + "ceps": 0.2, + }, + "axial_induction_model": "Madsen", + "deflection_model": {"name": "None"}, + "turbulence_model": {"name": "None"}, + "superposition_model": { + "ws_superposition": "Product", + "ti_superposition": "Squared", + }, + "rotor_averaging": { + "name": "none", + "background_averaging": "center", + "wake_averaging": "center", + "wind_speed_exponent_for_power": 3, + "wind_speed_exponent_for_ct": 2, + }, + "blockage_model": {"name": "None"}, + "wm_coupling": {"method": "PB", "wake_tool": "foxes"}, + "apm_grid": {"Lx": 1.0e5, "Ly": 1.0e5, "dx": 1.0e3}, +} + + +def test_xy_points_matches_xy_plane_foxes(): + """Foxes twin of the _xy_plane_points drift guard. The solve itself also + exercises the background_flow_direction shim: without it, wayve@e87780a's + broken e_spanwise import crashes every foxes-coupled solve.""" + pytest.importorskip("foxes") + import numpy as np + from wayve.apm import APM + from wayve.grid.grid import Stat2Dgrid + from wayve.momentum_flux_parametrizations import FrictionCoefficients + from wayve.solvers import FixedPointIteration + + from wifa.wayve_api import ( + _pressure_for_state, + _xy_plane_points, + flow_io_abl, + wf_setup, + ) + + system = _solver_system( + [0.0, 400.0, 800.0], [0.0, 0.0, 0.0], _scalar_wind_resource(270.0) + ) + farm_dat = system["wind_farm"] + resource = system["site"]["energy_resource"]["wind_resource"] + abl, rotation = flow_io_abl(resource, 0, zh=80.0, h1=160.0) + wind_farm, forcing, _, _ = wf_setup( + farm_dat, _FOXES_ANALYSIS, 1.0e3, rotation=rotation + ) + wake_model = wind_farm.coupling.wake_model + assert type(wake_model).__name__ == "FoxesWakeModel" + grid = Stat2Dgrid(1.0e5, 100, 1.0e5, 100) + model = APM(grid, forcing, abl, FrictionCoefficients(), _pressure_for_state(abl, 1)) + model.solve(method=FixedPointIteration(tol=5.0e-3, relax=0.7)) + coupling = wind_farm.coupling + u_bg_evaluator = coupling.set_up_u_bg_evaluator(abl) + apm_evaluator = coupling.apm_evaluator + xs = np.linspace(-2000.0, 4000.0, 13) + ys = np.linspace(-1500.0, 1500.0, 11) + ref = wake_model.xy_plane( + wind_farm, abl, u_bg_evaluator, apm_evaluator, xs, ys, 80.0 + ) + Xs, Ys = np.meshgrid(xs, ys, indexing="ij") + new = _xy_plane_points( + wake_model, wind_farm, abl, u_bg_evaluator, apm_evaluator, Xs, Ys, 80.0 + ) + for ref_field, new_field in zip(ref, new): + np.testing.assert_array_equal(new_field, ref_field) + + +def test_turbulence_profile_stress_rotation(): + """Explicit tau_x/tau_y stress profiles are earth-frame vectors: they must + rotate into the solver frame with the velocities (magnitude preserved).""" + import numpy as np + + from wifa.wayve_api import flow_io_abl + + zs = np.linspace(10.0, 1700.0, 18) + # Jet-like profile: wayve's ABL estimates the eddy viscosity from where + # the shear turns negative, so the speed must peak below the profile top + ws = 8.0 + 3.0 * np.minimum(zs, 1100.0) / 1100.0 + ws -= 0.5 * np.maximum(0.0, zs - 1100.0) / 600.0 + wd = np.full(zs.size, 200.0) + theta = 290.0 + 3.5e-3 * zs + # Stress magnitude decaying to zero at 800 m, at a fixed earth angle + tau_mag = 0.12 * np.maximum(0.0, 1.0 - zs / 800.0) ** 1.5 + phi_tau = np.deg2rad(30.0) + tau_x_in = tau_mag * np.cos(phi_tau) + tau_y_in = tau_mag * np.sin(phi_tau) + scalar = lambda val: {"data": [val]} # noqa: E731 + resource = { + "time": [0], + "height": list(zs), + "wind_speed": {"data": [list(ws)]}, + "wind_direction": {"data": [list(wd)]}, + "potential_temperature": {"data": [list(theta)]}, + "tau_x": {"data": [list(tau_x_in)]}, + "tau_y": {"data": [list(tau_y_in)]}, + "turbulence_intensity": scalar(0.08), + "z0": scalar(0.03), + "fc": scalar(1.0e-4), + "thermal_stratification": { + "capping_inversion": { + "ABL_height": scalar(1100.0), + "dH": scalar(200.0), + "dtheta": scalar(3.0), + "lapse_rate": scalar(3.5e-3), + } + }, + } + abl, rotation = flow_io_abl(resource, 0, zh=80.0, h1=160.0) + assert rotation == pytest.approx(np.deg2rad(270.0 - 200.0)) + c, s = np.cos(rotation), np.sin(rotation) + np.testing.assert_allclose( + abl.tauxs, c * tau_x_in + s * tau_y_in, rtol=1e-12, atol=1e-15 + ) + np.testing.assert_allclose( + abl.tauys, c * tau_y_in - s * tau_x_in, rtol=1e-12, atol=1e-15 + ) + np.testing.assert_allclose( + np.hypot(abl.tauxs, abl.tauys), tau_mag, rtol=1e-12, atol=1e-15 + ) diff --git a/wifa/foxes_api.py b/wifa/foxes_api.py index 911f6fb..af160f3 100644 --- a/wifa/foxes_api.py +++ b/wifa/foxes_api.py @@ -66,7 +66,9 @@ def run_foxes( idir = input_dir else: input_yaml = Path(input_yaml) - wio = load_yaml(input_yaml) + # Keep an included wind_resource.nc as numpy arrays (foxes' reader uses + # ndarrays directly), avoiding the dict-of-lists memory blow-up. + wio = load_yaml(input_yaml, nc_data="array") idir = input_yaml.parent idict, algo, odir = read_windio_dict(wio, verbosity=verbosity) diff --git a/wifa/pywake_api.py b/wifa/pywake_api.py index 91cf99a..20c8000 100644 --- a/wifa/pywake_api.py +++ b/wifa/pywake_api.py @@ -5,15 +5,19 @@ 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", @@ -33,7 +37,11 @@ "rotor_averaging": { "name": "Center", }, - "blockage_model": {"name": None, "ss_alpha": 0.8888888888888888}, + "blockage_model": { + "name": None, + "ss_alpha": 0.8888888888888888, + "ground_mirror": False, + }, } @@ -43,16 +51,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"): @@ -69,8 +86,11 @@ def load_and_validate_config(yaml_input, default_output_dir="output"): from windIO import validate as validate_yaml if not isinstance(yaml_input, dict): - validate_yaml(yaml_input, "plant/wind_energy_system") - system_dat = load_yaml(Path(yaml_input)) + # Keep an included wind_resource.nc as numpy arrays instead of exploding + # it into Python lists; validate structure-only so jsonschema does not + # require/iterate the bulk data. + validate_yaml(yaml_input, "plant/wind_energy_system", array_data=True) + system_dat = load_yaml(Path(yaml_input), nc_data="array") else: system_dat = yaml_input @@ -100,8 +120,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 @@ -109,9 +131,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 = [] @@ -150,11 +170,25 @@ 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) + 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"], 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, ) @@ -193,6 +227,23 @@ def dict_to_site(resource_dict): from windIO import dict_to_netcdf resource_ds = dict_to_netcdf(resource_dict) + + if "probability" in resource_ds: + # windIO histogram resource. Producers that follow the windkit/WAsP + # convention write ``probability`` as the per-sector *conditional* + # distribution over wind speed (sums to 1 per wind_direction slice) + # with the direction marginal in ``sector_probability``; pyWake's + # ``P`` is the *joint* distribution over (wd, ws). Multiply the two + # when both are present; a bare ``probability`` is taken as already + # joint (matching py_wake.utils.ieawind37_utils). + if "sector_probability" in resource_ds: + resource_ds["P"] = ( + resource_ds["probability"] * resource_ds["sector_probability"] + ) + resource_ds = resource_ds.drop_vars(["probability", "sector_probability"]) + else: + resource_ds = resource_ds.rename({"probability": "P"}) + rename_map = { "height": "h", "weibull_a": "Weibull_A", @@ -231,7 +282,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) @@ -267,68 +317,51 @@ 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 - + elif "probability" in wind_resource: + result = _construct_histogram_site(resource_dat, hub_heights, x_positions) 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. @@ -343,21 +376,32 @@ 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", {}) + ) + do_subset = False + if not times_run.get("all_occurences", True) and "subset" in times_run: + cases_idx = times_run["subset"] + do_subset = True heights = wind_resource.get("height") + # ``heights`` may be a list or a numpy array (array-backed resource), so use + # an explicit boolean rather than the ambiguous truth value of an array. + has_heights = heights is not None and len(heights) > 0 - # Helper to get data and dimensions safely + def _subset(arr): + # Only copy when a real time subset is requested. An all-True boolean + # index still copies the whole array, so skip it in the common case. + return arr[cases_idx] if do_subset else arr + + # Helper to get data and dimensions safely. np.asarray avoids copying when + # the data is already an ndarray (e.g. an array-backed windIO resource). def get_resource_data(var_name): data_obj = wind_resource[var_name] - vals = np.array(data_obj["data"]) + vals = np.asarray(data_obj["data"]) dims = data_obj.get("dims", ["time"]) return vals, dims @@ -366,8 +410,8 @@ def get_resource_data(var_name): wd_vals, wd_dims = get_resource_data("wind_direction") # Apply subsetting - ws_vals = ws_vals[cases_idx] - wd_vals = wd_vals[cases_idx] + ws_vals = _subset(ws_vals) + wd_vals = _subset(wd_vals) # Prepare reference arrays - average across turbines if turbine-specific if "wind_turbine" in ws_dims: @@ -386,7 +430,7 @@ def get_resource_data(var_name): # Handle operating status if "operating" in wind_resource: - operating = np.array(wind_resource["operating"]["data"])[cases_idx].T + operating = _subset(np.asarray(wind_resource["operating"]["data"])).T assert operating.shape[0] == len(x_positions) else: operating = np.ones((len(x_positions), len(cases_idx))) @@ -397,69 +441,118 @@ def get_resource_data(var_name): site = None if len(hub_heights) > 1: - # Multiple turbine types - need height interpolation - flow_field_spec = ( - system_dat["attributes"] - .get("model_outputs_specification", {}) - .get("flow_field", {}) - ) - if ( - "z_planes" in flow_field_spec - and flow_field_spec["z_planes"] != "hub_heights" - ): - additional_heights = flow_field_spec.get("z_planes", {}).get("z_list", []) + # Multiple turbine types with differing hub heights. + if ("wind_turbine" in ws_dims or "wind_turbine" in wd_dims) and not has_heights: + # Per-turbine ws/wd are given without a vertical profile (e.g. terrain + # speedups produced by a microscale model, already at each turbine's + # hub height). Preserve the per-turbine values via dict_to_site rather + # than averaging across turbines or interpolating a vertical profile. + site = dict_to_site(wind_resource) + if "turbulence_intensity" not in wind_resource: + TI = 0.02 + else: + TI = _subset(np.asarray(wind_resource["turbulence_intensity"]["data"])) + elif "wind_turbine" in ws_dims or "wind_turbine" in wd_dims: + # Both a vertical profile and per-turbine data is ambiguous — we + # cannot tell whether height or turbine indexes the inflow. + raise NotImplementedError( + "A wind resource with both a 'height' profile and a " + "'wind_turbine' dimension is not supported. Provide one or the " + "other for mixed hub heights." + ) + else: + # Vertical wind profile (a `height` dimension) with mixed hub + # heights: interpolate the profile to each distinct hub height (plus + # any requested flow-field z-planes) and build a height-indexed + # XRSite. pyWake then assigns every turbine its inflow at its own + # hub height. + flow_field_spec = ( + system_dat["attributes"] + .get("model_outputs_specification", {}) + .get("flow_field", {}) + ) + if ( + "z_planes" in flow_field_spec + and flow_field_spec["z_planes"] != "hub_heights" + ): + additional_heights = flow_field_spec.get("z_planes", {}).get( + "z_list", [] + ) - speeds, dirs, TIs, seen = [], [], [], [] - for hh in sorted(np.append(list(hub_heights.values()), additional_heights)): - if hh in seen: - continue - seen.append(hh) - ws_int, wd_int = _interpolate_wind_data(heights, ws, wd, hh) - speeds.append(ws_int) - dirs.append(wd_int) + h_levels = sorted(set(hub_heights.values()) | set(additional_heights)) - ws, wd = ws_int, wd_int + # Interpolate WS (linear) and WD (vector / circular) to each level. + speeds = [_interpolate_wind_data(heights, ws, wd, h)[0] for h in h_levels] + dirs = [_interpolate_wind_dir(heights, wd, h) for h in h_levels] + n_time = np.asarray(speeds[0]).shape[0] - # Handle TI interpolation - if "turbulence_intensity" not in wind_resource: - TI = 0.02 - else: - TI_data = np.array(wind_resource["turbulence_intensity"]["data"])[cases_idx] - for hh in sorted(np.append(list(hub_heights.values()), additional_heights)): - if hh in seen[len(speeds) :]: - continue - if heights: - ti_int = _interpolate_with_min(heights, TI_data, hh, min_val=0.02) - else: - ti_int = TI_data - TIs.append(ti_int) - TI = ti_int + # Interpolate TI to each level (floored), else a constant default. + if "turbulence_intensity" not in wind_resource: + ti_levels = [np.full(n_time, 0.02) for _ in h_levels] + else: + ti_data = _subset( + np.asarray(wind_resource["turbulence_intensity"]["data"]) + ) + ti_levels = [ + _interpolate_with_min(heights, ti_data, h, min_val=0.02) + for h in h_levels + ] data_vars = { "WS": (["h", "time"], np.array(speeds)), "WD": (["h", "time"], np.array(dirs)), - "TI": (["h", "time"], np.array(TIs)), - "P": 1, + "TI": (["h", "time"], np.array(ti_levels)), + "P": (["time"], np.ones(n_time) / n_time), } if "density" in wind_resource: - density_vals = np.array(wind_resource["density"]["data"])[cases_idx] + density_vals = _subset(np.asarray(wind_resource["density"]["data"])) density_dims = wind_resource["density"].get("dims", ["time"]) - if "wind_turbine" in density_dims: - density_vals = np.mean(density_vals, axis=1) - data_vars["Air_density"] = (["time"], density_vals) + if "height" in density_dims: + data_vars["Air_density"] = ( + ["h", "time"], + np.array( + [ + _interpolate_with_min(heights, density_vals, h, 0.0) + for h in h_levels + ] + ), + ) + elif "wind_turbine" in density_dims: + data_vars["Air_density"] = (["time"], np.mean(density_vals, axis=1)) + else: + data_vars["Air_density"] = (["time"], density_vals) + site = XRSite( xr.Dataset( data_vars=data_vars, - coords={"h": seen, "time": np.arange(len(times))}, + coords={"h": h_levels, "time": np.arange(n_time)}, ) ) + + # Reference inflow arrays (1-D over time). The height-indexed site + # is authoritative per turbine, so these only define the flow cases; + # use the turbine-averaged inflow for a representative reference. + turbine_types = system_dat["wind_farm"]["layouts"][0]["turbine_types"] + ordered_hh = list(hub_heights.values()) + level_of = {h: i for i, h in enumerate(h_levels)} + idx_per_turbine = [level_of[ordered_hh[t]] for t in turbine_types] + ws = np.mean([speeds[i] for i in idx_per_turbine], axis=0) + rads = np.deg2rad([dirs[i] for i in idx_per_turbine]) + wd = np.mod( + np.rad2deg( + np.arctan2( + np.mean(np.sin(rads), axis=0), np.mean(np.cos(rads), axis=0) + ) + ), + 360.0, + ) + TI = np.mean([ti_levels[i] for i in idx_per_turbine], axis=0) else: # Single turbine type - print(np.array(ws).shape, np.array(heights).shape) - if heights: + if has_heights: ws, wd = _interpolate_wind_data(heights, ws, wd, hh) - assert len(np.array(times)[cases_idx]) == len(ws) + assert len(_subset(np.asarray(times))) == len(ws) assert len(wd) == len(ws) if "wind_turbine" in ws_dims or "wind_turbine" in wd_dims: @@ -467,15 +560,15 @@ def get_resource_data(var_name): else: site = Hornsrev1Site() if "density" in wind_resource: - density_vals = np.array(wind_resource["density"]["data"])[cases_idx] + density_vals = _subset(np.asarray(wind_resource["density"]["data"])) site.ds["Air_density"] = (("time",), density_vals) # Handle TI if "turbulence_intensity" not in wind_resource: TI = 0.02 else: - TI = np.array(wind_resource["turbulence_intensity"]["data"])[cases_idx] - if heights: + TI = _subset(np.asarray(wind_resource["turbulence_intensity"]["data"])) + if has_heights: TI = interp1d(heights, TI, axis=1)(hh) return { @@ -490,27 +583,126 @@ def get_resource_data(var_name): } -def _construct_weibull_site(resource_dat, hub_heights, x_positions): - """Construct site from Weibull distribution data. +def _coord_values(coord): + """Extract a windIO coordinate as a float array. + + windIO permits coordinates either as plain arrays or as + ``{"data": ..., "dims": ...}`` mappings. + """ + if isinstance(coord, dict): + coord = coord["data"] + return np.asarray(coord, dtype=float) + + +def _subsector_wind_directions(wd_raw, n_subsector): + """Split each wind-direction sector into ``n_subsector`` sub-directions. + + Strips a trailing 360° wrap-around value first. Refinement requires at + least 4 sectors (equidistant, full circle); otherwise the sector centers + are returned unchanged. + """ + 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, + ) + return np.sort( + (wd_sectors[:, np.newaxis] + offsets[np.newaxis, :]).ravel() % 360 + ) + return wd_sectors + + +def _construct_histogram_site(resource_dat, hub_heights, x_positions, n_subsector=5): + """Construct site from a (wind_direction × wind_speed) histogram resource. Internal helper for construct_site(). + + The windIO resource carries ``probability`` on a shared wind-speed bin + grid; dict_to_site() builds the joint pyWake ``P`` from it (multiplying + in ``sector_probability`` when present). Flow cases: pyWake cannot + interpolate a ws-dependent ``P`` onto foreign wind speeds, so the + histogram bin centers are used verbatim; wind directions are refined to + sub-sectors (as in the Weibull path) — XRSite interpolates ``P`` + circularly over wd and rescales by the sub-sector width. + + 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. Default 5 + (matching the Weibull path). """ - from windIO import dict_to_netcdf + wind_resource = resource_dat["wind_resource"] + + ws = _coord_values(wind_resource["wind_speed"]) + wd = _subsector_wind_directions( + _coord_values(wind_resource["wind_direction"]), n_subsector + ) + + if "turbulence_intensity" in wind_resource: + # dict_to_site carries the resource TI into site.ds, which + # run_simulation prefers over this entry; kept for interface parity + # with the Weibull path. + TI = wind_resource["turbulence_intensity"]["data"] + else: + TI = 0.06 # default when TI is absent from wind resource + + return { + "site": dict_to_site(wind_resource), + "ws": ws, + "wd": wd, + "TI": TI, + "timeseries": False, + "operating": np.ones((len(x_positions), 1)), + "additional_heights": [], + "cases_idx": np.ones(1).astype(bool), + } + + +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). + """ 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 +715,59 @@ 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 + # Start at the first nonzero bin: a ws=0 reference case carries zero + # energy for every model, but it is a degenerate flow case for the + # WeightedSum superposition (Zong), whose convection-velocity iteration + # divides by the convection speed and is undefined at zero wind speed. + # Including ws=0 silently corrupts the WeightedSum AEP (collapsing the + # apparent wake loss); dropping it is harmless for all other models. + ws = np.arange(0.5, np.ceil(ws_max_ref) + 0.5, 0.5) + + # -- Wind direction sub-sectors --------------------------------------- + wd = _subsector_wind_directions(wd_raw, n_subsector) + 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: + # Reuse the dataset already materialized inside ``site`` instead of + # rebuilding it with a second dict_to_netcdf. dict_to_site renames + # turbulence_intensity -> TI and height -> h. + ti_da = site.ds["TI"] + if "x" in ti_da.dims: + interpolated_ti = ti_da.interp(x=x_positions, y=x_positions) + if "h" in interpolated_ti.dims: + interpolated_ti = interpolated_ti.interp(h=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, @@ -573,6 +803,33 @@ def _interpolate_wind_data(heights, ws, wd, target_height): return ws_int, wd_int +def _interpolate_wind_dir(heights, wd, target_height): + """Interpolate wind direction (degrees) to ``target_height``. + + Uses vector (sin/cos) interpolation so the 0/360 wrap-around is handled + correctly — plain linear interpolation of the raw degrees would average, + e.g., 350° and 10° to 180° instead of 0°. + """ + if heights is None: + return wd + rad = np.deg2rad(np.asarray(wd, dtype=float)) + try: + sin_i = interp1d(heights, np.sin(rad), axis=1, fill_value="extrapolate")( + target_height + ) + cos_i = interp1d(heights, np.cos(rad), axis=1, fill_value="extrapolate")( + target_height + ) + except ValueError: + sin_i = interp1d(heights, np.sin(rad).T, axis=1, fill_value="extrapolate")( + target_height + ) + cos_i = interp1d(heights, np.cos(rad).T, axis=1, fill_value="extrapolate")( + target_height + ) + return np.mod(np.rad2deg(np.arctan2(sin_i, cos_i)), 360.0) + + def _interpolate_with_min(heights, values, target_height, min_val=0.02): """Interpolate values to target height with minimum value clipping.""" try: @@ -589,40 +846,37 @@ def _interpolate_with_min(heights, values, target_height, min_val=0.02): ) -def configure_wake_model(system_dat, rotor_diameter, hub_height): +def configure_wake_model( + system_dat, + rotor_diameter, + hub_height, + resource_dat=None, + turbine_geometries=None, +): """Configure the wake model components based on system configuration. Args: system_dat: System data dictionary rotor_diameter: Rotor diameter for FUGA LUT generation hub_height: Hub height for FUGA LUT generation + resource_dat: Optional energy_resource dict; lets FUGA derive its LUT + roughness/inversion height from the site (z0 from TI, zi from + ABL_height) instead of using defaults. + turbine_geometries: Optional list of (rotor_diameter, hub_height) for + every turbine type; FUGA builds one LUT set per geometry so mixed + farms interpolate over d_h. Defaults to the single + (rotor_diameter, hub_height). Returns: dict with keys: wake_model_class, deficit_args, deflection_model, 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) @@ -630,35 +884,49 @@ 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, deficit_post_attrs = _configure_deficit_model( + wind_deficit_data, + analysis, + rotor_diameter, + hub_height, + resource_dat, + turbine_geometries, ) - - 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 + # WeightedSum/CumulativeWakeSum impose two constraints that PyWake otherwise + # enforces via bare AssertionErrors deep in a run. Fail fast with actionable + # messages; windIO cannot express these cross-field constraints. + from py_wake.deficit_models.deficit_model import ConvectionDeficitModel + from py_wake.rotor_avg_models.rotor_avg_model import NodeRotorAvgModel + from py_wake.superposition_models import CumulativeWakeSum, WeightedSum + + if isinstance(superposition_model, (WeightedSum, CumulativeWakeSum)): + # 1. Requires a node-based rotor-averaging model, or None (rotor centre, + # which PyWake accepts; an explicit RotorCenter is rejected). + if rotor_averaging is not None and not isinstance( + rotor_averaging, NodeRotorAvgModel + ): + raise ValueError( + "WeightedSum/CumulativeWakeSum superposition requires a node " + "rotor-averaging model (grid/eq_grid/gq_grid/cgi) or 'none'; " + "center, gaussian_overlap and area_overlap are not node models." + ) + # 2. Requires a convection-based deficit (carries a convective velocity). + if not issubclass(wake_model_class, ConvectionDeficitModel): + raise ValueError( + f"WeightedSum/CumulativeWakeSum superposition requires a " + f"ConvectionDeficitModel-based deficit (e.g. Zong2020); " + f"'{wind_deficit_data['name']}' " + f"({wake_model_class.__name__}) is not one. SuperGaussian, " + f"SuperGaussian2023 and GCL do not support it — use Linear." + ) + + # Blockage requires All2AllIterative solver solver_args = {} if blockage_model is not None: solver_class = All2AllIterative @@ -669,7 +937,8 @@ 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, + "deficit_post_attrs": deficit_post_attrs, + "wake_deficit_key": None, # Deprecated: kept for API compatibility "deflection_model": deflection_model, "turbulence_model": turbulence_model, "superposition_model": superposition_model, @@ -680,103 +949,554 @@ def configure_wake_model(system_dat, rotor_diameter, hub_height): } +# --- Fuga LUT generation ----------------------------------------------------- +# Fuga is a linearised-RANS wake model that reads a precomputed look-up table +# (LUT). Historically those came from a Windows GUI; pyfuga (conda-forge, pure +# Python) now generates them, so WIFA can build a LUT on the fly for any farm. +# +# Fuga has NO turbulence-intensity input: ambient turbulence enters implicitly +# through the roughness z0 (which sets the neutral shear/mixing) and the +# stability zeta0. We therefore derive z0 from the site's representative TI via +# the same inversion PyWake uses at runtime (z0 = zref * exp(-1/TI), neutral), +# unless the windIO config or the resource supplies z0 directly. + + +def _fuga_default_lut_dir(): + """Persistent, shared LUT cache dir (override with $WIFA_FUGA_LUT_DIR). + + LUTs are content-addressed by filename, so a single shared dir is safe and + lets the expensive preLUT stage be reused across runs and farms. + """ + return Path( + os.environ.get( + "WIFA_FUGA_LUT_DIR", Path.home() / ".cache" / "wifa" / "fuga_luts" + ) + ) + + +def _resource_field_array(resource_dat, key): + """Finite values of a windIO wind_resource field (dict-with-'data' or array). + + Returns a 1-D numpy array, or None if the field is absent/empty. + """ + if resource_dat is None: + return None + field = resource_dat.get("wind_resource", {}).get(key) + if field is None: + return None + data = np.asarray(field["data"] if isinstance(field, dict) else field, dtype=float) + data = data[np.isfinite(data)].ravel() + return data if data.size else None + + +def _mean_resource_field(resource_dat, key): + """Finite mean of a windIO wind_resource field, or None.""" + data = _resource_field_array(resource_dat, key) + return float(np.mean(data)) if data is not None else None + + +def _fuga_atmosphere(resource_dat, fuga_cfg, hub_height): + """Resolve (z0, zi, zeta0, ti) for a Fuga LUT from config + site resource. + + Precedence for z0/zi: explicit fuga config > site resource field > default. + z0 is derived from the site TI (Fuga has no TI knob) when not given. + """ + from py_wake.utils import fuga_utils + + zeta0 = float(fuga_cfg.get("zeta0", 0.0)) + + zi = fuga_cfg.get("zi") + if zi is None: + zi = _mean_resource_field(resource_dat, "ABL_height") + if zi is None: + zi = 500.0 + + ti = _mean_resource_field(resource_dat, "turbulence_intensity") + z0 = fuga_cfg.get("z0") + if isinstance(z0, (list, tuple)): + # An explicit z0 list is a sweep, handled by _fuga_z0_sweep; the single + # scalar here is only a fallback, so don't treat the list as scalar. + z0 = None + if z0 is None: + z0 = _mean_resource_field(resource_dat, "z0") + if z0 is None and ti is not None and ti > 0: + z0 = float(np.ravel(fuga_utils.z0(ti, hub_height, zeta0))[0]) + if z0 is None: + z0 = 0.03 # open-farmland fallback if neither z0 nor TI is available + return float(z0), float(zi), zeta0, ti + + +def _fuga_z0_sweep(resource_dat, fuga_cfg, hub_height, zeta0, z0_single): + """z0 values for a TI-faithful multi-LUT, spanning the site TI distribution. + + Fuga reads TI off the LUT roughness, so a single mean-TI LUT evaluates the + wake at loss(mean TI) and misses the low-TI tail that drives the deepest + wakes. A sweep of LUTs across z0 lets FugaDeficit interpolate z0 = z0(TI) + per flow case at run time, so the farm loss is integrated over the TI + distribution instead of taken at its mean (cf. the GCL free-stream-TI gap). + + Returns a sorted list of distinct z0. Falls back to ``[z0_single]`` when TI + data is unavailable, z0 is pinned in config, or n_z0 <= 1. Out-of-range TI + is handled by FugaDeficit's bounds='limit' (clamped to the nearest LUT). + """ + from py_wake.utils import fuga_utils + + if fuga_cfg.get("z0") is not None: + z0s = fuga_cfg["z0"] + z0s = z0s if isinstance(z0s, (list, tuple)) else [z0s] + return sorted({float(z) for z in z0s}) + + n = int(fuga_cfg.get("n_z0", 5)) + ti = _resource_field_array(resource_dat, "turbulence_intensity") + if n <= 1 or ti is None: + return [z0_single] + ti = ti[ti > 0] + if ti.size == 0: + return [z0_single] + lo, hi = np.quantile( + ti, [fuga_cfg.get("ti_qlo", 0.05), fuga_cfg.get("ti_qhi", 0.95)] + ) + # Clamp to a TI band that keeps the neutral-inversion z0 physical: the + # mapping z0 = zhub*exp(-1/TI) sends high TI to absurd roughness (TI 0.30 -> + # z0 ~2.8 m, well outside Fuga's linearisation). The low-TI tail drives the + # deepest, most TI-sensitive wakes, so cover it; high-TI cases saturate to + # shallow wakes and clamp to the roughest LUT via bounds='limit'. Band + # [0.03, 0.18] keeps z0 in ~[1e-5, 0.3] m. + ti_lo = float(fuga_cfg.get("ti_min", 0.03)) + ti_hi = float(fuga_cfg.get("ti_max", 0.18)) + lo, hi = float(np.clip(lo, ti_lo, ti_hi)), float(np.clip(hi, ti_lo, ti_hi)) + + def _z0(ti_val): + return round(float(np.ravel(fuga_utils.z0(ti_val, hub_height, zeta0))[0]), 8) + + if hi <= lo: + # Whole TI distribution sits at/over a clamp bound -> a single LUT at + # the clamped TI (still physical), not the unclamped mean-TI z0. + return [_z0(hi)] + return sorted({_z0(t) for t in np.linspace(lo, hi, n)}) + + +def _ensure_fuga_luts( + *, + folder, + zeta0, + nkz0, + nbeta, + geometries, + z0_list, + zi, + lut_vars, + nx, + ny, + zlow=None, + zhigh=None, + dx=None, + dy=None, + n_cpu=None, +): + """Generate/reuse a LUT for every (geometry, z0) pair; return the path list. + + All LUTs share the costly preLUT (which depends only on zeta0/nkz0/nbeta), + so extra z0 values and turbine geometries add only the cheap per-LUT stage. + FugaDeficit/FugaMultiLUTDeficit interpolate the resulting set over d_h + (turbine geometry) and z0 (per-flow-case TI). + + zlow/zhigh/dx/dy default to each geometry's own hub height and D/4, D/16. + Mixed-geometry layouts must pass a shared zlow/zhigh (spanning every hub + height) and a shared dx/dy so FugaMultiLUTDeficit can merge the LUTs onto + one z/x/y grid: merging single-height LUTs at different hub heights turns + the whole table NaN (xarray cannot interpolate a size-1 z axis), which + surfaced as zero power at every cross-type waked turbine. + """ + paths = [] + for diameter, zhub in geometries: + for z0 in z0_list: + paths.append( + _ensure_fuga_lut( + folder=folder, + zeta0=zeta0, + nkz0=nkz0, + nbeta=nbeta, + diameter=diameter, + zhub=zhub, + z0=z0, + zi=zi, + lut_vars=lut_vars, + nx=nx, + ny=ny, + zlow=zlow, + zhigh=zhigh, + dx=dx, + dy=dy, + n_cpu=n_cpu, + ) + ) + return paths + + +def _ensure_fuga_lut( + *, + folder, + zeta0, + nkz0, + nbeta, + diameter, + zhub, + z0, + zi, + lut_vars, + nx, + ny, + zlow=None, + zhigh=None, + dx=None, + dy=None, + n_cpu=None, +): + """Generate (or reuse a cached) Fuga LUT; return its path. + + The LUT filename encodes every physical/grid parameter, so an existing file + with the right name is a valid cache hit. pyfuga reuses the costly preLUT + stage (which depends only on zeta0/nkz0/nbeta) across geometries. + """ + from pyfuga import get_luts + from pyfuga.paths import get_luts_path + + folder = Path(folder) + folder.mkdir(parents=True, exist_ok=True) + # pyfuga's own defaults; pass explicitly so the cache-probe path built by + # get_luts_path matches the filename get_luts actually writes. + if dx is None: + dx = diameter / 4 + if dy is None: + dy = diameter / 16 + # zlow == zhigh == zhub -> single hub-height level (the cheap path). + if zlow is None: + zlow = zhub + if zhigh is None: + zhigh = zhub + lut_vars = list(lut_vars) + lut_path = get_luts_path( + folder, + zeta0, + nkz0, + nbeta, + diameter, + zhub, + z0, + zi, + zlow, + zhigh, + lut_vars, + nx, + ny, + dx, + dy, + ) + if not lut_path.exists(): + get_luts( + folder=folder, + zeta0=zeta0, + nkz0=nkz0, + nbeta=nbeta, + diameter=diameter, + zhub=zhub, + z0=z0, + zi=zi, + zlow=zlow, + zhigh=zhigh, + lut_vars=lut_vars, + nx=nx, + ny=ny, + dx=dx, + dy=dy, + n_cpu=n_cpu, + ) + return str(lut_path) + + def _configure_deficit_model( - wind_deficit_data, analysis, rotor_diameter, hub_height, deficit_args + wind_deficit_data, + analysis, + rotor_diameter, + hub_height, + resource_dat=None, + turbine_geometries=None, ): """Configure the wind deficit model. Returns: - tuple: (wake_model_class, deficit_args, wake_deficit_key) + tuple: (wake_model_class, deficit_args, deficit_post_attrs) where + deficit_post_attrs is a dict of attributes to set on the built + deficit after construction (e.g. TurbOPark's WS_key). """ 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) + + wind_deficit_cfg = analysis.get("wind_deficit_model", {}) + # Honor the windIO use_effective_ws flag (local vs free-stream inflow at the + # waking turbine); deficits that don't accept it pop it below (NOJDeficit). + deficit_args = {"use_effective_ws": wind_deficit_cfg.get("use_effective_ws", True)} + 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"} + # Deficits that expose a use_effective_ti param (TI-dependent expansion/width). + # NOJLocalDeficit (Jensen) accepts it too: with a=[k_a, k_b] it references + # effective TI, so honoring free_stream_ti lets a no-turbulence config use + # ambient TI. GCLDeficit also accepts it (GCLLocal sets use_effective_ti=True). + # Bastankhah2014, free-stream NOJDeficit and FUGA do not. + TI_CAPABLE = { + "jensen", + "nojlocaldeficit", + "niayifar2016", + "carbajofuertes2018", + "zong2020", + "turbopark", + "gcl", + "supergaussian", + "supergaussian2023", + "eddyviscosity", + } - if model_name == "Jensen": + # windIO convention: k = k_a + k_b * TI (k_a constant, k_b multiplies TI). + # PyWake's a-parametrized deficits compute k = a[0] * TI + a[1], so the + # windIO pair maps to a = [k_b, k_a]. Scalar-k deficits take k_a and cannot + # represent a nonzero k_b. + 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: + if "k_a" in wake_expansion or "k_b" in wake_expansion: + deficit_args["a"] = [ + wake_expansion.get("k_b", 0) or 0, + wake_expansion.get("k_a", 0) or 0, + ] + + elif normalized in ("jensen1983", "nojdeficit"): + wake_model_class = NOJDeficit + deficit_args.pop("use_effective_ws", None) + # NOJDeficit takes a scalar k. windIO's wake_expansion_coefficient has + # no scalar `k` field, so accept k_a (the constant) as well as `k`. + 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 "k_a" in wake_expansion: + deficit_args["k"] = wake_expansion["k_a"] + if wake_expansion.get("k_b"): + warnings.warn( + f"{model_name} takes a constant wake expansion k (= k_a); " + f"the TI coefficient k_b={wake_expansion['k_b']} is ignored." + ) - elif model_name == "SuperGaussian": + elif normalized in GAUSSIAN_MODELS: + wake_model_class = GAUSSIAN_MODELS[normalized] + if normalized in A_PARAM_MODELS: + # Niayifar, Zong, Carbajofuertes: k = k_a + k_b*TI -> a=[k_b, k_a] + if "k" in wake_expansion: + warnings.warn( + f"{model_name} uses k_a/k_b (k = k_a + k_b*TI) for wake " + f"expansion, not scalar k. Scalar 'k' is ignored." + ) + if "k_a" in wake_expansion or "k_b" in wake_expansion: + if "k_b" not in wake_expansion: + warnings.warn( + f"k_b not specified for {model_name}, defaulting to 0 " + f"(TI-independent wake expansion)" + ) + deficit_args["a"] = [ + wake_expansion.get("k_b", 0) or 0, + wake_expansion.get("k_a", 0) or 0, + ] + else: + # Bastankhah2014 uses k (scalar) = the windIO constant k_a + if "k_a" in wake_expansion: + deficit_args["k"] = wake_expansion["k_a"] + elif "k" in wake_expansion: + deficit_args["k"] = wake_expansion["k"] + if wake_expansion.get("k_b"): + warnings.warn( + f"{model_name} takes a constant wake expansion k (= k_a); " + f"the TI coefficient k_b={wake_expansion['k_b']} is ignored." + ) + # ceps maps to the deficit's near-wake epsilon coefficient. Bastankhah, + # Niayifar and Carbajofuertes name it `ceps`; Zong names it `eps_coeff`. + if "ceps" in wind_deficit_cfg: + if normalized == "zong2020": + deficit_args["eps_coeff"] = wind_deficit_cfg["ceps"] + else: + deficit_args["ceps"] = wind_deficit_cfg["ceps"] + + 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 + # Canonical Nygaard (2022) recipe (py_wake.literature.turbopark): a Mirror + # ground model and ctlim=0.96 as constructor args; the WS_key='WS_jlk' + # attribute (scale the deficit by the downstream turbine's ambient WS) is + # applied post-construction via deficit_post below. + from py_wake.ground_models.ground_models import Mirror + from py_wake.superposition_models import SquaredSum + + deficit_args["groundModel"] = Mirror(superpositionModel=SquaredSum()) + deficit_args["ctlim"] = 0.96 + + 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 == "eddyviscosity": + # Only on pyWake's unmerged EV branch (cj_add_eddy_viscosity_model); + # import lazily so the released pyWake keeps working without it. + from py_wake.deficit_models.eddy_viscosity import EddyViscosityDeficitModel + + wake_model_class = EddyViscosityDeficitModel + # EV (Ainslie 1988) references deficits to the free-stream wind speed + # and combines them with MaxSum (the WindFarmer convention), so unlike + # the analytic deficits use_effective_ws defaults to False. + deficit_args["use_effective_ws"] = wind_deficit_cfg.get( + "use_effective_ws", False + ) - elif model_name.upper() == "FUGA": - wake_model_class = FugaDeficit - from pyfuga import get_luts - - lut = get_luts( - folder="luts", - zeta0=0, - nkz0=8, - nbeta=32, - diameter=rotor_diameter, - zhub=hub_height, - z0=0.00001, - zi=500, - zlow=70, - zhigh=70, - lut_vars=["UL"], - nx=2048, - ny=512, - n_cpu=1, + elif normalized == "bastankhah2016": + raise NotImplementedError( + "Bastankhah2016 is not available in PyWake. Use flow_model 'foxes', " + "or choose Bastankhah2014/Zong2020 for PyWake." ) - deficit_args["LUT_path"] = ( - f"luts/LUTs_Zeta0=0.00e+00_8_32_D{rotor_diameter:.1f}_zhub{hub_height:.1f}" - f"_zi500_z0=0.00001000_z69.2-72.8_UL_nx2048_ny512_dx44.575_dy11.14375.nc" + + elif normalized == "fuga": + wake_model_class = FugaDeficit + # FugaDeficit reads a LUT instead of an analytic expansion; it takes no + # use_effective_ws (it always uses the free-stream-referenced deficit). + deficit_args.pop("use_effective_ws", None) + fuga_cfg = wind_deficit_cfg.get("fuga", {}) or {} + z0_single, zi, zeta0, _ti = _fuga_atmosphere(resource_dat, fuga_cfg, hub_height) + # A z0 sweep across the site TI distribution + a LUT per turbine + # geometry; FugaDeficit interpolates z0 (per-flow-case TI) and d_h at + # run time. Degenerates to a single LUT for one geometry + n_z0<=1. + z0_list = _fuga_z0_sweep(resource_dat, fuga_cfg, hub_height, zeta0, z0_single) + geometries = turbine_geometries or [(rotor_diameter, hub_height)] + # Dedupe: two turbine types with the same geometry share one LUT, and + # duplicate d_h coordinates would break FugaMultiLUTDeficit's merge. + geometries = list(dict.fromkeys((float(d), float(h)) for d, h in geometries)) + hub_heights = sorted({h for _, h in geometries}) + diameters = sorted({d for d, _ in geometries}) + # Mixed hub heights: every LUT must span all hub heights (a source + # turbine's wake is evaluated at each target's hub height), and mixed + # diameters need one shared x/y grid; otherwise FugaMultiLUTDeficit's + # merge yields NaN deficits -> zero power at cross-type waked turbines. + mixed_grid = {} + if len(hub_heights) > 1: + mixed_grid["zlow"] = hub_heights[0] + mixed_grid["zhigh"] = hub_heights[-1] + # Interpolate the merged LUTs at exactly the hub heights; pyfuga's + # log-spaced z levels differ per z0, so without this the z-union + # across a z0 sweep would reintroduce NaN edge cells. + deficit_args["z_lst"] = hub_heights + if len(diameters) > 1: + # Finest natural resolution; identical x/y coords across LUTs. + mixed_grid["dx"] = diameters[0] / 4 + mixed_grid["dy"] = diameters[0] / 16 + lut_paths = _ensure_fuga_luts( + folder=fuga_cfg.get("cache_dir", _fuga_default_lut_dir()), + zeta0=zeta0, + nkz0=int(fuga_cfg.get("nkz0", 8)), + nbeta=int(fuga_cfg.get("nbeta", 32)), + geometries=geometries, + z0_list=z0_list, + zi=zi, + lut_vars=fuga_cfg.get("lut_vars", ["UL"]), + nx=int(fuga_cfg.get("nx", 2048)), + ny=int(fuga_cfg.get("ny", 512)), + n_cpu=fuga_cfg.get("n_cpu"), + **mixed_grid, ) + # Single LUT -> plain path; multiple -> list (FugaDeficit globs/lists). + deficit_args["LUT_path"] = lut_paths[0] if len(lut_paths) == 1 else lut_paths 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] + # TI reference: windIO carries this as the nested wake_expansion_coefficient + # .free_stream_ti flag (foxes-compatible). PyWake's deficits expose the + # inverse use_effective_ti param (use_effective_ti = not free_stream_ti), + # but only the TI-dependent deficits accept it. + if normalized in TI_CAPABLE and "free_stream_ti" in wake_expansion: + deficit_args["use_effective_ti"] = not wake_expansion["free_stream_ti"] + + # Axial induction: windIO's axial_induction_model maps to PyWake's ct2a + # (1D -> ct2a_mom1d, Madsen -> ct2a_madsen). Honor it on every deficit that + # accepts a ct2a parameter; without this the deficit silently keeps its + # ct2a_madsen default, so a "1D" request was previously ignored on the + # pyWake path. + import inspect + + from py_wake.deficit_models.utils import ct2a_madsen, ct2a_mom1d + + axial = analysis.get("axial_induction_model") + if axial is not None: + ct2a_fn = {"1d": ct2a_mom1d, "madsen": ct2a_madsen}.get(_normalize_name(axial)) + if ( + ct2a_fn is not None + and "ct2a" in inspect.signature(wake_model_class.__init__).parameters + ): + deficit_args["ct2a"] = ct2a_fn + + # Attributes set on the deficit *after* construction (not constructor + # kwargs), applied by run_simulation. + deficit_post = {} + if normalized == "turbopark": + deficit_post["WS_key"] = "WS_jlk" - return wake_model_class, deficit_args, wake_deficit_key + return wake_model_class, deficit_args, deficit_post 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["name"].lower() - if name == "none": + name = deflection_data.get("name") + if name is None: return None - elif name == "jimenez": + + normalized = _normalize_name(name) + if normalized == "none": + return None + 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): @@ -786,67 +1506,209 @@ 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": + c = turbulence_data.get("c") + if c is not None: + # A paper's calibration (e.g. Niayifar 2016, Zong 2020): the + # literature CrespoHernandez uses 1D induction and SqrMaxSum + # added-turbulence summation alongside the calibrated coefficients. + from py_wake.deficit_models.utils import ct2a_mom1d + from py_wake.superposition_models import SqrMaxSum + + return CrespoHernandez( + c=list(c), + ct2a=ct2a_mom1d, + addedTurbulenceSuperpositionModel=SqrMaxSum(), + ) return CrespoHernandez() - else: - raise NotImplementedError( - f"Turbulence model '{turbulence_data['name']}' is not supported" + if normalized == "gcl": + return GCLTurbulence() + if normalized in ("quartonandainslie", "modifiedquartonandainslie"): + # Only on pyWake's unmerged EV branch; the Hassan (1992) modified + # variant is the one the EV bundle uses. It carries its own added-TI + # combination (SqrMaxSum default); windIO's ti_superposition is not + # consulted on the pyWake path. + from py_wake.turbulence_models.quarton_and_ainslie import ( + ModifiedQuartonAndAinslieTurbulenceModel, ) + return ModifiedQuartonAndAinslieTurbulenceModel() + 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: + 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.") + if normalized == "vector": raise NotImplementedError( - f"Superposition model '{superposition_data['ws_superposition']}' is not supported" + "Vector superposition is foxes-only; 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 ( + AreaOverlapAvgModel, + CGIRotorAvg, + EqGridRotorAvg, + GaussianOverlapAvgModel, + GQGridRotorAvg, + GridRotorAvg, + PolarGridRotorAvg, + RotorCenter, + ) + + name = rotor_avg_data["name"] + normalized = _normalize_name(name) - name = rotor_avg_data["name"].lower() - if name == "center": - print("Using Center Average") + if normalized == "none": + # No rotor-averaging model. PyWake's Weighted superposition accepts this + # (rotor centre) but rejects an explicit RotorCenter; the Zong (2020) + # literature model uses None. + return None + if normalized == "center": return RotorCenter() - elif name == "avg_deficit": + # "grid" is the canonical windIO name; "avgdeficit" is a deprecated alias + if normalized in ("grid", "avgdeficit"): return GridRotorAvg() - else: - raise NotImplementedError( - f"Rotor averaging model '{rotor_avg_data['name']}' is not supported" + if normalized == "gaussianoverlap": + return GaussianOverlapAvgModel() + if normalized == "areaoverlap": + return AreaOverlapAvgModel() + 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)) + if normalized == "simplifiedgaussian": + # Only on pyWake's unmerged EV branch: the LUT-based line average + # across the rotor that the EV bundle pairs with its deficit. + from py_wake.rotor_avg_models.simplified_gaussian_rotor_average_model import ( + SimplifiedGaussianRotorAverageModel, ) + return SimplifiedGaussianRotorAverageModel() + 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 + from py_wake.superposition_models import LinearSum 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": - return FugaDeficit(deficit_args["LUT_path"]) - else: - raise ValueError(f"Unknown blockage model: {name}") + + normalized = _normalize_name(name) + if normalized == "none": + return None + + # The analytic blockage models are calibrated without ground effects and + # allow flow through the ground; ground_mirror enforces the slip boundary + # condition with an image rotor (unlike wake models, which are calibrated + # including ground effects and must not be mirrored). + ground_model = None + if blockage_data.get("ground_mirror", False): + if normalized == "fuga": + warnings.warn( + "blockage_model.ground_mirror is ignored for FUGA: the Fuga " + "LUTs come from a linearized RANS solver that already includes " + "the ground." + ) + else: + from py_wake.ground_models.ground_models import Mirror + + ground_model = Mirror() + + # Models that take no constructor arguments + SIMPLE_BLOCKAGE_MODELS = { + "selfsimilaritydeficit": SelfSimilarityDeficit, + "rankinehalfbody": RankineHalfBody, + "rathmann": Rathmann, + "vortexcylinder": VortexCylinder, + "vortexdipole": VortexDipole, + "hybridinduction": HybridInduction, + } + + # Sum blockage deficits (across turbines, and real+image when mirrored) + # linearly: potential-flow induction superposes linearly, and pyWake + # otherwise falls back to the *wake* superposition model — SquaredSum + # (e.g. the TurbOPark recipe) asserts on the speed-up (negative deficit) + # regions every blockage model produces. + blockage_superposition = LinearSum() + + if normalized == "selfsimilaritydeficit2020": + return SelfSimilarityDeficit2020( + ss_alpha=blockage_data.get("ss_alpha", 0.8888888888888888), + groundModel=ground_model, + superpositionModel=blockage_superposition, + ) + if normalized in SIMPLE_BLOCKAGE_MODELS: + return SIMPLE_BLOCKAGE_MODELS[normalized]( + groundModel=ground_model, + superpositionModel=blockage_superposition, + ) + if normalized == "fuga": + return FugaDeficit(deficit_args["LUT_path"], z_lst=deficit_args.get("z_lst")) + raise NotImplementedError(f"Blockage model '{name}' is not supported") def run_simulation(site, turbine, wake_config, site_data, x, y, turbine_types): @@ -864,16 +1726,17 @@ def run_simulation(site, turbine, wake_config, site_data, x, y, turbine_types): Returns: dict with keys: sim_res, aep, aep_per_turbine """ - # Build deficit model - print("Running ", wake_config["wake_model_class"], wake_config["deficit_args"]) + # Build deficit model. groundModel comes from deficit_args when a model needs + # a specific one (e.g. TurbOPark's Mirror); otherwise the deficit's own + # default (None) applies. + deficit_args = dict(wake_config["deficit_args"]) + deficit_args.setdefault("groundModel", None) deficit_model = wake_config["wake_model_class"]( rotorAvgModel=wake_config["rotor_averaging"], - groundModel=None, - **wake_config["deficit_args"], + **deficit_args, ) - - if wake_config["wake_deficit_key"]: - deficit_model.WS_key = wake_config["wake_deficit_key"] + for attr, value in wake_config.get("deficit_post_attrs", {}).items(): + setattr(deficit_model, attr, value) # Build wind farm model wind_farm_model = wake_config["solver_class"]( @@ -905,20 +1768,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} @@ -938,9 +1791,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", {}) @@ -948,20 +1800,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", @@ -969,7 +1818,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) @@ -984,70 +1833,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" - ) - yaml_content = yaml_content.replace( - "INCLUDE_POWER_TABLE_PLACEHOLDER", "!include PowerTable.nc" - ) - yaml_content = yaml_content.replace( - "INCLUDE_FLOW_FIELD_PLACEHOLDER", "!include FarmFlow.nc" + # 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" ) - - 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"): @@ -1091,16 +1922,25 @@ def run_pywake(yaml_input, output_dir="output"): site = site_data["site"] # Step 4: Configure wake model - # Use first turbine's dimensions for FUGA LUT if needed - first_hh = list(hub_heights.values())[0] - # Get rotor diameter from farm data + # Collect every turbine geometry so FUGA can build a LUT set per type + # (mixed farms interpolate over d_h); the first one drives single-type paths. if "turbines" in farm_dat: - rd = farm_dat["turbines"]["rotor_diameter"] + geometries = [ + ( + farm_dat["turbines"]["rotor_diameter"], + farm_dat["turbines"]["hub_height"], + ) + ] else: - first_key = list(farm_dat["turbine_types"].keys())[0] - rd = farm_dat["turbine_types"][first_key]["rotor_diameter"] + geometries = [ + (t["rotor_diameter"], t["hub_height"]) + for t in farm_dat["turbine_types"].values() + ] + rd, first_hh = geometries[0] - wake_config = configure_wake_model(system_dat, rd, first_hh) + wake_config = configure_wake_model( + system_dat, rd, first_hh, resource_dat, geometries + ) # Step 5: Run simulation sim_results = run_simulation( diff --git a/wifa/wayve_api.py b/wifa/wayve_api.py index f731d19..fa49c6e 100644 --- a/wifa/wayve_api.py +++ b/wifa/wayve_api.py @@ -19,7 +19,6 @@ def run_wayve(yamlFile, output_dir="output", debug_mode=False): from wayve.apm import APM from wayve.grid.grid import Stat2Dgrid from wayve.momentum_flux_parametrizations import FrictionCoefficients - from wayve.pressure.gravity_waves.gravity_waves import NonUniform, Uniform from wayve.solvers import FixedPointIteration ##################### @@ -80,24 +79,41 @@ def run_wayve(yamlFile, output_dir="output", debug_mode=False): if h1 < h1_min: raise UserWarning("Lower layer height too low, please specify a higher value") + # Optional ABL-setup knobs (analysis.abl_setup), forwarded to flow_io_abl: + # geostrophic-wind mode and capping-inversion fit settings. + abl_setup = analysis_dat.get("abl_setup", {}) + abl_kwargs = {} + if "Gmode" in abl_setup: + abl_kwargs["gmode"] = abl_setup["Gmode"] + if "dh_max" in abl_setup: + abl_kwargs["dh_max"] = abl_setup["dh_max"] + if "serz" in abl_setup: + abl_kwargs["serz"] = bool(abl_setup["serz"]) + ################## # Other APM components ################## # Momentum flux parametrization mfp = FrictionCoefficients() - # Pressure feedback parametrization - pressure = Uniform(dynamic=True, rotating=False) + # Pressure feedback parametrization: any value but 1 selects the + # profile-resolving NonUniform gravity-wave closure (<= 0 means "one + # sublayer per profile level"). The object is built per state + # (_pressure_for_state), because whether the profile actually supports + # NonUniform depends on the per-state inversion fit. + n_fa_layers = 1 if "layers_description" in analysis_dat: - if "number_of_fa_layers" in analysis_dat["layers_description"]: - n_layers = analysis_dat["layers_description"]["number_of_fa_layers"] - if n_layers > 1: - pressure = NonUniform(n_layers=n_layers, order=1) + n_fa_layers = analysis_dat["layers_description"].get("number_of_fa_layers", 1) ###################### # Read output settings ###################### - # Select timestamps - times = resource_dat["wind_resource"]["time"] + # Select timestamps. `time_indices` are positions into the *full* wind + # resource arrays; they must be carried alongside the timestamp labels, + # because flow_io_abl() indexes the full arrays. Enumerating the subsetted + # timestamps instead would simulate rows 0..n-1 while labelling them with + # the requested timestamps. + all_times = resource_dat["wind_resource"]["time"] + time_indices = list(range(len(all_times))) run_config = system_dat["attributes"]["model_outputs_specification"][ "run_configuration" ] @@ -105,8 +121,8 @@ def run_wayve(yamlFile, output_dir="output", debug_mode=False): "all_occurences", True ): if "subset" in run_config["times_run"]: - subset = run_config["times_run"]["subset"] - times = [times[i] for i in subset] + time_indices = list(run_config["times_run"]["subset"]) + times = [all_times[i] for i in time_indices] # Get turbine variables to output turbine_nc_filename = "turbine_data.nc" turbine_output_variables = ["power", "rotor_effective_velocity"] @@ -116,7 +132,11 @@ def run_wayve(yamlFile, output_dir="output", debug_mode=False): ] if "turbine_nc_filename" in turb_out_dat: turbine_nc_filename = turb_out_dat["turbine_nc_filename"] - if "turbine_output_variables" in turb_out_dat: + # The schema (and every other WIFA engine) calls this `output_variables`; + # `turbine_output_variables` is kept as a fallback for legacy yamls. + if "output_variables" in turb_out_dat: + turbine_output_variables = turb_out_dat["output_variables"] + elif "turbine_output_variables" in turb_out_dat: turbine_output_variables = turb_out_dat["turbine_output_variables"] # Check flow field output specification flow_nc_filename = "flow_field.nc" @@ -168,19 +188,39 @@ def run_wayve(yamlFile, output_dir="output", debug_mode=False): ##################### # Perform model runs ##################### + # Validate the capping-inversion spelling once, before the loop: an + # ambiguous or incomplete block is a defect of the file, and the loop's + # ``except Exception`` would otherwise report it as every state crashing. + capping_inversion_spelling(resource_dat["wind_resource"]) # Initialize crash counter crashes = 0 + # NonUniform free-atmosphere tally (per-state fallback is warned once) + nonuniform_states = 0 # List of datasets ds_list = [] ds_ff_list = [] # Loop over timeseries - for time_index, time in enumerate(times): + for run_index, time_index in enumerate(time_indices): + time = times[run_index] if debug_mode: # Print timestep - print(f"time {time_index+1}/{len(times)}") + print(f"time {run_index + 1}/{len(times)}") try: - # Set up ABL - abl = flow_io_abl(resource_dat["wind_resource"], time_index, hh, h1) + # Set up ABL (solver frame: hub-height wind along +x) + abl, rotation = flow_io_abl( + resource_dat["wind_resource"], time_index, hh, h1, **abl_kwargs + ) + # Rebuild the wind farm in the solver frame: wayve's gravity-wave + # solve assumes westerly flow, so the layout turns with the wind. + wind_farm, forcing, wf_offset_x, wf_offset_y = wf_setup( + farm_dat, analysis_dat, L_filter, debug_mode, rotation=rotation + ) + coupling = wind_farm.coupling + wake_model = coupling.wake_model + # Pressure feedback for this state + pressure = _pressure_for_state(abl, n_fa_layers) + if type(pressure).__name__ == "NonUniform": + nonuniform_states += 1 # Set up APM from components model = APM(grid, forcing, abl, mfp, pressure) # Use a fixed-point iteration solver with a relaxation factor of 0.7 @@ -203,14 +243,20 @@ def run_wayve(yamlFile, output_dir="output", debug_mode=False): ["turbine"], wind_farm.coupling.St, ) + # State-level ABL diagnostics. The capping inversion and the + # geostrophic wind are often *fitted* here rather than given (see + # flow_io_abl), so without these the atmosphere a run actually used + # is unrecoverable from its outputs. Scalars, no turbine dim. + diagnostics = _abl_diagnostics(abl, pressure) + assert not set(diagnostics) & set(turb_out_dict) + turb_out_dict |= diagnostics # NC setup ds = xr.Dataset( turb_out_dict, coords={"states": time, "turbine": range(Nt)}, ) - # Add to output list - ds_list.append(ds) # Flow field outputs # + ds_ff = None if report_flow and not debug_mode: # Callables for flow evaluation u_bg_evaluator = coupling.set_up_u_bg_evaluator( @@ -220,22 +266,35 @@ def run_wayve(yamlFile, output_dir="output", debug_mode=False): # Output arrays wind_speed = np.zeros([len(x_ff), len(y_ff), len(z_ff)]) wind_dir = np.zeros([len(x_ff), len(y_ff), len(z_ff)]) + # Query points: the requested earth-frame grid, rotated into + # the solver frame (where the farm sits and the wind blows + # along +x). Evaluating the rotated points directly keeps the + # output on the requested grid with no regridding. + c, s = np.cos(rotation), np.sin(rotation) + x_e, y_e = np.meshgrid( + x_ff - wf_offset_x, y_ff - wf_offset_y, indexing="ij" + ) + x_q = c * x_e + s * y_e + y_q = c * y_e - s * x_e # Loop over z-planes for k, z_k in enumerate(z_ff): - # Get velocities - u_bg, v_bg, u_wm, v_wm = wake_model.xy_plane( + # Get velocities (solver frame) + u_bg, v_bg, u_wm, v_wm = _xy_plane_points( + wake_model, wind_farm, abl, u_bg_evaluator, apm_evaluator, - x_ff - wf_offset_x, - y_ff - wf_offset_y, + x_q, + y_q, z_k, ) - # Convert to speed and direction + # Rotate velocity vectors back to the earth frame + u_wm, v_wm = c * u_wm - s * v_wm, s * u_wm + c * v_wm + # Convert to speed and direction (wrapped to [0, 360)) wind_speed[:, :, k] = np.sqrt(np.square(u_wm) + np.square(v_wm)) - wind_dir[:, :, k] = np.rad2deg( - np.pi / 2 - (np.arctan2(v_wm, u_wm) + np.pi) + wind_dir[:, :, k] = ( + np.rad2deg(np.pi / 2 - (np.arctan2(v_wm, u_wm) + np.pi)) % 360.0 ) # Flow output dictionary flow_out_dict = {} @@ -248,7 +307,10 @@ def run_wayve(yamlFile, output_dir="output", debug_mode=False): flow_out_dict, coords={"states": time, "x": x_ff, "y": y_ff, "z": z_ff}, ) - # Add to output list + # Append outputs together, so a flow-field failure cannot leave + # turbine_data.nc and flow_field.nc with different states axes. + ds_list.append(ds) + if ds_ff is not None: ds_ff_list.append(ds_ff) except Exception as exc: @@ -258,6 +320,12 @@ def run_wayve(yamlFile, output_dir="output", debug_mode=False): continue if debug_mode: print(f"crashes: {crashes}/{len(times)}") + if n_fa_layers != 1: + # Make the per-state Uniform/NonUniform closure mixture visible + print( + f"NonUniform free atmosphere: {nonuniform_states}/{len(times)} " + "states (the rest fell back to Uniform)" + ) # Combine into total dataset output_dir = Path(output_dir) @@ -271,6 +339,253 @@ def run_wayve(yamlFile, output_dir="output", debug_mode=False): ds_ff_full.to_netcdf(output_fn) +def _abl_diagnostics(abl, pressure): + """Per-state record of the atmosphere the solver actually ran on. + + ``flow_io_abl`` fits the capping inversion when the wind resource does not + state one, derives the geostrophic wind through ``Gmode``, and closes the + stress profile from surface scalars; ``_pressure_for_state`` then picks the + free-atmosphere closure per state. All of that is invisible in the turbine + outputs, so record it alongside them: it is what makes a run reproducible + and its inversion checkable against an independent boundary-layer height. + + Names are prefixed ``abl_`` so they cannot collide with a turbine output: + ``air_density`` and ``turbulence_intensity`` are per-turbine fields in other + engines, and these share one dataset with them. + + Returned as ``{name: ([], value)}`` — scalars on the state axis, no turbine + dimension. + """ + gravity = 9.80665 # [m s-2], as in flow_io_abl + # gprime = g * dtheta / theta0, with theta0 the mixed-layer potential + # temperature at the inversion; invert it to report dtheta itself. + dtheta = abl.gprime * np.interp(abl.H, abl.zs, abl.ths) / gravity + values = { + "abl_height": (abl.H, "m"), + "abl_capping_inversion_thickness": (abl.inv_top - abl.inv_bottom, "m"), + "abl_capping_inversion_strength": (dtheta, "K"), + "abl_free_atmosphere_N": (abl.N, "s-1"), + "abl_geostrophic_wind_speed": (abl.S3, "m s-1"), + # Angle of the geostrophic wind in the solver frame, where the + # hub-height wind lies along +x: the cross-isobar angle, not a + # compass direction. + "abl_geostrophic_veer": (abl.WD3, "degree"), + # In the turbulence-profile branch wayve's ``utau`` carries the surface + # stress itself (m2 s-2), not its square root; everywhere else it is a + # velocity. Reported as wayve holds it. + "abl_friction_velocity": (abl.utau, "m s-1"), + "abl_air_density": (abl.rho, "kg m-3"), + "abl_turbulence_intensity": (abl.TI, "1"), + # 1 when this state resolved the free atmosphere on the profile, + # 0 when it fell back to the bulk closure. + "abl_nonuniform_free_atmosphere": ( + float(type(pressure).__name__ == "NonUniform"), + "1", + ), + } + # Reporting must not be able to cost a state its results: anything + # unexpected here (a None, a non-scalar) becomes NaN rather than an + # exception that the caller's `except Exception` would turn into a drop. + out = {} + for name, (value, units) in values.items(): + try: + out[name] = ([], float(value)) + except (TypeError, ValueError): + out[name] = ([], float("nan")) + out[name] = (*out[name], {"units": units}) + return out + + +def _pressure_for_state(abl, n_fa_layers): + """Build the gravity-wave pressure closure for one state. + + Anything other than exactly 1 requests wayve's NonUniform closure, which + resolves the free-atmosphere N(z) and wind shear by slicing the profile + between the inversion top and ``abl.h_strat`` into layers. A value <= 0 is + wayve's own "use the ABL's vertical grid" setting: the sublayers follow the + profile levels themselves, so a reanalysis column resolves the free + atmosphere at the resolution it actually has instead of an arbitrary count. + + Either way that needs actual profile points in that range: truncated or + synthetic profiles (e.g. the scalar branch's Nieuwstadt profile, which + stops at the inversion) fall back to the bulk Uniform closure with a + warning. The warning text is deliberately state-independent so Python's + warning dedup collapses it on long time series; run_wayve prints a per-run + NonUniform/Uniform tally. + """ + from wayve.pressure.gravity_waves.gravity_waves import NonUniform, Uniform + + if n_fa_layers != 1: + h_min = max(abl.H, abl.inv_top if abl.inv_top is not None else 0.0) + n_pts = int(np.sum((abl.zs > h_min) & (abl.zs < abl.h_strat))) + if n_pts >= 2: + return NonUniform(n_layers=n_fa_layers, order=1) + warnings.warn( + "number_of_fa_layers requested but too few profile points lie " + "between the inversion top and the free-atmosphere top; falling " + "back to the Uniform free atmosphere for such states" + ) + return Uniform(dynamic=True, rotating=False) + + +def _xy_plane_points( + wake_model, wind_farm, abl, u_bg_evaluator, apm_evaluator, Xs, Ys, z +): + """Evaluate the coupled flow at arbitrary (x, y) coordinate arrays. + + Mirror of wayve 2.0.0's ``xy_plane`` methods (Lanzilao and foxes wake + models) with the axis-aligned ``meshgrid(xs, ys)`` replaced by + caller-provided 2-D coordinate arrays, so a rotated (solver-frame) grid + can be evaluated directly — every operation downstream of the meshgrid is + pointwise in the coordinates. Returns ``(u_bg, v_bg, u_wm, v_wm)`` with + the shape of ``Xs``. Equivalence with ``xy_plane`` on axis-aligned grids + is pinned by regression tests for both wake-model paths. + """ + Nx, Ny = Xs.shape + if hasattr(wake_model, "_algo"): # foxes coupling + import foxes.variables as FV + from foxes import Engine + from foxes.utils import wd2uv + + locations = np.stack( + [np.ravel(Xs), np.ravel(Ys), np.full(Xs.size, float(z))], axis=1 + ) + with Engine.new(**wake_model._engine_pars): + point_results = wake_model._algo.calc_points( + wake_model._farm_results, + locations[None], + outputs=[FV.AMB_WS, FV.WS, FV.WD], + ) + amb_uv = wd2uv( + point_results[FV.WD].to_numpy()[0], point_results[FV.AMB_WS].to_numpy()[0] + ) + uv = wd2uv( + point_results[FV.WD].to_numpy()[0], point_results[FV.WS].to_numpy()[0] + ) + amb_uv = amb_uv.reshape(Nx, Ny, 2) + uv = uv.reshape(Nx, Ny, 2) + return amb_uv[..., 0], amb_uv[..., 1], uv[..., 0], uv[..., 1] + + # Lanzilao (UniDirectionalSelfSimilar) path + from wayve.forcing.wind_farms.wake_model_coupling.wake_models.lanzilao_merging import ( + array_of_matrices, + dot_matrix_vec_arrays, + ) + from wayve.forcing.wind_farms.wake_model_coupling.wake_models.wake_model_tools import ( + evaluate_TI, + ) + + Nz = 1 + zs = np.array([z]) + # Get the turbine thrust coefficients + _, Ct, _ = wake_model.get_St_Ct_et(wind_farm, abl, u_bg_evaluator, apm_evaluator) + # Get wind farm information + turbines = wind_farm.turbines + Nt = wind_farm.Nturb + xloc = np.array([turbines[k].x for k in range(Nt)]) + yloc = np.array([turbines[k].y for k in range(Nt)]) + # Ambient TI + TI_inf = abl.TI + # Get turbine direction (streamwise) + e_str, e_span = wake_model.background_flow_direction(wind_farm, abl) + theta_str = np.arctan2(e_str[1], e_str[0]) + # Sort turbines along wind direction + order = wake_model.sort_turbines(wind_farm, e_str) + # Turbine direction evaluation + if wake_model.wake_deflection: # Base turbine direction on APM velocity + u_1, v_1, h_1 = apm_evaluator(xloc, yloc) + theta_turb = np.array([np.arctan2(v_1[i], u_1[i]) for i in range(Nt)]) + else: + theta_turb = np.array([theta_str for _ in range(Nt)]) + # Vector normal (t) and parallel (p) to rotor + ets = np.array( + [np.array([np.cos(theta_turb[i]), np.sin(theta_turb[i])]) for i in order] + ) + eps = np.array( + [np.array([-np.sin(theta_turb[i]), np.cos(theta_turb[i])]) for i in order] + ) + # Sort turbines along wind direction - for TI + xloc_sort = np.array([xloc[i] for i in order]) + yloc_sort = np.array([yloc[i] for i in order]) + D_sort = np.array([turbines[i].D for i in order]) + Ct_sort = Ct[order] + theta_turb_sort = np.array([theta_turb[i] for i in order]) + et_sort = np.array([ets[i, :] for i in order]) + ep_sort = np.array([eps[i, :] for i in order]) + # Evaluate TI at turbine locations + TI = evaluate_TI( + Nt, + et_sort, + ep_sort, + order, + xloc_sort, + yloc_sort, + D_sort, + Ct_sort, + TI_inf, + theta_turb_sort, + wake_model.ka, + wake_model.kb, + ) + # Set up output + if wake_model.wake_deflection: + W_comb = np.zeros((Nx, Ny, Nz, 2)) + e_str[None, None, None, :] + else: + W_comb = np.ones((Nx, Ny, Nz)) + for turb in range(Nt): + if Ct[turb] != 0.0: + # Evaluate W_t + W_t = type(wake_model).wake_function( + Nx, + Ny, + Nz, + Xs, + Ys, + zs, + TI[turb], + Ct[turb], + xloc[turb], + yloc[turb], + turbines[turb].D, + turbines[turb].zh, + theta_turb[turb], + ka=wake_model.ka, + kb=wake_model.kb, + ind=wake_model.induction, + mirr=wake_model.mirrored, + ) + if wake_model.wake_deflection: + # Matrix of current turbine + et = ets[turb] + ep = eps[turb] + A_t = array_of_matrices(1 - W_t, np.outer(et, et)) + np.outer(ep, ep) + # Update wake deficit field + W_comb = dot_matrix_vec_arrays(A_t, W_comb) + else: + W_comb *= 1 - W_t + # Wake deficit in main flow direction + if wake_model.wake_deflection: + wake_deficit = np.inner(e_str, W_comb) + else: + wake_deficit = W_comb + # Evaluate background velocities + locations = np.stack( + [np.ravel(Xs), np.ravel(Ys), np.full(Xs.size, float(z))], axis=1 + ) + vel_bg = u_bg_evaluator(locations) + u_bg = np.reshape(vel_bg[:, 0], (Nx, Ny, Nz)) + v_bg = np.reshape(vel_bg[:, 1], (Nx, Ny, Nz)) + # Split up into streamwise and spanwise components + str_bg = u_bg * e_str[0] + v_bg * e_str[1] + span_bg = u_bg * e_span[0] + v_bg * e_span[1] + # Add wakes + str_wm = np.multiply(str_bg, wake_deficit) + # Convert to u and v components + u_wm = str_wm * e_str[0] + span_bg * e_span[0] + v_wm = str_wm * e_str[1] + span_bg * e_span[1] + return u_bg[:, :, 0], v_bg[:, :, 0], u_wm[:, :, 0], v_wm[:, :, 0] + + def nieuwstadt83_profiles(zh, v, wd, z0=1.0e-1, h=1.5e3, fc=1.0e-4, ust=0.666): """Set up the cubic analytical profile from Nieuwstadt (1983), based on hub height velocity information""" import mpmath @@ -523,7 +838,7 @@ def read_turbine_type(turb_dat): return hh, rd, ct_curve, cp_curve -def wf_setup(farm_dat, analysis_dat, L_filter=1.0e3, debug_mode=False): +def wf_setup(farm_dat, analysis_dat, L_filter=1.0e3, debug_mode=False, rotation=0.0): # WAYVE imports from wayve.forcing.apm_forcing import ForcingComposite from wayve.forcing.wind_farms.dispersive_stresses import DispersiveStresses @@ -533,14 +848,21 @@ def wf_setup(farm_dat, analysis_dat, L_filter=1.0e3, debug_mode=False): #################### # Set up WindFarm object #################### - # Get x and y positions - x = farm_dat["layouts"][0]["coordinates"]["x"] - y = farm_dat["layouts"][0]["coordinates"]["y"] + # Get x and y positions. Copy: this is called once per state now (the + # solver-frame rotation differs per state), so farm_dat must stay pristine. + x = np.asarray(farm_dat["layouts"][0]["coordinates"]["x"], dtype=float) + y = np.asarray(farm_dat["layouts"][0]["coordinates"]["y"], dtype=float) # Reposition to be at grid center wf_offset_x = np.mean(x) wf_offset_y = np.mean(y) - x -= wf_offset_x - y -= wf_offset_y + x = x - wf_offset_x + y = y - wf_offset_y + # Rotate the layout into the solver frame (see flow_io_abl: the ABL is + # built with the hub-height wind along +x, wayve's westerly convention, + # so the layout must turn with it to keep the wind-relative geometry). + if rotation != 0.0: + c, s = np.cos(rotation), np.sin(rotation) + x, y = c * x + s * y, c * y - s * x # Number of turbines Nt = len(x) # Get turbine types @@ -648,27 +970,32 @@ def wake_model_setup(analysis_dat, debug_mode=False): Lanzilao, ) - # WM tool - wake_tool = analysis_dat.get( - "wake_tool", "wayve" - ) # updated by Jonas -TODO update this according to updated schema + # WM tool. The windIO schema places `wake_tool` inside `wm_coupling`; a + # top-level `analysis.wake_tool` is rejected by validation (the schema sets + # additionalProperties: false). Older hand-written yamls put it at the + # analysis level, so keep reading that as a fallback. + wake_tool = analysis_dat.get("wm_coupling", {}).get( + "wake_tool", analysis_dat.get("wake_tool", "wayve") + ) if wake_tool == "wayve": # Read wake model settings # wm_dat = analysis_dat["wind_deficit_model"] k_dat = wm_dat["wake_expansion_coefficient"] - # k, k_a, k_b, ceps + # windIO convention: k = k_a + k_b * TI. wayve's Lanzilao computes + # kwake = ka * TI + kb, so the pair maps swapped: ka=k_b, kb=k_a. + # A scalar k is a constant expansion -> kb, with no TI term. if "k_a" in k_dat and "k_b" in k_dat and "ceps" in wm_dat: - k_a = k_dat["k_a"] - k_b = k_dat["k_b"] + ti_coef = k_dat["k_b"] + k_const = k_dat["k_a"] ceps = wm_dat["ceps"] elif "k" in k_dat and "ceps" in wm_dat: - k_a = k_dat["k"] - k_b = 0.0 + ti_coef = 0.0 + k_const = k_dat["k"] ceps = wm_dat["ceps"] else: raise ValueError("Wake spreading parameter not specified!") # Use wake merging method of Lanzilao and Meyers (2021) - wake_model = Lanzilao(ka=k_a, kb=k_b, eps_beta=ceps) + wake_model = Lanzilao(ka=ti_coef, kb=k_const, eps_beta=ceps) elif wake_tool == "foxes": require("foxes") from foxes import ModelBook @@ -678,26 +1005,159 @@ def wake_model_setup(analysis_dat, debug_mode=False): verbosity = 1 if debug_mode else 0 + # foxes' Dict takes its own label as `_name`; a plain `name=` kwarg is + # stored as a *data* key and would be forwarded into the foxes + # Algorithm constructor (TypeError: unexpected keyword argument 'name'). algo_dict = Dict( algo_type="Downwind", wake_models=[], verbosity=verbosity, - name="wayve.algorithm", + _name="wayve.algorithm", ) - ana_dict = Dict(analysis_dat, name="analysis") - idict = Dict(algorithm=algo_dict, name="wayve") + ana_dict = Dict(analysis_dat, _name="analysis") + idict = Dict(algorithm=algo_dict, _name="wayve") mbook = ModelBook() _read_analysis(ana_dict, idict, mbook=mbook, verbosity=verbosity) wake_model = FoxesWakeModel(mbook=mbook, **idict["algorithm"]) + # wayve@e87780a overrides background_flow_direction with logic + # identical to the base class but a broken lazy import (e_spanwise + # from wake_models.wake_model_tools instead of forcing_tools), which + # crashes every foxes-coupled solve. Rebind the base implementation + # while the override still carries the broken import; probing the + # override's own source (rather than the import target) keeps the + # shim inert once upstream fixes or rewrites the method. + import inspect + + try: + bfd_src = inspect.getsource(FoxesWakeModel.background_flow_direction) + except (OSError, TypeError): + bfd_src = "" + if "wake_model_tools import e_spanwise" in bfd_src: + import types + + from wayve.forcing.wind_farms.wake_model_coupling.wake_model_interface import ( + UniDirectionalSelfSimilar, + ) + + wake_model.background_flow_direction = types.MethodType( + UniDirectionalSelfSimilar.background_flow_direction, wake_model + ) else: raise NotImplementedError(f"Wake tool '{wake_tool}' not implemented!") return wake_model -def flow_io_abl(wind_resource_dat, time_index, zh, h1, dh_max=None, serz=True): +# Capping inversion: the windIO energy-resource schema spells the four scalars +# flat, and WIFA's own code_saturne adapter already reads them that way +# (cs_launch_modules.py). The nested +# ``thermal_stratification.capping_inversion`` block is the older spelling and +# is still accepted, so files written against it keep working; it cannot be +# expressed in a netCDF ``!include``, whose root group is flat. +# Ordered (h, dH, dtheta, lapse_rate) — the tuple read_capping_inversion returns. +_CI_FLAT_KEYS = ( + "ABL_height", + "capping_inversion_thickness", + "capping_inversion_strength", + "lapse_rate", +) +_CI_NESTED_KEYS = ("ABL_height", "dH", "dtheta", "lapse_rate") + + +def capping_inversion_spelling(wind_resource_dat): + """Which spelling this wind resource uses for the capping inversion. + + Parameters + ---------- + wind_resource_dat: dict + Wind resource data + + Returns + ------- + str or None + ``"flat"``, ``"nested"``, or None when no inversion is specified (the + caller then fits one from the potential-temperature profile, or falls + back to defaults). + + Raises + ------ + ValueError + If both spellings are present, or if either is incomplete. Both are + defects of the *file*, not of a single state, which is why this check + is separate from the per-state value read: ``run_wayve`` calls it once + before the state loop, whose ``except Exception`` would otherwise turn + a malformed file into a run where every state silently "crashed". + + A partial block is rejected rather than completed from defaults: of the + four scalars, ``capping_inversion_strength`` alone sets the amplitude of + the gravity-wave forcing, so silently pairing a real ``ABL_height`` with + an invented dtheta yields a result that looks fitted and is not. + """ + flat = [key for key in _CI_FLAT_KEYS if key in wind_resource_dat] + nested = wind_resource_dat.get("thermal_stratification", {}).get( + "capping_inversion" + ) + if flat and nested is not None: + raise ValueError( + "Wind resource specifies the capping inversion twice: flat keys " + f"{sorted(flat)} and a nested thermal_stratification." + "capping_inversion block. Keep one (the flat keys are the windIO " + "schema spelling)." + ) + if nested is not None: + missing = [key for key in _CI_NESTED_KEYS if key not in nested] + if missing: + raise ValueError( + "Incomplete thermal_stratification.capping_inversion block; " + f"missing {missing}. All of {list(_CI_NESTED_KEYS)} are required." + ) + return "nested" + if not flat: + return None + if len(flat) != len(_CI_FLAT_KEYS): + raise ValueError( + f"Incomplete capping inversion: found {sorted(flat)}, missing " + f"{sorted(set(_CI_FLAT_KEYS) - set(flat))}. All four are required " + "(they are not completed from defaults; see " + "capping_inversion_spelling). Omit all four to fit the inversion " + "from the potential-temperature profile instead." + ) + return "flat" + + +def read_capping_inversion(wind_resource_dat, time_index): + """The capping inversion of one state, in either windIO spelling. + + Parameters + ---------- + wind_resource_dat: dict + Wind resource data + time_index: int + Index of the timestamp to read + + Returns + ------- + tuple or None + ``(h, dh, dth, dthdz)`` — inversion centre height [m], thickness [m], + strength [K] and free-atmosphere lapse rate [K/m] — or None when the + resource specifies no inversion. + """ + spelling = capping_inversion_spelling(wind_resource_dat) + if spelling is None: + return None + if spelling == "nested": + block = wind_resource_dat["thermal_stratification"]["capping_inversion"] + return tuple(float(block[key]["data"][time_index]) for key in _CI_NESTED_KEYS) + return tuple( + float(wind_resource_dat[key]["data"][time_index]) for key in _CI_FLAT_KEYS + ) + + +def flow_io_abl( + wind_resource_dat, time_index, zh, h1, dh_max=None, serz=True, gmode="avg" +): """ Method to set up an ABL object based on FLOW IO @@ -715,6 +1175,29 @@ def flow_io_abl(wind_resource_dat, time_index, zh, h1, dh_max=None, serz=True): Maximum depth of the inversion layer used in the inversion curve fitting procedure (default: None) serz (optional): boolean Whether the surface-extended version of the RZ model is used (default: True) + gmode (optional): str + How the free-atmosphere (geostrophic) velocity is derived from + height-resolved wind resources with mesoscale surface scalars + (see wayve's ``abl_setup.mesoscale_based``): "h1" (at the inversion + center), "h2" (at the inversion top), "avg" (profile average between + the inversion and 5 km; needs profiles reaching 5 km), or "trop" + (average to the tropopause; needs profiles reaching the stratosphere). + Ignored by the scalar and turbulence-profile input paths. + (default: "avg") + + Returns + ------- + abl: wayve.abl.abl.ABL + Atmospheric state in the *solver frame*: the hub-height wind blows + along +x (wayve's westerly convention, which its gravity-wave + machinery and anisotropic grids assume). Vertical veer is preserved + as spanwise components. + rotation: float + Angle (radians, counterclockwise) of the earth-frame hub-height flow + vector measured from east. Rotating earth-frame coordinates by + ``-rotation`` maps them into the solver frame (``wf_setup`` does this + to the layout); rotating solver-frame vectors by ``+rotation`` maps + results back to earth. """ # Atmospheric state setup from wayve.abl.abl import ABL @@ -724,7 +1207,14 @@ def flow_io_abl(wind_resource_dat, time_index, zh, h1, dh_max=None, serz=True): kappa = 0.41 # Von Karman constant omega = 7.2921159e-5 # angular speed of the Earth [rad/s] # Basic atmospheric scalars # - air_density = 1.225 # Hard-coded for now + # Operating air density. The APM velocity solution does not depend on it + # (see wayve's power_turbines docstring); it linearly scales the reported + # turbine power — the Cp curve itself stays referenced to the standard + # 1.225 kg/m3 of the windIO power curve (read_turbine_type) — and is + # forwarded to foxes as FV.RHO by the foxes coupling. + air_density = 1.225 + if "air_density" in wind_resource_dat.keys(): + air_density = wind_resource_dat["air_density"]["data"][time_index] # Surface roughness z0 = 1.0e-1 if "z0" in wind_resource_dat.keys(): @@ -744,28 +1234,30 @@ def flow_io_abl(wind_resource_dat, time_index, zh, h1, dh_max=None, serz=True): # Wind speed and direction v = wind_resource_dat["wind_speed"]["data"][time_index] wd = wind_resource_dat["wind_direction"]["data"][time_index] + # Solver-frame normalization: build the profile as if the wind were + # westerly (the Ekman veer relative to the hub-height direction is + # unchanged) and report the rotation undone by doing so. + rotation = np.deg2rad(270.0 - wd) + wd = 270.0 # Friction velocity ust = 0.666 if "friction_velocity" in wind_resource_dat.keys(): ust = wind_resource_dat["friction_velocity"]["data"][time_index] - # Turbulence intensity + # Turbulence intensity. windIO carries TI as a fraction (as does the + # vertical-profile branch below, and the 0.04 default just above), so + # it is used as-is. Guard on the variable actually being read. TI = 0.04 - if "z0" in wind_resource_dat.keys(): - TI = wind_resource_dat["turbulence_intensity"]["data"][time_index] / 100.0 + if "turbulence_intensity" in wind_resource_dat.keys(): + TI = wind_resource_dat["turbulence_intensity"]["data"][time_index] # Capping inversion information h = 1.5e3 dh = 100.0 dth = 5.0 dthdz = 2.0e-3 th0 = 293.15 - if "thermal_stratification" in wind_resource_dat.keys(): - thermal_data = wind_resource_dat["thermal_stratification"] - if "capping_inversion" in thermal_data.keys(): - ci_data = thermal_data["capping_inversion"] - h = ci_data["ABL_height"]["data"][time_index] - dh = ci_data["dH"]["data"][time_index] - dth = ci_data["dtheta"]["data"][time_index] - dthdz = ci_data["lapse_rate"]["data"][time_index] + ci = read_capping_inversion(wind_resource_dat, time_index) + if ci is not None: + h, dh, dth, dthdz = ci inv_bottom, inv_top = h - dh / 2, h + dh / 2 # Nieuwstadt profiles for velocity and shear stress zs, us, vs, U3, V3, tauxs, tauys, nus = nieuwstadt83_profiles( @@ -776,63 +1268,201 @@ def flow_io_abl(wind_resource_dat, time_index, zh, h1, dh_max=None, serz=True): else: # Read out vertical profile zs = np.array(wind_resource_dat["height"]) + if not np.all(np.diff(zs) > 0): + # np.interp silently returns garbage on unsorted abscissae, and + # the hub interpolation below now derives the whole solver frame. + raise UserWarning("wind resource 'height' must be strictly ascending") vs = np.array(wind_resource_dat["wind_speed"]["data"][time_index]) wds = np.array(wind_resource_dat["wind_direction"]["data"][time_index]) ths = np.array(wind_resource_dat["potential_temperature"]["data"][time_index]) - TIs = np.array(wind_resource_dat["turbulence_intensity"]["data"][time_index]) - # Interpolate TI - TI = np.interp(zh, zs, TIs) + # Solver-frame normalization: shift the direction profile so the + # hub-height wind is westerly (flow along +x). Shifting the input + # keeps the veer and makes everything derived below (velocity + # components, Gmode geostrophic wind, momentum-flux alignment) land + # in the solver frame consistently — without wayve's ABL.rotate(), + # whose momentum-flux rotation is broken in wayve 2.0.0. Unwrap + # first so interpolation across the 0/360 seam is safe. + wds = np.rad2deg(np.unwrap(np.deg2rad(wds))) + wd_hub = np.interp(zh, zs, wds) + rotation = np.deg2rad(270.0 - wd_hub) + wds = wds + (270.0 - wd_hub) + # Turbulence intensity: height-resolved profile or one value per state + TIs = np.atleast_1d( + np.array(wind_resource_dat["turbulence_intensity"]["data"][time_index]) + ) + TI = np.interp(zh, zs, TIs) if TIs.size > 1 else float(TIs[0]) # Velocity components us = -vs * np.sin(np.deg2rad(wds)) vs = -vs * np.cos(np.deg2rad(wds)) # Check available inputs - if "k" in wind_resource_dat.keys(): # RANS-like inputs - tkes = np.array(wind_resource_dat["k"]["data"][time_index]) - eps = np.array(wind_resource_dat["epsilon"]["data"][time_index]) - # Eddy viscosity - C_mu = 0.09 # k-epsilon model value - nus = C_mu * np.divide( - np.square(tkes), eps, out=np.zeros_like(tkes), where=eps != 0 - ) - # Momentum fluxes - dudz = np.gradient(us, zs, edge_order=2) - dvdz = np.gradient(vs, zs, edge_order=2) - tauxs = nus * dudz - tauys = nus * dvdz - else: # Shear stress profile directly available - tauxs = np.array(wind_resource_dat["tau_x"]["data"][time_index]) - tauys = np.array(wind_resource_dat["tau_y"]["data"][time_index]) - nus = None - # Total momentum flux - taus = np.sqrt(np.square(tauxs) + np.square(tauys)) - # Friction velocity - ust = taus[0] # Assume friction velocity is not given explicitly - # Estimate boundary layer height based on momentum flux # - f_tau = interp1d(taus, zs) - blh = f_tau(0.01 * ust) - # Capping inversion information - if ( - "thermal_stratification" in wind_resource_dat.keys() - and "capping_inversion" - in wind_resource_dat["thermal_stratification"].keys() + if "k" in wind_resource_dat.keys() or "tau_x" in wind_resource_dat.keys(): + # Turbulence profiles provided directly (LES/RANS-like inputs) + if "k" in wind_resource_dat.keys(): # RANS-like inputs + tkes = np.array(wind_resource_dat["k"]["data"][time_index]) + eps = np.array(wind_resource_dat["epsilon"]["data"][time_index]) + # Eddy viscosity + C_mu = 0.09 # k-epsilon model value + nus = C_mu * np.divide( + np.square(tkes), eps, out=np.zeros_like(tkes), where=eps != 0 + ) + # Momentum fluxes + dudz = np.gradient(us, zs, edge_order=2) + dvdz = np.gradient(vs, zs, edge_order=2) + tauxs = nus * dudz + tauys = nus * dvdz + else: # Shear stress profile directly available + tauxs = np.array(wind_resource_dat["tau_x"]["data"][time_index]) + tauys = np.array(wind_resource_dat["tau_y"]["data"][time_index]) + # Explicit stress components are earth-frame vectors: rotate + # them into the solver frame along with the velocities. + c_r, s_r = np.cos(rotation), np.sin(rotation) + tauxs, tauys = c_r * tauxs + s_r * tauys, c_r * tauys - s_r * tauxs + nus = None + # Total momentum flux + taus = np.sqrt(np.square(tauxs) + np.square(tauys)) + # Friction velocity + ust = taus[0] # Assume friction velocity is not given explicitly + # Estimate boundary layer height based on momentum flux # + f_tau = interp1d(taus, zs) + blh = f_tau(0.01 * ust) + # Capping inversion information + ci = read_capping_inversion(wind_resource_dat, time_index) + if ci is not None: + h, dh, dth, dthdz = ci + # Mixed-layer temperature from the profile the resource + # supplies, as every other branch does. A fixed 293.15 K here + # made the buoyancy jump g*dtheta/th0 disagree with the ABL's + # own theta profile by up to ~4% at realistic temperatures. + th0 = np.interp(h, zs, ths) + inv_bottom, inv_top = h - dh / 2, h + dh / 2 + else: + inv_bottom, h, inv_top, th0, dth, dthdz = ci_fitting( + zs, ths, l_mo, blh, dh_max=dh_max, serz=serz + ) + # Geostrophic wind speed + z = np.linspace(h, 15.0e3, 1000) + _trapezoid = np.trapezoid if hasattr(np, "trapezoid") else np.trapz + U3 = _trapezoid(np.interp(z, zs, us), z) / (15.0e3 - h) + V3 = _trapezoid(np.interp(z, zs, vs), z) / (15.0e3 - h) + elif ( + "friction_velocity" in wind_resource_dat.keys() + and "boundary_layer_height" in wind_resource_dat.keys() ): - thermal_data = wind_resource_dat["thermal_stratification"] - ci_data = thermal_data["capping_inversion"] - th0 = 293.15 - h = ci_data["ABL_height"]["data"][time_index] - dh = ci_data["dH"]["data"][time_index] - dth = ci_data["dtheta"]["data"][time_index] - dthdz = ci_data["lapse_rate"]["data"][time_index] - inv_bottom, inv_top = h - dh / 2, h + dh / 2 + # Mesoscale/reanalysis-style inputs (e.g. ERA5): u/v/theta profiles + # plus surface scalars (ust, blh, L_MO) instead of turbulence + # profiles — wayve's abl_setup.mesoscale_based() territory. + ust = wind_resource_dat["friction_velocity"]["data"][time_index] + blh = wind_resource_dat["boundary_layer_height"]["data"][time_index] + if zs[-1] >= 12.0e3: + # Profiles reach the stratosphere: use wayve's full setup, + # including the tropopause two-line fit. Note the radiating + # top-layer state differs from the truncated branch below: + # here Uinf/Vinf are means above the fitted tropopause and + # Ninf is the fitted stratospheric N, while the truncated + # branch uses the top-of-profile point and the inversion + # fit's free lapse rate. + from wayve.abl.abl_setup import mesoscale_based + + lat = np.rad2deg(np.arcsin(fc / (2.0 * omega))) + # us/vs are already solver-frame (see above), so the ABL + # comes out westerly-aligned without further rotation. + return ( + mesoscale_based( + zs, + us, + vs, + ths, + ust, + blh, + l_mo, + lat, + h1, + z0=(z0 if "z0" in wind_resource_dat.keys() else None), + TI=TI, + rho=air_density, + dh_max=(300.0 if dh_max is None else dh_max), + Gmode=gmode, + serz=serz, + ), + rotation, + ) + # Truncated profiles (reanalysis subsets often stop below the + # tropopause, e.g. ERA5 levels 96-137 end near 6 km): run the same + # core setup but skip mesoscale_based's unconditional tropopause + # fit, which would feed a garbage stratosphere into the ABL. The + # free atmosphere is capped at the profile top instead (see the + # ABL constructor below), and Gmode is restricted to what the + # data supports. + stable = 0.0 < l_mo < 100 + # Capping inversion: an explicit windIO block wins over the fit + ci = read_capping_inversion(wind_resource_dat, time_index) + if ci is not None: + h, dh, dth, dthdz = ci + inv_bottom, inv_top = h - dh / 2, h + dh / 2 + th0 = np.interp(h, zs, ths) + else: + inv_bottom, h, inv_top, th0, dth, dthdz = ci_fitting( + zs, + ths, + l_mo, + blh, + dh_max=(300.0 if dh_max is None else dh_max), + serz=serz, + ) + # Momentum flux and eddy viscosity profiles from surface scalars + # (same closure as mesoscale_based: stable turbulence extends to + # blh, convective turbulence to the inversion) + tau = np.zeros(zs.shape) + nus = np.zeros(zs.shape) + if stable: + m = zs <= blh + tau[m] = ust**2 * (1 - zs[m] / blh) ** 1.5 + nus[m] = kappa * ust * zs[m] * (1 - zs[m] / blh) ** 2 + else: + m = zs <= h + tau[m] = ust**2 * (1 - zs[m] / h) + nus[m] = kappa * ust * zs[m] * (1 - zs[m] / h) ** 2 + # Geostrophic wind from the profile, restricted to the data range + if gmode == "trop": + raise UserWarning( + "Gmode 'trop' needs profiles reaching the stratosphere; " + f"these stop at {zs[-1]:.0f} m" + ) + if gmode == "avg" and zs[-1] < 5.0e3: + warnings.warn( + f"Gmode 'avg' averages the wind up to 5 km but the " + f"profiles stop at {zs[-1]:.0f} m; falling back to 'h2' " + "(wind at the inversion top)" + ) + gmode = "h2" + if gmode == "h1": + U3 = np.interp(h, zs, us) + V3 = np.interp(h, zs, vs) + elif gmode == "h2": + U3 = np.interp(inv_top, zs, us) + V3 = np.interp(inv_top, zs, vs) + elif gmode == "avg": + z = np.linspace(h, 5.0e3, 1000) + _trapezoid = np.trapezoid if hasattr(np, "trapezoid") else np.trapz + U3 = _trapezoid(np.interp(z, zs, us), z) / (5.0e3 - h) + V3 = _trapezoid(np.interp(z, zs, vs), z) / (5.0e3 - h) + else: + raise UserWarning(f"Gmode '{gmode}' unknown") + # Momentum flux aligned with the geostrophic wind (as in + # mesoscale_based) + tau_angle = np.arctan2(V3, U3) + tauxs = np.cos(tau_angle) * tau + tauys = np.sin(tau_angle) * tau + # Log-law surface roughness estimate when not provided + if "z0" not in wind_resource_dat.keys(): + z0 = zs[0] / np.exp(kappa * np.hypot(us[0], vs[0]) / ust) else: - inv_bottom, h, inv_top, th0, dth, dthdz = ci_fitting( - zs, ths, l_mo, blh, dh_max=dh_max, serz=serz + raise UserWarning( + "Vertical-profile wind resource needs either turbulence " + "profiles ('k'/'epsilon' or 'tau_x'/'tau_y') or mesoscale " + "surface scalars ('friction_velocity' and " + "'boundary_layer_height')" ) - # Geostrophic wind speed - z = np.linspace(h, 15.0e3, 1000) - _trapezoid = np.trapezoid if hasattr(np, "trapezoid") else np.trapz - U3 = _trapezoid(np.interp(z, zs, us), z) / (15.0e3 - h) - V3 = _trapezoid(np.interp(z, zs, vs), z) / (15.0e3 - h) # Upper layer thickness h2 = h - h1 if ( @@ -845,28 +1475,40 @@ def flow_io_abl(wind_resource_dat, time_index, zh, h1, dh_max=None, serz=True): # gprime and N gprime = gravity * dth / th0 N = np.sqrt(gravity * dthdz / th0) - # Set up ABL object - return ABL( - zs, - us, - vs, - ths, - tauxs, - tauys, - h1, - h2, - gprime, - N, - U3, - V3, - fc, - nus=nus, - rho=air_density, - TI=TI, - z0=z0, - ust=ust, - inv_bottom=inv_bottom, - inv_top=inv_top, + # Set up ABL object. The free atmosphere is capped at the top of the + # profile data: NonUniform's layer setup spans [inversion top, h_strat], + # and the default h_strat of 10 km would put most of its layers in a + # no-data region where the profile splines just hold constants. Above + # the cap sits the semi-infinite radiating layer with the top-of-profile + # state (Uinf/Vinf) and the fitted free-atmosphere stratification (Ninf). + return ( + ABL( + zs, + us, + vs, + ths, + tauxs, + tauys, + h1, + h2, + gprime, + N, + U3, + V3, + fc, + nus=nus, + rho=air_density, + TI=TI, + z0=z0, + ust=ust, + inv_bottom=inv_bottom, + inv_top=inv_top, + h_strat=min(10.0e3, zs[-1]), + Uinf=us[-1], + Vinf=vs[-1], + Ninf=N, + ), + rotation, )