KineticForces - FEATURE! - Suppress kinetic terms inside the near-axis validity limit and report validity profiles

[logan-nc]

Stacked on #408 (performance/decoupled-el-matrix-grid). Physics ruling from the issue #376 DIII-D kinetic investigation: do not chase drift-kinetic physics where the model has lost validity near the axis.

What

1. Near-axis kinetic validity suppression (axis_validity_suppression, default on; a single Bool — no tuning parameters). The zero-orbit-width drift-kinetic ordering fails where thermal ion orbit widths reach the local minor radius. The boundary is computed from the equilibrium + kinetic profiles at runtime: ψ_c = outermost ψ where max( potato width (q²ρ²R₀)^⅓, banana width qρ/√ε, poloidal gyroradius qρ/ε ) ≥ ⟨r⟩ (all coefficient 1). A C² quintic envelope zeroes the calculated kinetic increments below ψ_c (kernel evaluation is skipped there), rising to 1 at 2ψ_c; the kernel grid is augmented with knots across the band so coarse decks resolve the envelope. The same boundary and envelope apply to the NTV ψ torque quadrature — one source of truth.

2. KineticForces/Validity/ output group, written whenever kinetic profiles are used: rho_i, rho_banana, rho_theta, w_potato, r_minor, L_p, L_q, d_separatrix, psi_c, envelope, is_valid (true where max orbit width < ⟨r⟩, ρ_banana < L_p and L_q, and max orbit width < distance to separatrix). The far edge and steep-gradient regions are flagged, never suppressed — they can dominate the physical NTV; suppression is reserved for the region where the model both fails and poisons the numerics.

