FFS - REFACTOR - Extract local stability into standalone LocalStability module (refactor plan PR 2)

REFACTOR_PLAN.md

@@ -670,7 +670,11 @@ manual smoke: run the 4-line UX from the Context section in a REPL against
   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`
+- [ ] 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`

docs/development/architecture.md

@@ -12,7 +12,7 @@ This workflow is reflected in the modular structure and data flow.
 
 ## Module Structure
 
-GPEC consists of **seven main modules** organized in `src/`:
+GPEC consists of **eight main modules** organized in `src/`:
 
 ### Foundation Modules
 
@@ -67,18 +67,23 @@ GPEC consists of **seven main modules** organized in `src/`:
      - `Fourfit.jl` - Fourier fitting routines
      - `FixedBoundaryStability.jl` - Fixed boundary analysis
      - `Free.jl` - Free boundary stability
-     - `Ballooning.jl` - Local stability scan: Mercier D_I, resistive interchange D_R, and high-n ballooning Δ' (s–α). Replaces the former standalone `Mercier.jl`.
    - Status: Stable, core DCON functionality implemented
 
+6. **LocalStability** (`src/LocalStability/`) - Local high-n stability
+   - `Ballooning.jl` - Local stability scan: Mercier D_I, resistive interchange D_R, and high-n ballooning Δ' (s–α). Replaces the former standalone `Mercier.jl`.
+   - Depends only on Equilibrium (plus math libraries); carries no stability-solver state
+   - Main entry points: `compute_local_stability`, `ballooning_alpha_boundary`
+   - Status: Stable
+
 ### Perturbed Equilibrium Modules
 
-6. **ForcingTerms** (`src/ForcingTerms/`) - External field specification
+7. **ForcingTerms** (`src/ForcingTerms/`) - External field specification
    - Handles external magnetic field perturbations (coils, RMP, etc.)
    - Supports ASCII and HDF5 forcing data formats
    - `ForcingMode` data structure specifies amplitude and phase for each (m,n) component
    - Status: Complete and functional
 
-7. **PerturbedEquilibrium** (`src/PerturbedEquilibrium/`) - **GPEC-style plasma response**
+8. **PerturbedEquilibrium** (`src/PerturbedEquilibrium/`) - **GPEC-style plasma response**
    - Computes plasma response to external forcing
    - Calculates singular coupling metrics at rational surfaces
    - Key files:
@@ -191,6 +196,7 @@ GeneralizedPerturbedEquilibrium
 ├── Utilities (shared tools)
 │   └── FourierTransforms
 ├── Equilibrium (uses Splines)
+├── LocalStability (uses Equilibrium)
 ├── Vacuum (uses Splines, Equilibrium, Utilities)
 ├── ForcingTerms (data I/O)
 ├── ForceFreeStates (uses Equilibrium, Vacuum, Splines)

docs/src/ballooning.md

@@ -911,3 +911,9 @@ The local-stability output now stores ballooning ``\Delta'`` in the fourth
 `locstab_fs` entry. In the HDF5 output this is written as
 `LocalStability/ballooning_Delta_prime`, distinct from the tearing ``\Delta'`` outputs
 under `SingularSurfaces/` and `PerturbedEquilibrium/SingularCoupling/`.
+
+## API
+
+```@autodocs
+Modules = [GeneralizedPerturbedEquilibrium.LocalStability]
+```

docs/src/stability.md

@@ -270,7 +270,7 @@ The Galerkin Δ′ solver (`src/ForceFreeStates/Galerkin/`) is documented separa
 
 ```@autodocs
 Modules = [GeneralizedPerturbedEquilibrium.ForceFreeStates]
-Pages = ["ForceFreeStates.jl", "ForceFreeStatesStructs.jl", "Ballooning.jl", "Resist.jl", "EulerLagrange.jl", "Sing.jl", "Fourfit.jl", "Kinetic.jl", "FixedBoundaryStability.jl", "Utils.jl", "Free.jl", "Riccati.jl"]
+Pages = ["ForceFreeStates.jl", "ForceFreeStatesStructs.jl", "Resist.jl", "EulerLagrange.jl", "Sing.jl", "Fourfit.jl", "Kinetic.jl", "FixedBoundaryStability.jl", "Utils.jl", "Free.jl", "Riccati.jl"]
 ```
 
 ## Example usage

