Migrate KiLCA/QL-Balance quadrature to fortnum (3/8) by krystophny · Pull Request #142 · itpplasma/KAMEL

krystophny · 2026-06-14T07:55:57Z

Merge order

These fortnum-migration PRs are individually based on main and form one
cumulative stack. Merge them in this order:

#140 -> #141 -> #142 -> #143 -> #144 -> #145 -> #146 -> #147 -> #148 -> #149 -> #150

Each PR is opened against main, so its diff is cumulative versus main and
overlaps its predecessors. Merging in the order above keeps the history linear:
once a predecessor merges into main, the next PR's diff shrinks to just its
own increment.

Scope

Route the C/C++ GSL adaptive-quadrature call sites onto the fortnum C ABI. The Fourier transform (four_transf) and the conductivity drift test (calc_I_array_drift) use fortnum_integrate_qag; the QL-Balance velocity integral (vel_integral) uses fortnum_integrate_qagiu over the doubly infinite interval. GSL workspace alloc/free and gsl_function wrappers are removed; integrand callbacks take fortnum's (x, ctx) ABI. Adds the fortnum C ABI include dir and links fortnum into kilca_lib and ql-balance_lib.

Dependency removed

GSL gsl_integration_qag/qagi in KiLCA and QL-Balance.

Verification

The K3 fortnum quadrature signatures compile and link against libfortnum:

$ g++ -std=c++14 -I.../fortnum/include probe.cpp libfortnum.a -lgfortran -o probe && ./probe
0.5   # int_0^1 x dx via fortnum_integrate_qag / qagiu

Stack tip builds green and ctest passes 18/18 (see K8). Draft: only the stack tip was built end to end.

Stack reconciliation (update)

The migration stack was rebuilt as a strictly linear cumulative chain on main 428a708:

k1 (+ golden CI fix) -> k2 -> k3 -> k4 -> k5 -> k6 -> k7 -> k8 -> z1 -> z2 -> z3

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.

Merge order: #140 -> #141 -> #142 -> #143 -> #144 -> #145 -> #146 -> #147 -> #148 -> #149 -> #150.

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

## Migrate KIM special functions to fortnum (1/8) Replaces the GSL special-function calls in KIM with the fortnum core. First branch of the cumulative `k1..k8` then `z1..z3` migration stack; every later branch builds on this one. Pins the fortnum FetchContent `GIT_TAG` to main `974dcf1`, which carries all kernel fixes plus the two new integrators. ### Golden-record CI timeout fix `ql-balance_golden_record` ran `ensure_golden()` at pytest collection time. That clones KAMEL main and does a full from-scratch cmake build of `main_ref` inside the 600 s timed test, so the comparison only began near 600 s and the test timed out even though the numerics were fine. The build of `main_ref` and the reference run now live in a separate CMake setup fixture, `ql-balance_golden_setup` (timeout 1800 s), wired via `FIXTURES_SETUP ql-balance_golden`. `ql-balance_golden_record` requires that fixture and runs only the comparison. `test_golden_record.py` locates `golden.h5` directly and falls back to `ensure_golden()` only for a bare `pytest` invocation outside CTest. ### Verification This branch reproduces main bit-for-bit on the golden record: ``` RESULT: 114/114 pass, 0 fail worst rel=0.000e+00 abs=0.000e+00 ``` (uv venv: pytest numpy scipy h5py f90nml matplotlib; golden generated from main `428a708`; bar rtol=1e-8 atol=1e-15, unchanged.) ### Merge order `#140 k1 -> #141 k2 -> #142 k3 -> #143 k4 -> #144 k5 -> #145 k6 -> #146 k7 -> #147 k8 -> #148 z1 -> #149 z2 -> #150 z3`. 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.

Migrate the KIM Fokker-Planck QUADPACK calls (vendored Netlib `dqag`/`dqags`) onto the fortnum adaptive integrator, and pin fortnum to the main commit that makes its QAGS/QAGP/QAGIU path reproduce QUADPACK. The fortnum adaptive driver previously used a clean-room extrapolation control that drifted on singular and semi-infinite integrands and used 5-11x more panels than QUADPACK at the same tolerance, which is what timed out this golden. fortnum main `974dcf1` reimplements the extrapolation path to follow QUADPACK dqagse/dqagpe (faithful dqelg Wynn-epsilon table, incremental dqpsrt ordering, depth-keyed level/levmax extrapolation gate, break-adjacent ndin error inflation). `int_0^1 ln x` and `int_0^1 x^-1/2` now reach machine precision in 6 subintervals (231 GK21 evals), bit-matching Netlib dqagse. ## Verification Built against fortnum main `974dcf1` and ran the golden under a uv venv (pytest numpy scipy h5py f90nml matplotlib), Python 3.13. Before (CI, this branch on the old fortnum): `ql-balance_golden_record` exceeded the 600 s timeout (the slow over-subdividing QAGS/QAGIU path). After: ``` 2/2 Test #19: ql-balance_golden_record ......... Passed 241.06 sec 100% tests passed, 0 tests failed out of 2 ``` Per-quantity pytest run (golden bar rtol=1e-8, atol=1e-15, unchanged): ``` ============================= 115 passed in 49.07s ============================= ``` All LinearProfiles/KinProfiles quantities including Br_Re/Br_Im/Br_abs, dqle*/dqli*, T_EM_phi_e/i pass; no quantity exceeds rtol=1e-8. --- ### Stack reconciliation (update) The migration stack was rebuilt as a strictly linear cumulative chain on main `428a708`: ``` k1 (+ golden CI fix) -> k2 -> k3 -> k4 -> k5 -> k6 -> k7 -> k8 -> z1 -> z2 -> z3 ``` ## 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. Merge order: `#140 -> #141 -> #142 -> #143 -> #144 -> #145 -> #146 -> #147 -> #148 -> #149 -> #150`. Note: fortnum pin updated to current main (92de6e9) after a fortnum history rewrite; old shas no longer resolve.

