GPEC - NEW FEATURE - Per-stage runtime records in gpec.h5 under Info/Runtimes

benchmarks/benchmark_git_branches.jl

@@ -43,6 +43,7 @@ For each branch/commit, reports:
 - Eigenmode energy (et[1]) from gpec.h5
 - Integration steps from gpec.h5
 - Runtime (averaged over warm runs)
+- Per-stage runtime breakdown from `Info/Runtimes` (averaged over warm runs, when both refs record it)
 - Git commit hash
 
 # Example Output
@@ -67,6 +68,12 @@ Comparison:
   Eigenmode Δ:    +0.0067 (+0.4%)
   Steps Δ:        -775 (-85.1%)
   Runtime Δ:      -2.5 s (-1.7% faster)
+
+Per-stage runtime (avg over warm runs):
+  Stage                    Branch 1   Branch 2        Δ s       Δ %
+  equilibrium                 2.10s      2.08s     -0.02     -1.0%
+  force_free_states         100.40s     98.10s     -2.30     -2.3%
+  total                     145.30s    142.80s     -2.50     -1.7%
 ```
 """
 
@@ -156,6 +163,28 @@ function checkout_ref(branch, commit)
     end
 end
 
+# Read the per-stage wall-clock seconds a run recorded under Info/Runtimes.
+# Returns an empty Dict for refs that predate the feature, so those still benchmark.
+function read_stage_runtimes(gpec_path)
+    stages = Dict{String,Float64}()
+    isfile(gpec_path) || return stages
+    h5open(gpec_path, "r") do h5
+        haskey(h5, "Info/Runtimes") || return
+        group = h5["Info/Runtimes"]
+        for stage in keys(group)
+            stages[stage] = read(group[stage])
+        end
+    end
+    return stages
+end
+
+# Average per-stage runtimes over the warm runs, keeping only stages present in every run.
+function average_stage_runtimes(stage_times)
+    isempty(stage_times) && return Dict{String,Float64}()
+    shared = intersect(map(keys, stage_times)...)
+    return Dict(stage => sum(s[stage] for s in stage_times) / length(stage_times) for stage in shared)
+end
+
 # Run the example benchmark.
 # Each run is a fresh Julia subprocess so that git branch switches take effect
 # (using statement only fires once per process; re-checkout without subprocess restart
@@ -176,19 +205,22 @@ function run_example_benchmark(example_path, num_runs)
         println("\n[1/$(num_runs+1)] First run (JIT compilation)...")
         run(`julia --project=$project_root $tmpscript $abs_example_path`)
 
-        # Warm runs for timing
+        # Warm runs for timing. Each run overwrites gpec.h5, so the per-stage
+        # Info/Runtimes record must be collected inside the loop, not after it.
+        gpec_path = joinpath(abs_example_path, "gpec.h5")
         runtimes = Float64[]
+        stage_times = Vector{Dict{String,Float64}}()
         for i in 1:num_runs
             println("\n[$((i+1))/$(num_runs+1)] Warm run $i...")
             t_start = time()
             run(`julia --project=$project_root $tmpscript $abs_example_path`)
             runtime = time() - t_start
             push!(runtimes, runtime)
+            push!(stage_times, read_stage_runtimes(gpec_path))
             println("  Runtime: $(round(runtime, digits=2)) s")
         end
 
         # Extract metrics from gpec.h5
-        gpec_path = joinpath(abs_example_path, "gpec.h5")
         if !isfile(gpec_path)
             error("gpec.h5 not found at $gpec_path")
         end
@@ -203,13 +235,33 @@ function run_example_benchmark(example_path, num_runs)
         return (
             eigenvalue=real(et[1]),
             steps=nsteps,
-            runtime=avg_runtime
+            runtime=avg_runtime,
+            stages=average_stage_runtimes(stage_times)
         )
     finally
         rm(tmpscript; force=true)
     end
 end
 
+# Per-stage runtime breakdown (Info/Runtimes), averaged over the warm runs of each branch.
+# Silently skipped when either side lacks the record, so comparisons against older refs work.
+function print_stage_comparison(r1, r2)
+    shared = intersect(keys(r1.metrics.stages), keys(r2.metrics.stages))
+    isempty(shared) && return nothing
+    ordered = sort(collect(shared); by=stage -> (stage == "total", stage))
+
+    println("\nPer-stage runtime (avg over warm runs):")
+    @printf("  %-22s %10s %10s %10s %9s\n", "Stage", "Branch 1", "Branch 2", "Δ s", "Δ %")
+    for stage in ordered
+        t1 = r1.metrics.stages[stage]
+        t2 = r2.metrics.stages[stage]
+        # A stage that ran in no measurable time has no meaningful percentage baseline.
+        pct = t1 == 0 ? "—" : @sprintf("%+.1f%%", 100 * (t2 - t1) / t1)
+        @printf("  %-22s %9.2fs %9.2fs %+9.2f %9s\n", stage, t1, t2, t2 - t1, pct)
+    end
+    return nothing
+end
+
 # Main benchmarking function
 function benchmark_branches(options)
     println("="^60)
@@ -293,6 +345,8 @@ function benchmark_branches(options)
     speedup_desc = time_delta < 0 ? "faster" : "slower"
     @printf("  Runtime Δ:      %+.2f s (%+.1f%% %s)\n", time_delta, abs(time_pct), speedup_desc)
 
+    print_stage_comparison(r1, r2)
+
     println("\n" * "="^60)
 
     # Write to output file if requested

docs/development/hdf5-conventions.md

@@ -37,7 +37,7 @@ Top level (10 groups):
 
 | Group | Contents |
 |---|---|
-| `Info/` | Run metadata: `git_version`, mode-number ranges (`mpert`, `mlow`, …, `mn_index`), `psilim`, `qlim` |
+| `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/` |

docs/src/workflow.md

@@ -174,7 +174,7 @@ All results are written to a single HDF5 file (default: `gpec.h5`). The top-leve
 
 | Group | Contents |
 |---|---|
-| `Info/` | Run metadata: git version, mode-number ranges, ψ limit |
+| `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/` |

src/GeneralizedPerturbedEquilibrium.jl

@@ -192,6 +192,11 @@ function main_from_inputs(
     preloaded_coil_sets::Union{Nothing,Vector{ForcingTerms.CoilSet}}=nothing
 )
     total_start = time()
+    # Per-stage wall-clock seconds, written to Info/Runtimes at the end of the run.
+    runtimes = Vector{Pair{String,Float64}}()
+    # Output file this run actually wrote, last writer wins (matching where the SLAYER stage
+    # appends). `nothing` means no stage produced a file, so there is nothing to stamp.
+    written_h5 = nothing
 
     # ----------------------------------------------------------------
     # Equilibrium
@@ -212,7 +217,9 @@ function main_from_inputs(
     kf_ctrl, kinetic_profiles, kf_species = load_kinetic_context(inputs, intr, ctrl, equil)
     equil = maybe_reform_equilibrium(equil, eq_config, additional_input, intr, ctrl, kinetic_profiles)
 
-    @info "Equilibrium construction completed in $(@sprintf("%.3f", time() - equil_start)) s"
+    equil_dt = time() - equil_start
+    push!(runtimes, "equilibrium" => equil_dt)
+    @info "Equilibrium construction completed in $(@sprintf("%.3f", equil_dt)) s"
 
     # Early exit if user only requested equilibrium setup
     if equil.config.force_termination
@@ -242,7 +249,7 @@ function main_from_inputs(
 
     locstab, ballooning_boundary = run_local_stability(ctrl, equil)
     metric, ffit = prepare_force_free_states!(intr, ctrl, equil, kf_ctrl, kinetic_profiles; species=kf_species)
-    ffs_result = run_force_free_states(ctrl, equil, ffit, intr, metric)
+    ffs_result = run_force_free_states(ctrl, equil, ffit, intr, metric; runtimes=runtimes)
 
     if ctrl.write_outputs_to_HDF5
         write_outputs_to_HDF5(
@@ -253,36 +260,72 @@ function main_from_inputs(
             locstab=locstab,
             ballooning_boundary=ballooning_boundary
         )
+        written_h5 = ctrl.HDF5_filename
         @info "Results written to $(ctrl.HDF5_filename)"
     end
 
-    @info "Force-Free States completed in $(@sprintf("%.3f", time() - ffs_start)) s"
+    ffs_dt = time() - ffs_start
+    push!(runtimes, "force_free_states" => ffs_dt)
+    @info "Force-Free States completed in $(@sprintf("%.3f", ffs_dt)) s"
 
     # Early exit if user only requested force-free states (SLAYER still runs).
     if ctrl.force_termination
-        slayer_result = run_slayer_stage(ffs_result, inputs, nothing)
-        @info "\n$_BANNER\n  GPEC completed successfully in $(@sprintf("%.3f", time() - total_start)) s\n$_BANNER"
+        slayer_result = run_slayer_stage(ffs_result, inputs, nothing; runtimes=runtimes)
+        # A non-nothing result means the Tearing/ group was appended inside the guarded stage.
+        slayer_result === nothing || (written_h5 = ctrl.HDF5_filename)
+        total_dt = time() - total_start
+        push!(runtimes, "total" => total_dt)
+        if written_h5 !== nothing
+            _write_runtimes!(joinpath(intr.dir_path, written_h5), runtimes)
+        end
+        @info "\n$_BANNER\n  GPEC completed successfully in $(@sprintf("%.3f", total_dt)) s\n$_BANNER"
         return (; ffs=ffs_result, pe=nothing, slayer=slayer_result)
     end
 
-    pe_state = run_perturbed_equilibrium(ffs_result, inputs, forcing_modes_snapshot, preloaded_coil_sets)
-
-    run_kinetic_forces(inputs, ffs_result, pe_state, kf_ctrl, kinetic_profiles, kf_species)
-
-    # SLAYER runs after PE so it appends to the PE output file; it falls back to the
-    # ForceFreeStates file when PE did not run.
+    # PE (and SLAYER, which appends after it) writes to PE's own output filename when set,
+    # falling back to the ForceFreeStates file.
     pe_file = if "PerturbedEquilibrium" in keys(inputs)
         pe_out = get(inputs["PerturbedEquilibrium"], "output_filename", "")
         isempty(pe_out) ? ctrl.HDF5_filename : pe_out
     else
         ctrl.HDF5_filename
     end
-    slayer_result = run_slayer_stage(ffs_result, inputs, pe_file)
+
+    pe_start = time()
+    pe_state = run_perturbed_equilibrium(ffs_result, inputs, forcing_modes_snapshot, preloaded_coil_sets)
+    # Record only when the stage ran; an absent [PerturbedEquilibrium] section would otherwise
+    # stamp a ~0 s entry for work that never happened.
+    if "PerturbedEquilibrium" in keys(inputs)
+        push!(runtimes, "perturbed_equilibrium" => time() - pe_start)
+        # Mirrors the write gate inside `perturbed_equilibrium` (write flag defaults to true).
+        if get(inputs["PerturbedEquilibrium"], "write_outputs_to_HDF5", true)
+            written_h5 = pe_file
+        end
+    end
+
+    kf_start = time()
+    run_kinetic_forces(inputs, ffs_result, pe_state, kf_ctrl, kinetic_profiles, kf_species)
+    if "KineticForces" in keys(inputs)
+        push!(runtimes, "kinetic_forces" => time() - kf_start)
+        # Mirrors the write gate inside `run_kinetic_forces` (needs a PE state to contract against).
+        if pe_state !== nothing && kf_ctrl.write_outputs_to_HDF5
+            written_h5 = kf_ctrl.HDF5_filename
+        end
+    end
+
+    slayer_result = run_slayer_stage(ffs_result, inputs, pe_file; runtimes=runtimes)
+    # A non-nothing result means the Tearing/ group was appended inside the guarded stage.
+    slayer_result === nothing || (written_h5 = pe_file)
 
     # ----------------------------------------------------------------
     # Done
     # ----------------------------------------------------------------
-    @info "\n$_BANNER\n  GPEC completed successfully in $(@sprintf("%.3f", time() - total_start)) s\n$_BANNER"
+    total_dt = time() - total_start
+    push!(runtimes, "total" => total_dt)
+    if written_h5 !== nothing
+        _write_runtimes!(joinpath(intr.dir_path, written_h5), runtimes)
+    end
+    @info "\n$_BANNER\n  GPEC completed successfully in $(@sprintf("%.3f", total_dt)) s\n$_BANNER"
 
     # TODO: Do not allow perturbed equilibrium calculations if zero crossings are found
 
@@ -602,14 +645,16 @@ end
 
 Run the formalism selected by `ctrl.integrator` — the standalone Galerkin solve, or the
 Euler-Lagrange sweep with its free-boundary energies and Δ′ BVP — and publish its products as a
-`ForceFreeStatesResult`.
+`ForceFreeStatesResult`. A `runtimes` collector, when given, receives the Galerkin solve's
+wall-clock seconds as a `"galerkin" => dt` pair for the `Info/Runtimes` record.
 """
 function run_force_free_states(
     ctrl::ForceFreeStatesControl,
     equil::Equilibrium.PlasmaEquilibrium,
     ffit,
     intr::ForceFreeStatesInternal,
-    metric
+    metric;
+    runtimes::Union{Nothing,Vector{Pair{String,Float64}}}=nothing
 )
     nstring = _mode_range_label(intr)
 
@@ -626,7 +671,9 @@ function run_force_free_states(
         gal_start = time()
         wv = ctrl.vac_flag ? first(ForceFreeStates.compute_scaled_wv(ctrl, equil, intr)) : nothing
         gal_data, gal_dp = galerkin_solve(ctrl, equil, ffit, intr; wv=wv)
-        @info "Galerkin solve completed in $(@sprintf("%.3f", time() - gal_start)) s"
+        gal_dt = time() - gal_start
+        runtimes === nothing || push!(runtimes, "galerkin" => gal_dt)
+        @info "Galerkin solve completed in $(@sprintf("%.3f", gal_dt)) s"
     else
         # Integrate Euler-Lagrange Equation
         if ctrl.verbose
@@ -988,9 +1035,12 @@ end
 
 Run the SLAYER tearing-mode analysis off the force-free-states `result`, appending its group to
 `pe_file` (or the force-free-states output when PE did not run). Needs only the result, so it
-runs in both the `force_termination = true` path and the full pipeline.
+runs in both the `force_termination = true` path and the full pipeline. A `runtimes` collector,
+when given, receives the solver's wall-clock seconds as a `"tearing" => dt` pair for the
+`Info/Runtimes` record.
 """
-function run_slayer_stage(result::ForceFreeStatesResult, inputs::Dict{String,Any}, pe_file::Union{String,Nothing})
+function run_slayer_stage(result::ForceFreeStatesResult, inputs::Dict{String,Any}, pe_file::Union{String,Nothing};
+    runtimes::Union{Nothing,Vector{Pair{String,Float64}}}=nothing)
     ("SLAYER" in keys(inputs)) || return nothing
     # SLAYER is a post-processing diagnostic. A failure here must not
     # discard the equilibrium / stability / PE results already computed,
@@ -1003,7 +1053,9 @@ function run_slayer_stage(result::ForceFreeStatesResult, inputs::Dict{String,Any
         slayer_start = time()
         slayer_result = Runner.run_slayer(result, slayer_ctrl;
             dir_path=result.dir_path)
-        @info "SLAYER completed in $(@sprintf("%.3f", time() - slayer_start)) s"
+        slayer_dt = time() - slayer_start
+        runtimes === nothing || push!(runtimes, "tearing" => slayer_dt)
+        @info "SLAYER completed in $(@sprintf("%.3f", slayer_dt)) s"
         h5_filename = pe_file === nothing ? result.control.HDF5_filename : pe_file
         h5_path = joinpath(result.dir_path, h5_filename)
         # Append the Tearing/ group; create the file if no prior stage wrote
@@ -1379,6 +1431,29 @@ function _write_coil_snapshot!(h5_path::String, coil_sets::Vector{ForcingTerms.C
     return nothing
 end
 
+"""
+    _write_runtimes!(h5_path::String, runtimes)
+
+Write the per-stage wall-clock seconds collected during a run into `Info/Runtimes/` of an
+existing gpec.h5 file. `runtimes` iterates `stage => seconds` pairs; only the stages that ran
+are recorded. Metadata comes from `RUNTIME_H5_ANNOTATIONS`, which skips the absent stages.
+These timings are informational only — machine- and load-dependent, never a regression quantity.
+
+Call it only when this run produced the file; the `isfile` guard alone would also stamp a
+stale gpec.h5 left over from an earlier run.
+"""
+function _write_runtimes!(h5_path::String, runtimes)
+    isfile(h5_path) || return nothing
+    h5open(h5_path, "r+") do out_h5
+        haskey(out_h5, "Info/Runtimes") && HDF5.delete_object(out_h5, "Info/Runtimes")
+        for (stage, dt) in runtimes
+            out_h5["Info/Runtimes/$stage"] = dt
+        end
+        Utilities.HDF5Annotations.annotate!(out_h5, RUNTIME_H5_ANNOTATIONS)
+    end
+    return nothing
+end
+
 """
     write_imas(dd, result)
 

src/HDF5Schema.jl

@@ -246,6 +246,19 @@ const MAIN_H5_ANNOTATIONS = [
     "SurfaceGeometries/Wall/z" => (; long_name="Cartesian z of wall point cloud", units="m"),
 ]
 
+# Per-stage wall-clock records written by `_write_runtimes!` after every stage has run.
+# Stages that did not run leave their path absent and are skipped by `annotate!`.
+const RUNTIME_H5_ANNOTATIONS = [
+    "Info/Runtimes/equilibrium" => (; long_name="wall-clock time of the Equilibrium construction stage", units="s"),
+    "Info/Runtimes/galerkin" => (; long_name="wall-clock time of the Galerkin outer-region solve (measured inside force_free_states, not additional to it)", units="s"),
+    "Info/Runtimes/force_free_states" => (; long_name="wall-clock time of the ForceFreeStates stability stage", units="s"),
+    "Info/Runtimes/tearing" => (; long_name="wall-clock time of the Tearing stage (SLAYER inner-layer solve)", units="s"),
+    "Info/Runtimes/perturbed_equilibrium" => (; long_name="wall-clock time of the PerturbedEquilibrium stage", units="s"),
+    "Info/Runtimes/kinetic_forces" => (; long_name="wall-clock time of the KineticForces (NTV) stage", units="s"),
+    "Info/Runtimes/total" =>
+        (; long_name="wall-clock time of the full GPEC run (the per-stage values nest rather than partition it, so they do not sum to this)", units="s"),
+]
+
 # Euler-Lagrange operator matrices: same wording per letter, Ideal/ and Kinetic/ variants.
 const _ELM_IDEAL_LETTERS = [
     ("A", "Euler-Lagrange primitive coefficient matrix A"),

test/runtests_fullruns.jl

@@ -28,6 +28,9 @@ using HDF5
             et = read(h5["ForceFreeStates/FreeBoundaryStability/eigenmode_energies"])
             @test isfinite(real(et[1]))
             @test real(et[1]) > 0
+            # Per-stage wall-clock records (informational, not regression quantities).
+            @test haskey(h5, "Info/Runtimes/total") && read(h5["Info/Runtimes/total"]) > 0
+            @test haskey(h5, "Info/Runtimes/force_free_states")
         end
         rm(joinpath(ex3, "gpec.h5"); force=true)
         true

test/runtests_rerun_from_h5.jl

@@ -7,10 +7,11 @@ using TOML
 
 # Collect every leaf dataset path under an open HDF5 file, skipping the groups/paths that
 # legitimately differ between a source run and its replay (`Input/` is re-emitted with the
-# rerun's own filename/TOML blob; `Info/git_version` reflects the running commit).
+# rerun's own filename/TOML blob; `Info/git_version` reflects the running commit;
+# `Info/Runtimes/` records wall-clock seconds, which never repeat).
 function _rerun_leaf_paths(h5)
     skip_toplevel = Set(["Input"])
-    skip_paths = Set(["Info/git_version"])
+    skip_paths = Set(["Info/git_version", "Info/Runtimes"])
     paths = String[]
     function walk(node, prefix)
         for k in keys(node)