KineticForces - BUGFIX! - Fit the parallel-velocity spline periodically

[logan-nc]

Release note

• Audience: users
• Numerical impact: Kinetic eigenvalues shift sub-ppm (~3.5e-6 absolute on O(1) values, 0.00% relative) and the ODE step count moves by one; diiid_n1 and solovev_kinetic_ntv are unchanged (harness @ 228a8e9)
• Migration: Re-baseline only if you compare results bitwise; no config or API change.

The parallel-velocity spline used to locate bounce points is now fitted periodically, matching the periodicity of the quantity it samples. The previous non-periodic fit was only C⁰ at the θ=0/1 seam and could manufacture false near-seam extrema and bounce roots; upstream this produced non-finite bounce/precession frequencies and an integrator abort on an ITER case.

---

What

Fit the parallel-velocity (bounce-point) spline with periodic rather than non-periodic endpoint conditions.

Julia analog of Fortran GPEC PR #284, which made the same spline_fit(vspl, "extrap") → "periodic" correction in pentrc/torque.F90.

Why

The spline samples 1 − (λ/B₀)·B(θ), which is periodic on the closed poloidal interval. The adjacent equilibrium field spline (tspl) was already periodic, but this derived one was not — so the fit is only C⁰ at the θ=0/1 seam and can manufacture false near-seam extrema and bounce-root pairs.

On a synthetic shaped field (B = 1 + 0.30cos2πθ + 0.04cos6πθ, global max at the seam), the one-sided derivative mismatch across the seam drops ~5×:

fit	seam derivative mismatch
non-periodic endpoint	-1.376e-04
periodic	2.606e-05

Upstream, the same defect produced non-finite omega_b/omega_D and an LSODE ISTATE=-5 failure on an ITER TC-24 case.

Fixed at both construction sites — the tpsi! ψ-quadrature path (Torque.jl) and the _setup_surface_state kinetic-matrix path — since they build the same object for the same purpose.

Naming

B_extrap → B_vpar and _vpar_from_extrap → _vpar_from_spline: after this change the old names describe an extrapolating fit that no longer exists. Mechanical rename, 23 references, confined to Torque.jl and BounceAveraging.jl.

Scope note vs. PR #341

#341 (from a fork) proposes this same periodic fix plus a ~90-line rewrite of bounce-root enumeration that splits each cubic knot interval at its analytic stationary points and reaches into spline internals (.cache.x, .y, .z).

This PR deliberately carries only the periodic fit, because that is the root cause: #341's own rationale for the root rewrite is that the non-periodic fit "can create false near-seam extrema and root pairs" — which this fix removes at source. I could not reproduce a case where the current Roots.find_zeros misses roots once the fit is periodic (identical root counts across λ ∈ [0.80, 1.28] on the two-well test field above).

Per the repo's minimal-change discipline, the root-enumeration rewrite should stand or fall on its own evidence — a case where find_zeros demonstrably misses a root with a periodic fit. Happy to take it as a follow-up if such a case exists.

Validation

• runtests_kinetic.jl: 277/277. Note this is not sufficient on its own — those unit tests do not exercise the pitch-integration / compute_bounce_data paths; the kinetic regression cases do.
• Regression harness: pending — draft until solovev_kinetic_ntv, solovev_kinetic_calculated, solovev_kinetic_nuzero are run against develop. Expect small movements in the bounce-averaged quantities (the seam handling changes), which is the point of the fix; I'll post the table before marking ready.

🤖 Generated with Claude Code

[logan-nc]

Ready for review; not yet ready to merge.

The physics/code argument is complete and self-contained, so review can start now. The regression-harness table is still outstanding and is a merge gate per CLAUDE.md — I'll post it here before this is merged.

Note the machine currently has other regression runs in flight, so the kinetic cases are queued behind the PR #339 re-baseline rather than run in parallel.

[logan-nc]

Regression harness — develop (b3abe074) vs this branch

regress --cases diiid_n1,solovev_kinetic_ntv,solovev_kinetic_nuzero,solovev_kinetic_calculated --refs develop,bugfix/kf-periodic-vpar-spline

Case	Result
diiid_n1	48 unchanged (NTV torque FGAR and dW FGAR bit-identical)
solovev_kinetic_ntv	6 unchanged (ψ-quadrature NTV bit-identical)
solovev_kinetic_nuzero	6 changed, all at 0.00%
solovev_kinetic_calculated	6 changed, all at 0.00%

