ql-balance/KiLCA golden record: run-to-run nondeterminism in the linear solve (root cause + fixes)

[krystophny]

Summary

The ql-balance / KiLCA golden-record comparison is unstable because the KiLCA linear solve is not reproducible run to run for identical inputs — even single-threaded. Two independent memory bugs in code shared by all backends are the cause; both are pre-existing on main and unrelated to the fortnum migration. This issue records the root-cause analysis and links the fixes.

Case used throughout: golden f_6_2 (m=6, n=2). Primary observable: flre/linear-data/m_6_n_2_flab_[1,0]/zone_0_poy_test_err.dat and EB.dat.

Root causes

1. Uninitialized interpolation index (dominant) — fixed in #170

interp_profile (in QL-Balance/src/base/wave_code_data_64bit.f90 and KiLCA/interface/wave_code_data_64bit.f90) passes a local integer ind to eval_neville_polynom without initializing it. ind is the in/out warm-start window index; on the first call it is read uninitialized (guard in KiLCA/interp/interp.h only catches out-of-range garbage; in-range garbage flows through, out-of-range garbage triggers an OOB read of the grid array). Because this interpolates the initial background profiles onto the balance grid, it seeds a stack-garbage-dependent perturbation that propagates through the whole solve.

• valgrind memcheck: 22 "uninitialised value" reads in eval_neville_polynom, origin interp_profile's ind.
• After fixing: those reads -> 0; run-to-run EB.dat goes from ~7e-4 relative spread to byte-identical across runs.

2. Buffer over-read on the run path — fixed in #171

calc_wave_code_data_, calc_wave_code_data_for_mode_, calc_conductivity_matrices_for_mode_ in KiLCA/interface/wave_code_interface.cpp strcpy a non-NUL-terminated Fortran character buffer of length *pathlength. strcpy over-reads past the buffer (and can overflow the fixed 1024-byte destination). Fixed by memcpy of the known length.

• valgrind: "Invalid read of size 1" in calc_wave_code_data_for_mode_ -> 0.

3. Residual nondeterminism (open)

After #170 and #171, a much smaller run-to-run effect remains: the solve settles into one of two discrete states differing by ~1e-9 in EB.dat (~1 ULP in the resonant radius r_res). It is valgrind-clean (not uninitialized memory), single-threaded, and appears in both GSL and fortnum builds. Cause not yet identified (candidate: address/ordering-dependent branch). This is the last blocker to a byte-exact golden.

poy_test_err is a diagnostic, not a physics output

zone_0_poy_test_err.dat is the Poynting energy-balance consistency residual |div(PoyntingFlux) - (absorbed_power + j.E)| — the solution's own numerical error (max ~1e-3), a difference of nearly-equal terms. It amplifies last-ULP input differences by ~1e8 (a ~1e-9 field perturbation becomes ~50% per-node). It should not be compared byte/rtol-exact across runs or backends. Recommend excluding it from the golden or flooring it at its own scale (atol ~1e-3), documented as a diagnostic.

The fortnum migration is transparent

Relevant to #144 and #164: with the solve made deterministic (threads pinned), GSL (main) and fortnum builds produce byte-identical output across the whole linear-data tree, including r_res. The apparent GSL-vs-fortnum differences reported earlier (sort tie order; r_res 1-ULP) are the run-to-run bimodality of item 3, present in the unmodified GSL code too — not a fortnum defect. #164's own control experiment (perturbing GSL's own root reproduces the failure) is consistent with this. Once items 1-3 land, the ql-balance step exclusions added for #164 should be revisited and the fortnum-regression framing corrected.

Reproduction

Build ql-balance.x, run it twice on the same f_6_2 runfolder with OMP_NUM_THREADS=1 OPENBLAS_NUM_THREADS=1, and diff flre/linear-data/m_6_n_2_flab_[1,0]/EB.dat. On main the two runs differ; with #170 applied they are byte-identical (modulo item 3). valgrind --tool=memcheck --track-origins=yes on the same run surfaces both memory errors.

Status

• Uninitialized interp index -> Initialize interpolation index in interp_profile (fixes nondeterministic solve) #170
• Run-path buffer over-read -> Copy KiLCA run path by known length instead of strcpy (fixes buffer over-read) #171
• Residual ~1e-9 bimodal nondeterminism (item 3)
• poy_test_err golden handling (exclude/floor as diagnostic)
• Revisit ql-balance golden record: near-resonant root location is ULP-sensitive (pre-existing, not a fortnum regression) #164 step exclusions and fortnum framing once the above land

@marjohma for visibility — this reframes the #144 golden failures as pre-existing KiLCA nondeterminism, not a fortnum regression.

[krystophny]

Item 3 root cause confirmed: OpenBLAS run-to-run FP nondeterminism

Ran the KiLCA field solve (ql-balance.x, f_6_2, OMP_NUM_THREADS=1) under three BLAS conditions:

• OpenBLAS: EB.dat usually one state, rarely (~1 in 40) a second state ~1.8e-8 apart.
• Reference Netlib BLAS/LAPACK (LD_PRELOAD=/usr/lib/liblapack.so.3 /usr/lib/libblas.so.3): 8/8 valid runs bit-identical (deterministic).
• OpenBLAS vs reference: EB differs by max_rel = 7.5e-8 — same order as the run-to-run flip.

So BLAS is in the numeric path, reference BLAS is reproducible, and OpenBLAS is not. Every uninitialized-memory hypothesis was excluded first:

Hypothesis	Test	Result
uninit heap	MALLOC_PERTURB_=0/1/255	identical -> not heap
ASLR	setarch -R	still flips -> not ASLR
uninit C++ local	-ftrivial-auto-var-init=pattern	unchanged -> not C++
uninit Fortran int	-finit-integer	unchanged -> not int
uninit Fortran real (arith)	-finit-real=snan	no NaN in EB -> not real-in-arith
uninit real in solve	-finit-real=snan -ffpe-trap=invalid	ran through solve, no trap
module/static vars	(BSS zero-filled)	deterministic
memcheck (post #170/#171)	valgrind	clean

Conclusion: item 3 is OpenBLAS SIMD-kernel/alignment-dependent rounding, run-to-run, ~1e-8. It is below the golden's EB rtol (1e-7), so EB/mode_data/profiles pass; only poy_test_err (which amplifies last-ULP noise ~1e8x) is broken. Not a KAMEL code defect.

Recommendation

• Run the golden with a deterministic BLAS (reference Netlib, or a pinned/deterministic OpenBLAS build) so every run is bit-reproducible; and/or
• Floor/exclude poy_test_err as a diagnostic (warranted regardless of BLAS, given the ~1e8 amplification).

Neither is a KAMEL source change. The read-before-set found during this trace is fixed in #173.

Status

• Uninitialized interp index -> Initialize interpolation index in interp_profile (fixes nondeterministic solve) #170
• Run-path buffer over-read -> Copy KiLCA run path by known length instead of strcpy (fixes buffer over-read) #171
• Item 3 root-caused: OpenBLAS FP nondeterminism (env/CI fix: deterministic BLAS, or floor poy_test_err)
• Read-before-set constant_time_step -> Initialize constant_time_step and set_constant_time_step (read-before-set) #173
• poy_test_err golden handling (floor/exclude as diagnostic) — libneo harness
• Revisit ql-balance golden record: near-resonant root location is ULP-sensitive (pre-existing, not a fortnum regression) #164 step exclusions once the above land