Repository navigation
Performance: Why do large equilibrium splines slow down the code? #376
Description
Activity
- added 7 commits that reference this issue
on Aug 14, 2026 Answer: spline evaluation is not the bottleneck — adaptive step count is
I measured this end to end on
examples/DIIID-like_ideal_example(ideal path,kinetic_factor = 0),
single-threaded, Vern9,eulerlagrange_tolerance = 1e-10unless stated otherwise. Short version:The hint-based spline lookup already does exactly what the issue hoped — evaluation cost is
flat in knot count. There is nothing left to win there. The mpsi slowdown is entirely
adaptive step-count growth: the integrator's accepted step size is slaved to the local knot
spacing, sonstepscales with the number of knots in the regions where steps concentrate.Two things worth saying up front, because they set what is actionable. Refining the grid genuinely
does make the integrand rougher — the coefficient matrices carry a node-error floor that grid
refinement amplifies rather than reduces (§6) — but when I repaired that directly it recovered
only 11–20% of the steps (§7). The dominant lever is much simpler: the two-pass auto grid reaches
the same answer with 558 knots and a third of the steps of a fixed mpsi=1024 grid.1. Spline evaluation is knot-count-independent
CubicSeriesInterpolant, 676ComplexF64series (mirrorsfmats/kmats/gmatsatnumpert = 26),
non-uniformVectorgrid, in-place eval, persistentRef(1)hint, 16k queries:npts mono ns/ev stage ns/ev no-hint mono range grid 3-sep stage B/eval 33 880 928 890 889 883 0 129 898 913 915 925 893 0 513 898 906 930 890 896 0 1025 887 889 937 991 1020 0 2049 892 909 948 901 935 0 Flat from 33 to 2049 knots, zero allocations. The hint buys only ~5% here because the 676-series
payload dominates the lookup entirely. A single-series control (theq_splineanalogue) runs
~2 ns/eval hinted and 5→13 ns unhinted across the same range — the expected O(log n), and
negligible. Soprecompute_transpose!, extra hints, and mergedSeriesInterpolants were measured
and dropped: there is no knot-length sensitivity left in the evaluation path to remove.2. What actually grows: the number of ODE steps
mpsi ladder, warm runs (second
main()call in one process, JIT excluded):mpsi knots EL integration (s) nstep ms/step ballooning (s) et[1] 128 129 9.3 1190 7.8 1.7 1.391 256 257 10.0 1753 5.7 3.7 0.376 512 513 18.5 3010 6.2 8.2 0.783 1024 1025 32.4 5683 5.7 18.7 0.803 auto 558 16.9 2542 6.6 17.4 0.801 EL time ∝ nstep, and ms/step is flat — per-step cost is mpsi-independent (Npert-sized LAPACK
dominates it, not interpolation). The whole slowdown is the step count.3. The steps pile up where the knots pack, not where the physics is
Accepted steps binned by ψ:
region 128 256 512 1024 auto ψ < 0.1 (axis packing) 553 964 1891 3806 581 0.985–0.9935 (edge pack) 117 157 217 418 227 q = 2 surface ± 63 62 64 67 112 q = 3 surface ± 59 60 58 59 87 Singular-surface neighbourhoods are mpsi-flat — those steps are physics-controlled, as they
should be. All the growth is in the near-axis region where the log-family grid packs hardest.Within ψ < 0.1, knot counts 73/145/289 (mpsi 128/256/512) give 553/964/1891 steps, and the ratio
dψ/knotΔ has p25/median/p75 = 0.06/0.12/0.21 at every mpsi — about 6.5 accepted steps per knot
interval, invariant. If steps were tolerance- or physics-limited, dψ would be mpsi-independent;
instead it halves whenever the knot spacing halves.Two discriminators confirm the mechanism at mpsi = 512:
grid_type = "uniform"(un-packs the axis): near-axis steps 1891 → 559, total nstep
3010 → 2237. Steps follow knot density, not axis physics. (This is a mechanism probe, not a
config recommendation — the uniform grid is physically under-resolved near axis, et[1] = 0.993.)eulerlagrange_tolerance = 1e-8(100× looser, same grid): accepted steps 3977 → 2116
with identical physics (et[1] unchanged to 8 digits). The step size is genuinely error-controlled
— but at an effective order well short of Vern9's, see §4.
So the step size is genuinely error-controlled — but the error magnitude per unit ψ is set by
knot-scale roughness in the interpolated F/K/G coefficients, not by the physics. §6 is about where
that roughness comes from, because "it's the same smooth function, just sampled more finely" is
the natural expectation and it turns out to be false.4. Tolerance sweep — how tight does
eulerlagrange_toleranceactually need to be?Five runs at mpsi = 512 (513 knots), everything but the tolerance held fixed. Δ' deviations are
relative to the 1e-12 point;Δ'diagis the worst per-surface relative error on the
singular/delta_prime_matrixdiagonal.tol accepted steps saved EL (s) et[1] Δ'diag rel 1e-6 1159 929 10.1 0.78255904 1.7e-02 1e-7 1565 1240 12.2 0.78255903 3.3e-03 1e-8 2116 1656 14.7 0.78255906 6.9e-04 1e-10 3977 3018 18.9 0.78255905 6.2e-06 1e-12 9620 7297 32.9 0.78255905 0 (ref) (
integration/nstep_totalis accepted steps;integration/nstepis saved snapshots at
save_interval = 3. Rejected steps are not recorded anywhere.)Three things worth knowing:
- et[1] is not a discriminator — identical to 8 significant digits across six decades of
tolerance. Δ' is what actually degrades, and it degrades fast below 1e-8 (the worst surface is
the q = 5 row: −1460 at 1e-6, −1332 at 1e-7, −1321 at 1e-8, −1319.5 at 1e-10). 1e-8 — the
ForceFreeStatesControldefault — is the sweet spot: 1.8× fewer steps than 1e-10 for 7e-4
relative Δ' error. - Vern9 never achieves its formal order. Per-decade exponents of h ∝ tol^α: 0.130, 0.131,
0.137, then 0.192 for the last two decades. A 9th-order method on a smooth integrand should
give α = 1/(p+1) = 0.100, or 0.111 under error-per-unit-step. Measured α is outside both
everywhere and degrades sharply at tight tolerance — an error floor that is not truncation
error. Vern7 over the same range gives α = 0.131, comfortably inside its own formal band
(0.125–0.143): the 7th-order method behaves as advertised, the 9th-order one does not. §6
explains why. - EL wall time is sub-linear in step count: 1e-10 → 1e-8 cuts accepted steps 1.88× but EL time only
1.29×. There is a ~7 s tolerance-independent floor in that phase on this deck. And the whole
warm run is ~220 s regardless of tolerance — at mpsi = 512 with the NTV section active, EL
integration is under 10% of the run. Tolerance relaxation is real but second-order here; it
matters where EL dominates (large mpsi, no kinetic post-processing).
5. Integrator order — Vern9 is the wrong tool here
If steps cannot span a knot interval anyway, Vern9's 16 stages per step are wasted. I swapped all
sixVern9()call sites insideForceFreeStatesforVern7()(10 stages), left the
equilibrium-side ones alone, and reran at mpsi = 512 (edits reverted afterwards):method tol accepted steps EL (s) ms/step et[1] Δ'diag rel Vern9 1e-10 3977 18.9 4.75 0.782559048 6.2e-06 Vern7 1e-10 5355 15.3 2.86 0.782559048 2.5e-04 Vern9 1e-8 2116 14.7 6.95 0.782559057 6.9e-04 Vern7 1e-8 2927 10.9 3.72 0.782559050 2.5e-04 Vern7 takes ~40% more steps but each step is 1.7× cheaper, so EL wall time drops ~20% at 1e-10 and
~26% at 1e-8. Vern7 at 1e-8 dominates Vern9 at 1e-8 on both axes — faster and closer to the
converged Δ' — and is 42% faster than the deck's current Vern9 @ 1e-10.One caveat before anyone acts on this: Vern7's Δ' deviation is 2.47e-4 at 1e-10 and 2.50e-4 at
1e-8, i.e. tolerance-independent, so it is not step-truncation error. Something in the
propagator/matching path carries a method-dependent offset that Vern9 does not. It is small, but
it should be understood rather than absorbed.6. Refining the grid makes the integrand measurably rougher
The integrand is the same smooth function only down to the accuracy of the node data. Below that,
adding knots does not add information; it just interpolates the error more finely.Diagnostic: second divided differences of
matrices/ideal/F,K,G(the actual EL integrand splines)
at their own knots. Two statistics —r1, the knot-to-knot autocorrelation of f'' (+1 = smooth,
−2/3 = white noise), and |f''|/|f|.q(ψ)on the same grid is the control.region quantity mpsi=256 mpsi=512 mpsi=1024 ψ<0.1 K: r1 / |f''|/|f| −0.42 / 1.01e7 −0.45 / 4.53e7 −0.56 / 1.62e8 ψ<0.1 F: r1 / |f''|/|f| +0.91 / 4.48e3 +0.67 / 4.10e3 −0.32 / 4.65e3 0.3–0.7 F: r1 / |f''|/|f| +0.91 / 1.730 +0.96 / 1.712 +0.97 / 1.705 0.3–0.7 G: r1 / |f''|/|f| +0.83 / 3.66 +0.86 / 3.87 −0.18 / 3.99 all q control r1 +0.88…+0.94 +0.94…+0.97 +0.97 Near-axis K grows 4.5× then 3.6× per doubling — the Δ⁻² of a fixed-amplitude node error — while
r1 marches toward the white-noise value. Meanwhile mid-plasma F converges beautifully
(1.730 → 1.712 → 1.705): where the grid is still coarse relative to the noise, the physics
dominates and your intuition holds exactly. The noise floor simply spreads outward with mpsi —
axis first, then edge, then mid-plasma.Three candidate sources tested and ruled out, all at fixed mpsi = 512:
knob range tested effect on accepted steps etol(field-line integration tolerance)1e-8 → 1e-12 4038 → 3930 (2.7%) mtheta(poloidal resolution)256 → 1024 3977 → 3983 (0.3%) input equilibrium 257×257 EFIT vs analytic Solovev Solovev diverges the same, or worse So it is not the ODE tolerance, not the poloidal quadrature, and not the finite resolution of the
input reconstruction. It is generated inside the equilibrium → coefficient pipeline.Localised. A local degree-4 fit residual over a 9-knot window measures the node-error amplitude
directly (smooth data → falls like Δ⁵; a floor → flat). Mid-plasma:matrix built from resid @512 resid @1024 r1 @512 r1 @1024 A, B, D g22, g23, g33 ~6e-08 ~1.3e-08 +0.972 +0.987 C + g31, jtheta·imat 8.90e-06 1.21e-05 +0.694 −0.401 E + g31, jtheta·imat, q1 8.74e-06 1.17e-05 +0.547 −0.400 H + g31, jtheta·imat 1.11e-05 1.41e-05 +0.837 +0.022 F = F̃ − D†A⁻¹D clean inputs 2.48e-07 1.96e-07 +0.957 +0.973 K = E − K†A⁻¹C dirty inputs 1.22e-05 1.56e-05 +0.916 +0.593 G = H − C†A⁻¹C dirty inputs 1.01e-05 1.70e-05 +0.856 −0.176 The floor enters at C, E, H at the ~1e-5 relative level. A/B/D are clean at 1e-8 and still
converging, exactly as smooth data should. The derived-matrix algebra is not the culprit either —
F goes through the sameA⁻¹path and stays clean at 2e-7 because its inputs are clean; K and G
are dirty only because E, C and H are.Traced to its source: the dirty matrices are the ones using g11/g12/g31, and per
Fourfit.jl:84-90those are built from ψ-derivatives (partials[2]) of the per-surface
geometry, while the clean matrices use only θ-derivatives (partials[3]) on the fixed
mthetagrid. The rawsplines/rzphinode data shows why that matters —nu,offsetand
rcoordssit on a ψ-direction floor that is flat in the grid spacing, whilejacon the same
grid converges:quantity resid @256 resid @512 resid @1024 nu (ψ<0.1) 2.41e-05 2.11e-05 2.12e-05 offset (ψ<0.1) 6.39e-04 5.37e-04 5.60e-04 jac (0.3–0.7) 2.22e-07 4.29e-08 3.26e-09 Each flux surface is traced independently, so its geometry carries an error that refining ψ does
not reduce; differentiating that across surfaces gives ε/Δψ. The 1D profiles are not involved
— q, mu0p, 2piF and dVdpsi all converge, clearingq1,p1andjtheta.7. How much of the slowdown does this actually explain? Only 11–20%
Rather than assume the roughness drives the step count, I repaired the amplification and measured
it: replace only the ψ-derivative channel with derivatives of a fit over a coarser ψ subgrid,
leaving values, θ-derivatives,jacand the coordinate mapping untouched.The repair does exactly what it was designed to do — mid-plasma at mpsi=1024, C's residual goes
1.21e-05 → 5.30e-07 and its r1 −0.401 → +0.915; G goes 1.70e-05 → 9.38e-07 and −0.176 →
+0.942; near-axis K's |f''|/|f| drops 1.62e8 → 2.10e7. Every autocorrelation flips positive,
everywhere. The matrices are genuinely smooth afterwards.The step count barely moves:
grid accepted steps, stock with repair change mpsi=512 3977 3551 −11% mpsi=1024 7638 6081 −20% So this partial repair is worth 11–20% — it caps the ψ-derivatives of data that is still
noisy, rather than removing the noise. The next section measures what removing it is worth.Two candidate origins of ε are also ruled out, so this is not ODE error control: the
etolsweep
above, and the hard-codedabstol=1e-8atDirectEquilibrium.jl:294(onlyreltolis
configurable, andu0starts at zeros with its components becoming ν/offset/r², so it should
bind wherereltolcannot) — tightening it to 1e-14 gives 3977 → 4071 at mpsi=512 and
7638 → 7474 at mpsi=1024. Null.Where the rest lives — measured. Partial repairs recover little because they never remove ε
itself. Five ladders varying input, fill site and grid one at a time show what removing it is
worth:case input construction geometry ε across the ladder step ratios DIII-D efitEFIT g-file field-line trace flat 1.8e-6 → 1.3e-6 1.72 / 1.92 DIII-D efit_by_inversionEFIT + contours inversion rising 1.3e-6 → 5.6e-6 1.76 / 1.96 Solovev solanalytic field-line trace flat 0.19, r1 = −0.65 1.72 / 1.83 LAR tj_analyticanalytic inversion converging 1.2e-7 → 5.2e-9 1.14 / 1.10 LAR + DIII-D packed grid analytic inversion converging 5.9e-7 → 2.9e-8 1.10 / 1.15 Clean geometry ⟺ flat step count, across two constructions, two inputs and two grids. Where the
surface data converges, quadrupling the knots costs ~25% more steps; where it sits on a floor, 3.3×.Two results worth highlighting. Solovev is analytic — perfectly smooth input — yet the standard
field-line construction hands the metricnunode data that is flat at 19% relative with
r1 = −0.65, i.e. pure white noise. Each surface is traced independently, splined on that surface's
own solver-chosen abscissae and resampled to the common θ grid, so the remap error is uncorrelated
between neighbours; that is where ε is manufactured, and it explains why neitheretolnor
abstoltouches it. And giving the clean LAR case the DIII-D axis-packed grid keeps it flat, so
grid packing is not the driver — those near-axis steps are noise-chasing, not physics.So the premise "same smooth function, just sampled more finely" is false of the data in a real and
measurable way — but only partly responsible. The practical consequence is in the recommendations:
the largest lever is not making the coefficients smoother, it is not placing the knots there.Recommendations
-
Use the two-pass auto grid (
mpsi = 0withpsi_accuracy) — already the example default,
and by some distance the biggest measured lever. It reaches et[1] = 0.801 with 558 knots and
2542 steps; a fixed mpsi=1024 log grid needs 7638 steps to reach 0.803. 3× the work for a
0.3% difference — that is the bulk of the "unnecessary" steps this issue is about, and it is
fixable today with one config line. -
Revisit the per-deck
eulerlagrange_tolerance. The struct default is1e-8; the DIII-D
decks pin1e-10and the LAR decks1e-12. On the evidence of §4 those tighter values buy
Δ' accuracy this case does not need, at 1.8–4× the step count. This is a per-deck edit on its
own branch with a regression-harness run, not a change to the struct default. -
Switch the ForceFreeStates integrator to Vern7 — measured in §5, ~20–26% off the EL phase,
with better Δ' than Vern9 at the same tolerance. Gated on explaining the tolerance-independent
2.5e-4 Δ' offset first; then its own branch and a regression-harness run. -
Make the surface-to-surface θ parametrization consistent — this is the main lever. Each
surface is currently traced independently and splined on its own solver-chosen SFL-angle
nodes before being resampled to the common θ grid (DirectEquilibrium.jl:500-525), so the
remap error is uncorrelated between neighbours. Sampling every surface at the common abscissae
(dense output plus a root-solve per node) would make that error smooth in ψ instead of white.
The ladders above bound the prize: clean geometry gives ~1.1× step growth per mpsi doubling
instead of ~1.9×.I tried the tidier-looking alternative first — deriving the ψ-derivatives from metric
identities — and measured it: it is worth ~2%, not 20%. ∇ψ = (e_θ × e_ζ)/J needs only
θ-derivatives, and duality gives e_ψ·∇ψ = 1 exactly (verified: mean 1.0000000092, pointwise
spread 0.9978–1.0037, so the two routes disagree by up to 0.4%). Projecting e_ψ onto that
constraint is exact and unbiased but gives 7638 → 7495 steps. Capping ∂ν/∂ψ alone — which
cleans C/E/H/K completely — gives 7524. Only cleaning g11 as well reaches 6081. The step count
cares about G via g11, and g11's residual noise sits in the tangential component of e_ψ,
i.e. the θ-parametrization slip — coordinate freedom that no metric identity can constrain. -
Open question worth more than any of the above: what accounts for the other ~80% of the
step growth (§7).
One bookkeeping note: the mpsi ladder (§2–§3) was measured a few merges earlier than the tolerance sweep (§4); the shared
mpsi = 512 / 1e-10 point agrees to 0.3% in nstep between the two, but the tables should not be
mixed element-by-element.Reacted by Min-Gu Yoo- added 5 commits that reference this issue
on Aug 15, 2026 Disposition update: the certified adaptive kinetic grid (#410) is closed without merging — it failed its pre-registered gate (global-scale certificate blind to the dynamically dominant structure; certified raw increments instead of the post-Schur matrices) and convergence evidence across five grids removed its premise. Its role is filled by resonance-aware pinning in #414. The revival path, the bisection map diagnostic (branch
experiment/kinetic-bisection-refine), and the instruction to certify against the torque cross-check rather than eigenvalues are recorded inhandoff/issue376/RESULTS.md§28–§35 and #423.A follow-up to §3 and §7, from #480 and #379. All of the ladders above ran the ForceFreeStates (FFS) forward solve with only
reltolset, soabstolwas the OrdinaryDiffEq default, 1e-6. The error norm is then element-wise: each entry of U is scaled by1e-6 + reltol·|u_i|. Small entries of a large solution column are held to their own relative accuracy, and those are the entries that follow knot-scale roughness in the coefficients (§6). Fortran DCON'sode_stepinstead scalesatolto each column,max|u(:,j)|·tol, so error is measured against the column's size. #480 does the same.DIIID-like_ideal_example, auto grid (573 knots), ForceFreeStates only:Error control Accepted steps Steps per knot interval element-wise (current) 4738 8.3 per-column (#480) 1464 2.6 per-column + knots as d_discontinuities1301 2.3 The least-stable energy is unchanged in all three.
- §3 knot slaving: on the EQUIL - BUGFIX - Fix ideal Δ′ non-convergence on the auto psi-grid #350-era grid, per-column control cuts the steps per knot interval from about 8 to about 2.6. Steps still grow with knots, 1.7× for 2.9× the knots, but more slowly. That is the knot dependence the rest of §7 ("the other ~80 %") is about. Stepping to each knot (
d_discontinuities) removes only another 11 %. - §4 and §5 need re-measuring under per-column control. Under element-wise control, the α floor and "Vern9 never achieves its formal order" came from small entries tracking the coefficient noise. In our tests under per-column control, Vern7 took more steps than Vern9 and gave worse NTV torque (ForceFreeStates - BUGFIX! - Scale the ODE absolute tolerance to each solution column #480). So recommendation 3 (switch to Vern7) doesn't carry over, and the tolerance choice in recommendation 2 should be re-measured.
- Recommendation 4 (consistent θ parametrization across surfaces) is untouched by this. It removes the noise itself rather than how it is weighed, so the two should stack.
Reacted by Nikolas Logan- §3 knot slaving: on the EQUIL - BUGFIX - Fix ideal Δ′ non-convergence on the auto psi-grid #350-era grid, per-column control cuts the steps per knot interval from about 8 to about 2.6. Steps still grow with knots, 1.7× for 2.9× the knots, but more slowly. That is the knot dependence the rest of §7 ("the other ~80 %") is about. Stepping to each knot (
I confirmed the column normalization and theta alignment stack well. @jhalpern30 let's close this when those two are merged
Empirically, we all know larger
mpsislows down the ForceFreeStates.The hand waving reason is the 10s of thousands of spline interpolations done during the ODE integration are a bit slower. But why would this be the case when we are using hints to cleanly provide the knot locations of interest and these are only slowly (and monotonically) incrementing? It seems like there would be effectively no search difference with knot number in this case. Is there some further optimization of the spline evaluation calls that would make this loop insensitive to the spline knot length? It seems like there should be.
@mgyoo86 might have some insights.
Possible topic for David Smitth ( @drsmith48 ) and UW hackathon