Skip to content

Make KiLCA background ODE reproduce the GSL-recorded golden bit-for-bit#153

Closed
krystophny wants to merge 10 commits into
mainfrom
fortnum/k4-imhd-rk8pd-fix
Closed

Make KiLCA background ODE reproduce the GSL-recorded golden bit-for-bit#153
krystophny wants to merge 10 commits into
mainfrom
fortnum/k4-imhd-rk8pd-fix

Conversation

@krystophny

@krystophny krystophny commented Jun 15, 2026

Copy link
Copy Markdown
Member

Summary

The QL-Balance golden record (ql-balance_golden_record) compares the
fortnum-backed build against the GSL-era golden at rtol=1e-8/atol=1e-15. The
KiLCA linear-response solve amplifies a ~1-2 ULP background-field difference by
~1e7, so dqle*/dqli*/Br_* diverged from the golden by ~2.7e-6 and the test
failed at every time step.

Root cause: two fortnum kernels did not reproduce their GSL predecessors.

  1. The background equilibrium ODE in KiLCA/background/calc_back.cpp is driven
    through the old gsl_odeiv API (gsl_odeiv_evolve_apply +
    gsl_odeiv_control_y_new). fortnum's rk8pd evolve loop modelled the odeiv2
    semantics: it scaled the failure ratio by the pre-step solution and rejected
    every step with rmax > 1.1 unconditionally. The old API updates y in place
    before the controller runs and keeps an otherwise-rejected step when the
    suggested decrease would not move t by at least one ULP. That mismatch
    produced a different accepted-step mesh and a ULP-level background profile.

  2. fortnum's deriv_central used a different 5-point stencil than
    gsl_deriv_central and skipped GSL's optimal-stepsize second pass, diverging
    by up to ~1e-10 on smooth functions.

Both are fixed in fortnum; this PR only bumps the fortnum pin. No golden
regeneration, no tolerance change, no skipped cases.

Verification

Before (fortnum at the previous pin):

test_golden_record.py::test_golden_record FAILED
[FAIL] /f_6_2/LinearProfiles/0/Br_Re   max rel 2.722e-06
[FAIL] /f_6_2/LinearProfiles/0/dqle11  max rel 2.705e-06
... 100+ failing quantities across all time steps

After (fortnum at the bumped pin):

============================= 115 passed in 33.57s =============================
worst_abs=0.000e+00 worst_rel=0.000e+00   (cur output bit-identical to golden)

Clean-room check, old gsl_odeiv rk8pd + control_y_new(1e-16,1e-16) vs
fortnum on a stiff scalar ODE carried across a grid: max relative gap 1e-10
before, 0.0 after.

gsl_deriv_central vs fortnum_deriv_central on smooth functions: max
relative gap 2e-10 before, 0.0 after.

Remaining KAMEL ctest cases (special functions, QUADPACK QAGI, integration
methods, KIM adapter, sparse, ampere matrices): 19/19 pass. fortnum's own suite:
74/74 pass.

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

calc_back's background equilibrium ODE was tuned against GSL
gsl_odeiv_step_rk8pd under gsl_odeiv_control_y_new, run as one continuous
adaptive evolve carrying the step across the radial grid. The per-segment
dop853 path reset the integrator at every grid point and uses a different 8(7)
error norm. Replace it with fortnum's GSL-faithful rk8pd C ABI
(fortnum_rk8pd_create / integrate_to / destroy), evolving continuously with the
golden initial step h0 = x[1]-x[0]; the background profile is now bit-identical
to the golden's GSL rk8pd output. incompressible.cpp / compressible_flow.cpp
stay on dop853.

Re-pin fortnum to the main commit carrying rk8pd and its C ABI.
The incompressible and compressible-flow imhd zone basis solvers built their
radial grid from the original gsl_odeiv_evolve_apply loop under
gsl_odeiv_control_y_new(eps_abs, eps_rel): one continuous rk8pd evolve that
stored every accepted step. The dop853 replacement reset per segment, used a
different 8(7) error norm, and produced a different accepted-step mesh.

Evolve continuously with fortnum's rk8pd C ABI, recording each accepted step
via fortnum_rk8pd_step_to (one gsl_odeiv_evolve_apply per call), starting from
h0 = 1e-6 in the integration direction as the original did. Re-pin fortnum to
the main commit carrying the single-step rk8pd ABI and the GSL std_control /
yerr fidelity fixes.
Pin fortnum to the commit whose rk8pd evolve mirrors the old gsl_odeiv
gsl_odeiv_evolve_apply control flow and whose deriv_central matches
gsl_deriv_central. The KiLCA background equilibrium ODE now reproduces the
GSL-recorded background bit-for-bit, so ql-balance_golden_record passes
115/115 at rtol=1e-8/atol=1e-15 with the cur output bit-identical to the
golden.
ql-balance golden stays bit-identical (115/115, worst rel/abs 0) at the
unchanged rtol=1e-8, atol=1e-15.
@krystophny

Copy link
Copy Markdown
Member Author

Superseded by #143. The winning state from this branch (calc_back evolved as one continuous fortnum rk8pd run via the C ABI, plus the imhd incompressible and compressible-flow zones moved from dop853 to the same continuous rk8pd) has been consolidated into the single k4 PR #143. #143 now carries:

  • Replace GSL rk8pd ODE evolve in KiLCA with fortnum dop853 (the base step)
  • Evolve KiLCA background ODE continuously with fortnum rk8pd (calc_back)
  • Evolve KiLCA imhd zones with GSL-faithful rk8pd instead of dop853

The downstream stack (#144..#150) is rebased on top of this consolidated k4, so the continuous-rk8pd background and imhd state flows through the whole stack. The k4 tip reproduces the ql-balance golden record 115/115 (rtol=1e-8, atol=1e-15). fortnum is pinned at 92de6e9 on every branch. Closing in favor of #143.

@krystophny krystophny closed this Jun 15, 2026
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.

1 participant