From 260a38a47409a6554e4fdc54c022225f705bb9ea Mon Sep 17 00:00:00 2001 From: logan-nc Date: Thu, 13 Aug 2026 09:56:29 -0400 Subject: [PATCH 1/2] FFS - BUGFIX - Emit zero-extent ca_left/ca_right when unpopulated MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit SingularSurfaces/ca_left and ca_right were written to gpec.h5 unconditionally, but only ideal singular-surface crossings populate odet.ca_l/ca_r — kinetic and galerkin-matched runs dumped uninitialized heap memory (non-reproducible NaN counts between identical runs; subnormal garbage values). - Allocate ca_l/ca_r with zeros instead of undef (kills in-memory nondeterminism at the root). - New OdeState.ca_populated flag (mirrors du_store_populated), set by the two ideal crossing routines (EulerLagrange + Riccati) and carried through the Riccati dense-xi save/restore; galerkin-matched OdeStates keep the default false. - The writer emits rank-4 zero-extent sentinels when the flag is false — datasets stay always-present (no reader KeyErrors, metadata annotations unchanged), matching the established empty-sentinel idiom. - Guard the two readers that index ca: Analysis.plot_delta_prime returns a placeholder on empty ca, and the gal_{epsilon,beta}_scan benchmarks record NaN for the ca-jump diagnostic (they were consuming garbage on galerkin runs already); also fix a leftover legacy haskey(f["singular"], ...) in both. - Zero-extent not-computed sentinel codified in hdf5-conventions.md; stale stability.md claim about SingularCoupling reading ca_l/ca_r corrected. - Tests: kinetic fullrun asserts empty ca datasets; ideal schema run asserts populated + finite. New solovev_n1 harness quantity checksums ca_left to lock bitwise reproducibility. Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_0129rSTCmYJDBbcH9khHqYnz --- benchmarks/gal_validation/gal_beta_scan.jl | 4 ++-- benchmarks/gal_validation/gal_epsilon_scan.jl | 4 ++-- docs/development/hdf5-conventions.md | 1 + docs/src/stability.md | 6 ++++-- regression-harness/cases/solovev_n1.toml | 10 ++++++++++ src/Analysis/ForceFreeStates.jl | 2 ++ src/ForceFreeStates/EulerLagrange.jl | 1 + src/ForceFreeStates/ForceFreeStatesStructs.jl | 9 +++++++-- src/ForceFreeStates/Riccati.jl | 3 +++ src/GeneralizedPerturbedEquilibrium.jl | 6 ++++-- src/HDF5Schema.jl | 6 ++++-- test/runtests_fullruns.jl | 4 ++++ test/runtests_h5_schema.jl | 7 +++++++ 13 files changed, 51 insertions(+), 12 deletions(-) diff --git a/benchmarks/gal_validation/gal_beta_scan.jl b/benchmarks/gal_validation/gal_beta_scan.jl index 11b8fd69e..b9f8f0e4d 100644 --- a/benchmarks/gal_validation/gal_beta_scan.jl +++ b/benchmarks/gal_validation/gal_beta_scan.jl @@ -26,7 +26,7 @@ function extract(h5) # surface m-values (STRIDE side) m_s = vec(read(f, "SingularSurfaces/m")); n_s = vec(read(f, "SingularSurfaces/n")); msing = length(m_s) # STRIDE delta_prime_matrix - dpm = haskey(f["singular"], "delta_prime_matrix") ? read(f, "SingularSurfaces/delta_prime_matrix") : nothing + dpm = haskey(f, "SingularSurfaces/delta_prime_matrix") ? read(f, "SingularSurfaces/delta_prime_matrix") : nothing cal = read(f, "SingularSurfaces/ca_left"); car = read(f, "SingularSurfaces/ca_right") # gal gm = haskey(f, "SingularSurfaces/GalerkinDeltaPrime/sing_m") ? vec(read(f, "SingularSurfaces/GalerkinDeltaPrime/sing_m")) : Int[] @@ -35,7 +35,7 @@ function extract(h5) key = "m$(m)" # ca-jump ipr = 1 + m - mlow + (n_s[s] - nlow) * mpert - out["cajump_$key"] = (car[ipr, ipr, 2, s] - cal[ipr, ipr, 2, s]) / denom + out["cajump_$key"] = isempty(cal) ? NaN + NaN * im : (car[ipr, ipr, 2, s] - cal[ipr, ipr, 2, s]) / denom # zero-extent when not computed # stride out["stride_$key"] = dpm !== nothing && s <= size(dpm, 1) ? dpm[s, s] : NaN + NaN*im # gal (match by m) diff --git a/benchmarks/gal_validation/gal_epsilon_scan.jl b/benchmarks/gal_validation/gal_epsilon_scan.jl index 04fe00f13..7ff7653fe 100644 --- a/benchmarks/gal_validation/gal_epsilon_scan.jl +++ b/benchmarks/gal_validation/gal_epsilon_scan.jl @@ -26,7 +26,7 @@ function extract(h5) # surface m-values (STRIDE side) m_s = vec(read(f, "SingularSurfaces/m")); n_s = vec(read(f, "SingularSurfaces/n")); msing = length(m_s) # STRIDE delta_prime_matrix - dpm = haskey(f["singular"], "delta_prime_matrix") ? read(f, "SingularSurfaces/delta_prime_matrix") : nothing + dpm = haskey(f, "SingularSurfaces/delta_prime_matrix") ? read(f, "SingularSurfaces/delta_prime_matrix") : nothing cal = read(f, "SingularSurfaces/ca_left"); car = read(f, "SingularSurfaces/ca_right") # gal gm = haskey(f, "SingularSurfaces/GalerkinDeltaPrime/sing_m") ? vec(read(f, "SingularSurfaces/GalerkinDeltaPrime/sing_m")) : Int[] @@ -35,7 +35,7 @@ function extract(h5) key = "m$(m)" # ca-jump ipr = 1 + m - mlow + (n_s[s] - nlow) * mpert - out["cajump_$key"] = (car[ipr, ipr, 2, s] - cal[ipr, ipr, 2, s]) / denom + out["cajump_$key"] = isempty(cal) ? NaN + NaN * im : (car[ipr, ipr, 2, s] - cal[ipr, ipr, 2, s]) / denom # zero-extent when not computed # stride out["stride_$key"] = dpm !== nothing && s <= size(dpm, 1) ? dpm[s, s] : NaN + NaN*im # gal (match by m) diff --git a/docs/development/hdf5-conventions.md b/docs/development/hdf5-conventions.md index b24610d2a..fd3ff96ba 100644 --- a/docs/development/hdf5-conventions.md +++ b/docs/development/hdf5-conventions.md @@ -57,6 +57,7 @@ Mechanism: writers stay table-driven — each writer keeps a `path => (; long_na - Complex numbers are stored as the native HDF5.jl compound type (readable by h5py as a compound dtype). - `NaN` is the not-computed sentinel in numeric datasets (e.g. auto-derived settings, rootless growth-rate entries). +- A **zero-extent array** is the not-computed sentinel for whole datasets that a given run never produces (e.g. `SingularSurfaces/ca_left`/`ca_right` on kinetic or galerkin-matched runs, the free-boundary energies when `vac_flag=false`, the on-demand derivative stores). Never write unpopulated (`undef`) memory. - Ragged (variable-length) data uses the flat-plus-`offsets` companion pattern (`offsets[k+1] - offsets[k]` = length of row `k`) rather than HDF5 VLEN types, e.g. `KineticForces//EnergyIntegrals/` and `Tearing/Diagnostics/*`. ## Back-compatibility policy diff --git a/docs/src/stability.md b/docs/src/stability.md index ccd061460..b6f48e6ff 100644 --- a/docs/src/stability.md +++ b/docs/src/stability.md @@ -342,8 +342,10 @@ end ## Notes -- The standard path does not populate `delta_prime`; use `PerturbedEquilibrium.SingularCoupling` - for Δ' on the standard path (it reads `ca_l`/`ca_r` directly). +- The standard path does not populate `delta_prime`; the canonical Δ' is the STRIDE BVP + `SingularSurfaces/delta_prime_matrix` from the parallel FM path. `ca_l`/`ca_r` are filled + only by ideal surface crossings (kinetic and galerkin-matched runs emit zero-extent + `ca_left`/`ca_right` sentinels). - The Riccati and parallel FM paths compute Δ' inline at each crossing, using the direct diagonal formula (no GR permutation). The result in `delta_prime_col[ipert_res, i]` equals `delta_prime[i]` to machine precision. diff --git a/regression-harness/cases/solovev_n1.toml b/regression-harness/cases/solovev_n1.toml index d6fab13e3..0930b0dad 100644 --- a/regression-harness/cases/solovev_n1.toml +++ b/regression-harness/cases/solovev_n1.toml @@ -182,6 +182,16 @@ label = "pressure profile (checksum)" noise_threshold = 0 order = 71 +# Bitwise reproducibility guard for the ideal asymptotic coefficients (issue: kinetic +# runs used to dump uninitialized memory here; ideal values must stay byte-stable). +[quantities.ca_left] +h5path = "SingularSurfaces/ca_left" +type = "complex_matrix" +extract = "checksum" +label = "ca_left (checksum)" +noise_threshold = 0 +order = 72 + # Runtime (special: not from H5) [quantities.runtime] h5path = "" diff --git a/src/Analysis/ForceFreeStates.jl b/src/Analysis/ForceFreeStates.jl index c92d6e698..14062f68e 100644 --- a/src/Analysis/ForceFreeStates.jl +++ b/src/Analysis/ForceFreeStates.jl @@ -327,6 +327,8 @@ function plot_delta_prime(h5path; save_path=nothing) end msing == 0 && return plot(; title="No singular surfaces found", legend=false) + # ca_left/ca_right are zero-extent sentinels on kinetic/galerkin-matched runs (never computed there). + isempty(ca_l) && return plot(; title="No asymptotic coefficients — ca_left/ca_right not computed for this run", legend=false) numpert_total = size(ca_l, 1) chi1 = 2π * psio diff --git a/src/ForceFreeStates/EulerLagrange.jl b/src/ForceFreeStates/EulerLagrange.jl index af73d37da..4e7558807 100644 --- a/src/ForceFreeStates/EulerLagrange.jl +++ b/src/ForceFreeStates/EulerLagrange.jl @@ -630,6 +630,7 @@ function cross_ideal_singular_surf!( end # Get asymptotic coefficients after crossing rational surface odet.ca_r[:, :, :, ising] .= sing_get_ca(odet.u, ua, intr) + odet.ca_populated = true # Δ' is NOT computed for the standard path. The physical Δ' requires the solution # columns to be in the Riccati gauge (U₂=I), maintained only by Riccati renormalization. diff --git a/src/ForceFreeStates/ForceFreeStatesStructs.jl b/src/ForceFreeStates/ForceFreeStatesStructs.jl index b32b59df7..c9a03dd8b 100644 --- a/src/ForceFreeStates/ForceFreeStatesStructs.jl +++ b/src/ForceFreeStates/ForceFreeStatesStructs.jl @@ -541,6 +541,10 @@ and a small set of temporary matrices and factors used to compute singular-layer - `ca_l::Array{ComplexF64,4}` - Asymptotic coefficients just to the left of each singular surface with shape `(numpert_total, numpert_total, 2, msing)`. + - `ca_populated::Bool` - True once an ideal singular-surface crossing has filled `ca_l`/`ca_r`; kinetic and + galerkin-matched runs never populate them and leave this false, and the HDF5 writer then emits zero-extent + `ca_left`/`ca_right` datasets instead of unpopulated arrays. + - `edge_scan::EdgeScanState` - Edge dW scan state and results. Initialized as a disabled sentinel (N_edge=0) and replaced by `findmax_dW_edge!` when a scan runs. - `psifac::Float64` - Current normalized flux coordinate for the integrator. @@ -604,8 +608,9 @@ and a small set of temporary matrices and factors used to compute singular-layer xi_s_store::Array{ComplexF64,3} = Array{ComplexF64}(undef, numpert_total, numpert_total, numsteps_init) du_store_populated::Bool = false crit_store::Vector{Float64} = Vector{Float64}(undef, numsteps_init) - 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) + ca_r::Array{ComplexF64,4} = zeros(ComplexF64, numpert_total, numpert_total, 2, msing) + ca_l::Array{ComplexF64,4} = zeros(ComplexF64, numpert_total, numpert_total, 2, msing) + ca_populated::Bool = false # Edge dW scan state and results (disabled sentinel when psiedge >= psilim, i.e. no edge scan) edge_scan::EdgeScanState = EdgeScanState(numpert_total, 0) diff --git a/src/ForceFreeStates/Riccati.jl b/src/ForceFreeStates/Riccati.jl index dbdd4440b..c57c8648a 100644 --- a/src/ForceFreeStates/Riccati.jl +++ b/src/ForceFreeStates/Riccati.jl @@ -1271,6 +1271,7 @@ function _capture_right_crossing_data!(odet::OdeState, singp::SingType, sing_asy end end odet.ca_r[:, :, :, ising] .= sing_get_ca(odet.u, ua, intr) + odet.ca_populated = true end # STUB: per-surface ca-based Δ' (not physically valid; see SingType.delta_prime docstring). @@ -1947,6 +1948,7 @@ function _populate_dense_xi_via_serial_el!( qlim = intr.qlim, ca_l = copy(odet.ca_l), ca_r = copy(odet.ca_r), + ca_populated = odet.ca_populated, sing_state = [( delta_prime = copy(intr.sing[s].delta_prime), delta_prime_col = copy(intr.sing[s].delta_prime_col), @@ -1992,6 +1994,7 @@ function _populate_dense_xi_via_serial_el!( # written against the (S, I) Riccati convention. fresh_odet.ca_l .= saved.ca_l fresh_odet.ca_r .= saved.ca_r + fresh_odet.ca_populated = saved.ca_populated # Return the fresh serial-EL odet (self-consistent for ξ-function # storage in axis basis; `ca_l`/`ca_r` carry the parallel-BVP diff --git a/src/GeneralizedPerturbedEquilibrium.jl b/src/GeneralizedPerturbedEquilibrium.jl index 656316b75..573eaf204 100755 --- a/src/GeneralizedPerturbedEquilibrium.jl +++ b/src/GeneralizedPerturbedEquilibrium.jl @@ -845,8 +845,10 @@ function write_outputs_to_HDF5( out_h5["SingularSurfaces/psi"] = [sing.psifac for sing in intr.sing] out_h5["SingularSurfaces/q"] = [sing.q for sing in intr.sing] out_h5["SingularSurfaces/q1"] = [sing.q1 for sing in intr.sing] - out_h5["SingularSurfaces/ca_left"] = odet.ca_l - out_h5["SingularSurfaces/ca_right"] = odet.ca_r + # Kinetic and galerkin-matched runs never populate ca_l/ca_r (only ideal surface + # crossings do); emit zero-extent sentinels instead of unpopulated arrays. + out_h5["SingularSurfaces/ca_left"] = odet.ca_populated ? odet.ca_l : zeros(ComplexF64, 0, 0, 0, 0) + out_h5["SingularSurfaces/ca_right"] = odet.ca_populated ? odet.ca_r : zeros(ComplexF64, 0, 0, 0, 0) if intr.msing > 0 # Mode numbers at each surface (jagged — pad with 0 to max_modes width) diff --git a/src/HDF5Schema.jl b/src/HDF5Schema.jl index ace08b783..49fd3989c 100644 --- a/src/HDF5Schema.jl +++ b/src/HDF5Schema.jl @@ -121,9 +121,11 @@ const MAIN_H5_ANNOTATIONS = [ "SingularSurfaces/n" => (; long_name="resonant toroidal mode numbers per surface (0-padded)", dims=("surface", "mode")), "SingularSurfaces/di0" => (; long_name="Mercier D_I evaluated at each rational surface", dims=("surface",)), "SingularSurfaces/ca_left" => - (; long_name="asymptotic large/small-solution coefficient matrices just left of each surface", dims=("mode", "solution", "large_small", "surface")), + (; long_name="asymptotic large/small-solution coefficient matrices just left of each surface (zero-extent when not computed — ideal crossings only)", + dims=("mode", "solution", "large_small", "surface")), "SingularSurfaces/ca_right" => - (; long_name="asymptotic large/small-solution coefficient matrices just right of each surface", dims=("mode", "solution", "large_small", "surface")), + (; long_name="asymptotic large/small-solution coefficient matrices just right of each surface (zero-extent when not computed — ideal crossings only)", + dims=("mode", "solution", "large_small", "surface")), "SingularSurfaces/E" => (; long_name="Glasser-Greene-Johnson coefficient E per surface", dims=("surface",)), "SingularSurfaces/F" => (; long_name="Glasser-Greene-Johnson coefficient F per surface", dims=("surface",)), "SingularSurfaces/G" => (; long_name="Glasser-Greene-Johnson coefficient G per surface", dims=("surface",)), diff --git a/test/runtests_fullruns.jl b/test/runtests_fullruns.jl index ce6470385..a76515bf5 100644 --- a/test/runtests_fullruns.jl +++ b/test/runtests_fullruns.jl @@ -28,6 +28,10 @@ using HDF5 et = read(h5["ForceFreeStates/FreeBoundaryStability/eigenmode_energies"]) @test isfinite(real(et[1])) @test real(et[1]) > 0 + # Kinetic runs never populate the asymptotic ca coefficients; the writer must + # emit deterministic zero-extent sentinels, not uninitialized memory. + @test isempty(read(h5["SingularSurfaces/ca_left"])) + @test isempty(read(h5["SingularSurfaces/ca_right"])) end rm(joinpath(ex3, "gpec.h5"); force=true) true diff --git a/test/runtests_h5_schema.jl b/test/runtests_h5_schema.jl index 3e11f978a..b2ec0816e 100644 --- a/test/runtests_h5_schema.jl +++ b/test/runtests_h5_schema.jl @@ -104,6 +104,13 @@ end fwd = "ForceFreeStates/Solutions/ForwardIntegration" @test HDF5.API.h5ds_is_scale(h5["$fwd/psi"]) @test HDF5.API.h5ds_is_attached(h5["$fwd/q"], h5["$fwd/psi"], 0) + + # Ideal run with rational surfaces: the asymptotic ca coefficients are + # populated and finite (the kinetic/galerkin not-computed case emits + # zero-extent sentinels instead — asserted in runtests_fullruns.jl). + ca_l = read(h5["SingularSurfaces/ca_left"]) + @test !isempty(ca_l) + @test all(isfinite, ca_l) end end end From 809366db34c50e0fe158f4bdfb1ba6daea43f4f1 Mon Sep 17 00:00:00 2001 From: logan-nc Date: Sat, 15 Aug 2026 09:57:22 -0400 Subject: [PATCH 2/2] DOCS - FIX - Capitalize the Delta_prime_matrix HDF5 path in stability.md The note added by this PR cited the dataset as SingularSurfaces/delta_prime_matrix; the literature-capitalization rule introduced with the schema overhaul writes it as Delta_prime_matrix. Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_01VGmFuAw5JdYrAyBSXCVssR --- docs/src/stability.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/stability.md b/docs/src/stability.md index 82f636364..0139384f6 100644 --- a/docs/src/stability.md +++ b/docs/src/stability.md @@ -343,7 +343,7 @@ end ## Notes - The standard path does not populate `delta_prime`; the canonical Δ' is the STRIDE BVP - `SingularSurfaces/delta_prime_matrix` from the parallel FM path. `ca_l`/`ca_r` are filled + `SingularSurfaces/Delta_prime_matrix` from the parallel FM path. `ca_l`/`ca_r` are filled only by ideal surface crossings (kinetic and galerkin-matched runs emit zero-extent `ca_left`/`ca_right` sentinels). - The Riccati and parallel FM paths compute Δ' inline at each crossing, using the