`TurbinePointCloud` silently scrambles per-turbine data when turbines have different hub heights

[ews-ffarella]

First of all, great library!

I think I have found a bug.

Summary

TurbinePointCloud states return permuted/duplicated ambient data across turbines as soon as the wind farm contains turbines with different hub heights. With uniform hub heights the results are correct. The corruption is silent — no warning or error — and affects all downstream results (ambient rotor variables, wakes, power, AEP).

Environment

• foxes 1.8.3 (reproduced on current main, commit d2449ce)
• numpy 2.x, xarray 2024.07, single-chunk engine (also reproduces with default engine)

Minimal reproducer

import numpy as np, pandas as pd, xarray as xr
import foxes, foxes.variables as FV, foxes.constants as FC
from foxes.core import Engine

ws = np.array([[10.0, 10.5, 9.5]] * 2)  # distinct value per turbine
sdata = xr.Dataset(
    coords={FC.STATE: pd.date_range("2000-01-01", periods=2, freq="1h"),
            FC.TURBINE: np.arange(3)},
    data_vars={"ws": ((FC.STATE, FC.TURBINE), ws),
               "wd": ((FC.STATE, FC.TURBINE), np.full((2, 3), 270.0))},
)

for hubs in ([100.0, 100.0, 100.0], [100.0, 100.0, 90.0]):
    states = foxes.input.states.TurbinePointCloud(
        data_source=sdata,
        output_vars=[FV.WS, FV.WD, FV.TI, FV.RHO],
        var2ncvar={FV.WS: "ws", FV.WD: "wd"},
        fixed_vars={FV.TI: 0.06, FV.RHO: 1.225},
    )
    farm = foxes.WindFarm()
    for i, h in enumerate(hubs):
        farm.add_turbine(
            foxes.Turbine(xy=np.array([i * 500.0, 0.0]), H=h,
                          turbine_models=["null_type"]),
            verbosity=0,
        )
    algo = foxes.algorithms.Downwind(
        farm, states, wake_models=[], rotor_model="centre",
        mbook=foxes.models.ModelBook(), verbosity=0,
    )
    with Engine.new("single", verbosity=0):
        got = algo.calc_farm()[FV.AMB_REWS].to_numpy()[0]
    print(f"H={hubs} -> AMB_REWS={got}, expected {ws[0]}")

Output:

H=[100.0, 100.0, 100.0] -> AMB_REWS=[10.  10.5  9.5], expected [10.  10.5  9.5]   # OK
H=[100.0, 100.0, 90.0]  -> AMB_REWS=[10.5 10.5 10. ], expected [10.  10.5  9.5]   # scrambled

With mixed hub heights the per-turbine data comes back as [T1, T1, T0] instead of [T0, T1, T2].

Root cause

In DatasetStates.calculate → _analyze_points (foxes/input/states/dataset_states.py), the turbine-point-cloud branch (3D FC.TURBINE hcoords) sets

_points_data["up"] = points          # L1257
_points_data["points_vary"] = False  # L1258

but leaves heights_vary unset. The subsequent height analysis then runs: with per-turbine differing z, pmax[2] - pmin[2] > 1e-4 holds, so heights_vary=True and uh, uh2h = np.unique(points[:, :, 2], return_inverse=True) are computed (L1271–1274).

Later, the height reconstruction intended for gridded (x, y, h) states executes (L1382):

d = d[:, points_data["uh2h"], :].reshape(shp)

For a turbine point cloud, axis 1 of d is the turbine axis (the data was already returned per turbine, verbatim, by TurbinePointCloud.interpolate_data's exact-match branch). Re-indexing it with uh2h (inverse indices into the unique-heights array) permutes/duplicates turbines by their hub-height rank: H=[100, 100, 90] → unique [90, 100] → uh2h=[1, 1, 0] → data served as [T1, T1, T0].

Suggested fix

The turbine-cloud branch carries per-turbine coordinates verbatim and needs no height reconstruction — mark heights as non-varying there:

                 ):
                     _points_data["up"] = points
                     _points_data["points_vary"] = False
+                    # turbine point clouds carry per-turbine coordinates verbatim;
+                    # the uh2h height reconstruction must not re-index the turbine axis
+                    _points_data["heights_vary"] = False
+                    _points_data["uh"] = None
+                    _points_data["uh2h"] = None
                 elif np.max(pmax - pmin) > 1e-4:

With this patch the reproducer returns correct data for uniform and mixed hub heights (tested [100,100,100], [100,100,90], [120,100,90]); the existing uniform-height behavior is unchanged. Happy to open a PR.