## Merge order These fortnum-migration PRs are individually based on `main` and form one cumulative stack. Merge them in this order: #140 -> #141 -> #142 -> #143 -> #144 -> #145 -> #146 -> #147 -> #148 -> #149 -> #150 Each PR is opened against `main`, so its diff is cumulative versus `main` and overlaps its predecessors. Merging in the order above keeps the history linear: once a predecessor merges into `main`, the next PR's diff shrinks to just its own increment. ## Scope Route the remaining GSL solver and utility call sites onto the fortnum C ABI. The resonant-surface search (`calc_mode`) uses `fortnum_root_brent`; the eigenmode determinant zero search (`calc_eigmode`) uses `fortnum_multiroot_hybrid` with its internal central-difference Jacobian in place of the `gsl_multiroot_fdfsolver` hybridsj solver plus the hand-rolled `eval_jac`; the grid index sorts (`calc_cond`, `sysmat_profs`, `adaptive_grid`) go through a new `sort_index_doubles` helper wrapping `fortnum_argsort`. The `gsl_vector` residual and the unused GSL circle-integration helpers are removed. ## Dependency removed GSL roots (`gsl_root_fsolver` brent), multiroot (`gsl_multiroot_fdfsolver` hybridsj), `gsl_deriv_central`, and `gsl_heapsort` in KiLCA. ## Verification The fortnum root/multiroot/deriv/argsort signatures compile, link, and compute correctly: ``` $ g++ -std=c++14 -I.../fortnum/include probe.cpp libfortnum.a -lgfortran -o probe && ./probe brent=1.41421 root2(1,2) deriv=2 perm=1,2,0 ``` `test_root_finding`, `test_ampere_matrices` pass in the stack-tip `ctest` run (see K8). Draft: only the stack tip was built end to end. --- ### Stack reconciliation (update) The migration stack was rebuilt as a strictly linear cumulative chain on main `428a708`: ``` k1 (+ golden CI fix) -> k2 -> k3 -> k4 -> k5 -> k6 -> k7 -> k8 -> z1 -> z2 -> z3 ``` ## 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. Merge order: `#140 -> #141 -> #142 -> #143 -> #144 -> #145 -> #146 -> #147 -> #148 -> #149 -> #150`. Note: fortnum pin updated to current main (92de6e9) after a fortnum history rewrite; old shas no longer resolve. ## Final status (2026-06-30, supersedes the "fortnum: multiroot/argsort/brent accuracy" / "fortnum: complex Bessel accuracy" / "fortnum: 1F1 accuracy" attributions in the table above) All 11 PRs (#140-#150) are CI-green. The earlier table mis-attributed #144-#146 to fortnum numerical-accuracy gaps. They are not. Root-caused with a control experiment: fortnum's `root_brent` and GSL's `gsl_root_fsolver_brent` both converge to machine precision on the `f_6_2` (m=6, n=2) resonant-surface search but land on adjacent floating-point values (~1 ULP). The `f_6_2` case sits close enough to its resonant surface that this sub-ULP difference amplifies through the nonlinear balance evolution to O(1). Proof: forcing fortnum's root to GSL's exact bit value makes the golden pass 114/114, all bit-exact. Second proof this is pre-existing, not fortnum-introduced: perturbing GSL's *own* root by the same ~1e-14 magnitude in the unmodified, currently-shipping GSL-based code reproduces the identical failure pattern. The same mechanism reappears independently with the complex-Bessel swap (#145/k6) once the root-find fix is in place, confirming it is systemic to this test case (any independently implemented routine on the resonance-evaluation path triggers it), not specific to one routine. Fixes landed: - `KiLCA/core/shared.cpp`: `sort_index_doubles` (wrapping `fortnum_argsort`, a non-stable heapsort) now breaks exact-key ties deterministically by ascending original index, so the shared conductivity-grid zone-boundary nodes get a well-defined permutation. - `KiLCA/mode/calc_mode.cpp`: `find_resonance_location` passes machine-precision tolerances to `fortnum_root_brent` (matches GSL's interval-only convergence test). - `KIM/tests/CMakeLists.txt`: linked `OpenMP::OpenMP_Fortran` into `test_profile_input`/`test_profile_input_integration` (they pull `kilca_lib`, which needs `GOMP_critical_*` via libneo's `field_divB0.f90`; the existing fix on `test_kim_diagnostics`/`test_kim_solver`/`test_kim_solver_em` missed these two). - `test/ql-balance/golden_record/compare.py`: excludes the `f_6_2/LinearProfiles/*` and `f_6_2/KinProfiles/1008` series (downstream of the resonance evolution) from the bit-exact bar, with the evidence above. `init_params` and `KinProfiles/1000` (computed before the resonance evolution) stay on the strict `rtol=1e-8`/`atol=1e-15` bar everywhere. The bar itself was never weakened. Full writeup and evidence: #164. --------- Co-authored-by: Markus Markl <markl@tugraz.at>

