OpenFUSIONToolkit · matt-pharr · Aug 17, 2026 · Aug 15, 2026 · Aug 15, 2026 · Aug 15, 2026
diff --git a/.pre-commit-config.yaml b/.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|reform_eq_with_psilim|use_riccati|use_parallel|parallel_threads|populate_dense_xi|power_bp|power_b|power_r|power_rc)\s*='
+        entry: '^(mer_flag|force_wv_symmetry|ode_flag|cyl_flag|mat_flag|reform_eq_with_psilim|use_riccati|use_parallel|parallel_threads|populate_dense_xi|gal_flag|power_bp|power_b|power_r|power_rc)\s*='
         files: ^(examples/.*\.toml|test/test_data/.*\.toml)$
diff --git a/Project.toml b/Project.toml
@@ -6,6 +6,7 @@ version = "0.1.0"
 
 [deps]
 AdaptiveArrayPools = "4f381ef7-9af0-4cbe-99d4-cf36d7b0f233"
+CommonSolve = "38540f10-b2f7-11e9-35d8-d573e4eb0ff2"
 Contour = "d38c429a-6771-53c6-b99e-75d170b6e991"
 Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
 DelaunayTriangulation = "927a84f5-c5f4-47a5-9785-b46e178433df"
@@ -38,6 +39,7 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [compat]
 AdaptiveArrayPools = "0.3.5"
+CommonSolve = "0.2"
 Contour = "0.6.3"
 DelaunayTriangulation = "1.6.6"
 DelimitedFiles = "1.9.1"

diff --git a/REFACTOR_PLAN.md b/REFACTOR_PLAN.md
diff --git a/benchmarks/benchmark_diiid_ideal_ntv_torque.jl b/benchmarks/benchmark_diiid_ideal_ntv_torque.jl
@@ -31,7 +31,6 @@ using Plots
 # Load GPEC
 using GeneralizedPerturbedEquilibrium
 const GPE = GeneralizedPerturbedEquilibrium
-const FFS = GPE.ForceFreeStates
 const KF = GPE.KineticForces
 const Eq = GPE.Equilibrium
 const PE = GPE.PerturbedEquilibrium
@@ -237,15 +236,12 @@ function run_benchmark(fortran_dir::String=default_fortran_dir())
     _p("\n--- Equilibrium + ForceFreeStates (via main()) ---")
     t0 = time()
     result = GPE.main([tomldir])
-    equil = result.equil
-    intr = result.intr
-    ctrl = result.ctrl
+    ffs = result.ffs
+    equil = ffs.equil
+    metric = ffs.metric
 
     _pf("  FFS completed in %.1f s\n", time() - t0)
-    _pf("  mpert=%d, mlow=%d, mhigh=%d\n", intr.mpert, intr.mlow, intr.mhigh)
-
-    # Build metric (needed for JBB deweighting)
-    metric = FFS.make_metric(equil, intr.mpert)
+    _pf("  mpert=%d, mlow=%d, mhigh=%d\n", ffs.mpert, ffs.mlow, ffs.mhigh)
 
     # Load Fortran xclebsch data
     _p("\n--- Load Fortran xclebsch ---")
@@ -257,8 +253,8 @@ function run_benchmark(fortran_dir::String=default_fortran_dir())
         npsi_f, mpert_f, mlow_f)
     _pf("  ψ range: [%.6f, %.6f]\n", psi_grid_f[1], psi_grid_f[end])
 
-    if mpert_f != intr.mpert || mlow_f != intr.mlow
-        @warn "Mode ranges differ: Fortran mpert=$mpert_f,mlow=$mlow_f vs Julia mpert=$(intr.mpert),mlow=$(intr.mlow)"
+    if mpert_f != ffs.mpert || mlow_f != ffs.mlow
+        @warn "Mode ranges differ: Fortran mpert=$mpert_f,mlow=$mlow_f vs Julia mpert=$(ffs.mpert),mlow=$(ffs.mlow)"
     end
 
     # Build PE state from Fortran data and run JBB deweighting
@@ -291,7 +287,7 @@ function run_benchmark(fortran_dir::String=default_fortran_dir())
     kf_intr = KF.KineticForcesInternal(equil; verbose=false)
 
     # Run set_perturbation_data! — builds dbob_m, divx_m, xs_m via JBB deweighting
