inaccuracy: legacy single-precision physics constants (c/pi/sol/kB) in KiLCA + KIM

[krystophny]

Category: inaccuracy (latent, pre-existing; encoded in the golden record).

Several legacy Fortran files declare real(8)/real(dp) physics constants from
literals that are parsed at single precision (no _dp/d0 suffix, or an e
exponent), so a ~7-significant-figure value is stored in a double and then used
in double-precision physics. Example: 29979245800.0 (bare) rounds to
29979246592.0 — off by ~792, i.e. ~2.6e-8 relative.

These are not bugs introduced by the C++->Fortran port; they predate it and
the KiLCA golden record was generated with them, so they must not be "fixed"
in isolation — doing so shifts the golden outputs. (The port already had to
reproduce one of these: the KiLCA constants module keeps its single-precision
c for the legacy Fortran physics, while the ported-from-C++ modules carry a
local double c mirroring the oracle's constants.h; see PR #161.)

Affected (verified legacy, 0 lines changed vs the pre-port base)

• KiLCA/hom_medium/hmedium.f90:21,23,106,108,351,352,408,409 — c = 29979245800.0, pi = 3.14159... (single).
• KiLCA/imhd/imhd.f90:137,138,221,222 — pi, c (single).
• KiLCA/math/W/wfunc.f90:11, wmfunc.f90:11, wm_asympt.f90:11 — pi = 3.141592653589793238462643383279502884197 (single).
• KIM/src/util/constants_mod.f90:8 — sol = 29979245800.0 (single; used across kernels/electromagnetic_solver/species/calculate_equil).
• KIM/src/util/constants_mod.f90:13 — kB = 1.380649e-16 (e exponent -> single: 1.3806490132968306e-16 vs 1.380649e-16).

Fix (deferred, requires golden regeneration)

Add explicit double-precision suffixes (_dp) / d-exponents to the literals,
consolidate onto a single shared constants module where possible, then
regenerate the KiLCA and KIM golden records and confirm the strict bar holds.
Coordinate with any golden exclusions.

Environment

Branch fortran-port-all (PR #161); files unchanged since the pre-port base.

[krystophny]

Automated triage pass (Claude), grounded in the repository at 2722090f. Verify before acting.

Still valid on main. Every line cited in the body is unchanged at 2722090f, and the sweep below turns up two more affected files that the issue does not list — the shared KiLCA constants module, which is where most of the KiLCA Fortran physics actually gets c from. One listed item (kB) turns out to be dead weight and is a free fix. The rest need the golden re-baseline the issue anticipates.

Purpose

KAMEL is CGS/Gaussian throughout, so c = 2.99792458e10 cm/s enters nearly every electromagnetic expression — 4*pi/c in the Maxwell stitching matrices (KiLCA/hom_medium/hmedium.f90:353, KiLCA/imhd/imhd.f90:139), 4.0d0*pi/sol in the Ampère block of KIM's EM system (KIM/src/electromagnetic/electromagnetic_solver.f90:186,192), omega_c = ZeB/(mc) in KIM/src/background_equilibrium/species_mod.f90:490-491, and om_E = -c*k_s*E_r/B_0 in KIM/src/background_equilibrium/calculate_equil.f90:148. A constant silently truncated to single precision injects a 2.6e-8 relative error — one to two orders of magnitude above the golden-record tolerance (rtol 1e-7, test/golden/bin/gr_compare.sh:2) and far above double round-off. That is a defined, non-physical bias sitting underneath every result the code produces, and it defeats any attempt to attribute residual golden drift to a real numerical cause. It also makes the source a trap: real(dp), parameter :: c = 29979245800.0 reads as exact but is not.

Reproduced with gfortran (all four rows re-run during this triage pass; see the recipe below):

literal as written	stored double	relative error
29979245800.0	2.99792465920000000E+10	+2.6418e-08
1.380649e-16	1.38064901329683061E-16	+9.6309e-09
3.141592653589793238462643383279502884197 (bare)	3.14159274101257324E+00	+2.7828e-08
3.1415926535897932384626 (bare)	3.14159274101257324E+00	+2.7828e-08

Note the last row: the number of digits written is irrelevant. Without a d exponent or _dp suffix, Fortran's default real kind applies and everything past ~7 significant figures is discarded before the widening conversion to real(8).

Current state

Confirmed affected, in the build, at 2722090f (each verified by reading the file, not by grep alone):

1. KiLCA/core/constants_m_64bit.f90:35 and KiLCA/core/constants_m_32bit.f90:35 — not listed in the issue body, and the highest-impact item:

   real(dp), parameter :: c  = 29979245800.0   !speed of light in vacuum

   Every other real constant in that module is correctly kinded (pi = 3.14159...197_dp at :20, boltz = 1.60216428d-12 at :25, mp/me/e at :36-38) — c is the lone offender. The module is used by 59 files under KiLCA/.
2. KiLCA/hom_medium/hmedium.f90:21,23 and :106,108 — c = 29979245800.0, pi = 3.1415926535897932384626; :351,352 and :408,409 — pi, c, feeding fpc = 4.0d0*pi/c at :353.
   Caveat the issue does not make: lines 21 and 23 live in eval_fields_in_hom_media (:5), and rg -n 'eval_fields_in_hom_media' -g '!build/**' . finds no caller anywhere in the repo — that subroutine is dead code inside a live file. The output-moving lines are :106,108 (eval_basis_in_hom_media, called from hmedium_zone.cpp:107,115), :351,352 (stitching_equations_hommed_hommed, called from calc_mode.cpp:433) and :408,409 (calc_coeffs_vacuum_antenna_vacuum).
3. KiLCA/imhd/imhd.f90:137,138 and :221,222 — same pi/c pair, with fpc at :139.
4. KiLCA/math/W/wfunc.f90:11, wmfunc.f90:11, wm_asympt.f90:11 — pi = 3.14...197 bare; wfunc.f90:12 then derives srpi = sqrt(pi) from it. All three are in the build (KiLCA/CMakeLists.txt:236-238).
5. KIM/src/util/constants_mod.f90:8 — sol = 29979245800.0. Eleven real consumers: asymptotics/FLR2_asymptotics.f90, background_equilibrium/{calculate_equil,profile_input_m,species_mod}.f90, diagnostics/kim_qldiff_mod.f90, electromagnetic/electromagnetic_solver.f90, electrostatic_poisson/fields_mod.f90, kernels/{FP_kernel_plasma_prefacs,kernel_adaptive,kernel,Krook_kernel_plasma_prefacs}.f90. (KIM/src/math/RKF45.f90 also matches \bsol\b but that is a local intent(out) dummy named sol, not the constant — do not touch it.)
6. KIM/src/util/constants_mod.f90:13 — kB = 1.380649e-16. rg -n '\bkB\b' KIM/ returns this declaration and nothing else — the constant has no consumer. This one is golden-neutral and can be fixed immediately.

Not affected, do not touch:

• KiLCA/core/constants.h:18 and QL-Balance/src/base/constants.h:14 (const double c = 29979245800.0;) — in C/C++ an unsuffixed floating literal is already double. Correct as written.
• python/KIMpy/constants.py:8, python/profile_processor/profile_processor.py:26 — Python floats are doubles.

Dead, out of scope (present in the tree, absent from every CMakeLists.txt): KiLCA/math/W/wfunc_old.f90:14, KiLCA/math/det/test_det.f90:17,18, KiLCA/math/hyper/test_hyper.f90:12, KiLCA/math/W/test_wfunc.f90:7, KiLCA/math/igamma/test_igamma.f90:7, KiLCA/flre/conductivity/tests/Wmn_NaN.f90:35,544, .../tests/test_mods.f90:35, KiLCA/interface/ql-balance_mod.f90:34. (That last one is distinct from the built QL-Balance/src/base/ql-balance_mod.f90 referenced at QL-Balance/CMakeLists.txt:71 — do not confuse them.) Fix opportunistically or leave them; they cannot move any output.

Prior art you must read before touching KiLCA/core/constants_m_*.f90: commit ee1fc521 "fix(kilca): mirror oracle's split speed-of-light precision", which exists only on origin/fortran-port-all, not on main (git branch -a --contains ee1fc521). Its message records that a blanket _dp fix on this exact constant was already tried and reverted, because it made the legacy Fortran path (conduct_parameters → Vth_m, j0th_m) diverge from the golden record while the C++ path (dPhi0, j0th) did not. Its resolution was to make the truncation explicit rather than silent:

real(dp), parameter :: c  = real(29979245800.0_sp, dp)   !speed of light in vacuum; legacy ...

and to give three ported-from-C++ modules a local c = 29979245800.0_dp: kilca_background_data_m (KiLCA/background/background_data_m.f90), kilca_flre_zone_m (KiLCA/flre/flre_zone_m.f90), kilca_transforms_m (KiLCA/mode/transforms_m.f90), each dropping c from its use constants, only: list. That is the documentation fix. This issue is the correction fix, and it supersedes the split.

Reproduction / evidence

Demonstrate the defect standalone (no build needed):

cat > /tmp/prec.f90 <<'EOF'
program prec
    implicit none
    real(8), parameter :: c_leg  = 29979245800.0
    real(8), parameter :: c_true = 29979245800.0_8
    print '(a,es24.17)', 'as written ', c_leg
    print '(a,es24.17)', 'intended   ', c_true
    print '(a,es12.4)',  'rel. error ', (c_leg - c_true)/c_true
end program
EOF
gfortran -o /tmp/prec /tmp/prec.f90 && /tmp/prec

Observed: rel. error 2.6418E-08. Expected after the fix: 0.0000E+00.

Demonstrate the golden impact (this is what gates the PR):

cd <KAMEL>
test/golden/run_golden.sh golden-baseline HEAD

The harness stores no reference data files — it builds golden-baseline and HEAD in the same environment and diffs the two output trees at rtol 1e-7 / atol 1e-12 (test/golden/run_golden.sh:1-27, test/golden/bin/gr_numcompare.py:20-22). The golden-baseline tag is annotated; git rev-parse golden-baseline^{commit} resolves to eb511831 ("Migrate KiLCA 1F1 to fortnum (7/8) (#146)"). Consequently "regenerate the golden record" here means move the golden-baseline tag, not rewrite committed .dat files. Corpus: test/golden/kilca/cases/flre_m6n2/, test/golden/kim/cases/electromagnetic/, test/golden/qlbalance/cases/default/.

Two facts that shape the plan, both established by reading the case fixtures and the zone dispatch:

• test/golden/kilca/cases/flre_m6n2/ contains zone_1_flre.in, zone_2_vacuum.in, zone_3_vacuum.in (rg -n 'plasma model' test/golden/kilca/cases/flre_m6n2/*.in confirms the three model fields are flre, vacuum, vacuum). KiLCA/mode/calc_mode.cpp:146-150 instantiates an hmedium_zone for PLASMA_MODEL_VACUUM — and hmedium_zone.cpp:107,115 plus calc_mode.cpp:314,319,433,517,522 call into the hmedium.f90 subroutines. So hmedium.f90 IS exercised by the golden case and item 2 will move outputs.
• No case in the corpus declares an imhd zone (the only imhd strings under test/golden/ are inside the explanatory comment on the model line), and imhd_zone is only reachable via PLASMA_MODEL_IMHD (calc_mode.cpp:158-162). So item 3 is golden-neutral for the current corpus — but it does change physics for imhd users and must not be smuggled in silently.

Behaviour (BDD)

Scenario: kB is fixed with zero golden impact
  Given KIM/src/util/constants_mod.f90:13 reads "kB       = 1.380649e-16"
  And "rg -n '\bkB\b' KIM/" reports exactly 1 line, the declaration itself
  When kB is changed to 1.380649e-16_dp
  And "test/golden/run_golden.sh golden-baseline HEAD" is run
  Then the kim/electromagnetic case reports 0 files over rtol=1e-7
  And the overall line reads "golden-record overall: PASS"
  And "ctest" passes all 15 registered KIM tests

Scenario: the speed of light is exact in every built Fortran source
  Given the repository at HEAD after the fix
  When "rg -nP '29979245800\.0(?![_d0-9])' --glob '!build/**' -g '*.f90'" is run
  Then the only remaining hits are the 8 files listed under "Dead, out of scope"
  And no hit is in a file listed in KiLCA/CMakeLists.txt or under KIM/src/
  And KiLCA/core/constants_m_64bit.f90:35 declares "c  = 2.99792458e10_dp"
  And KiLCA/core/constants_m_32bit.f90:35 declares the same literal
  And no built source declares a local c, sol, or pi shadowing the shared module

Scenario: the golden shift is bounded and attributed
  Given the constants fix applied on a branch
  When "test/golden/run_golden.sh golden-baseline HEAD" is run
  Then every reported per-file max relative deviation is below 1e-5
  And the PR body tabulates, per output file, the observed max relative deviation
  And no file is newly NaN, newly infinite, or newly absent versus the baseline

Scenario: the new baseline is self-consistent
  Given the golden-baseline tag has been moved to the merge commit of the constants fix
  When "test/golden/run_golden.sh golden-baseline HEAD" is run on the next commit to main
  Then all three codes report PASS with 0 files over rtol=1e-7

Implementation plan

Split into two PRs. PR-A is golden-neutral and mergeable today; PR-B is the shifting change and needs the re-baseline.

PR-A — golden-neutral hygiene

1. KIM/src/util/constants_mod.f90:13 — change kB = 1.380649e-16 to kB = 1.380649e-16_dp. (Deleting the line is also defensible since it has no consumer; prefer fixing, so the module stays a complete CGS reference.)
   Do not touch :7 (pi = 3.14159265358979d0) in this PR. That literal is correctly parsed as double but truncated at 15 digits; measured relative error vs 4*atan(1.0_dp) is -9.9e-16 (~4 ulp), which is below rtol 1e-7 per-value but can still flip last bits through the kernels. Defer it to PR-B where a golden shift is already expected.
2. KIM/tests/CMakeLists.txt — add a new test following the existing pattern. That file currently holds 15 add_executable/add_test pairs, each of the shape:

   add_executable(test_constants ${CMAKE_SOURCE_DIR}/KIM/tests/test_constants.f90)
   set_target_properties(test_constants PROPERTIES
                           OUTPUT_NAME test_constants.x
                           RUNTIME_OUTPUT_DIRECTORY "${CMAKE_BINARY_DIR}/tests/")
   target_link_libraries(test_constants KIM_lib)
   add_test(NAME test_constants
            COMMAND ${CMAKE_BINARY_DIR}/tests/test_constants.x)

   Write KIM/tests/test_constants.f90 to error stop unless sol == 2.99792458e10_dp exactly, kB == 1.380649e-16_dp exactly, and abs(pi - 4.0_dp*atan(1.0_dp)) <= 2.0_dp*epsilon(1.0_dp)*pi. This test is what stops the regression coming back.
3. KiLCA/CMakeLists.txt — no change. The dead files in §"Dead, out of scope" may be corrected in the same PR for tidiness; they cannot affect any test.

PR-B — the correction (shifts the golden record)

4. KiLCA/core/constants_m_64bit.f90:35 and KiLCA/core/constants_m_32bit.f90:35 — replace the bare literal with real(dp), parameter :: c = 2.99792458e10_dp. The two files are currently identical except line 7 (pp = 8 vs pp = 4); keep them that way.
5. KiLCA/background/background_data_m.f90, KiLCA/flre/flre_zone_m.f90, KiLCA/mode/transforms_m.f90 — only if the branch you work from already carries ee1fc521. On main at 2722090f those local c declarations do not exist, so this step is a no-op there. If ee1fc521/PR Port remaining KiLCA C++ to Fortran (fortnum-backed, single PR) #161 has landed, delete the three local real(dp), parameter :: c = 29979245800.0_dp declarations and restore c to each module's use constants, only: list. With the module constant now exact, the split ee1fc521 deliberately created is no longer needed and keeping it would leave two names for one quantity. Ordering note: PR-B must rebase onto whatever state of Port remaining KiLCA C++ to Fortran (fortnum-backed, single PR) #161 is on main.
6. KiLCA/hom_medium/hmedium.f90 lines 21, 23, 106, 108, 351, 352, 408, 409 — delete the eight local parameter declarations and add use constants, only: pi, c at the top of each of the four affected subroutines (these are free-standing subroutines, not module procedures, so the use goes immediately after the subroutine line and before implicit none). Keep the derived fpc = 4.0d0*pi/c where it is; it will pick up the exact values automatically. Consolidating rather than re-suffixing in place is the point: eight independent copies of c is how this defect propagated.
7. KiLCA/imhd/imhd.f90 lines 137, 138, 221, 222 — same consolidation onto use constants, only: pi, c.
8. KiLCA/math/W/wfunc.f90:11, wmfunc.f90:11, wm_asympt.f90:11 — replace the local pi with use constants, only: pi. In wfunc.f90 the module already exports srpi = sqrt(pi) (constants_m_64bit.f90:27), so use use constants, only: pi, srpi and delete the local srpi = sqrt(pi) at :12.
9. KIM/src/util/constants_mod.f90:8 — sol = 2.99792458e10_dp. Optionally also :7 — pi = 3.141592653589793238462643383279502884197_dp (deferred here from PR-A).
10. Run test/golden/run_golden.sh golden-baseline HEAD locally. Capture the full per-file max-relative-deviation table from gr_numcompare.py and paste it into the PR body. Any file deviating by more than 1e-5 is a red flag that something other than the constant changed — investigate before proceeding.
11. After PR-B merges: retag the baseline on the merge commit (git tag -f golden-baseline <sha> && git push -f origin golden-baseline). Also check the pre-fortnum diagnostic job in .github/workflows/golden.yml, which runs test/golden/run_golden.sh "$PRE_FORTNUM_BASELINE" ${{ github.sha }} at :82 against PRE_FORTNUM_BASELINE: 0fb4f296... (:10) — it will legitimately show the constants shift. gr_compare.sh honors a GOLDEN_RECORD_SKIP_CASES env var (gr_compare.sh:3), but note that it skips whole cases, not individual files, so either skip the affected case for that job or retire the diagnostic job if the fortnum migration is considered closed.
12. docs/plans/ — extend docs/plans/2026-06-30-golden-record-ci.md (it exists) with a short note recording which baseline shift corresponds to which physics change, so a future reader can tell the constants re-baseline apart from a genuine regression.

No sibling-repo change is required. libneo declares C = 2.99792458e10_dp correctly in src/util/math_constants.f90:26, but KAMEL links only LIBNEO::collision_freqs, polylag, species, transport (QL-Balance/CMakeLists.txt:146-149) and LIBNEO::magfie (common/equil/CMakeLists.txt:22) — pulling libneo into KiLCA solely for constants is a heavier coupling than this fix warrants. Adopting math_constants as the single org-wide source is a reasonable follow-up issue, not part of this one.

Acceptance criteria

• rg -nP '29979245800\.0(?![_d0-9])' --glob '!build/**' -g '*.f90' -g '*.F90' returns no hit in any source listed in a CMakeLists.txt (the -P flag is required — rg's default engine rejects look-ahead).
• rg -n 'parameter\s*::.*=\s*[0-9]+\.[0-9]{8,}\s*[;!]?\s*$' --glob '!build/**' -g '*.f90' | grep -v '_dp\|d[+-]\?[0-9]' returns no hit in a built source. (At 2722090f this returns 16 hits; after PR-B only the dead-file hits should remain.)
• make all && ctest passes; test_constants is registered and green, bringing the KIM test count from 15 to 16.
• PR-A: test/golden/run_golden.sh golden-baseline HEAD prints golden-record overall: PASS with zero files over rtol=1e-7.
• PR-B: the golden run's per-file max relative deviations are all < 1e-5, tabulated in the PR body, with no new NaN/Inf and no missing output file.
• After the re-tag, test/golden/run_golden.sh golden-baseline HEAD on the following main commit is clean.
• KiLCA/core/constants_m_32bit.f90 and constants_m_64bit.f90 still differ only at line 7 (pp); verify with diff returning exactly that one hunk.

References

The exact value c = 2.99792458e10 cm/s is definitional: the metre was redefined in terms of a fixed speed of light at the 17th CGPM (1983, Resolution 1). k_B = 1.380649e-16 erg/K is definitional under the revised SI adopted at the 26th CGPM (2018, Resolution 1), effective 20 May 2019. Both are exact, so there is no accuracy/precision trade-off to argue about.

The _dp-suffix discipline for literals is standard Fortran practice; Metcalf, Reid & Cohen, Modern Fortran Explained (Oxford University Press) covers kind-parameterised literals and the default-real trap.

Otherwise: no specific literature reference applies; this is an implementation detail internal to the code. The in-repo model to imitate is libneo's src/util/math_constants.f90, and the golden/regression methodology is defined by docs/plans/2026-06-30-golden-record-ci.md, which is the authority here.

Note for the implementer: QL-Balance/src/base/ql-balance.f90:4 attributes the balance equations to "Heyn et al. NF2014", but the header gives no volume or article number, and this triage pass could not verify the full bibliographic record. Do not cite a volume/page for it in the PR body without checking the actual reference first.

Related issues

• ql-balance golden record: near-resonant root location is ULP-sensitive (pre-existing, not a fortnum regression) #164 and ql-balance/KiLCA golden record: run-to-run nondeterminism in the linear solve (root cause + fixes) #172 — ULP-sensitivity and nondeterminism in the golden records. Both are referenced directly in the harness source (test/golden/bin/gr_numcompare.py:29, test/golden/bin/test_gr_numcompare.py:78,105), which documents that zone_*_poy_test_err.dat amplifies last-bit differences ~1e8x. Both will be easier to reason about once a 2.6e-8 systematic bias is off the table; do not merge PR-B concurrently with work on those, or the attributions will be entangled.
• Port remaining KiLCA C++ to Fortran (fortnum-backed, single PR) #161 — the fortran-port-all PR the issue body itself names, carrying ee1fc521. See step 5 for the ordering constraint.
• KiLCA flre_m6n2 golden: FLRE 1F1m is ULP-sensitive near resonance (pre-existing, not a port regression) #165 and [math-06] Verify equilibrium and plasma-background formulas #198 were referenced in an earlier draft of this triage. Neither is referenced anywhere in the checkout, so their content could not be verified here — check them on the tracker before relying on any claimed overlap.
• Worth flagging separately (not part of this issue): KiLCA/core/constants_m_64bit.f90:25 declares boltz = 1.60216428d-12 ! erg/ev, correctly double-parsed but numerically off from the exact 1.602176634e-12 by 8e-6 relative — three orders of magnitude worse than the precision defect this issue is about. That is a wrong-value problem, not a wrong-kind problem, and deserves its own issue with its own golden re-baseline. The same value appears in KiLCA/core/constants.h:16 and QL-Balance/src/base/constants.h:12.

Scope and risk

Two PRs as scoped above: PR-A is a handful of lines plus one new test and carries no risk. PR-B touches six built source files and will move golden outputs by ~1e-8 relative, so the real risk is procedural — merging it without moving the golden-baseline tag turns CI red for every subsequent PR, and merging it alongside #164/#172 work makes all three undebuggable. Land PR-B alone, on a quiet main, with the deviation table in the PR body.