Unification of the three integrators - PR #1

.pre-commit-config.yaml

@@ -69,5 +69,5 @@ repos:
       - id: toml-no-deprecated-keys
         name: 'TOML conventions: no deprecated config keys'
         language: pygrep
-        entry: '^(mer_flag|force_wv_symmetry|ode_flag|cyl_flag|mat_flag|power_bp|power_b|power_r|power_rc)\s*='
+        entry: '^(mer_flag|force_wv_symmetry|ode_flag|cyl_flag|mat_flag|use_riccati|use_parallel|parallel_threads|populate_dense_xi|power_bp|power_b|power_r|power_rc)\s*='
         files: ^(examples/.*\.toml|test/test_data/.*\.toml)$

REFACTOR_PLAN.md

@@ -8,11 +8,17 @@
 > changes ride with it) so everyone with open PRs can see what is coming and where it
 > will touch their work. Key coordination points:
 >
-> - The PR sequence below assumes **nothing else merges into `develop` before it
->   finishes**. If your PR must land mid-sequence, talk to Matthew first.
-> - The HDF5 schema PRs (#363/#364) are treated as guidance on the final HDF5 shape;
->   writers refactored here keep today's dataset paths and will be re-targeted by
->   those PRs afterward. #367 (immutable control structs) merges after this sequence.
+> - The PR sequence below assumes **nothing else merges into `develop` mid-sequence**.
+>   If your PR must land before it finishes, talk to Matthew first.
+> - **Amendment (post-PR-1):** the module-mirroring HDF5 schema (#363) and the vacuum
+>   surface-inductance migration (#345) merged into develop before PR 1 branched, so
+>   the sequence is built on the NEW schema (`ForceFreeStates/…`, `SingularSurfaces/…`,
+>   `LocalStability/…`, `Equilibrium/…`). Writer refactors in PR 3/4 therefore target
+>   that schema directly; dataset-path references below have been updated accordingly.
+>   Only #364 (self-describing metadata) remains as a later re-target, and #367
+>   (immutable control structs) still merges after this sequence. Neither merged PR
+>   changes any decision: #345 left `calc_surface_inductance` in PerturbedEquilibrium
+>   (only its Vacuum support types moved), so the PE consumer map stands.
 > - The regression harness will be re-baselined during this work; a fresh
 >   Fortran-agreement comparison is the final validation gate for the whole sequence.
 > - The "Decision record" and "Verified code facts" sections are settled — please do
@@ -53,8 +59,8 @@ pe  = perturbed_equilibrium(ffs, rmp)
 | D7 | Public API via **CommonSolve.jl**: `solve(eq::PlasmaEquilibrium, alg; kwargs...)`. `PlasmaEquilibrium(path; kwargs...)` constructor. Module names unchanged. |
 | D8 | Galerkin standalone computes its own vacuum `wv` (no ODE state needed — verified); its result has `free_boundary = nothing`. |
 | D9 | TOML: new `integrator = "forward"|"riccati"|"galerkin"` key. Old keys (`use_riccati`, `use_parallel`, `parallel_threads`, `populate_dense_xi`, later `gal_flag`) go to the `_DEPRECATED_FFS_KEYS` warn-and-ignore list AND the `toml-no-deprecated-keys` pre-commit hook pattern. |
-| D10 | No back-compat burden; examples/fixtures updated freely; regression re-baselining accepted (HDF5 schema PRs #363/#364 churn it anyway). Final validation = fresh Fortran comparison after the sequence. Assume nothing else merges first. |
-| D11 | New structs immutable from day one (eases the later #367 merge). HDF5 writers become functions on result structs (schema itself unchanged here; #363/#364 re-target them later). |
+| D10 | No back-compat burden; examples/fixtures updated freely; regression re-baselining accepted. Final validation = fresh Fortran comparison after the sequence. Nothing else merges mid-sequence without coordination (#363/#345 landed before PR 1 and are absorbed — see header amendment). |
+| D11 | New structs immutable from day one (eases the later #367 merge). HDF5 writers become functions on result structs, keeping the merged #363 schema paths unchanged; #364 (metadata) re-targets them later. |
 | D12 | Analysis module reads HDF5 files, not live structs — untouched except where dataset names would change (they don't in this plan). |
 
 ### Verified code facts the workers must not re-derive
@@ -227,7 +233,7 @@ All code must be JuliaFormatter-clean per `.JuliaFormatter.toml` before commit.
     `integrator="riccati"`) so the canonical Δ′-matrix fixture survives the
     DIIID-like_ideal switch to forward. Add a matching regression case
     `regression-harness/cases/diiid_n1_riccati.toml` tracking
-    `singular/delta_prime_matrix`-derived quantities (mirror the Δ′ entries of the
+    `SingularSurfaces/delta_prime_matrix`-derived quantities (mirror the Δ′ entries of the
     existing `diiid_n1` case; ξ/PE quantities stay on `diiid_n1`).
 - `benchmarks/benchmark_threads.jl`, `benchmarks/benchmark_delta_prime_methods.jl`:
   update flag names (`use_riccati`/`parallel_threads` → `integrator`/`nchunks`);
@@ -323,7 +329,7 @@ All code must be JuliaFormatter-clean per `.JuliaFormatter.toml` before commit.
 
 Full test suite; targeted `runtests_resist_eval.jl`, `runtests_fullruns.jl`
 (exercises `local_stability_flag=true` decks); docs build (missing-docs gate);
-harness `--cases diiid_n1 --refs develop,local` (locstab datasets must be identical).
+harness `--cases diiid_n1 --refs develop,local` (`LocalStability/*` datasets must be identical; note the #363 group name already matches the new module name).
 
 ---
 
@@ -370,8 +376,8 @@ end
 Contract of `solution`/`solution_basis`:
 - Forward → dense EL-basis odet, `:el_axis`. PE-usable.
 - Riccati → its odet IS carried (psi/q/crit/edge_scan are valid, u_store is chunk
-  snapshots), `:riccati`. NOT PE-usable; HDF5 `integration/psi|crit|EdgeScan` still
-  written from it, `integration/xi_*` written empty.
+  snapshots), `:riccati`. NOT PE-usable; HDF5 `ForceFreeStates/Solutions/ForwardIntegration/psi|crit` and `EdgeScan/*` still
+  written from it, `.../ForwardIntegration/xi_*` written empty.
 - Galerkin with `gal_match_flag` → `gal_matched_odestate(...)`, `:gal_native`
   (PE-usable; `du_store_populated=true` analytic derivatives). Without match →
   `nothing`, `:none`.
@@ -454,7 +460,7 @@ crosses module boundaries after `build_result`.
 ### 5.4 Verification
 
 Full suite; `runtests_fullruns.jl` (all decks — forward decks produce identical
-HDF5 vs pre-PR, riccati decks now emit empty `integration/xi_*` + PE-skip warnings);
+HDF5 vs pre-PR, riccati decks now emit empty `ForwardIntegration/xi_*` + PE-skip warnings);
 harness `--cases diiid_n1,solovev_n1 --refs develop,local` (forward-deck tracked
 quantities unchanged); docs build.
 
@@ -512,14 +518,14 @@ logic stays a pre-FFS stage but is owned by the FFS-facing function
   disappear from the HDF5 — accepted per D10; gal datasets identical because
   `galerkin_solve` inputs are unchanged). `gal_*` sub-knobs stay (they become
   `Galerkin(...)` fields in PR 5).
-- FFS writer: tolerate `solution === nothing` (write empty `integration/*` datasets —
+- FFS writer: tolerate `solution === nothing` (write empty `ForwardIntegration/*` datasets —
   extend the existing empty-fallback pattern).
 
 ### 6.3 Tests / verification
 
 - Update `runtests_fullruns.jl` gal decks' expectations (gal-only HDF5).
 - New testset (in `runtests_fullruns.jl` or the gal tests): `integrator="galerkin"`
-  on `LAR_ideal_match_test` produces `galerkin/delta` identical to the PR-3 additive
+  on `LAR_ideal_match_test` produces the Galerkin `delta` dataset identical to the PR-3 additive
   run (same `wv` by construction — assert against a stored reference or a paired
   riccati+gal_flag run on the pre-PR commit during development).
 - Full suite; harness (gal cases if present, plus diiid/solovev); docs.
@@ -660,8 +666,15 @@ manual smoke: run the 4-line UX from the Context section in a REPL against
 
 ## 10. Progress
 
-- [ ] PR 1 — `refactor/riccati-unification`
-- [ ] PR 2 — `refactor/local-stability-module`
+- [ ] PR 1 — `refactor/riccati-unification` — **implemented, in review.** Two deltas
+  from the §3 spec, both improvements: the new Δ′ example references the DIIID geqdsk
+  by relative path instead of copying it, and the TOML sweep covered six regression
+  fixtures (two more had landed on develop since the plan was written), all `forward`.
+- [ ] PR 2 — `refactor/local-stability-module` — **implemented, in review.** One delta
+  from the §4 spec: the signature change also required updating two call-site groups the
+  section did not list — `examples/DIIID-like_ideal_example/analyze_example.jl` (five
+  ballooning entry points) and two docstring cross-references in
+  `src/Analysis/ForceFreeStates.jl`.
 - [ ] PR 3 — `refactor/forcefreestates-result`
 - [ ] PR 4 — `refactor/staged-main`
 - [ ] PR 5 — `feature/solve-api`

benchmarks/benchmark_delta_prime_methods.jl

@@ -23,7 +23,7 @@ function setup_and_run_solovev()
     ex = joinpath(@__DIR__, "..", "test", "test_data", "regression_solovev_ideal_example")
     inputs = TOML.parsefile(joinpath(ex, "gpec.toml"))
     inputs["ForceFreeStates"]["verbose"] = false
-    inputs["ForceFreeStates"]["use_riccati"] = true
+    inputs["ForceFreeStates"]["integrator"] = "riccati"
     intr = FFS.ForceFreeStatesInternal(; dir_path=ex)
     ctrl = FFS.ForceFreeStatesControl(;
         (Symbol(k) => v for (k, v) in inputs["ForceFreeStates"])...)
@@ -42,7 +42,7 @@ function setup_and_run_solovev()
     intr.numpert_total = intr.mpert * intr.npert
     metric = FFS.make_metric(equil, intr.mpert)
     ffit = FFS.make_matrix(equil, intr, metric)
-    odet = FFS.riccati_eulerlagrange_integration(ctrl, equil, ffit, intr)
+    odet, _, _, _ = FFS.riccati_eulerlagrange_integration(ctrl, equil, ffit, intr)
     return ctrl, equil, ffit, intr, odet
 end
 

benchmarks/benchmark_threads.jl

@@ -1,17 +1,17 @@
-# Thread-scaling benchmark for the bidirectional parallel FM integration.
-# Runs the Solovev (N=8) and DIIID-like (N=26) examples with use_parallel=true
-# across 1, 2, 4, 8 threads and compares against the serial Riccati path.
+# Thread-scaling benchmark for the Riccati chunked propagator integration.
+# Runs the Solovev (N=8) and DIIID-like (N=26) examples with integrator="riccati"
+# across 1, 2, 4, 8 threads and compares against the forward path. The chunk list is
+# thread-independent, so et[1] must not move with the thread count — only wall-clock.
 #
 # Usage (from JPEC_main root):
 #   for t in 1 2 4 8; do julia -t $t --project=. benchmarks/benchmark_threads.jl; done
 
 using GeneralizedPerturbedEquilibrium, TOML, Printf, Statistics
 
-function run_ffs(ex; use_parallel, use_riccati=false)
+function run_ffs(ex; integrator)
     inputs = TOML.parsefile(joinpath(ex, "gpec.toml"))
     inputs["ForceFreeStates"]["verbose"] = false
-    inputs["ForceFreeStates"]["use_parallel"] = use_parallel
-    inputs["ForceFreeStates"]["use_riccati"] = use_riccati
+    inputs["ForceFreeStates"]["integrator"] = integrator
     inputs["ForceFreeStates"]["write_outputs_to_HDF5"] = false
     intr = GeneralizedPerturbedEquilibrium.ForceFreeStates.ForceFreeStatesInternal(; dir_path=ex)
     ctrl = GeneralizedPerturbedEquilibrium.ForceFreeStates.ForceFreeStatesControl(;
@@ -36,17 +36,17 @@ function run_ffs(ex; use_parallel, use_riccati=false)
     return real(vac.et[1]), intr.numpert_total
 end
 
-function timed_run(ex; use_parallel, use_riccati=false, nwarm=1, nrep=2)
+function timed_run(ex; integrator, nwarm=1, nrep=2)
     # Warmup
     for _ in 1:nwarm
-        run_ffs(ex; use_parallel, use_riccati)
+        run_ffs(ex; integrator)
     end
     # Timed runs
     times = Float64[]
     local et1, N
     for _ in 1:nrep
         t0 = time()
-        et1, N = run_ffs(ex; use_parallel, use_riccati)
+        et1, N = run_ffs(ex; integrator)
         push!(times, time() - t0)
     end
     return mean(times), et1, N
@@ -60,18 +60,14 @@ diiid_ex = joinpath(root, "examples", "DIIID-like_ideal_example")
 println("\n=== Thread-scaling benchmark ($(nthreads) thread(s)) ===\n")
 
 for (label, ex) in [("Solovev", sol_ex), ("DIIID-like", diiid_ex)]
-    t_std, et_std, N = timed_run(ex; use_parallel=false, use_riccati=false)
-    t_ric, et_ric, _ = timed_run(ex; use_parallel=false, use_riccati=true)
-    t_par, et_par, _ = timed_run(ex; use_parallel=true, use_riccati=false)
+    t_fwd, et_fwd, N = timed_run(ex; integrator="forward")
+    t_ric, et_ric, _ = timed_run(ex; integrator="riccati")
 
-    err_ric = abs(et_ric - et_std) / abs(et_std) * 100
-    err_par = abs(et_par - et_std) / abs(et_std) * 100
+    err_ric = abs(et_ric - et_fwd) / abs(et_fwd) * 100
 
     println("$label (N=$N, nthreads=$nthreads)")
-    @printf("  standard   et[1]=%.5f  t=%.2fs  speedup=1.00×\n", et_std, t_std)
+    @printf("  forward    et[1]=%.5f  t=%.2fs  speedup=1.00×\n", et_fwd, t_fwd)
     @printf("  riccati    et[1]=%.5f  t=%.2fs  speedup=%.2f×  err=%.4f%%\n",
-        et_ric, t_ric, t_std/t_ric, err_ric)
-    @printf("  parallel   et[1]=%.5f  t=%.2fs  speedup=%.2f×  err=%.4f%%\n",
-        et_par, t_par, t_std/t_par, err_par)
+        et_ric, t_ric, t_fwd/t_ric, err_ric)
     println()
 end

benchmarks/compare_gal_vs_el.jl

@@ -5,7 +5,7 @@
 # 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).
 #
-# Needs ONE gpec.h5 from a run with populate_dense_xi=true (EL dense u_store), gal_match_flag=true,
+# Needs ONE gpec.h5 from a run with integrator="forward" (EL dense u_store), gal_match_flag=true,
 # gal_ideal_flag=true (gal ideal matched), vac_flag=true (energy operator).
 # Usage: julia --project=. benchmarks/compare_gal_vs_el.jl [gpec.h5] [out.png] [mode]
 #   mode = "highest" (default, most stable), "lowest" (most unstable), or an integer eigenmode index.

benchmarks/compare_jbgradpsi_m2.jl

@@ -1,16 +1,16 @@
 # 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)
+#   (2) FORWARD ξ                  (gal_match_flag=false)
 #
 # PE writes no ψ grid, so it's reconstructed: gal-ideal → ForceFreeStates/Solutions/GalerkinIntegration/Solution/psi minus issing points;
-# shooting → ForceFreeStates/Solutions/ForwardIntegration/psi.
+# forward → ForceFreeStates/Solutions/ForwardIntegration/psi.
 # Usage: julia --project=. benchmarks/compare_jbgradpsi_m2.jl [gal_h5] [shoot_h5] [out.png] [m]
 
 using HDF5, Plots, Printf
 
 gal_h5 = length(ARGS) >= 1 ? ARGS[1] : "/tmp/gal_ideal_test/gpec.h5"
-sh_h5 = length(ARGS) >= 2 ? ARGS[2] : "/tmp/shooting_test/gpec.h5"
+sh_h5 = length(ARGS) >= 2 ? ARGS[2] : "/tmp/forward_test/gpec.h5"
 outpng = length(ARGS) >= 3 ? ARGS[3] : joinpath(@__DIR__, "cmp_jbgradpsi_m2.png")
 mtarget = length(ARGS) >= 4 ? parse(Int, ARGS[4]) : 2
 
@@ -23,7 +23,7 @@ pa_g, psi_g, mlow, sing_psi, sing_m = h5open(gal_h5) do f
     (pa, read(f["ForceFreeStates/Solutions/GalerkinIntegration/Solution/psi"])[.!iss], read(f["Info/mlow"]),
         read(f["SingularSurfaces/GalerkinDeltaPrime/rational_psi"]), read(f["SingularSurfaces/GalerkinDeltaPrime/rational_m"]))
 end
-# shooting run: PE grid = ForceFreeStates/Solutions/ForwardIntegration/psi
+# forward run: PE grid = ForceFreeStates/Solutions/ForwardIntegration/psi
 pa_s, psi_s = h5open(sh_h5) do f
     (to_c(read(f["PerturbedEquilibrium/Response/psi_area"])), read(f["ForceFreeStates/Solutions/ForwardIntegration/psi"]))
 end
@@ -43,15 +43,15 @@ if !isnan(psi_res)
     ig = argmin(abs.(psi_g .- psi_res))
     is = argmin(abs.(psi_s .- psi_res))
     @printf("m=%d resonant surface ψ=%.4f\n", mtarget, psi_res)
-    @printf("  nearest |b^ψ_area|:  gal-ideal=%.4e (ψ=%.4f)   shooting=%.4e (ψ=%.4f)   ratio=%.3f\n",
+    @printf("  nearest |b^ψ_area|:  gal-ideal=%.4e (ψ=%.4f)   forward=%.4e (ψ=%.4f)   ratio=%.3f\n",
         abs(g[ig]), psi_g[ig], abs(s[is]), psi_s[is], abs(g[ig]) / abs(s[is]))
 end
 
 plt = plot(; size=(1000, 620), xlabel="ψ_N", ylabel="|b^ψ / ⟨J·|∇ψ|⟩|  (area-normalized)",
-    title="m=$mtarget perturbed normal field — gal-ideal vs shooting PE", legend=:topleft,
+    title="m=$mtarget perturbed normal field — gal-ideal vs forward PE", legend=:topleft,
     left_margin=13Plots.mm, bottom_margin=5Plots.mm, right_margin=4Plots.mm)
 plot!(plt, psi_g, abs.(g); lw=2.2, color=1, label="gal-ideal ξ → PE")
-plot!(plt, psi_s, abs.(s); lw=1.8, color=2, ls=:dash, label="shooting ξ → PE")
+plot!(plt, psi_s, abs.(s); lw=1.8, color=2, ls=:dash, label="forward ξ → PE")
 for (j, ps) in enumerate(sing_psi)
     ism2 = sing_m[j] == mtarget
     vline!(plt, [ps]; color=(ism2 ? :red : :black), ls=:dot, lw=(ism2 ? 1.8 : 1),

benchmarks/plot_xi_eigenmode.jl

@@ -5,7 +5,7 @@
 # largest Re(eigenvalue) and reconstructs its radial profile by projecting the EL fundamental matrix
 # (ForceFreeStates/Solutions/ForwardIntegration/xi_psi) onto that edge eigenvector:  c = U_edge \ w,  ξ(ψ) = U(ψ)·c.
 #
-# Requires a run with populate_dense_xi=true so ForceFreeStates/Solutions/ForwardIntegration/xi_psi is the dense axis-basis fundamental
+# Requires a run with integrator="forward" 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]
 

benchmarks/scan_resistivity_m2.jl

@@ -1,6 +1,6 @@
 # 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
+# Overlays the forward (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).
 # Usage: julia --project=. benchmarks/scan_resistivity_m2.jl [out.png] [m]
 
@@ -46,10 +46,10 @@ for (i, (d, e)) in enumerate(zip(scandirs, etas))
     psi, v = read_m2(joinpath(d, "gpec.h5"); gal=true)
     plot!(plt, psi, abs.(v); color=cols[i], lw=2, label=@sprintf("η = %g×η_ref", round(e / eta_ref; sigdigits=2)))
 end
-# shooting (ideal, η→0) reference
-if isfile("/tmp/shooting_test/gpec.h5")
-    psis, vs = read_m2("/tmp/shooting_test/gpec.h5"; gal=false)
-    plot!(plt, psis, abs.(vs); color=:black, ls=:dash, lw=2, label="shooting (ideal, η→0)")
+# forward (ideal, η→0) reference
+if isfile("/tmp/forward_test/gpec.h5")
+    psis, vs = read_m2("/tmp/forward_test/gpec.h5"; gal=false)
+    plot!(plt, psis, abs.(vs); color=:black, ls=:dash, lw=2, label="forward (ideal, η→0)")
 end
 isnan(psi_res) || vline!(plt, [psi_res]; color=:red, ls=:dot, lw=1.6, label="q=$mtarget surface")
 

benchmarks/scan_rotation_m2.jl

@@ -1,7 +1,7 @@
 # Plot the m=2 area-normalized b^ψ (PerturbedEquilibrium/Response/psi_area) across a ROTATION scan
 # of the resistive gal matched PE runs (gal_match_flag=true, gal_ideal_flag=false), fixed η=8e-8,
 # rotation f = 1,2,4,8,16 Hz (forced eigenvalue γ_s = 2πi·n·f). One curve per rotation; overlays the
-# shooting (ideal) reference. Scan dirs produced by the bash loop over /tmp/rotscan_<f>.
+# forward (ideal) reference. Scan dirs produced by the bash loop over /tmp/rotscan_<f>.
 # Usage: julia --project=. benchmarks/scan_rotation_m2.jl [out.png] [m]
 
 using HDF5, Plots, Printf, TOML
@@ -41,9 +41,9 @@ for (i, (d, r)) in enumerate(zip(scandirs, rots))
     psi, v = read_m2(joinpath(d, "gpec.h5"); gal=true)
     plot!(plt, psi, abs.(v); color=cols[i], lw=2, label=@sprintf("f = %g Hz", r))
 end
-if isfile("/tmp/shooting_test/gpec.h5")
-    psis, vs = read_m2("/tmp/shooting_test/gpec.h5"; gal=false)
-    plot!(plt, psis, abs.(vs); color=:black, ls=:dash, lw=2, label="shooting (ideal)")
+if isfile("/tmp/forward_test/gpec.h5")
+    psis, vs = read_m2("/tmp/forward_test/gpec.h5"; gal=false)
+    plot!(plt, psis, abs.(vs); color=:black, ls=:dash, lw=2, label="forward (ideal)")
 end
 isnan(psi_res) || vline!(plt, [psi_res]; color=:red, ls=:dot, lw=1.6, label="q=$mtarget surface")