-    KF.set_perturbation_data!(kf_intr, pe_state, intr, equil, metric)
+    KF.set_perturbation_data!(kf_intr, pe_state, ffs, equil, metric)
 
     _pf("  JBB deweighting completed in %.1f s\n", time() - t1)
 

diff --git a/benchmarks/compare_gal_vs_el.jl b/benchmarks/compare_gal_vs_el.jl
@@ -18,13 +18,13 @@ 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
+et, wt, u1, psiE, gxi, psiG, mlow, sing_psi = h5open(h5path) do f
     (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/is_rational"])), read(f["Info/mlow"]),
-        read(f["SingularSurfaces/GalerkinDeltaPrime/rational_psi"]))
+        to_c(read(f["ForceFreeStates/Solutions/GalerkinIntegration/xi_psi"])), read(f["ForceFreeStates/Solutions/GalerkinIntegration/psi"]),
+        read(f["Info/mlow"]),
+        read(f["ForceFreeStates/Solutions/GalerkinIntegration/rational_psi"]))
 end
 
 mpert = size(u1, 1)
@@ -38,10 +38,10 @@ w = wt[:, k]
 cEL = u1[:, :, end] \ w
 xiE = reduce(hcat, (u1[:, :, ip] * cEL for ip in 1:size(u1, 3)))   # (mpert, nE)
 
-# gal ideal profile (identity-at-edge ⇒ coefficient = w); drop on-surface points
-keep = .!issing
-psiGk = psiG[keep]
-xiG = reduce(hcat, (gxi[:, ip, :] * w for ip in findall(keep)))    # (mpert, nGk)
+# gal ideal profile (identity-at-edge ⇒ coefficient = w); the closed-profile grid
+# already excludes on-surface points
+psiGk = psiG
+xiG = reduce(hcat, (gxi[:, :, ip] * w for ip in eachindex(psiG)))    # (mpert, nGk)
 
 ms = mlow .+ (0:mpert-1)
 peak = [maximum(abs, @view xiE[i, :]) for i in 1:mpert]

diff --git a/benchmarks/compare_jbgradpsi_m2.jl b/benchmarks/compare_jbgradpsi_m2.jl
@@ -3,7 +3,7 @@
 #   (1) IDEAL galerkin matched ξ   (gal_match_flag=true, gal_ideal_flag=true)
 #   (2) FORWARD ξ                  (gal_match_flag=false)
 #
-# PE writes no ψ grid, so it's reconstructed: gal-ideal → ForceFreeStates/Solutions/GalerkinIntegration/Solution/psi minus issing points;
+# PE writes no ψ grid, so it's reconstructed: gal-ideal → ForceFreeStates/Solutions/GalerkinIntegration/psi (already excludes on-surface points);
 # forward → ForceFreeStates/Solutions/ForwardIntegration/psi.
 # Usage: julia --project=. benchmarks/compare_jbgradpsi_m2.jl [gal_h5] [shoot_h5] [out.png] [m]
 
@@ -19,9 +19,8 @@ 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["PerturbedEquilibrium/Response/psi_area"]))   # [npsi, mpert]
-    iss = Bool.(read(f["ForceFreeStates/Solutions/GalerkinIntegration/Solution/is_rational"]))
-    (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"]))
+    (pa, read(f["ForceFreeStates/Solutions/GalerkinIntegration/psi"]), read(f["Info/mlow"]),
+        read(f["ForceFreeStates/Solutions/GalerkinIntegration/rational_psi"]), read(f["ForceFreeStates/Solutions/GalerkinIntegration/rational_m"]))
 end
 # forward run: PE grid = ForceFreeStates/Solutions/ForwardIntegration/psi
 pa_s, psi_s = h5open(sh_h5) do f

diff --git a/benchmarks/plot_xi_eigenmode.jl b/benchmarks/plot_xi_eigenmode.jl
@@ -22,7 +22,7 @@ et, wt, u1, psi, mlow, sing_psi = h5open(h5path) do f
         to_c(read(f["ForceFreeStates/Solutions/ForwardIntegration/xi_psi"])),
         read(f["ForceFreeStates/Solutions/ForwardIntegration/psi"]),
         read(f["Info/mlow"]),
-        haskey(f, "SingularSurfaces/GalerkinDeltaPrime/rational_psi") ? read(f["SingularSurfaces/GalerkinDeltaPrime/rational_psi"]) : Float64[])
+        haskey(f, "ForceFreeStates/Solutions/GalerkinIntegration/rational_psi") ? read(f["ForceFreeStates/Solutions/GalerkinIntegration/rational_psi"]) : Float64[])
 end
 
 mpert, _, nstep = size(u1)