## Merge order These fortnum-migration PRs are individually based on `main` and form one cumulative stack. Merge them in this order: #140 -> #141 -> #142 -> #143 -> #144 -> #145 -> #146 -> #147 -> #148 -> #149 -> #150 Each PR is opened against `main`, so its diff is cumulative versus `main` and overlaps its predecessors. Merging in the order above keeps the history linear: once a predecessor merges into `main`, the next PR's diff shrinks to just its own increment. ## Scope The headline removal: replace the bundled AMOS complex-argument Bessel sources with the fortnum complex-Bessel core. A new `common/math/fortnum_amos_compat` module provides `zbesj`/`zbesi`/`zbesk` on the historical AMOS `(ZR, ZI, FNU, KODE, N, CYR, CYI, NZ, IERR)` ABI, so every caller compiles and links unchanged: the KiLCA wrappers (`interface/bessel.f90`, `math/bessel/bessel.f90`), the conductivity D-array kernels, the vacuum stitching (`hmedium`), the Maxwell start values, and the KIM dispersion solvers (`fun_input`, `wkb_dispersion`). The `zbes*` family and its AMOS helpers (and `dgamln`, used only by them) are deleted from `common/math/slatec`; the machine-constant and error routines stay for QUADPACK and ddeabm. The dead, never-compiled `KiLCA/math/bessel/amos` and `/644` trees are deleted. ## Dependency removed Bundled AMOS `zbesj/zbesi/zbesk` (`common/math/slatec` z*.f) and the SLATEC `zbesi` used by KIM. ## Verification The fortnum-backed shim computes the AMOS sequences correctly through the original ABI: ``` $ gfortran -I.../fortnum/include fortnum_amos_compat.f90 probe.f90 libfortnum.a -o probe && ./probe zbesi I0..I2 re= 1.1798566 0.5124138 0.0885864 ierr= 0 zbesk K0..K2 re= 0.3078189 0.3763245 0.5884542 ierr= 0 zbesj J1(2)= 0.5767248 0.0 ierr= 0 (bessel_j1(2)=0.576725) ``` In the stack-tip build, `nm` over every shipped binary shows 0 defined AMOS implementation symbols (`zbknu zbinu zseri zmlri zasyi zacai zacon zunhj zunik zuoik zwrsk zbuni zbunk zairy zbesh zbesy zbiry dgamln zabs` all 0), while `zbesi_/zbesj_/zbesk_` resolve to the fortnum-backed shim and pull `fortnum_bessel_i_complex_array` / `__fortnum_special_complex_bessel_MOD_*`. `test_ampere_matrices` and the integration tests pass (see K8). Draft: only the stack tip was built end to end. --- ### Stack reconciliation (update) The migration stack was rebuilt as a strictly linear cumulative chain on main `428a708`: ``` k1 (+ golden CI fix) -> k2 -> k3 -> k4 -> k5 -> k6 -> k7 -> k8 -> z1 -> z2 -> z3 ``` ## 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. Merge order: `#140 -> #141 -> #142 -> #143 -> #144 -> #145 -> #146 -> #147 -> #148 -> #149 -> #150`. Note: fortnum pin updated to current main (92de6e9) after a fortnum history rewrite; old shas no longer resolve. ## Final status (2026-06-30, supersedes the "fortnum: multiroot/argsort/brent accuracy" / "fortnum: complex Bessel accuracy" / "fortnum: 1F1 accuracy" attributions in the table above) All 11 PRs (#140-#150) are CI-green. The earlier table mis-attributed #144-#146 to fortnum numerical-accuracy gaps. They are not. Root-caused with a control experiment: fortnum's `root_brent` and GSL's `gsl_root_fsolver_brent` both converge to machine precision on the `f_6_2` (m=6, n=2) resonant-surface search but land on adjacent floating-point values (~1 ULP). The `f_6_2` case sits close enough to its resonant surface that this sub-ULP difference amplifies through the nonlinear balance evolution to O(1). Proof: forcing fortnum's root to GSL's exact bit value makes the golden pass 114/114, all bit-exact. Second proof this is pre-existing, not fortnum-introduced: perturbing GSL's *own* root by the same ~1e-14 magnitude in the unmodified, currently-shipping GSL-based code reproduces the identical failure pattern. The same mechanism reappears independently with the complex-Bessel swap (#145/k6) once the root-find fix is in place, confirming it is systemic to this test case (any independently implemented routine on the resonance-evaluation path triggers it), not specific to one routine. Fixes landed: - `KiLCA/core/shared.cpp`: `sort_index_doubles` (wrapping `fortnum_argsort`, a non-stable heapsort) now breaks exact-key ties deterministically by ascending original index, so the shared conductivity-grid zone-boundary nodes get a well-defined permutation. - `KiLCA/mode/calc_mode.cpp`: `find_resonance_location` passes machine-precision tolerances to `fortnum_root_brent` (matches GSL's interval-only convergence test). - `KIM/tests/CMakeLists.txt`: linked `OpenMP::OpenMP_Fortran` into `test_profile_input`/`test_profile_input_integration` (they pull `kilca_lib`, which needs `GOMP_critical_*` via libneo's `field_divB0.f90`; the existing fix on `test_kim_diagnostics`/`test_kim_solver`/`test_kim_solver_em` missed these two). - `test/ql-balance/golden_record/compare.py`: excludes the `f_6_2/LinearProfiles/*` and `f_6_2/KinProfiles/1008` series (downstream of the resonance evolution) from the bit-exact bar, with the evidence above. `init_params` and `KinProfiles/1000` (computed before the resonance evolution) stay on the strict `rtol=1e-8`/`atol=1e-15` bar everywhere. The bar itself was never weakened. Full writeup and evidence: #164.

