perf: reuse triangulation/point-location in grid interpolation (drop Dask) - #53
Merged
Conversation
…Dask)
griddata rebuilt the Delaunay triangulation of the SOURCE points AND re-ran the
point-location of the (fixed) TARGET points on every call, so interpolating the
~60 vertical levels of the 3D CO2 field -- and the per-field 2D emissions/albedo
loops -- repeated that work N times. The 3D CO2 step wrapped this in a
dask.distributed LocalCluster whose process startup + pathological ~150 s
teardown dwarfed the ~seconds of actual interpolation.
Add ocp_tool/interp_utils.ReusableGridInterp: builds the Delaunay + KD-tree once
and the target-point barycentric weights once, so each subsequent field is just
a gather + weighted sum. Use it in:
- co2_interpolation: horizontal loop (was per-level griddata + Dask); the
vertical interp1d loop is likewise collapsed to a single vectorised axis-0
interp1d (after the horizontal nearest-fill there are no NaNs, so every
column shares the full level set). Dask removed entirely.
- field_interpolation: the 2D CO2-emissions (linear) and albedo (nearest)
per-field loops reuse one triangulation / KD-tree across fields.
Results are bit-identical to griddata to interpolation precision (barycentric
weights match LinearNDInterpolator; max abs diff ~1e-19), and bit-for-bit after
GRIB packing (grib_compare: identical). Measured on TCO95/feomdyn: the CO2
interpolation drops ~101 s -> ~6.7 s (~15x), and a full OCP-Tool run
211.6 s -> 138.6 s. The win grows with resolution (the triangulation build,
previously repeated, is amortised once).
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Problem
scipy.interpolate.griddatarebuilds the Delaunay triangulation of the source points and re-runs the point-location of the (fixed) target points on every call. Interpolating the ~60 vertical levels of the 3D CO2 field — and the per-field 2D emissions/albedo loops — repeats that identical work N times. The 3D CO2 step then wrapped it in adask.distributedLocalClusterwhose process startup + a pathological ~150 s worker teardown dwarfed the few seconds of actual interpolation.Profiling (TCO95/feomdyn) showed the CO2 interpolation dominated by exactly this:
LinearNDInterpolator.__call__re-doingfind_simplexat ~2.7 s × 60 levels.Fix
New
ocp_tool/interp_utils.ReusableGridInterp: builds the Delaunay + KD-tree once, and the target-point barycentric weights once — so each subsequent field is just a gather + weighted sum. Applied in:co2_interpolation— horizontal loop (was per-levelgriddata+ Dask). The verticalinterp1dper-point loop is collapsed to a single vectorisedinterp1d(axis=0)(after the horizontal nearest-fill there are no NaNs, so every column shares the full level set). Dask removed entirely.field_interpolation— the 2D CO2-emissions (linear,fill_value=0) and albedo (nearest) per-field loops reuse one triangulation / KD-tree across fields.Correctness
Bit-identical to
griddatato interpolation precision (the barycentric weights are the same onesLinearNDInterpolatoruses; max abs diff ~1e-19), and bit-for-bit after GRIB packing — a full-pipeline A/B (old Dask path vs new) givesgrib_compare: IDENTICAL.Impact (TCO95/feomdyn)
run_ocp_tool.py: 211.6 s → 138.6 sdask/distributeddependency in the interpolation path.