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). The remaining golden gap is a different fortnum-vs-GSL component, addressed elsewhere. No golden regeneration, no tolerance change.

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.

[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.