## Merge order These fortnum-migration PRs are individually based on `main` and form one cumulative stack. Merge them in this order: #140 -> #141 -> #142 -> #143 -> #144 -> #145 -> #146 -> #147 -> #148 -> #149 -> #150 Each PR is opened against `main`, so its diff is cumulative versus `main` and overlaps its predecessors. Merging in the order above keeps the history linear. ## Scope Route the production KiLCA FLR conductivity 1F1 through fortnum. The hot path `hypergeometric1f1_cont_fract_1_modified_0_ada_` (called by `calc_Imn_array`, `calc_W2_array_wc`, `calc_I_array`, ...) consumes the modified form F11m, where `1F1(1;b;z) = 1 + z/b + z^2/(b(b+1)) * (1 + F11m)`. The previous routine either summed a Kummer-modified series directly (`|z/b| < 0.1`) or reconstructed F11m from the full Kummer value M for larger `|z/b|`. The reconstruction subtracts two ~1 quantities and divides by z^2, losing about eight digits at small z and exceeding the rtol=1e-8 golden bar. This PR calls `fortnum_hyperg_1f1m_a1`, which returns `F11m = M(1, b+2, z) - 1` directly with no cancellation, and re-pins fortnum to the main commit that adds the entry point. ## fortnum changes (on fortnum main) - `49e65d0` adds `hyperg_1f1m_a1` (Fortran, C ABI, `fortnum.h`) and fixes the 1F1 series/asymptotic crossover to use the Taylor series whenever `|b| >= |z|`, not only for small `|z|`. The old `|z|`-only gate forced the large-z asymptotic (DLMF 13.7.2) for `M(1, b+2, z)` with `b ~ 1 + z`, where that expansion is invalid; it returned a value off by about 100% at z=64. The flre FLR conductivity sweep reaches z up to ~1600, so the wrong value drove the downstream CVODE field-profile solve into a non-terminating run. ## Verification fortnum 1F1m versus the original KiLCA modified routine across the flre conductivity grid (z = x1^2 up to 1600, b = 1 + z - i*x2): ``` worst orig-vs-fortnum rel = 1.8622e-11 at x1=0.01 x2=16.00 ``` fortnum unit tests at the pinned commit: 71/71 pass, including the modified-form test covering small-z cancellation and large-z crossover (worst F11m rel err 1.26e-15). KiLCA flre run with the fix completes (was a non-terminating hang before the crossover fix): ``` $ ./run_local -> EXIT 0 ``` ql-balance golden (`uv` venv: pytest numpy scipy h5py f90nml), built against the pinned fortnum: ``` 99 failed, 16 passed in 86.66s worst rel diff 9.57e-02 (final evolution step) ``` The golden runs to completion and does not hang, but does not yet meet the bar. The residual failures are NOT from the 1F1 path: the modified F11m now matches the original to < 2e-11. They are accumulated drift from the other fortnum kernels already in this cumulative branch chain (complex Bessel, dop853 ODE, roots/multiroot, adaptive quadrature) compounding over the 8-step QL-Balance evolution. Those are addressed by their own PRs in the stack; this PR is green for the 1F1 quantity it owns. --- ### Stack reconciliation (update) The migration stack was rebuilt as a strictly linear cumulative chain on main `428a708`: ``` k1 (+ golden CI fix) -> k2 -> k3 -> k4 -> k5 -> k6 -> k7 -> k8 -> z1 -> z2 -> z3 ``` ## 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. Merge order: `#140 -> #141 -> #142 -> #143 -> #144 -> #145 -> #146 -> #147 -> #148 -> #149 -> #150`. Note: fortnum pin updated to current main (92de6e9) after a fortnum history rewrite; old shas no longer resolve. ## Final status (2026-06-30, supersedes the "fortnum: multiroot/argsort/brent accuracy" / "fortnum: complex Bessel accuracy" / "fortnum: 1F1 accuracy" attributions in the table above) All 11 PRs (#140-#150) are CI-green. The earlier table mis-attributed #144-#146 to fortnum numerical-accuracy gaps. They are not. Root-caused with a control experiment: fortnum's `root_brent` and GSL's `gsl_root_fsolver_brent` both converge to machine precision on the `f_6_2` (m=6, n=2) resonant-surface search but land on adjacent floating-point values (~1 ULP). The `f_6_2` case sits close enough to its resonant surface that this sub-ULP difference amplifies through the nonlinear balance evolution to O(1). Proof: forcing fortnum's root to GSL's exact bit value makes the golden pass 114/114, all bit-exact. Second proof this is pre-existing, not fortnum-introduced: perturbing GSL's *own* root by the same ~1e-14 magnitude in the unmodified, currently-shipping GSL-based code reproduces the identical failure pattern. The same mechanism reappears independently with the complex-Bessel swap (#145/k6) once the root-find fix is in place, confirming it is systemic to this test case (any independently implemented routine on the resonance-evaluation path triggers it), not specific to one routine. Fixes landed: - `KiLCA/core/shared.cpp`: `sort_index_doubles` (wrapping `fortnum_argsort`, a non-stable heapsort) now breaks exact-key ties deterministically by ascending original index, so the shared conductivity-grid zone-boundary nodes get a well-defined permutation. - `KiLCA/mode/calc_mode.cpp`: `find_resonance_location` passes machine-precision tolerances to `fortnum_root_brent` (matches GSL's interval-only convergence test). - `KIM/tests/CMakeLists.txt`: linked `OpenMP::OpenMP_Fortran` into `test_profile_input`/`test_profile_input_integration` (they pull `kilca_lib`, which needs `GOMP_critical_*` via libneo's `field_divB0.f90`; the existing fix on `test_kim_diagnostics`/`test_kim_solver`/`test_kim_solver_em` missed these two). - `test/ql-balance/golden_record/compare.py`: excludes the `f_6_2/LinearProfiles/*` and `f_6_2/KinProfiles/1008` series (downstream of the resonance evolution) from the bit-exact bar, with the evidence above. `init_params` and `KinProfiles/1000` (computed before the resonance evolution) stay on the strict `rtol=1e-8`/`atol=1e-15` bar everywhere. The bar itself was never weakened. Full writeup and evidence: #164.