[SchmJo]

Thanks for the issue! Can you please pull the latest dev branch and try your code there? The above minimal example returns

H=[100.0, 100.0, 100.0] -> AMB_REWS=[10.  10.5  9.5], expected [10.  10.5  9.5]
H=[100.0, 100.0, 90.0] -> AMB_REWS=[10.  10.5  9.5], expected [10.  10.5  9.5]

for me. Do you agree that the issue is solved on dev, i.e., in the upcoming release?

Thanks

[ews-ffarella]

I will check it out next week. Thanks for the reply

Btw, I compared py-wake and your repo and get machine precision parity, with very complex sector management rules and everything. Great work

If you want, I have introduced a caching class as well, to resolve on a grid. Basically rews is cached and interpolated from grid corners. Speeds up a lot ;) need to test it on blockage scenarios, because my throughput is barely 1 state per second 😀 ( I have 25 years time series and 10k plus states )

[ews-ffarella]

Thanks for pointing me at dev — I tested it. Short version: dev does fix this issue, but I can't move to it yet, because in fixing it two other things appear in the TurbinePointCloud path. Both are reproducible with the script at the end (self-contained, no data files).

Builds compared:

	upstream main d2449ce	dev 63c965c	main + the 2-line patch from this issue
A. per-turbine data served verbatim, 2-D layout, uniform hub heights, WD varying	OK	FAIL	OK
B. per-turbine data served verbatim, mixed hub heights (this issue)	FAIL	OK	OK
C. scaling, 20 turbines + Jensen, uniform H, 2-D layout	linear	quadratic, OOM	linear

---

B — this issue is fixed on dev

Confirmed for H = [100,100,90], [120,100,90], [90,100,120] and uniform. main returns [10.5, 10.5, 10.0] for [100,100,90]; dev returns the input verbatim. Presumably from 6442fb8 "fixing issues with point cloud states". 👍

---

A — on dev, a 2-D layout corrupts TurbinePointCloud data (no mixed hub heights needed)

This is the one I'd flag first, because it is the same class of failure as this issue — silent, plausible-looking wrong numbers — but with a much wider trigger.

Invariant under test: with no wake models, AMB_REWS[s, t] must equal the ws supplied for turbine t in state s, exactly. The support points are the evaluation points, so nothing should be interpolated, and the ambient wind certainly cannot depend on where the turbines sit.

dev, 6 turbines, 50 states, uniform H = 100 m, no wakes:

OK    row,  WD constant 270              max|AMB_REWS - input| = 0
OK    row,  WD varies per state only     max|AMB_REWS - input| = 3.6e-15
OK    row,  WD varies per state+turbine  max|AMB_REWS - input| = 3.6e-15
OK    grid, WD constant 270              max|AMB_REWS - input| = 0
FAIL  grid, WD varies per state only     max|AMB_REWS - input| = 5.443     <-- m/s
FAIL  grid, WD varies per state+turbine  max|AMB_REWS - input| = 6.480     <-- m/s

