FFS - BUG - Riccati Δ' carries a singfac_min matching-radius bias that dominates small matrix elements

[d-burg]

Symptom

The Riccati-path Δ' matrix carries a matching-radius (singfac_min) truncation error that the asymptotic Frobenius series does not absorb. It is an absolute error, so it is invisible on large diagonal elements (~1%) and dominant on small ones (~15%). Found while benchmarking against Fortran STRIDE, TJ, and our own Galerkin path on a circular TJ equilibrium (ε = 0.1499, q_edge = 3.6, n = 1; ldp mpsi=256 unless noted).

Δ'_31 (small element, consensus ≈ 1.05–1.07):

singfac_min	Δ'_31	Δ'_21
1e-4 (default)	0.9087	10.8678
3e-5	1.0217	10.8709
1e-5	1.0464	10.8767
3e-6	1.0545	10.8743
1e-6	1.0552 (plateau)	10.8750

Cross-checks at the same point: julia Galerkin path (gal_flag=true) = 1.0567; Fortran STRIDE = 1.05 ± 0.02 (mpsi 128/192/256, at the same singfac_min=1e-4, sing_order=6); TJ (converted to the ψ_N reference) = 1.067. So the default-singfac_min Riccati value is ~14 % low and everything else agrees.

The bias is grid-converged: identical (±0.2 %) across ldp mpsi 128/512/1024 and the auto grid at τ = 1e-3/3e-4 — resolution scans cannot detect it. It is also insensitive to qhigh truncation (0.001 per 0.02 in q_edge) and to eulerlagrange_tolerance (1e-8 → 1e-10 changes nothing).

What it is NOT (measured/audited)

• Not the asymptotic series order. sing_order = 2, 6, 10 give identical Δ' to 4+ digits — expected, since the expansion parameter √dpsi ≈ 5e-3 makes terms beyond k=4 numerically irrelevant at these radii. Raising sing_order cannot fix this.
• Not a port defect in the series or the BVP. A line-by-line comparison against Fortran STRIDE (stride/sing.F, stride/ode.F) found the vmat recursion, sing_get_ua power restoration, matching-point formula (singfac_min/|n·q'|), BVP assembly, boundary blocks, and PEST3 four-term combination all term-for-term equivalent.
• Not spline-knot pollution of the series coefficients: the same bias appears on the auto grid, which brackets rationals away from knots.

Mechanism status

The residual decays as ~dpsi^1.2 — a non-half-integer power outside the span of the Frobenius series, consistent with the one error family the series provably cannot cancel (big-solution contamination of the extracted small-solution coefficient, expected to scale as dpsi^(2α); here 2α = 1.00 at both surfaces, vs 1.22 measured — right family, imperfect exponent match, so mechanism is indicated, not proven). Open question: Fortran shows only −1.4 % over the same singfac decade at identical settings, and the audits found no convention difference that explains its immunity. Remaining unexamined candidates: per-column atol policy (see below), left/right alpha_override sharing, axis-midpoint placement.

Two concrete defects found along the way

1. ForceFreeStatesStructs.jl:272 — singfac_min::Float64 = 1e-4 # Matches Fortran STRIDE — the comment is wrong: Fortran STRIDE's code default is 1e-5 (stride/ode.F:39).
2. Riccati.jl:1411/1421 (and the other solve(...) calls in the FM path) never pass abstol, so DifferentialEquations' uniform default abstol=1e-6 applies to columns spanning ~10⁻³ to ~10³. Fortran scales atol per column to the column's own magnitude (ode.F:411–431). Independent fidelity gap; tightening it is safe (it strengthens, not relaxes, the tolerance).

Proposed fixes — NOT applied, need maintainer sign-off (published Δ' values move)

• P1 (safe, do first): per-column-scaled abstol in integrate_fm_with_ua_ic and the FM solves, mirroring ode.F:422-431.
• P2 (the main fix, needs decision): lower the default matching radius — the measured plateau is at singfac_min ≈ 3e-6–1e-6 (runtime unchanged in tests) — or decouple the chunk-layout radius (EulerLagrange.jl:533) from the ua-matching radius so the layout stays at 1e-4. This moves Δ' (+16 % on the benchmark's 3/1) and every downstream tearing quantity; regression harness + human review required.
• P3: fix the misleading comment either way.

Repro

Scripts + logs: julia_deltaprime_harness/outputs/develop_benchmark/{riccati31_sens.jl,riccati31_sens2.jl,galerkin31_probe.jl,dprime31_probes.jl,stride31_probe.py} in the maintainer's benchmark repo (equilibrium TJ_epsilon_scan_0.1499.geqdsk, settings in-script). Measured on develop @ 9491f89.

🤖 Generated with Claude Code

[d-burg]

Root cause found and fixed — draft PR #413

The matching-radius sensitivity is error control, not convention or series truncation: the shooting solves batched all columns under the ODE library's uniform default abstol = 1e-6, so the resonant big column (~dpsi^−α) set the absolute error floor of the small column (~dpsi^+α), producing the dpsi-dependent contamination. Verified chain:

• sing_order 2/6/10 identical to 4+ digits (expected: √dpsi ≈ 5e-3 makes higher terms numerically irrelevant — raising the order can never fix this).
• Structural confinement of the ua IC to the surface-adjacent chunk (the reference implementation's layout) was implemented and tested: only half-fixed the benchmark and regressed the small-radius plateau (1.055 → 1.011) — reverted.
• Per-column abstol alone (11 lines, mirroring the reference implementation's per-column atol policy): Δ'_31 at the default radius 0.9087 → 1.0544, radius-independent to 0.07 %/decade, in agreement with Galerkin (1.0567), Fortran (1.05 ± 0.02), TJ (1.067). Shaped DIII-D case unchanged ≤ 0.05 %; harness vs merge base: diiid_slayer_n1 15/15 unchanged, riccati/gal cases 0.00–0.01 %.

Also closed en route: Galerkin robustness verified (gal_dx0 over 20× range → 1.052–1.061), and note the default-radius workaround previously suggested here (singfac_min = 1e-6) is NOT safe globally — it degrades high-μ surfaces (shipped-case q=5, μ = 0.83: −2476 → −2110 at 1e-6) — the PR fixes the mechanism instead, making the workaround unnecessary.

Remaining open (minor, not blocking): the misleading singfac_min comment ("Matches Fortran STRIDE" — actual reference default is 1e-5) still stands in the struct docs; and the can_reintegrate == false fallback for uShootL composes forward with a surface-side T_init, which looks direction-inconsistent (latent, fires only without ctrl/equil/ffit) — both worth folding into #413 review or a follow-up.