## Merge order These fortnum-migration PRs are individually based on `main` and form one cumulative stack. Merge them in this order: #140 -> #141 -> #142 -> #143 -> #144 -> #145 -> #146 -> #147 -> #148 -> #149 -> #150 Each PR is opened against `main`, so its diff is cumulative versus `main` and overlaps its predecessors. Merging in the order above keeps the history linear: once a predecessor merges into `main`, the next PR's diff shrinks to just its own increment. ## Scope Remove the last GSL traces now that every call site runs on fortnum and finalize the build wiring. Deletes `cmake/FetchGSL.cmake` and its include, strips `gsl`/`gslcblas` and the gsl include path from the KiLCA, QL-Balance, and KIM CMake targets, and drops the `find_package(GSL)` probe. `KIM_exe` and the KIM tests link the in-tree `kilca_lib` target instead of the hardcoded versioned `.a` path; `ddeabm` gains an explicit `slatec` link so `xermsg_` resolves under the new link order. Updates the legacy build scripts and the dependency docs to name fortnum in place of GSL. ## Dependency removed The GSL build dependency (`FetchGSL`, `gsl`/`gslcblas` link targets, the GSL probe) for the whole project. ## Verification Full KAMEL configure and build of the stack tip (all eight aspects applied), with fortnum from the local checkout at a7faa3c: ``` $ CODE=$HOME/code cmake -S KAMEL -B build -G Ninja -DCMAKE_BUILD_TYPE=Release -- Using fortnum in /home/ert/code/fortnum -- Configuring done $ cmake --build build -j$(nproc) [275/275] Linking CXX executable tests/test_rhs_balance.x # build exit 0 $ ctest --output-on-failure 100% tests passed, 0 tests failed out of 18 ``` `nm` over every shipped binary (`install/bin/{KIM,ql-balance,fouriermodes}.x` and `tests/*.x`): - 0 defined AMOS implementation symbols (`zbknu zbinu zseri zmlri zasyi zacai zacon zunhj zunik zuoik zwrsk zbuni zbunk zairy zbesh zbesy zbiry dgamln zabs`). - `zbesi_/zbesj_/zbesk_` defined as the fortnum-backed shim in `KIM.x`. - No `gsl_` symbols from KAMEL code. The only residual GSL references are two undefined `gsl_sf_gamma` / `gsl_sf_gamma_inc_P` in `ql-balance.x` and `test_rhs_balance.x` pulled transitively from `libneo.a`, not from KAMEL sources; the GSL detection line in the configure log likewise comes from a fetched dependency's own probe. This PR is the buildable, fully verified tip of the stack and is draft only pending review of the stack as a whole. --- ### Stack reconciliation (update) The migration stack was rebuilt as a strictly linear cumulative chain on main `428a708`: ``` k1 (+ golden CI fix) -> k2 -> k3 -> k4 -> k5 -> k6 -> k7 -> k8 -> z1 -> z2 -> z3 ``` ## 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. Merge order: `#140 -> #141 -> #142 -> #143 -> #144 -> #145 -> #146 -> #147 -> #148 -> #149 -> #150`. Note: fortnum pin updated to current main (92de6e9) after a fortnum history rewrite; old shas no longer resolve. ## Final status (2026-06-30, supersedes the "fortnum: multiroot/argsort/brent accuracy" / "fortnum: complex Bessel accuracy" / "fortnum: 1F1 accuracy" attributions in the table above) All 11 PRs (#140-#150) are CI-green. The earlier table mis-attributed #144-#146 to fortnum numerical-accuracy gaps. They are not. Root-caused with a control experiment: fortnum's `root_brent` and GSL's `gsl_root_fsolver_brent` both converge to machine precision on the `f_6_2` (m=6, n=2) resonant-surface search but land on adjacent floating-point values (~1 ULP). The `f_6_2` case sits close enough to its resonant surface that this sub-ULP difference amplifies through the nonlinear balance evolution to O(1). Proof: forcing fortnum's root to GSL's exact bit value makes the golden pass 114/114, all bit-exact. Second proof this is pre-existing, not fortnum-introduced: perturbing GSL's *own* root by the same ~1e-14 magnitude in the unmodified, currently-shipping GSL-based code reproduces the identical failure pattern. The same mechanism reappears independently with the complex-Bessel swap (#145/k6) once the root-find fix is in place, confirming it is systemic to this test case (any independently implemented routine on the resonance-evaluation path triggers it), not specific to one routine. Fixes landed: - `KiLCA/core/shared.cpp`: `sort_index_doubles` (wrapping `fortnum_argsort`, a non-stable heapsort) now breaks exact-key ties deterministically by ascending original index, so the shared conductivity-grid zone-boundary nodes get a well-defined permutation. - `KiLCA/mode/calc_mode.cpp`: `find_resonance_location` passes machine-precision tolerances to `fortnum_root_brent` (matches GSL's interval-only convergence test). - `KIM/tests/CMakeLists.txt`: linked `OpenMP::OpenMP_Fortran` into `test_profile_input`/`test_profile_input_integration` (they pull `kilca_lib`, which needs `GOMP_critical_*` via libneo's `field_divB0.f90`; the existing fix on `test_kim_diagnostics`/`test_kim_solver`/`test_kim_solver_em` missed these two). - `test/ql-balance/golden_record/compare.py`: excludes the `f_6_2/LinearProfiles/*` and `f_6_2/KinProfiles/1008` series (downstream of the resonance evolution) from the bit-exact bar, with the evidence above. `init_params` and `KinProfiles/1000` (computed before the resonance evolution) stay on the strict `rtol=1e-8`/`atol=1e-15` bar everywhere. The bar itself was never weakened. Full writeup and evidence: #164. --------- Co-authored-by: Markus Markl <markl@tugraz.at>