diff --git a/benchmarks/scan_resistivity_m2.jl b/benchmarks/scan_resistivity_m2.jl
@@ -16,7 +16,7 @@ function read_m2(h5; gal::Bool)
     h5open(h5) do f
         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/is_rational"]))] :
+        psi = gal ? read(f["ForceFreeStates/Solutions/GalerkinIntegration/psi"]) :
               read(f["ForceFreeStates/Solutions/ForwardIntegration/psi"])
         (psi, pa[:, col])
     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["SingularSurfaces/GalerkinDeltaPrime/rational_psi"]), read(f["SingularSurfaces/GalerkinDeltaPrime/rational_m"]))
+    (read(f["ForceFreeStates/Solutions/GalerkinIntegration/rational_psi"]), read(f["ForceFreeStates/Solutions/GalerkinIntegration/rational_m"]))
 end
 psi_res = mtarget in sing_m ? sing_psi[findfirst(==(mtarget), sing_m)] : NaN
 

diff --git a/benchmarks/scan_rotation_m2.jl b/benchmarks/scan_rotation_m2.jl
@@ -15,7 +15,7 @@ function read_m2(h5; gal::Bool)
     h5open(h5) do f
         pa = to_c(read(f["PerturbedEquilibrium/Response/psi_area"]))
         col = mtarget - read(f["Info/mlow"]) + 1
-        psi = gal ? read(f["ForceFreeStates/Solutions/GalerkinIntegration/Solution/psi"])[.!Bool.(read(f["ForceFreeStates/Solutions/GalerkinIntegration/Solution/is_rational"]))] :
+        psi = gal ? read(f["ForceFreeStates/Solutions/GalerkinIntegration/psi"]) :
               read(f["ForceFreeStates/Solutions/ForwardIntegration/psi"])
         (psi, pa[:, col])
     end
@@ -29,7 +29,7 @@ scandirs, rots = scandirs[ord], rots[ord]
 @printf("%d scan runs: rotation f = %s Hz  (η fixed = 8e-8)\n", length(rots), join((@sprintf("%g", r) for r in rots), ", "))
 
 sing_psi, sing_m = h5open(joinpath(scandirs[1], "gpec.h5")) do f
-    (read(f["SingularSurfaces/GalerkinDeltaPrime/rational_psi"]), read(f["SingularSurfaces/GalerkinDeltaPrime/rational_m"]))
+    (read(f["ForceFreeStates/Solutions/GalerkinIntegration/rational_psi"]), read(f["ForceFreeStates/Solutions/GalerkinIntegration/rational_m"]))
 end
 psi_res = mtarget in sing_m ? sing_psi[findfirst(==(mtarget), sing_m)] : NaN
 

diff --git a/benchmarks/verify_gal_ideal.jl b/benchmarks/verify_gal_ideal.jl
@@ -2,27 +2,30 @@
 #   1. cout / cin / deltar all zero (no resistive plasma combination, no inner layer)
 #   2. matched ξ_j == the gal coil column sols(:,:,2·msing+j)   (and ξ′ likewise)
 # Usage: julia --project=. benchmarks/verify_gal_ideal.jl [path/to/gpec.h5]
+# Requires [DEBUG] gal_basis_output = true in the deck (compares against the raw-basis dump).
 using HDF5, Printf, LinearAlgebra
 
 h5path = length(ARGS) >= 1 ? ARGS[1] : "/tmp/gal_ideal_test/gpec.h5"
 to_c(a) = eltype(a) <: Complex ? ComplexF64.(a) : map(x -> ComplexF64(x.re, x.im), a)
 
