EQUIL/FFS - IMPROVEMENT - Consistent SFL surface sampling and a core knot-density cap for the EL coefficient splines - #398
Conversation
…fieldline angles Each surface was traced independently and then splined on that surface's OWN solver-chosen abscissae before being resampled onto the common theta grid, so the resample error was uncorrelated between neighbouring surfaces -- white noise in psi that grid refinement amplifies rather than reduces. Neither reltol nor abstol touched it, because it is remap interpolation error rather than integration error. The trace now returns its dense solution, and equilibrium_solver root-solves (Brent, bracketed by the monotone jac-weighted flux integral) for the angle at which the normalised straight-fieldline angle reaches each target node, evaluating there. Every surface is sampled at identical abscissae and the resample error at the output nodes is zero. The arclength tracer returns nothing for the solution and keeps the previous path. Measured on DIII-D stripped decks at eulerlagrange_tolerance 1e-10, accepted Euler-Lagrange steps fall 2309/3977/7638 -> 1881/2768/4403 for mpsi 256/512/1024, a 42% reduction at mpsi=1024, and the per-doubling growth drops from 1.92x to 1.59x. Surface geometry residuals now converge with refinement instead of sitting on a floor, and every knot-to-knot correlation turns positive. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
…ed as PR #398 The A3 kill-switch fires: route (a) geometry already agrees with a near-exact trace to 5e-10..2e-8, i.e. at the integration tolerance, so the planned SFL reparametrisation would fix an error that is not there. Also records that nstep is hypersensitive -- a 5.3e-11 geometry perturbation moves it 1.2% -- so the leftover 15-20% gaps are not reliable signal; that A1 places the residual in the traced construction rather than the EFIT input (analytic input, traced 1.47x vs inversion 1.18x); and that tightening the trace tolerance 1000x on that case is null. Notes the harness baseline was refreshed from a stale local develop, and that the first etol test was a no-op because the deck has no etol key. No src changes on this branch. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
The original PR evidence used forward-integrator cases, which emit no BVP Delta-prime at all, so route (a) had never been checked against the observable that section 19 showed can break silently. Ran diiid_n1_riccati and gal_resistive_diiid. Tearing-consumed quantities are unchanged: PEST3 Delta diagonal 0.01%, Delta-prime matrix norm 0.03%, inner-layer Delta 0.00%. The headline 16.73% on the raw BVP diagonal is a single element (q=5) that runs 1e5 -> -1319 -> -2369 and changes sign across the mpsi ladder on BOTH versions, so it cannot discriminate between them. Records that route (a) does not improve Delta-prime convergence either, flagging the unconverged q=5 element as a pre-existing issue worth its own investigation. Wall time improves on all three cases. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
…plines The Euler-Lagrange coefficient splines inherited every knot of the equilibrium grid. A cubic spline's third-derivative jumps at knots scale as (node error)/dpsi^3, so the equilibrium's near-axis packing (dpsi ~ 1e-6) amplifies even tolerance-level node error into huge C2 kinks, and the adaptive integrator's step size becomes slaved to the knot spacing -- measured directly: core jump magnitudes grow ~34x per mpsi doubling while a (tol/J)^(1/4) step model reproduces the observed step-count ladder. The coefficients are near-cylindrical in the core and do not need that packing. Build their splines on a subset of the equilibrium grid with core density capped at dpsi >= 0.05*psi below psi = 0.1; node values are unchanged, only knot density. Measured on DIII-D stripped decks (route-a base, mpsi 512/1024): accepted EL steps 2768 -> 2188 and 4403 -> 3034, per-doubling growth 1.59x -> 1.39x, warm run -19% at mpsi=512, with et[1] unchanged to 3e-8 relative and the Riccati BVP Delta-prime diagonal unchanged to 4e-7 on all five surfaces. Interim fixed-cap form; a matrix-curvature knot selection is planned to replace the fixed rule, with this commit as the fallback. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
…es knots On the production two-pass auto grid the cap is a no-op (its core spacing already satisfies the density rule), so the unconditional message was noise. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
|
Report with visuals for workflow and results can be found here: https://claude.ai/code/artifact/873c8376-3f75-4fe4-934b-4d39f23e2376 |
… protect rationals Reframes the cap as log-uniform sampling of the Frobenius region: near the axis every component is a power law in psi, and cubic interpolation of psi^p on a log-uniform grid errs by ~(p*c)^4/384, so c = 0.05 resolves even the steepest spectrum component (p = mmax/2) to ~2e-4 while the physics responds far below that. Two generalization guards, motivated by cross-equilibrium testing: the capped region now ends at the innermost rational surface when that sits inside psi = 0.1, and no knot inside a rational's RATIONAL_RES_RADIUS window is ever removed -- preserving the Delta'-stencil structure for decks (e.g. higher n) whose rationals reach the core. Both guards are no-ops on every current case, verified: DIII-D m512 and Solovev m512 reproduce the previous cap's step counts and knot sets exactly. Cross-equilibrium check of the rule itself: tj_analytic_direct m1024 (analytic, traced) 1470 -> 959 steps at et[1] 1.3e-9; LAR m1024 (inversion path, clean geometry) 951 -> 875 at et[1] identical to 8 digits; Solovev ldp m512/m1024 ~unchanged steps at ~4.5e-7 absolute et[1] shift (a +-10.4 cancellation amplifies this to 3e-5 relative). Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
|
Review package (visual companion to this PR — mechanism diagram, the Solovev et[1] grid-convergence before/after, step ladders, geometry-residual evidence, and the Δ′ referee table): 📦 https://claude.ai/code/artifact/873c8376-3f75-4fe4-934b-4d39f23e2376 Self-contained page; complements rather than duplicates the description and diff. If the link does not resolve, ask @logan-nc to enable sharing on it. |
… vs #398; deltas inherited Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
|
This pull request is missing an assignee and a reviewer. If you are not ready to name them, mark this pull request as a draft. |
…III-D step explosion is a #398 regression Co-Authored-By: Claude Fable 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01LzbLFQKyuRE5DYZmLokKmk
Port the kinetic guards onto develop's structures rather than fighting the reorg:
- MatrixSplines replaces FourFitVars: kinetic matrices read mats.ideal and return
a fresh MatrixSplines, so the ideal-restore workaround is gone (develop's
_compute_fkg_matrices is idempotent by construction).
- Arbitrary-psi kernel evaluation returns as a type-stable psis::Vector{Float64}
(empty = metric.xs), no Union sentinel.
- Multi-ion: the near-axis boundary is now the widest-orbit species' psi_c.
- Sub-threshold cond(F-bar) reporting follows find_kinetic_singular_surfaces!
into Surfaces/Finding.jl.
- main keeps develop's shape: resonance pinning, the validity boundary, the
regularization policy and the Validity output move into named helpers in the
owning modules (KineticForces.resonance_grid_nodes / axis_validity_boundary /
write_validity!, and kinetic_regularization_kwargs beside perturbed_equilibrium).
Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01LzbLFQKyuRE5DYZmLokKmk
Summary
Each flux surface was traced independently and then splined on that surface's own solver-chosen abscissae before being resampled onto the common θ grid (
DirectEquilibrium.jl). The resample error was therefore uncorrelated between neighbouring surfaces — white noise in ψ that grid refinement amplifies rather than reduces. Neitherreltolnorabstoltouched it, because it is remap interpolation error, not integration error.The trace now returns its dense solution, and
equilibrium_solverroot-solves (Brent, bracketed by the monotone jac-weighted flux integral) for the angle at which the normalised straight-fieldline angle reaches each target node, evaluating there. Every surface is sampled at identical abscissae and the resample error at the output nodes is zero. The arclength tracer returnsnothingand keeps the previous path.Investigated under issue #376.
The headline: the answer stops depending on the grid
Solovev, same commit, only this change differing — the free-boundary energy was not converging before and is now converged to 6 significant figures:
et[1]developet[1]this PRdevelop's
et[1]drifts 4.5× across the ladder and is still moving; the plasma energy drifts with it (−10.4049 → −10.3144). With this change both are grid-invariant, and the step count is nearly flat (1.04×, 1.11× per doubling versus 1.45×, 1.73×).That is the property we want: once the equilibrium splines resolve the equilibrium, adding knots should change neither the answer nor the work.
Regression harness
regress --cases diiid_n1,solovev_n1 --refs develop,local, baseline develop @ 9491f89(branch is current with origin/develop; 1 ahead, 0 behind).
diiid_n1 — physics moves only in the 3rd–4th digit, cost halves:
q0,q95,beta_t,beta_nand the singular-surface locations/count are unchanged to 0.00%.solovev_n1 —
et[1]moves 98%. This is the fix working, not a regression: as the table above shows, the develop value is grid-dependent and non-convergent. Solovev's ν node data sat at ~19% relative white noise under the old resample, so it is the case this change corrects most.et[1]is also a near-cancellation (ep ≈ −10.4, ev ≈ +10.4), so a 7% plasma-energy correction dominates it.On
solovev_n1,q0,q95, the singular-surface count/locations/q-values andmpertare allexactly unchanged;
beta_t/beta_nmove 0.02%/0.08%.Runtime 210.6s → 184.5s.
Kinetic cases — large moves, all of which enter at this PR
regress --cases solovev_kinetic_calculated,solovev_kinetic_ntv,solovev_kinetic_nuzero --refs develop,local. This PR moves the kinetic harness values, some of them dramatically — that is owned here, not hidden:Attribution: the stacked follow-up branches (EL grid cap + certified kinetic grid, knob off) reproduce these local values bit-identically (34/34 tracked quantities unchanged vs this branch's head), so every kinetic delta vs develop enters at this commit — none at the follow-ups.
Interpretation: the old per-surface resample injected ~1e-6-level white-in-ψ geometry error (measured; see the mechanism section), and resonance-dominated kinetic quantities amplify exactly that kind of noise — the Solovev cases are the most sensitive in the suite. The quadrature-cost drop (840 → 60 evaluations for the same tolerance) is direct evidence the develop-side torque integrand carried noise structure the quadrature was chasing. The eigenvalue moves (2–8%, 23% for Im at ν→0) are consistent with the ideal-case finding that develop's Solovev values were grid-dependent while this branch's are grid-convergent.
Second change, folded in from #408: cap the core knot density of the EL coefficient splines
(#408 was folded here on 2026-09-04 — both are grid-side fixes aimed at the same symptom, the cap is a no-op on production auto grids so its measurements only mean anything alongside the step numbers above, and the two were benchmarked as a pair throughout. One file,
src/ForceFreeStates/Fourfit.jl.)Summary
The Euler-Lagrange coefficient splines (
fmats/kmats/gmatsand the primitives) inherited every knot of the equilibrium grid. A cubic spline's third-derivative jumps at knots scale as (node error)/Δψ³, so the equilibrium's near-axis packing (Δψ ~ 1e-6 at high mpsi) amplifies even tolerance-level node error (~1e-9, measured) into huge C² kinks — and the adaptive integrator's step size becomes slaved to the knot spacing. Measured directly: core jump magnitudes grow ~34× per mpsi doubling, and an (tol/J)^¼ step model reproduces the observed step-count ladder.The coefficients are near-cylindrical in the core and do not need that packing. This PR builds their splines on a subset of the equilibrium grid with core density capped at Δψ ≥ 0.05·ψ below ψ = 0.1. Node values are unchanged — only interpolation density — which is why the physics moves at the 1e-7–1e-8 level.
Stacked on #398 (route (a), same investigation); diff shows only the cap once #398 lands.
Measured (DIII-D stripped decks, route-(a) base, mpsi 512/1024)
Not a tuned hack: the rule's basis and its generalization
Form: near the axis every component is a Frobenius power law in ψ; power laws are scale-free, so log-uniform sampling (Δψ ≥ c·ψ) resolves them at constant relative accuracy. Constant: cubic interpolation of ψ^p on a log-uniform grid errs by ~(p·c)⁴/384, so c = 0.05 resolves even the steepest spectrum component (p = m_max/2 = 11) to ~2e-4 — and the physics responds far below that because the steep components carry vanishing solution amplitude. Region: ends at min(0.1, innermost rational −
RATIONAL_RES_RADIUS), and no knot inside a rational's resolution window is ever removed — guards for decks (e.g. higher n) whose rationals reach the core; verified no-ops on every current case.Cross-equilibrium check (four equilibria, two construction paths, three grid families):
The LAR row is the anti-over-fit witness: clean geometry with no noise to exploit, and the cap is still harmless.
Two properties reviewers should know
diiid_n1, solovev_n1, diiid_n1_riccati, gal_resistive_diiid) reproduces EQUIL/FFS - IMPROVEMENT - Consistent SFL surface sampling and a core knot-density cap for the EL coefficient splines #398's numbers exactly. The cap engages only on explicitly packed fine grids (largempsi), which is precisely the issue Performance: Why do large equilibrium splines slow down the code? #376 scenario.handoff/issue376/RESULTS.md§22–§23 (experiment branch).Verification
runtests_eulerlagrange,runtests_riccati,runtests_sing— all pass.handoff/issue376/RESULTS.md§22.Rebased onto develop (349a0c2)
Merged develop on 2026-09-04. The only conflict was
Fourfit.jl, where develop's MatrixSplinesrefactor now builds an immutable
IdealMatrices: the cap logic moved into a named helper,core_capped_knots(xs, rationals), and the constructor builds on the capped subset. Capbehaviour is unchanged — only its plumbing follows develop.
Verification after the port:
runtests_equil279/279,runtests_grid_refinement59/59.Harness vs develop reproduces the headline effect —
diiid_n1ODE steps 4572 → 1974 (−56.8%),with energies moving 0.02–0.28% (Re(et[1]) 0.10%) and q0/q95/β/singular-surface locations at 0.00%.
Bisected cleanly (same deck, worktrees per commit): the DIII-D
kinetic_source="calculated"forward integration takes 223,271 steps on this branch vs 7,972 on develop (FFS wall 1079 s vs 427 s), with physics agreeing to ~0.1% (et[1] 1.005545 vs 1.006441). This PR alone reproduces the full-stack number bit-for-bit; the stacked follow-ups add nothing. The proximate signature is an integrator grind at ψ ≈ 0.0012–0.0015 (median step 9e-9 vs develop's 2e-6); the deep mechanism is unresolved at the h5 level — the kinetic matrices, their spline third-derivative jumps, condition numbers,crit, integration start point, and solution norms are all identical at every stored digit between the two runs (open follow-up: the unstored F-matrix spline / step-control internals). No harness case covers this configuration — a DIII-D kinetic-calculated FFS-only case should be added. The stacked PR #414 (near-axis validity suppression) removes the sensitive region on physics grounds and cures the regression outright (4,246 steps — below develop's own baseline).Caveat, stated plainly: a three-orders-of-magnitude move in the ntv-case torque means the develop baseline for that tracked value was noise-dominated, and the new value has no independent reference yet. This needs a physics reviewer's judgment, not just the attribution argument. If this PR is accepted, the kinetic harness baselines must be re-pinned on the merge commit.
Tests
test/runtests_equil.jl— 279 pass (includes the arclength path; its return-type annotation was updated to match)test/runtests_grid_refinement.jl— passtest/runtests_tj_analytic.jl— pass (exercisestj_analytic_direct, the same code path)Scope / follow-ups
eq_type ∈ {efit, efit_arclength, imas, sol, tj_analytic_direct}— everything reachingequilibrium_solverinDirectEquilibrium.jl.InverseEquilibrium.jl(chease / lar / tj_analytic / efit_by_inversion), which does its own two-stage SFL resample. Follow-up.🤖 Generated with Claude Code
Follow-up verification: Δ′ and wall time (requested in review)
The original evidence used
diiid_n1/solovev_n1, which runintegrator = "forward"and thereforedo not emit the BVP Δ′ matrix at all. Re-ran the Δ′-tracking cases.
diiid_n1_riccatigal_resistive_diiid— the quantities the tearing/matching path consumesReading the 16.73%
It is one element. Δ′ BVP diagonal across an mpsi ladder on the riccati deck:
The 4th entry (q = 5 surface) runs 1e5 → −1319 → −2369 and changes sign: it is not converged in
mpsi on either version, so it cannot discriminate between them. The other four diagonal entries
agree between versions at every grid, and version-to-version agreement on the full diagonal is
0.25% / 2.0% / 0.03% at mpsi 256 / 512 / 1024.
Caveat worth stating plainly: this PR does not improve Δ′ convergence either (drift
174.7%/83.1% versus develop's 174.0%/79.5%). The unconverged q = 5 Δ′ element is a pre-existing
issue that this change neither causes nor fixes, and it deserves separate attention.
Wall time improves on every case measured: forward
diiid_n1210.6 → 184.5 s,diiid_n1_riccati202.5 → 192.8 s,
gal_resistive_diiid216.8 → 213.6 s.