3. Resonance-aware auto grid (was EQUIL - NEW FEATURE - Pin located kinetic-resonance surfaces into the auto psi grid #422, combined here): the two-pass auto grid's criterion is ideal-driven and knows nothing about kinetic resonance locations, so when a run builds calculated kinetic matrices the located Ω_ℓ = 0 surfaces (same kinetic_resonance_psi_nodes locator the NTV quadrature panels at) are pinned into the grid as plain knots via merge_mandatory_nodes — knot-at-node, no cleared zone, inserted before rational bracketing so the Δ′ clean-interval treatment wins where a resonance sits inside a rational's bracket. Nodes below ψ_c are suppressed anyway and not pinned. DIII-D: 6 nodes located, +2 net knots, et[1] unchanged to 2e-6, EL steps −14% (the ψ=0.17 resonance was previously unknown to the grid). Stress-tested with toroidal_rotation_factor = 0.2 (parks resonances at ψ = 0.578, 0.812) against a 1024-knot gold: naive auto grid already agrees to 1.24e-3 and pinning reproduces it to 1e-6 — cheap insurance plus the step win, and the same margin held collisionless (nufac = 0.02).

4. Regularization off in self-consistent kinetic runs.reg_spot smooths the ideal 1/(m−nq) divergence of ξ^ψ′ and ξ^α before they drive the NTV integrand. The self-consistent kinetic Euler–Lagrange operator has no such divergence — Park & Logan (Phys. Plasmas 24, 032505 (2017) §III D) decompose F_k = Q F̄_k Q − P_l†Q − Q P_u + R₁ with R₁ ≠ 0 at Q = 0, and with finite torque det F̄ goes complex, removing the singularity from the solution and the torque integral. Kinetic runs now force reg_spot = 0 and log the override; ideal runs are untouched.

   Verified in this run's own operator: σ_min(F̄) at the rational is 2.6e-17 ideal vs 4.2e-3 kinetic (never below 3.3e-3 across the window) — which is also why ksing_find correctly reports no kinetic singular surfaces. And in the displacements: ideal unregularized |ξ^α| peaks at 643.8 while kinetic unregularized peaks at 0.059, indistinguishable from the regularized kinetic value (0.055).

   configuration	max |ξ^α|	NTV torque [N·m]	EL dissipation [N·m]
   ideal, reg_spot = 0	643.8	6074.3	—
   ideal, reg_spot = 0.05	0.058	0.554	—
   kinetic, reg_spot = 0.05	0.055	0.1655	0.1322
   kinetic, reg_spot = 0	0.059	0.1324	0.1322

   The two independent torque calculations agree to 0.15% with it off and 20% with it on. Full derivation, equations and figures in the review package; a docs section landed in docs/src/kinetic_forces.md.

5. Near-singular F̄ guards. The cond(F̄) scan that locates kinetic singular surfaces (Park & Logan Eq. 70 assembly, so shifted/split resonances are included) already swept 2000 points into SingularSurfaces/Kinetic/scan_cond, but only reported peaks above the 1e8 singular threshold. Sub-threshold peaks are now reported, and where a peak's FWHM contains fewer than three matrix knots the kernel is evaluated at a few targeted ψ and spliced in — existing values reused, so one kernel call per added knot rather than a re-formation — respecting MIN_KNOT_SPACING, the validity band, and a 24-knot cap.

   case	qualifying peaks	knots inside FWHM	inserted	result
   nominal	none	—	0	bit-identical (et[1], 3665 steps, torque unchanged)
   low-collisionality	none	—	0	bit-identical
   collisionless + rot 0.2	ψ=0.5296, 0.5072	0 and 0	6	et[1] Δ3e-6, PE-vs-NTV 3.66% → 3.57%

   One of the two unresolved peaks coincides with a top NTV torque-density peak. A resolved grid inserts nothing, so this is invisible on healthy cases by construction.

Why (measured, DIII-D kinetic-calculated on the production auto grid)

The kinetic increments diverge toward the axis (~ψ^−0.7; rtol-invariant, identical on develop and the stack — real model output, evaluated outside its validity domain, no clamp anywhere in the kernel), and on the stack 94% of 223,271 EL integration steps land at ψ < 0.01 while contributing nothing physical:

	grid / kernel evals	et[1]	EL steps	wall (FFS, 32 threads)
develop c42d558 (ad-hoc auto grid, no adjustments)	554 / 554	1.006441 − 0.259620i	7,972	427 s
stack unsuppressed	288 / 288	1.005545 − 0.259371i	223,271	905 s
stack suppressed (ψ_c = 0.04)	288 / ~258	1.005298 − 0.259472i (Δ 2.5e-4)	4,246	365 s

Honest framing (bisected): the 223k-step explosion is a stack regression introduced by #398 (that branch alone reproduces it bit-for-bit; develop does not explode — see #398's disclosure). This PR cures it on physics grounds and finishes below the develop baseline on every axis: fewer steps (4,246 vs 7,972), less than half the kernel evaluations (~258 vs 554 — develop spends 45 of its evals below ψ_c on invalid physics), and et[1] within 1.1e-3 of develop / 2.5e-4 of the stack's converged unsuppressed reference.

A 53× step collapse (also eliminating multi-GB solution dumps) for a 2.5e-4 eigenvalue change, measured against a doubly-converged reference (kernel tolerance ×10: Δ 2.7e-7; denser equilibrium grid: Δ 1.6e-4) and cross-checked against develop (Δ 1.1e-3).

Fortran precedent

ktanh_flag (dcon/fourfit.F:1117) suppressed the same region with four hand-tuned user knobs (ktc, ktw, kinfac1/2) and no physics setting the location. Adopted as supporting evidence for the decision; rejected as a design — here the boundary is profile-derived with zero user parameters.

Runtimes

run	Force-Free States	KineticForces	total	insertions
nominal, before	206.3 s	36.9 s	275.7 s	—
nominal, with guards	205.1 s	37.6 s	274.5 s	0
collisionless + rot 0.2, before	211.0 s	63.8 s	306.6 s	—
collisionless + rot 0.2, with guards	214.1 s	60.4 s	305.8 s	6

The cond(F̄) scan sits below the run-to-run noise floor (1000 evaluations of four spline lookups plus a 35×35 condition number); each inserted knot costs ~0.5 s at 28 threads, so the 24-knot cap bounds the worst case near +12 s on a ~275 s run.

Regression harness — and an honest coverage gap

regress --cases diiid_n1,solovev_kinetic_calculated,solovev_kinetic_ntv,solovev_kinetic_nuzero --refs performance/decoupled-el-matrix-grid,local:

• diiid_n1 (ideal): 47/47 unchanged — bit-identical.
• Solovev kinetic cases: energies move ≤ 0.06%, step counts ≤ 1.1%, NTV torque 0.12% (its quadrature domain now starts at ψ_c). These are the axis-validity deltas already reported; the regularization change and the F̄ guards add nothing to them, for the reason below.
• DIII-D kinetic-calculated, nominal: bit-identical with all guards active (et[1] = 1.005300−0.259473i, 3665 steps), and a deck that omits reg_spot entirely reproduces the same numbers — verifying the override picks up the struct default, not only explicit settings.

The gap, stated plainly: no current harness case exercises the regularization change. solovev_kinetic_ntv runs kinetic_factor = 0 (not self-consistent), and Solovev_kinetic_calculated_example has no [PerturbedEquilibrium] section at all — its PE stage runs in 0.000 s. The only coverage is the DIII-D scratch runs above. #407 (unmerged) adds [ForcingTerms]/[PerturbedEquilibrium] to that example and would give this change real harness coverage; reviewers of the two together should note the interaction. A DIII-D kinetic-calculated FFS-only harness case remains the outstanding test-infrastructure item (also flagged in #398 and #423).

Rebased onto develop (349a0c2) — and what the port changed

Merged develop on 2026-09-04, adopting its structures rather than working around them. Three of
this branch's changes became deletions, because develop had solved the same problems
independently: _compute_fkg_matrices is now idempotent (the ideal-restore workaround is gone),
KineticForces/Output.jl already opens-or-creates its group, and the :static kinetic-threading
fix landed in 9fe8a60c9. What remains is ported onto MatrixSplines, with the arbitrary-ψ kernel
argument type-stable (psis::Vector, empty meaning metric.xs) and the orchestration layer gaining
two lines — the logic lives in KineticForces.resonance_grid_nodes / axis_validity_boundary /
write_validity! and PerturbedEquilibrium.kinetic_regularization_kwargs.

Multi-ion ψ_c (recorded for the historical record)

develop's multi-ion NTV means ψ_c is now the widest-orbit species' boundary. Measured per
species on the DIII-D-like H-mode equilibrium:

species	ψ_c	ratio to main ion
deuterium (main, z=1 m=2)	0.040	1.00
carbon impurity (z=6 m=12)	0.010	0.25
hydrogen (z=1 m=1)	0.030	0.75
helium (z=2 m=4)	0.030	0.75
electrons	0.000	—

Because ρ ∝ √(mT)/(Z·e·B₀), impurity orbits are narrower than the main ion's — carbon's ψ_c is
4× smaller — so adding impurities cannot widen the suppressed region on this case; the main ion
sets it. A species with larger √m/Z does move it: the Solovev D-T deck's tritium pushes ψ_c from
0.10 to 0.12, which is what exposed the bug below.

Bug found and fixed during the port

With tritium setting ψ_c = 0.12, the envelope's transition band landed on that equilibrium's
rational surface at ψ = 0.1221. Multiplying near-singular kinetic increments by a rapidly varying,
near-zero envelope is not representable by the matrix splines, and the overshoot propagated NaNs
into the stability solve — solovev_kinetic_calculated failed while develop passed. The boundary is
now moved clear of any rational whose window the band would cut through (ψ_c → 0.127 there), on the
physical grounds that where orbit widths already reach ⟨r⟩ the resonance is not trustworthy anyway,
so suppressing it wholly is more honest than half-suppressing it. The boundary is computed once
per run and threaded to the kernel, the band knots and the Validity output so the three cannot
disagree.

Verification after the port

• Unit: runtests_kinetic 277/277, runtests_multiion 52/52, runtests_sing 76/76,
  runtests_h5_schema 22/22.
• Harness vs develop: all four cases run; solovev_kinetic_calculated Re(et[1]) 0.94%,
  Im(et[1]) 6.7%, steps 737 → 792; solovev_kinetic_ntv and _nuzero similar magnitude. These are
  the deliberate physics of the suppression and the regularization change.
• DIII-D headline set re-measured on the ported stack (scratch decks, measured before develop's KineticForces - BUGFIX! - Drop the spurious major radius from the omega_D prefactor #391/KineticForces - BUGFIX! - Give the energy integration its own tolerances #392 kinetic bugfixes landed on 2026-09-04): nominal 3665 steps, PE 0.1322 / NTV 0.1324
  (0.2%); suppression off 222,489 steps; collisionless + slow rotation 3031 steps — all
  reproducing the pre-port values.

Release note

• Audience: users
• Numerical impact: kinetic runs move — near-axis kinetic terms are suppressed inside the orbit-width validity limit and reg_spot is forced to 0 in self-consistent kinetic runs (harness @ 14254f5)
• Migration: kinetic results move, so re-pin local baselines; set axis_validity_suppression = false in [KineticForces] to recover the previous behaviour for debugging

Calculated-kinetic runs no longer evaluate the drift-kinetic model where its zero-orbit-width
ordering has failed: the code computes the validity boundary from the profiles themselves and
smoothly suppresses the kinetic terms inside it, which both removes unphysical near-axis structure
and cuts the Euler–Lagrange work that chasing it cost. Every run that uses KineticForces now writes
an orbit-width and gradient-length Validity profile so the trustworthy region is visible rather
than assumed, and self-consistent kinetic runs no longer regularize the displacement — the kinetic
singularity is shifted and split, not divergent, so the ideal-only regularization was inconsistent.

Regression report

regress --cases diiid_n1,solovev_n1,solovev_kinetic_calculated,solovev_kinetic_ntv,solovev_kinetic_nuzero --refs develop,local,
baseline develop @ 349a0c2, working tree at 14254f5. Deltas below bundle this PR with
its base #398 (this branch is stacked on it); #398's own table is in that PR.

diiid_n1 develop local diff
total energy Re(et[1]) 8.012318e-01 8.004564e-01 0.10%
PE toroidal torque 5.087465e-02 5.086417e-02 0.02%
NTV torque FGAR [N.m] 5.296762e-01 5.285453e-01 OK (within tolerance)
ODE steps (total) 4572 1974 -56.8%
singular surfaces / psi / q unchanged; q0, q95, beta_t, beta_n 0.00%
solovev_kinetic_calculated develop local diff
total energy Re(et[1]) 1.873745e+00 1.828561e+00 2.41%
total energy Im(et[1]) -1.434288e+00 -1.325305e+00 7.60%
ODE steps (total) 734 789 +7.5%
solovev_kinetic_nuzero develop local diff
total energy Re(et[1]) 2.243336e+00 2.290643e+00 2.11%
total energy Im(et[1]) -2.026025e+00 -1.561707e+00 22.92%
solovev_kinetic_ntv develop local diff
NTV torque fgar (see #398) (see #398) dominated by #398
NTV psi quadrature evals 840 60 -92.9%

Attribution — almost all of the Solovev kinetic movement belongs to the base, not to this PR.
Running #398 alone against the same develop gives calculated Re 1.829090 / Im −1.322557 and
nuzero Re 2.290779 / Im −1.557856; adding this PR moves them to Re 1.828561 / Im −1.325305 and
Re 2.290643 / Im −1.561707. So this PR's own increment on Solovev is 0.03% / 0.21% (calculated)
and 0.006% / 0.25% (nuzero), with ODE steps 759 → 789 and 754 → 794. solovev_kinetic_ntv's
torque move is likewise inherited from #398. That the increment is this small is the expected
result, not a disappointment: Solovev's ψ_c is small and its et[1] barely sees the suppressed
region, matching the 2.5e-4 et[1] change measured on DIII-D. The physics this PR is really being
judged on — the 61× step reduction and the PE-vs-NTV torque agreement — lives on the DIII-D decks
in the review package, which the harness cannot yet reach (see the caveat below).

Both sides were measured against a develop that includes its own kinetic bugfixes (#391 ω_D
prefactor, #392 energy-integration tolerances, both merged 2026-09-04).

Caveat on coverage: the Solovev calculated example has no [PerturbedEquilibrium] section and
solovev_kinetic_ntv runs at kinetic_factor = 0, so the harness cannot yet exercise the
regularization change; #407 adds the missing sections and would give it real coverage. The
DIII-D evidence for it is in the review package.

⚠️REQUIRES THIRD-PARTY HUMAN REVIEW BEFORE MERGING — NON-NEGOTIABLE. DO NOT MERGE WITHOUT AN APPROVING HUMAN REVIEW.

🤖 Generated with Claude Code

https://claude.ai/code/session_01LzbLFQKyuRE5DYZmLokKmk

[logan-nc]

Reviewer package (figures: validity-scale crossings defining ψ_c, before/after step distribution, flow diagram, verification record): https://claude.ai/code/artifact/09745676-ff9b-49c6-a14b-74f093427d00 — flip to shared for reviewers as with #398/#408.

[logan-nc]

Review package updated with the regularization analysis (same link): https://claude.ai/code/artifact/09745676-ff9b-49c6-a14b-74f093427d00 — adds the operator σ_min evidence, the ideal-vs-kinetic ξ^α overplot, and the torque table. Verified on the shipped DIII-D deck (which sets reg_spot = 0.05): the override logs, PE/fgar agree to 0.2%, and the EL solve is bit-identical (et[1] and 3665 steps unchanged) — reg_spot only ever entered the post-processing that builds the NTV input.

[logan-nc]

Note for anyone benchmarking from this branch before it is reconciled with develop: this stack forked at 9491f893d (2026-08-16), which predates develop's :static kinetic-threading fix (9fe8a60c9). Runs launched from feature/kinetic-axis-validity as it stands are therefore nondeterministic at rare-task-migration frequency — identical decks were measured at et[1] = 1.0053 / 1.1355 / 0.9895. Reconciliation with develop brings the fix; until then, treat single runs off this branch as indicative rather than exact. All numbers reported in this PR were taken on a stack that carried the fix locally.

[github-actions]

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.
docs/development/contributors.md suggests lead developers to ask.
Merging is not blocked here, but no pull request may be merged without human review.

[logan-nc]

Review package (self-contained, private artifact): https://claude.ai/code/artifact/09745676-ff9b-49c6-a14b-74f093427d00

Covers the validity boundary and where it comes from, what the invalid region was costing, the resonance-aware grid, the near-singular F̄ guards, and the regularization-off argument — plus, new since the port onto develop, the per-species ψ_c table (deuterium sets it; carbon's is 4× smaller) and the rational-inside-the-transition-band bug the port exposed and how it is now guarded.