From 1885a20cde3d35c1f37d5f1c58fc6a479f30a52f Mon Sep 17 00:00:00 2001
From: lejeunemax <max-lej@hotmail.com>
Date: Tue, 16 Dec 2025 15:47:01 +0100
Subject: [PATCH 1/7] Adding a coupling with Floris

---
 wifa/floris_api.py | 409 +++++++++++++++++++++++++++++++++++++++++++++
 wifa/main_api.py   |  11 +-
 2 files changed, 418 insertions(+), 2 deletions(-)
 create mode 100644 wifa/floris_api.py

diff --git a/wifa/floris_api.py b/wifa/floris_api.py
new file mode 100644
index 0000000..9a0ad93
--- /dev/null
+++ b/wifa/floris_api.py
@@ -0,0 +1,409 @@
+from windIO import load_yaml
+from typing import Dict, Any, List
+import numpy as np
+import windIO
+from pathlib import Path
+
+from typing import TYPE_CHECKING
+if TYPE_CHECKING:
+    from floris import FlorisModel
+
+def run_floris(yaml_input):
+
+    from floris import FlorisModel
+    from floris.read_windio import TrackedDict
+
+    if isinstance(yaml_input, (str, Path)):
+        windio_dict = windIO.load_yaml(yaml_input)
+    elif isinstance(yaml_input, dict):
+        windio_dict = yaml_input
+    else:
+        raise TypeError(f"yaml_input must be a file path or dict, got {type(yaml_input)}")
+
+    fmodel = FlorisModel.from_windio(windio_dict)
+    fmodel.run()
+
+    windio_dict['_context'] = ''
+    
+    with TrackedDict.from_parent(windio_dict, 'attributes') as attrs:
+        attrs.untrack('analysis')
+        attrs.untrack('flow_model')
+
+        out_specs = attrs['model_outputs_specification']
+        run_config = out_specs['run_configuration']
+
+        slice_selection = None
+        times_run: Dict = run_config.get('times_run', None)
+        wind_speeds_run: Dict = run_config.get('wind_speeds_run', None)
+        directions_run: Dict = run_config.get('directions_run', None)
+
+        if times_run:
+            all_occurences = times_run.get("all_occurences", False)
+            
+            if all_occurences:
+                slice_selection = slice(None)
+            else:
+                slice_selection = times_run.get("subset", None)
+
+        elif wind_speeds_run and directions_run:
+            
+            if not wind_speeds_run.get("all_values", False):
+                raise NotImplementedError(
+                    f"Wind speed and direction subsets are not yet supported, got {wind_speeds_run.name} {wind_speeds_run}"
+                )
+            if not directions_run.get("all_values", False):
+                raise NotImplementedError(
+                    f"Wind speed and direction subsets are not yet supported, got {directions_run.name} {directions_run}"
+                )
+        
+        output_dir = out_specs.get("output_folder", ".")
+
+        if ("turbine_outputs" in out_specs):
+            _read_turbines(fmodel, output_dir, out_specs, slice_selection)
+
+        if ("flow_field" in out_specs):
+            _read_fields(fmodel, output_dir, out_specs, slice_selection)
+        
+        pass
+        
+def _read_turbines(fmodel: "FlorisModel", output_dir, out_specs, slice_selection):
+    """
+    Extract turbine output data from FLORIS model and write to NetCDF file.
+    
+    Args:
+        fmodel: FlorisModel instance that has been run
+        output_dir: Output directory path
+        out_specs: TrackedDict containing turbine_outputs configuration from out_specs["turbine_outputs"]
+        slice_selection: Time slice selection for filtering output
+    
+    Returns:
+        xarray.Dataset containing turbine output data
+    """
+    import xarray as xr 
+
+    turbine_outputs = out_specs["turbine_outputs"]
+
+    output_vars: List[str] = turbine_outputs["output_variables"]
+    turbine_nc_filename = turbine_outputs.get("turbine_nc_filename", "turbine_outputs.nc")
+        
+    # Make a copy to avoid modifying the original list
+    requested_vars = output_vars.copy()
+    data_vars = {}
+    coords = {}
+    
+    # Determine the wind data type to understand output structure
+    wind_data_type = None
+    if fmodel.wind_data is not None:
+        from floris.wind_data import WindRose, WindRoseWRG, WindTIRose
+        if isinstance(fmodel.wind_data, WindTIRose):
+            wind_data_type = 'wind_ti_rose'
+        elif isinstance(fmodel.wind_data, (WindRose, WindRoseWRG)):
+            wind_data_type = 'wind_rose'
+    
+    # Define handlers for each variable type
+    # Format: variable_name -> (extractor_func, is_coord)
+    variable_handlers = {
+        'time': (lambda: fmodel.get_metadata('wind_resource', {}).get('time'), True),
+        'x': (lambda: fmodel.layout_x, False),
+        'y': (lambda: fmodel.layout_y, False),
+        'z': (lambda: fmodel.get_turbine_layout(z=True)[2], False),
+        'power': (lambda: fmodel.get_turbine_powers(), False),
+        'thrust_coefficient': (lambda: fmodel.get_turbine_thrust_coefficients(), False),
+        'axial_induction': (lambda: fmodel.get_turbine_ais(), False),
+        'turbulence_intensity': (lambda: fmodel.get_turbine_TIs(), False),
+    }
+    
+    # Process each requested variable
+    unhandled_vars = []
+    for var_name in requested_vars:
+        if var_name in variable_handlers:
+            try:
+                extractor, is_coord = variable_handlers[var_name]
+                data = extractor()
+                
+                if data is not None:
+                    data = np.asarray(data)
+                    
+                    # Determine dimensions based on data shape
+                    if is_coord:
+                        # Coordinate variables are 1D
+                        if data.ndim == 0:
+                            data = np.array([data])
+                        coords[var_name] = data
+                    else:
+                        # Data variables - infer dimensions based on shape and wind data type
+                        if data.ndim == 1:
+                            # 1D: turbine dimension only
+                            data_vars[var_name] = (['turbine'], data)
+                        elif data.ndim == 2:
+                            # 2D: (time, turbine) - TimeSeries or flattened
+                            data_vars[var_name] = (['time', 'turbine'], data)
+                        elif data.ndim == 3:
+                            # 3D: (wind_direction, wind_speed, turbine) - WindRose
+                            data_vars[var_name] = (['wind_direction', 'wind_speed', 'turbine'], data)
+                        elif data.ndim == 4:
+                            # 4D: (wind_direction, wind_speed, turbulence_intensity, turbine) - WindTIRose
+                            data_vars[var_name] = (['wind_direction', 'wind_speed', 'turbulence_intensity', 'turbine'], data)
+                        else:
+                            print(f"Warning: Unexpected data shape for '{var_name}': {data.shape}; skipping.")
+                            continue
+                else:
+                    print(f"Warning: '{var_name}' data not available; skipping.")
+            except Exception as e:
+                print(f"Warning: Failed to extract '{var_name}': {e}; skipping.")
+        else:
+            unhandled_vars.append(var_name)
+    
+    # Warn about unhandled variables
+    if unhandled_vars:
+        print(f"Warning: Unhandled output variables: {unhandled_vars}")
+    
+    # Build coordinates based on data structure
+    final_coords = {}
+    
+    # Add turbine coordinate
+    n_turbines = len(fmodel.layout_x)
+    final_coords['turbine'] = np.arange(n_turbines)
+    
+    # Determine primary dimension and add appropriate coordinates
+    if 'time' in coords and coords['time'] is not None:
+        # TimeSeries case - use time coordinate
+        final_coords['time'] = coords['time']
+        # Rename findex to time in data_vars
+        renamed_data_vars = {}
+        for var_name, (dims, data) in data_vars.items():
+            new_dims = tuple('time' if d == 'findex' else d for d in dims)
+            renamed_data_vars[var_name] = (new_dims, data)
+        data_vars = renamed_data_vars
+    elif wind_data_type in ['wind_rose', 'wind_ti_rose']:
+        # WindRose case - add wind direction and speed coordinates
+        if fmodel.wind_data is not None:
+            final_coords['wind_direction'] = fmodel.wind_data.wind_directions
+            final_coords['wind_speed'] = fmodel.wind_data.wind_speeds
+            if wind_data_type == 'wind_ti_rose':
+                final_coords['turbulence_intensity'] = fmodel.wind_data.turbulence_intensities
+    else:
+        # Fallback: infer findex dimension
+        n_findex = None
+        for var_name, (dims, data) in data_vars.items():
+            if 'findex' in dims:
+                n_findex = data.shape[dims.index('findex')]
+                break
+        if n_findex is not None:
+            final_coords['findex'] = np.arange(n_findex)
+    
+    # Create dataset
+    ds = xr.Dataset(data_vars=data_vars, coords=final_coords)
+    
+    # Add metadata attributes
+    ds.attrs['description'] = 'FLORIS turbine output data'
+    ds.attrs['source'] = 'FLORIS wind farm wake model'
+    
+    # Add units and descriptions for known variables
+    var_attrs = {
+        'x': {'units': 'm', 'long_name': 'turbine x-coordinate'},
+        'y': {'units': 'm', 'long_name': 'turbine y-coordinate'},
+        'z': {'units': 'm', 'long_name': 'turbine hub height'},
+        'power': {'units': 'W', 'long_name': 'turbine power output'},
+        'thrust_coefficient': {'units': '-', 'long_name': 'turbine thrust coefficient'},
+        'axial_induction': {'units': '-', 'long_name': 'turbine axial induction factor'},
+        'turbulence_intensity': {'units': '-', 'long_name': 'turbulence intensity at turbine'},
+    }
+    
+    for var_name, attrs in var_attrs.items():
+        if var_name in ds.data_vars:
+            ds[var_name].attrs.update(attrs)
+    
+    # Apply time selection if specified
+    if slice_selection is not None:
+        # Determine which dimension to slice
+        if 'time' in ds.dims:
+            ds = ds.isel({'time': slice_selection})
+        elif 'findex' in ds.dims:
+            ds = ds.isel({'findex': slice_selection})
+        elif 'wind_direction' in ds.dims:
+            # For WindRose, could slice wind_direction or wind_speed
+            # This would need more sophisticated handling based on user intent
+            print("Warning: slice_selection not applied for WindRose data")
+    
+    # Write to NetCDF
+    ds.to_netcdf(f"{output_dir}/{turbine_nc_filename}")
+    
+    return ds
+    
+
+def _read_fields(fmodel: "FlorisModel", output_dir, out_specs, slice_selection):
+    """
+    Extract flow field data from FLORIS model and write to NetCDF file.
+    
+    Args:
+        fmodel: FlorisModel instance that has been run
+        output_dir: Output directory path
+        out_specs: TrackedDict containing flow_field configuration from out_specs["flow_field"]
+        slice_selection: Time slice selection for filtering output
+    
+    Returns:
+        xarray.Dataset containing flow field data
+    """
+    import xarray as xr
+
+    flow_field = out_specs["flow_field"]
+
+    # Check if flow field reporting is enabled
+    report = flow_field.get("report", True)
+    
+    if not report:
+        return None
+    
+    # Get output configuration
+    output_vars: List[str] = flow_field["output_variables"]
+    flow_nc_filename = flow_field.get("flow_nc_filename", "flow_field.nc")
+    
+    # Get z_planes configuration
+    z_planes = flow_field["z_planes"]
+    z_sampling = z_planes["z_sampling"]
+    xy_sampling = z_planes.get("xy_sampling", "grid")
+    
+    # Determine z values to sample
+    if z_sampling == "hub_heights":
+        z_list = fmodel.core.farm.hub_heights
+        # Get unique hub heights
+        z_list = np.unique(z_list)
+    elif z_sampling == "plane_list":
+        z_list = z_planes["z_list"]
+    else:
+        raise ValueError(f"Unknown z_sampling type: {z_sampling}")
+    
+    # Get xy sampling parameters
+    if xy_sampling == "grid":
+        x_bounds = z_planes["x_bounds"]
+        y_bounds = z_planes["y_bounds"]
+        
+        # Support either Nx/Ny or dx/dy specification
+        if "Nx" in z_planes and "Ny" in z_planes:
+            Nx = z_planes["Nx"]
+            Ny = z_planes["Ny"]
+        elif "dx" in z_planes and "dy" in z_planes:
+            dx = z_planes["dx"]
+            dy = z_planes["dy"]
+            Nx = int((x_bounds[1] - x_bounds[0]) / dx) + 1
+            Ny = int((y_bounds[1] - y_bounds[0]) / dy) + 1
+        else:
+            # Default resolution
+            Nx = 200
+            Ny = 200
+
+    elif xy_sampling == "original_grid":
+        # Use the grid from the model
+        x_bounds = None
+        y_bounds = None
+        Nx = None
+        Ny = None
+    else:
+        raise ValueError(f"Unknown xy_sampling type: {xy_sampling}")
+
+    # Calculate horizontal planes for each z level
+    data_vars = {}
+    
+    # Prepare coordinate arrays
+    n_findex = fmodel.core.flow_field.n_findex
+    if slice_selection is not None:
+        if isinstance(slice_selection, slice):
+            findex_range = range(n_findex)[slice_selection]
+        else:
+            findex_range = slice_selection
+    else:
+        findex_range = range(n_findex)
+    
+    # For each z level, calculate the horizontal plane
+    for z_idx, z_height in enumerate(z_list):
+        # Calculate horizontal plane at this height
+        # Use keep_inertial_frame=True to avoid rotation for WIFA output
+        horizontal_plane = fmodel.calculate_horizontal_plane(
+            height=z_height,
+            x_resolution=Nx,
+            y_resolution=Ny,
+            x_bounds=x_bounds,
+            y_bounds=y_bounds,
+            findex_for_viz=0,  # Will loop over findices if needed
+            keep_inertial_frame=True,
+        )
+        
+        # Extract coordinates from first plane
+        if z_idx == 0:
+            # Get unique x and y coordinates
+            x_coords = np.unique(horizontal_plane.df.x1.values)
+            y_coords = np.unique(horizontal_plane.df.x2.values)
+            
+            # Initialize data arrays
+            if "u" in output_vars:
+                u_data = np.zeros((len(findex_range), len(z_list), len(y_coords), len(x_coords)))
+            if "v" in output_vars:
+                v_data = np.zeros((len(findex_range), len(z_list), len(y_coords), len(x_coords)))
+            if "w" in output_vars:
+                w_data = np.zeros((len(findex_range), len(z_list), len(y_coords), len(x_coords)))
+        
+        # Reshape plane data to grid
+        # The dataframe has x1, x2, x3, u, v, w columns
+        for t_idx, findex in enumerate(findex_range):
+            # Recalculate for this findex if multiple conditions
+            if findex > 0:
+                horizontal_plane = fmodel.calculate_horizontal_plane(
+                    height=z_height,
+                    x_resolution=Nx,
+                    y_resolution=Ny,
+                    x_bounds=x_bounds,
+                    y_bounds=y_bounds,
+                    findex_for_viz=findex,
+                    keep_inertial_frame=True,
+                )
+            
+            # Reshape to 2D grid (y, x)
+            u_grid = horizontal_plane.df.u.values.reshape(len(y_coords), len(x_coords))
+            v_grid = horizontal_plane.df.v.values.reshape(len(y_coords), len(x_coords))
+            w_grid = horizontal_plane.df.w.values.reshape(len(y_coords), len(x_coords))
+            
+            if "u" in output_vars:
+                u_data[t_idx, z_idx, :, :] = u_grid
+            if "v" in output_vars:
+                v_data[t_idx, z_idx, :, :] = v_grid
+            if "w" in output_vars:
+                w_data[t_idx, z_idx, :, :] = w_grid
+    
+    # Create xarray dataset
+    coords = {
+        "time": list(findex_range),
+        "z": z_list,
+        "y": y_coords,
+        "x": x_coords,
+    }
+    
+    if "u" in output_vars:
+        data_vars["u"] = (("time", "z", "y", "x"), u_data)
+    if "v" in output_vars:
+        data_vars["v"] = (("time", "z", "y", "x"), v_data)
+    if "w" in output_vars:
+        data_vars["w"] = (("time", "z", "y", "x"), w_data)
+    
+    ds = xr.Dataset(data_vars=data_vars, coords=coords)
+    
+    # Add metadata
+    ds.attrs["description"] = "FLORIS flow field output"
+    ds.x.attrs["units"] = "m"
+    ds.y.attrs["units"] = "m"
+    ds.z.attrs["units"] = "m"
+    
+    if "u" in output_vars:
+        ds.u.attrs["units"] = "m/s"
+        ds.u.attrs["long_name"] = "x-component of velocity"
+    if "v" in output_vars:
+        ds.v.attrs["units"] = "m/s"
+        ds.v.attrs["long_name"] = "y-component of velocity"
+    if "w" in output_vars:
+        ds.w.attrs["units"] = "m/s"
+        ds.w.attrs["long_name"] = "z-component of velocity"
+    
+    # Write to NetCDF
+    ds.to_netcdf(f"{output_dir}/{flow_nc_filename}")
+    
+    return ds
\ No newline at end of file
diff --git a/wifa/main_api.py b/wifa/main_api.py
index 0ddef48..8ee8ed3 100644
--- a/wifa/main_api.py
+++ b/wifa/main_api.py
@@ -11,16 +11,20 @@
 from .foxes_api import run_foxes
 from .pywake_api import run_pywake
 from .wayve_api import run_wayve
