Skip to content

singular/ca_left and ca_right are written as uninitialized memory on kinetic runs #366

Description

@logan-nc

Summary

singular/ca_left and singular/ca_right are written to gpec.h5 unconditionally, but
nothing populates them on kinetic runs. The datasets therefore contain uninitialized
memory
whenever kinetic_factor > 0.

Filing as a sub-issue of #226 because that effort already rewrites every singular/ write
site (singular/SingularSurfaces/), so the guard belongs in the same pass rather than
as a separate drive-by change.

Mechanism

  • The arrays are allocated with undef and never zeroed:
    src/ForceFreeStates/ForceFreeStatesStructs.jl:615-616
    ca_r::Array{ComplexF64,4} = Array{ComplexF64}(undef, numpert_total, numpert_total, 2, msing)
    ca_l::Array{ComplexF64,4} = Array{ComplexF64}(undef, numpert_total, numpert_total, 2, msing)
  • They are only ever filled by the ideal surface crossings —
    cross_ideal_singular_surf! (EulerLagrange.jl:594, 631) and the Riccati ideal crossing
    (Riccati.jl:1226, 1267).
  • cross_kinetic_singular_surf! (EulerLagrange.jl) does no asymptotic analysis and never
    writes them, so on a kinetic run they keep whatever was on the heap.
  • They are written out regardless: src/GeneralizedPerturbedEquilibrium.jl:839-840
    out_h5["singular/ca_left"]  = odet.ca_l
    out_h5["singular/ca_right"] = odet.ca_r

Evidence

Running examples/Solovev_kinetic_calculated_example twice with identical code and
diffing singular/ca_left (4096 elements) gives a different number of NaNs each time —
238 in one run, 100 in the other. One ca_right element came back as 7.026e-310, a
subnormal double, which is the classic signature of reading uninitialized memory rather
than a computed value.

This surfaced while validating #365 (on-demand solution derivatives): a develop-vs-branch
comparison flagged these two datasets as differing even though the change is provably
numerics-neutral everywhere else. It is pre-existing on develop and unrelated to that PR.

Impact

Low severity, but genuinely bad behaviour:

  • Output is non-reproducible run to run for a deterministic calculation, which
    undermines any bitwise output comparison (regression harness, benchmark diffing, rerun
    from gpec.h5).
  • Consumers cannot distinguish "not computed" from "computed and happens to be small".
  • NaNs propagate silently if anything downstream ever starts reading these.

No current consumer is affected: the only reader is compute_delta_prime_from_ca!
(EulerLagrange.jl:16), documented as a stub that is not physically valid and not
called from any integration driver.

Suggested fix

Pick whichever fits the new schema best:

  1. Allocate with zeros instead of undef — one-line, and makes "not computed" read as
    zero rather than garbage. Cheapest, and safe since these arrays are small (msing
    surfaces).
  2. Skip the datasets entirely when they were never populated (kinetic runs), matching how
    ForceFreeStates - IMPROVEMENT - Compute solution derivatives on demand #365 now handles the derivative stores — empty dataset rather than meaningless values.
  3. If they are to stay always-present, add an explicit computed flag / attribute so
    readers can tell.

Option 1 plus a note in the dataset docs is probably enough; option 2 is more consistent
with how the rest of the schema is heading.

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Type

    No type

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions