FFS - BUGFIX - Zero-extent ca_left/ca_right when unpopulated (issue #366)

docs/development/hdf5-conventions.md

@@ -67,6 +67,7 @@ Mechanism: writers stay table-driven — each writer keeps a `path => (; long_na
 
 - Complex quantities are stored as the native HDF5.jl compound type (readable by h5py as a compound dtype) — **never split into `*_real`/`*_imag` dataset pairs**. Sole sanctioned exception: `Input/RawInputs/ForcingTerms/amplitude_{real,imag}`, which mirrors the external forcing ingest-file format and keeps pre-existing snapshots replayable.
 - `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/<method>/EnergyIntegrals/` and `Tearing/Diagnostics/*`.
 
 ## Back-compatibility policy

docs/development/regression-harness.md

@@ -98,14 +98,12 @@ regress --cases solovev_n1 --ref-range develop~10..develop
 - `--force` — re-run even if cached
 - `--verbose` — print GPEC subprocess output
 - `--no-instantiate` — skip `Pkg.instantiate()` (faster if deps are already resolved)
-<<<<<<< HEAD
-
-GPEC subprocesses run with `-t auto` (all cores) so GPEC's threaded kernels are active; set `GPEC_REGRESS_THREADS=1` to force single-threaded runs. Tracked quantities are thread-count independent, and the count each run actually used is recorded in its environment fingerprint (shown in the report's `env:` lines). Thread count is deliberately not part of the cache key, so `Runtime (s)` rows cached from single-threaded runs are not comparable to threaded ones — re-baseline with `--force` if runtime tracking matters.
-=======
 - `--no-pin-manifest` — let each ref resolve its own package set (see below)
 - `--allow-env-mismatch` — reuse cached results produced in a different environment
 - `--fail-on-change` — exit non-zero when any tracked quantity changed
 
+GPEC subprocesses run with `-t auto` (all cores) so GPEC's threaded kernels are active; set `GPEC_REGRESS_THREADS=1` to force single-threaded runs. Tracked quantities are thread-count independent, and the count each run actually used is recorded in its environment fingerprint (shown in the report's `env:` lines). Thread count is deliberately not part of the cache key, so `Runtime (s)` rows cached from single-threaded runs are not comparable to threaded ones — re-baseline with `--force` if runtime tracking matters.
+
 ## Making source code the only variable
 
 `Manifest.toml` is untracked, so a worktree checked out at an old commit used to resolve whatever
@@ -140,11 +138,9 @@ than leaving you to infer it from the numbers.
 Results cached before environment fingerprinting existed carry no environment and are therefore
 re-run once — those are exactly the entries whose provenance cannot be established.
 
-Thread counts are recorded but **not** forced: the harness does not silently change how your runs
-execute. If the two refs in a comparison ran under different thread counts, the report flags it.
+If the two refs in a comparison ran under different thread counts, the report flags it.
 
 ## Exit status
 
 - `0` — every run completed (and, with `--fail-on-change`, nothing changed)
 - `1` — a run failed, or a quantity changed under `--fail-on-change`
->>>>>>> performance/regression-harness-worktree-reuse

docs/src/stability.md

@@ -363,8 +363,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.

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 = ""

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

src/ForceFreeStates/EulerLagrange.jl

@@ -641,6 +641,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.

src/ForceFreeStates/ForceFreeStatesStructs.jl

@@ -532,6 +532,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.
@@ -591,8 +595,9 @@ and a small set of temporary matrices and factors used to compute singular-layer
     u_store_el_basis::Bool = true
     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)

src/ForceFreeStates/Riccati.jl

@@ -1265,6 +1265,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).

src/GeneralizedPerturbedEquilibrium.jl

@@ -1205,8 +1205,10 @@ function write_outputs_to_HDF5(
         out_h5["SingularSurfaces/rational_psi"] = [sing.psifac for sing in result.surfaces]
         out_h5["SingularSurfaces/rational_q"] = [sing.q for sing in result.surfaces]
         out_h5["SingularSurfaces/dqdpsi"] = [sing.q1 for sing in result.surfaces]
-        out_h5["SingularSurfaces/ca_left"] = diag !== nothing ? diag.ca_l : ComplexF64[]
-        out_h5["SingularSurfaces/ca_right"] = diag !== nothing ? diag.ca_r : ComplexF64[]
+        # 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"] = (diag !== nothing && diag.ca_populated) ? diag.ca_l : zeros(ComplexF64, 0, 0, 0, 0)
+        out_h5["SingularSurfaces/ca_right"] = (diag !== nothing && diag.ca_populated) ? diag.ca_r : zeros(ComplexF64, 0, 0, 0, 0)
 
         if msing > 0
             # Mode numbers at each surface (jagged — pad with 0 to max_modes width)

src/HDF5Schema.jl

@@ -181,9 +181,11 @@ const MAIN_H5_ANNOTATIONS = [
     "SingularSurfaces/D_I" =>
         (; long_name="Mercier D_I evaluated at each rational surface", dims=("surface",), attach=(1 => "SingularSurfaces/rational_psi", 1 => "SingularSurfaces/rational_q")),
     "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",), attach=(1 => "SingularSurfaces/rational_psi", 1 => "SingularSurfaces/rational_q")),
     "SingularSurfaces/F" =>

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

test/runtests_h5_schema.jl

@@ -95,6 +95,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