+from .floris_api import run_floris
 
 sys.path.append(windIO.__path__[0])
 
 
 def run_api(yaml_input):
     # validate input
-    validate_yaml(yaml_input, windIO.__path__[0] + "/plant/wind_energy_system.yaml")
+    validate_yaml(yaml_input, "/root/code/libs/windIO_flow/windIO/schemas/plant/wind_energy_system")
 
     # get number of turbines
-    yaml_dat = load_yaml(yaml_input)
+    if isinstance(yaml_input, dict):
+        yaml_dat = yaml_input
+    else:
+        yaml_dat = load_yaml(yaml_input)
 
     model_name = yaml_dat["attributes"]["flow_model"]["name"]
 
@@ -30,6 +34,9 @@ def run_api(yaml_input):
     elif model_name.lower() == "foxes":
         foxes_aep = run_foxes(yaml_input)
 
+    elif model_name.lower() == "floris":
+        floris_aep = run_floris(yaml_input)
+
     elif model_name.lower() == "wayve":
         # Output directory
         # yaml_input_no_ext = os.path.splitext(yaml_input)[0]  # Remove the file extension

From 2aa1f9d3c8023a5e92891b0ca0b1b14c953526d4 Mon Sep 17 00:00:00 2001
From: lejeunemax <max-lej@hotmail.com>
Date: Tue, 16 Dec 2025 15:51:24 +0100
Subject: [PATCH 2/7] Add Flois coupling

---
 tests/test_floris.py | 132 +++++++++++++++++++++++++++++++++++++++++++
 wifa/floris_api.py   |  27 +++++++--
 2 files changed, 155 insertions(+), 4 deletions(-)
 create mode 100644 tests/test_floris.py

diff --git a/tests/test_floris.py b/tests/test_floris.py
new file mode 100644
index 0000000..b802535
--- /dev/null
+++ b/tests/test_floris.py
@@ -0,0 +1,132 @@
+import os
+from pathlib import Path
+
+from windIO import __path__ as wiop
+from windIO import validate as validate_yaml
+from windIO import load_yaml
+
+from wifa.floris_api import run_floris
+
+test_path = Path(os.path.dirname(__file__))
+windIO_path = Path(wiop[0])
+
+
+def _run_floris(wes_dir, expected_error=None):
+    """
+    Helper function to run FLORIS on all system.yaml files in a directory.
+    
+    Args:
+        wes_dir: Path to wind energy system directory containing system.yaml files
+        expected_error: If provided, the test expects this error type/message
+    """
+    assert wes_dir.is_dir(), f"{wes_dir} is not a directory"
+
+    for yaml_path in wes_dir.glob("system.yaml"):
+        print("\nRUNNING FLORIS ON", yaml_path, "\n")
+        
+        # Validate the YAML file against windIO schema
+        validate_yaml(yaml_path, Path("plant/wind_energy_system"))
+        
+        # Load the YAML file as a dictionary
+        yaml_input = load_yaml(yaml_path)
+        yaml_input['attributes']['flow_model']['name'] = 'floris'
+        
+        # Create output directory
+        output_dir_name = Path("output_test_floris")
+        output_dir_name.mkdir(parents=True, exist_ok=True)
+        
+        # Run FLORIS
+        result = run_floris(yaml_input)
+
+
+def test_floris_4wts():
+    """Test FLORIS with 4 turbines configuration"""
+    wes_dir = test_path / "../examples/cases/windio_4turbines/wind_energy_system/"
+    _run_floris(wes_dir)
+
+
+def test_floris_simple_wind_rose():
+    """Test FLORIS with simple wind rose configuration"""
+    wes_dir = test_path / "../examples/cases/simple_wind_rose/wind_energy_system/"
+    _run_floris(wes_dir)
+
+
+def test_floris_timeseries():
+    """Test FLORIS with timeseries and operating flag"""
+    wes_dir = test_path / "../examples/cases/timeseries_with_operating_flag/wind_energy_system/"
+    _run_floris(wes_dir)
+
+
+def test_floris_heterogeneous_wind_rose():
+    """Test FLORIS with heterogeneous wind rose map"""
+    wes_dir = test_path / "../examples/cases/heterogeneous_wind_rose_map/wind_energy_system/"
+    _run_floris(wes_dir)
+
+
+def test_floris_heterogeneous_at_turbines():
+    """Test FLORIS with heterogeneous wind rose at turbines"""
+    wes_dir = test_path / "../examples/cases/heterogeneous_wind_rose_at_turbines/wind_energy_system/"
+    _run_floris(wes_dir)
+
+
+def test_floris_multiple_turbines():
+    """Test FLORIS with multiple turbine types"""
+    wes_dir = test_path / "../examples/cases/windio_4turbines_multipleTurbines/wind_energy_system/"
+    _run_floris(wes_dir)
+
+
+if __name__ == "__main__":
+    """Run tests when executed directly"""
+    import sys
+    import traceback
+    
+    print("=" * 80)
+    print("Running FLORIS tests...")
+    print("=" * 80)
+    
+    tests = [
+        ("4 Turbines", test_floris_4wts),
+        ("Simple Wind Rose", test_floris_simple_wind_rose),
+        ("Timeseries", test_floris_timeseries),
+        ("Multiple Turbine Types", test_floris_multiple_turbines),
+    ]
+    
+    failed = []
+    passed = []
+    known_issues = []
+    
+    for test_name, test_func in tests:
+        print(f"\n{'='*80}")
+        print(f"Testing: {test_name}")
+        print(f"{'='*80}")
+        try:
+            test_func()
+            passed.append((test_name, ""))
+            print(f"✓ {test_name} PASSED")
+        except Exception as e:
+            failed.append((test_name, str(e)))
+            print(f"✗ {test_name} FAILED: {e}")
+            if "--verbose" in sys.argv or "-v" in sys.argv:
+                traceback.print_exc()
+    
+    # Summary
+    print(f"\n{'='*80}")
+    print("Test Summary")
+    print(f"{'='*80}")
+    print(f"Passed: {len(passed)}/{len(tests)}")
+    print(f"Known Issues: {len(known_issues)}/{len(tests)}")
+    print(f"Failed: {len(failed)}/{len(tests)}")
+    
+    if known_issues:
+        print("\nKnown issues (NotImplementedError):")
+        for name, error in known_issues:
+            print(f"  WARNING: {name}: {error}")
+    
+    if failed:
+        print("\nFailed tests:")
+        for name, error in failed:
+            print(f"  ✗ {name}: {error}")
+        sys.exit(1)
+    else:
+        print("\n✓ All tests passed!")
+        sys.exit(0)
diff --git a/wifa/floris_api.py b/wifa/floris_api.py
index 9a0ad93..7bca640 100644
--- a/wifa/floris_api.py
+++ b/wifa/floris_api.py
@@ -1,5 +1,5 @@
 from windIO import load_yaml
-from typing import Dict, Any, List
+from typing import Dict, List
 import numpy as np
 import windIO
 from pathlib import Path
@@ -9,6 +9,26 @@
     from floris import FlorisModel
 
 def run_floris(yaml_input):
