Skip to content

Integrate the equilibrium ODE with fortnum ddeabm instead of SLATEC ddeabm (Z2)#149

Merged
marjohma merged 3 commits into
mainfrom
fortnum/z2-ddeabm-dop853
Jul 10, 2026
Merged

Integrate the equilibrium ODE with fortnum ddeabm instead of SLATEC ddeabm (Z2)#149
marjohma merged 3 commits into
mainfrom
fortnum/z2-ddeabm-dop853

Conversation

@krystophny

@krystophny krystophny commented Jun 14, 2026

Copy link
Copy Markdown
Member

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.

@krystophny krystophny force-pushed the fortnum/z1-zeal-region-roots branch from 756729d to 7af37fd Compare June 14, 2026 14:29
@krystophny krystophny force-pushed the fortnum/z2-ddeabm-dop853 branch from 67351be to e605baa Compare June 14, 2026 14:29
@krystophny krystophny changed the base branch from fortnum/z1-zeal-region-roots to main June 14, 2026 18:26
@krystophny krystophny force-pushed the fortnum/z2-ddeabm-dop853 branch from e605baa to 124517b Compare June 14, 2026 18:32
@krystophny krystophny force-pushed the fortnum/z2-ddeabm-dop853 branch 2 times, most recently from 1a8f62f to 72c8dab Compare June 15, 2026 19:28
@krystophny krystophny changed the base branch from main to fortnum/z1-zeal-region-roots June 15, 2026 19:29
@krystophny krystophny requested a review from marjohma June 30, 2026 19:21
@marjohma marjohma force-pushed the fortnum/z1-zeal-region-roots branch from ad14235 to c6f465f Compare July 1, 2026 11:26
@marjohma marjohma force-pushed the fortnum/z2-ddeabm-dop853 branch from c1c3c44 to 3a116d3 Compare July 1, 2026 11:26
marjohma pushed a commit that referenced this pull request Jul 1, 2026
## 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.
@marjohma marjohma force-pushed the fortnum/z1-zeal-region-roots branch from c6f465f to 7b04350 Compare July 1, 2026 11:38
@marjohma marjohma force-pushed the fortnum/z2-ddeabm-dop853 branch from 3a116d3 to e4b430d Compare July 1, 2026 11:38
marjohma pushed a commit that referenced this pull request Jul 1, 2026
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.
@marjohma marjohma force-pushed the fortnum/z1-zeal-region-roots branch from 7b04350 to be942c9 Compare July 1, 2026 11:48
@marjohma marjohma force-pushed the fortnum/z2-ddeabm-dop853 branch from e4b430d to 96bf656 Compare July 1, 2026 11:48
marjohma pushed a commit that referenced this pull request Jul 1, 2026
## 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.
@marjohma marjohma force-pushed the fortnum/z1-zeal-region-roots branch from be942c9 to 9127aa9 Compare July 1, 2026 11:59
@marjohma marjohma force-pushed the fortnum/z2-ddeabm-dop853 branch from 96bf656 to 23f1387 Compare July 1, 2026 11:59
marjohma pushed a commit that referenced this pull request Jul 1, 2026
… (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.
@marjohma marjohma force-pushed the fortnum/z1-zeal-region-roots branch from 9127aa9 to ea914dd Compare July 1, 2026 12:10
@marjohma marjohma force-pushed the fortnum/z2-ddeabm-dop853 branch from 23f1387 to c02f6c0 Compare July 1, 2026 12:10
marjohma added a commit that referenced this pull request Jul 2, 2026
## 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>
krystophny added a commit that referenced this pull request Jul 2, 2026
## 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.
krystophny added a commit that referenced this pull request Jul 2, 2026
## 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.

@slopqueue slopqueue Bot left a comment

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

ORCHESTRATE_REVIEW_STATUS: request_changes

Here's my complete review of PR #149.


PR #149 Review: Replace SLATEC ddeabm with fortnum ddeabm in KIM

Build & Test Results

Check Result
Full build (macOS, gfortran) Pass
test_equil_ode_dop853 Pass (rel err ~4.5e-16)
Full test suite (26 tests) 25 pass, 1 fail
Failed test: test_rhs_balance Pre-existing QL-Balance bug (forces_nl stale value, documented in AGENTS.md) — NOT related to this PR

Production Code: calculate_equil.f90 — Correct

The swap from SLATEC ddeabm to fortnum ddeabm is correctly implemented:

  • ddeabm_init at r_grid(1) replaces SLATEC INFO(1)=0 (fresh start)
  • ddeabm_integrate_to with carried ode_state replaces SLATEC INFO(1)=1 (continuation)
  • tstop=rstop replaces SLATEC INFO(4)=1, RWORK(1)=tstop (boundary bound)
  • u_seg(1) extracts solution at target — replaces old u0 (which SLATEC updated in-place)
  • Error handling improved: on failure, holds u(i-1) instead of using potentially stale u0
  • dudr interface correctly changed from scalar to array to match ode_rhs_t

Blocking Issue: Test Does Not Test Production Code

test_equil_ode_dop853.f90 uses DOP853; production uses ddeabm. The test:

  • Calls ode_solve_dop (DOP853 integrator from fortnum_ode_dop853)
  • Does NOT call calculate_equil or use ddeabm_init/ddeabm_integrate_to
  • Reimplements the ODE with a different integrator on a simplified problem (constant q, dpress=0)
  • Test comment says "ddeabm -> fortnum DOP853 swap" but production actually swaps to fortnum ddeabm, not DOP853

The test passes but validates the wrong integrator. The changed production code path (fortnum ddeabm) has no test coverage.

Minor Issues

  1. Stale comment (common/math/CMakeLists.txt:6): "stay for QUADPACK and ddeabm" — ddeabm is now from fortnum. Pre-existing (introduced by earlier merged PR #145), not introduced here.

  2. KiLCA ddeabm references: calc_back_slatec.cpp and calc_I_array_drift.f90 call ddeabm_ but are dead code (not in CMakeLists.txt source list). Not a breakage.

Verdict: Request Changes

The production code change is correct and well-implemented. The test must be fixed to exercise the actual ddeabm path — either call calculate_equil directly, use ddeabm_init+ddeabm_integrate_to, or at minimum test the same integrator as production. Once the test matches the production code path, this PR is approvable.

@krystophny krystophny force-pushed the fortnum/z1-zeal-region-roots branch from 5685753 to 9fa2299 Compare July 3, 2026 05:22
@krystophny krystophny force-pushed the fortnum/z2-ddeabm-dop853 branch from dfbb4e7 to 7eb83ff Compare July 3, 2026 05:23
@krystophny krystophny changed the base branch from fortnum/z1-zeal-region-roots to main July 3, 2026 05:32
@krystophny krystophny force-pushed the fortnum/z2-ddeabm-dop853 branch 2 times, most recently from ddc31c9 to f2997df Compare July 3, 2026 07:07
marjohma added a commit that referenced this pull request Jul 10, 2026
## 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>
marjohma added a commit that referenced this pull request Jul 10, 2026
…#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>
ddeabm was replaced by the fortnum DOP853 integrator; the library no
longer exists, so linking -lddeabm dangles. Remove it from the KIM
test targets, matching the other migrated targets.
@marjohma marjohma force-pushed the fortnum/z2-ddeabm-dop853 branch from 85b73a5 to 1fc92d0 Compare July 10, 2026 08:14
@marjohma marjohma changed the title Integrate the equilibrium ODE with fortnum DOP853 instead of ddeabm (Z2) Integrate the equilibrium ODE with fortnum ddeabm instead of SLATEC ddeabm (Z2) Jul 10, 2026
@marjohma marjohma merged commit eb2ed22 into main Jul 10, 2026
3 checks passed
@marjohma marjohma deleted the fortnum/z2-ddeabm-dop853 branch July 10, 2026 08:34
marjohma added a commit that referenced this pull request Jul 10, 2026
## 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>
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants