HDF5Schema - API! - ⚠️ Make the Euler-Lagrange matrix output optional and off by default

docs/development/hdf5-conventions.md

@@ -40,7 +40,7 @@ Top level (11 groups):
 | `Info/` | Run metadata: `git_version`, mode-number ranges (`mpert`, `mlow`, …, `mn_index`), `psilim`, `qlim`, `Runtimes/` (per-stage wall-clock seconds) |
 | `Input/` | Rerun snapshot: `gpec_toml_raw`, `RawInputs/{Equilibrium, ForcingTerms, Coils/<name>}` |
 | `Equilibrium/` | Scalars (`beta_N`, `q_axis`, `q_95`, `I_p`, …) plus `Profiles/` (1-D on `psi`: 2piF, mu0p, dVdpsi, q) and `Geometry/` (2-D on `psi`×`theta`: rcoords, offset, nu, jac) |
-| `ForceFreeStates/` | `Solutions/ForwardIntegration/` (u-solutions), `Solutions/GalerkinIntegration/` (closed ξ profiles in the shared layout, `Match/` diagnostics, the gal surface list, debug-gated `Basis/`), `EulerLagrangeMatrices/{Ideal,Kinetic}`, `FreeBoundaryStability/`, `EdgeScan/` |
+| `ForceFreeStates/` | `Solutions/ForwardIntegration/` (u-solutions), `Solutions/GalerkinIntegration/` (closed ξ profiles in the shared layout, `Match/` diagnostics, the gal surface list, debug-gated `Basis/`), `write_el_matrices`-gated `EulerLagrangeMatrices/{Ideal,Kinetic}`, `FreeBoundaryStability/`, `EdgeScan/` |
 | `LocalStability/` | Mercier `D_I`, resistive interchange `D_R`, `ballooning_Delta_prime` on `psi`; the ballooning α boundary on `ballooning_psi` |
 | `SingularSurfaces/` | Per-rational-surface data: `rational_psi`/`rational_q`/`rational_m`/`rational_n`, GGJ coefficients, `Delta_prime_matrix`/`Delta_prime_raw`/`Delta_coil`/`pest3_A`/`pest3_B`/`pest3_Gamma` (Riccati or Galerkin alike), `Kinetic/` |
 | `PerturbedEquilibrium/` | `ForcingModes/`, `Response/`, `ResponseMatrices/`, `SingularCoupling/`, `Energies/`, control-surface spectra |

docs/src/stability.md

@@ -280,6 +280,7 @@ local_stability_flag = false  # scan Mercier D_I, resistive D_R, and ballooning
 # Output
 verbose              = true
 write_outputs_to_HDF5 = true