+    """
+    Run FLORIS model based on windIO YAML input and extract specified outputs.
+    
+    Parameters
+    ----------
+    yaml_input: str, Path, or dict
+        Path to the windIO YAML input file or a dictionary representing the input data.
+        
+    Notes
+    -----
+    WindIO is not yet suppoeted by FLOIS main branch. 
+    Currently refer to https://github.com/lejeunemax/floris/tree/windIO.
+    
+    Some features may not be fully supported, these include:
+    - Subsetting of wind conditions (times_run, wind_speeds_run, directions_run)
+    - Only "power", "thrust_coefficient", "axial_induction", and "turbulence_intensity"
+      are currently supported as turbine output variables.
+    - Some advanced windIO properties may not be fully supported. These include:
+      blockage_model, rotor_averaging, and axial_induction_model.
+    """
 
     from floris import FlorisModel
     from floris.read_windio import TrackedDict
@@ -57,15 +77,14 @@ def run_floris(yaml_input):
                 )
         
         output_dir = out_specs.get("output_folder", ".")
+        Path(output_dir).mkdir(parents=True, exist_ok=True)
 
         if ("turbine_outputs" in out_specs):
             _read_turbines(fmodel, output_dir, out_specs, slice_selection)
 
         if ("flow_field" in out_specs):
             _read_fields(fmodel, output_dir, out_specs, slice_selection)
-        
-        pass
-        
+
 def _read_turbines(fmodel: "FlorisModel", output_dir, out_specs, slice_selection):
     """
     Extract turbine output data from FLORIS model and write to NetCDF file.

From d498b974ff9ac4faf03064d4a7a787528170b030 Mon Sep 17 00:00:00 2001
From: lejeunemax <max-lej@hotmail.com>
Date: Thu, 18 Dec 2025 11:31:58 +0100
Subject: [PATCH 3/7] removing local path

---
 wifa/main_api.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/wifa/main_api.py b/wifa/main_api.py
index 8ee8ed3..4a18637 100644
--- a/wifa/main_api.py
+++ b/wifa/main_api.py
@@ -18,8 +18,8 @@
 
 def run_api(yaml_input):
     # validate input
-    validate_yaml(yaml_input, "/root/code/libs/windIO_flow/windIO/schemas/plant/wind_energy_system")
-
+    validate_yaml(yaml_input, windIO.__path__[0] + "/plant/wind_energy_system.yaml")
+    
     # get number of turbines
     if isinstance(yaml_input, dict):
         yaml_dat = yaml_input

From 9a044ba763717e7f3e59ae7714f36771e0b48287 Mon Sep 17 00:00:00 2001
From: lejeunemax <max-lej@hotmail.com>
Date: Wed, 7 Jan 2026 13:20:54 +0100
Subject: [PATCH 4/7] Make DTU10MW def consistent accros files

---
 .../plant_energy_turbine/DTU_10MW_turbine.yaml            | 8 ++++----
 .../plant_energy_turbine/DTU_10MW_turbine.yaml            | 8 ++++----
 2 files changed, 8 insertions(+), 8 deletions(-)

diff --git a/examples/cases/timeseries_with_operating_flag/plant_energy_turbine/DTU_10MW_turbine.yaml b/examples/cases/timeseries_with_operating_flag/plant_energy_turbine/DTU_10MW_turbine.yaml
index 0320701..8ac1d6d 100644
--- a/examples/cases/timeseries_with_operating_flag/plant_energy_turbine/DTU_10MW_turbine.yaml
+++ b/examples/cases/timeseries_with_operating_flag/plant_energy_turbine/DTU_10MW_turbine.yaml
@@ -1,10 +1,10 @@
 name: DTU 10MW Offshore Reference Turbine
 performance:
   power_curve:
-    power_values: [263388., 751154., 1440738., 2355734., 3506858., 4993092., 6849310., 9116402., 10000754., 10009590., 10000942., 10042678., 10003480., 10001600., 10001506., 10013632., 10007428., 10005360., 10002728., 10001130., 10004984., 9997558.]
-    power_wind_speeds: [4.,5.,6.,7.,8.,9.,10.,11.,12.,13.,14.,15.,16.,17.,18.,19.,20.,21.,22.,23.,24.,25.]
+    power_values: [0, 263388., 751154., 1440738., 2355734., 3506858., 4993092., 6849310., 9116402., 10000754., 10009590., 10000942., 10042678., 10003480., 10001600., 10001506., 10013632., 10007428., 10005360., 10002728., 10001130., 10004984., 9997558.]
+    power_wind_speeds: [3, 4.,5.,6.,7.,8.,9.,10.,11.,12.,13.,14.,15.,16.,17.,18.,19.,20.,21.,22.,23.,24.,25.]
   Ct_curve:
-    Ct_values: [0.923,0.919,0.904,0.858,0.814,0.814,0.814,0.814,0.577,0.419,0.323,0.259,0.211,0.175,0.148,0.126,0.109,0.095,0.084,0.074,0.066,0.059]
-    Ct_wind_speeds: [4.,5.,6.,7.,8.,9.,10.,11.,12.,13.,14.,15.,16.,17.,18.,19.,20.,21.,22.,23.,24.,25.]
+    Ct_values: [0.923, 0.923,0.919,0.904,0.858,0.814,0.814,0.814,0.814,0.577,0.419,0.323,0.259,0.211,0.175,0.148,0.126,0.109,0.095,0.084,0.074,0.066,0.059]
+    Ct_wind_speeds: [3, 4.,5.,6.,7.,8.,9.,10.,11.,12.,13.,14.,15.,16.,17.,18.,19.,20.,21.,22.,23.,24.,25.]
 hub_height: 119.0
 rotor_diameter: 178.3
diff --git a/examples/cases/windio_4turbines/plant_energy_turbine/DTU_10MW_turbine.yaml b/examples/cases/windio_4turbines/plant_energy_turbine/DTU_10MW_turbine.yaml
index 0320701..8ac1d6d 100644
--- a/examples/cases/windio_4turbines/plant_energy_turbine/DTU_10MW_turbine.yaml
+++ b/examples/cases/windio_4turbines/plant_energy_turbine/DTU_10MW_turbine.yaml
@@ -1,10 +1,10 @@
 name: DTU 10MW Offshore Reference Turbine
 performance:
   power_curve:
-    power_values: [263388., 751154., 1440738., 2355734., 3506858., 4993092., 6849310., 9116402., 10000754., 10009590., 10000942., 10042678., 10003480., 10001600., 10001506., 10013632., 10007428., 10005360., 10002728., 10001130., 10004984., 9997558.]
-    power_wind_speeds: [4.,5.,6.,7.,8.,9.,10.,11.,12.,13.,14.,15.,16.,17.,18.,19.,20.,21.,22.,23.,24.,25.]
+    power_values: [0, 263388., 751154., 1440738., 2355734., 3506858., 4993092., 6849310., 9116402., 10000754., 10009590., 10000942., 10042678., 10003480., 10001600., 10001506., 10013632., 10007428., 10005360., 10002728., 10001130., 10004984., 9997558.]
+    power_wind_speeds: [3, 4.,5.,6.,7.,8.,9.,10.,11.,12.,13.,14.,15.,16.,17.,18.,19.,20.,21.,22.,23.,24.,25.]
   Ct_curve:
-    Ct_values: [0.923,0.919,0.904,0.858,0.814,0.814,0.814,0.814,0.577,0.419,0.323,0.259,0.211,0.175,0.148,0.126,0.109,0.095,0.084,0.074,0.066,0.059]
-    Ct_wind_speeds: [4.,5.,6.,7.,8.,9.,10.,11.,12.,13.,14.,15.,16.,17.,18.,19.,20.,21.,22.,23.,24.,25.]
+    Ct_values: [0.923, 0.923,0.919,0.904,0.858,0.814,0.814,0.814,0.814,0.577,0.419,0.323,0.259,0.211,0.175,0.148,0.126,0.109,0.095,0.084,0.074,0.066,0.059]
+    Ct_wind_speeds: [3, 4.,5.,6.,7.,8.,9.,10.,11.,12.,13.,14.,15.,16.,17.,18.,19.,20.,21.,22.,23.,24.,25.]
 hub_height: 119.0
 rotor_diameter: 178.3

From 965092369b0d4a7de7bcb2812daf19bb5f4728f6 Mon Sep 17 00:00:00 2001
From: lejeunemax <max-lej@hotmail.com>
Date: Wed, 7 Jan 2026 13:22:13 +0100
Subject: [PATCH 5/7] printing output files

---
 wifa/floris_api.py | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/wifa/floris_api.py b/wifa/floris_api.py
index 7bca640..8cbf0bb 100644
--- a/wifa/floris_api.py
+++ b/wifa/floris_api.py
@@ -247,6 +247,7 @@ def _read_turbines(fmodel: "FlorisModel", output_dir, out_specs, slice_selection
     
     # Write to NetCDF
     ds.to_netcdf(f"{output_dir}/{turbine_nc_filename}")
+    print(f"Turbine output data written to {output_dir}/{turbine_nc_filename}")
     
     return ds
     
@@ -424,5 +425,6 @@ def _read_fields(fmodel: "FlorisModel", output_dir, out_specs, slice_selection):
     
     # Write to NetCDF
     ds.to_netcdf(f"{output_dir}/{flow_nc_filename}")
+    print(f"Flow field data written to {output_dir}/{flow_nc_filename}")
     
     return ds
\ No newline at end of file

From 98f521c427a7b6751fc7fea86a7d73f13f1691f5 Mon Sep 17 00:00:00 2001
From: lejeunemax <max-lej@hotmail.com>
Date: Wed, 7 Jan 2026 13:22:55 +0100
Subject: [PATCH 6/7] Improved testing

---
 tests/test_floris.py | 408 ++++++++++++++++++++++++++++++++-----------
 1 file changed, 308 insertions(+), 100 deletions(-)

diff --git a/tests/test_floris.py b/tests/test_floris.py
index b802535..8a5f847 100644
--- a/tests/test_floris.py
+++ b/tests/test_floris.py
@@ -1,132 +1,340 @@
+"""
+FLORIS Testing Suite for WIFA Integration
+
+Validates WIFA-generated results against direct FLORIS calculations.
+Tests single/multiple turbine farms, wind roses, and operating flags.
+
+Run with: pytest tests/test_floris.py -v
+"""
+
 import os
+import numpy as np
+import shutil
 from pathlib import Path
+import pytest
 
 from windIO import __path__ as wiop
 from windIO import validate as validate_yaml
 from windIO import load_yaml
 
+import floris
+from floris.turbine_library import build_cosine_loss_turbine_dict
 from wifa.floris_api import run_floris
+import xarray as xr
+
+# ============================================================================ #
+# CONFIGURATION
+# ============================================================================ #
 
+CWD = Path(__file__).parent.resolve()
+save_dir = CWD / ".output_test_floris"
 test_path = Path(os.path.dirname(__file__))
 windIO_path = Path(wiop[0])
 
+# ============================================================================ #
+# PYTEST FIXTURES
+# ============================================================================ #
+
+@pytest.fixture(scope="session")
+def floris_config():
+    """Default FLORIS configuration for all tests."""
+    return get_default_floris_dict()
+
+@pytest.fixture
+def cleanup_temp_files():
+    """Clean up temporary files after each test."""
+    yield
+    # Cleanup after test
+    if save_dir.exists():
+        shutil.rmtree(save_dir)
+
+# ============================================================================ #
+# UTILITY FUNCTIONS
+# ============================================================================ #
+
+def recursive_dict_update(original, updates):
+    """Recursively update dictionary with nested structure preservation."""
+    for key, value in updates.items():
+        if isinstance(value, dict) and key in original:
+            recursive_dict_update(original[key], value)
+        else:
+            original[key] = value
+    return original
+
+def get_default_floris_dict():
+    """Load default FLORIS configuration."""
+    floris_dict_default_yaml = Path(floris.__path__[0]) / 'default_inputs.yaml'
+    return load_yaml(floris_dict_default_yaml)
+
+def create_output_dir():
+    """Create test output directory."""
+    output_dir_name = Path("output_test_floris")
+    output_dir_name.mkdir(parents=True, exist_ok=True)
+    return output_dir_name
+
+def _run_floris(yaml_path):
+    """Run FLORIS via WIFA, load outputs, and clean up."""
+    print(f"\nRUNNING FLORIS ON {yaml_path}\n")
+    
+    # Validate and load configuration
+    validate_yaml(yaml_path, Path("plant/wind_energy_system"))
+    yaml_input = load_yaml(yaml_path)
+    yaml_input['attributes']['model_outputs_specification']['output_folder'] = save_dir
+    yaml_input['attributes']['flow_model']['name'] = 'floris'
+    
+    # Run simulation
+    save_dir.mkdir(parents=True, exist_ok=True)
+    run_floris(yaml_input)
+    
+    # Load NetCDF outputs
+    output_data = {}
+    if save_dir.exists():
+        for file_path in save_dir.iterdir():
+            if file_path.is_file() and file_path.suffix == '.nc':
+                output_data[file_path.name] = xr.load_dataset(file_path)
+    
+    # Clean up
+    if save_dir.exists():
+        shutil.rmtree(save_dir)
+    
+    return output_data
+
+# ============================================================================ #
+# TEST ENVIRONMENT SETUP
+# ============================================================================ #
 
-def _run_floris(wes_dir, expected_error=None):
-    """
-    Helper function to run FLORIS on all system.yaml files in a directory.
+def create_turbine(turbine_name="dtu_10mw"):
+    """Create FLORIS turbine dictionary from YAML. Supports DTU and IEA turbines."""
+    # Determine YAML path
+    if turbine_name == "dtu_10mw":
+        turbine_yaml_path = test_path / "../examples/cases/windio_4turbines/plant_energy_turbine/DTU_10MW_turbine.yaml"
+    elif turbine_name in ["IEA_10MW", "IEA_15MW", "IEA_22MW"]:
+        turbine_yaml_path = test_path / f"../examples/cases/windio_4turbines_multipleTurbines/plant_energy_turbine/{turbine_name}_turbine.yaml"
+    else:
+        raise ValueError(f"Unknown turbine: {turbine_name}. Options: 'dtu_10mw', 'IEA_10MW', 'IEA_15MW', 'IEA_22MW'")
+    
+    turbine_data = load_yaml(turbine_yaml_path)
     
