ForceFreeStates - Investigate residual Galerkin and Riccati Delta-prime disagreement #603

Track the residual disagreement between the ForceFreeStates (FFS) Galerkin and Riccati outer-region Δ′ calculations separately from the fixes in #513. The cause remains unresolved; Riccati is a comparison reference, not an independently established exact answer.

Established result

Galerkin directly interpolated the reduced matrix F̄, while Riccati and the singular-surface expansions represented it using the interpolated Cholesky factor L: F̄=LL†. Aligning Galerkin with that representation reduces the DIII-D relative Δ′ gap from about 10.98% to 0.152% at nx=256. The residual exists with both LU and Cholesky. Across the expanded battery, their full Δ′ matrices agree within 2.53e-9 relative norm.

The gap is ||Δ′_Galerkin−Δ′_Riccati||_F/||Δ′_Riccati||_F for the full complex PEST-3 matrix, including off-diagonal coupling. The large outer-surface DIII-D entry dominates this norm; retain per-surface diagnostics.

Conditions and results

Examples: DIIID-like_gal_resistive_example and LAR_ideal_match_test. Each comparison uses the same equilibrium, Fourier band, rational-surface order and fixed-displacement inner/outer boundary conditions. Riccati uses extended-precision BVP assembly/solve. Galerkin inner-layer matching is disabled. This concerns outer-region Δ′, not matched inner-layer solutions or free-boundary energies.

At nx=1024 DIII-D diagonal relative errors at q=2,3,4,5 are approximately 0.01045%, 0.00561%, 0.00136%, and 0.14512%. The q=5 absolute complex diagonal error is 4.406; Galerkin's imaginary part there is −0.007812. LAR diagonal errors are about 0.00209% and 0.00215%. Radial refinement has not demonstrated agreement of the residual.

Investigations completed

Eleven Galerkin and six Riccati variants per equilibrium tested nx through 1024, jointly scaled Galerkin widths ¼×/½×/2×, expansion orders 4/8/10, quadrature order 10 with tighter resonant tolerance, Riccati tolerances 1e-10/1e-12, half Riccati crossing distance, and combined refinements.

• Expansion order 4–10 has negligible effect.
• DIII-D's baseline Riccati tolerance is already 1e-10. Tightening to 1e-12 changes its matrix by 5.80e-6 relative; halving crossing distance changes it by 9.12e-6. Both are much smaller than the residual (~1.45e-3).
• LAR's baseline tolerance is already 1e-12; its 1e-10 test is a loosening.
• Joint width changes and higher quadrature do not consistently improve both examples. Quarter widths reduce DIII-D's gap to 0.1263% but worsen LAR to 0.00443%.

A second eight-variant study at nx=512 varied only gal_dx0 down to ⅛× while holding element widths fixed, raised gal_gnstep from 5000 to 200000, and split resonant integrals at coefficient-spline knots.

• Δ′ barely changes: DIII-D remains about 0.1471%, LAR about 0.0023029%, against a fresh Riccati reference at tolerance 1e-12.
• All 96 resonant integrations exhausted their budgets without meeting estimated error targets, even with the larger budget. Cancellation/roundoff in the near-zero residual integrand is a hypothesis, not an established explanation of the Δ′ gap.
• DIII-D resonant intervals contain no interior coefficient-spline knots, so splitting there was a no-op. LAR had one affected interval. Ordinary and extension cells crossing knots remain untested.

Next investigation

1. Test knot-aware quadrature in ordinary and extension cells, retaining the original Hermite element coordinates and basis.
2. Check weak-form and surface-term consistency with the factorized operator and singular basis, including cancellation-sensitive terms.
3. Investigate scale-aware stopping/error diagnostics for near-zero residual integrals. Increasing evaluation budgets alone is ineffective; do not weaken tolerances solely to suppress exhausted-budget reports.
4. Establish an error budget per method and per surface before accepting the remaining gap; use an independent analytic/high-precision reference if necessary.

Reproduction notes

Measurements used Julia 1.12.6, 4 Julia threads / 1 BLAS thread, and the CI-pinned Julia 1.12 manifest. Source base was b75eeb3a8a6e6770daeafd9c86d6b1cb49f55e5b plus #513's coefficient-alignment change, uncommitted at measurement time: evaluate F_spline_lower and reconstruct F=QLL†Q in gal_get_fkg. Do not reproduce from the base SHA alone.

Build each equilibrium and coefficient splines once per case. Assemble Galerkin for each control variant, capture its matrix/RHS before factorization, solve copies with LU and auto, and compare PEST-3 blocks. Compute fresh Riccati references from the same equilibrium/splines and copied integration state with fixed-edge BVP conditions and extended_precision_bvp=true. Check surface count/order, retain full complex matrices, and report matrix norm, per-diagonal absolute/relative error and off-diagonal errors.

These diagnostics are separate from #513's final regression harness. This issue records deferred investigation, not a claim of full Galerkin/Riccati agreement.