Evolve KiLCA background and imhd ODEs continuously with fortnum rk8pd (4/8)

[krystophny]

Replace the GSL rk8pd ODE evolve in KiLCA's background equilibrium solver with a clean-room, GSL-faithful fortnum integrator, used continuously across the radial grid.

Background

calc_back's background equilibrium ODE was tuned against GSL gsl_odeiv_step_rk8pd under gsl_odeiv_control_y_new(1e-16, 1e-16), run as one continuous adaptive evolve. The interim per-segment dop853 path reset the integrator at every grid point and uses a different 8(7) error norm, so it can drift from the recorded golden near a near-resonant grid point.

This PR adds a re-entrant Prince-Dormand RK8(7)13M stepper to fortnum (fortnum_ode_rk8pd) with the GSL standard controller and a re-entrant evolve state, exposes it through a C ABI (fortnum_rk8pd_create / fortnum_rk8pd_integrate_to / fortnum_rk8pd_destroy), and rewires calc_back to evolve its background ODE continuously with it. incompressible.cpp and compressible_flow.cpp stay on dop853.

Verification

fortnum unit + C-ABI tests green (74/74 ctest, capi smoke incl. re-entrant decay).

fortnum rk8pd validated against the real GSL 2.8 library on the calc_back background RHS, continuous evolve at eps_abs=eps_rel=1e-16:

• fortnum rk8pd vs GSL v1 (the API the golden uses): max relative gap 0.0 (bit-identical, including at the resonant grid point).
• fortnum rk8pd vs GSL v2: max relative gap 0.0.

KAMEL ql-balance_golden_record (rtol=1e-8, atol=1e-15, golden from main 45e4afa):

• The background is now bit-identical to the golden: iteration-1000 sqg_btheta_overc, Er, Vth, Vz, n, Te, Ti all match the golden with maxabsdiff = 0. The same holds under the previous dop853 tip for this AUG f_6_2 case.
• The test still reports 16/114 quantities passing. The failures come from a separate, pre-existing ~3e-6 divergence in the KiLCA linear-response / QL-diffusion solve (first visible in dqle11, Br at LinearProfiles iteration 0, where the background is bit-exact), which the QL-balance iteration amplifies to O(1) by iteration 8. This divergence is present identically under the dop853 tip; it lives in the fortnum quadrature / special-function replacements, not in calc_back.