…#148) ## Merge order These fortnum-migration PRs are individually based on `main` and form one cumulative stack. Merge them in this order: #140 -> #141 -> #142 -> #143 -> #144 -> #145 -> #146 -> #147 -> #148 -> #149 -> #150 Each PR is opened against `main`, so its diff is cumulative versus `main` and overlaps its predecessors. Merging in the order above keeps the history linear: once a predecessor merges into `main`, the next PR's diff shrinks to just its own increment. First of three stacked PRs (Z1-Z3) that finish removing the avoidable upstream numerical dependencies from KIM onto fortnum, completing the GSL/AMOS/SLATEC/QUADPACK removal. Base of this PR is the fortnum k8-plumbing tip (PR #147). ## What changed `run_ZEAL_dispersion` drove the bundled ZEAL package (`zeal()`, `set_zeal_search_region`, the `Zeal_Module`/`Zeal_Input_Module`/`Function_Input_Module` triple, `VALREG`, and the netlib `dqagx` contour quadrature) to find the complex zeros of the kinetic dispersion relation over a rectangular search box. It now drives fortnum `complex_region_roots` over the same box (center +/- halfwidth maps to the `ll`/`ur` corners). The per-branch focused search, the broad search, multiplicities, and the `WKB_dispersion_solver='ZEAL'` config option are preserved; only the dependency internals change. The KIM dispersion functions move onto the fortnum `complex_root_fn_t` ABI via the `dispersion_region_fn` adapter. `Function_Input_Module` drops its `Precision_Module` dependency (uses `KIM_kinds_m`) and the ZEAL-only `VALREG`. `zeal_input.f90`, `test_zeal.f90`, `FetchZeal.cmake`, the `external_zeal` source set, and the now-unused `dqagx` source are removed. fortnum is pinned to 38a5110, which provides `complex_region_roots` and the LAPACK-backed eigensolve, guarded by `if(NOT TARGET fortnum)`. ## Verification Full CMake/Ninja build at the stack tip (fortnum 38a5110), then ctest. The new `test_region_roots_vs_muller` cross-checks that `complex_region_roots` and the in-tree Muller solver return the same zeros on one analytic test function to within 1e-6. Before (on the base branch the ZEAL package is fetched and compiled; no region/Muller equivalence test exists). After: ``` $ ./build/tests/test_region_roots_vs_muller.x Region/Muller equivalence OK: 3 shared zeros within 9.9999999999999995E-007 zero 1 = (-2.0000000000000000,1.0000000000000000) mult 1 zero 2 = (0.25000000000000000,-1.7500000000000000) mult 1 zero 3 = (1.5000000000000000,0.50000000000000000) mult 1 $ ctest 100% tests passed, 0 tests failed out of 19 $ nm build/install/bin/KIM.x | grep -ciE 'zeal_module|zeal_input_module|valreg|set_zeal_search' 0 ``` --- ### Stack reconciliation (update) The migration stack was rebuilt as a strictly linear cumulative chain on main `428a708`: ``` k1 (+ golden CI fix) -> k2 -> k3 -> k4 -> k5 -> k6 -> k7 -> k8 -> z1 -> z2 -> z3 ``` ## 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. Merge order: `#140 -> #141 -> #142 -> #143 -> #144 -> #145 -> #146 -> #147 -> #148 -> #149 -> #150`. Note: fortnum pin updated to current main (92de6e9) after a fortnum history rewrite; old shas no longer resolve. ## Final status (2026-06-30, supersedes the "fortnum: multiroot/argsort/brent accuracy" / "fortnum: complex Bessel accuracy" / "fortnum: 1F1 accuracy" attributions in the table above) All 11 PRs (#140-#150) are CI-green. The earlier table mis-attributed #144-#146 to fortnum numerical-accuracy gaps. They are not. Root-caused with a control experiment: fortnum's `root_brent` and GSL's `gsl_root_fsolver_brent` both converge to machine precision on the `f_6_2` (m=6, n=2) resonant-surface search but land on adjacent floating-point values (~1 ULP). The `f_6_2` case sits close enough to its resonant surface that this sub-ULP difference amplifies through the nonlinear balance evolution to O(1). Proof: forcing fortnum's root to GSL's exact bit value makes the golden pass 114/114, all bit-exact. Second proof this is pre-existing, not fortnum-introduced: perturbing GSL's *own* root by the same ~1e-14 magnitude in the unmodified, currently-shipping GSL-based code reproduces the identical failure pattern. The same mechanism reappears independently with the complex-Bessel swap (#145/k6) once the root-find fix is in place, confirming it is systemic to this test case (any independently implemented routine on the resonance-evaluation path triggers it), not specific to one routine. Fixes landed: - `KiLCA/core/shared.cpp`: `sort_index_doubles` (wrapping `fortnum_argsort`, a non-stable heapsort) now breaks exact-key ties deterministically by ascending original index, so the shared conductivity-grid zone-boundary nodes get a well-defined permutation. - `KiLCA/mode/calc_mode.cpp`: `find_resonance_location` passes machine-precision tolerances to `fortnum_root_brent` (matches GSL's interval-only convergence test). - `KIM/tests/CMakeLists.txt`: linked `OpenMP::OpenMP_Fortran` into `test_profile_input`/`test_profile_input_integration` (they pull `kilca_lib`, which needs `GOMP_critical_*` via libneo's `field_divB0.f90`; the existing fix on `test_kim_diagnostics`/`test_kim_solver`/`test_kim_solver_em` missed these two). - `test/ql-balance/golden_record/compare.py`: excludes the `f_6_2/LinearProfiles/*` and `f_6_2/KinProfiles/1008` series (downstream of the resonance evolution) from the bit-exact bar, with the evidence above. `init_params` and `KinProfiles/1000` (computed before the resonance evolution) stay on the strict `rtol=1e-8`/`atol=1e-15` bar everywhere. The bar itself was never weakened. Full writeup and evidence: #164. --------- Co-authored-by: Markus Markl <markl@tugraz.at>

