matt-pharr · Oct 6, 2026 · Oct 6, 2026 · Oct 6, 2026 · Oct 6, 2026 · Oct 6, 2026
diff --git a/REFACTOR_PLAN.md b/REFACTOR_PLAN.md
diff --git a/benchmarks/scan_resistivity_m2.jl b/benchmarks/scan_resistivity_m2.jl
diff --git a/benchmarks/scan_rotation_m2.jl b/benchmarks/scan_rotation_m2.jl
diff --git a/benchmarks/verify_gal_ideal.jl b/benchmarks/verify_gal_ideal.jl
@@ -9,7 +9,7 @@ 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, 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["SingularSurfaces/Match/cout"])), to_c(read(f["SingularSurfaces/Match/Delta_r"])),
         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"])),

diff --git a/benchmarks/verify_gal_match.jl b/benchmarks/verify_gal_match.jl
@@ -1,5 +1,5 @@
 # Piece 2 verification: RPEC outer↔inner matched solution (closed profiles under
-# GalerkinIntegration/xi_psi + matching diagnostics under GalerkinIntegration/Match/*).
+# GalerkinIntegration/xi_psi + matching diagnostics under SingularSurfaces/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),
@@ -13,9 +13,9 @@ h5path = length(ARGS) >= 1 ? ARGS[1] : "examples/DIIID-like_gal_resistive_exampl
 
 xi, dxi, cout, cin, deltar, eig, resid, sing_psi = h5open(h5path) do f
     (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["ForceFreeStates/Solutions/GalerkinIntegration/rational_psi"]))