row = 6 turbines collinear on the x-axis; grid = the same 6 turbines at (i%5, i//5) * 500 m. Same input data in every row of that table.

Trigger, empirically: a non-collinear layout AND a wind direction that varies across states. Both main and main + the 2-line patch pass all six cases. Since "2-D layout with a varying wind direction" describes essentially every real wind farm and every real time series, this isn't an edge case on dev.

My guess at the mechanism (not verified in a debugger) is that it shares a root cause with C below: the per-state downwind reordering of the turbine axis makes the point set look state-dependent, and the reconstruction then puts values back on the wrong turbines. A collinear row survives that permutation; a 2-D layout doesn't.

---

C — O(n_states²) time and memory in the wake stage on dev

20 turbines, uniform hub heights (so this is isolated from this issue), 2-D layout, Jensen_linear_k007, single core:

n_states	main d2449ce	dev 63c965c	main + patch
500	0.08 s / 0.16 GB	0.53 s / 0.24 GB	0.08 s / 0.16 GB
1000	0.12 s / 0.16 GB	2.30 s / 0.48 GB	0.12 s / 0.16 GB
2000	0.24 s / 0.17 GB	8.09 s / 1.42 GB	0.23 s / 0.17 GB
4000	0.48 s / 0.19 GB	27.72 s / 5.19 GB	0.42 s / 0.19 GB
8000	1.01 s / 0.23 GB	MemoryError	1.08 s / 0.23 GB

~4x time and 2x memory per doubling. main is linear here, so this looks like a regression rather than a pre-existing cost. My production case is a 20-year hourly series (175,320 states), where dev requests a single 22.9 GiB array:

File ".../foxes/input/states/dataset_states.py", line 1641, in calculate
    d = d[:, points_data["up2p"], :].reshape(shp)
numpy._core._exceptions._ArrayMemoryError:
    Unable to allocate 22.9 GiB for an array with shape (175320, 8766, 2)

Where I think it comes from

_analyze_points(has_p, has_h, hcoords=None) still has the turbine-point-cloud short-circuit:

if (hcoords is not None and FC.TURBINE in hcoords
        and len(hcoords[FC.TURBINE].shape) == 3):
    _points_data["up"] = points
    _points_data["points_vary"] = False          # cheap path
elif np.any(pmax - pmin > 1e-4):
    _points_data["up"], _points_data["up2p"] = np.unique(
        points.reshape(n_states * n_pts, 3), axis=0, return_inverse=True)
    _points_data["points_vary"] = True           # gridded-states path

but 6442fb8 changed the FC.TURBINE call site from _analyze_points(has_p, has_h, hcoords) to _analyze_points(has_p, has_h). On dev no caller passes hcoords — all three call sites (≈ lines 1600, 1609, 1613) are two-argument — so that first branch is now unreachable, and the related guard at ≈ line 1578 is commented out.

Adding a print to _analyze_points on dev, a 20-turbine / 8,766-state farm calculation hits it twice:

[DBG] points.shape=(8766, 20, 3) spread=0     points_vary=False   <- ambient rotor stage
[DBG] points.shape=(8766, 20, 3) spread=2016  points_vary=True    <- wake stage

spread=2016 is exactly the farm's x-extent: in the wake stage the turbine axis is in per-state downwind order, so points[:, i] is a different turbine in different states and the point set looks state-dependent. That selects the np.unique(... n_states * n_pts ...) reconstruction, whose intermediate is (n_states, n_states * n_turbines, n_vars).

If that reading is right, then restoring an explicit "these points are the turbines, in whatever order the algorithm has them" signal — rather than inferring point variability from coordinate spread across states — would fix A and C together. I'm happy to open a PR if you tell me which shape you'd prefer that signal to take; I didn't want to guess at the design of the reworked point_cloud_data.py.

---

Reproducer

foxes_dev_branch_findings.py — attached / pasted below. No data files, no network. Run it against any checkout:

PYTHONPATH=/path/to/foxes python foxes_dev_branch_findings.py

It prints A, B and C for whichever build it imports, so the three columns of the first table are three runs of the same script.

foxes_dev_branch_findings.py

[SchmJo]

Thank you very much for this interesting issue and the detailed analysis. I just merged a mypy-related branch into dev, including a fix for this issue. On my laptop your script now gives the following result:

A — verbatim ambient data (no wakes, uniform H=100 m)
  OK   row , WD constant 270                max|AMB-input| = 0
  OK   row , WD varies per state            max|AMB-input| = 3.553e-15
  OK   row , WD varies per state+turbine    max|AMB-input| = 3.553e-15
  OK   grid, WD constant 270                max|AMB-input| = 0
  OK   grid, WD varies per state            max|AMB-input| = 3.553e-15
  OK   grid, WD varies per state+turbine    max|AMB-input| = 3.553e-15
  => 2-D-layout corruption: ABSENT on this build

B — issue #63: mixed hub heights (no wakes)
  OK   H=[100, 100, 100] -> [10.  10.5  9.5]  (expect [10.  10.5  9.5])
  OK   H=[100, 100, 90] -> [10.  10.5  9.5]  (expect [10.  10.5  9.5])
  OK   H=[120, 100, 90] -> [10.  10.5  9.5]  (expect [10.  10.5  9.5])
  OK   H=[90, 100, 120] -> [10.  10.5  9.5]  (expect [10.  10.5  9.5])
  => issue #63: FIXED on this build

C — scaling: 20 turbines, UNIFORM H, 2-D layout, Jensen_linear_k007
  n_states=500        0.04 s   peak RSS   0.16 GB
  n_states=1000       0.05 s   peak RSS   0.17 GB
  n_states=2000       0.09 s   peak RSS   0.18 GB
  n_states=4000       0.17 s   peak RSS   0.20 GB
  n_states=8000       0.35 s   peak RSS   0.23 GB

Would you say that the issue can now be closed? Cheers