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. (This was claimed but not true until the review round: the quadrature recomputed ψ_c per species and skipped the rational clearing. See Review round below.)

2. KineticForces/Validity/ output group, written whenever kinetic profiles are used (from the ForceFreeStates output stage — it was previously gated behind an NTV stage that calculated-kinetic runs never reach, so FFS-only runs got no group at all): 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 kinetic runs (harness @ 2abea51)
• 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.

Review round (2026-09-19) — what changed since @ebursch's review

Reconciled onto the base branch and develop (#403), then addressed the review plus an
independent pass (fortran-physics-reviewer, clean-code-reviewer, and a manual read).

Reviewer's three comments — all addressed. GridRefinement.jl's pinned-knot block and the
kinetic_forces.md regularization section are trimmed. The third site
(KineticForcesStructs.jl tolerances) turned out not to be this branch's code — it arrived on
develop with #392 and only appears in this diff because the merge base makes develop's commits
look like additions. Trimmed anyway.

Six substantive defects the review did not reach (2dcbbf170, 2abea51d2):

#	defect	consequence
1	ψ_c computed three different ways	Validity/psi_c in gpec.h5 described a different domain than the torque was integrated over; the standalone torque kept near-axis contributions the EL build had already zeroed — the exact cross-check the PE-vs-NTV agreement rests on
2	clear_rational_windows cascaded unbounded	the scan window widens as ψ_c grows; rationals at 0.15/0.25/0.45 walk ψ_c from 0.10 to 0.455
3	boundary scan had no contiguity requirement	a non-monotonic q (edge rise, reversed shear) re-triggers the criterion far from the axis and drags the band outward
4	Validity never written on FFS-only kinetic runs	the group was missing on exactly the deck this PR flags as uncovered
5	Validity datasets carried no long_name/units	its table sat in MAIN_H5_ANNOTATIONS, applied at the end of write_outputs_to_HDF5; write_validity! runs after that, so the entries were dead
6	zero unit coverage for ~230 lines of boundary code	—

Fixes: axis_validity_boundary is now the only entry point, computed once per run and threaded
into the quadrature; the clearing is capped at VALIDITY_CLEAR_MAX_FACTOR× the orbit-width
boundary with a warning naming the two ways out; the scan stops at the first valid surface;
write_validity! moved to the ForceFreeStates output stage and annotates its own group (matching
HDF5Schema's stated rule that sub-writers keep their tables next to their writers); and
test/runtests_kinetic_validity.jl covers the envelope's C² join, the cascade cap, the contiguity
guard and the reg_spot override (27 tests).

Not changed: scoping reg_spot = 0 to kinetic_source == "calculated" was proposed and
rejected — the fixed kinetic terms also remove the resonant singularity, so both sources
correctly force it to 0.

Regression report

regress --cases diiid_n1,solovev_n1,solovev_kinetic_calculated,solovev_kinetic_ntv,solovev_kinetic_nuzero,solovev_kinetic_multiion --refs develop,<pre-fix>,<branch>,
baseline develop @ f36ecef, branch at 2abea51 (HEAD). Run with three refs so this
round's fixes are separable from the PR as a whole; <pre-fix> is cf7a40cc4, the reconciled tree
before the review fixes. Deltas vs develop still bundle this PR with its base #398.

The three-ref run was taken at 0f8dd766f; it was then repeated at HEAD (2abea51d2, the HDF5
metadata fix) against develop, and every tracked value is identical — that commit changes
attribute strings and a warning message only, with no numerical path.

This round's fixes moved only the NTV quadrature domain, and nothing else:

case                        pre-fix -> branch
diiid_n1 (ideal)            48 unchanged  -- bit-identical
solovev_n1 (ideal)          22 unchanged
solovev_kinetic_calculated  14 unchanged
solovev_kinetic_nuzero      14 unchanged
solovev_kinetic_ntv          NTV torque fgar        0.10%
                             NTV psi quad evals     60 -> 45  (-25%)
solovev_kinetic_multiion     NTV total (D+T+imp+e)  0.68%
                             NTV Deuterium          0.10%
                             NTV Tritium            0.03%
                             NTV electron           2.62%

That the two ideal cases are bit-identical confirms the grid-pinning rework did not leak into
ideal runs. The electron moving most is the clearest evidence defect #1 was real: the electron's
own ψ_c is ~0, so it previously integrated from psilow while the D/T kinetic matrices had
already been zeroed below their boundary. It now starts where the widest-orbit species says the
model is valid.

Against develop (bundling #398):

diiid_n1                    develop         branch          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_n1                  develop         branch          diff
total energy Re(et[1])      6.773398e-01    1.073540e-02    (dominated by #398)
ODE steps (total)           618             495             -19.9%
q0, q95 unchanged

solovev_kinetic_calculated  develop         branch          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         branch          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         branch          diff
NTV torque fgar             (see #398)      (see #398)      dominated by #398
NTV psi quadrature evals    840             45              -94.6%

solovev_kinetic_multiion    develop         branch          diff
NTV total (D+T+imp+e)       1.315589e-04    3.049754e-01    dominated by #398

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). 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
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).

Unit tests at this state: runtests_kinetic 277/277, runtests_kinetic_validity 27/27
(new), runtests_kinetic_profiles 9/9, runtests_multiion 52/52, runtests_sing 76/76,
runtests_grid_refinement 59/59, runtests_h5_schema 16/16 + 6/6, runtests_fullruns 19/19.

Caveat on coverage: the Solovev calculated example has no [PerturbedEquilibrium] section and
solovev_kinetic_ntv runs at kinetic_factor = 0, so the harness still cannot 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, measured before #391/#392 and before this round's
fixes. The boundary change only alters multi-species or rational-adjacent runs (see the regression
table), so the single-species DIII-D headline is not expected to move.

Merge order: this PR is stacked on #398, which is still open. #414 cannot merge before it.

⚠️ 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]

Review round pushed (2abea51d2)

@ebursch — your three comments are addressed and replied to inline. One is worth repeating here: the tolerance block you flagged in KineticForcesStructs.jl is not this branch's code. It came in on develop with #392 and only shows in this diff because the merge base makes develop's commits look like additions. I trimmed it anyway, but the diff you were reading is wider than the branch — worth knowing for the rest of the review too (the BounceAveraging.jl ω_D change and the Asymptotics.jl/HDF5Schema.jl formatter churn are also develop's, not mine).

I also ran an independent pass while I was in here, and it found six things your review did not reach — the full table is in the PR body under Review round. The one that matters most: the suppression boundary was computed three different ways. The NTV torque quadrature recomputed ψ_c per-species and never cleared rationals, while the EL matrices and the Validity output used the widest-species cleared boundary. So Validity/psi_c in gpec.h5 described a different domain than the torque was actually integrated over — and the PE-vs-NTV torque agreement, which is this PR's headline evidence, was being compared across two different suppression domains. The body's "one source of truth" claim was simply false until now.

The harness confirms the fix does what it should and nothing else: both ideal cases are bit-identical against the pre-fix tree, the two calculated-kinetic cases are unchanged, and the only movement is the NTV quadrature domain. On solovev_kinetic_multiion the electron term moves most (2.62%) — which is exactly the signature of the bug, since the electron's own ψ_c is ~0 and it was integrating from psilow while the D/T matrices had already been zeroed below their boundary.

Also added test/runtests_kinetic_validity.jl (27 tests) — the validity machinery had shipped with none.

---

On your PENTRC question ("could be good to have someone who has used pentrc before to look it over. Not sure who that would be in our group") — docs/development/contributors.md lists @jaesun57 (Sunjae Lee) as the only entry whose focus is KineticForces. I have not added them as a reviewer; that is your call and @logan-nc's.

Merge order is unchanged: this is stacked on #398, still open, so #414 cannot merge before it does — and

⚠️ NO MERGE WITHOUT AN APPROVING THIRD-PARTY HUMAN REVIEW. NON-NEGOTIABLE.

🤖 Generated with Claude Code

https://claude.ai/code/session_01VdRsGGg4YRG6E9fgLGbpdF