+        read(f["SingularSurfaces/Match/cout"]), read(f["SingularSurfaces/Match/cin"]),
+        read(f["SingularSurfaces/Match/Delta_r"]), read(f["SingularSurfaces/Match/rpec_eig"]),
+        read(f["SingularSurfaces/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)

diff --git a/docs/development/architecture.md b/docs/development/architecture.md
@@ -63,7 +63,7 @@ Splines are provided by the external `FastInterpolations` package rather than by
      - `Surfaces/` - Singular-surface finding, Frobenius asymptotics, and GGJ coefficients
      - `Riccati/` - Chunked fundamental-matrix (STRIDE) driver and Δ' boundary-value problem
      - `Galerkin/` - RDCON outer-region singular Galerkin Δ' solver
-     - `Matching/` - Outer↔inner resistive matching (`DeltaPrimeData`, `resonant_match_rpec`)
+     - `Matching/` - Outer↔inner resistive matching (`DeltaPrimeData`, `MatchProblem`, `MatchResult`)
      - `Fourfit.jl` - Fourier fitting routines (`MatrixSplines`)
      - `FixedBoundaryStability.jl` - Fixed boundary analysis
      - `Free.jl` - Free boundary stability

diff --git a/docs/development/hdf5-conventions.md b/docs/development/hdf5-conventions.md
@@ -40,9 +40,9 @@ 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, the gal surface list, debug-gated `Basis/`), `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/` |
+| `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), the inner-layer `Match/` diagnostics, `Kinetic/` |
 | `PerturbedEquilibrium/` | `ForcingModes/`, `Response/`, `ResponseMatrices/`, `SingularCoupling/`, `Energies/`, control-surface spectra |
 | `KineticForces/` | `<method>/` (torque/energy profiles, `EnergyIntegrals/`, `KineticMatrices/`); multi-ion runs add `PerSpecies/<species>/<method>/` with the same per-method layout, summing to the top-level total |
 | `ErrorFields/` | `CoilSensitivities/` (per-coil-set control-surface spectra and their rigid shift/tilt derivatives, `DominantMode/` full-window projection); `MonteCarlo/` (intrinsic and corrected `\|δ\|` histograms over the sampled tolerances, per batch and averaged); `Risk/` (threshold density, `P(lock\|δ)`, locking probabilities, `ToleranceScan/`); `NTV/` (correction-coil overlap and NTV torque couplings per kAt) |
@@ -53,7 +53,7 @@ Reserved (documented, not yet written): `ForceFreeStates/Solutions/RiccatiIntegr
 
 ## Metadata contract (self-describing datasets)
 
-Every dataset outside `Input/` (raw snapshot) and `GalerkinIntegration/Match/` (debug-only) must answer "what is this, in what units, plotted against what" without opening the source — enforced by `test/runtests_h5_schema.jl`:
+Every dataset outside `Input/` (raw snapshot) and `SingularSurfaces/Match/` (debug-only) must answer "what is this, in what units, plotted against what" without opening the source — enforced by `test/runtests_h5_schema.jl`:
 
 - **`long_name`** — plain-text physics description.
 - **`units`** — SI string (`"T"`, `"Wb/rad"`, `"A"`, `"m"`, `"J"`, `"N*m"`, `"Hz"`, `"Ohm*m"`); `"1"` for dimensionless (CF convention). Normalized quantities state the normalization in `long_name` (e.g. the power-normalized stability energies are per unit ⟨|ξ|²⟩, not joules).

diff --git a/docs/src/api.md b/docs/src/api.md
@@ -78,18 +78,56 @@ Which products each formalism can supply differs; a result carries `nothing` in
 its integrator does not produce and consumers warn and skip rather than erroring. See the
 [Stability Analysis](stability.md) page for the result struct and its capability gates.
 
-Inner-layer matching is requested with the integrator-agnostic `match` keyword, which closes
-the basis with a resistive layer solution instead of the ideal jump condition:
+## Inner-layer matching
+
+Inner-layer matching is its own problem, posed on a finished solve: a `MatchProblem` holds
+the outer Δ′ the solve published plus the per-surface layer parameters, and the inner-layer
+model passed to `solve` computes the layer response at the prescribed rotation. Solving it
+returns a new result with the closure changed from `:ideal` to `:matched` and the
+eigenfunctions replaced — the expensive outer solve is reused, so layer-parameter scans
+cost one cheap match solve per point:
+
+```julia
+ffs = solve(eq, Galerkin(; rpec_flag=true, cut_solution=true); nn=1)
+
+matched = solve(MatchProblem(ffs; eta=[1e-6, 2e-6], rho=[1e-7, 1e-7], rotation=[0.0, 0.0]), GGJModel())
+@assert matched.closure === :matched
+
+# A rotation scan reuses the one outer solve:
+bpens = [solve(MatchProblem(ffs; eta=[1e-6, 2e-6], rho=[1e-7, 1e-7], rotation=[f, f]), GGJModel()).bpen
+         for f in 0.0:50.0:500.0]
+```
+
+The per-surface η/ρ/rotation can also be derived from the kinetic profiles attached to the
+equilibrium (`layer_parameters`), with the explicit vectors as overrides. The problem needs
+a Δ′ payload with coil-response columns, so it accepts Galerkin (`rpec_flag=true`) and
+Riccati (`vac_flag=true`) results; a Riccati-fed match fills `bpen` and the resonant data but keeps
+`solution === nothing` (no outer basis is retained). Only a closure-capable model is
+accepted — `GGJModel()` today; `SLAYERModel()` is slab-only and drives the free-eigenvalue tearing
+solve instead.
+
+## Tearing stability
+
+The free-eigenvalue tearing solve is the second flavor of inner-layer matching: instead of
+prescribing the layer rotation, a `TearingProblem` holds the outer Δ′ fixed and root-finds
+the growth rate where the inner-layer response matches it. The same model slot applies —
+`SLAYERModel()` is the slab layer that exists for exactly this problem. `GGJModel()` also runs
+through the scan, but GGJ growth-rate extraction is not implemented yet: its γ are placeholders.
 
 ```julia
-ffs = solve(eq, Galerkin(); nn=1,
-    match=ResistiveMatch(; eta=[1e-6, 2e-6], rho=[1e-7, 1e-7], rotation=[0.0, 0.0]))
-@assert ffs.closure === :matched
+attach_kinetic_profiles!(eq, "kin.h5")                      # n, T, ω for the layer parameters
+ffs  = solve(eq, Riccati(); nn=1, vac_flag=true)            # Δ′ matrix for the dispersion
+tear = solve(TearingProblem(ffs; coupling_mode=:coupled), SLAYERModel())
+tear.gamma_Hz, tear.rational_q                               # root-found rates per surface
 ```
 
-Only the Galerkin formalism implements the match today; requesting one from `Forward` or
-`Riccati` errors. Kinetic runs (`kinetic_factor > 0`) need the `[KineticForces]` profiles and
-remain TOML-driven.
+Keyword arguments of `TearingProblem` are the `[SLAYER]` deck section's procedure knobs
+(scan mode and Q-domain, coupling mode, critical-Δ convention, extraction filters);
+kinetic profiles come from `profile_file` or, when none is named, from the profiles
+attached to the equilibrium.
+
+Kinetic runs (`kinetic_factor > 0`) need the `[KineticForces]` profiles and remain
+TOML-driven.
 
 ## Entry points
 

diff --git a/docs/src/galerkin.md b/docs/src/galerkin.md
@@ -17,8 +17,7 @@ whichever formalism produced it.  These are the outer-region inputs to resistive
 Select it with `integrator = "galerkin"` in `[ForceFreeStates]`. It replaces the radial ODE
 integration rather than supplementing it: the run computes its own vacuum response at the
 control surface (when `vac_flag`) and produces no free-boundary energies or ODE trace.
-Setting `gal_match_flag` additionally matches the inner layer, giving a driven ξ solution that
-`PerturbedEquilibrium` consumes in place of a forward solution.
+Setting `gal_match_flag` additionally matches the inner layer, giving a driven ξ solution.
 
 The implementation lives in `src/ForceFreeStates/Galerkin/`:
 
@@ -28,14 +27,13 @@ The implementation lives in `src/ForceFreeStates/Galerkin/`:
 | `GalerkinGrid.jl` | Packed grid construction and local→global DOF mapping |
 | `GalerkinAssembly.jl` | Element-level assembly: Hermite basis, Gauss-Lobatto stiffness, resonant and extension cells, boundary conditions |
 | `GalerkinSolution.jl` | Reconstruct ξ(ψ) and analytic ξ′(ψ) on the gal-native grid |
-| `GalerkinMatch.jl` | DRIVEN/RPEC outer↔inner asymptotic matching, whose matched solution PerturbedEquilibrium consumes |
 | `GalerkinSolve.jl` | Top-level driver `galerkin_solve`, banded solve, Δ′ extraction, PEST-3 blocks, HDF5 output |
 
 ## API Reference
 
 ```@autodocs
 Modules = [GeneralizedPerturbedEquilibrium.ForceFreeStates]
-Pages = ["Galerkin/GalerkinStructs.jl", "Galerkin/GalerkinGrid.jl", "Galerkin/GalerkinAssembly.jl", "Galerkin/GalerkinSolution.jl", "Galerkin/GalerkinMatch.jl", "Galerkin/GalerkinSolve.jl"]
+Pages = ["Galerkin/GalerkinStructs.jl", "Galerkin/GalerkinGrid.jl", "Galerkin/GalerkinAssembly.jl", "Galerkin/GalerkinSolution.jl", "Galerkin/GalerkinSolve.jl"]
 ```
 
 ## See also

diff --git a/docs/src/stability.md b/docs/src/stability.md
@@ -144,8 +144,8 @@ zeroing vs GR), not from ODE tolerance; it is present at every thread count.
 the outer region is discretized on packed Hermite-cubic elements and solved as one global banded
 system, giving the RDCON resistive ``\Delta'`` matrix and the PEST-3 matching blocks.  It computes its own vacuum response and
 returns no free-boundary energies, no ODE trace, and no fixed-boundary `crit` scan.  With
