ALL - REFACTOR - Module-mirroring CamelCase HDF5 schema (issue #226, PR A)

CLAUDE.md

@@ -146,6 +146,10 @@ Additional file hygiene (enforced by pre-commit hooks):
 - Files must end with exactly one newline
 - LF line endings only (no CRLF)
 
+### HDF5 Output Conventions
+
+The `gpec.h5` schema follows one physics-first convention (CamelCase groups at all levels, snake_case datasets, data-driven tokens verbatim, inputs only under `Input/`, five named top-level physics-topic exceptions). **Do not invent new group names or echo inputs into output groups** — read **[`docs/development/hdf5-conventions.md`](docs/development/hdf5-conventions.md)** before adding or moving any HDF5 output; renames are clean breaks (update writer, readers, and harness case TOMLs together — there is no legacy-path shim).
+
 ### TOML Annotation Conventions
 
 Config-style TOML files (`examples/*/gpec.toml`, `examples/*/sol.toml`, `test/test_data/*` fixtures, `regression-harness/cases/*.toml`) follow one shared annotation style (header comment block, inline `# description` on every variable line sourced from the matching config struct's docstring, no Fortran references, no deprecated variables). **Do not invent a new convention** — read **[`docs/development/toml-conventions.md`](docs/development/toml-conventions.md)** in full before adding or editing one of these files.

benchmarks/benchmark_against_fortran_run.jl

@@ -312,17 +312,17 @@ function load_julia_outputs(h5_path::String)
 
     julia = Dict{String,Any}()
     h5open(h5_path, "r") do f
-        julia["psilim"] = read(f, "info/psilim")
-        julia["qlim"] = read(f, "info/qlim")
-        julia["et"] = read(f, "vacuum/et")
-        julia["psi_q"] = read(f, "splines/profiles/xs")
-        julia["q"] = read(f, "splines/profiles/q")
-        julia["di"] = haskey(f, "locstab/di") ? read(f, "locstab/di") : Float64[]
-        julia["dr"] = haskey(f, "locstab/dr") ? read(f, "locstab/dr") : Float64[]
+        julia["psilim"] = read(f, "Info/psilim")
+        julia["qlim"] = read(f, "Info/qlim")
+        julia["et"] = read(f, "ForceFreeStates/FreeBoundaryStability/eigenmode_energies")
+        julia["psi_q"] = read(f, "Equilibrium/Profiles/xs")
+        julia["q"] = read(f, "Equilibrium/Profiles/q")
+        julia["di"] = haskey(f, "LocalStability/di") ? read(f, "LocalStability/di") : Float64[]
+        julia["dr"] = haskey(f, "LocalStability/dr") ? read(f, "LocalStability/dr") : Float64[]
 
-        julia["psio"] = haskey(f, "equil/psio") ? read(f, "equil/psio") : NaN
+        julia["psio"] = haskey(f, "Equilibrium/psio") ? read(f, "Equilibrium/psio") : NaN
 
-        sc = "perturbed_equilibrium/singular_coupling"
+        sc = "PerturbedEquilibrium/SingularCoupling"
         julia["rational_psi"] = haskey(f, "$sc/rational_psi") ? read(f, "$sc/rational_psi") : Float64[]
         julia["rational_q"] = haskey(f, "$sc/rational_q") ? read(f, "$sc/rational_q") : Float64[]
         julia["rational_n"] = haskey(f, "$sc/rational_n") ? read(f, "$sc/rational_n") : Int[]
@@ -333,30 +333,31 @@ function load_julia_outputs(h5_path::String)
         julia["chirikov_parameter"] = haskey(f, "$sc/chirikov_parameter") ? read(f, "$sc/chirikov_parameter") : Float64[]
         julia["delta_prime"] = haskey(f, "$sc/delta_prime") ? read(f, "$sc/delta_prime") : ComplexF64[]
 
-        pe = "perturbed_equilibrium"
+        pe = "PerturbedEquilibrium"
         # Fortran Phi_x/Phi_tot are the area-weighted field b̄ (tesla), matching forcing/response_b_area directly.
         julia["forcing_vec"] = haskey(f, "$pe/forcing_b_area") ? read(f, "$pe/forcing_b_area") : ComplexF64[]   # Phi_x   = b̄_x
         julia["response_vec"] = haskey(f, "$pe/response_b_area") ? read(f, "$pe/response_b_area") : ComplexF64[]   # Phi_tot = b̄_tot
-        julia["b_n"] = haskey(f, "$pe/response/b_n") ? read(f, "$pe/response/b_n") : Matrix{ComplexF64}(undef, 0, 0)
-        julia["Jbgradpsi"] = haskey(f, "$pe/response/b_psi_area_weighted") ? read(f, "$pe/response/b_psi_area_weighted") : Matrix{ComplexF64}(undef, 0, 0)
-        julia["xi_psi"] = haskey(f, "$pe/response/xi_psi") ? read(f, "$pe/response/xi_psi") : Matrix{ComplexF64}(undef, 0, 0)
-        julia["xi_n"] = haskey(f, "$pe/response/xi_n") ? read(f, "$pe/response/xi_n") : Matrix{ComplexF64}(undef, 0, 0)
-        julia["clebsch_psi1"] = haskey(f, "$pe/response/clebsch_psi1") ? read(f, "$pe/response/clebsch_psi1") : Matrix{ComplexF64}(undef, 0, 0)
-        julia["clebsch_alpha"] = haskey(f, "$pe/response/clebsch_alpha") ? read(f, "$pe/response/clebsch_alpha") : Matrix{ComplexF64}(undef, 0, 0)
-        julia["psi_grid"] = haskey(f, "integration/psi") ? read(f, "integration/psi") : Float64[]
+        julia["b_n"] = haskey(f, "$pe/Response/b_n") ? read(f, "$pe/Response/b_n") : Matrix{ComplexF64}(undef, 0, 0)
+        julia["Jbgradpsi"] = haskey(f, "$pe/Response/b_psi_area_weighted") ? read(f, "$pe/Response/b_psi_area_weighted") : Matrix{ComplexF64}(undef, 0, 0)
+        julia["xi_psi"] = haskey(f, "$pe/Response/xi_psi") ? read(f, "$pe/Response/xi_psi") : Matrix{ComplexF64}(undef, 0, 0)
+        julia["xi_n"] = haskey(f, "$pe/Response/xi_n") ? read(f, "$pe/Response/xi_n") : Matrix{ComplexF64}(undef, 0, 0)
+        julia["clebsch_psi1"] = haskey(f, "$pe/Response/clebsch_psi1") ? read(f, "$pe/Response/clebsch_psi1") : Matrix{ComplexF64}(undef, 0, 0)
+        julia["clebsch_alpha"] = haskey(f, "$pe/Response/clebsch_alpha") ? read(f, "$pe/Response/clebsch_alpha") : Matrix{ComplexF64}(undef, 0, 0)
+        julia["psi_grid"] = haskey(f, "ForceFreeStates/Solutions/ForwardIntegration/psi") ? read(f, "ForceFreeStates/Solutions/ForwardIntegration/psi") : Float64[]
 
         # R,Z,φ: loaded via modes_to_theta helper below (not raw modes)
         julia["h5_path"] = h5_path  # stash for modes_to_theta
         # mn_index[:, 1] = m values, mn_index[:, 2] = n values for each mode index
-        julia["m_modes"] = haskey(f, "info/mn_index") ? Int.(read(f, "info/mn_index")[:, 1]) : Int[]
+        julia["m_modes"] = haskey(f, "Info/mn_index") ? Int.(read(f, "Info/mn_index")[:, 1]) : Int[]
 
         # Control surface matrices
-        rm = "$pe/response_matrices"
+        rm = "$pe/ResponseMatrices"
         julia["permeability"] = haskey(f, "$rm/permeability") ? read(f, "$rm/permeability") : Matrix{ComplexF64}(undef, 0, 0)
         julia["plasma_inductance"] = haskey(f, "$rm/plasma_inductance") ? read(f, "$rm/plasma_inductance") : Matrix{ComplexF64}(undef, 0, 0)
         julia["surface_inductance"] = haskey(f, "$rm/surface_inductance") ? read(f, "$rm/surface_inductance") : Matrix{ComplexF64}(undef, 0, 0)
         julia["reluctance"] = haskey(f, "$rm/reluctance") ? read(f, "$rm/reluctance") : Matrix{ComplexF64}(undef, 0, 0)
-        julia["wt0"] = haskey(f, "vacuum/wt0") ? read(f, "vacuum/wt0") : Matrix{ComplexF64}(undef, 0, 0)
+        fbs = "ForceFreeStates/FreeBoundaryStability"
+        julia["wt0"] = haskey(f, "$fbs/W_freeboundary") ? read(f, "$fbs/W_freeboundary") : Matrix{ComplexF64}(undef, 0, 0)
     end
     return julia
 end
@@ -1194,12 +1195,12 @@ function generate_plots(fort, julia, bench_dir, nn)
     h5_path = julia["h5_path"]
     rzphi_panels = []
     for (comp_label, f_key, h5_var) in [
-        ("ξ_R", "xi_r_fun", "perturbed_equilibrium/response/xi_R"),
-        ("ξ_Z", "xi_z_fun", "perturbed_equilibrium/response/xi_Z"),
-        ("ξ_φ", "xi_phi_fun", "perturbed_equilibrium/response/xi_phi"),
-        ("b_R", "b_r_fun", "perturbed_equilibrium/response/b_R"),
-        ("b_Z", "b_z_fun", "perturbed_equilibrium/response/b_Z"),
-        ("b_φ", "b_phi_fun", "perturbed_equilibrium/response/b_phi")]
+        ("ξ_R", "xi_r_fun", "PerturbedEquilibrium/Response/xi_R"),
+        ("ξ_Z", "xi_z_fun", "PerturbedEquilibrium/Response/xi_Z"),
+        ("ξ_φ", "xi_phi_fun", "PerturbedEquilibrium/Response/xi_phi"),
+        ("b_R", "b_r_fun", "PerturbedEquilibrium/Response/b_R"),
+        ("b_Z", "b_z_fun", "PerturbedEquilibrium/Response/b_Z"),
+        ("b_φ", "b_phi_fun", "PerturbedEquilibrium/Response/b_phi")]
         f_fun = get(fort, f_key, Matrix{ComplexF64}(undef, 0, 0))
 
         # Reconstruct theta-space from mode-space via modes_to_theta (applies ν phase + helicity)

benchmarks/benchmark_diiid_kinetic_stability.jl

@@ -8,7 +8,7 @@ Fortran GPEC's kinetic DCON reference.
 Runs `GPE.main()` with `kinetic_source="calculated"` and `kinetic_factor=1.0`
 against the EFIT g-file and `.kin` profile taken from a Fortran GPEC
 kinetic example directory, then compares the least-stable total-energy
-eigenvalue `vacuum/et[1]` against `W_t_eigenvalue[:, 0]` in the Fortran
+eigenvalue `ForceFreeStates/FreeBoundaryStability/eigenmode_energies[1]` against `W_t_eigenvalue[:, 0]` in the Fortran
 `dcon_output_n1.nc`.
 
 No inputs are duplicated into this repo — everything is read from the
@@ -202,9 +202,9 @@ function run_benchmark(fortran_dir::String=default_fortran_dir())
     isfile(h5path) || error("Expected Julia output not found: $h5path")
 
     et = h5open(h5path, "r") do h5
-        read(h5["vacuum/et"])
+        read(h5["ForceFreeStates/FreeBoundaryStability/eigenmode_energies"])
     end
-    isempty(et) && error("vacuum/et is empty in $h5path")
+    isempty(et) && error("ForceFreeStates/FreeBoundaryStability/eigenmode_energies is empty in $h5path")
     # et is stored as a length-2*N real array (re,im interleaved) by HDF5.jl
     # when the underlying Julia array is ComplexF64. NCDatasets and HDF5 give
     # us a ComplexF64 array directly here.

benchmarks/benchmark_git_branches.jl

@@ -194,8 +194,8 @@ function run_example_benchmark(example_path, num_runs)
         end
 
         h5 = h5open(gpec_path, "r")
-        et = read(h5["FreeBoundaryStability/eigenmode_energies"])
-        nsteps = read(h5["integration/nstep"])
+        et = read(h5["ForceFreeStates/FreeBoundaryStability/eigenmode_energies"])
+        nsteps = read(h5["ForceFreeStates/Solutions/ForwardIntegration/nstep"])
         close(h5)
 
         avg_runtime = sum(runtimes) / length(runtimes)

benchmarks/benchmark_solovev_kinetic_stability.jl

@@ -155,7 +155,7 @@ function run_julia_reference()
     GPE.main([rundir])
     wall = time() - t0
     et = h5open(joinpath(rundir, "gpec.h5"), "r") do h5
-        read(h5["FreeBoundaryStability/eigenmode_energies"])
+        read(h5["ForceFreeStates/FreeBoundaryStability/eigenmode_energies"])
     end
     return real(et[1]), imag(et[1]), wall
 end

benchmarks/compare_gal_vs_el.jl

@@ -1,7 +1,7 @@
 # Overlay the IDEAL gal matched ξ(ψ) against the EL total-energy eigenmode ξ(ψ), for the same edge
 # eigenvector w. The ideal gal solution should reproduce the EL (DCON) ideal solution.
 #
-# EL :  ξ_EL(ψ)  = U_EL(ψ) · (U_EL_edge \ w)          (fundamental matrix, integration/xi_psi)
+# EL :  ξ_EL(ψ)  = U_EL(ψ) · (U_EL_edge \ w)          (fundamental matrix, ForceFreeStates/Solutions/ForwardIntegration/xi_psi)
 # gal:  ξ_gal(ψ) = U_gal(ψ) · w                       (identity-at-edge ⇒ coefficient is w itself)
 # w = eigenvector of the total energy operator W = W_plasma + W_vacuum (FreeBoundaryStability).
 #
@@ -19,12 +19,12 @@ ksel = length(ARGS) >= 3 ? ARGS[3] : "highest"
 to_c(a) = eltype(a) <: Complex ? ComplexF64.(a) : map(x -> ComplexF64(x.re, x.im), a)
 
 et, wt, u1, psiE, gxi, psiG, issing, mlow, sing_psi = h5open(h5path) do f
-    (to_c(read(f["FreeBoundaryStability/eigenmode_energies"])),
-        to_c(read(f["FreeBoundaryStability/W_freeboundary_eigenmodes"])),
-        to_c(read(f["integration/xi_psi"])), read(f["integration/psi"]),
-        to_c(read(f["galerkin/match/xi"])), read(f["galerkin/solution/psi"]),
-        Bool.(read(f["galerkin/solution/issing"])), read(f["info/mlow"]),
-        read(f["galerkin/sing_psi"]))
+    (to_c(read(f["ForceFreeStates/FreeBoundaryStability/eigenmode_energies"])),
+        to_c(read(f["ForceFreeStates/FreeBoundaryStability/W_freeboundary_eigenmodes"])),
+        to_c(read(f["ForceFreeStates/Solutions/ForwardIntegration/xi_psi"])), read(f["ForceFreeStates/Solutions/ForwardIntegration/psi"]),
+        to_c(read(f["ForceFreeStates/Solutions/GalerkinIntegration/Match/xi"])), read(f["ForceFreeStates/Solutions/GalerkinIntegration/Solution/psi"]),
+        Bool.(read(f["ForceFreeStates/Solutions/GalerkinIntegration/Solution/issing"])), read(f["Info/mlow"]),
+        read(f["SingularSurfaces/GalerkinDeltaPrime/sing_psi"]))
 end
 
 mpert = size(u1, 1)

benchmarks/compare_jbgradpsi_m2.jl

@@ -1,10 +1,10 @@
-# Compare the area-normalized b^ψ (perturbed_equilibrium/response/psi_area = b^ψ/⟨J·|∇ψ|⟩_θ) for one
+# Compare the area-normalized b^ψ (PerturbedEquilibrium/Response/psi_area = b^ψ/⟨J·|∇ψ|⟩_θ) for one
 # poloidal harmonic between two GPEC runs that are identical except for which ξ feeds PerturbedEquilibrium:
 #   (1) IDEAL galerkin matched ξ   (gal_match_flag=true, gal_ideal_flag=true)
 #   (2) SHOOTING ξ                  (gal_match_flag=false)
 #
-# PE writes no ψ grid, so it's reconstructed: gal-ideal → galerkin/solution/psi minus issing points;
-# shooting → integration/psi.
+# PE writes no ψ grid, so it's reconstructed: gal-ideal → ForceFreeStates/Solutions/GalerkinIntegration/Solution/psi minus issing points;
+# shooting → ForceFreeStates/Solutions/ForwardIntegration/psi.
 # Usage: julia --project=. benchmarks/compare_jbgradpsi_m2.jl [gal_h5] [shoot_h5] [out.png] [m]
 
 using HDF5, Plots, Printf
@@ -18,14 +18,14 @@ to_c(a) = eltype(a) <: Complex ? ComplexF64.(a) : map(x -> ComplexF64(x.re, x.im
 
 # gal-ideal run: PE grid = gal solution grid with the on-surface (issing) points dropped
 pa_g, psi_g, mlow, sing_psi, sing_m = h5open(gal_h5) do f
-    pa = to_c(read(f["perturbed_equilibrium/response/psi_area"]))   # [npsi, mpert]
-    iss = Bool.(read(f["galerkin/solution/issing"]))
-    (pa, read(f["galerkin/solution/psi"])[.!iss], read(f["info/mlow"]),
-        read(f["galerkin/sing_psi"]), read(f["galerkin/sing_m"]))
+    pa = to_c(read(f["PerturbedEquilibrium/Response/psi_area"]))   # [npsi, mpert]
+    iss = Bool.(read(f["ForceFreeStates/Solutions/GalerkinIntegration/Solution/issing"]))
+    (pa, read(f["ForceFreeStates/Solutions/GalerkinIntegration/Solution/psi"])[.!iss], read(f["Info/mlow"]),
+        read(f["SingularSurfaces/GalerkinDeltaPrime/sing_psi"]), read(f["SingularSurfaces/GalerkinDeltaPrime/sing_m"]))
 end
-# shooting run: PE grid = integration/psi
+# shooting run: PE grid = ForceFreeStates/Solutions/ForwardIntegration/psi
 pa_s, psi_s = h5open(sh_h5) do f
-    (to_c(read(f["perturbed_equilibrium/response/psi_area"])), read(f["integration/psi"]))
+    (to_c(read(f["PerturbedEquilibrium/Response/psi_area"])), read(f["ForceFreeStates/Solutions/ForwardIntegration/psi"]))
 end
 
 size(pa_g, 1) == length(psi_g) || error("gal grid mismatch: npsi=$(size(pa_g,1)) vs grid=$(length(psi_g))")

benchmarks/equil_psihigh_scan.jl

@@ -69,7 +69,7 @@ function run_ffs_et1(config_path::String, eq_type::String, psihigh::Float64)::Fl
         try
             GeneralizedPerturbedEquilibrium.main([tmpdir])
             h5open(joinpath(tmpdir, "gpec.h5"), "r") do h5
-                et = read(h5["FreeBoundaryStability/eigenmode_energies"])
+                et = read(h5["ForceFreeStates/FreeBoundaryStability/eigenmode_energies"])
                 return real(et[1])
             end
         catch

benchmarks/plot_xi_eigenmode.jl

@@ -3,9 +3,9 @@
 # The total energy operator W = W_plasma + W_vacuum (free_run, Free.jl); its eigenvectors are the
 # free-boundary edge displacement patterns and the eigenvalues are δW. This picks the eigenmode with the
 # largest Re(eigenvalue) and reconstructs its radial profile by projecting the EL fundamental matrix
-# (integration/xi_psi) onto that edge eigenvector:  c = U_edge \ w,  ξ(ψ) = U(ψ)·c.
+# (ForceFreeStates/Solutions/ForwardIntegration/xi_psi) onto that edge eigenvector:  c = U_edge \ w,  ξ(ψ) = U(ψ)·c.
 #
-# Requires a run with populate_dense_xi=true so integration/xi_psi is the dense axis-basis fundamental
+# Requires a run with populate_dense_xi=true so ForceFreeStates/Solutions/ForwardIntegration/xi_psi is the dense axis-basis fundamental
 # matrix (not the Riccati S-matrices).
 # Usage: julia --project=. benchmarks/plot_xi_eigenmode.jl [path/to/gpec.h5] [out.png]
 
@@ -17,12 +17,12 @@ outpng = length(ARGS) >= 2 ? ARGS[2] : joinpath(@__DIR__, "xi_eigenmode.png")
 to_c(a) = eltype(a) <: Complex ? ComplexF64.(a) : map(x -> ComplexF64(x.re, x.im), a)
 
 et, wt, u1, psi, mlow, sing_psi = h5open(h5path) do f
-    (to_c(read(f["FreeBoundaryStability/eigenmode_energies"])),
-        to_c(read(f["FreeBoundaryStability/W_freeboundary_eigenmodes"])),
-        to_c(read(f["integration/xi_psi"])),
-        read(f["integration/psi"]),
-        read(f["info/mlow"]),
-        haskey(f, "galerkin/sing_psi") ? read(f["galerkin/sing_psi"]) : Float64[])
+    (to_c(read(f["ForceFreeStates/FreeBoundaryStability/eigenmode_energies"])),
+        to_c(read(f["ForceFreeStates/FreeBoundaryStability/W_freeboundary_eigenmodes"])),
+        to_c(read(f["ForceFreeStates/Solutions/ForwardIntegration/xi_psi"])),
+        read(f["ForceFreeStates/Solutions/ForwardIntegration/psi"]),
+        read(f["Info/mlow"]),
+        haskey(f, "SingularSurfaces/GalerkinDeltaPrime/sing_psi") ? read(f["SingularSurfaces/GalerkinDeltaPrime/sing_psi"]) : Float64[])
 end
 
 mpert, _, nstep = size(u1)

benchmarks/scan_resistivity_m2.jl

@@ -1,4 +1,4 @@
-# Plot the m=2 area-normalized b^ψ (perturbed_equilibrium/response/psi_area) across a resistivity scan
+# Plot the m=2 area-normalized b^ψ (PerturbedEquilibrium/Response/psi_area) across a resistivity scan
 # of the RESISTIVE gal matched PE runs (gal_match_flag=true, gal_ideal_flag=false), one curve per η.
 # Overlays the shooting (ideal, η→0) reference. The η-scan dirs are produced by the bash loop over
 # /tmp/etascan_<factor> (each a copy of the 0.993 config with gal_eta scaled).
@@ -14,10 +14,10 @@ to_c(a) = eltype(a) <: Complex ? ComplexF64.(a) : map(x -> ComplexF64(x.re, x.im
 # read m=target area-normalized b^ψ on the run's PE grid
 function read_m2(h5; gal::Bool)
     h5open(h5) do f
-        pa = to_c(read(f["perturbed_equilibrium/response/psi_area"]))   # [npsi, mpert]
-        col = mtarget - read(f["info/mlow"]) + 1
-        psi = gal ? read(f["galerkin/solution/psi"])[.!Bool.(read(f["galerkin/solution/issing"]))] :
-              read(f["integration/psi"])
+        pa = to_c(read(f["PerturbedEquilibrium/Response/psi_area"]))   # [npsi, mpert]
+        col = mtarget - read(f["Info/mlow"]) + 1
+        psi = gal ? read(f["ForceFreeStates/Solutions/GalerkinIntegration/Solution/psi"])[.!Bool.(read(f["ForceFreeStates/Solutions/GalerkinIntegration/Solution/issing"]))] :
+              read(f["ForceFreeStates/Solutions/ForwardIntegration/psi"])
         (psi, pa[:, col])
     end
 end
@@ -34,7 +34,7 @@ eta_ref = 8e-8
 
 # rational surface for m=target
 sing_psi, sing_m = h5open(joinpath(scandirs[1], "gpec.h5")) do f
-    (read(f["galerkin/sing_psi"]), read(f["galerkin/sing_m"]))
+    (read(f["SingularSurfaces/GalerkinDeltaPrime/sing_psi"]), read(f["SingularSurfaces/GalerkinDeltaPrime/sing_m"]))
 end
 psi_res = mtarget in sing_m ? sing_psi[findfirst(==(mtarget), sing_m)] : NaN