-    Args:
-        wes_dir: Path to wind energy system directory containing system.yaml files
-        expected_error: If provided, the test expects this error type/message
-    """
-    assert wes_dir.is_dir(), f"{wes_dir} is not a directory"
-
-    for yaml_path in wes_dir.glob("system.yaml"):
-        print("\nRUNNING FLORIS ON", yaml_path, "\n")
+    # Handle different data formats
+    if 'power_curve' in turbine_data['performance']:
+        # Direct power curve (DTU)
+        wind_speeds = turbine_data['performance']['power_curve']['power_wind_speeds']
+        power_values = turbine_data['performance']['power_curve']['power_values']
+        ct_values = turbine_data['performance']['Ct_curve']['Ct_values']
         
-        # Validate the YAML file against windIO schema
-        validate_yaml(yaml_path, Path("plant/wind_energy_system"))
+        # Auto-convert W to kW
+        power_values_kw = [p / 1000.0 for p in power_values] if max(power_values) > 100000 else power_values
+            
+    elif 'Cp_curve' in turbine_data['performance']:
+        # Calculate from Cp curve (IEA)
+        wind_speeds = turbine_data['performance']['Cp_curve']['Cp_wind_speeds']
+        cp_values = turbine_data['performance']['Cp_curve']['Cp_values']
+        ct_values = turbine_data['performance']['Ct_curve']['Ct_values']
         
-        # Load the YAML file as a dictionary
-        yaml_input = load_yaml(yaml_path)
-        yaml_input['attributes']['flow_model']['name'] = 'floris'
+        # P = 0.5 * ρ * A * V³ * Cp
+        air_density = 1.225
+        rotor_area = np.pi * (turbine_data['rotor_diameter'] / 2) ** 2
+        power_values_kw = [0.5 * air_density * rotor_area * (ws ** 3) * cp / 1000 
+                          for ws, cp in zip(wind_speeds, cp_values)]
+    else:
+        raise ValueError("YAML must contain 'power_curve' or 'Cp_curve'")
+    
+    power_thrust_table = {
+        'power': power_values_kw,
+        'thrust_coefficient': ct_values,
+        'wind_speed': wind_speeds,
+    }
+
+    return build_cosine_loss_turbine_dict(
+        power_thrust_table, turbine_name.lower(), generator_efficiency=1.0,
+        hub_height=turbine_data["hub_height"], rotor_diameter=turbine_data["rotor_diameter"], TSR=7.0
+    )
+
+def create_four_turbine_farm(layout_type: str):
+    """Create 4-turbine farm. 'same' = all DTU, 'multiple' = mixed IEA types."""
+    if layout_type == "same":
+        turbine_dict = create_turbine("dtu_10mw")
+        turbine_dicts = [turbine_dict] * 4
+    elif layout_type == "multiple":
+        turbine_types = ["IEA_10MW", "IEA_15MW", "IEA_15MW", "IEA_22MW"]
+        turbine_dicts = [create_turbine(t_type) for t_type in turbine_types]
+    else:
+        raise ValueError(f"Unknown layout_type: {layout_type}. Use 'same' or 'multiple'.")
+    
+    return {
+        "layout_x": [0, 1248.1, 2496.2, 3744.3],
+        "layout_y": [0.0, 0.0, 0.0, 0.0],
+        "turbine_type": turbine_dicts,
+    }
+
+def create_timeseries():
+    """Create 3-step wind time series for testing."""
+    return {
+        'wind_speeds': [10.091022491455078, 8.797999382019043, 10.307791709899902],
+        'wind_directions': [271.8246154785156, 268.6852111816406, 271.0501403808594],
+        'turbulence_intensities': [2.6189446449279785, 1.6507437229156494, 3.7473068237304688],
+    }
+
+def create_timeseries_with_operating_flag():
+    """Create time series with turbine availability flags."""
+    return {
+        'wind_directions': [271.8246, 266.2015, 268.6852, 273.6164, 263.4558, 271.0501],
+        'wind_speeds': [10.09102, 10.23302, 8.797999, 9.662098, 9.78371, 10.30779],
+        'turbulence_intensities': [2.618945, 1.537051, 1.650744, 3.149817, 0.5826571, 3.747307],
+        'disable_turbines': np.logical_not(np.array([
+            [0, 1, 1, 1], [0, 1, 1, 1], [0, 1, 1, 1],  # T0 offline
+            [0, 1, 1, 1], [1, 1, 1, 1], [1, 1, 1, 1]   # All online
+        ]))
+    }
+
+def create_simple_wind_rose_energy_resource():
+    """Create wind rose with Weibull distributions for 12 sectors."""
+    # Sector probabilities and Weibull parameters
+    sector_probability = [0.0359715203597152, 0.0394868203948682, 0.0516739505167395, 0.0700015407000154,
+                         0.0836454708364547, 0.064348500643485, 0.0864319408643194, 0.117705101177051,
+                         0.151575701515757, 0.147379201473792, 0.100120501001205, 0.0516597505165975, 0.0359715203597152]
+    
+    weibull_a = [9.176929, 9.782334, 9.531809, 9.909545, 10.04269, 9.593921, 9.584007, 10.51499, 
+                11.39895, 11.68746, 11.63732, 10.08803, 9.176929]
+    
+    weibull_k = [2.392578, 2.447266, 2.412109, 2.591797, 2.755859, 2.595703, 2.583984, 2.548828,
+                2.470703, 2.607422, 2.626953, 2.326172, 2.392578]
+
+    wind_direction = [0, 30, 60, 90, 120, 150, 180, 210, 240, 270, 300, 330, 360]
+    wind_speed = [3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25]
+    
+    # Normalize probabilities and calculate frequency table
+    sector_probability = np.array(sector_probability) / np.sum(sector_probability)
+    wind_speeds = np.array(wind_speed)
+    freq_table = np.zeros((len(wind_direction), len(wind_speed)))
+    
+    for i_dir in range(len(wind_direction)):
+        A, k = weibull_a[i_dir], weibull_k[i_dir]
+        wind_speed_edges = np.arange(wind_speeds[0] - 0.5, wind_speeds[-1] + 1.5, 1.0)
         
-        # Create output directory
-        output_dir_name = Path("output_test_floris")
-        output_dir_name.mkdir(parents=True, exist_ok=True)
+        # Weibull CDF: F(x) = 1 - exp(-(x/A)^k)
+        cdf_edges = 1.0 - np.exp(-((wind_speed_edges / A) ** k))
+        cdf_edges[wind_speed_edges < 0] = 0.0
         
-        # Run FLORIS
-        result = run_floris(yaml_input)
+        freq = np.diff(cdf_edges)
+        freq_table[i_dir, :] = (freq / freq.sum()) * sector_probability[i_dir]
 
+    wr = floris.WindRose(wind_directions=np.array(wind_direction), wind_speeds=wind_speeds, 
+                        ti_table=0.1, freq_table=freq_table)
+    return {"wind_data": wr}
 
-def test_floris_4wts():
-    """Test FLORIS with 4 turbines configuration"""
-    wes_dir = test_path / "../examples/cases/windio_4turbines/wind_energy_system/"
-    _run_floris(wes_dir)
+def create_wake_model_config():
+    """Standard wake model configuration for testing."""
+    return {
+        'enable_active_wake_mixing': False,
+        'enable_secondary_steering': False,
+        'enable_transverse_velocities': False,
+        'enable_yaw_added_recovery': False,
+        'model_strings': {
+            'combination_model': 'fls',
+            'deflection_model': 'none',
+            'velocity_model': 'gauss',
+            'turbulence_model': 'crespo_hernandez',
+        },
+        'wake_velocity_parameters': {'gauss': {'ka': 0.0, 'kb': 0.04}},
+    }
 
+# ============================================================================ #
+# TEST FUNCTIONS
+# ============================================================================ #
 
-def test_floris_simple_wind_rose():
-    """Test FLORIS with simple wind rose configuration"""
-    wes_dir = test_path / "../examples/cases/simple_wind_rose/wind_energy_system/"
-    _run_floris(wes_dir)
+def test_floris_4wts(floris_config):
+    """Test homogeneous 4-turbine farm. Max diff < 1e-6."""
+    wes_dir = test_path / "../examples/cases/windio_4turbines/wind_energy_system/system.yaml"
+    output_data = _run_floris(wes_dir)
+    powers_wifa = output_data['turbine_data.nc']['power'].values
 
+    floris_dict = floris_config.copy()
+    recursive_dict_update(floris_dict['farm'], create_four_turbine_farm("same"))
+    recursive_dict_update(floris_dict['flow_field'], create_timeseries())
+    recursive_dict_update(floris_dict['wake'], create_wake_model_config())
 
-def test_floris_timeseries():
-    """Test FLORIS with timeseries and operating flag"""
-    wes_dir = test_path / "../examples/cases/timeseries_with_operating_flag/wind_energy_system/"
-    _run_floris(wes_dir)
+    fmodel = floris.FlorisModel(floris_dict)
+    fmodel.run()
+    powers_floris = fmodel.get_turbine_powers()
 
+    # Use pytest.approx for better numerical comparison
+    assert powers_wifa == pytest.approx(powers_floris, rel=1e-6), f"Power outputs don't match within 1e-6 tolerance"
 
-def test_floris_heterogeneous_wind_rose():
-    """Test FLORIS with heterogeneous wind rose map"""
-    wes_dir = test_path / "../examples/cases/heterogeneous_wind_rose_map/wind_energy_system/"
-    _run_floris(wes_dir)
+def test_floris_simple_wind_rose(floris_config):
+    """Test wind rose analysis. Max diff < 1e-2."""
+    wes_dir = test_path / "../examples/cases/simple_wind_rose/wind_energy_system/system.yaml"
+    output_data = _run_floris(wes_dir)
+    powers_wifa = output_data['PowerTable.nc']['power'].values
 
+    floris_dict = floris_config.copy()
+    recursive_dict_update(floris_dict['farm'], create_four_turbine_farm("same"))
+    recursive_dict_update(floris_dict['wake'], create_wake_model_config())
 
-def test_floris_heterogeneous_at_turbines():
-    """Test FLORIS with heterogeneous wind rose at turbines"""
-    wes_dir = test_path / "../examples/cases/heterogeneous_wind_rose_at_turbines/wind_energy_system/"
-    _run_floris(wes_dir)
+    fmodel = floris.FlorisModel(floris_dict)
+    fmodel.set(**create_simple_wind_rose_energy_resource())
+    fmodel.run()
+    powers_floris = fmodel.get_turbine_powers()
 
+    assert powers_wifa == pytest.approx(powers_floris, rel=1e-2), f"Wind rose power outputs don't match within 1e-2 tolerance"
 
-def test_floris_multiple_turbines():
-    """Test FLORIS with multiple turbine types"""
-    wes_dir = test_path / "../examples/cases/windio_4turbines_multipleTurbines/wind_energy_system/"
-    _run_floris(wes_dir)
+def test_floris_timeseries_with_operating_flag(floris_config):
+    """Test turbine availability flags. Max diff < 1e-2."""
+    wes_dir = test_path / "../examples/cases/timeseries_with_operating_flag/wind_energy_system/system.yaml"
+    output_data = _run_floris(wes_dir)
+    powers_wifa = output_data['turbine_data.nc']['power'].values
 
-
-if __name__ == "__main__":
-    """Run tests when executed directly"""
-    import sys
-    import traceback
-    
-    print("=" * 80)
-    print("Running FLORIS tests...")
-    print("=" * 80)
-    
-    tests = [
-        ("4 Turbines", test_floris_4wts),
-        ("Simple Wind Rose", test_floris_simple_wind_rose),
-        ("Timeseries", test_floris_timeseries),
-        ("Multiple Turbine Types", test_floris_multiple_turbines),
-    ]
-    
-    failed = []
-    passed = []
-    known_issues = []
-    
-    for test_name, test_func in tests:
-        print(f"\n{'='*80}")
-        print(f"Testing: {test_name}")
-        print(f"{'='*80}")
-        try:
-            test_func()
-            passed.append((test_name, ""))
-            print(f"✓ {test_name} PASSED")
-        except Exception as e:
-            failed.append((test_name, str(e)))
-            print(f"✗ {test_name} FAILED: {e}")
-            if "--verbose" in sys.argv or "-v" in sys.argv:
-                traceback.print_exc()
-    
-    # Summary
-    print(f"\n{'='*80}")
-    print("Test Summary")
-    print(f"{'='*80}")
-    print(f"Passed: {len(passed)}/{len(tests)}")
-    print(f"Known Issues: {len(known_issues)}/{len(tests)}")
-    print(f"Failed: {len(failed)}/{len(tests)}")
+    floris_dict = floris_config.copy()
+    recursive_dict_update(floris_dict['farm'], create_four_turbine_farm("same"))
+    recursive_dict_update(floris_dict['wake'], create_wake_model_config())
     
-    if known_issues:
-        print("\nKnown issues (NotImplementedError):")
-        for name, error in known_issues:
-            print(f"  WARNING: {name}: {error}")
-    
-    if failed:
-        print("\nFailed tests:")
-        for name, error in failed:
-            print(f"  ✗ {name}: {error}")
-        sys.exit(1)
-    else:
-        print("\n✓ All tests passed!")
-        sys.exit(0)
+    fmodel = floris.FlorisModel(floris_dict)
+    fmodel.set_operation_model("mixed")
+    fmodel.set(**create_timeseries_with_operating_flag())
+    fmodel.run()
+    powers_floris = fmodel.get_turbine_powers()
+
+    assert powers_wifa == pytest.approx(powers_floris, rel=1e-2), f"Operating flag power outputs don't match within 1e-2 tolerance"
+
+def test_floris_multiple_turbines(floris_config):
+    """Test mixed turbine types. Max diff < 0.1. Marked as slow test."""
+    wes_dir = test_path / "../examples/cases/windio_4turbines_multipleTurbines/wind_energy_system/system.yaml"
+    output_data = _run_floris(wes_dir)
+    powers_wifa = output_data['turbine_data.nc']['power'].values
+
+    floris_dict = floris_config.copy()
+    recursive_dict_update(floris_dict['farm'], create_four_turbine_farm("multiple"))
+    recursive_dict_update(floris_dict['flow_field'], create_timeseries())
+    recursive_dict_update(floris_dict['wake'], create_wake_model_config())
+
+    # Set reference height to mean hub height
+    mean_hub_height = np.mean([t["hub_height"] for t in floris_dict["farm"]["turbine_type"]])
+    floris_dict['flow_field']['reference_wind_height'] = mean_hub_height
+
+    fmodel = floris.FlorisModel(floris_dict)
+    fmodel.run()
+    powers_floris = fmodel.get_turbine_powers()
+
+    max_diff = np.max(np.abs(powers_wifa - powers_floris)) / np.max(np.abs(powers_wifa))
+    print(f"Max relative difference: {max_diff:.6f}")
+    assert powers_wifa == pytest.approx(powers_floris, rel=0.1), f"Mixed turbine power outputs don't match within 0.1 tolerance"
+
+
+# ============================================================================ #
+# PYTEST FIXTURES
+# ============================================================================ #
+
+@pytest.fixture(scope="session")  
+def test_output_dir():
+    """Shared test output directory fixture with automatic cleanup."""
+    output_dir = Path("output_test_floris")
+    output_dir.mkdir(parents=True, exist_ok=True)
+    yield output_dir
+    if output_dir.exists():
+        shutil.rmtree(output_dir)
+
+@pytest.fixture
+def floris_config():
+    """Default FLORIS configuration for all tests."""
+    return get_default_floris_dict()
+

From 679152314982772d06fa79c08d482261af4b6efb Mon Sep 17 00:00:00 2001
From: lejeunemax <max-lej@hotmail.com>
Date: Mon, 19 Jan 2026 13:59:11 +0100
Subject: [PATCH 7/7] precommit check

---
 tests/test_floris.py | 361 ++++++++++++++++++++++++++++++-------------
 wifa/floris_api.py   | 261 +++++++++++++++++--------------
 2 files changed, 397 insertions(+), 225 deletions(-)

diff --git a/tests/test_floris.py b/tests/test_floris.py
index 8a5f847..5bed6e5 100644
--- a/tests/test_floris.py
+++ b/tests/test_floris.py
@@ -8,19 +8,19 @@
 """
 
 import os