-`gal_match_flag` it also matches the inner layer, producing a driven ``\xi`` solution that
-`PerturbedEquilibrium` consumes.  Kinetic runs are not supported.  See
+`gal_match_flag` it also matches the inner layer, producing a driven ``\xi`` solution.
+Kinetic runs are not supported.  See
 `docs/src/galerkin.md` for the solver and its `gal_*` knobs.
 
 Enable with:
@@ -292,7 +292,7 @@ The Galerkin Δ′ solver (`src/ForceFreeStates/Galerkin/`) is documented separa
 
 ```@autodocs
 Modules = [GeneralizedPerturbedEquilibrium.ForceFreeStates]
-Pages = ["ForceFreeStates.jl", "CoreTypes.jl", "Surfaces/Types.jl", "Riccati/Types.jl", "Matching/DeltaPrime.jl", "Result.jl", "Surfaces/Resist.jl", "Surfaces/ResistEval.jl", "Matching/ResonantMatch.jl", "EulerLagrange.jl", "Surfaces/Finding.jl", "Surfaces/Asymptotics.jl", "Fourfit.jl", "Kinetic.jl", "FixedBoundaryStability.jl", "Utils.jl", "Free.jl", "Riccati/Propagators.jl", "Riccati/Crossings.jl", "Riccati/DeltaPrimeBVP.jl", "Riccati/Driver.jl"]
+Pages = ["ForceFreeStates.jl", "CoreTypes.jl", "Surfaces/Types.jl", "Riccati/Types.jl", "Matching/DeltaPrime.jl", "Result.jl", "Surfaces/ResistEval.jl", "Matching/ResonantMatch.jl", "Matching/LayerParameters.jl", "Matching/MatchProblem.jl", "EulerLagrange.jl", "Surfaces/Finding.jl", "Surfaces/Asymptotics.jl", "Fourfit.jl", "Kinetic.jl", "FixedBoundaryStability.jl", "Utils.jl", "Free.jl", "Riccati/Propagators.jl", "Riccati/Crossings.jl", "Riccati/DeltaPrimeBVP.jl", "Riccati/Driver.jl"]
 ```
 
 ## Example usage