diff --git a/docs/development/hdf5-conventions.md b/docs/development/hdf5-conventions.md index c4631c984..a08f7efd4 100644 --- a/docs/development/hdf5-conventions.md +++ b/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/}` | | `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 | diff --git a/docs/src/stability.md b/docs/src/stability.md index c61ca716e..2dba52cfa 100644 --- a/docs/src/stability.md +++ b/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` diff --git a/docs/src/workflow.md b/docs/src/workflow.md index 86a7ae0df..79978d7b9 100644 --- a/docs/src/workflow.md +++ b/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/` | diff --git a/src/ForceFreeStates/CoreTypes.jl b/src/ForceFreeStates/CoreTypes.jl index 5544b0e60..49a1d96a4 100644 --- a/src/ForceFreeStates/CoreTypes.jl +++ b/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 diff --git a/src/GeneralizedPerturbedEquilibrium.jl b/src/GeneralizedPerturbedEquilibrium.jl index 23de76dad..e288a4a19 100755 --- a/src/GeneralizedPerturbedEquilibrium.jl +++ b/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). diff --git a/test/runtests_fullruns.jl b/test/runtests_fullruns.jl index d652bb35b..297044339 100644 --- a/test/runtests_fullruns.jl +++ b/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 diff --git a/test/runtests_h5_schema.jl b/test/runtests_h5_schema.jl index c743d6a7d..2a9967906 100644 --- a/test/runtests_h5_schema.jl +++ b/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")