KineticForces - BUGFIX - Decide cubic and quadratic root counts by conditioned tests

[logan-nc]

Release note

• Audience: developers
• Numerical impact: none expected — every tracked quantity unchanged across the five kinetic cases (harness @ be686c3)
• Migration: not required; internal helpers only, no interface or configuration changes

The closed-form root solvers behind the NTV bounce-point search now decide how many real roots a cubic or quadratic has by a conditioned test rather than by the sign of a discriminant that cancels to rounding noise when two roots nearly coincide. Previously such a case could silently drop real roots; the bounce-point pipeline fell back to its adaptive scan and gave the right answer, but the helper is generic and would not have in another caller.

What this fixes

Review of #456 found, and reproduced with an out-of-tree probe, that _real_cubic_roots misclassifies a cubic with three genuine real roots when two are about 1e-8 apart: disc = (q/2)^2 + (p/3)^3 computes to +8.8e-20 for a true negative, lands just outside the repeated-root tolerance, and the one-real-root branch returns only the far root. Independent review of the merged code found the same cancellation at two more sites.

site	failure mode	consequence before this fix
_real_cubic_roots, depressed-cubic discriminant	near-double root reports one root, drops two	NaN in the monotone window, lambda emptied, find_zeros fallback — correct answer, pre-#456 cost
_real_cubic_roots, quadratic branch for a degenerate leading coefficient	near-double root reports no roots	same fallback
_surface_b_field, stationary-point quadratic	near-double stationary point reports none	an interval that is not monotone, with no fallback — a root pair lost silently

The third is physically a near-flat inflection, so the lost pair is nearly coincident and the impact small; but the discipline has to be consistent across all three sites.

The fix

• _real_cubic_roots decides the count with the conditioned ratio: for p < 0, three real roots iff |q/2| <= r^3 with r = sqrt(-p/3), tested as the ratio against 1, where nothing cancels. A ratio a rounding step past 1 is a repeated root, which the clamped trigonometric form returns directly — so the explicit repeated-root branch and its ad-hoc tolerance are deleted. One real root otherwise (Cardano, radicals do not cancel there).
• New _quadratic_real_roots treats a discriminant within rounding of zero as a double root rather than none, and serves both the degenerate-cubic branch and the stationary-point search.
• One shared constant, ROOT_COUNT_TOL = 64 eps. Deliberately generous: a wrongly admitted repeated root is a near-touch at a stationary point, which cannot sit strictly inside the monotone window _cell_level_root filters on, so it can only lose the residual ranking there, never win it. The measured margin for the reported case is 6 eps, so the earlier 8 eps would have been one FMA reordering away from failing on Julia 1.12.

Tests

5073 assertions in runtests_kinetic.jl (from 1844). New coverage:

• the near-double cubic, built from its roots (0.37, 0.37000001, 0.91) and two more triples so the polynomial has them exactly — all three return, residual < 1e-12;
• near-double and exact-double quadratics, same construction, plus the degenerate linear and no-real-root forms;
• the previously untested path where both stationary points bounding an interval fall inside one spline cell — found deterministically on an 11-knot grid (2 + 0.45 cos 2πθ + 0.20 cos 6πθ: six stationary points, a pair 0.077 apart sharing a cell), swept across the full B-range between them against Roots.find_zeros: counts agree, the crossing lies strictly between the two stationary points, and |v_par| < 1e-14 at every root.

A fixture that drives the stationary-point discriminant itself to the rounding boundary through a spline is not attempted — that is fiddly and non-deterministic across Julia versions; the helper test is the unit check and the existing stationary-point-count test (against a dense deriv1 sign scan) is the integration check.

Verification

Real DIII-D geometry, 40 surfaces, every trapped lambda on the production grid plus lambdas injected to within 1e-13 of both trapped/passing boundaries, against Roots.find_zeros on the identical spline:

check	#456 as merged	this fix
root-set mismatches on the production lambda grid	0 of 2,560	0 of 2,560
falls back to find_zeros	2	0
cases where the direct solve found fewer roots than find_zeros	2	0
cases where it found more (injected near-tangent lambda only)	78	80
max abs(v_par) at our roots, production grid	4.44e-16	4.44e-16
max abs(v_par) at our roots, injected lambda within 1e-13 of a boundary	4.44e-16 (via fallback)	5.5e-14

The two previous fallbacks were exactly the misclassification firing on injected lambda; they now solve directly. The one figure that moved is the residual at those same two lambdas: they are genuinely near-double roots, where the single Newton polish in _cell_level_root converges linearly rather than quadratically, so the direct path leaves 5.5e-14 where find_zeros left 4.4e-16. That is the conditioning limit of a near-double root, at a point where the root location itself is only determined to about 1e-7; the production grid never gets within 1e-4 of a boundary and stays at machine precision. A multiplicity-aware iteration would fix it and is not worth the machinery for a regime the lambda grid cannot reach.

Regression report

case	result
diiid_n1	53 unchanged
solovev_kinetic_ntv	6 unchanged
solovev_kinetic_multiion	6 unchanged
solovev_kinetic_calculated	14 unchanged
solovev_kinetic_nuzero	14 unchanged

Baseline develop@56eafc2a5, branch @be686c388 (develop merged in; comment-only commit since the last run); no stale-baseline warning. Every tracked quantity unchanged, as predicted below — including solovev_kinetic_calculated and solovev_kinetic_nuzero, which moved at the 1e-8 level under #456.

Regression Report (5 cases)
Ref 1: develop  @ 56eafc2a5
Ref 2: bugfix/kf-cubic-root-classification  @ be686c388

NTV torque FGAR [N·m]                         5.336823e-01     5.336823e-01                         0.0e+00    OK    
NTV kinetic energy dW FGAR [J]                7.997363e-02     7.997363e-02                         0.0e+00    OK    
NTV total (D+T+imp+e) [N·m]  1.317860e-04   1.317860e-04                         0.0e+00  OK    
NTV Deuterium [N·m]          7.794788e-05   7.794788e-05                         0.0e+00  OK    
NTV Tritium [N·m]            8.206022e-05   8.206022e-05                         0.0e+00  OK    
NTV electron [N·m]           -2.822212e-05  -2.822212e-05                        0.0e+00  OK    

diiid_n1                    Summary: 53 unchanged
solovev_kinetic_ntv         Summary: 6 unchanged
solovev_kinetic_multiion    Summary: 6 unchanged
solovev_kinetic_calculated  Summary: 14 unchanged
solovev_kinetic_nuzero      Summary: 14 unchanged

Notes for reviewers

Prediction, stated before the harness ran: every row unchanged, including the two Mode A cases that moved under #456. The real-geometry sweep for #456 showed zero fallbacks on the production lambda grid, so the NaN path never fired there; where three real roots genuinely exist the trigonometric branch is the identical formula; and the stationary-point change only acts when a discriminant is within rounding of zero, which clean geometry never hits. Any movement would have been a red flag to investigate, not a reason for a !.

No review package for this one — a small bug fix, below the threshold for it.

Credit. The cubic misclassification and the untested same-cell path were both found by @adrianaghiozzi's review of #456, probe and all; the two further sites with the same cancellation came from re-reviewing the merged code with her example in hand.

🤖 Generated with Claude Code

https://claude.ai/code/session_0186CLfjeJeSusNTva5RMfGU

[adrianaghiozzi]

Reviewed with Claude Code. Verified independently in an isolated worktree: the exact #456 misclassification case, and a full sweep of the old failure band (1e-7 to 1e-17, both signs), now resolve correctly. test/runtests_kinetic.jl passes 5073/5073. Triple-root and monotone-cubic edge cases the old explicit branches handled are still correct.