ForceFreeStates - BUGFIX - Decouple the ideal-crossing jump width from singfac_min #481

Rewritten 2026-10-13. The original body blamed the trapezoid step and proposed re-projection, ODE sub-steps or extended precision; none of those addresses the actual source. The original text is in the edit history.

Problem

At each ideal singular surface, both integrators (forward and Riccati) stop at ψ_s − d and jump to ψ_s + d, where d = singfac_min/(n·|q′|). The jump zeroes the resonant component of the regular solutions, but that component only vanishes like x at the surface. Zeroing it at x = d instead of 0 leaves an error in W, δW and the energies that is first order in singfac_min:

• 100× smaller jump → 100× smaller change in et₁ (DIIID 2.69e-3 → 3.0e-5; Solovev 3.40e-3 → 3.4e-5), so first order, as expected.
• Riccati et₁ drifts by 0.7% (n=1) to 10% (n=6) on DIIID-like at the default.
• The non-Hermiticity of W after each crossing (the original symptom here) is the same truncation error.

Why lowering singfac_min isn't the fix

singfac_min sets four things at once:

1. the jump half-width (wants to be small);
2. the Riccati Δ′ matching distance (wants to be large: precision degrades as tol·x^(−2α); DIIID Riccati breaks at 1e-7);
3. the chunk ends and kinetic chunking;
4. the knot-free equilibrium grid bracket around each rational (GridRefinement).

After #520 the forward path stays accurate at small singfac_min (Solovev to 1e-8), so forward-only runs can lower it as a workaround, at the cost of moving the grid, the per-surface ca_left/ca_right stubs and kinetic chunking. Riccati cannot.

Proposed fix

Give the jump its own distance and leave singfac_min for everything else. Two candidate designs; prototype both and pick one:

A. ODE legs (prototyped). Add δ_jump (default 1e-6/(min(n)·|q′|)). Inside both crossings (cross_ideal_singular_surf!, riccati_cross_ideal_singular_surf!), integrate the same ODE from ψ_s − d to ψ_s − δ_jump, jump across ±δ_jump, and integrate back out to ψ_s + d. Every handoff point stays at d.

B. Crossing in asymptotic-coefficient space (not prototyped). Write u = V c_A (Glasser 2016 Eq. 50), set the large resonant coefficient to exactly zero, carry the non-resonant coefficients across (Glasser–Kolemen 2018 Eq. 36) and rebuild with ua on the right. This removes the O(d) error at today's 1e-4 with no legs and no new key, and also the d^(2α) residual that resonant-row elimination leaves for α < ½. Open risk: extracting coefficients through an ill-conditioned ua.

Prototype evidence (design A, δ_jump = 1e-6)

• DIIID-like n=1–6: ≤ 1.1e-3 error at δ_jump = 1e-6; forward and Riccati extrapolations in δ_jump agree within 2e-5.
• Riccati Δ′ matrix is bit-identical in every case (the STRIDE BVP reads only ua at the matching distance and the chunk propagators). Kinetic is identical.
• Cost: +5–10% runtime at 1e-6.
• Developer page with figures and the draft diff: https://claude.ai/artifact/MaMkwYpbRxALD99cT6GAWB

Caveat: the prototype predates #520's final trigger (cond(U₁) on the ucrit key, U₂ / lu(U₁)); re-run the table on #520's head. The earlier "DIIID n=3 forward 2.7× wrong" was mostly forward basis collapse, which #520 fixes, not this jump error (Riccati at the default is within ~1%).

Required corrections before a PR (design A)

• Capture ca_left/ca_right and the Δca stub at d, outside the legs. They go to HDF5 and SLAYER's fallback, and ua is far worse conditioned at δ_jump.
• Riccati: the column zeroing assumes (S, I) form, and the left leg ends without a renormalization, so the large solution is effectively removed at d. Renormalize at ψ_s − δ_jump before zeroing, or reuse ForceFreeStates - BUGFIX! - ⚠️ Remove the resonant solution at ideal crossings regardless of reduction timing #489's resonant-row elimination. This is the likely cause of the remaining 7e-5 forward/Riccati gap at DIIID n=3.
• Scale δ_jump with min(n) in multi-n.
• Attribute non-Hermiticity flags with a fixed axis start (frobenius_psi_max = 0); the default Frobenius start adds its own m=0 floor (ForceFreeStates - Repair or retire the Frobenius axis start #464).
• The kinetic crossing (cross_kinetic_singular_surf!) has the same structure at singfac_min; its convergence is untested. Leave it unchanged and say so, or open a follow-up issue.
• Check reverse-shear double rationals (ForceFreeStates - BUGFIX! - ⚠️ Keep both crossings of a rational surface on reverse-shear profiles #492) where two surfaces sit within d.

Impact

• Every free-boundary energy moves toward convergence by roughly singfac_min relative (DIIID et₁ −0.34%), so the PR takes !. Shipped Fortran inputs use 1e-4, so results move away from Fortran benchmarks by about that amount.
• Re-pin the harness and add a small-jump Solovev case.

Dependencies

• ForceFreeStates - BUGFIX! - ⚠️ Remove the resonant solution at ideal crossings regardless of reduction timing #489 (crossing pivots on the resonant row) → ForceFreeStates - BUGFIX! - Trigger the Gaussian reduction on the conditioning of U₁ and factor U₁ in the plasma response matrix #520 (cond(U₁) reduction trigger, resizable reduction storage, stable W). Both open; XB stacks on ForceFreeStates - BUGFIX! - Trigger the Gaussian reduction on the conditioning of U₁ and factor U₁ in the plasma response matrix #520. Without ForceFreeStates - BUGFIX! - Trigger the Gaussian reduction on the conditioning of U₁ and factor U₁ in the plasma response matrix #520, a small jump collapses the forward basis after high-α surfaces.
• Independent of the Δ′-input work (Equilibrium - BUGFIX! - Build F and P of g-file and IMAS equilibria from the tabulated FF′ and p′ #506, ForceFreeStates - BUGFIX! - Use trace-aware Frobenius exponents at resonant surfaces #517, field-line and grid PRs), since Δ′ is unchanged by construction. Those PRs do move the equilibrium, so take the final before/after numbers on the base this lands on. If design B is chosen, land ForceFreeStates - BUGFIX! - Use trace-aware Frobenius exponents at resonant surfaces #517 (trace-aware exponents) first, since B relies on ua.
• The grid-bracket PR must keep its bracket ≥ d so the matching points and legs stay in one cubic cell.

Related

• Large non-hermitian component of plasma response matrix when calculating crit #93: the original report, closed by ForceFreeStates - BUGFIX! - Scale the ODE absolute tolerance to each solution column #480.
• ForceFreeStates - Gaussian reduction trigger misses near-parallel solutions #490: forward basis collapse undetected by ucrit; fixed by ForceFreeStates - BUGFIX! - Trigger the Gaussian reduction on the conditioning of U₁ and factor U₁ in the plasma response matrix #520 (open).
• ForceFreeStates - Repair or retire the Frobenius axis start #464: Frobenius axis-start floor; separate from this.