-cout, deltar, mxi, mdxi, sols, sols_d, sing_psi = h5open(h5path) do f
+cout, deltar, mxi, mdxi, sols, sols_d, iss, sing_psi = h5open(h5path) do f
     (to_c(read(f["ForceFreeStates/Solutions/GalerkinIntegration/Match/cout"])), to_c(read(f["ForceFreeStates/Solutions/GalerkinIntegration/Match/Delta_r"])),
-        to_c(read(f["ForceFreeStates/Solutions/GalerkinIntegration/Match/xi"])), to_c(read(f["ForceFreeStates/Solutions/GalerkinIntegration/Match/dxidpsi"])),
-        to_c(read(f["ForceFreeStates/Solutions/GalerkinIntegration/Solution/xi_psi"])), to_c(read(f["ForceFreeStates/Solutions/GalerkinIntegration/Solution/dxi_psidpsi"])),
-        read(f["SingularSurfaces/GalerkinDeltaPrime/rational_psi"]))
+        to_c(read(f["ForceFreeStates/Solutions/GalerkinIntegration/xi_psi"])), to_c(read(f["ForceFreeStates/Solutions/GalerkinIntegration/dxi_psidpsi"])),
+        to_c(read(f["ForceFreeStates/Solutions/GalerkinIntegration/Basis/xi_psi"])), to_c(read(f["ForceFreeStates/Solutions/GalerkinIntegration/Basis/dxi_psidpsi"])),
+        Bool.(read(f["ForceFreeStates/Solutions/GalerkinIntegration/Basis/is_rational"])),
+        read(f["ForceFreeStates/Solutions/GalerkinIntegration/rational_psi"]))
 end
 msing = length(sing_psi)
-mpert, ngrid, mcoil = size(mxi)
+mpert, mcoil, ngrid = size(mxi)
+keep = .!iss
 
 @printf("[1] ‖cout‖ = %.2e, ‖deltar‖ = %.2e   %s\n", norm(cout), norm(deltar),
     (norm(cout) == 0 && norm(deltar) == 0) ? "✓ no resistive combination (ideal)" : "✗")
 
 # matched ξ_j should equal the coil column sols(:,:,2msing+j)
 err = maximum(1:mcoil) do j
     csol = 2msing + j
-    max(norm(mxi[:, :, j] - sols[:, :, csol]), norm(mdxi[:, :, j] - sols_d[:, :, csol]))
+    max(norm(mxi[:, j, :] - sols[:, csol, keep]), norm(mdxi[:, j, :] - sols_d[:, csol, keep]))
 end
 @printf("[2] max‖ξ_matched − coil column‖ = %.2e   %s\n", err,
     err < 1e-12 ? "✓ matched ξ == bare ideal coil column" : "✗")

diff --git a/benchmarks/verify_gal_match.jl b/benchmarks/verify_gal_match.jl
@@ -1,4 +1,5 @@
-# Piece 2 verification: RPEC outer↔inner matched solution (ForceFreeStates/Solutions/GalerkinIntegration/Match/*).
+# Piece 2 verification: RPEC outer↔inner matched solution (closed profiles under
+# GalerkinIntegration/xi_psi + matching diagnostics under GalerkinIntegration/Match/*).
 #   1. linear-solve residual ‖mat·cof − rmat‖/‖rmat‖
 #   2. matched ξ / ξ′ finiteness
 #   3. edge column == identity basis: each coil drive j must give ξ_edge = e_j (the j-th harmonic),
@@ -11,16 +12,16 @@ h5path = length(ARGS) >= 1 ? ARGS[1] : "examples/DIIID-like_gal_resistive_exampl
 @info "Reading $h5path"
 
 xi, dxi, cout, cin, deltar, eig, resid, sing_psi = h5open(h5path) do f