examples/DIIID-like_ideal_example/analyze_example.jl

@@ -2,7 +2,7 @@ using Pkg;
 Pkg.activate(joinpath(@__DIR__, "../.."))
 using GeneralizedPerturbedEquilibrium, Plots
 using GeneralizedPerturbedEquilibrium: Analysis
-using GeneralizedPerturbedEquilibrium: Equilibrium, ForceFreeStates
+using GeneralizedPerturbedEquilibrium: Equilibrium, ForceFreeStates, LocalStability
 using LaTeXStrings, TOML
 isinteractive() ? plotlyjs() : gr()
 
@@ -29,7 +29,7 @@ Plots.savefig(p_baloo, baloo_path)
 # ----------------------------------------------------------------------
 # Local stability: s-alpha profiles, D_I / ballooning Δ', and 2D s-alpha
 # maps (ported from the former Bal_salpha_delta_di_summary notebook).
-# These recompute directly from the equilibrium via the Ballooning.jl helpers
+# These recompute directly from the equilibrium via the LocalStability helpers
 # rather than reading gpec.h5, so the s-alpha scan can perturb (p', q').
 # ----------------------------------------------------------------------
 
@@ -66,7 +66,7 @@ for geqdsk_file in geqdsk_cases
     alpha_profile = fill(NaN, length(psi_norm))
     for i in eachindex(psi_norm)
         try
-            ref = ForceFreeStates.salpha_reference(i, equil)
+            ref = LocalStability.salpha_reference(i, equil)
             s_profile[i] = ref.s_ref
             alpha_profile[i] = ref.alpha_ref
         catch err
@@ -107,8 +107,7 @@ for geqdsk_file in geqdsk_cases
 
     # D_I (Mercier) and ballooning Δ' profiles
     locstab_fs = zeros(length(psi_norm), 5)
-    ctrl = ForceFreeStates.ForceFreeStatesControl(; verbose=false)
-    ForceFreeStates.compute_ballooning_stability!(ctrl, locstab_fs, equil; compute_delta_prime=true)
+    LocalStability.compute_ballooning_stability!(locstab_fs, equil; compute_delta_prime=true)
 
     delta_prime = Vector(locstab_fs[:, 4])
     di_profile = fill(NaN, length(psi_norm))
@@ -152,7 +151,7 @@ for geqdsk_file in geqdsk_cases
     psi_idx_scan = argmin(abs.(psi_norm .- psi_target))
     s_scales = collect(range(-5.0, 5.0; length=30))
     alpha_scales = collect(range(-5.0, 5.0; length=30))
-    scan = ForceFreeStates.scan_delta_prime_map(psi_idx_scan, equil; theta_k=0.0, s_scales=s_scales, alpha_scales=alpha_scales)
+    scan = LocalStability.scan_delta_prime_map(psi_idx_scan, equil; theta_k=0.0, s_scales=s_scales, alpha_scales=alpha_scales)
 
     scan_qprime = scan.dqdpsi_ref .* scan.s_scales
     scan_pprime = scan.pprime_ref .* scan.alpha_scales
@@ -210,8 +209,8 @@ for geqdsk_file in geqdsk_cases
     Plots.savefig(p_zero, zero_path)
 
     # BALOO-style α vs ψ stability boundaries (1st and 2nd) over the Δ'(ψ, α) map