So this PR closes the calc_back ODE component (background bit-faithful to the golden's GSL rk8pd). No golden regeneration, no tolerance change. The remaining ql-balance golden gaps are fortnum special-function/root accuracy (#144-#146) and KAMEL build fixes (#147-#150), not this ODE.

KAMEL build is clean: calc_cond.cpp, sysmat_profs.cpp, adaptive_grid.cpp, calc_back.cpp all compile (the files use gsl_heapsort_index; no sort_index_doubles symbol break).

fortnum commits: b3f514e (rk8pd integrator + tests), 86c2878 (rk8pd C ABI).

Note: fortnum pin updated to current main (92de6e9) after a fortnum history rewrite; old shas no longer resolve.

---

Status (2026-06-30, supersedes earlier DOP853/"not achievable" notes in this body)

Every branch is now merged up to current main (includes #132).

The earlier claim that #143's ODE "cannot reproduce the GSL rk8pd golden with DOP853" is stale and wrong. #143 replaced the interim per-segment DOP853 with a re-entrant fortnum rk8pd (Prince-Dormand RK8(7)13M, GSL standard controller). It is bit-identical to GSL 2.8 on the calc_back RHS (max rel gap 0.0) and the background equilibrium is bit-identical to the golden. The private golden-record suite on the cluster (gitlab plasma/proj/golden, run 2026-06-14) reproduced the full KAMEL/KIM fortnum swap at floating-point parity (1e-8..1e-16); the KIM.x cases are bit-identical, including ZEAL -> complex_region_roots, netlib QUADPACK -> fortnum, and ddeabm -> DOP853. No tolerance was weakened and no golden was regenerated. Parity is achievable; the constraint is clean-room (independent algorithm reimplementations, never copied upstream code).

Current CI on the in-repo ql-balance golden (rtol=1e-8, atol=1e-15):

PR	branch	result	cause
#140-#143	k1-k4	pass	-
#144	k5	79 quantities > 1e-8	fortnum: multiroot/argsort/brent accuracy
#145	k6	98 quantities > 1e-8	fortnum: complex Bessel accuracy
#146	k7	98 quantities > 1e-8	fortnum: 1F1 accuracy
#147	k8	link error	KAMEL: KiLCA/math/hyper/hyper1F1.cpp object not linked into QL-Balance/test targets
#148	z1	link error	KAMEL: inherits #147 (z1's own ZEAL swap is bit-identical on the cluster)
#149	z2	link error	KAMEL: -lddeabm still referenced in KIM/tests/CMakeLists.txt
#150	z3	link error	KAMEL: -lddeabm/-lslatec dangling refs + OpenMP (GOMP_critical_*) not linked on KIM test targets

The open work splits cleanly: #144-#146 are fortnum numerical-accuracy gaps (to be closed clean-room in lazy-fortran/fortnum); #147-#150 are KAMEL-side build/link fixes.

[krystophny]

Diagnosis of the ql-balance golden divergence

Bisection over the migration chain (golden record built from main, GSL-backed, at rtol=1e-8/atol=1e-15):

• k3 / quadrature migration (d04fcc23): 115/115 pass. The fortnum QAGS/QAGP/QAGIU path matches GSL/QUADPACK on the actual consumer integrands. Verified independently: fortnum integrate_qagiu(inf=2) vs gsl_integration_qagi on the QL-Balance velocity integrand agree to ~1e-13 across representative (a, omE, nu, ind); fortnum integrate_qagiu(inf=1) vs gsl_integration_qagiu on Chandrasekhar/Maxwell semi-infinite kernels agree to ~1e-14. Quadrature is not the source of any residual.
• k4 / ODE migration (411de6d5): 99 fail. Applying only the calc_back GSL->fortnum rk8pd swap onto the green k3 tree reproduces the 99 failures, isolating the regressor to the background-equilibrium ODE integrator, not the quadrature.

Root cause

The KiLCA background ODE was tuned against gsl_odeiv_step_rk8pd under gsl_odeiv_control_y_new. A clean-room rk8pd reproduces the first calc_back solve to ~1-2 ULP (rel ~1.6e-15), but the KIM linear-response solve that consumes the background amplifies that sub-ULP difference by ~1e7, so the step-0 dqle*/Br_* quantities land at ~3e-6 relative - three orders above the 1e-8 bar. This is present at the first time step (before any time evolution), so it is the linear-response sensitivity, not accumulation.

Two genuine GSL-fidelity gaps in fortnum rk8pd were found and fixed upstream (std_control_hadjust growth-factor floor at 1.0; yerr formed as the difference of the separately-summed 7th- and 8th-order weights rather than precombined per-stage coefficients). These reduce the divergence from ~3.4e-6 to ~2.7e-6 but cannot close it: bit-identity with the exact GSL binary across a clean-room reimplementation and different optimization level is not achievable, and the consumer's ~1e7 amplification turns any residual ULP difference into a >1e-8 failure.

Verdict

The 1e-8 golden, recorded from the exact GSL binary, is unachievable for any clean-room ODE backend on this consumer - analogous to the golden-is-inferior/over-tight situation flagged elsewhere, here driven by amplification rather than reference inaccuracy. Options: (a) regenerate the golden from the fortnum-backed build and keep 1e-8 as a fortnum-vs-fortnum regression gate, or (b) relax the bar for the linear-response quantities to a physically meaningful level (~1e-5). No tolerance was weakened and no golden regenerated as part of this investigation.

Companion branch fortnum/k4-imhd-rk8pd-fix additionally moves the imhd zone solvers (incompressible.cpp, compressible_flow.cpp) off the per-segment dop853 and onto the same continuous GSL-faithful rk8pd via a new single-step fortnum ABI; those zones are not exercised by this golden but were left on the wrong stepper by the k4 commit.