-    (read(f["ForceFreeStates/Solutions/GalerkinIntegration/Match/xi"]), read(f["ForceFreeStates/Solutions/GalerkinIntegration/Match/dxidpsi"]),
+    (read(f["ForceFreeStates/Solutions/GalerkinIntegration/xi_psi"]), read(f["ForceFreeStates/Solutions/GalerkinIntegration/dxi_psidpsi"]),
         read(f["ForceFreeStates/Solutions/GalerkinIntegration/Match/cout"]), read(f["ForceFreeStates/Solutions/GalerkinIntegration/Match/cin"]),
         read(f["ForceFreeStates/Solutions/GalerkinIntegration/Match/Delta_r"]), read(f["ForceFreeStates/Solutions/GalerkinIntegration/Match/rpec_eig"]),
-        read(f["ForceFreeStates/Solutions/GalerkinIntegration/Match/residual"]), read(f["SingularSurfaces/GalerkinDeltaPrime/rational_psi"]))
+        read(f["ForceFreeStates/Solutions/GalerkinIntegration/Match/residual"]), read(f["ForceFreeStates/Solutions/GalerkinIntegration/rational_psi"]))
 end
 # HDF5 stores ComplexF64 as a compound (re,im); convert if needed
 to_c(a) = eltype(a) <: Complex ? a : map(x -> ComplexF64(x.re, x.im), a)
 xi = to_c(xi); dxi = to_c(dxi); cout = to_c(cout); cin = to_c(cin); deltar = to_c(deltar); eig = to_c(eig)
 
-mpert, ngrid, mcoil = size(xi)
+mpert, mcoil, ngrid = size(xi)
 msing = size(deltar, 1)
 @printf("matched solution: mpert=%d  ngrid=%d  mcoil=%d   msing=%d\n", mpert, ngrid, mcoil, msing)
 
@@ -32,7 +33,7 @@ nbad = count(!isfinite, xi) + count(!isfinite, dxi)
     nbad, maximum(abs, xi), maximum(abs, dxi), nbad == 0 ? "✓" : "✗")
 
 # [3] edge column == identity (last grid point = psihigh edge)
-edge = xi[:, ngrid, :]                       # (mpert mode, mcoil drive)
+edge = xi[:, :, ngrid]                       # (mpert mode, mcoil drive)
 id_err = norm(edge - Matrix{ComplexF64}(I, mpert, mcoil)) / sqrt(mpert)
 @printf("[3] edge basis: ‖ξ(edge) − I‖/√mpert = %.3e   %s\n",
     id_err, id_err < 1e-8 ? "✓ identity-at-edge" : "✗")

diff --git a/benchmarks/verify_gal_solution.jl b/benchmarks/verify_gal_solution.jl
@@ -1,5 +1,6 @@
 # Piece 1 verification: reconstructed gal ξ(ψ) and analytic ξ′(ψ).
-#   1. shapes / finiteness sanity of ForceFreeStates/Solutions/GalerkinIntegration/Solution arrays
+#   1. shapes / finiteness sanity of the raw-basis dump (GalerkinIntegration/Basis; requires
+#      [DEBUG] gal_basis_output = true in the deck)
 #   2. analytic ξ′ vs centered finite-difference of ξ — agree in the smooth interior, diverge at the
 #      packed edge (the spline-endpoint-derivative artifact we deliberately avoid)
 # Usage: julia --project=. verify_gal_solution.jl [path/to/gpec.h5]
@@ -9,20 +10,20 @@ h5path = length(ARGS) >= 1 ? ARGS[1] : "examples/DIIID-like_gal_resistive_exampl
 @info "Reading $h5path"
 
 psi, issing, xi, dxi, sing_psi = h5open(h5path) do f
-    (read(f["ForceFreeStates/Solutions/GalerkinIntegration/Solution/psi"]),
-        read(f["ForceFreeStates/Solutions/GalerkinIntegration/Solution/is_rational"]), read(f["ForceFreeStates/Solutions/GalerkinIntegration/Solution/xi_psi"]),
-        read(f["ForceFreeStates/Solutions/GalerkinIntegration/Solution/dxi_psidpsi"]), read(f["SingularSurfaces/GalerkinDeltaPrime/rational_psi"]))
+    (read(f["ForceFreeStates/Solutions/GalerkinIntegration/Basis/psi"]),
+        read(f["ForceFreeStates/Solutions/GalerkinIntegration/Basis/is_rational"]), read(f["ForceFreeStates/Solutions/GalerkinIntegration/Basis/xi_psi"]),
+        read(f["ForceFreeStates/Solutions/GalerkinIntegration/Basis/dxi_psidpsi"]), read(f["ForceFreeStates/Solutions/GalerkinIntegration/rational_psi"]))
 end
 issing = Bool.(issing)