-    bnd = ForceFreeStates.ballooning_alpha_boundaries(ctrl, equil)
-    dpmap = ForceFreeStates.ballooning_delta_prime_map(ctrl, equil)
+    bnd = LocalStability.ballooning_alpha_boundaries(equil)
+    dpmap = LocalStability.ballooning_delta_prime_map(equil)
     p_alpha_bnd = Analysis.ForceFreeStates.plot_ballooning_alpha_boundaries(bnd, dpmap)
     display(p_alpha_bnd)
     alpha_bnd_path = joinpath(@__DIR__, "alpha_stability_boundary_$(case_label).png")

src/Analysis/ForceFreeStates.jl

@@ -584,7 +584,7 @@ arrays leave natural gaps over always-stable surfaces without requiring explicit
 ### Arguments
 
   - `bnd`: NamedTuple with fields `psi`, `alpha`, `alpha_critical1`, `alpha_critical2`
-    (as returned by `ForceFreeStates.ballooning_alpha_boundaries`)
+    (as returned by `LocalStability.ballooning_alpha_boundaries`)
 
 ### Keyword arguments
 
@@ -617,7 +617,7 @@ end
     plot_ballooning_alpha_boundaries(bnd, dpmap; save_path=nothing, psi_min=0.0)
 
 Same diagram drawn over a heatmap of the signed Δ' from
-`ForceFreeStates.ballooning_delta_prime_map`: each surface's Δ'(α) is oriented by the
+`LocalStability.ballooning_delta_prime_map`: each surface's Δ'(α) is oriented by the
 sign of its α=0 (stable) value so that positive is stable everywhere, regridded from
 its native physical α = α_ref*scale onto a shared uniform α axis, and shown with the
 Δ'=0 contour, the extracted boundaries, and the scan cap `max_alpha_scale*α_exp`.

src/ForceFreeStates/ForceFreeStates.jl

@@ -24,7 +24,6 @@ import StaticArrays: @MMatrix
 
 # Include all necessary files
 include("ForceFreeStatesStructs.jl")
-include("Ballooning.jl")
 include("Resist.jl")
 include("EulerLagrange.jl")
 include("Sing.jl")

src/GeneralizedPerturbedEquilibrium.jl

@@ -20,6 +20,11 @@ include("Equilibrium/Equilibrium.jl")
 import .Equilibrium as Equilibrium
 export Equilibrium
 
+# Local high-n stability (Mercier, resistive interchange, ballooning Δ'); depends only on Equilibrium.
+include("LocalStability/LocalStability.jl")
+import .LocalStability as LocalStability
+export LocalStability
+
 include("Vacuum/Vacuum.jl")
 import .Vacuum as Vacuum
 export Vacuum
@@ -64,7 +69,6 @@ include("Rerun.jl")
 # Import ForceFreeStates types and functions needed for main
 using .ForceFreeStates: ForceFreeStatesInternal, ForceFreeStatesControl, DebugSettings, FreeBoundaryResult, OdeState, FourFitVars
 using .ForceFreeStates: sing_lim!, sing_min!, sing_find!, resist_eval_all!, resist_geometry, ResistGeometry
-using .ForceFreeStates: compute_local_stability, compute_ballooning_stability!, ballooning_alpha_boundary, ballooning_alpha_boundaries
 using .ForceFreeStates: make_metric, make_matrix, make_kinetic_matrix
 using .ForceFreeStates: find_kinetic_singular_surfaces!
 using .ForceFreeStates: eulerlagrange_integration, free_run, normalize_eigenfunctions!
@@ -336,9 +340,9 @@ function main_from_inputs(
     locstab = nothing
     ballooning_boundary = (psi=Float64[], alpha=Float64[], alpha_critical=Float64[])
     if ctrl.local_stability_flag
-        locstab = compute_local_stability(ctrl, equil)
+        locstab = LocalStability.compute_local_stability(equil; verbose=ctrl.verbose)
         # First ballooning stability boundary (α vs ψ_N) for BALOO-style diagnostics.
-        ballooning_boundary = ballooning_alpha_boundary(ctrl, equil)
+        ballooning_boundary = LocalStability.ballooning_alpha_boundary(equil; verbose=ctrl.verbose)
     end
 
     # Find all singular surfaces in the equilibrium