LocalStability - API! - Solve a LocalStabilityProblem from its own deck section

.pre-commit-config.yaml

@@ -77,5 +77,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|gal_flag|nstep|diagnose_ca|diagnose|ksing|power_bp|power_b|power_r|power_rc|use_galgrid)\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|nstep|diagnose_ca|diagnose|ksing|power_bp|power_b|power_r|power_rc|use_galgrid|local_stability_flag)\s*='
         files: ^(examples/.*\.toml|test/test_data/.*\.toml)$

benchmarks/benchmark_against_fortran_run.jl

@@ -409,8 +409,9 @@ function write_gpec_toml_coil(
         println(io, "[Wall]")
         println(io, "shape = \"nowall\"")
         println(io)
+        println(io, "[LocalStability]")
+        println(io)
         println(io, "[ForceFreeStates]")
-        println(io, "local_stability_flag = true")
         println(io, "vac_flag = true")
         @printf(io, "qlow     = %.4f\n", p.qlow)
         # Match the Fortran sas_flag truncation: integration stops at q = qhigh
@@ -472,8 +473,9 @@ function write_gpec_toml_file(
         println(io, "[Wall]")
         println(io, "shape = \"nowall\"")
         println(io)
+        println(io, "[LocalStability]")
+        println(io)
         println(io, "[ForceFreeStates]")
-        println(io, "local_stability_flag = true")
         println(io, "vac_flag = true")
         @printf(io, "qlow     = %.4f\n", p.qlow)
         @printf(io, "dW_edge_scan_start = %.4f\n", p.psiedge)

benchmarks/benchmark_diiid_ideal_ntv_torque.jl

@@ -128,8 +128,9 @@ flux_surface_rtol = 1e-7
 [Wall]
 shape = "nowall"
 
+[LocalStability]
+
 [ForceFreeStates]
-local_stability_flag = true
 vac_flag = true
 force_termination = false
 

benchmarks/benchmark_diiid_kinetic_stability.jl

@@ -135,8 +135,9 @@ flux_surface_rtol = 1e-7
 [Wall]
 shape = "nowall"
 
+[LocalStability]
+
 [ForceFreeStates]
-local_stability_flag = true
 vac_flag = true
 force_termination = true       # Skip PE+KF post-processing — we only need FFS eigenvalues
 

docs/development/architecture.md

@@ -72,7 +72,8 @@ Splines are provided by the external `FastInterpolations` package rather than by
 5. **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`
+   - `Problem.jl` - `solve(LocalStabilityProblem(equil), Ballooning())` → `LocalStabilityResult`; a deck's `[LocalStability]` section runs the same solve
+   - Kernels: `compute_ballooning_stability!`, `ballooning_alpha_boundary`
    - Status: Stable
 
 ### Perturbed Equilibrium Modules

docs/src/api.md

@@ -125,6 +125,20 @@ kinetic profiles come from `profile_file`, relative to the run directory.
 Kinetic runs (`kinetic_factor > 0`) take `Forward()` and the `"fixed"` kinetic source;
 the `"calculated"` source needs the `[KineticForces]` profiles and remains TOML-driven.
 
+## Local stability
+
+The Mercier, resistive interchange and ballooning profiles are their own stage, posed on the
+equilibrium alone. A [`LocalStabilityProblem`](@ref) holds the equilibrium, the
+[`Ballooning`](@ref) algorithm holds the scan knobs, and `solve` returns a
+[`LocalStabilityResult`](@ref) with the profiles as named fields. A deck runs the same solve
+through its `[LocalStability]` section, whose keys are the `Ballooning` keywords.
+
+```julia
+ls = solve(LocalStabilityProblem(eq), Ballooning(; n_scan=64))
+ls.psi, ls.D_I, ls.D_R, ls.ballooning_Delta_prime     # profiles on the equilibrium ψ grid
+ls.ballooning_psi, ls.alpha, ls.alpha_critical        # first ballooning boundary, NaN outside its scan window
+```
+
 ## Entry points
 
 ```@docs

docs/src/stability.md

@@ -156,8 +156,11 @@ integrator = "galerkin"
 
 ## Local stability: Mercier and ballooning (s–α)
 
-Setting `local_stability_flag = true` in `[ForceFreeStates]` runs a local high-``n``
-stability scan over every flux surface, in addition to the global ideal analysis above.
+A `[LocalStability]` section in the deck runs a local high-``n`` stability scan over every
+flux surface, as its own stage before the global ideal analysis above. The section's keys are
+the [`Ballooning`](@ref) keywords (the ballooning angle, ODE tolerance and the α-boundary scan
+knobs), and a script runs the same stage with `solve(LocalStabilityProblem(eq), Ballooning())`
+(see the [scripting API](api.md)).
 For the derivation and implementation details behind these diagnostics, see
 [Ballooning and Mercier Local Stability](ballooning.md).
 Three diagnostics are produced and stored under the `LocalStability/` HDF5 group, each a profile
@@ -176,6 +179,10 @@ in normalized poloidal flux ``\psi``:
   taking the jump in the logarithmic derivative of the solution between the two asymptotic
   ends.
 
+With `alpha_boundary_scan = true` the stage also writes the critical ballooning
+``\alpha`` boundary: `LocalStability/ballooning_psi`, `LocalStability/alpha` and
+`LocalStability/alpha_critical` (see the [scripting API](api.md)).
+
 !!! note "Two different Δ' quantities"
     `LocalStability/ballooning_Delta_prime` is the **local high-``n`` ballooning** index and is
     distinct from the **resistive tearing** ``\Delta'`` described in the next section, which
@@ -379,9 +386,6 @@ numsteps_init     = 200    # initial step budget per chunk
 numunorms_init    = 50     # renorm checkpoint budget
 reltol            = 1e-6   # ODE relative tolerance
 
-# Local stability
-local_stability_flag = false  # scan Mercier D_I, resistive D_R, and ballooning Δ' over ψ
-
 # Output
 verbose              = true
 write_outputs_to_HDF5 = true

examples/DIIID-like_SLAYER_example/gpec.toml

@@ -28,9 +28,19 @@ triangularity = 0.5            # Triangularity parameter for wall shapes
 corner_sharpness = 0.05        # Sharpness of wall corners (try 0.05 as initial value)
 equal_arc_wall = true          # Equal arc length distribution of nodes on wall
 
+# Presence runs the LocalStability stage (Mercier D_I, resistive D_R, ballooning Δ', first α boundary); keys are the Ballooning keywords
+[LocalStability]
+theta_k = 0.0                 # Ballooning angle θ_k, used by both the profile scan and the α-boundary scan
+ode_tol = 1e-8                # Relative and absolute tolerance of the ballooning ODE integration behind Δ'
+alpha_boundary_scan = true    # Run the first ballooning α-boundary scan (a root find per ψ surface)
+n_scan = 24                   # Resolution of the α scan for the first boundary
+max_alpha_scale = 8.0         # Upper end of the α search, as a multiple of the experimental α
+alpha_tol = 1e-3              # Bisection tolerance of the critical α, in the same multiples
+psi_window = [0.1, 0.99]      # ψ_N range of the α scan; surfaces outside it are NaN
+verbose = true                # Print progress and per-surface α-scan failures
+
 [ForceFreeStates]
 force_termination = true       # Run FFS + SLAYER, skip PerturbedEquilibrium
-local_stability_flag = true    # Perform local stability analysis (Mercier and ballooning) across the ψ profile
 vac_flag = true                # Compute plasma, vacuum, and total energies for free-boundary modes
 
 dW_edge_scan_start = 0.99      # Edge dW(ψ) diagnostic scan band [dW_edge_scan_start, psilim]; set ≥ psilim to disable

examples/DIIID-like_error_field_example/gpec.toml

@@ -32,8 +32,18 @@ triangularity = 0.5            # Triangularity parameter for wall shapes
 corner_sharpness = 0.05        # Sharpness of wall corners (try 0.05 as initial value)
 equal_arc_wall = true          # Equal arc length distribution of nodes on wall
 
+# Presence runs the LocalStability stage (Mercier D_I, resistive D_R, ballooning Δ', first α boundary); keys are the Ballooning keywords
+[LocalStability]
+theta_k = 0.0                 # Ballooning angle θ_k, used by both the profile scan and the α-boundary scan
+ode_tol = 1e-8                # Relative and absolute tolerance of the ballooning ODE integration behind Δ'
+alpha_boundary_scan = true    # Run the first ballooning α-boundary scan (a root find per ψ surface)
+n_scan = 24                   # Resolution of the α scan for the first boundary
+max_alpha_scale = 8.0         # Upper end of the α search, as a multiple of the experimental α
+alpha_tol = 1e-3              # Bisection tolerance of the critical α, in the same multiples
+psi_window = [0.1, 0.99]      # ψ_N range of the α scan; surfaces outside it are NaN
+verbose = true                # Print progress and per-surface α-scan failures
+
 [ForceFreeStates]
-local_stability_flag = true    # Perform local stability analysis (Mercier and ballooning) across the ψ profile
 vac_flag = true                # Compute plasma, vacuum, and total energies for free-boundary modes
 
 dW_edge_scan_start = 0.99      # Edge dW(ψ) diagnostic scan band [dW_edge_scan_start, psilim]; set ≥ psilim to disable

examples/DIIID-like_gal_resistive_example/gpec.toml

@@ -28,8 +28,18 @@ triangularity = 0.5            # Triangularity parameter for wall shapes
 corner_sharpness = 0.05        # Sharpness of wall corners (try 0.05 as initial value)
 equal_arc_wall = true          # Equal arc length distribution of nodes on wall
 
+# Presence runs the LocalStability stage (Mercier D_I, resistive D_R, ballooning Δ', first α boundary); keys are the Ballooning keywords
+[LocalStability]
+theta_k = 0.0                 # Ballooning angle θ_k, used by both the profile scan and the α-boundary scan
+ode_tol = 1e-8                # Relative and absolute tolerance of the ballooning ODE integration behind Δ'
+alpha_boundary_scan = true    # Run the first ballooning α-boundary scan (a root find per ψ surface)
+n_scan = 24                   # Resolution of the α scan for the first boundary
+max_alpha_scale = 8.0         # Upper end of the α search, as a multiple of the experimental α
+alpha_tol = 1e-3              # Bisection tolerance of the critical α, in the same multiples
+psi_window = [0.1, 0.99]      # ψ_N range of the α scan; surfaces outside it are NaN
+verbose = true                # Print progress and per-surface α-scan failures
+
 [ForceFreeStates]
-local_stability_flag = true    # Perform local stability analysis (Mercier and ballooning) across the ψ profile
 vac_flag = true                # Compute plasma, vacuum, and total energies for free-boundary modes
 
 dW_edge_scan_start = 0.99      # Edge dW(ψ) diagnostic scan band [dW_edge_scan_start, psilim]; set ≥ psilim to disable

examples/DIIID-like_gal_resistive_pe_example/gpec.toml

@@ -29,8 +29,18 @@ triangularity = 0.5            # Triangularity parameter for wall shapes
 corner_sharpness = 0.05        # Sharpness of wall corners (try 0.05 as initial value)
 equal_arc_wall = true          # Equal arc length distribution of nodes on wall
 
+# Presence runs the LocalStability stage (Mercier D_I, resistive D_R, ballooning Δ', first α boundary); keys are the Ballooning keywords
+[LocalStability]
+theta_k = 0.0                 # Ballooning angle θ_k, used by both the profile scan and the α-boundary scan
+ode_tol = 1e-8                # Relative and absolute tolerance of the ballooning ODE integration behind Δ'
+alpha_boundary_scan = true    # Run the first ballooning α-boundary scan (a root find per ψ surface)
+n_scan = 24                   # Resolution of the α scan for the first boundary
+max_alpha_scale = 8.0         # Upper end of the α search, as a multiple of the experimental α
+alpha_tol = 1e-3              # Bisection tolerance of the critical α, in the same multiples
+psi_window = [0.1, 0.99]      # ψ_N range of the α scan; surfaces outside it are NaN
+verbose = true                # Print progress and per-surface α-scan failures
+
 [ForceFreeStates]
-local_stability_flag = true    # Perform local stability analysis (Mercier and ballooning) across the ψ profile
 vac_flag = true                # Compute plasma, vacuum, and total energies for free-boundary modes
 
 dW_edge_scan_start = 0.99      # Edge dW(ψ) diagnostic scan band [dW_edge_scan_start, psilim]; set ≥ psilim to disable

examples/DIIID-like_ideal_example/gpec.toml

@@ -27,8 +27,18 @@ triangularity = 0.5            # Triangularity parameter for wall shapes
 corner_sharpness = 0.05        # Sharpness of wall corners (try 0.05 as initial value)
 equal_arc_wall = true          # Equal arc length distribution of nodes on wall
 
+# Presence runs the LocalStability stage (Mercier D_I, resistive D_R, ballooning Δ', first α boundary); keys are the Ballooning keywords
+[LocalStability]
+theta_k = 0.0                 # Ballooning angle θ_k, used by both the profile scan and the α-boundary scan
+ode_tol = 1e-8                # Relative and absolute tolerance of the ballooning ODE integration behind Δ'
+alpha_boundary_scan = true    # Run the first ballooning α-boundary scan (a root find per ψ surface)
+n_scan = 24                   # Resolution of the α scan for the first boundary
+max_alpha_scale = 8.0         # Upper end of the α search, as a multiple of the experimental α
+alpha_tol = 1e-3              # Bisection tolerance of the critical α, in the same multiples
+psi_window = [0.1, 0.99]      # ψ_N range of the α scan; surfaces outside it are NaN
+verbose = true                # Print progress and per-surface α-scan failures
+
 [ForceFreeStates]
-local_stability_flag = true    # Perform local stability analysis (Mercier and ballooning) across the ψ profile
 vac_flag = true                # Compute plasma, vacuum, and total energies for free-boundary modes
 
 dW_edge_scan_start = 0.99      # Edge dW(ψ) diagnostic scan band [dW_edge_scan_start, psilim]; set ≥ psilim to disable

examples/DIIID-like_ideal_example_IMAS/gpec.toml

@@ -24,7 +24,6 @@ shape = "nowall"               # Wall shape (nowall, conformal, elliptical, dee,
 nn_low               = 1       # Smallest toroidal mode number to include
 nn_high              = 1       # Largest toroidal mode number to include
 vac_flag             = true    # Compute plasma, vacuum, and total energies for free-boundary modes
-local_stability_flag = false   # Perform local stability analysis (Mercier and ballooning) across the ψ profile
 write_outputs_to_HDF5 = false  # Write outputs to the HDF5 file
 verbose              = false   # Enable verbose logging
 delta_mhigh          = 16      # Expands upper bound of Fourier harmonics

examples/DIIID-like_ideal_example_IMAS/run_imas_example.jl

@@ -29,7 +29,6 @@ const SHARED_CONFIG = Dict(
         "nn_low" => 1,
         "nn_high" => 1,
         "vac_flag" => true,
-        "local_stability_flag" => false,
         "write_outputs_to_HDF5" => false,
         "verbose" => false,
         "delta_mhigh" => 16,

examples/DIIID-like_ideal_multi_n_example/gpec.toml

@@ -28,8 +28,18 @@ triangularity = 0.5            # Triangularity parameter for wall shapes
 corner_sharpness = 0.05        # Sharpness of wall corners (try 0.05 as initial value)
 equal_arc_wall = true          # Equal arc length distribution of nodes on wall
 
+# Presence runs the LocalStability stage (Mercier D_I, resistive D_R, ballooning Δ', first α boundary); keys are the Ballooning keywords
+[LocalStability]
+theta_k = 0.0                 # Ballooning angle θ_k, used by both the profile scan and the α-boundary scan
+ode_tol = 1e-8                # Relative and absolute tolerance of the ballooning ODE integration behind Δ'
+alpha_boundary_scan = true    # Run the first ballooning α-boundary scan (a root find per ψ surface)
+n_scan = 24                   # Resolution of the α scan for the first boundary
+max_alpha_scale = 8.0         # Upper end of the α search, as a multiple of the experimental α
+alpha_tol = 1e-3              # Bisection tolerance of the critical α, in the same multiples
+psi_window = [0.1, 0.99]      # ψ_N range of the α scan; surfaces outside it are NaN
+verbose = true                # Print progress and per-surface α-scan failures
+
 [ForceFreeStates]
-local_stability_flag = true    # Perform local stability analysis (Mercier and ballooning) across the ψ profile
 vac_flag = true                # Compute plasma, vacuum, and total energies for free-boundary modes
 
 dW_edge_scan_start = 0.99      # Edge dW(ψ) diagnostic scan band [dW_edge_scan_start, psilim]; set ≥ psilim to disable

examples/DIIID-like_riccati_deltaprime_example/gpec.toml

@@ -28,8 +28,18 @@ triangularity = 0.5            # Triangularity parameter for wall shapes
 corner_sharpness = 0.05        # Sharpness of wall corners (try 0.05 as initial value)
 equal_arc_wall = true          # Equal arc length distribution of nodes on wall
 
+# Presence runs the LocalStability stage (Mercier D_I, resistive D_R, ballooning Δ', first α boundary); keys are the Ballooning keywords
+[LocalStability]
+theta_k = 0.0                 # Ballooning angle θ_k, used by both the profile scan and the α-boundary scan
+ode_tol = 1e-8                # Relative and absolute tolerance of the ballooning ODE integration behind Δ'
+alpha_boundary_scan = true    # Run the first ballooning α-boundary scan (a root find per ψ surface)
+n_scan = 24                   # Resolution of the α scan for the first boundary
+max_alpha_scale = 8.0         # Upper end of the α search, as a multiple of the experimental α
+alpha_tol = 1e-3              # Bisection tolerance of the critical α, in the same multiples
+psi_window = [0.1, 0.99]      # ψ_N range of the α scan; surfaces outside it are NaN
+verbose = true                # Print progress and per-surface α-scan failures
+
 [ForceFreeStates]
-local_stability_flag = true    # Perform local stability analysis (Mercier and ballooning) across the ψ profile
 vac_flag = true                # Compute plasma, vacuum, and total energies for free-boundary modes
 force_termination = true       # Terminate after force-free states (skip perturbed equilibrium calculations)
 

examples/DIIID-like_truncation_example/gpec.toml

@@ -21,7 +21,6 @@ flux_surface_rtol = 1e-10      # Relative tolerance of the ODE solves that trace
 force_termination = false      # Terminate after equilibrium setup (skip stability calculations)
 
 [ForceFreeStates]
-local_stability_flag = false   # Perform local stability analysis (Mercier and ballooning) across the ψ profile
 vac_flag = true                # Compute plasma, vacuum, and total energies for free-boundary modes
 force_termination = true       # Terminate after force-free states (skip perturbed equilibrium calculations)
 

examples/LAR_beta_scan/gpec.toml

@@ -33,8 +33,18 @@ mtheta           = 128         # TJ-analytic internal poloidal grid resolution (
 shape = "conformal"            # Wall shape (nowall, conformal, elliptical, dee, mod_dee, filepath)
 standoff_fraction = 20         # Wall standoff in plasma minor radii (conformal, dee) (20× → effectively no wall)
 
+# Presence runs the LocalStability stage (Mercier D_I, resistive D_R, ballooning Δ', first α boundary); keys are the Ballooning keywords
+[LocalStability]
+theta_k = 0.0                 # Ballooning angle θ_k, used by both the profile scan and the α-boundary scan
+ode_tol = 1e-8                # Relative and absolute tolerance of the ballooning ODE integration behind Δ'
+alpha_boundary_scan = true    # Run the first ballooning α-boundary scan (a root find per ψ surface)
+n_scan = 24                   # Resolution of the α scan for the first boundary
+max_alpha_scale = 8.0         # Upper end of the α search, as a multiple of the experimental α
+alpha_tol = 1e-3              # Bisection tolerance of the critical α, in the same multiples
+psi_window = [0.1, 0.99]      # ψ_N range of the α scan; surfaces outside it are NaN
+verbose = true                # Print progress and per-surface α-scan failures
+
 [ForceFreeStates]
-local_stability_flag = true    # Perform local stability analysis (Mercier and ballooning) across the ψ profile
 vac_flag = true                # Compute plasma, vacuum, and total energies for free-boundary modes
 
 qlow         = 1.02            # Integration initiated at q determined by min(q0, qlow)