+write_el_matrices    = false  # also write ForceFreeStates/EulerLagrangeMatrices (mpert² × npsi per matrix)
 ```
 
 The number of Julia threads is controlled at startup via `-t N` or the `JULIA_NUM_THREADS`

docs/src/workflow.md

@@ -178,7 +178,7 @@ All results are written to a single HDF5 file (default: `gpec.h5`). The top-leve
 | `Info/` | Run metadata: git version, mode-number ranges, ψ limit, `Runtimes/` (per-stage wall-clock seconds) |
 | `Input/` | Self-contained rerun snapshot: merged TOML blob, raw equilibrium/forcing/coil inputs |
 | `Equilibrium/` | Equilibrium scalars (`beta_N`, `q_axis`, `q_95`, …), 1-D profiles (`Profiles/`), 2-D geometry (`Geometry/`) |
-| `ForceFreeStates/` | Stability solve: `Solutions/{ForwardIntegration,GalerkinIntegration}`, `EulerLagrangeMatrices/`, `FreeBoundaryStability/`, `EdgeScan/` |
+| `ForceFreeStates/` | Stability solve: `Solutions/{ForwardIntegration,GalerkinIntegration}`, `EulerLagrangeMatrices/` (opt-in: `write_el_matrices = true`), `FreeBoundaryStability/`, `EdgeScan/` |
 | `LocalStability/` | Mercier D_I, resistive interchange D_R, ballooning Δ' profiles |
 | `SingularSurfaces/` | Per-rational-surface data: ψ_s, q, m/n, GGJ coefficients, Δ'/PEST-3 matching matrices, kinetic surfaces (`Kinetic/`) |
 | `PerturbedEquilibrium/` | Plasma response: `ForcingModes/`, `Response/`, `ResponseMatrices/`, `SingularCoupling/`, `Energies/` |

src/ForceFreeStates/CoreTypes.jl

@@ -147,6 +147,7 @@ gpec.toml.
   - `diagnose::Bool` - Enable diagnostic output (not yet implemented)
   - `diagnose_ca::Bool` - Enable asymptotic coefficient diagnostics (not yet implemented)
   - `write_outputs_to_HDF5::Bool` - Write results to HDF5 format
+  - `write_el_matrices::Bool` - Write the Euler-Lagrange matrices (`ForceFreeStates/EulerLagrangeMatrices`: the ideal A-K and, on a kinetic run, the kinetic set) sampled on the ψ grid. Each is `mpert × mpert × npsi` complex, so the group scales as mpert²·npsi and dominates the file size. Default `false`.
   - `HDF5_filename::String` - Name of HDF5 output file
   - `save_interval::Int` - Save every Nth ODE step (1=all). Always saves near rational surfaces. Default `1`: PerturbedEquilibrium and KineticForces interpolate ξ(ψ) between saved steps, so their accuracy follows the saved-step density.
   - `force_termination::Bool` - Terminate after force-free states (skip perturbed equilibrium calculations)
@@ -184,6 +185,7 @@ gpec.toml.
     diagnose::Bool = false
     diagnose_ca::Bool = false
     write_outputs_to_HDF5::Bool = true
+    write_el_matrices::Bool = false
     HDF5_filename::String = "gpec.h5"
     save_interval::Int = 1
     force_termination::Bool = false

src/GeneralizedPerturbedEquilibrium.jl

@@ -1647,45 +1647,47 @@ function write_outputs_to_HDF5(
         out_h5["SurfaceGeometries/Wall/y"] = free_energies !== nothing ? free_energies.wall_pts[:, 2] : Float64[]
         out_h5["SurfaceGeometries/Wall/z"] = free_energies !== nothing ? free_energies.wall_pts[:, 3] : Float64[]
 
-        # Write fundamental matrices on the ψ grid
-        xs = equil.rzphi_xs
-        npsi = length(xs)
-        np = result.numpert_total
-
-        # Helper: evaluate a matrix spline on the psi grid → (npsi, np, np) array
-        function _eval_mat_spline(spline)
-            arr = zeros(ComplexF64, npsi, np, np)
-            hint = Ref(1)
-            for i in 1:npsi
-                arr[i, :, :] .= reshape(spline(xs[i]; hint=hint), np, np)
+        # Euler-Lagrange matrices on the ψ grid; opt-in since the group scales as mpert²·npsi.
+        if ctrl.write_el_matrices
+            xs = equil.rzphi_xs
+            npsi = length(xs)
+            np = result.numpert_total
+
+            # Helper: evaluate a matrix spline on the psi grid → (npsi, np, np) array
+            function _eval_mat_spline(spline)
+                arr = zeros(ComplexF64, npsi, np, np)
+                hint = Ref(1)
+                for i in 1:npsi
+                    arr[i, :, :] .= reshape(spline(xs[i]; hint=hint), np, np)
+                end
+                return arr
             end
-            return arr
-        end
 
-        elm = "ForceFreeStates/EulerLagrangeMatrices"
-        out_h5["$elm/psi"] = xs
-        # Ideal primitive matrices (A, B, C, D, E, H)
-        out_h5["$elm/Ideal/A"] = _eval_mat_spline(mats.ideal.A_spline)
-        out_h5["$elm/Ideal/B"] = _eval_mat_spline(mats.ideal.B_spline)
-        out_h5["$elm/Ideal/C"] = _eval_mat_spline(mats.ideal.C_spline)
-        out_h5["$elm/Ideal/D"] = _eval_mat_spline(mats.ideal.D_spline_prim)
-        out_h5["$elm/Ideal/E"] = _eval_mat_spline(mats.ideal.E_spline_prim)
-        out_h5["$elm/Ideal/H"] = _eval_mat_spline(mats.ideal.H_spline)
-
-        # Ideal derived matrices (F, K, G)
-        out_h5["$elm/Ideal/F"] = _eval_mat_spline(mats.ideal.F_spline_lower)
-        out_h5["$elm/Ideal/K"] = _eval_mat_spline(mats.ideal.K_spline)
-        out_h5["$elm/Ideal/G"] = _eval_mat_spline(mats.ideal.G_spline)
-
-        # Kinetic-modified matrices
-        kin = mats.kinetic
-        if kin !== nothing
-            out_h5["$elm/Kinetic/A"] = _eval_mat_spline(kin.A_spline)
-            out_h5["$elm/Kinetic/B"] = _eval_mat_spline(kin.B_spline)
-            out_h5["$elm/Kinetic/C"] = _eval_mat_spline(kin.C_spline)
-            out_h5["$elm/Kinetic/f0"] = _eval_mat_spline(kin.F0_spline)
-            out_h5["$elm/Kinetic/K"] = _eval_mat_spline(kin.Kk_spline)
-            out_h5["$elm/Kinetic/G"] = _eval_mat_spline(kin.G_spline_adj)
+            elm = "ForceFreeStates/EulerLagrangeMatrices"
+            out_h5["$elm/psi"] = xs
+            # Ideal primitive matrices (A, B, C, D, E, H)
+            out_h5["$elm/Ideal/A"] = _eval_mat_spline(mats.ideal.A_spline)
+            out_h5["$elm/Ideal/B"] = _eval_mat_spline(mats.ideal.B_spline)
+            out_h5["$elm/Ideal/C"] = _eval_mat_spline(mats.ideal.C_spline)
+            out_h5["$elm/Ideal/D"] = _eval_mat_spline(mats.ideal.D_spline_prim)
+            out_h5["$elm/Ideal/E"] = _eval_mat_spline(mats.ideal.E_spline_prim)
+            out_h5["$elm/Ideal/H"] = _eval_mat_spline(mats.ideal.H_spline)
+
+            # Ideal derived matrices (F, K, G)
+            out_h5["$elm/Ideal/F"] = _eval_mat_spline(mats.ideal.F_spline_lower)
+            out_h5["$elm/Ideal/K"] = _eval_mat_spline(mats.ideal.K_spline)
+            out_h5["$elm/Ideal/G"] = _eval_mat_spline(mats.ideal.G_spline)
+
+            # Kinetic-modified matrices
+            kin = mats.kinetic
+            if kin !== nothing
+                out_h5["$elm/Kinetic/A"] = _eval_mat_spline(kin.A_spline)
+                out_h5["$elm/Kinetic/B"] = _eval_mat_spline(kin.B_spline)
+                out_h5["$elm/Kinetic/C"] = _eval_mat_spline(kin.C_spline)
+                out_h5["$elm/Kinetic/f0"] = _eval_mat_spline(kin.F0_spline)
+                out_h5["$elm/Kinetic/K"] = _eval_mat_spline(kin.Kk_spline)
+                out_h5["$elm/Kinetic/G"] = _eval_mat_spline(kin.G_spline_adj)
+            end
         end
 
         # Self-describing metadata pass (long_name/units/dims + dimension scales).

test/runtests_fullruns.jl

@@ -35,6 +35,8 @@ using HDF5
             # emit deterministic zero-extent sentinels, not uninitialized memory.
             @test isempty(read(h5["SingularSurfaces/ca_left"]))
             @test isempty(read(h5["SingularSurfaces/ca_right"]))
+            # The Euler-Lagrange matrices are opt-in; a default deck writes no group.
+            @test !haskey(h5, "ForceFreeStates/EulerLagrangeMatrices")
         end
         rm(joinpath(ex3, "gpec.h5"); force=true)
         true

test/runtests_h5_schema.jl

@@ -58,7 +58,8 @@ end
             cp(joinpath(template_dir, name), joinpath(run_dir, name))
         end
         toml_path = joinpath(run_dir, "gpec.toml")
-        write(toml_path, replace(read(toml_path, String), "write_outputs_to_HDF5 = false" => "write_outputs_to_HDF5 = true"))
+        # Opt into the Euler-Lagrange matrices so the metadata pass below covers that group too.
+        write(toml_path, replace(read(toml_path, String), "write_outputs_to_HDF5 = false" => "write_el_matrices = true\nwrite_outputs_to_HDF5 = true"))
 
         GeneralizedPerturbedEquilibrium.main([run_dir])
         h5_path = joinpath(run_dir, "gpec.h5")
@@ -74,6 +75,9 @@ end
             @test !haskey(h5, "FreeBoundaryStability")
             @test !haskey(h5, "EdgeScan")
 
+            # The opt-in above writes the Euler-Lagrange matrix group.
+            @test haskey(h5, "ForceFreeStates/EulerLagrangeMatrices/Ideal/A")
+
             # Inputs live only under Input/; spot-check the rerun-critical paths.
             @test haskey(h5, "Input/gpec_toml_raw")
             @test haskey(h5, "Info/git_version")