-import numpy as np
 import shutil
 from pathlib import Path
-import pytest
 
+import floris
+import numpy as np
+import pytest
+import xarray as xr
+from floris.turbine_library import build_cosine_loss_turbine_dict
 from windIO import __path__ as wiop
-from windIO import validate as validate_yaml
 from windIO import load_yaml
+from windIO import validate as validate_yaml
 
-import floris
-from floris.turbine_library import build_cosine_loss_turbine_dict
 from wifa.floris_api import run_floris
-import xarray as xr
 
 # ============================================================================ #
 # CONFIGURATION
@@ -35,11 +35,13 @@
 # PYTEST FIXTURES
 # ============================================================================ #
 
+
 @pytest.fixture(scope="session")
 def floris_config():
     """Default FLORIS configuration for all tests."""
     return get_default_floris_dict()
 
+
 @pytest.fixture
 def cleanup_temp_files():
     """Clean up temporary files after each test."""
@@ -48,10 +50,12 @@ def cleanup_temp_files():
     if save_dir.exists():
         shutil.rmtree(save_dir)
 
+
 # ============================================================================ #
 # UTILITY FUNCTIONS
 # ============================================================================ #
 
+
 def recursive_dict_update(original, updates):
     """Recursively update dictionary with nested structure preservation."""
     for key, value in updates.items():
@@ -61,95 +65,119 @@ def recursive_dict_update(original, updates):
             original[key] = value
     return original
 
+
 def get_default_floris_dict():
     """Load default FLORIS configuration."""
-    floris_dict_default_yaml = Path(floris.__path__[0]) / 'default_inputs.yaml'
+    floris_dict_default_yaml = Path(floris.__path__[0]) / "default_inputs.yaml"
     return load_yaml(floris_dict_default_yaml)
 
+
 def create_output_dir():
     """Create test output directory."""
     output_dir_name = Path("output_test_floris")
     output_dir_name.mkdir(parents=True, exist_ok=True)
     return output_dir_name
 
+
 def _run_floris(yaml_path):
     """Run FLORIS via WIFA, load outputs, and clean up."""
     print(f"\nRUNNING FLORIS ON {yaml_path}\n")
-    
+
     # Validate and load configuration
     validate_yaml(yaml_path, Path("plant/wind_energy_system"))
     yaml_input = load_yaml(yaml_path)
-    yaml_input['attributes']['model_outputs_specification']['output_folder'] = save_dir
-    yaml_input['attributes']['flow_model']['name'] = 'floris'
-    
+    yaml_input["attributes"]["model_outputs_specification"]["output_folder"] = save_dir
+    yaml_input["attributes"]["flow_model"]["name"] = "floris"
+
     # Run simulation
     save_dir.mkdir(parents=True, exist_ok=True)
     run_floris(yaml_input)
-    
+
     # Load NetCDF outputs
     output_data = {}
     if save_dir.exists():
         for file_path in save_dir.iterdir():
-            if file_path.is_file() and file_path.suffix == '.nc':
+            if file_path.is_file() and file_path.suffix == ".nc":
                 output_data[file_path.name] = xr.load_dataset(file_path)
-    
+
     # Clean up
     if save_dir.exists():
         shutil.rmtree(save_dir)
-    
+
     return output_data
 
+
 # ============================================================================ #
 # TEST ENVIRONMENT SETUP
 # ============================================================================ #
 
+
 def create_turbine(turbine_name="dtu_10mw"):
     """Create FLORIS turbine dictionary from YAML. Supports DTU and IEA turbines."""
     # Determine YAML path
     if turbine_name == "dtu_10mw":
-        turbine_yaml_path = test_path / "../examples/cases/windio_4turbines/plant_energy_turbine/DTU_10MW_turbine.yaml"
+        turbine_yaml_path = (
+            test_path
+            / "../examples/cases/windio_4turbines/plant_energy_turbine/DTU_10MW_turbine.yaml"
+        )
     elif turbine_name in ["IEA_10MW", "IEA_15MW", "IEA_22MW"]:
-        turbine_yaml_path = test_path / f"../examples/cases/windio_4turbines_multipleTurbines/plant_energy_turbine/{turbine_name}_turbine.yaml"
+        turbine_yaml_path = (
+            test_path
+            / f"../examples/cases/windio_4turbines_multipleTurbines/plant_energy_turbine/{turbine_name}_turbine.yaml"
+        )
     else:
-        raise ValueError(f"Unknown turbine: {turbine_name}. Options: 'dtu_10mw', 'IEA_10MW', 'IEA_15MW', 'IEA_22MW'")
-    
+        raise ValueError(
+            f"Unknown turbine: {turbine_name}. Options: 'dtu_10mw', 'IEA_10MW', 'IEA_15MW', 'IEA_22MW'"
+        )
+
     turbine_data = load_yaml(turbine_yaml_path)
-    
+
     # Handle different data formats
-    if 'power_curve' in turbine_data['performance']:
+    if "power_curve" in turbine_data["performance"]:
         # Direct power curve (DTU)
-        wind_speeds = turbine_data['performance']['power_curve']['power_wind_speeds']
-        power_values = turbine_data['performance']['power_curve']['power_values']
-        ct_values = turbine_data['performance']['Ct_curve']['Ct_values']
-        
+        wind_speeds = turbine_data["performance"]["power_curve"]["power_wind_speeds"]
+        power_values = turbine_data["performance"]["power_curve"]["power_values"]
+        ct_values = turbine_data["performance"]["Ct_curve"]["Ct_values"]
+
         # Auto-convert W to kW
-        power_values_kw = [p / 1000.0 for p in power_values] if max(power_values) > 100000 else power_values
-            
-    elif 'Cp_curve' in turbine_data['performance']:
+        power_values_kw = (
+            [p / 1000.0 for p in power_values]
+            if max(power_values) > 100000
+            else power_values
+        )
+
+    elif "Cp_curve" in turbine_data["performance"]:
         # Calculate from Cp curve (IEA)
-        wind_speeds = turbine_data['performance']['Cp_curve']['Cp_wind_speeds']
-        cp_values = turbine_data['performance']['Cp_curve']['Cp_values']
-        ct_values = turbine_data['performance']['Ct_curve']['Ct_values']
-        
+        wind_speeds = turbine_data["performance"]["Cp_curve"]["Cp_wind_speeds"]
+        cp_values = turbine_data["performance"]["Cp_curve"]["Cp_values"]
+        ct_values = turbine_data["performance"]["Ct_curve"]["Ct_values"]
+
         # P = 0.5 * ρ * A * V³ * Cp
         air_density = 1.225
-        rotor_area = np.pi * (turbine_data['rotor_diameter'] / 2) ** 2
-        power_values_kw = [0.5 * air_density * rotor_area * (ws ** 3) * cp / 1000 
-                          for ws, cp in zip(wind_speeds, cp_values)]
+        rotor_area = np.pi * (turbine_data["rotor_diameter"] / 2) ** 2
+        power_values_kw = [
+            0.5 * air_density * rotor_area * (ws**3) * cp / 1000
+            for ws, cp in zip(wind_speeds, cp_values)
+        ]
     else:
         raise ValueError("YAML must contain 'power_curve' or 'Cp_curve'")