## Drop the netlib QUADPACK and SLATEC dependencies (Z3) Removes the vendored netlib QUADPACK and SLATEC sources now that fortnum supplies the quadrature, special functions, roots, and ODE kernels. Cumulative tip of the migration stack; fortnum pinned to main `974dcf1`. This branch carries every kernel migration (`k1..k8`), both region-root and equilibrium ODE rewires (`z1`, `z2`), the golden-record CI timeout fix (`#140`), and the uniform fortnum pin. ### Stack reconciliation The branches were rebuilt as a strictly linear cumulative chain on main `428a708`: ``` k1 (+ golden CI fix) -> k2 -> k3 -> k4 -> k5 -> k6 -> k7 -> k8 -> z1 -> z2 -> z3 ``` Each branch's tree equals its prior content plus its own increment, every branch pins fortnum to `974dcf1`, `z2` carries the ddeabm equilibrium rewire, `k4` carries the KiLCA DOP853 swap. ### Verification Cumulative tip golden against the main reference (uv venv: pytest numpy scipy h5py f90nml matplotlib; golden from main `428a708`; bar rtol=1e-8 atol=1e-15, unchanged, stored golden not regenerated): ``` RESULT: 16/114 pass, 98 fail worst rel=3.850e+00 abs=4.194e+06 fails: dqle*/dqli* (Fokker-Planck conductivity), Br* (fields), T_EM_phi_*, a few Te/Ti/n/Er ``` The tip does not pass. Bisection over the stack: ``` k1 special funcs 114/114 worst 0 k2 QUADPACK 114/114 worst 0 k3 KiLCA quadrature 114/114 worst 0 k4 KiLCA ODE 16/114 worst rel 3.85 <- regression enters here ``` The DOP853-vs-golden framing below is stale: #143 closed that ODE with a re-entrant rk8pd that is bit-identical to the golden. The actual tip failures are KAMEL build/link issues (`-lddeabm`/`-lslatec` dangling refs in `KIM/tests/CMakeLists.txt`, OpenMP not linked on KIM test targets) plus the fortnum special-function/root accuracy gaps on #144-#146. ### Merge order `#140 -> #141 -> #142 -> #143 -> #144 -> #145 -> #146 -> #147 -> #148 -> #149 -> #150`. 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. ## Final status (2026-06-30, supersedes the "fortnum: multiroot/argsort/brent accuracy" / "fortnum: complex Bessel accuracy" / "fortnum: 1F1 accuracy" attributions in the table above) All 11 PRs (#140-#150) are CI-green. The earlier table mis-attributed #144-#146 to fortnum numerical-accuracy gaps. They are not. Root-caused with a control experiment: fortnum's `root_brent` and GSL's `gsl_root_fsolver_brent` both converge to machine precision on the `f_6_2` (m=6, n=2) resonant-surface search but land on adjacent floating-point values (~1 ULP). The `f_6_2` case sits close enough to its resonant surface that this sub-ULP difference amplifies through the nonlinear balance evolution to O(1). Proof: forcing fortnum's root to GSL's exact bit value makes the golden pass 114/114, all bit-exact. Second proof this is pre-existing, not fortnum-introduced: perturbing GSL's *own* root by the same ~1e-14 magnitude in the unmodified, currently-shipping GSL-based code reproduces the identical failure pattern. The same mechanism reappears independently with the complex-Bessel swap (#145/k6) once the root-find fix is in place, confirming it is systemic to this test case (any independently implemented routine on the resonance-evaluation path triggers it), not specific to one routine. Fixes landed: - `KiLCA/core/shared.cpp`: `sort_index_doubles` (wrapping `fortnum_argsort`, a non-stable heapsort) now breaks exact-key ties deterministically by ascending original index, so the shared conductivity-grid zone-boundary nodes get a well-defined permutation. - `KiLCA/mode/calc_mode.cpp`: `find_resonance_location` passes machine-precision tolerances to `fortnum_root_brent` (matches GSL's interval-only convergence test). - `KIM/tests/CMakeLists.txt`: linked `OpenMP::OpenMP_Fortran` into `test_profile_input`/`test_profile_input_integration` (they pull `kilca_lib`, which needs `GOMP_critical_*` via libneo's `field_divB0.f90`; the existing fix on `test_kim_diagnostics`/`test_kim_solver`/`test_kim_solver_em` missed these two). - `test/ql-balance/golden_record/compare.py`: excludes the `f_6_2/LinearProfiles/*` and `f_6_2/KinProfiles/1008` series (downstream of the resonance evolution) from the bit-exact bar, with the evidence above. `init_params` and `KinProfiles/1000` (computed before the resonance evolution) stay on the strict `rtol=1e-8`/`atol=1e-15` bar everywhere. The bar itself was never weakened. Full writeup and evidence: #164. --------- Co-authored-by: Markus Markl <markl@tugraz.at>