Every movement is sub-ppm — the largest relative shift in a physical quantity is ~4e-7:

solovev_kinetic_calculated develop branch Diff
total energy Re(et[1]) 1.845905e+00 1.845909e+00 3.310e-06 (0.00%)
total energy Im(et[1]) -1.421059e+00 -1.421063e+00 3.554e-06 (0.00%)
plasma energy Re(ep[1]) -8.593964e+00 -8.593958e+00 5.730e-06 (0.00%)
ODE steps (total) 1800 1803 3 (0.17%)
solovev_kinetic_nuzero
total energy Re(et[1]) 2.245439e+00 2.245439e+00 4.455e-07 (0.00%)
total energy Im(et[1]) -2.034869e+00 -2.034866e+00 2.558e-06 (0.00%)
ODE steps (total) 1837 1836 1 (0.05%)

Reading this

The near-zero movement is the expected and desired outcome, and it is consistent with what I found while evaluating the scope question above: on well-resolved, well-behaved equilibria the two fits give the same bounce roots (identical root counts across λ ∈ [0.80, 1.28] on a synthetic two-well field), so the seam correction only perturbs the last few digits.

The justification for the fix is therefore not a change in these numbers — it is that the sampled quantity is periodic and the fit must be too. The failure mode it removes is the pathological one: upstream, the non-periodic seam produced false near-seam extrema, non-finite omega_b/omega_D, and an LSODE ISTATE=-5 abort on an ITER TC-24 case. This table shows the correction is safe on the cases we track; it does not, and cannot, demonstrate the bug it prevents.

The two ψ-quadrature NTV cases being bit-identical is the useful negative control: the torque path is untouched on these equilibria.

[logan-nc]

Correction to the table above

I described solovev_kinetic_ntv and the diiid_n1 NTV rows as "bit-identical". That overstates what the harness shows — those quantities are guarded by absolutenoise_thresholds that are coarse relative to their magnitudes:

Quantity	value	noise_threshold	effective relative sensitivity
solovev_kinetic_ntv NTV torque fgar	~1.3e-4 N·m	1e-8	~1e-4
diiid_n1 NTV torque FGAR	5.59e-1 N·m	1e-2	~1.8%
diiid_n1 NTV kinetic energy dW FGAR	6.9e-2 J	1e-3	~1.4%

So "unchanged" for those rows means "changed by less than the noise floor", not "identical". Given the sub-ppm movements this fix produces in the quantities that are finely guarded (the kinetic eigenvalues, at 1e-10), a same-sized change in the NTV torque would sit far below its threshold and report as unchanged either way.

The conclusion in the comment above is unaffected — the movements are sub-ppm wherever we can resolve them — but the NTV rows are not the tight negative control I called them.

[logan-nc]

Regression re-run after re-syncing with develop

PR #339 merged, so this branch was updated to develop@c42d558e and the harness re-run at the new head. Baseline and branch now both contain the multi-ion work.

diiid_n1	Summary: 47 unchanged
solovev_kinetic_ntv	Summary: 6 unchanged
solovev_kinetic_nuzero	Summary: 7 changed, 7 unchanged
solovev_kinetic_calculated	Summary: 7 changed, 7 unchanged

Every changed row on the two Solovev kinetic cases is at 0.00% relative (absolute diffs ~1e-6 to 1e-8 on O(1) eigenvalues, plus one-step ODE count moves) — they cross the harness's tight 1e-10 eigenvalue threshold but are not physically meaningful.

[logan-nc]

Note on the harness stamp. The branch was updated to develop@4e430f170 after the regression run, so the _(harness @ …)_ stamp now predates the head and the conventions check will flag it as stale.

That delta is #424, and it cannot change any tracked quantity. Its entire src/ footprint is three PerturbedEquilibrium files, and the only non-comment line in the whole diff is an HDF5 long_nameattribute string:

- "Energies/toroidal_torque" => (; long_name="net toroidal torque on the plasma", units="N*m")
+ "Energies/toroidal_torque" => (;
+ long_name="boundary-response toroidal torque ...", units="N*m")

Everything else is docstrings and comments. Re-running the harness would reproduce the posted table exactly, so the numbers above stand. Title and release note is not a required check, so this does not block the merge; Documentation and Tests are the gates and are running on the updated head.