-    
+
     power_thrust_table = {
-        'power': power_values_kw,
-        'thrust_coefficient': ct_values,
-        'wind_speed': wind_speeds,
+        "power": power_values_kw,
+        "thrust_coefficient": ct_values,
+        "wind_speed": wind_speeds,
     }
 
     return build_cosine_loss_turbine_dict(
-        power_thrust_table, turbine_name.lower(), generator_efficiency=1.0,
-        hub_height=turbine_data["hub_height"], rotor_diameter=turbine_data["rotor_diameter"], TSR=7.0
+        power_thrust_table,
+        turbine_name.lower(),
+        generator_efficiency=1.0,
+        hub_height=turbine_data["hub_height"],
+        rotor_diameter=turbine_data["rotor_diameter"],
+        TSR=7.0,
     )
 
+
 def create_four_turbine_farm(layout_type: str):
     """Create 4-turbine farm. 'same' = all DTU, 'multiple' = mixed IEA types."""
     if layout_type == "same":
@@ -159,157 +187,271 @@ def create_four_turbine_farm(layout_type: str):
         turbine_types = ["IEA_10MW", "IEA_15MW", "IEA_15MW", "IEA_22MW"]
         turbine_dicts = [create_turbine(t_type) for t_type in turbine_types]
     else:
-        raise ValueError(f"Unknown layout_type: {layout_type}. Use 'same' or 'multiple'.")
-    
+        raise ValueError(
+            f"Unknown layout_type: {layout_type}. Use 'same' or 'multiple'."
+        )
+
     return {
         "layout_x": [0, 1248.1, 2496.2, 3744.3],
         "layout_y": [0.0, 0.0, 0.0, 0.0],
         "turbine_type": turbine_dicts,
     }
 
+
 def create_timeseries():
     """Create 3-step wind time series for testing."""
     return {
-        'wind_speeds': [10.091022491455078, 8.797999382019043, 10.307791709899902],
-        'wind_directions': [271.8246154785156, 268.6852111816406, 271.0501403808594],
-        'turbulence_intensities': [2.6189446449279785, 1.6507437229156494, 3.7473068237304688],
+        "wind_speeds": [10.091022491455078, 8.797999382019043, 10.307791709899902],
+        "wind_directions": [271.8246154785156, 268.6852111816406, 271.0501403808594],
+        "turbulence_intensities": [
+            2.6189446449279785,
+            1.6507437229156494,
+            3.7473068237304688,
+        ],
     }
 
+
 def create_timeseries_with_operating_flag():
     """Create time series with turbine availability flags."""
     return {
-        'wind_directions': [271.8246, 266.2015, 268.6852, 273.6164, 263.4558, 271.0501],
-        'wind_speeds': [10.09102, 10.23302, 8.797999, 9.662098, 9.78371, 10.30779],
-        'turbulence_intensities': [2.618945, 1.537051, 1.650744, 3.149817, 0.5826571, 3.747307],
-        'disable_turbines': np.logical_not(np.array([
-            [0, 1, 1, 1], [0, 1, 1, 1], [0, 1, 1, 1],  # T0 offline
-            [0, 1, 1, 1], [1, 1, 1, 1], [1, 1, 1, 1]   # All online
-        ]))
+        "wind_directions": [271.8246, 266.2015, 268.6852, 273.6164, 263.4558, 271.0501],
+        "wind_speeds": [10.09102, 10.23302, 8.797999, 9.662098, 9.78371, 10.30779],
+        "turbulence_intensities": [
+            2.618945,
+            1.537051,
+            1.650744,
+            3.149817,
+            0.5826571,
+            3.747307,
+        ],
+        "disable_turbines": np.logical_not(
+            np.array(
+                [
+                    [0, 1, 1, 1],
+                    [0, 1, 1, 1],
+                    [0, 1, 1, 1],  # T0 offline
+                    [0, 1, 1, 1],
+                    [1, 1, 1, 1],
+                    [1, 1, 1, 1],  # All online
+                ]
+            )
+        ),
     }
 
+
 def create_simple_wind_rose_energy_resource():
     """Create wind rose with Weibull distributions for 12 sectors."""
     # Sector probabilities and Weibull parameters
-    sector_probability = [0.0359715203597152, 0.0394868203948682, 0.0516739505167395, 0.0700015407000154,
-                         0.0836454708364547, 0.064348500643485, 0.0864319408643194, 0.117705101177051,
-                         0.151575701515757, 0.147379201473792, 0.100120501001205, 0.0516597505165975, 0.0359715203597152]
-    
-    weibull_a = [9.176929, 9.782334, 9.531809, 9.909545, 10.04269, 9.593921, 9.584007, 10.51499, 
-                11.39895, 11.68746, 11.63732, 10.08803, 9.176929]
-    
-    weibull_k = [2.392578, 2.447266, 2.412109, 2.591797, 2.755859, 2.595703, 2.583984, 2.548828,
-                2.470703, 2.607422, 2.626953, 2.326172, 2.392578]
+    sector_probability = [
+        0.0359715203597152,
+        0.0394868203948682,
+        0.0516739505167395,
+        0.0700015407000154,
+        0.0836454708364547,
+        0.064348500643485,
+        0.0864319408643194,
+        0.117705101177051,
+        0.151575701515757,
+        0.147379201473792,
+        0.100120501001205,
+        0.0516597505165975,
+        0.0359715203597152,
+    ]
+
+    weibull_a = [
+        9.176929,
+        9.782334,
+        9.531809,
+        9.909545,
+        10.04269,
+        9.593921,
+        9.584007,
+        10.51499,
+        11.39895,
+        11.68746,
+        11.63732,
+        10.08803,
+        9.176929,
+    ]
+
+    weibull_k = [
+        2.392578,
+        2.447266,
+        2.412109,
+        2.591797,
+        2.755859,
+        2.595703,
+        2.583984,
+        2.548828,
+        2.470703,
+        2.607422,
+        2.626953,
+        2.326172,
+        2.392578,
+    ]
 
     wind_direction = [0, 30, 60, 90, 120, 150, 180, 210, 240, 270, 300, 330, 360]
-    wind_speed = [3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25]
-    
+    wind_speed = [
+        3,
+        4,
+        5,
+        6,
+        7,
+        8,
+        9,
+        10,
+        11,
+        12,
+        13,
+        14,
+        15,
+        16,
+        17,
+        18,
+        19,
+        20,
+        21,
+        22,
+        23,
+        24,
+        25,
+    ]
+
     # Normalize probabilities and calculate frequency table
     sector_probability = np.array(sector_probability) / np.sum(sector_probability)
     wind_speeds = np.array(wind_speed)
     freq_table = np.zeros((len(wind_direction), len(wind_speed)))
-    
+
     for i_dir in range(len(wind_direction)):
         A, k = weibull_a[i_dir], weibull_k[i_dir]
         wind_speed_edges = np.arange(wind_speeds[0] - 0.5, wind_speeds[-1] + 1.5, 1.0)
-        
+
         # Weibull CDF: F(x) = 1 - exp(-(x/A)^k)
         cdf_edges = 1.0 - np.exp(-((wind_speed_edges / A) ** k))
         cdf_edges[wind_speed_edges < 0] = 0.0
-        
+
         freq = np.diff(cdf_edges)
         freq_table[i_dir, :] = (freq / freq.sum()) * sector_probability[i_dir]
 
-    wr = floris.WindRose(wind_directions=np.array(wind_direction), wind_speeds=wind_speeds, 
-                        ti_table=0.1, freq_table=freq_table)
+    wr = floris.WindRose(
+        wind_directions=np.array(wind_direction),
+        wind_speeds=wind_speeds,
+        ti_table=0.1,
+        freq_table=freq_table,
+    )
     return {"wind_data": wr}
 