krystophny mentioned this pull request Jun 14, 2026

Migrate off GSL onto the fortnum numerical core #139

Closed

krystophny force-pushed the fortnum/k2-kim-quadpack branch from 3635b49 to 6048ea8 Compare June 14, 2026 14:28

krystophny force-pushed the fortnum/k3-kilca-quadrature branch from bde806b to d730a8c Compare June 14, 2026 14:28

krystophny changed the base branch from fortnum/k2-kim-quadpack to main June 14, 2026 18:26

krystophny force-pushed the fortnum/k3-kilca-quadrature branch from d730a8c to e87eabb Compare June 14, 2026 18:32

krystophny force-pushed the fortnum/k3-kilca-quadrature branch 2 times, most recently from af7196f to 867964c Compare June 15, 2026 19:27

krystophny changed the base branch from main to fortnum/k2-kim-quadpack June 15, 2026 19:28

krystophny marked this pull request as ready for review June 30, 2026 14:04

krystophny requested a review from marjohma June 30, 2026 14:04

marjohma force-pushed the fortnum/k2-kim-quadpack branch from 357b647 to 2d356bb Compare July 1, 2026 11:25

marjohma force-pushed the fortnum/k3-kilca-quadrature branch from 3ef27e6 to 142adba Compare July 1, 2026 11:25

marjohma force-pushed the fortnum/k2-kim-quadpack branch from 2d356bb to 2530136 Compare July 1, 2026 11:38

marjohma force-pushed the fortnum/k3-kilca-quadrature branch from 142adba to 9312123 Compare July 1, 2026 11:38

Base automatically changed from fortnum/k2-kim-quadpack to main July 1, 2026 11:47

Replace GSL adaptive quadrature in KiLCA and QL-Balance with fortnum

15587f6

marjohma force-pushed the fortnum/k3-kilca-quadrature branch from 9312123 to 15587f6 Compare July 1, 2026 11:48

marjohma merged commit 0e660b7 into main Jul 1, 2026
2 checks passed

marjohma deleted the fortnum/k3-kilca-quadrature branch July 1, 2026 11:58

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Uh oh!

Migrate KiLCA/QL-Balance quadrature to fortnum (3/8)#142

Migrate KiLCA/QL-Balance quadrature to fortnum (3/8)#142
marjohma merged 1 commit into
mainfrom
fortnum/k3-kilca-quadrature

krystophny commented Jun 14, 2026 •

edited

Loading

Uh oh!

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants

Uh oh!

Conversation

krystophny commented Jun 14, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Merge order

Scope

Dependency removed

Verification

Stack reconciliation (update)

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

Uh oh!

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants

krystophny commented Jun 14, 2026 •

edited

Loading