Drop GSL and finalize fortnum wiring (8/8)#147
Conversation
|
Stacked on this branch, PRs #148-#150 finish the upstream numerical dependency removal that this PR (AMOS) was part of, completing the GSL/AMOS/SLATEC/QUADPACK drop. Only BLAS/LAPACK/OpenMP/MPI remain.
nm on every shipped binary: 0 dqag/quadpack-netlib, 0 d1mach/i1mach/xermsg/slatec, 0 zeal, 0 ddeabm. ctest green at each step. The only residual undefined GSL symbol is |
b41221c to
eba278e
Compare
7a3574a to
7a8cbba
Compare
7a8cbba to
38393a4
Compare
2ac023b to
a17d191
Compare
0e5ce8d to
a08c0c4
Compare
701e07c to
0df5ceb
Compare
## 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.
a08c0c4 to
080be51
Compare
0df5ceb to
e19a5c4
Compare
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.
080be51 to
f3c6dca
Compare
e19a5c4 to
3503265
Compare
## 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.
f3c6dca to
a0570a5
Compare
3503265 to
7c25073
Compare
… (4/8) (#143) 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.
a0570a5 to
b1716ee
Compare
7c25073 to
eef5757
Compare
## 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.
7ca1846 to
c65d98e
Compare
## 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.
These targets link kilca_lib, which pulls in libneo's field_divB0.f90 (OMP critical sections). CI link error: undefined reference to GOMP_critical_start/end. test_kim_diagnostics/solver/solver_em already link OpenMP::OpenMP_Fortran for the same reason; these two were missed.
a8ce0cf to
3c4945b
Compare
The GOMP_critical_* symbols come from libneo's field_divB0.f90, which is pulled into any executable that links kilca_lib. Until now each consumer that happened to drag in that object had to link OpenMP::OpenMP_Fortran itself; the KIM test targets were fixed in two separate passes as the missing ones surfaced, and five other kilca_lib consumers still omit it (they only build because they don't currently reference that object). Make OpenMP part of kilca_lib's PUBLIC link interface so every consumer (KIM_exe, ql-balance_lib, all KIM tests) inherits it automatically, which removes the per-target annotations and the latent breakage if one of the un-annotated targets later pulls the OMP object. To reference the imported target from kilca_lib's interface it must exist on all platforms, so hoist find_package(OpenMP) out of the APPLE-only branch to top scope (it was previously never resolved on Linux, where the top level only appends -fopenmp to the compile flags). The macOS libomp discovery and per-platform flag logic are unchanged.
There was a problem hiding this comment.
ORCHESTRATE_REVIEW_STATUS: approve
This PR cleanly removes the GSL build dependency and finalizes the fortnum/OpenMP build wiring across the whole KAMEL project. The changes are minimal, correct, and well-scoped. CI is green: build-and-test passes 24/24 tests, golden-record (main) passes (58 files, 0 over rtol=1e-7), and golden-record (pre-fortnum diagnostic) fails as expected (diagnostic-only, continue-on-error). I verified the build log: no undefined references, no link errors, no GSL symbols in the final binaries.
Key correctness points verified:
find_package(OpenMP REQUIRED COMPONENTS C Fortran)moved to top scope (CMakeLists.txt:60) —OpenMP::OpenMP_Fortranis now visible to all subdirectories on both macOS and Linux. Previously it was only called insideif(APPLE), so Linux builds had no imported target.kilca_libnow PUBLIC-linksOpenMP::OpenMP_Fortran(KiLCA/CMakeLists.txt:256), carryingGOMP_critical_*resolution transitively to KIM_exe, ql-balance_lib, and all KIM tests. No per-target OpenMP linkage needed.- Hardcoded
.apath (libKiLCA_Lib_V_2.4.2_..._64bit.a) replaced with thekilca_libCMake target in KIM_exe and all 10 KIM test targets. This fixes the link-order fragility and thehyper1F1.cppnot-linked error noted in the PR body. FetchGSL.cmakedeleted;include(FetchGSL)removed fromDependencies.cmake; allfind_package(GSL),GSL::gsl/GSL::gslcblas, andGSL_INCLUDE_DIRSremoved from CMake files. No remaining GSL references in any CMake build file (verified by grep).ddeabmretains itstarget_link_libraries(ddeabm PUBLIC slatec)with an updated comment —xermsg_resolution is preserved.add_subdirectoryorder (KiLCA before KIM before QL-Balance) ensureskilca_libis defined before use.
Low-severity finding (pre-existing, not in diff): .github/workflows/ci.yml:40 and .github/workflows/golden.yml:32,67 still install libgsl-dev via apt. The build no longer uses GSL, so this package is unnecessary. Not a blocker — the build passes without it — but it's a leftover that should be cleaned up in a follow-up for consistency with utility_scripts/setup_debian.sh (which this PR correctly cleans).
{"findings":[{"file":".github/workflows/ci.yml","line":40,"severity":"low","message":"CI workflow still installs libgsl-dev; no longer needed after GSL removal. Pre-existing, not introduced by this PR."},{"file":".github/workflows/golden.yml","line":32,"severity":"low","message":"CI workflow still installs libgsl-dev; no longer needed after GSL removal. Pre-existing, not introduced by this PR."}]}The build no longer links GSL after the fortnum migration, so the apt install of libgsl-dev in ci.yml and golden.yml is dead weight.
The golden-record job rebuilds the golden-baseline (and pre-fortnum) reference trees, which predate the GSL removal and still call find_package(GSL REQUIRED) in KIM/src/CMakeLists.txt. Dropping libgsl-dev in d1c9b5d broke the reference build's CMake configure: Could NOT find GSL (missing: GSL_INCLUDE_DIR ...) at KIM/src/CMakeLists.txt HEAD no longer needs GSL, so ci.yml correctly drops it, but golden.yml must keep it to compile the frozen references it diffs against.
…#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>
…deabm (Z2) (#149) ## Integrate the equilibrium ODE with fortnum ddeabm (Z2) `KIM/src/background_equilibrium/calculate_equil.f90` now integrates the scalar force-balance ODE `du/dr` with the fortnum variable-order Adams integrator (`fortnum_ode_ddeabm`), a clean-room equivalent of SLATEC ddeabm. This replaces the earlier DOP853 substitution and restores the original SLATEC continuation semantics. One `ddeabm_init` seeds the state at `r_grid(1)`. Each outer grid value `u(i)` continues the same re-entrant state to `r_grid(i)` (the SLATEC `INFO(1)=1` restart), with the step barred from overshooting `r_grid(end)` (the `RWORK(1)` / `INFO(4)=1` tstop bound). `idid < 1` maps to the existing `log_warning` path via `ode_status%code /= FORTNUM_OK`. The original scalar RHS `dudr(r, u, du)` is wrapped as the array-valued `ode_rhs_t` callback `dudr(r, u, du, ctx)` with `neq = 1`; the q-profile and pressure data ride on host association, so `ctx` is unused. Cumulative on `#148 z1`; fortnum pinned to main `974dcf1`. ### Verification This rewire compiles and runs. Its own change is sound: on the cluster the KIM.x `electromagnetic`/`quadpack` cases (which exercise the ddeabm -> integrator path) are bit-identical. The CI failure on this branch is a KAMEL build issue - `-lddeabm` is still referenced in `KIM/tests/CMakeLists.txt` after the library was dropped. Tolerances unchanged (rtol=1e-12, atol=1e-12 for this ODE; golden bar rtol=1e-8 atol=1e-15). 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.
## 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>
Merge order
These fortnum-migration PRs are individually based on
mainand form onecumulative 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 versusmainandoverlaps 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 itsown increment.
Scope
Remove the last GSL traces now that every call site runs on fortnum and finalize the build wiring. Deletes
cmake/FetchGSL.cmakeand its include, stripsgsl/gslcblasand the gsl include path from the KiLCA, QL-Balance, and KIM CMake targets, and drops thefind_package(GSL)probe.KIM_exeand the KIM tests link the in-treekilca_libtarget instead of the hardcoded versioned.apath;ddeabmgains an explicitslateclink soxermsg_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/gslcblaslink 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:
nmover every shipped binary (install/bin/{KIM,ql-balance,fouriermodes}.xandtests/*.x):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 inKIM.x.gsl_symbols from KAMEL code. The only residual GSL references are two undefinedgsl_sf_gamma/gsl_sf_gamma_inc_Pinql-balance.xandtest_rhs_balance.xpulled transitively fromlibneo.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: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):
KiLCA/math/hyper/hyper1F1.cppobject not linked into QL-Balance/test targets-lddeabmstill referenced inKIM/tests/CMakeLists.txt-lddeabm/-lslatecdangling refs + OpenMP (GOMP_critical_*) not linked on KIM test targetsThe 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_brentand GSL'sgsl_root_fsolver_brentboth converge to machine precision on thef_6_2(m=6, n=2) resonant-surface search but land on adjacent floating-point values (~1 ULP). Thef_6_2case 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(wrappingfortnum_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_locationpasses machine-precision tolerances tofortnum_root_brent(matches GSL's interval-only convergence test).KIM/tests/CMakeLists.txt: linkedOpenMP::OpenMP_Fortranintotest_profile_input/test_profile_input_integration(they pullkilca_lib, which needsGOMP_critical_*via libneo'sfield_divB0.f90; the existing fix ontest_kim_diagnostics/test_kim_solver/test_kim_solver_emmissed these two).test/ql-balance/golden_record/compare.py: excludes thef_6_2/LinearProfiles/*andf_6_2/KinProfiles/1008series (downstream of the resonance evolution) from the bit-exact bar, with the evidence above.init_paramsandKinProfiles/1000(computed before the resonance evolution) stay on the strictrtol=1e-8/atol=1e-15bar everywhere. The bar itself was never weakened.Full writeup and evidence: #164.