From 3c36f2a95f3a8c705dd3dc23cde463b335c6c141 Mon Sep 17 00:00:00 2001 From: Matthew Pharr Date: Fri, 23 Oct 2026 14:40:33 +0200 Subject: [PATCH 1/7] LocalStability - API - Add LocalStabilityProblem, Ballooning and LocalStabilityResult --- docs/src/api.md | 14 ++++ src/GeneralizedPerturbedEquilibrium.jl | 2 + src/LocalStability/LocalStability.jl | 3 + src/LocalStability/Problem.jl | 90 ++++++++++++++++++++++++++ test/runtests.jl | 1 + test/runtests_localstability.jl | 54 ++++++++++++++++ 6 files changed, 164 insertions(+) create mode 100644 src/LocalStability/Problem.jl create mode 100644 test/runtests_localstability.jl diff --git a/docs/src/api.md b/docs/src/api.md index 1f80de680..ee901de94 100644 --- a/docs/src/api.md +++ b/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 diff --git a/src/GeneralizedPerturbedEquilibrium.jl b/src/GeneralizedPerturbedEquilibrium.jl index c03a52542..e4f91c59b 100755 --- a/src/GeneralizedPerturbedEquilibrium.jl +++ b/src/GeneralizedPerturbedEquilibrium.jl @@ -89,6 +89,7 @@ using .ForceFreeStates: AbstractIntegrator, Forward, Riccati, Galerkin using .ForceFreeStates: MatchProblem using .InnerLayer: GGJModel, SLAYERModel using .Tearing.Runner: TearingProblem +using .LocalStability: LocalStabilityProblem, Ballooning, LocalStabilityResult using .Equilibrium: PlasmaEquilibrium using .ForcingTerms: RMPField @@ -1897,6 +1898,7 @@ export main, write_imas export solve, perturbed_equilibrium export PlasmaEquilibrium, EulerLagrangeProblem, Forward, Riccati, Galerkin, ForceFreeStatesResult, RMPField export MatchProblem, TearingProblem, GGJModel, SLAYERModel +export LocalStabilityProblem, Ballooning, LocalStabilityResult include("Precompile.jl") diff --git a/src/LocalStability/LocalStability.jl b/src/LocalStability/LocalStability.jl index eb96780e2..e68c8029a 100644 --- a/src/LocalStability/LocalStability.jl +++ b/src/LocalStability/LocalStability.jl @@ -7,11 +7,14 @@ using FFTW using OrdinaryDiffEq using FastInterpolations using StaticArrays: SVector +import CommonSolve import ..Equilibrium include("Ballooning.jl") +include("Problem.jl") +export LocalStabilityProblem, Ballooning, LocalStabilityResult, mercier_D_I export compute_local_stability, compute_ballooning_stability!, ballooning_alpha_boundary, ballooning_alpha_boundaries end diff --git a/src/LocalStability/Problem.jl b/src/LocalStability/Problem.jl new file mode 100644 index 000000000..6ee45fa3d --- /dev/null +++ b/src/LocalStability/Problem.jl @@ -0,0 +1,90 @@ +# Problem.jl +# +# The LocalStability stage posed as solve(problem, algorithm): Mercier D_I, resistive D_R, +# ballooning Δ' and the first ballooning α boundary over the equilibrium ψ grid. + +""" + LocalStabilityProblem(equil) + +Local high-n stability of an equilibrium: the Mercier and resistive interchange criteria and +the ballooning Δ' on every flux surface, and the first ballooning α boundary. Solve with +`solve(prob, Ballooning())` to get a [`LocalStabilityResult`](@ref). A deck poses the same +problem through its `[LocalStability]` section. + +## Fields + + - `equil::Equilibrium.PlasmaEquilibrium` - The equilibrium the profiles are computed on. +""" +struct LocalStabilityProblem + equil::Equilibrium.PlasmaEquilibrium +end + +""" + Ballooning(; theta_k=0.0, n_scan=24, alpha_boundary_scan=true, verbose=true) + +The local-stability algorithm: the ballooning coefficient system per flux surface, plus the +BALOO-style α-boundary scan. Its keywords are the `[LocalStability]` deck keys. + +## Fields + + - `theta_k::Float64` - Ballooning angle θ_k, used by both the profile scan and the α-boundary scan. + - `n_scan::Int` - Resolution of the α scan for the first boundary (see [`ballooning_alpha_boundary`](@ref)). + - `alpha_boundary_scan::Bool` - Run the α-boundary scan, a root find per ψ surface; `false` leaves the α fields empty. + - `verbose::Bool` - Print progress and per-surface α-scan failures. +""" +@kwdef struct Ballooning + theta_k::Float64 = 0.0 + n_scan::Int = 24 + alpha_boundary_scan::Bool = true + verbose::Bool = true +end + +""" + LocalStabilityResult + +The local-stability profiles returned by `solve(::LocalStabilityProblem, ::Ballooning)`. + +## Fields + + - `psi::Vector{Float64}` - Normalized flux grid of the profiles (the equilibrium `profiles.xs`). + - `D_I::Vector{Float64}` - Mercier interchange criterion `det(d0bar)`; stable where negative. + - `D_R::Vector{Float64}` - Resistive interchange criterion `D_I + (H - 1/2)²` (see [`resistive_interchange_h`](@ref)). + - `ballooning_Delta_prime::Vector{Float64}` - High-n ballooning Δ', zero for ψ > 1. + - `ballooning_psi::Vector{Float64}` - Grid of the α-boundary scan; empty when the scan is off. + - `alpha::Vector{Float64}` - Experimental α on `ballooning_psi`, NaN outside the scan window. + - `alpha_critical::Vector{Float64}` - First-boundary critical α on `ballooning_psi`, NaN outside the scan window or with no boundary. + - `di_psi` - Cubic interpolant of `D_I·ψ` over `psi`, for evaluating `D_I` off the grid (see [`mercier_D_I`](@ref)). +""" +struct LocalStabilityResult{S} + psi::Vector{Float64} + D_I::Vector{Float64} + D_R::Vector{Float64} + ballooning_Delta_prime::Vector{Float64} + ballooning_psi::Vector{Float64} + alpha::Vector{Float64} + alpha_critical::Vector{Float64} + di_psi::S +end + +""" + mercier_D_I(result::LocalStabilityResult, psi) -> Float64 + +Mercier `D_I` at `psi`, interpolated through the smooth `D_I·ψ` profile rather than `D_I` itself. +""" +mercier_D_I(result::LocalStabilityResult, psi::Real) = result.di_psi(psi)[1] / psi + +function CommonSolve.solve(prob::LocalStabilityProblem, alg::Ballooning) + equil = prob.equil + xs = equil.profiles.xs + locstab_fs = zeros(Float64, length(xs), 5) + compute_ballooning_stability!(locstab_fs, equil; theta_k=alg.theta_k, verbose=alg.verbose) + boundary = if alg.alpha_boundary_scan + ballooning_alpha_boundary(equil; theta_k=alg.theta_k, n_scan=alg.n_scan, verbose=alg.verbose) + else + (psi=Float64[], alpha=Float64[], alpha_critical=Float64[]) + end + # A Series spline, not a plain-vector one: the two differ in the last bits off the knots. + di_psi = cubic_interp(xs, Series(locstab_fs[:, 1:1]); extrap=ExtendExtrap()) + return LocalStabilityResult(collect(xs), locstab_fs[:, 1] ./ xs, locstab_fs[:, 2] ./ xs, locstab_fs[:, 4], + boundary.psi, boundary.alpha, boundary.alpha_critical, di_psi) +end diff --git a/test/runtests.jl b/test/runtests.jl index ad3563c91..b3891e63c 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -49,6 +49,7 @@ else include("./runtests_tj_analytic.jl") include("./runtests_kinetic_profiles.jl") include("./runtests_resist_eval.jl") + include("./runtests_localstability.jl") include("./runtests_slayer_params.jl") include("./runtests_slayer_riccati.jl") include("./runtests_slayer_inputs.jl") diff --git a/test/runtests_localstability.jl b/test/runtests_localstability.jl new file mode 100644 index 000000000..ef937cbe9 --- /dev/null +++ b/test/runtests_localstability.jl @@ -0,0 +1,54 @@ +using TOML + +# The LocalStability stage on the coarse Solovev fixture deck (mpsi=16): the named profiles of +# solve(LocalStabilityProblem, Ballooning()) against their kernels. +@testset "LocalStability" begin + GPEC = GeneralizedPerturbedEquilibrium + LS = GPEC.LocalStability + template = joinpath(@__DIR__, "test_data", "regression_solovev_ideal_example") + deck = TOML.parsefile(joinpath(template, "gpec.toml")) + equil = GPEC.Equilibrium.setup_equilibrium( + GPEC.Equilibrium.EquilibriumConfig(deck["Equilibrium"], template), + GPEC.Equilibrium.SolovevConfig(deck["SOL_INPUT"]) + ) + xs = equil.profiles.xs + res = solve(LocalStabilityProblem(equil), Ballooning(; verbose=false)) + + @testset "profiles match their kernels" begin + @test res isa LocalStabilityResult + @test res.psi == xs + for i in eachindex(xs) + di = LS.prepare_ballooning_coefficients(i, equil).di + h = LS.resistive_interchange_h(i, equil) + @test res.D_I[i] ≈ di rtol = 1e-14 + @test res.D_R[i] ≈ di + (h - 0.5)^2 rtol = 1e-12 + @test LS.mercier_D_I(res, xs[i]) ≈ res.D_I[i] rtol = 1e-12 + end + @test all(isfinite, res.ballooning_Delta_prime) + end + + @testset "α boundary is NaN outside its scan window" begin + @test res.ballooning_psi == xs + inside = [LS._in_ballooning_scan_window(p, xs[end]) for p in xs] + @test any(inside) && !all(inside) + @test all(isfinite, res.alpha[inside]) + @test all(isnan, res.alpha[.!inside]) + @test all(isnan, res.alpha_critical[.!inside]) + end + + @testset "the α scan can be switched off" begin + quick = solve(LocalStabilityProblem(equil), Ballooning(; alpha_boundary_scan=false, verbose=false)) + @test isempty(quick.ballooning_psi) && isempty(quick.alpha) && isempty(quick.alpha_critical) + @test quick.D_I == res.D_I + @test quick.ballooning_Delta_prime == res.ballooning_Delta_prime + end + + @testset "solve writes no file" begin + mktempdir() do dir + cd(dir) do + solve(LocalStabilityProblem(equil), Ballooning(; alpha_boundary_scan=false, verbose=false)) + end + @test isempty(readdir(dir)) + end + end +end From b3eb9aed5d2ff0a9982193c06be4cb172367499d Mon Sep 17 00:00:00 2001 From: Matthew Pharr Date: Fri, 23 Oct 2026 15:06:38 +0200 Subject: [PATCH 2/7] =?UTF-8?q?LocalStability=20-=20FEATURE=20-=20Expose?= =?UTF-8?q?=20the=20=CE=B1-scan=20window,=20search=20bound,=20bisection=20?= =?UTF-8?q?and=20ODE=20tolerances=20on=20Ballooning?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- src/LocalStability/Ballooning.jl | 40 +++++++++++++++++++------------- src/LocalStability/Problem.jl | 20 ++++++++++++---- test/runtests_localstability.jl | 13 +++++++++++ 3 files changed, 52 insertions(+), 21 deletions(-) diff --git a/src/LocalStability/Ballooning.jl b/src/LocalStability/Ballooning.jl index 83ad76192..9979b913b 100644 --- a/src/LocalStability/Ballooning.jl +++ b/src/LocalStability/Ballooning.jl @@ -10,9 +10,9 @@ const MATCHING_POINT = 1e-3 # mid-radius and pedestal; the packed axis and far-edge surfaces only add scan cost. const BALLOONING_SCAN_PSI_WINDOW = (0.1, 0.99) -# Effective window [0.1, min(0.99, ψ_edge)]; the grid never extends past psihigh/psi_edge. -_in_ballooning_scan_window(psi::Float64, psi_edge::Float64) = - BALLOONING_SCAN_PSI_WINDOW[1] <= psi <= min(BALLOONING_SCAN_PSI_WINDOW[2], psi_edge) +# Effective window [lo, min(hi, ψ_edge)]; the grid never extends past psihigh/psi_edge. +_in_ballooning_scan_window(psi::Float64, psi_edge::Float64, window=BALLOONING_SCAN_PSI_WINDOW) = + window[1] <= psi <= min(window[2], psi_edge) """ compute_ballooning_stability!(locstab_fs, plasma_eq) @@ -26,6 +26,7 @@ to compute Delta Prime. - `locstab_fs::Matrix{Float64}`: Local stability matrix to store results (modified in place). - `plasma_eq::Equilibrium.PlasmaEquilibrium`: Plasma equilibrium data. + - `ode_tol::Float64`: Relative and absolute tolerance of the ballooning ODE integration. - `verbose::Bool`: Print progress messages. This function modifies `locstab_fs` in place with: @@ -40,6 +41,7 @@ function compute_ballooning_stability!( plasma_eq::Equilibrium.PlasmaEquilibrium; theta_k::Float64=0.0, compute_delta_prime::Bool=true, + ode_tol::Float64=1e-8, verbose::Bool=false ) @@ -64,7 +66,8 @@ function compute_ballooning_stability!( if compute_delta_prime && plasma_eq.profiles.xs[flux_surface_index] <= 1.0 result = integrate_ballooning_ode( coeff_data.ode_coefficient_spline; - theta_k=theta_k + theta_k=theta_k, + tol=ode_tol ) locstab_fs[flux_surface_index, 4] = result.value end @@ -839,7 +842,8 @@ function critical_ballooning_alpha( end """ - ballooning_alpha_boundary(plasma_eq; theta_k=0.0, n_scan=24, verbose=false) + ballooning_alpha_boundary(plasma_eq; theta_k=0.0, n_scan=24, max_alpha_scale=8.0, tol=1e-3, + psi_window=BALLOONING_SCAN_PSI_WINDOW, verbose=false) Profile driver for the BALOO-style ballooning stability diagram. Loops over flux surfaces returning the experimental pressure gradient `alpha` (from @@ -852,8 +856,9 @@ first crossing, so cost scales with the boundary location rather than the full Use this (not [`ballooning_alpha_boundaries`](@ref)) when only the 1st boundary is needed, e.g. a pedestal-stability constraint. `n_scan` sets the scan resolution; raise it (e.g. 64) on edge surfaces where Δ' poles can shift a coarse first crossing. +`max_alpha_scale` and `tol` bound and resolve the search, as in [`critical_ballooning_alpha`](@ref). -The scan is evaluated only inside `BALLOONING_SCAN_PSI_WINDOW` — ballooning boundaries +The scan is evaluated only inside `psi_window` (default `BALLOONING_SCAN_PSI_WINDOW`) — ballooning boundaries are physically relevant in the mid-radius and pedestal, and the tightly packed axis and far-edge surfaces dominate the scan cost. Per-surface failures, skipped surfaces, and surfaces with no boundary within range are returned as `NaN`. @@ -862,6 +867,9 @@ function ballooning_alpha_boundary( plasma_eq::Equilibrium.PlasmaEquilibrium; theta_k::Float64=0.0, n_scan::Int=24, + max_alpha_scale::Float64=8.0, + tol::Float64=1e-3, + psi_window::Tuple{Float64,Float64}=BALLOONING_SCAN_PSI_WINDOW, verbose::Bool=false ) xs = plasma_eq.profiles.xs @@ -871,10 +879,11 @@ function ballooning_alpha_boundary( alpha_critical = fill(NaN, npsi) Threads.@threads :greedy for i in 1:npsi - _in_ballooning_scan_window(xs[i], xs[end]) || continue + _in_ballooning_scan_window(xs[i], xs[end], psi_window) || continue try alpha[i] = salpha_reference(i, plasma_eq).alpha_ref - alpha_critical[i] = critical_ballooning_alpha(i, plasma_eq; theta_k=theta_k, n_scan=n_scan).alpha_crit + alpha_critical[i] = critical_ballooning_alpha(i, plasma_eq; + theta_k=theta_k, n_scan=n_scan, max_alpha_scale=max_alpha_scale, tol=tol).alpha_crit catch err if verbose @warn "ballooning alpha boundary failed" psi_idx = i exception = err @@ -1140,17 +1149,16 @@ end # Integrates the ideal marginal ballooning equations # ====================================================================== """ - integrate_ballooning_ode(ode_coefficient_spline; theta_k=0.0) + integrate_ballooning_ode(ode_coefficient_spline; theta_k=0.0, tol=1e-8) Integrates the ballooning ODE from `-Infinity` to `-1e-3` and `+Infinity` to `+1e-3`. -Computes Delta' = (y'/y)_right - (y'/y)_left. +Computes Delta' = (y'/y)_right - (y'/y)_left. `tol` is the relative and absolute ODE tolerance. ## Returns - `delta_prime`: The stability index Δ'. """ -function integrate_ballooning_ode(ode_coefficient_spline; theta_k::Float64=0.0) - TOLERANCE = 1e-8 +function integrate_ballooning_ode(ode_coefficient_spline; theta_k::Float64=0.0, tol::Float64=1e-8) MINIMUM_STEP = 1e-10 theta_max = max(ode_coefficient_spline.theta_max, 10.0 * MATCHING_POINT) @@ -1179,8 +1187,8 @@ function integrate_ballooning_ode(ode_coefficient_spline; theta_k::Float64=0.0) sol_left = solve( problem_left, DP5(); - reltol=TOLERANCE, - abstol=TOLERANCE, + reltol=tol, + abstol=tol, dtmin=MINIMUM_STEP, adaptive=true, save_everystep=false, @@ -1200,8 +1208,8 @@ function integrate_ballooning_ode(ode_coefficient_spline; theta_k::Float64=0.0) sol_right = solve( problem_right, DP5(); - reltol=TOLERANCE, - abstol=TOLERANCE, + reltol=tol, + abstol=tol, dtmin=MINIMUM_STEP, adaptive=true, save_everystep=false, diff --git a/src/LocalStability/Problem.jl b/src/LocalStability/Problem.jl index 6ee45fa3d..5c9d97516 100644 --- a/src/LocalStability/Problem.jl +++ b/src/LocalStability/Problem.jl @@ -20,7 +20,8 @@ struct LocalStabilityProblem end """ - Ballooning(; theta_k=0.0, n_scan=24, alpha_boundary_scan=true, verbose=true) + Ballooning(; theta_k=0.0, ode_tol=1e-8, alpha_boundary_scan=true, n_scan=24, max_alpha_scale=8.0, + alpha_tol=1e-3, psi_window=(0.1, 0.99), verbose=true) The local-stability algorithm: the ballooning coefficient system per flux surface, plus the BALOO-style α-boundary scan. Its keywords are the `[LocalStability]` deck keys. @@ -28,14 +29,22 @@ BALOO-style α-boundary scan. Its keywords are the `[LocalStability]` deck keys. ## Fields - `theta_k::Float64` - Ballooning angle θ_k, used by both the profile scan and the α-boundary scan. - - `n_scan::Int` - Resolution of the α scan for the first boundary (see [`ballooning_alpha_boundary`](@ref)). + - `ode_tol::Float64` - Relative and absolute tolerance of the ballooning ODE integration behind Δ'. - `alpha_boundary_scan::Bool` - Run the α-boundary scan, a root find per ψ surface; `false` leaves the α fields empty. + - `n_scan::Int` - Resolution of the α scan for the first boundary (see [`ballooning_alpha_boundary`](@ref)). + - `max_alpha_scale::Float64` - Upper end of the α search, as a multiple of the experimental α; no crossing below it gives NaN. + - `alpha_tol::Float64` - Bisection tolerance of the critical α, in the same multiples. + - `psi_window::Tuple{Float64,Float64}` - ψ_N range of the α scan; surfaces outside it are NaN. - `verbose::Bool` - Print progress and per-surface α-scan failures. """ @kwdef struct Ballooning theta_k::Float64 = 0.0 - n_scan::Int = 24 + ode_tol::Float64 = 1e-8 alpha_boundary_scan::Bool = true + n_scan::Int = 24 + max_alpha_scale::Float64 = 8.0 + alpha_tol::Float64 = 1e-3 + psi_window::Tuple{Float64,Float64} = BALLOONING_SCAN_PSI_WINDOW verbose::Bool = true end @@ -77,9 +86,10 @@ function CommonSolve.solve(prob::LocalStabilityProblem, alg::Ballooning) equil = prob.equil xs = equil.profiles.xs locstab_fs = zeros(Float64, length(xs), 5) - compute_ballooning_stability!(locstab_fs, equil; theta_k=alg.theta_k, verbose=alg.verbose) + compute_ballooning_stability!(locstab_fs, equil; theta_k=alg.theta_k, ode_tol=alg.ode_tol, verbose=alg.verbose) boundary = if alg.alpha_boundary_scan - ballooning_alpha_boundary(equil; theta_k=alg.theta_k, n_scan=alg.n_scan, verbose=alg.verbose) + ballooning_alpha_boundary(equil; theta_k=alg.theta_k, n_scan=alg.n_scan, max_alpha_scale=alg.max_alpha_scale, + tol=alg.alpha_tol, psi_window=alg.psi_window, verbose=alg.verbose) else (psi=Float64[], alpha=Float64[], alpha_critical=Float64[]) end diff --git a/test/runtests_localstability.jl b/test/runtests_localstability.jl index ef937cbe9..c52769b92 100644 --- a/test/runtests_localstability.jl +++ b/test/runtests_localstability.jl @@ -43,6 +43,19 @@ using TOML @test quick.ballooning_Delta_prime == res.ballooning_Delta_prime end + @testset "numerical knobs reach the kernels" begin + narrow = solve(LocalStabilityProblem(equil), Ballooning(; psi_window=(0.3, 0.7), verbose=false)) + inside = [0.3 <= p <= 0.7 for p in xs] + @test all(isfinite, narrow.alpha[inside]) + @test all(isnan, narrow.alpha[.!inside]) + @test narrow.alpha[inside] == res.alpha[inside] + + loose = solve(LocalStabilityProblem(equil), Ballooning(; ode_tol=1e-6, alpha_boundary_scan=false, verbose=false)) + @test loose.D_I == res.D_I + @test loose.ballooning_Delta_prime ≈ res.ballooning_Delta_prime rtol = 1e-3 + @test loose.ballooning_Delta_prime != res.ballooning_Delta_prime + end + @testset "solve writes no file" begin mktempdir() do dir cd(dir) do From 85faf7657b97d272da022dc747ee145295b0ef91 Mon Sep 17 00:00:00 2001 From: Matthew Pharr Date: Fri, 23 Oct 2026 16:03:49 +0200 Subject: [PATCH 3/7] LocalStability - API! - Run LocalStability from its own [LocalStability] deck section --- .pre-commit-config.yaml | 2 +- src/Analysis/ForceFreeStates.jl | 5 +- src/ForceFreeStates/CoreTypes.jl | 2 - src/GeneralizedPerturbedEquilibrium.jl | 98 +++++++++---------- src/HDF5Schema.jl | 1 + src/LocalStability/Ballooning.jl | 13 --- src/LocalStability/LocalStability.jl | 2 +- test/runtests_localstability.jl | 57 +++++++++++ test/runtests_rerun_from_h5.jl | 8 ++ test/runtests_solve_api.jl | 3 +- .../gpec.toml | 5 +- .../gpec.toml | 5 +- .../gpec.toml | 5 +- .../gpec.toml | 5 +- .../gpec.toml | 5 +- .../gpec.toml | 5 +- 16 files changed, 144 insertions(+), 77 deletions(-) diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 5efde01b0..0fd39daad 100644 --- a/.pre-commit-config.yaml +++ b/.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)$ diff --git a/src/Analysis/ForceFreeStates.jl b/src/Analysis/ForceFreeStates.jl index 03462531f..304b9ba54 100644 --- a/src/Analysis/ForceFreeStates.jl +++ b/src/Analysis/ForceFreeStates.jl @@ -377,8 +377,7 @@ Plot the BALOO-style infinite-n ballooning stability diagram: the experimental pressure gradient α (solid) and the first stability boundary α_crit (dashed) versus normalized poloidal flux ψ_N. Surfaces where the experimental α lies above the boundary are ballooning-unstable. Reads `LocalStability/ballooning_psi`, `LocalStability/alpha`, and -`LocalStability/alpha_critical` (populated when ForceFreeStates runs with -`local_stability_flag = true`). +`LocalStability/alpha_critical` (populated when the deck has a `[LocalStability]` section). ### Arguments @@ -400,7 +399,7 @@ function plot_ballooning_alpha_boundary(h5path; save_path=nothing, psi_min=0.0) read(fid["LocalStability/ballooning_psi"]), read(fid["LocalStability/alpha"]), read(fid["LocalStability/alpha_critical"]) end - isempty(alpha) && return plot(; title="No local stability data (set local_stability_flag)", legend=false) + isempty(alpha) && return plot(; title="No local stability data (add a [LocalStability] section)", legend=false) p = plot(; xlims=(psi_min, 1), diff --git a/src/ForceFreeStates/CoreTypes.jl b/src/ForceFreeStates/CoreTypes.jl index b783c6f50..75c05ea8f 100644 --- a/src/ForceFreeStates/CoreTypes.jl +++ b/src/ForceFreeStates/CoreTypes.jl @@ -118,7 +118,6 @@ gpec.toml. ## Fields - `verbose::Bool` - Enable verbose output - - `local_stability_flag::Bool` - Enable local stability analysis (`D_I` and ballooning) - `vac_flag::Bool` - Enable vacuum region calculation - `mtheta_vacuum::Int` - Number of vacuum poloidal grid points (corresponds to `mtheta` in VacuumInput) - `nzeta_vacuum::Int` - Number of vacuum toroidal grid points (corresponds to `nzeta` in VacuumInput3D) @@ -153,7 +152,6 @@ gpec.toml. """ @kwdef struct ForceFreeStatesControl verbose::Bool = true - local_stability_flag::Bool = false vac_flag::Bool = false mtheta_vacuum::Int = 480 nzeta_vacuum::Int = 1 diff --git a/src/GeneralizedPerturbedEquilibrium.jl b/src/GeneralizedPerturbedEquilibrium.jl index e4f91c59b..816fa3781 100755 --- a/src/GeneralizedPerturbedEquilibrium.jl +++ b/src/GeneralizedPerturbedEquilibrium.jl @@ -110,6 +110,24 @@ function _drop_deprecated_keys!(table, deprecated_keys, section::String) return table end +""" + local_stability_algorithm(inputs) -> Union{LocalStability.Ballooning,Nothing} + +The [`Ballooning`](@ref) algorithm a deck's `[LocalStability]` section poses, its keys being the +`Ballooning` keywords, or `nothing` when the deck has no such section. Unknown keys are an error. +""" +function local_stability_algorithm(inputs) + haskey(inputs, "LocalStability") || return nothing + table = inputs["LocalStability"] + allowed = string.(fieldnames(LocalStability.Ballooning)) + for k in keys(table) + k in allowed || throw(ArgumentError("unknown key \"$k\" in [LocalStability] (allowed: $(join(allowed, ", ")))")) + end + kwargs = Dict{Symbol,Any}(Symbol(k) => v for (k, v) in table) + haskey(kwargs, :psi_window) && (kwargs[:psi_window] = Tuple{Float64,Float64}(kwargs[:psi_window])) + return LocalStability.Ballooning(; kwargs...) +end + function main(args::Vector{String}=String[]; dd::Union{IMASdd.dd,Nothing}=nothing) # Every input source builds a ready `(inputs, eq_config, additional_input)` and hands it to # `main_from_inputs`: a gpec.toml working directory, an IMAS `dd`, or a gpec.h5 snapshot. @@ -184,9 +202,10 @@ forcing file it writes from that snapshot. against the current equilibrium) without the original `.dat`/`.h5` files. The coil geometry actually used by the run is always written back into `Input/RawInputs/Coils/`. -Returns `(; ffs, pe, slayer)`: the `ForceFreeStates.ForceFreeStatesResult`, the -`PerturbedEquilibriumState` (`nothing` when that stage did not run) and the SLAYER result -(`nothing` when that stage did not run or failed). An equilibrium-only run +Returns `(; ffs, pe, slayer, local_stability, ...)`: the `ForceFreeStates.ForceFreeStatesResult`, the +`PerturbedEquilibriumState` (`nothing` when that stage did not run), the SLAYER result +(`nothing` when that stage did not run or failed) and the `LocalStabilityResult` (`nothing` +without a `[LocalStability]` section). An equilibrium-only run (`force_termination` in `[Equilibrium]`) returns `nothing` — it never reaches the solve. """ function main_from_inputs( @@ -215,6 +234,7 @@ function main_from_inputs( ffs_table = inputs["ForceFreeStates"] _drop_deprecated_keys!(ffs_table, _DEPRECATED_FFS_KEYS, "ForceFreeStates") ctrl = ForceFreeStatesControl(; (Symbol(k) => v for (k, v) in ffs_table)...) + local_stability_alg = local_stability_algorithm(inputs) resolve_mode_space!(intr, ctrl) @@ -248,13 +268,22 @@ function main_from_inputs( forcing_snapshot = snapshot_forcing_modes(inputs, path, ctrl) + # ---------------------------------------------------------------- + # Local Stability + # ---------------------------------------------------------------- + local_stability = nothing + if local_stability_alg !== nothing + local_stability = _timed_stage(runtimes, "local_stability", "Local Stability") do + solve(LocalStabilityProblem(equil), local_stability_alg) + end + end + # ---------------------------------------------------------------- # Force-Free States # ---------------------------------------------------------------- @info "\n Force-Free States\n$_SECTION" ffs_start = time() - locstab, ballooning_boundary = run_local_stability(ctrl, equil) metric, mats, layer_overlap = prepare_force_free_states!(intr, ctrl, equil, kf_ctrl, kinetic_profiles; species=kf_species, overlap_inputs=_overlap_slayer_inputs(inputs, intr.dir_path)) @@ -266,8 +295,7 @@ function main_from_inputs( git_version=git_version, inputs=inputs, forcing_snapshot=forcing_snapshot, - locstab=locstab, - ballooning_boundary=ballooning_boundary, + local_stability=local_stability, layer_overlap=layer_overlap ) written_h5 = ctrl.HDF5_filename @@ -356,7 +384,7 @@ function main_from_inputs( # TODO: Do not allow perturbed equilibrium calculations if zero crossings are found - return (; ffs=ffs_result, pe=pe_state, slayer=slayer_result, coil_sensitivities, monte_carlo, locking_risk, efc_couplings) + return (; ffs=ffs_result, pe=pe_state, slayer=slayer_result, local_stability, coil_sensitivities, monte_carlo, locking_risk, efc_couplings) end @@ -537,24 +565,6 @@ function snapshot_forcing_modes(inputs::Dict{String,Any}, path::String, ctrl::Fo return modes, normalization end -""" - run_local_stability(ctrl, equil) -> (locstab, ballooning_boundary) - -Run the LocalStability stage when `local_stability_flag` is set. `locstab` holds `D_I` from the -ballooning coefficient system and the local ballooning result, `ballooning_boundary` the first -α-vs-ψ_N stability boundary; both are empty placeholders when the stage is off. -""" -function run_local_stability(ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium) - locstab = nothing - ballooning_boundary = (psi=Float64[], alpha=Float64[], alpha_critical=Float64[]) - if ctrl.local_stability_flag - locstab = LocalStability.compute_local_stability(equil; verbose=ctrl.verbose) - # First ballooning stability boundary (α vs ψ_N) for BALOO-style diagnostics. - ballooning_boundary = LocalStability.ballooning_alpha_boundary(equil; verbose=ctrl.verbose) - end - return locstab, ballooning_boundary -end - # Kinetic profiles for the overlap scan come from the [SLAYER] section, through the same control # and loader the SLAYER analysis itself uses, so the scan sees the same profiles, the same HDF5 # group and the same chi(psi). Returns `nothing` when there is no usable [SLAYER] profile_file. @@ -900,6 +910,8 @@ function EulerLagrangeProblem( kwargs... ) ctrl_kwargs = Dict{Symbol,Any}(kwargs) + haskey(ctrl_kwargs, :local_stability_flag) && + error("local stability is its own stage now; drop local_stability_flag and call solve(LocalStabilityProblem(equil), Ballooning())") (haskey(ctrl_kwargs, :nn_low) || haskey(ctrl_kwargs, :nn_high)) && error("the toroidal mode range comes from the `nn` keyword; drop nn_low/nn_high") ctrl_kwargs[:nn_low] = first(nn) @@ -960,12 +972,11 @@ function solve(prob::EulerLagrangeProblem, alg::ForceFreeStates.AbstractIntegrat equil = maybe_reform_equilibrium(equil, equil.config, nothing, intr, ctrl, nothing) end - locstab, ballooning_boundary = run_local_stability(ctrl, equil) metric, mats, _ = prepare_force_free_states!(intr, ctrl, equil, kf_ctrl, nothing) result = run_force_free_states(ctrl, equil, mats, intr, metric; cut_solution=alg isa ForceFreeStates.Galerkin && alg.cut_solution) if ctrl.write_outputs_to_HDF5 - write_outputs_to_HDF5(result; locstab=locstab, ballooning_boundary=ballooning_boundary) + write_outputs_to_HDF5(result) @info "Results written to $(ctrl.HDF5_filename)" end @@ -1453,7 +1464,7 @@ end """ write_outputs_to_HDF5(result::ForceFreeStatesResult; git_version, inputs, forcing_snapshot, - locstab, ballooning_boundary) + local_stability) Write the HDF5 output file with the run, equilibrium and stability products carried by a force-free-states `result`. This combines the functionality of several pieces of the Fortran @@ -1469,8 +1480,8 @@ instead), and ψ, q, the step counts, `crit`, the asymptotic coefficients and th from `result.diagnostics`. Either may be absent, and is then written empty. `forcing_snapshot` is the `(modes, normalization)` pair from `snapshot_forcing_modes`, -written to `Input/RawInputs/ForcingTerms/`. `locstab` and `ballooning_boundary` come from the -LocalStability stage, which is not part of the force-free-states solve. +written to `Input/RawInputs/ForcingTerms/`. `local_stability` is the `LocalStabilityResult` +of the separate LocalStability stage, or `nothing` when it did not run. ### TODOs @@ -1481,8 +1492,7 @@ function write_outputs_to_HDF5( git_version::String="unknown", inputs::Union{Nothing,Dict{String,Any}}=nothing, forcing_snapshot::Union{Nothing,Tuple{Vector{ForcingTerms.ForcingMode},String}}=nothing, - locstab::Union{FastInterpolations.CubicSeriesInterpolant,Nothing}=nothing, - ballooning_boundary=(psi=Float64[], alpha=Float64[], alpha_critical=Float64[]), + local_stability::Union{LocalStabilityResult,Nothing}=nothing, layer_overlap=nothing ) @@ -1585,25 +1595,13 @@ function write_outputs_to_HDF5( # LocalStability/D_I = Mercier D_I (det(d0bar)); LocalStability/D_R = resistive interchange D_R; # LocalStability/ballooning_Delta_prime = high-n ballooning Δ' (distinct from the Riccati # tearing Δ' under PerturbedEquilibrium/SingularCoupling/Delta_prime). - if locstab !== nothing - locstab_xs = locstab.cache.x - out_h5["LocalStability/psi"] = collect(locstab_xs) # cached vector → dense for HDF5 - out_h5["LocalStability/D_I"] = locstab.y[:, 1] ./ locstab_xs - out_h5["LocalStability/D_R"] = locstab.y[:, 2] ./ locstab_xs - else - out_h5["LocalStability/psi"] = Float64[] - out_h5["LocalStability/D_I"] = Float64[] - out_h5["LocalStability/D_R"] = Float64[] + # The α boundary (ballooning_psi, alpha, alpha_critical) is on its own scan grid, distinct from LocalStability/psi. + ls = local_stability + for name in (:psi, :D_I, :D_R, :ballooning_Delta_prime, :ballooning_psi, :alpha, :alpha_critical) + out_h5["LocalStability/$name"] = ls === nothing ? Float64[] : getfield(ls, name) end - out_h5["SingularSurfaces/D_I"] = (locstab !== nothing && !isempty(result.surfaces)) ? - [locstab(sing.psifac)[1] / sing.psifac for sing in result.surfaces] : Float64[] - out_h5["LocalStability/ballooning_Delta_prime"] = locstab !== nothing ? locstab.y[:, 4] : Float64[] - - # First ballooning stability boundary: experimental α vs critical α (BALOO-style). - # Its own scan grid, distinct from the LocalStability/psi profile grid above. - out_h5["LocalStability/ballooning_psi"] = ballooning_boundary.psi - out_h5["LocalStability/alpha"] = ballooning_boundary.alpha - out_h5["LocalStability/alpha_critical"] = ballooning_boundary.alpha_critical + out_h5["SingularSurfaces/D_I"] = (ls !== nothing && !isempty(result.surfaces)) ? + [LocalStability.mercier_D_I(ls, sing.psifac) for sing in result.surfaces] : Float64[] # Resistive-layer overlap scan. Written whenever the scan ran, whether or not it # constrained the domain, so a run always shows where layer physics would have cut. diff --git a/src/HDF5Schema.jl b/src/HDF5Schema.jl index e32b94b9a..4bfefee26 100644 --- a/src/HDF5Schema.jl +++ b/src/HDF5Schema.jl @@ -281,6 +281,7 @@ const MAIN_H5_ANNOTATIONS = [ # Stages that did not run leave their path absent and are skipped by `annotate!`. const RUNTIME_H5_ANNOTATIONS = [ "Info/Runtimes/equilibrium" => (; long_name="wall-clock time of the Equilibrium construction stage", units="s"), + "Info/Runtimes/local_stability" => (; long_name="wall-clock time of the LocalStability stage (Mercier, resistive interchange and ballooning scans)", units="s"), "Info/Runtimes/force_free_states" => (; long_name="wall-clock time of the ForceFreeStates stability stage", units="s"), "Info/Runtimes/tearing" => (; long_name="wall-clock time of the Tearing stage (SLAYER inner-layer solve)", units="s"), "Info/Runtimes/perturbed_equilibrium" => (; long_name="wall-clock time of the PerturbedEquilibrium stage", units="s"), diff --git a/src/LocalStability/Ballooning.jl b/src/LocalStability/Ballooning.jl index 9979b913b..2c33a2d9a 100644 --- a/src/LocalStability/Ballooning.jl +++ b/src/LocalStability/Ballooning.jl @@ -79,19 +79,6 @@ function compute_ballooning_stability!( end -""" - compute_local_stability(plasma_eq; verbose=false) -> CubicSeriesInterpolant - -Local stability profile spline over `plasma_eq.profiles.xs`, with the columns filled by -[`compute_ballooning_stability!`](@ref): 1 = `D_I·ψ`, 2 = `D_R·ψ`, 4 = ballooning `Δ'`. -""" -function compute_local_stability(plasma_eq::Equilibrium.PlasmaEquilibrium; verbose::Bool=false) - xs = plasma_eq.profiles.xs - locstab_fs = zeros(Float64, length(xs), 5) - compute_ballooning_stability!(locstab_fs, plasma_eq; verbose=verbose) - return cubic_interp(xs, Series(locstab_fs); extrap=ExtendExtrap()) -end - """ resistive_interchange_h(flux_surface_index, plasma_eq) diff --git a/src/LocalStability/LocalStability.jl b/src/LocalStability/LocalStability.jl index e68c8029a..6aa492c38 100644 --- a/src/LocalStability/LocalStability.jl +++ b/src/LocalStability/LocalStability.jl @@ -15,6 +15,6 @@ include("Ballooning.jl") include("Problem.jl") export LocalStabilityProblem, Ballooning, LocalStabilityResult, mercier_D_I -export compute_local_stability, compute_ballooning_stability!, ballooning_alpha_boundary, ballooning_alpha_boundaries +export compute_ballooning_stability!, ballooning_alpha_boundary, ballooning_alpha_boundaries end diff --git a/test/runtests_localstability.jl b/test/runtests_localstability.jl index c52769b92..dd05b59ad 100644 --- a/test/runtests_localstability.jl +++ b/test/runtests_localstability.jl @@ -1,3 +1,4 @@ +using HDF5 using TOML # The LocalStability stage on the coarse Solovev fixture deck (mpsi=16): the named profiles of @@ -64,4 +65,60 @@ using TOML @test isempty(readdir(dir)) end end + + # Copy the fixture deck into `dir` with HDF5 output on, no [LocalStability] section, and the + # given `[ForceFreeStates]` keys and extra sections. + function _deck(dir; ffs_extra=Dict{String,Any}(), sections=Dict{String,Any}()) + for name in readdir(template) + cp(joinpath(template, name), joinpath(dir, name)) + end + toml_path = joinpath(dir, "gpec.toml") + inputs = TOML.parsefile(toml_path) + delete!(inputs, "LocalStability") + inputs["ForceFreeStates"]["write_outputs_to_HDF5"] = true + merge!(inputs["ForceFreeStates"], ffs_extra) + merge!(inputs, sections) + open(io -> TOML.print(io, inputs), toml_path, "w") + return dir + end + _h5_local_stability(dir) = h5open(joinpath(dir, "gpec.h5"), "r") do f + Dict(k => read(f, "LocalStability/$k") for k in keys(f["LocalStability"])) + end + quiet_section = Dict{String,Any}("LocalStability" => Dict{String,Any}("verbose" => false)) + + @testset "deck section and solve give identical results" begin + mktempdir() do dir + run = GPEC.main([_deck(dir; sections=quiet_section)]) + for name in fieldnames(LocalStabilityResult) + name === :di_psi && continue + @test isequal(getfield(run.local_stability, name), getfield(res, name)) + end + h5 = _h5_local_stability(dir) + @test h5["D_I"] == res.D_I + @test isequal(h5["alpha_critical"], res.alpha_critical) + sing_di = h5open(f -> read(f, "SingularSurfaces/D_I"), joinpath(dir, "gpec.h5")) + @test !isempty(sing_di) + @test sing_di == [LS.mercier_D_I(res, s.psifac) for s in run.ffs.surfaces] + end + end + + @testset "no section, no stage" begin + mktempdir() do dir + run = GPEC.main([_deck(dir)]) + @test run.local_stability === nothing + @test all(isempty, values(_h5_local_stability(dir))) + end + end + + @testset "deck errors" begin + mktempdir() do dir + @test_throws MethodError GPEC.main([_deck(dir; ffs_extra=Dict{String,Any}("local_stability_flag" => true))]) + end + @test_throws ArgumentError GPEC.local_stability_algorithm(Dict{String,Any}("LocalStability" => Dict{String,Any}("n_scna" => 32))) + alg = GPEC.local_stability_algorithm(Dict{String,Any}("LocalStability" => Dict{String,Any}("psi_window" => [0.2, 0.8], "n_scan" => 32))) + @test alg.psi_window === (0.2, 0.8) + @test alg.n_scan == 32 + @test GPEC.local_stability_algorithm(Dict{String,Any}()) === nothing + @test_throws ErrorException EulerLagrangeProblem(equil; nn=1, local_stability_flag=true) + end end diff --git a/test/runtests_rerun_from_h5.jl b/test/runtests_rerun_from_h5.jl index 7aea77bab..d3415d396 100644 --- a/test/runtests_rerun_from_h5.jl +++ b/test/runtests_rerun_from_h5.jl @@ -126,6 +126,14 @@ end end end + @testset "the LocalStability section replays" begin + mktempdir() do out_dir + inputs, = GeneralizedPerturbedEquilibrium.build_inputs_from_h5([source_h5, "--output-dir", out_dir]) + @test haskey(inputs, "LocalStability") + @test !isempty(h5open(f -> read(f, "LocalStability/D_I"), source_h5)) + end + end + @testset "rerun refuses to overwrite its own source" begin # Default output is `_rerun.h5`, so clobbering only happens if the user # routes the rerun back onto the source file; exercise that guard. diff --git a/test/runtests_solve_api.jl b/test/runtests_solve_api.jl index 7d8272e5d..a28073ce6 100644 --- a/test/runtests_solve_api.jl +++ b/test/runtests_solve_api.jl @@ -31,6 +31,7 @@ using HDF5 inputs = TOML.parsefile(toml_path) merge!(inputs["Equilibrium"], eq_overrides) merge!(inputs["ForceFreeStates"], ffs_overrides) + delete!(inputs, "LocalStability") merge!(inputs, extra_sections) open(io -> TOML.print(io, inputs), toml_path, "w") return dir @@ -373,7 +374,7 @@ using HDF5 @testset "analytic equilibrium on the two-pass auto grid" begin # A coarse psi_accuracy and no local-stability scan keep the 1024-interval pass-1 grid cheap. eq_overrides = Dict{String,Any}("grid_type" => "auto", "mpsi" => 0, "psi_accuracy" => 3e-2) - ffs_overrides = Dict{String,Any}("local_stability_flag" => false, "vac_flag" => false) + ffs_overrides = Dict{String,Any}("vac_flag" => false) mktempdir() do dir reference = GPEC.main([_stage_deck(dir; eq_overrides=eq_overrides, ffs_overrides=ffs_overrides)]).ffs pass1 = GPEC.Equilibrium.setup_equilibrium( diff --git a/test/test_data/regression_solovev_ideal_example/gpec.toml b/test/test_data/regression_solovev_ideal_example/gpec.toml index 2e1607f7a..5496d83f0 100644 --- a/test/test_data/regression_solovev_ideal_example/gpec.toml +++ b/test/test_data/regression_solovev_ideal_example/gpec.toml @@ -24,8 +24,11 @@ triangularity = 0.5 # Triangularity parameter for wall shape 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] +verbose = false # 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 diff --git a/test/test_data/regression_solovev_ideal_example_multi_n/gpec.toml b/test/test_data/regression_solovev_ideal_example_multi_n/gpec.toml index 9ba608c55..517826893 100644 --- a/test/test_data/regression_solovev_ideal_example_multi_n/gpec.toml +++ b/test/test_data/regression_solovev_ideal_example_multi_n/gpec.toml @@ -24,8 +24,11 @@ triangularity = 0.5 # Triangularity parameter for wall shape 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] +verbose = false # 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 diff --git a/test/test_data/regression_solovev_kinetic_calculated/gpec.toml b/test/test_data/regression_solovev_kinetic_calculated/gpec.toml index a05caa035..c2b0dcaa3 100644 --- a/test/test_data/regression_solovev_kinetic_calculated/gpec.toml +++ b/test/test_data/regression_solovev_kinetic_calculated/gpec.toml @@ -25,8 +25,11 @@ triangularity = 0.5 # Triangularity parameter for wall shape 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] +verbose = false # 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 diff --git a/test/test_data/regression_solovev_kinetic_example/gpec.toml b/test/test_data/regression_solovev_kinetic_example/gpec.toml index 92684ed25..e528e50cb 100644 --- a/test/test_data/regression_solovev_kinetic_example/gpec.toml +++ b/test/test_data/regression_solovev_kinetic_example/gpec.toml @@ -24,8 +24,11 @@ triangularity = 0.5 # Triangularity parameter for wall shape 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] +verbose = false # 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 diff --git a/test/test_data/regression_solovev_kinetic_multi_n/gpec.toml b/test/test_data/regression_solovev_kinetic_multi_n/gpec.toml index e70b6cccc..0d6fde149 100644 --- a/test/test_data/regression_solovev_kinetic_multi_n/gpec.toml +++ b/test/test_data/regression_solovev_kinetic_multi_n/gpec.toml @@ -24,8 +24,11 @@ triangularity = 0.5 # Triangularity parameter for wall shape 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] +verbose = false # 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 diff --git a/test/test_data/regression_solovev_kinetic_nuzero/gpec.toml b/test/test_data/regression_solovev_kinetic_nuzero/gpec.toml index 5df1debce..7fd01216f 100644 --- a/test/test_data/regression_solovev_kinetic_nuzero/gpec.toml +++ b/test/test_data/regression_solovev_kinetic_nuzero/gpec.toml @@ -25,8 +25,11 @@ triangularity = 0.5 # Triangularity parameter for wall shape 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] +verbose = false # 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 From 50b515adf70fd61ab4dc26f5c74d1a4764ca77d2 Mon Sep 17 00:00:00 2001 From: Matthew Pharr Date: Fri, 23 Oct 2026 16:06:35 +0200 Subject: [PATCH 4/7] Examples - MINOR - Move every deck's local_stability_flag to a [LocalStability] section --- benchmarks/benchmark_against_fortran_run.jl | 6 ++++-- benchmarks/benchmark_diiid_ideal_ntv_torque.jl | 3 ++- benchmarks/benchmark_diiid_kinetic_stability.jl | 3 ++- docs/development/architecture.md | 3 ++- docs/src/stability.md | 11 ++++++----- examples/DIIID-like_SLAYER_example/gpec.toml | 12 +++++++++++- examples/DIIID-like_error_field_example/gpec.toml | 12 +++++++++++- examples/DIIID-like_gal_resistive_example/gpec.toml | 12 +++++++++++- .../DIIID-like_gal_resistive_pe_example/gpec.toml | 12 +++++++++++- examples/DIIID-like_ideal_example/gpec.toml | 12 +++++++++++- examples/DIIID-like_ideal_example_IMAS/gpec.toml | 1 - .../run_imas_example.jl | 1 - examples/DIIID-like_ideal_multi_n_example/gpec.toml | 12 +++++++++++- .../DIIID-like_riccati_deltaprime_example/gpec.toml | 12 +++++++++++- examples/DIIID-like_truncation_example/gpec.toml | 1 - examples/LAR_beta_scan/gpec.toml | 12 +++++++++++- examples/LAR_epsilon_scan/gpec.toml | 12 +++++++++++- examples/LAR_ideal_match_test/gpec.toml | 12 +++++++++++- examples/LAR_resistive_match_test/gpec.toml | 12 +++++++++++- examples/Solovev_ideal_example/gpec.toml | 12 +++++++++++- examples/Solovev_ideal_example_3D/gpec.toml | 12 +++++++++++- examples/Solovev_ideal_example_multi_n/gpec.toml | 12 +++++++++++- .../single_n_1/gpec.toml | 12 +++++++++++- .../single_n_2/gpec.toml | 12 +++++++++++- examples/Solovev_kinetic_NTV_example/gpec.toml | 12 +++++++++++- .../Solovev_kinetic_calculated_example/gpec.toml | 12 +++++++++++- examples/Solovev_kinetic_multiion_example/gpec.toml | 12 +++++++++++- examples/a10_kinetic_example/gpec.toml | 12 +++++++++++- regression-harness/cases/diiid_n1.toml | 2 +- .../regression_solovev_ideal_example/gpec.toml | 7 +++++++ .../gpec.toml | 7 +++++++ .../regression_solovev_kinetic_calculated/gpec.toml | 7 +++++++ .../regression_solovev_kinetic_example/gpec.toml | 7 +++++++ .../regression_solovev_kinetic_multi_n/gpec.toml | 7 +++++++ .../regression_solovev_kinetic_nuzero/gpec.toml | 7 +++++++ 35 files changed, 279 insertions(+), 34 deletions(-) diff --git a/benchmarks/benchmark_against_fortran_run.jl b/benchmarks/benchmark_against_fortran_run.jl index 841995301..7adfffc22 100644 --- a/benchmarks/benchmark_against_fortran_run.jl +++ b/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) diff --git a/benchmarks/benchmark_diiid_ideal_ntv_torque.jl b/benchmarks/benchmark_diiid_ideal_ntv_torque.jl index 430c353cc..ade40501f 100644 --- a/benchmarks/benchmark_diiid_ideal_ntv_torque.jl +++ b/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 diff --git a/benchmarks/benchmark_diiid_kinetic_stability.jl b/benchmarks/benchmark_diiid_kinetic_stability.jl index 1c87d9814..21db14304 100644 --- a/benchmarks/benchmark_diiid_kinetic_stability.jl +++ b/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 diff --git a/docs/development/architecture.md b/docs/development/architecture.md index 059539a8b..ba99748f5 100644 --- a/docs/development/architecture.md +++ b/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 diff --git a/docs/src/stability.md b/docs/src/stability.md index b056753ab..8d845b5db 100644 --- a/docs/src/stability.md +++ b/docs/src/stability.md @@ -156,8 +156,12 @@ 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)). The section replaces the removed `local_stability_flag` +of `[ForceFreeStates]`. 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 @@ -379,9 +383,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 diff --git a/examples/DIIID-like_SLAYER_example/gpec.toml b/examples/DIIID-like_SLAYER_example/gpec.toml index 3732f8076..c32ca67c4 100644 --- a/examples/DIIID-like_SLAYER_example/gpec.toml +++ b/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 diff --git a/examples/DIIID-like_error_field_example/gpec.toml b/examples/DIIID-like_error_field_example/gpec.toml index 2abb9ef53..958067a38 100644 --- a/examples/DIIID-like_error_field_example/gpec.toml +++ b/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 diff --git a/examples/DIIID-like_gal_resistive_example/gpec.toml b/examples/DIIID-like_gal_resistive_example/gpec.toml index dac772545..d67b7bcb2 100644 --- a/examples/DIIID-like_gal_resistive_example/gpec.toml +++ b/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 diff --git a/examples/DIIID-like_gal_resistive_pe_example/gpec.toml b/examples/DIIID-like_gal_resistive_pe_example/gpec.toml index 5398fd0d6..f418ccf57 100644 --- a/examples/DIIID-like_gal_resistive_pe_example/gpec.toml +++ b/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 diff --git a/examples/DIIID-like_ideal_example/gpec.toml b/examples/DIIID-like_ideal_example/gpec.toml index 261667ba0..78bdaf445 100644 --- a/examples/DIIID-like_ideal_example/gpec.toml +++ b/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 diff --git a/examples/DIIID-like_ideal_example_IMAS/gpec.toml b/examples/DIIID-like_ideal_example_IMAS/gpec.toml index 0893460a2..b1e8c9119 100644 --- a/examples/DIIID-like_ideal_example_IMAS/gpec.toml +++ b/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 diff --git a/examples/DIIID-like_ideal_example_IMAS/run_imas_example.jl b/examples/DIIID-like_ideal_example_IMAS/run_imas_example.jl index c4fe1637f..ef01d67eb 100644 --- a/examples/DIIID-like_ideal_example_IMAS/run_imas_example.jl +++ b/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, diff --git a/examples/DIIID-like_ideal_multi_n_example/gpec.toml b/examples/DIIID-like_ideal_multi_n_example/gpec.toml index 325a24854..6b6c8b541 100644 --- a/examples/DIIID-like_ideal_multi_n_example/gpec.toml +++ b/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 diff --git a/examples/DIIID-like_riccati_deltaprime_example/gpec.toml b/examples/DIIID-like_riccati_deltaprime_example/gpec.toml index 2c369b253..8ed214bd9 100644 --- a/examples/DIIID-like_riccati_deltaprime_example/gpec.toml +++ b/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) diff --git a/examples/DIIID-like_truncation_example/gpec.toml b/examples/DIIID-like_truncation_example/gpec.toml index eb2180e03..341c01263 100644 --- a/examples/DIIID-like_truncation_example/gpec.toml +++ b/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) diff --git a/examples/LAR_beta_scan/gpec.toml b/examples/LAR_beta_scan/gpec.toml index da0e40fe1..bfd5ed0e5 100644 --- a/examples/LAR_beta_scan/gpec.toml +++ b/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) diff --git a/examples/LAR_epsilon_scan/gpec.toml b/examples/LAR_epsilon_scan/gpec.toml index ee6e3709c..e7ace3ad5 100644 --- a/examples/LAR_epsilon_scan/gpec.toml +++ b/examples/LAR_epsilon_scan/gpec.toml @@ -34,8 +34,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) diff --git a/examples/LAR_ideal_match_test/gpec.toml b/examples/LAR_ideal_match_test/gpec.toml index b24604a11..7f4378750 100644 --- a/examples/LAR_ideal_match_test/gpec.toml +++ b/examples/LAR_ideal_match_test/gpec.toml @@ -30,9 +30,19 @@ mtheta = 128 # Poloidal grid points shape = "conformal" # Wall shape (nowall, conformal, elliptical, dee, mod_dee, filepath) standoff_fraction = 20 # Wall standoff in plasma minor radii (conformal, dee) +# 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] HDF5_filename = "gpec_lar_ideal.h5" # Output filename for the stability-stage results -local_stability_flag = true # Perform local stability analysis (Mercier and ballooning) across the psi 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) diff --git a/examples/LAR_resistive_match_test/gpec.toml b/examples/LAR_resistive_match_test/gpec.toml index 1e2bb7ed6..23bf53ba5 100644 --- a/examples/LAR_resistive_match_test/gpec.toml +++ b/examples/LAR_resistive_match_test/gpec.toml @@ -31,9 +31,19 @@ mtheta = 128 # Poloidal grid points shape = "conformal" # Wall shape (nowall, conformal, elliptical, dee, mod_dee, filepath) standoff_fraction = 20 # Wall standoff in plasma minor radii (conformal, dee) +# 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] HDF5_filename = "gpec_lar_resistive.h5" # Output filename for the stability-stage results -local_stability_flag = true # Perform local stability analysis (Mercier and ballooning) across the psi 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) diff --git a/examples/Solovev_ideal_example/gpec.toml b/examples/Solovev_ideal_example/gpec.toml index a42fe4cb0..bed91d59c 100644 --- a/examples/Solovev_ideal_example/gpec.toml +++ b/examples/Solovev_ideal_example/gpec.toml @@ -49,8 +49,18 @@ regularization_width = 0.05 # Regularization width for singular surf # The SLAYER example and regression case use the DIII-D-like geqdsk instead # (see examples/DIIID-like_SLAYER_example). +# 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 diff --git a/examples/Solovev_ideal_example_3D/gpec.toml b/examples/Solovev_ideal_example_3D/gpec.toml index ab1909e15..fa8d52730 100644 --- a/examples/Solovev_ideal_example_3D/gpec.toml +++ b/examples/Solovev_ideal_example_3D/gpec.toml @@ -14,8 +14,18 @@ newq0 = 0 # Target on-axis safety factor q(0) (0 = flux_surface_rtol = 1e-7 # Relative tolerance of the ODE solves that trace each flux surface force_termination = false # Terminate after equilibrium setup (skip stability calculations) +# 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 = 1.0 # Edge dW(ψ) diagnostic scan band [dW_edge_scan_start, psilim]; set ≥ psilim to disable diff --git a/examples/Solovev_ideal_example_multi_n/gpec.toml b/examples/Solovev_ideal_example_multi_n/gpec.toml index 0150e1b73..949579bf6 100644 --- a/examples/Solovev_ideal_example_multi_n/gpec.toml +++ b/examples/Solovev_ideal_example_multi_n/gpec.toml @@ -26,8 +26,18 @@ triangularity = 0.5 # Triangularity parameter for wall shape 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 diff --git a/examples/Solovev_ideal_example_multi_n/single_n_1/gpec.toml b/examples/Solovev_ideal_example_multi_n/single_n_1/gpec.toml index 472546634..19a2cfeb7 100644 --- a/examples/Solovev_ideal_example_multi_n/single_n_1/gpec.toml +++ b/examples/Solovev_ideal_example_multi_n/single_n_1/gpec.toml @@ -25,8 +25,18 @@ triangularity = 0.5 # Triangularity parameter for wall shape 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 diff --git a/examples/Solovev_ideal_example_multi_n/single_n_2/gpec.toml b/examples/Solovev_ideal_example_multi_n/single_n_2/gpec.toml index 97f108a5a..c8a2708d0 100644 --- a/examples/Solovev_ideal_example_multi_n/single_n_2/gpec.toml +++ b/examples/Solovev_ideal_example_multi_n/single_n_2/gpec.toml @@ -25,8 +25,18 @@ triangularity = 0.5 # Triangularity parameter for wall shape 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 diff --git a/examples/Solovev_kinetic_NTV_example/gpec.toml b/examples/Solovev_kinetic_NTV_example/gpec.toml index 8ff181558..eb30f6c2d 100644 --- a/examples/Solovev_kinetic_NTV_example/gpec.toml +++ b/examples/Solovev_kinetic_NTV_example/gpec.toml @@ -46,8 +46,18 @@ verbose = true # Enable verbose logging write_outputs_to_HDF5 = true # Write perturbed equilibrium outputs to HDF5 regularization_width = 0.05 # Regularization width for singular surfaces (0 = disabled) +# 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 diff --git a/examples/Solovev_kinetic_calculated_example/gpec.toml b/examples/Solovev_kinetic_calculated_example/gpec.toml index c663ab844..8c3524319 100644 --- a/examples/Solovev_kinetic_calculated_example/gpec.toml +++ b/examples/Solovev_kinetic_calculated_example/gpec.toml @@ -43,8 +43,18 @@ verbose = false # Enable verbose logging write_outputs_to_HDF5 = true # Write perturbed equilibrium outputs to HDF5 regularization_width = 0.05 # Regularization width for singular surfaces (0 = disabled) +# 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 = false # 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 diff --git a/examples/Solovev_kinetic_multiion_example/gpec.toml b/examples/Solovev_kinetic_multiion_example/gpec.toml index bee59220e..37e9cc924 100644 --- a/examples/Solovev_kinetic_multiion_example/gpec.toml +++ b/examples/Solovev_kinetic_multiion_example/gpec.toml @@ -46,9 +46,19 @@ verbose = true # Enable verbose logging write_outputs_to_HDF5 = true # Write perturbed equilibrium outputs to HDF5 regularization_width = 0.05 # Regularization width for singular surfaces (0 = disabled) +# 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] integrator = "forward" # Dense ξ profiles required by PerturbedEquilibrium and kinetic torque diagnostics -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 diff --git a/examples/a10_kinetic_example/gpec.toml b/examples/a10_kinetic_example/gpec.toml index b56ff6bbf..fe8544c9a 100644 --- a/examples/a10_kinetic_example/gpec.toml +++ b/examples/a10_kinetic_example/gpec.toml @@ -21,8 +21,18 @@ flux_surface_rtol = 1e-7 # Relative tolerance of the ODE solves that trace [Wall] shape = "nowall" # Wall shape (nowall, conformal, elliptical, dee, mod_dee, filepath) +# 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 truncate_at_rational_offset = false # FALSE for limited/analytical equilibria — rationals sparse, rational_offset_fraction would chop too much edge diff --git a/regression-harness/cases/diiid_n1.toml b/regression-harness/cases/diiid_n1.toml index 1e033e5ca..72ccadc89 100644 --- a/regression-harness/cases/diiid_n1.toml +++ b/regression-harness/cases/diiid_n1.toml @@ -396,7 +396,7 @@ label = "pressure profile (checksum)" noise_threshold = 0 order = 71 -# Local stability profiles (local_stability_flag = true). Checksums hash the full +# Local stability profiles (the [LocalStability] section). Checksums hash the full # profile so they are robust to NaN edge-surface entries and catch any change. [quantities.locstab_di] h5path = "LocalStability/D_I" diff --git a/test/test_data/regression_solovev_ideal_example/gpec.toml b/test/test_data/regression_solovev_ideal_example/gpec.toml index 5496d83f0..2dc2c057c 100644 --- a/test/test_data/regression_solovev_ideal_example/gpec.toml +++ b/test/test_data/regression_solovev_ideal_example/gpec.toml @@ -26,6 +26,13 @@ equal_arc_wall = true # Equal arc length distribution of nodes # 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 = false # Print progress and per-surface α-scan failures [ForceFreeStates] diff --git a/test/test_data/regression_solovev_ideal_example_multi_n/gpec.toml b/test/test_data/regression_solovev_ideal_example_multi_n/gpec.toml index 517826893..f87c07b1b 100644 --- a/test/test_data/regression_solovev_ideal_example_multi_n/gpec.toml +++ b/test/test_data/regression_solovev_ideal_example_multi_n/gpec.toml @@ -26,6 +26,13 @@ equal_arc_wall = true # Equal arc length distribution of nodes # 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 = false # Print progress and per-surface α-scan failures [ForceFreeStates] diff --git a/test/test_data/regression_solovev_kinetic_calculated/gpec.toml b/test/test_data/regression_solovev_kinetic_calculated/gpec.toml index c2b0dcaa3..30fd5f854 100644 --- a/test/test_data/regression_solovev_kinetic_calculated/gpec.toml +++ b/test/test_data/regression_solovev_kinetic_calculated/gpec.toml @@ -27,6 +27,13 @@ equal_arc_wall = true # Equal arc length distribution of nodes # 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 = false # Print progress and per-surface α-scan failures [ForceFreeStates] diff --git a/test/test_data/regression_solovev_kinetic_example/gpec.toml b/test/test_data/regression_solovev_kinetic_example/gpec.toml index e528e50cb..af9bf78f8 100644 --- a/test/test_data/regression_solovev_kinetic_example/gpec.toml +++ b/test/test_data/regression_solovev_kinetic_example/gpec.toml @@ -26,6 +26,13 @@ equal_arc_wall = true # Equal arc length distribution of nodes # 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 = false # Print progress and per-surface α-scan failures [ForceFreeStates] diff --git a/test/test_data/regression_solovev_kinetic_multi_n/gpec.toml b/test/test_data/regression_solovev_kinetic_multi_n/gpec.toml index 0d6fde149..bd31bef3e 100644 --- a/test/test_data/regression_solovev_kinetic_multi_n/gpec.toml +++ b/test/test_data/regression_solovev_kinetic_multi_n/gpec.toml @@ -26,6 +26,13 @@ equal_arc_wall = true # Equal arc length distribution of nodes # 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 = false # Print progress and per-surface α-scan failures [ForceFreeStates] diff --git a/test/test_data/regression_solovev_kinetic_nuzero/gpec.toml b/test/test_data/regression_solovev_kinetic_nuzero/gpec.toml index 7fd01216f..0404c6025 100644 --- a/test/test_data/regression_solovev_kinetic_nuzero/gpec.toml +++ b/test/test_data/regression_solovev_kinetic_nuzero/gpec.toml @@ -27,6 +27,13 @@ equal_arc_wall = true # Equal arc length distribution of nodes # 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 = false # Print progress and per-surface α-scan failures [ForceFreeStates] From 871ffff93ce6c5f8ef0223ffdc7f8a49bfa38ead Mon Sep 17 00:00:00 2001 From: Matthew Pharr Date: Fri, 23 Oct 2026 16:25:51 +0200 Subject: [PATCH 5/7] LocalStability - MINOR - Drop the local_stability_flag checks and let the old key fail as an unknown keyword No backward-compatibility code: remove the EulerLagrangeProblem error for local_stability_flag, its deprecated-key pre-commit entry, the two tests asserting it errors, and the docs sentence naming it. Co-Authored-By: Claude Opus 5.5 --- .pre-commit-config.yaml | 2 +- docs/src/stability.md | 3 +-- src/GeneralizedPerturbedEquilibrium.jl | 2 -- test/runtests_localstability.jl | 6 +----- 4 files changed, 3 insertions(+), 10 deletions(-) diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 0fd39daad..5efde01b0 100644 --- a/.pre-commit-config.yaml +++ b/.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|local_stability_flag)\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)\s*=' files: ^(examples/.*\.toml|test/test_data/.*\.toml)$ diff --git a/docs/src/stability.md b/docs/src/stability.md index 8d845b5db..3f5701379 100644 --- a/docs/src/stability.md +++ b/docs/src/stability.md @@ -160,8 +160,7 @@ A `[LocalStability]` section in the deck runs a local high-``n`` stability scan 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)). The section replaces the removed `local_stability_flag` -of `[ForceFreeStates]`. +(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 diff --git a/src/GeneralizedPerturbedEquilibrium.jl b/src/GeneralizedPerturbedEquilibrium.jl index 816fa3781..ab02e8e7c 100755 --- a/src/GeneralizedPerturbedEquilibrium.jl +++ b/src/GeneralizedPerturbedEquilibrium.jl @@ -910,8 +910,6 @@ function EulerLagrangeProblem( kwargs... ) ctrl_kwargs = Dict{Symbol,Any}(kwargs) - haskey(ctrl_kwargs, :local_stability_flag) && - error("local stability is its own stage now; drop local_stability_flag and call solve(LocalStabilityProblem(equil), Ballooning())") (haskey(ctrl_kwargs, :nn_low) || haskey(ctrl_kwargs, :nn_high)) && error("the toroidal mode range comes from the `nn` keyword; drop nn_low/nn_high") ctrl_kwargs[:nn_low] = first(nn) diff --git a/test/runtests_localstability.jl b/test/runtests_localstability.jl index dd05b59ad..2d2866bbb 100644 --- a/test/runtests_localstability.jl +++ b/test/runtests_localstability.jl @@ -110,15 +110,11 @@ using TOML end end - @testset "deck errors" begin - mktempdir() do dir - @test_throws MethodError GPEC.main([_deck(dir; ffs_extra=Dict{String,Any}("local_stability_flag" => true))]) - end + @testset "deck section parsing" begin @test_throws ArgumentError GPEC.local_stability_algorithm(Dict{String,Any}("LocalStability" => Dict{String,Any}("n_scna" => 32))) alg = GPEC.local_stability_algorithm(Dict{String,Any}("LocalStability" => Dict{String,Any}("psi_window" => [0.2, 0.8], "n_scan" => 32))) @test alg.psi_window === (0.2, 0.8) @test alg.n_scan == 32 @test GPEC.local_stability_algorithm(Dict{String,Any}()) === nothing - @test_throws ErrorException EulerLagrangeProblem(equil; nn=1, local_stability_flag=true) end end From 03cba5d50941fd3451a06ebcbaa24d5580595f43 Mon Sep 17 00:00:00 2001 From: logan-nc Date: Sun, 25 Oct 2026 20:38:31 -0400 Subject: [PATCH 6/7] Docs - DOCS - List the alpha-boundary datasets under the local stability stage Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_0186CLfjeJeSusNTva5RMfGU --- docs/src/stability.md | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/docs/src/stability.md b/docs/src/stability.md index 3f5701379..39576f927 100644 --- a/docs/src/stability.md +++ b/docs/src/stability.md @@ -179,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 From e018bce14b4c632184475c450d8f30ccdc37b79c Mon Sep 17 00:00:00 2001 From: Matthew Pharr Date: Mon, 26 Oct 2026 11:03:51 +0100 Subject: [PATCH 7/7] ForceFreeStates - MINOR - List local_stability_flag as a deprecated key --- .pre-commit-config.yaml | 2 +- src/GeneralizedPerturbedEquilibrium.jl | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 5efde01b0..0fd39daad 100644 --- a/.pre-commit-config.yaml +++ b/.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)$ diff --git a/src/GeneralizedPerturbedEquilibrium.jl b/src/GeneralizedPerturbedEquilibrium.jl index ab02e8e7c..5ed395dec 100755 --- a/src/GeneralizedPerturbedEquilibrium.jl +++ b/src/GeneralizedPerturbedEquilibrium.jl @@ -95,7 +95,7 @@ using .ForcingTerms: RMPField const _DEPRECATED_FFS_KEYS = ("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") + "gal_flag", "nstep", "diagnose_ca", "diagnose", "ksing", "local_stability_flag") const _DEPRECATED_EQUIL_KEYS = ("power_bp", "power_b", "power_r", "power_rc", "use_galgrid") # Drop deprecated keys from a parsed gpec.toml section so legacy files keep parsing