-mpert, ngrid, nsol = size(xi)
+mpert, nsol, ngrid = size(xi)
 @printf("grid: mpert=%d  ngrid=%d  nsol=%d   psi∈[%.4f, %.4f]\n", mpert, ngrid, nsol, psi[1], psi[end])
 @printf("singular surfaces at psi = %s\n", join((@sprintf("%.4f", p) for p in sing_psi), ", "))
 
 # --- finiteness (skip on-surface points, which are intentionally left zero) ---
 good = .!issing
-nbad = count(!isfinite, xi[:, good, :]) + count(!isfinite, dxi[:, good, :])
+nbad = count(!isfinite, xi[:, :, good]) + count(!isfinite, dxi[:, :, good])
 @printf("non-finite entries (off-surface): %d   |xi|max=%.3e  |dxi|max=%.3e\n",
-    nbad, maximum(abs, xi[:, good, :]), maximum(abs, dxi[:, good, :]))
+    nbad, maximum(abs, xi[:, :, good]), maximum(abs, dxi[:, :, good]))
 
 # --- analytic ξ′ vs centered finite difference of ξ ---
 # For each column, centered diff at interior grid points (using off-surface neighbours), compared to the
@@ -40,8 +41,8 @@ for isol in 1:nsol
         (h1 <= 0 || h2 <= 0) && continue
         for m in 1:mpert
             # nonuniform centered difference
-            fd = (xi[m, ip+1, isol] - xi[m, ip-1, isol]) / (h2 + h1)
-            an = dxi[m, ip, isol]
+            fd = (xi[m, isol, ip+1] - xi[m, isol, ip-1]) / (h2 + h1)
+            an = dxi[m, isol, ip]
             scale = max(abs(an), abs(fd), 1e-30)
             push!(relerr, abs(fd - an) / scale)
             push!(dsurf, dist_sing(psi[ip]))

diff --git a/docs/development/architecture.md b/docs/development/architecture.md
@@ -62,6 +62,7 @@ GPEC consists of **eight main modules** organized in `src/`:
    - Identifies singular surfaces where ξ·∇ψ = 0
    - Key files:
      - `ForceFreeStatesStructs.jl` - Core data structures
+     - `Result.jl` - `ForceFreeStatesResult`, the published solve product every downstream stage reads
      - `Ode.jl` - ODE solver for Euler-Lagrange equations
      - `Sing.jl` - Singular point handling and layer analysis
      - `Fourfit.jl` - Fourier fitting routines
@@ -146,7 +147,7 @@ The complete GPEC analysis pipeline:
    - Compute Δ' at each singular surface
    - Calculate potential and kinetic energies
    - Check Mercier and ballooning stability criteria
-   - Outputs: Eigenmode structure ξ(ψ,θ)
+   - Outputs: `ForceFreeStatesResult` carrying the eigenmode structure ξ(ψ,θ) and the per-integrator products
 
 4. **Perturbed Equilibrium** (GPEC-style):
    - Load external forcing data (coil fields, RMP configuration)
@@ -176,8 +177,14 @@ The complete GPEC analysis pipeline:
 ### Stability
 - `SingType` - Singular surface data including:
   - Rational surface location (ψ, ρ, q = m/n, dq/dψ)
-  - Δ' (tearing stability parameter) — **stub**; the valid Δ' is `ForceFreeStatesInternal.delta_prime_matrix`
+  - Δ' (tearing stability parameter) — **stub**; the valid Δ' is `ForceFreeStatesResult.delta_prime.matrix`
   - Asymptotic solution bases at the inner-layer boundaries
+- `ForceFreeStatesResult` - Published product of a solve: mode space, metric/matrix fits, singular
+  surfaces, and the per-integrator products (ξ solution and its basis, free-boundary energies,
+  STRIDE Δ', Galerkin solve). Optional products are `nothing` when the integrator that ran cannot
+  supply them, and consumers warn-and-skip via `require` / `require_solution`.
+- `ForceFreeStatesInternal` - Solve-time scratch; does not cross a module boundary once the result
+  is built
 
 ### Perturbed Equilibrium
 - `PerturbedEquilibriumControl` - User-facing TOML configuration parameters