+
 def create_wake_model_config():
     """Standard wake model configuration for testing."""
     return {
-        'enable_active_wake_mixing': False,
-        'enable_secondary_steering': False,
-        'enable_transverse_velocities': False,
-        'enable_yaw_added_recovery': False,
-        'model_strings': {
-            'combination_model': 'fls',
-            'deflection_model': 'none',
-            'velocity_model': 'gauss',
-            'turbulence_model': 'crespo_hernandez',
+        "enable_active_wake_mixing": False,
+        "enable_secondary_steering": False,
+        "enable_transverse_velocities": False,
+        "enable_yaw_added_recovery": False,
+        "model_strings": {
+            "combination_model": "fls",
+            "deflection_model": "none",
+            "velocity_model": "gauss",
+            "turbulence_model": "crespo_hernandez",
         },
-        'wake_velocity_parameters': {'gauss': {'ka': 0.0, 'kb': 0.04}},
+        "wake_velocity_parameters": {"gauss": {"ka": 0.0, "kb": 0.04}},
     }
 
+
 # ============================================================================ #
 # TEST FUNCTIONS
 # ============================================================================ #
 
+
 def test_floris_4wts(floris_config):
     """Test homogeneous 4-turbine farm. Max diff < 1e-6."""
-    wes_dir = test_path / "../examples/cases/windio_4turbines/wind_energy_system/system.yaml"
+    wes_dir = (
+        test_path / "../examples/cases/windio_4turbines/wind_energy_system/system.yaml"
+    )
     output_data = _run_floris(wes_dir)
-    powers_wifa = output_data['turbine_data.nc']['power'].values
+    powers_wifa = output_data["turbine_data.nc"]["power"].values
 
     floris_dict = floris_config.copy()
-    recursive_dict_update(floris_dict['farm'], create_four_turbine_farm("same"))
-    recursive_dict_update(floris_dict['flow_field'], create_timeseries())
-    recursive_dict_update(floris_dict['wake'], create_wake_model_config())
+    recursive_dict_update(floris_dict["farm"], create_four_turbine_farm("same"))
+    recursive_dict_update(floris_dict["flow_field"], create_timeseries())
+    recursive_dict_update(floris_dict["wake"], create_wake_model_config())
 
     fmodel = floris.FlorisModel(floris_dict)
     fmodel.run()
     powers_floris = fmodel.get_turbine_powers()
 
     # Use pytest.approx for better numerical comparison
-    assert powers_wifa == pytest.approx(powers_floris, rel=1e-6), f"Power outputs don't match within 1e-6 tolerance"
+    assert powers_wifa == pytest.approx(
+        powers_floris, rel=1e-6
+    ), f"Power outputs don't match within 1e-6 tolerance"
+
 
 def test_floris_simple_wind_rose(floris_config):
     """Test wind rose analysis. Max diff < 1e-2."""
-    wes_dir = test_path / "../examples/cases/simple_wind_rose/wind_energy_system/system.yaml"
+    wes_dir = (
+        test_path / "../examples/cases/simple_wind_rose/wind_energy_system/system.yaml"
+    )
     output_data = _run_floris(wes_dir)
-    powers_wifa = output_data['PowerTable.nc']['power'].values
+    powers_wifa = output_data["PowerTable.nc"]["power"].values
 
     floris_dict = floris_config.copy()
-    recursive_dict_update(floris_dict['farm'], create_four_turbine_farm("same"))
-    recursive_dict_update(floris_dict['wake'], create_wake_model_config())
+    recursive_dict_update(floris_dict["farm"], create_four_turbine_farm("same"))
+    recursive_dict_update(floris_dict["wake"], create_wake_model_config())
 
     fmodel = floris.FlorisModel(floris_dict)
     fmodel.set(**create_simple_wind_rose_energy_resource())
     fmodel.run()
     powers_floris = fmodel.get_turbine_powers()
 
-    assert powers_wifa == pytest.approx(powers_floris, rel=1e-2), f"Wind rose power outputs don't match within 1e-2 tolerance"
+    assert powers_wifa == pytest.approx(
+        powers_floris, rel=1e-2
+    ), f"Wind rose power outputs don't match within 1e-2 tolerance"
+
 
 def test_floris_timeseries_with_operating_flag(floris_config):
     """Test turbine availability flags. Max diff < 1e-2."""
-    wes_dir = test_path / "../examples/cases/timeseries_with_operating_flag/wind_energy_system/system.yaml"
+    wes_dir = (
+        test_path
+        / "../examples/cases/timeseries_with_operating_flag/wind_energy_system/system.yaml"
+    )
     output_data = _run_floris(wes_dir)
-    powers_wifa = output_data['turbine_data.nc']['power'].values
+    powers_wifa = output_data["turbine_data.nc"]["power"].values
 
     floris_dict = floris_config.copy()
-    recursive_dict_update(floris_dict['farm'], create_four_turbine_farm("same"))
-    recursive_dict_update(floris_dict['wake'], create_wake_model_config())
-    
+    recursive_dict_update(floris_dict["farm"], create_four_turbine_farm("same"))
+    recursive_dict_update(floris_dict["wake"], create_wake_model_config())
+
     fmodel = floris.FlorisModel(floris_dict)
     fmodel.set_operation_model("mixed")
     fmodel.set(**create_timeseries_with_operating_flag())
     fmodel.run()
     powers_floris = fmodel.get_turbine_powers()
 
-    assert powers_wifa == pytest.approx(powers_floris, rel=1e-2), f"Operating flag power outputs don't match within 1e-2 tolerance"
+    assert powers_wifa == pytest.approx(
+        powers_floris, rel=1e-2
+    ), f"Operating flag power outputs don't match within 1e-2 tolerance"
+
 
 def test_floris_multiple_turbines(floris_config):
     """Test mixed turbine types. Max diff < 0.1. Marked as slow test."""
-    wes_dir = test_path / "../examples/cases/windio_4turbines_multipleTurbines/wind_energy_system/system.yaml"
+    wes_dir = (
+        test_path
+        / "../examples/cases/windio_4turbines_multipleTurbines/wind_energy_system/system.yaml"
+    )
     output_data = _run_floris(wes_dir)
-    powers_wifa = output_data['turbine_data.nc']['power'].values
+    powers_wifa = output_data["turbine_data.nc"]["power"].values
 
     floris_dict = floris_config.copy()
-    recursive_dict_update(floris_dict['farm'], create_four_turbine_farm("multiple"))
-    recursive_dict_update(floris_dict['flow_field'], create_timeseries())
-    recursive_dict_update(floris_dict['wake'], create_wake_model_config())
+    recursive_dict_update(floris_dict["farm"], create_four_turbine_farm("multiple"))
+    recursive_dict_update(floris_dict["flow_field"], create_timeseries())
+    recursive_dict_update(floris_dict["wake"], create_wake_model_config())
 
     # Set reference height to mean hub height
-    mean_hub_height = np.mean([t["hub_height"] for t in floris_dict["farm"]["turbine_type"]])
-    floris_dict['flow_field']['reference_wind_height'] = mean_hub_height
+    mean_hub_height = np.mean(
+        [t["hub_height"] for t in floris_dict["farm"]["turbine_type"]]
+    )
+    floris_dict["flow_field"]["reference_wind_height"] = mean_hub_height
 
     fmodel = floris.FlorisModel(floris_dict)
     fmodel.run()
@@ -317,14 +459,17 @@ def test_floris_multiple_turbines(floris_config):
 
     max_diff = np.max(np.abs(powers_wifa - powers_floris)) / np.max(np.abs(powers_wifa))
     print(f"Max relative difference: {max_diff:.6f}")
-    assert powers_wifa == pytest.approx(powers_floris, rel=0.1), f"Mixed turbine power outputs don't match within 0.1 tolerance"
+    assert powers_wifa == pytest.approx(
+        powers_floris, rel=0.1
+    ), f"Mixed turbine power outputs don't match within 0.1 tolerance"
 
 
 # ============================================================================ #
 # PYTEST FIXTURES
 # ============================================================================ #
 
-@pytest.fixture(scope="session")  
+
+@pytest.fixture(scope="session")
 def test_output_dir():
     """Shared test output directory fixture with automatic cleanup."""
     output_dir = Path("output_test_floris")
@@ -332,9 +477,3 @@ def test_output_dir():
     yield output_dir
     if output_dir.exists():
         shutil.rmtree(output_dir)
-
-@pytest.fixture
-def floris_config():
-    """Default FLORIS configuration for all tests."""
-    return get_default_floris_dict()
-
diff --git a/wifa/floris_api.py b/wifa/floris_api.py
index 8cbf0bb..b5d5812 100644
--- a/wifa/floris_api.py
+++ b/wifa/floris_api.py
@@ -1,27 +1,28 @@
-from windIO import load_yaml
-from typing import Dict, List
+from pathlib import Path
+from typing import TYPE_CHECKING, Any, Dict, List
+
 import numpy as np
 import windIO
-from pathlib import Path
+from windIO import load_yaml
 
-from typing import TYPE_CHECKING
 if TYPE_CHECKING:
     from floris import FlorisModel
 
+
 def run_floris(yaml_input):
     """
     Run FLORIS model based on windIO YAML input and extract specified outputs.
-    
+
     Parameters
     ----------
     yaml_input: str, Path, or dict
         Path to the windIO YAML input file or a dictionary representing the input data.
-        
+
     Notes
     -----
-    WindIO is not yet suppoeted by FLOIS main branch. 
+    WindIO is not yet suppoeted by FLOIS main branch.
     Currently refer to https://github.com/lejeunemax/floris/tree/windIO.
-    
+
     Some features may not be fully supported, these include:
     - Subsetting of wind conditions (times_run, wind_speeds_run, directions_run)
     - Only "power", "thrust_coefficient", "axial_induction", and "turbulence_intensity"
@@ -38,35 +39,36 @@ def run_floris(yaml_input):
     elif isinstance(yaml_input, dict):
         windio_dict = yaml_input
     else:
-        raise TypeError(f"yaml_input must be a file path or dict, got {type(yaml_input)}")
+        raise TypeError(
+            f"yaml_input must be a file path or dict, got {type(yaml_input)}"
+        )
 
     fmodel = FlorisModel.from_windio(windio_dict)
     fmodel.run()
 
-    windio_dict['_context'] = ''
-    
-    with TrackedDict.from_parent(windio_dict, 'attributes') as attrs:
-        attrs.untrack('analysis')
-        attrs.untrack('flow_model')
+    windio_dict["_context"] = ""
 
-        out_specs = attrs['model_outputs_specification']
-        run_config = out_specs['run_configuration']
+    with TrackedDict.from_parent(windio_dict, "attributes") as attrs:
+        attrs.untrack("analysis")
+        attrs.untrack("flow_model")
+
+        out_specs = attrs["model_outputs_specification"]
+        run_config = out_specs["run_configuration"]
 
         slice_selection = None
-        times_run: Dict = run_config.get('times_run', None)
-        wind_speeds_run: Dict = run_config.get('wind_speeds_run', None)
-        directions_run: Dict = run_config.get('directions_run', None)
+        times_run: Dict = run_config.get("times_run", None)
+        wind_speeds_run: Dict = run_config.get("wind_speeds_run", None)
+        directions_run: Dict = run_config.get("directions_run", None)
 
         if times_run:
             all_occurences = times_run.get("all_occurences", False)
-            
+
             if all_occurences:
                 slice_selection = slice(None)
             else:
                 slice_selection = times_run.get("subset", None)
 
         elif wind_speeds_run and directions_run:
-            
             if not wind_speeds_run.get("all_values", False):
                 raise NotImplementedError(
                     f"Wind speed and direction subsets are not yet supported, got {wind_speeds_run.name} {wind_speeds_run}"
@@ -75,63 +77,67 @@ def run_floris(yaml_input):
                 raise NotImplementedError(
                     f"Wind speed and direction subsets are not yet supported, got {directions_run.name} {directions_run}"
                 )
-        
+
         output_dir = out_specs.get("output_folder", ".")
         Path(output_dir).mkdir(parents=True, exist_ok=True)
 
-        if ("turbine_outputs" in out_specs):
+        if "turbine_outputs" in out_specs:
             _read_turbines(fmodel, output_dir, out_specs, slice_selection)
 
-        if ("flow_field" in out_specs):
+        if "flow_field" in out_specs:
             _read_fields(fmodel, output_dir, out_specs, slice_selection)
 
+
 def _read_turbines(fmodel: "FlorisModel", output_dir, out_specs, slice_selection):
     """
     Extract turbine output data from FLORIS model and write to NetCDF file.
-    
+
     Args:
         fmodel: FlorisModel instance that has been run
         output_dir: Output directory path
         out_specs: TrackedDict containing turbine_outputs configuration from out_specs["turbine_outputs"]
         slice_selection: Time slice selection for filtering output
-    
+
     Returns:
         xarray.Dataset containing turbine output data
     """
-    import xarray as xr 
+    import xarray as xr
 
     turbine_outputs = out_specs["turbine_outputs"]
 
     output_vars: List[str] = turbine_outputs["output_variables"]
-    turbine_nc_filename = turbine_outputs.get("turbine_nc_filename", "turbine_outputs.nc")
-        
+    turbine_nc_filename = turbine_outputs.get(
+        "turbine_nc_filename", "turbine_outputs.nc"
+    )
+
     # Make a copy to avoid modifying the original list
     requested_vars = output_vars.copy()
-    data_vars = {}
+    data_vars: Dict[str, tuple[List[str], Any]] = {}
     coords = {}
-    
+
     # Determine the wind data type to understand output structure
     wind_data_type = None
     if fmodel.wind_data is not None:
         from floris.wind_data import WindRose, WindRoseWRG, WindTIRose
+
         if isinstance(fmodel.wind_data, WindTIRose):
-            wind_data_type = 'wind_ti_rose'
+            wind_data_type = "wind_ti_rose"
         elif isinstance(fmodel.wind_data, (WindRose, WindRoseWRG)):
-            wind_data_type = 'wind_rose'
-    
+            wind_data_type = "wind_rose"
+
     # Define handlers for each variable type
     # Format: variable_name -> (extractor_func, is_coord)
     variable_handlers = {
-        'time': (lambda: fmodel.get_metadata('wind_resource', {}).get('time'), True),
-        'x': (lambda: fmodel.layout_x, False),
-        'y': (lambda: fmodel.layout_y, False),
-        'z': (lambda: fmodel.get_turbine_layout(z=True)[2], False),
-        'power': (lambda: fmodel.get_turbine_powers(), False),
-        'thrust_coefficient': (lambda: fmodel.get_turbine_thrust_coefficients(), False),
-        'axial_induction': (lambda: fmodel.get_turbine_ais(), False),
-        'turbulence_intensity': (lambda: fmodel.get_turbine_TIs(), False),
+        "time": (lambda: fmodel.get_metadata("wind_resource", {}).get("time"), True),
+        "x": (lambda: fmodel.layout_x, False),
+        "y": (lambda: fmodel.layout_y, False),
+        "z": (lambda: fmodel.get_turbine_layout(z=True)[2], False),
+        "power": (lambda: fmodel.get_turbine_powers(), False),
+        "thrust_coefficient": (lambda: fmodel.get_turbine_thrust_coefficients(), False),
+        "axial_induction": (lambda: fmodel.get_turbine_ais(), False),
+        "turbulence_intensity": (lambda: fmodel.get_turbine_TIs(), False),
     }
-    
+
     # Process each requested variable
     unhandled_vars = []
     for var_name in requested_vars:
@@ -139,10 +145,10 @@ def _read_turbines(fmodel: "FlorisModel", output_dir, out_specs, slice_selection
             try:
                 extractor, is_coord = variable_handlers[var_name]
                 data = extractor()
-                
+
                 if data is not None:
                     data = np.asarray(data)
-                    
+
                     # Determine dimensions based on data shape
                     if is_coord:
                         # Coordinate variables are 1D
@@ -153,18 +159,31 @@ def _read_turbines(fmodel: "FlorisModel", output_dir, out_specs, slice_selection
                         # Data variables - infer dimensions based on shape and wind data type
                         if data.ndim == 1:
                             # 1D: turbine dimension only
-                            data_vars[var_name] = (['turbine'], data)
+                            data_vars[var_name] = (["turbine"], data)
                         elif data.ndim == 2:
                             # 2D: (time, turbine) - TimeSeries or flattened
-                            data_vars[var_name] = (['time', 'turbine'], data)
+                            data_vars[var_name] = (["time", "turbine"], data)
                         elif data.ndim == 3:
                             # 3D: (wind_direction, wind_speed, turbine) - WindRose
-                            data_vars[var_name] = (['wind_direction', 'wind_speed', 'turbine'], data)
+                            data_vars[var_name] = (
+                                ["wind_direction", "wind_speed", "turbine"],
+                                data,
+                            )
                         elif data.ndim == 4:
                             # 4D: (wind_direction, wind_speed, turbulence_intensity, turbine) - WindTIRose
-                            data_vars[var_name] = (['wind_direction', 'wind_speed', 'turbulence_intensity', 'turbine'], data)
+                            data_vars[var_name] = (
+                                [
+                                    "wind_direction",
+                                    "wind_speed",
+                                    "turbulence_intensity",
+                                    "turbine",
+                                ],
+                                data,
+                            )
                         else:
-                            print(f"Warning: Unexpected data shape for '{var_name}': {data.shape}; skipping.")
+                            print(
+                                f"Warning: Unexpected data shape for '{var_name}': {data.shape}; skipping."
+                            )
                             continue
                 else:
                     print(f"Warning: '{var_name}' data not available; skipping.")
@@ -172,96 +191,104 @@ def _read_turbines(fmodel: "FlorisModel", output_dir, out_specs, slice_selection
                 print(f"Warning: Failed to extract '{var_name}': {e}; skipping.")
         else:
             unhandled_vars.append(var_name)
-    
+
     # Warn about unhandled variables
     if unhandled_vars:
         print(f"Warning: Unhandled output variables: {unhandled_vars}")
-    
+
     # Build coordinates based on data structure
     final_coords = {}
-    
+
     # Add turbine coordinate
     n_turbines = len(fmodel.layout_x)
-    final_coords['turbine'] = np.arange(n_turbines)
-    
+    final_coords["turbine"] = np.arange(n_turbines)
+
     # Determine primary dimension and add appropriate coordinates
-    if 'time' in coords and coords['time'] is not None:
+    if "time" in coords and coords["time"] is not None:
         # TimeSeries case - use time coordinate
-        final_coords['time'] = coords['time']
+        final_coords["time"] = coords["time"]
         # Rename findex to time in data_vars
         renamed_data_vars = {}
         for var_name, (dims, data) in data_vars.items():
-            new_dims = tuple('time' if d == 'findex' else d for d in dims)
+            new_dims = ["time" if d == "findex" else d for d in dims]
             renamed_data_vars[var_name] = (new_dims, data)
         data_vars = renamed_data_vars
-    elif wind_data_type in ['wind_rose', 'wind_ti_rose']:
+    elif wind_data_type in ["wind_rose", "wind_ti_rose"]:
         # WindRose case - add wind direction and speed coordinates
         if fmodel.wind_data is not None:
-            final_coords['wind_direction'] = fmodel.wind_data.wind_directions
-            final_coords['wind_speed'] = fmodel.wind_data.wind_speeds
-            if wind_data_type == 'wind_ti_rose':
-                final_coords['turbulence_intensity'] = fmodel.wind_data.turbulence_intensities
+            final_coords["wind_direction"] = fmodel.wind_data.wind_directions
+            final_coords["wind_speed"] = fmodel.wind_data.wind_speeds
+            if wind_data_type == "wind_ti_rose":
+                final_coords[
+                    "turbulence_intensity"
+                ] = fmodel.wind_data.turbulence_intensities
     else:
         # Fallback: infer findex dimension
         n_findex = None
         for var_name, (dims, data) in data_vars.items():
-            if 'findex' in dims:
-                n_findex = data.shape[dims.index('findex')]
+            if "findex" in dims:
+                n_findex = data.shape[dims.index("findex")]
                 break
         if n_findex is not None:
-            final_coords['findex'] = np.arange(n_findex)
-    
+            final_coords["findex"] = np.arange(n_findex)
+
     # Create dataset
     ds = xr.Dataset(data_vars=data_vars, coords=final_coords)
-    
+
     # Add metadata attributes
-    ds.attrs['description'] = 'FLORIS turbine output data'
-    ds.attrs['source'] = 'FLORIS wind farm wake model'
-    
+    ds.attrs["description"] = "FLORIS turbine output data"
+    ds.attrs["source"] = "FLORIS wind farm wake model"
+
     # Add units and descriptions for known variables
     var_attrs = {
-        'x': {'units': 'm', 'long_name': 'turbine x-coordinate'},
-        'y': {'units': 'm', 'long_name': 'turbine y-coordinate'},
-        'z': {'units': 'm', 'long_name': 'turbine hub height'},
-        'power': {'units': 'W', 'long_name': 'turbine power output'},
-        'thrust_coefficient': {'units': '-', 'long_name': 'turbine thrust coefficient'},
-        'axial_induction': {'units': '-', 'long_name': 'turbine axial induction factor'},
-        'turbulence_intensity': {'units': '-', 'long_name': 'turbulence intensity at turbine'},
+        "x": {"units": "m", "long_name": "turbine x-coordinate"},
+        "y": {"units": "m", "long_name": "turbine y-coordinate"},
+        "z": {"units": "m", "long_name": "turbine hub height"},
+        "power": {"units": "W", "long_name": "turbine power output"},
+        "thrust_coefficient": {"units": "-", "long_name": "turbine thrust coefficient"},
+        "axial_induction": {
+            "units": "-",
+            "long_name": "turbine axial induction factor",
+        },
+        "turbulence_intensity": {
+            "units": "-",
+            "long_name": "turbulence intensity at turbine",
+        },
     }
-    
+
     for var_name, attrs in var_attrs.items():
         if var_name in ds.data_vars:
             ds[var_name].attrs.update(attrs)
-    
+
     # Apply time selection if specified
     if slice_selection is not None:
         # Determine which dimension to slice
-        if 'time' in ds.dims:
-            ds = ds.isel({'time': slice_selection})
-        elif 'findex' in ds.dims:
-            ds = ds.isel({'findex': slice_selection})
-        elif 'wind_direction' in ds.dims:
+        if "time" in ds.dims:
+            ds = ds.isel({"time": slice_selection})
+        elif "findex" in ds.dims:
+            ds = ds.isel({"findex": slice_selection})
+        elif "wind_direction" in ds.dims:
             # For WindRose, could slice wind_direction or wind_speed
             # This would need more sophisticated handling based on user intent
             print("Warning: slice_selection not applied for WindRose data")
-    
+
     # Write to NetCDF
     ds.to_netcdf(f"{output_dir}/{turbine_nc_filename}")
     print(f"Turbine output data written to {output_dir}/{turbine_nc_filename}")
-    
+
     return ds
-    
+
 
 def _read_fields(fmodel: "FlorisModel", output_dir, out_specs, slice_selection):
     """
     Extract flow field data from FLORIS model and write to NetCDF file.
-    
+
     Args:
         fmodel: FlorisModel instance that has been run
         output_dir: Output directory path
         out_specs: TrackedDict containing flow_field configuration from out_specs["flow_field"]
         slice_selection: Time slice selection for filtering output
-    
+
     Returns:
         xarray.Dataset containing flow field data
     """
@@ -271,19 +298,19 @@ def _read_fields(fmodel: "FlorisModel", output_dir, out_specs, slice_selection):
 
     # Check if flow field reporting is enabled
     report = flow_field.get("report", True)
-    
+
     if not report:
         return None
-    
+
     # Get output configuration
     output_vars: List[str] = flow_field["output_variables"]
     flow_nc_filename = flow_field.get("flow_nc_filename", "flow_field.nc")
-    
+
     # Get z_planes configuration
     z_planes = flow_field["z_planes"]
     z_sampling = z_planes["z_sampling"]
     xy_sampling = z_planes.get("xy_sampling", "grid")
-    
+
     # Determine z values to sample
     if z_sampling == "hub_heights":
         z_list = fmodel.core.farm.hub_heights
@@ -293,12 +320,12 @@ def _read_fields(fmodel: "FlorisModel", output_dir, out_specs, slice_selection):
         z_list = z_planes["z_list"]
     else:
         raise ValueError(f"Unknown z_sampling type: {z_sampling}")
-    
+
     # Get xy sampling parameters
     if xy_sampling == "grid":
         x_bounds = z_planes["x_bounds"]
         y_bounds = z_planes["y_bounds"]
-        
+
         # Support either Nx/Ny or dx/dy specification
         if "Nx" in z_planes and "Ny" in z_planes:
             Nx = z_planes["Nx"]
@@ -324,7 +351,7 @@ def _read_fields(fmodel: "FlorisModel", output_dir, out_specs, slice_selection):
 
     # Calculate horizontal planes for each z level
     data_vars = {}
-    
+
     # Prepare coordinate arrays
     n_findex = fmodel.core.flow_field.n_findex
     if slice_selection is not None:
@@ -334,7 +361,7 @@ def _read_fields(fmodel: "FlorisModel", output_dir, out_specs, slice_selection):
             findex_range = slice_selection
     else:
         findex_range = range(n_findex)
-    
+
     # For each z level, calculate the horizontal plane
     for z_idx, z_height in enumerate(z_list):
         # Calculate horizontal plane at this height
@@ -348,21 +375,27 @@ def _read_fields(fmodel: "FlorisModel", output_dir, out_specs, slice_selection):
             findex_for_viz=0,  # Will loop over findices if needed
             keep_inertial_frame=True,
         )
-        
+
         # Extract coordinates from first plane
         if z_idx == 0:
             # Get unique x and y coordinates
             x_coords = np.unique(horizontal_plane.df.x1.values)
             y_coords = np.unique(horizontal_plane.df.x2.values)
-            
+
             # Initialize data arrays
             if "u" in output_vars:
-                u_data = np.zeros((len(findex_range), len(z_list), len(y_coords), len(x_coords)))
+                u_data = np.zeros(
+                    (len(findex_range), len(z_list), len(y_coords), len(x_coords))
+                )
             if "v" in output_vars:
-                v_data = np.zeros((len(findex_range), len(z_list), len(y_coords), len(x_coords)))
+                v_data = np.zeros(
+                    (len(findex_range), len(z_list), len(y_coords), len(x_coords))
+                )
             if "w" in output_vars:
-                w_data = np.zeros((len(findex_range), len(z_list), len(y_coords), len(x_coords)))
-        
+                w_data = np.zeros(
+                    (len(findex_range), len(z_list), len(y_coords), len(x_coords))
+                )
+
         # Reshape plane data to grid
         # The dataframe has x1, x2, x3, u, v, w columns
         for t_idx, findex in enumerate(findex_range):
@@ -377,19 +410,19 @@ def _read_fields(fmodel: "FlorisModel", output_dir, out_specs, slice_selection):
                     findex_for_viz=findex,
                     keep_inertial_frame=True,
                 )
-            
+
             # Reshape to 2D grid (y, x)
             u_grid = horizontal_plane.df.u.values.reshape(len(y_coords), len(x_coords))
             v_grid = horizontal_plane.df.v.values.reshape(len(y_coords), len(x_coords))
             w_grid = horizontal_plane.df.w.values.reshape(len(y_coords), len(x_coords))
-            
+
             if "u" in output_vars:
                 u_data[t_idx, z_idx, :, :] = u_grid
             if "v" in output_vars:
                 v_data[t_idx, z_idx, :, :] = v_grid
             if "w" in output_vars:
                 w_data[t_idx, z_idx, :, :] = w_grid
-    
+
     # Create xarray dataset
     coords = {
         "time": list(findex_range),
@@ -397,22 +430,22 @@ def _read_fields(fmodel: "FlorisModel", output_dir, out_specs, slice_selection):
         "y": y_coords,
         "x": x_coords,
     }
-    
+
     if "u" in output_vars:
         data_vars["u"] = (("time", "z", "y", "x"), u_data)
     if "v" in output_vars:
         data_vars["v"] = (("time", "z", "y", "x"), v_data)
     if "w" in output_vars:
         data_vars["w"] = (("time", "z", "y", "x"), w_data)
-    
+
     ds = xr.Dataset(data_vars=data_vars, coords=coords)
-    
+
     # Add metadata
     ds.attrs["description"] = "FLORIS flow field output"
     ds.x.attrs["units"] = "m"
     ds.y.attrs["units"] = "m"
     ds.z.attrs["units"] = "m"
-    
+
     if "u" in output_vars:
         ds.u.attrs["units"] = "m/s"
         ds.u.attrs["long_name"] = "x-component of velocity"
@@ -422,9 +455,9 @@ def _read_fields(fmodel: "FlorisModel", output_dir, out_specs, slice_selection):
     if "w" in output_vars:
         ds.w.attrs["units"] = "m/s"
         ds.w.attrs["long_name"] = "z-component of velocity"
-    
+
     # Write to NetCDF
     ds.to_netcdf(f"{output_dir}/{flow_nc_filename}")
     print(f"Flow field data written to {output_dir}/{flow_nc_filename}")
-    
-    return ds
\ No newline at end of file
+
+    return ds