From f8996d4f0d9a7261bc48a9fe2389028e703009dd Mon Sep 17 00:00:00 2001 From: Matthew Pharr Date: Fri, 14 Aug 2026 23:08:10 -0400 Subject: [PATCH] FFS - REFACTOR - Extract local stability into standalone LocalStability module --- REFACTOR_PLAN.md | 6 +- docs/development/architecture.md | 14 ++-- docs/src/ballooning.md | 6 ++ docs/src/stability.md | 2 +- .../analyze_example.jl | 15 ++-- src/Analysis/ForceFreeStates.jl | 4 +- src/ForceFreeStates/ForceFreeStates.jl | 1 - src/GeneralizedPerturbedEquilibrium.jl | 10 ++- .../Ballooning.jl | 70 +++++++++---------- src/LocalStability/LocalStability.jl | 17 +++++ test/runtests.jl | 1 + test/runtests_resist_eval.jl | 3 +- 12 files changed, 92 insertions(+), 57 deletions(-) rename src/{ForceFreeStates => LocalStability}/Ballooning.jl (96%) create mode 100644 src/LocalStability/LocalStability.jl diff --git a/REFACTOR_PLAN.md b/REFACTOR_PLAN.md index 0f5bdc00b..de545c427 100644 --- a/REFACTOR_PLAN.md +++ b/REFACTOR_PLAN.md @@ -670,7 +670,11 @@ manual smoke: run the 4-line UX from the Context section in a REPL against from the §3 spec, both improvements: the new Δ′ example references the DIIID geqdsk by relative path instead of copying it, and the TOML sweep covered six regression fixtures (two more had landed on develop since the plan was written), all `forward`. -- [ ] PR 2 — `refactor/local-stability-module` +- [ ] PR 2 — `refactor/local-stability-module` — **implemented, in review.** One delta + from the §4 spec: the signature change also required updating two call-site groups the + section did not list — `examples/DIIID-like_ideal_example/analyze_example.jl` (five + ballooning entry points) and two docstring cross-references in + `src/Analysis/ForceFreeStates.jl`. - [ ] PR 3 — `refactor/forcefreestates-result` - [ ] PR 4 — `refactor/staged-main` - [ ] PR 5 — `feature/solve-api` diff --git a/docs/development/architecture.md b/docs/development/architecture.md index ef756c9c0..a6ddef8d9 100644 --- a/docs/development/architecture.md +++ b/docs/development/architecture.md @@ -12,7 +12,7 @@ This workflow is reflected in the modular structure and data flow. ## Module Structure -GPEC consists of **seven main modules** organized in `src/`: +GPEC consists of **eight main modules** organized in `src/`: ### Foundation Modules @@ -67,18 +67,23 @@ GPEC consists of **seven main modules** organized in `src/`: - `Fourfit.jl` - Fourier fitting routines - `FixedBoundaryStability.jl` - Fixed boundary analysis - `Free.jl` - Free boundary stability - - `Ballooning.jl` - Local stability scan: Mercier D_I, resistive interchange D_R, and high-n ballooning Δ' (s–α). Replaces the former standalone `Mercier.jl`. - Status: Stable, core DCON functionality implemented +6. **LocalStability** (`src/LocalStability/`) - Local high-n stability + - `Ballooning.jl` - Local stability scan: Mercier D_I, resistive interchange D_R, and high-n ballooning Δ' (s–α). Replaces the former standalone `Mercier.jl`. + - Depends only on Equilibrium (plus math libraries); carries no stability-solver state + - Main entry points: `compute_local_stability`, `ballooning_alpha_boundary` + - Status: Stable + ### Perturbed Equilibrium Modules -6. **ForcingTerms** (`src/ForcingTerms/`) - External field specification +7. **ForcingTerms** (`src/ForcingTerms/`) - External field specification - Handles external magnetic field perturbations (coils, RMP, etc.) - Supports ASCII and HDF5 forcing data formats - `ForcingMode` data structure specifies amplitude and phase for each (m,n) component - Status: Complete and functional -7. **PerturbedEquilibrium** (`src/PerturbedEquilibrium/`) - **GPEC-style plasma response** +8. **PerturbedEquilibrium** (`src/PerturbedEquilibrium/`) - **GPEC-style plasma response** - Computes plasma response to external forcing - Calculates singular coupling metrics at rational surfaces - Key files: @@ -191,6 +196,7 @@ GeneralizedPerturbedEquilibrium ├── Utilities (shared tools) │ └── FourierTransforms ├── Equilibrium (uses Splines) +├── LocalStability (uses Equilibrium) ├── Vacuum (uses Splines, Equilibrium, Utilities) ├── ForcingTerms (data I/O) ├── ForceFreeStates (uses Equilibrium, Vacuum, Splines) diff --git a/docs/src/ballooning.md b/docs/src/ballooning.md index 72697da20..8ef6cd10b 100644 --- a/docs/src/ballooning.md +++ b/docs/src/ballooning.md @@ -911,3 +911,9 @@ The local-stability output now stores ballooning ``\Delta'`` in the fourth `locstab_fs` entry. In the HDF5 output this is written as `LocalStability/ballooning_Delta_prime`, distinct from the tearing ``\Delta'`` outputs under `SingularSurfaces/` and `PerturbedEquilibrium/SingularCoupling/`. + +## API + +```@autodocs +Modules = [GeneralizedPerturbedEquilibrium.LocalStability] +``` diff --git a/docs/src/stability.md b/docs/src/stability.md index 0e8b72880..dadf13762 100644 --- a/docs/src/stability.md +++ b/docs/src/stability.md @@ -270,7 +270,7 @@ The Galerkin Δ′ solver (`src/ForceFreeStates/Galerkin/`) is documented separa ```@autodocs Modules = [GeneralizedPerturbedEquilibrium.ForceFreeStates] -Pages = ["ForceFreeStates.jl", "ForceFreeStatesStructs.jl", "Ballooning.jl", "Resist.jl", "EulerLagrange.jl", "Sing.jl", "Fourfit.jl", "Kinetic.jl", "FixedBoundaryStability.jl", "Utils.jl", "Free.jl", "Riccati.jl"] +Pages = ["ForceFreeStates.jl", "ForceFreeStatesStructs.jl", "Resist.jl", "EulerLagrange.jl", "Sing.jl", "Fourfit.jl", "Kinetic.jl", "FixedBoundaryStability.jl", "Utils.jl", "Free.jl", "Riccati.jl"] ``` ## Example usage diff --git a/examples/DIIID-like_ideal_example/analyze_example.jl b/examples/DIIID-like_ideal_example/analyze_example.jl index 319a7063d..715e3ec76 100644 --- a/examples/DIIID-like_ideal_example/analyze_example.jl +++ b/examples/DIIID-like_ideal_example/analyze_example.jl @@ -2,7 +2,7 @@ using Pkg; Pkg.activate(joinpath(@__DIR__, "../..")) using GeneralizedPerturbedEquilibrium, Plots using GeneralizedPerturbedEquilibrium: Analysis -using GeneralizedPerturbedEquilibrium: Equilibrium, ForceFreeStates +using GeneralizedPerturbedEquilibrium: Equilibrium, ForceFreeStates, LocalStability using LaTeXStrings, TOML isinteractive() ? plotlyjs() : gr() @@ -29,7 +29,7 @@ Plots.savefig(p_baloo, baloo_path) # ---------------------------------------------------------------------- # Local stability: s-alpha profiles, D_I / ballooning Δ', and 2D s-alpha # maps (ported from the former Bal_salpha_delta_di_summary notebook). -# These recompute directly from the equilibrium via the Ballooning.jl helpers +# These recompute directly from the equilibrium via the LocalStability helpers # rather than reading gpec.h5, so the s-alpha scan can perturb (p', q'). # ---------------------------------------------------------------------- @@ -66,7 +66,7 @@ for geqdsk_file in geqdsk_cases alpha_profile = fill(NaN, length(psi_norm)) for i in eachindex(psi_norm) try - ref = ForceFreeStates.salpha_reference(i, equil) + ref = LocalStability.salpha_reference(i, equil) s_profile[i] = ref.s_ref alpha_profile[i] = ref.alpha_ref catch err @@ -107,8 +107,7 @@ for geqdsk_file in geqdsk_cases # D_I (Mercier) and ballooning Δ' profiles locstab_fs = zeros(length(psi_norm), 5) - ctrl = ForceFreeStates.ForceFreeStatesControl(; verbose=false) - ForceFreeStates.compute_ballooning_stability!(ctrl, locstab_fs, equil; compute_delta_prime=true) + LocalStability.compute_ballooning_stability!(locstab_fs, equil; compute_delta_prime=true) delta_prime = Vector(locstab_fs[:, 4]) di_profile = fill(NaN, length(psi_norm)) @@ -152,7 +151,7 @@ for geqdsk_file in geqdsk_cases psi_idx_scan = argmin(abs.(psi_norm .- psi_target)) s_scales = collect(range(-5.0, 5.0; length=30)) alpha_scales = collect(range(-5.0, 5.0; length=30)) - scan = ForceFreeStates.scan_delta_prime_map(psi_idx_scan, equil; theta_k=0.0, s_scales=s_scales, alpha_scales=alpha_scales) + scan = LocalStability.scan_delta_prime_map(psi_idx_scan, equil; theta_k=0.0, s_scales=s_scales, alpha_scales=alpha_scales) scan_qprime = scan.dqdpsi_ref .* scan.s_scales scan_pprime = scan.pprime_ref .* scan.alpha_scales @@ -210,8 +209,8 @@ for geqdsk_file in geqdsk_cases Plots.savefig(p_zero, zero_path) # BALOO-style α vs ψ stability boundaries (1st and 2nd) over the Δ'(ψ, α) map - bnd = ForceFreeStates.ballooning_alpha_boundaries(ctrl, equil) - dpmap = ForceFreeStates.ballooning_delta_prime_map(ctrl, equil) + bnd = LocalStability.ballooning_alpha_boundaries(equil) + dpmap = LocalStability.ballooning_delta_prime_map(equil) p_alpha_bnd = Analysis.ForceFreeStates.plot_ballooning_alpha_boundaries(bnd, dpmap) display(p_alpha_bnd) alpha_bnd_path = joinpath(@__DIR__, "alpha_stability_boundary_$(case_label).png") diff --git a/src/Analysis/ForceFreeStates.jl b/src/Analysis/ForceFreeStates.jl index c92d6e698..e4ce95103 100644 --- a/src/Analysis/ForceFreeStates.jl +++ b/src/Analysis/ForceFreeStates.jl @@ -584,7 +584,7 @@ arrays leave natural gaps over always-stable surfaces without requiring explicit ### Arguments - `bnd`: NamedTuple with fields `psi`, `alpha`, `alpha_critical1`, `alpha_critical2` - (as returned by `ForceFreeStates.ballooning_alpha_boundaries`) + (as returned by `LocalStability.ballooning_alpha_boundaries`) ### Keyword arguments @@ -617,7 +617,7 @@ end plot_ballooning_alpha_boundaries(bnd, dpmap; save_path=nothing, psi_min=0.0) Same diagram drawn over a heatmap of the signed Δ' from -`ForceFreeStates.ballooning_delta_prime_map`: each surface's Δ'(α) is oriented by the +`LocalStability.ballooning_delta_prime_map`: each surface's Δ'(α) is oriented by the sign of its α=0 (stable) value so that positive is stable everywhere, regridded from its native physical α = α_ref*scale onto a shared uniform α axis, and shown with the Δ'=0 contour, the extracted boundaries, and the scan cap `max_alpha_scale*α_exp`. diff --git a/src/ForceFreeStates/ForceFreeStates.jl b/src/ForceFreeStates/ForceFreeStates.jl index 4e6d3c0f1..a88f9e9b8 100644 --- a/src/ForceFreeStates/ForceFreeStates.jl +++ b/src/ForceFreeStates/ForceFreeStates.jl @@ -24,7 +24,6 @@ import StaticArrays: @MMatrix # Include all necessary files include("ForceFreeStatesStructs.jl") -include("Ballooning.jl") include("Resist.jl") include("EulerLagrange.jl") include("Sing.jl") diff --git a/src/GeneralizedPerturbedEquilibrium.jl b/src/GeneralizedPerturbedEquilibrium.jl index 27ba49f43..cd7497dc0 100755 --- a/src/GeneralizedPerturbedEquilibrium.jl +++ b/src/GeneralizedPerturbedEquilibrium.jl @@ -20,6 +20,11 @@ include("Equilibrium/Equilibrium.jl") import .Equilibrium as Equilibrium export Equilibrium +# Local high-n stability (Mercier, resistive interchange, ballooning Δ'); depends only on Equilibrium. +include("LocalStability/LocalStability.jl") +import .LocalStability as LocalStability +export LocalStability + include("Vacuum/Vacuum.jl") import .Vacuum as Vacuum export Vacuum @@ -64,7 +69,6 @@ include("Rerun.jl") # Import ForceFreeStates types and functions needed for main using .ForceFreeStates: ForceFreeStatesInternal, ForceFreeStatesControl, DebugSettings, FreeBoundaryResult, OdeState, FourFitVars using .ForceFreeStates: sing_lim!, sing_min!, sing_find!, resist_eval_all!, resist_geometry, ResistGeometry -using .ForceFreeStates: compute_local_stability, compute_ballooning_stability!, ballooning_alpha_boundary, ballooning_alpha_boundaries using .ForceFreeStates: make_metric, make_matrix, make_kinetic_matrix using .ForceFreeStates: find_kinetic_singular_surfaces! using .ForceFreeStates: eulerlagrange_integration, free_run, normalize_eigenfunctions! @@ -336,9 +340,9 @@ function main_from_inputs( locstab = nothing ballooning_boundary = (psi=Float64[], alpha=Float64[], alpha_critical=Float64[]) if ctrl.local_stability_flag - locstab = compute_local_stability(ctrl, equil) + locstab = LocalStability.compute_local_stability(equil; verbose=ctrl.verbose) # First ballooning stability boundary (α vs ψ_N) for BALOO-style diagnostics. - ballooning_boundary = ballooning_alpha_boundary(ctrl, equil) + ballooning_boundary = LocalStability.ballooning_alpha_boundary(equil; verbose=ctrl.verbose) end # Find all singular surfaces in the equilibrium diff --git a/src/ForceFreeStates/Ballooning.jl b/src/LocalStability/Ballooning.jl similarity index 96% rename from src/ForceFreeStates/Ballooning.jl rename to src/LocalStability/Ballooning.jl index 51da221be..95e7c0a08 100644 --- a/src/ForceFreeStates/Ballooning.jl +++ b/src/LocalStability/Ballooning.jl @@ -2,8 +2,6 @@ # Main Driver for Ballooning Stability Analysis # Computes ballooning stability criterion over all flux surfaces # ====================================================================== -using LinearAlgebra -using StaticArrays: SVector const BALLOONING_THETA_MAX_CAP = 16.5 const BALLOONING_THETA_SCALE_MULTIPLIER = 10.0 @@ -17,7 +15,7 @@ _in_ballooning_scan_window(psi::Float64, psi_edge::Float64) = BALLOONING_SCAN_PSI_WINDOW[1] <= psi <= min(BALLOONING_SCAN_PSI_WINDOW[2], psi_edge) """ - compute_ballooning_stability!(ctrl, locstab_fs, plasma_eq) + compute_ballooning_stability!(locstab_fs, plasma_eq) Main driver routine for local high-n stability analysis. Iterates over all magnetic flux surfaces, prepares ballooning coefficients, stores `det(d0bar)` @@ -26,9 +24,9 @@ to compute Delta Prime. ## Arguments - - `ctrl::ForceFreeStatesControl`: Control parameters for the analysis. - `locstab_fs::Matrix{Float64}`: Local stability matrix to store results (modified in place). - `plasma_eq::Equilibrium.PlasmaEquilibrium`: Plasma equilibrium data. + - `verbose::Bool`: Print progress messages. This function modifies `locstab_fs` in place with: @@ -38,14 +36,14 @@ This function modifies `locstab_fs` in place with: - Column 4: Delta Prime (Δ') """ function compute_ballooning_stability!( - ctrl::ForceFreeStatesControl, locstab_fs::Matrix{Float64}, plasma_eq::Equilibrium.PlasmaEquilibrium; theta_k::Float64=0.0, - compute_delta_prime::Bool=true + compute_delta_prime::Bool=true, + verbose::Bool=false ) - if ctrl.verbose + if verbose println("Evaluating local high-n ballooning stability...") end @@ -72,22 +70,22 @@ function compute_ballooning_stability!( end end - if ctrl.verbose + if verbose println("Ballooning analysis complete.") end end """ - compute_local_stability(ctrl, plasma_eq) -> CubicSeriesInterpolant + 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(ctrl::ForceFreeStatesControl, plasma_eq::Equilibrium.PlasmaEquilibrium) +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!(ctrl, locstab_fs, plasma_eq) + compute_ballooning_stability!(locstab_fs, plasma_eq; verbose=verbose) return cubic_interp(xs, Series(locstab_fs); extrap=ExtendExtrap()) end @@ -669,10 +667,10 @@ ballooning evaluator, and returns delta-prime and `det(d0bar)` `Di` maps. function scan_delta_prime_map( psi_idx::Int, plasma_eq::Equilibrium.PlasmaEquilibrium; - ctrl::ForceFreeStatesControl=ForceFreeStatesControl(; verbose=false), theta_k::Float64=0.0, s_scales::AbstractVector{<:Real}, - alpha_scales::AbstractVector{<:Real} + alpha_scales::AbstractVector{<:Real}, + verbose::Bool=false ) ref = salpha_reference(psi_idx, plasma_eq) @@ -697,7 +695,7 @@ function scan_delta_prime_map( theta_k=theta_k ) catch err - if ctrl.verbose + if verbose @warn "s-alpha scan point failed" psi_idx is ia err end nothing @@ -841,7 +839,7 @@ function critical_ballooning_alpha( end """ - ballooning_alpha_boundary(ctrl, plasma_eq; theta_k=0.0, n_scan=24) + ballooning_alpha_boundary(plasma_eq; theta_k=0.0, n_scan=24, verbose=false) Profile driver for the BALOO-style ballooning stability diagram. Loops over flux surfaces returning the experimental pressure gradient `alpha` (from @@ -861,10 +859,10 @@ far-edge surfaces dominate the scan cost. Per-surface failures, skipped surfaces surfaces with no boundary within range are returned as `NaN`. """ function ballooning_alpha_boundary( - ctrl::ForceFreeStatesControl, plasma_eq::Equilibrium.PlasmaEquilibrium; theta_k::Float64=0.0, - n_scan::Int=24 + n_scan::Int=24, + verbose::Bool=false ) xs = plasma_eq.profiles.xs npsi = length(xs) @@ -878,7 +876,7 @@ function ballooning_alpha_boundary( 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 catch err - if ctrl.verbose + if verbose @warn "ballooning alpha boundary failed" psi_idx = i exception = err end end @@ -910,7 +908,7 @@ function second_critical_ballooning_alpha( end """ - ballooning_alpha_boundaries(ctrl, plasma_eq; theta_k=0.0, max_alpha_scale=8.0, n_scan=24) + ballooning_alpha_boundaries(plasma_eq; theta_k=0.0, max_alpha_scale=8.0, n_scan=24, verbose=false) Profile driver returning the experimental pressure gradient `alpha`, the first stability boundary `alpha_critical1` (lowest Δ' zero), and the second stability boundary @@ -924,11 +922,11 @@ smooth boundary curves. Surfaces outside `BALLOONING_SCAN_PSI_WINDOW` are skippe (returned as `NaN`), like per-surface failures. """ function ballooning_alpha_boundaries( - ctrl::ForceFreeStatesControl, plasma_eq::Equilibrium.PlasmaEquilibrium; theta_k::Float64=0.0, max_alpha_scale::Float64=8.0, - n_scan::Int=24 + n_scan::Int=24, + verbose::Bool=false ) xs = plasma_eq.profiles.xs npsi = length(xs) @@ -946,7 +944,7 @@ function ballooning_alpha_boundaries( alpha_critical1[i] = cr.alphas[1] length(cr.alphas) >= 2 && (alpha_critical2[i] = cr.alphas[2]) catch err - if ctrl.verbose + if verbose @warn "ballooning alpha boundaries failed" psi_idx=i exception=err end end @@ -1002,7 +1000,7 @@ function ballooning_qprime_crossings( end """ - ballooning_qprime_boundaries(ctrl, plasma_eq; theta_k=0.0, min_qprime_scale=-2.0, max_qprime_scale=4.0, n_scan=24) + ballooning_qprime_boundaries(plasma_eq; theta_k=0.0, min_qprime_scale=-2.0, max_qprime_scale=4.0, n_scan=24, verbose=false) Profile driver returning the experimental shear profile `qprime` (`dq/dpsi_norm`), the critical shear `qprime_critical1` below which the surface is ballooning unstable at the @@ -1012,12 +1010,12 @@ normalized flux `psi`, from [`ballooning_qprime_crossings`](@ref) at each surfac Arrays contain `NaN` where no boundary exists. """ function ballooning_qprime_boundaries( - ctrl::ForceFreeStatesControl, plasma_eq::Equilibrium.PlasmaEquilibrium; theta_k::Float64=0.0, min_qprime_scale::Float64=-2.0, max_qprime_scale::Float64=4.0, - n_scan::Int=24 + n_scan::Int=24, + verbose::Bool=false ) xs = plasma_eq.profiles.xs npsi = length(xs) @@ -1037,7 +1035,7 @@ function ballooning_qprime_boundaries( qprime_critical1[i] = cr.qprimes[1] length(cr.qprimes) >= 2 && (qprime_critical2[i] = cr.qprimes[2]) catch err - if ctrl.verbose + if verbose @warn "ballooning qprime boundaries failed" psi_idx=i exception=err end end @@ -1047,7 +1045,7 @@ function ballooning_qprime_boundaries( end """ - ballooning_delta_prime_map(ctrl, plasma_eq; theta_k=0.0, max_alpha_scale=8.0, n_alpha=61, max_surfaces=40) + ballooning_delta_prime_map(plasma_eq; theta_k=0.0, max_alpha_scale=8.0, n_alpha=61, max_surfaces=40, verbose=false) Map the ballooning Δ' over the (ψ_N, α) plane at fixed magnetic shear: at each flux surface, scan the pressure-gradient scaling from 0 to `max_alpha_scale` times the @@ -1058,12 +1056,12 @@ Returns a NamedTuple `(psi, alpha_scales, alpha_ref, delta_prime)` where `delta_ is `[n_surf × n_alpha]` with `NaN` on failed surfaces; physical `α = alpha_ref * scale`. """ function ballooning_delta_prime_map( - ctrl::ForceFreeStatesControl, plasma_eq::Equilibrium.PlasmaEquilibrium; theta_k::Float64=0.0, max_alpha_scale::Float64=8.0, n_alpha::Int=61, - max_surfaces::Int=40 + max_surfaces::Int=40, + verbose::Bool=false ) xs = plasma_eq.profiles.xs idx_all = [i for i in eachindex(xs) if xs[i] <= 1.0] @@ -1077,11 +1075,11 @@ function ballooning_delta_prime_map( for (k, i) in enumerate(idx) try - res = scan_delta_prime_map(i, plasma_eq; ctrl=ctrl, theta_k=theta_k, s_scales=[1.0], alpha_scales=alpha_scales) + res = scan_delta_prime_map(i, plasma_eq; theta_k=theta_k, s_scales=[1.0], alpha_scales=alpha_scales, verbose=verbose) alpha_ref[k] = res.reference.alpha_ref delta_prime[k, :] .= vec(res.delta_prime) catch err - if ctrl.verbose + if verbose @warn "ballooning delta-prime map failed" psi_idx=i exception=err end end @@ -1091,7 +1089,7 @@ function ballooning_delta_prime_map( end """ - ballooning_qprime_delta_prime_map(ctrl, plasma_eq; theta_k=0.0, min_qprime_scale=-2.0, max_qprime_scale=4.0, n_qprime=61, max_surfaces=40) + ballooning_qprime_delta_prime_map(plasma_eq; theta_k=0.0, min_qprime_scale=-2.0, max_qprime_scale=4.0, n_qprime=61, max_surfaces=40, verbose=false) Map the ballooning Δ' over the (ψ_N, q') plane at the fixed experimental pressure gradient: at each flux surface, scan the magnetic shear scaling from `min_qprime_scale` @@ -1104,13 +1102,13 @@ Returns a NamedTuple `(psi, qprime_scales, qprime_ref, delta_prime)` where `q' = qprime_ref * scale` in `dq/dpsi_norm` units. """ function ballooning_qprime_delta_prime_map( - ctrl::ForceFreeStatesControl, plasma_eq::Equilibrium.PlasmaEquilibrium; theta_k::Float64=0.0, min_qprime_scale::Float64=-2.0, max_qprime_scale::Float64=4.0, n_qprime::Int=61, - max_surfaces::Int=40 + max_surfaces::Int=40, + verbose::Bool=false ) xs = plasma_eq.profiles.xs idx_all = [i for i in eachindex(xs) if xs[i] <= 1.0] @@ -1124,11 +1122,11 @@ function ballooning_qprime_delta_prime_map( for (k, i) in enumerate(idx) try - res = scan_delta_prime_map(i, plasma_eq; ctrl=ctrl, theta_k=theta_k, s_scales=qprime_scales, alpha_scales=[1.0]) + res = scan_delta_prime_map(i, plasma_eq; theta_k=theta_k, s_scales=qprime_scales, alpha_scales=[1.0], verbose=verbose) qprime_ref[k] = res.reference.qprime_norm_ref delta_prime[k, :] .= vec(res.delta_prime) catch err - if ctrl.verbose + if verbose @warn "ballooning qprime delta-prime map failed" psi_idx=i exception=err end end diff --git a/src/LocalStability/LocalStability.jl b/src/LocalStability/LocalStability.jl new file mode 100644 index 000000000..eb96780e2 --- /dev/null +++ b/src/LocalStability/LocalStability.jl @@ -0,0 +1,17 @@ +module LocalStability + +# Local high-n stability: Mercier D_I, resistive interchange D_R, and the ballooning Δ' +# scans. Depends only on Equilibrium (plus math libraries) — no stability-solver state. +using LinearAlgebra +using FFTW +using OrdinaryDiffEq +using FastInterpolations +using StaticArrays: SVector + +import ..Equilibrium + +include("Ballooning.jl") + +export compute_local_stability, compute_ballooning_stability!, ballooning_alpha_boundary, ballooning_alpha_boundaries + +end diff --git a/test/runtests.jl b/test/runtests.jl index 76eca69d1..7e0e46d5d 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -3,6 +3,7 @@ using Pkg using GeneralizedPerturbedEquilibrium.Vacuum using GeneralizedPerturbedEquilibrium.Equilibrium using GeneralizedPerturbedEquilibrium.ForceFreeStates +using GeneralizedPerturbedEquilibrium.LocalStability using GeneralizedPerturbedEquilibrium.ForcingTerms using GeneralizedPerturbedEquilibrium.PerturbedEquilibrium using GeneralizedPerturbedEquilibrium.Utilities diff --git a/test/runtests_resist_eval.jl b/test/runtests_resist_eval.jl index 58aee3d9b..8e334a67a 100644 --- a/test/runtests_resist_eval.jl +++ b/test/runtests_resist_eval.jl @@ -3,6 +3,7 @@ using GeneralizedPerturbedEquilibrium.Equilibrium using GeneralizedPerturbedEquilibrium.ForceFreeStates using GeneralizedPerturbedEquilibrium.ForceFreeStates: SingType, ResistGeometry + using GeneralizedPerturbedEquilibrium.LocalStability using GeneralizedPerturbedEquilibrium.Utilities using GeneralizedPerturbedEquilibrium.InnerLayer using GeneralizedPerturbedEquilibrium.Tearing: build_ggj_inputs @@ -44,7 +45,7 @@ # scan stores in locstab[:,1]. Build it on the radial grid, interpolate # to a few surface ψ values, and check against the GGJ reconstruction. xs = equil.profiles.xs - di_ref = Float64[ForceFreeStates.prepare_ballooning_coefficients(i, equil).di for i in eachindex(xs)] + di_ref = Float64[LocalStability.prepare_ballooning_coefficients(i, equil).di for i in eachindex(xs)] di_spline = cubic_interp(xs, di_ref) dq = deriv_view(equil.profiles.q_spline, 1)