From fd3580e61bc5ff28bb4b2161d8a280a86e47ded1 Mon Sep 17 00:00:00 2001 From: d-burg Date: Sat, 7 Nov 2026 11:18:43 -0600 Subject: [PATCH 01/12] ForceFreeStates - FEATURE - Add perpendicular inertia and a growth-rate solve MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit ELITE's inertia model in the force-balance equations: the perpendicular kinetic energy K = ½∫ρ|ξ_⊥|² dV, with parallel inertia dropped, is algebraic in the normal displacement Ξ = ξ·∇ψ and the in-surface variable ξ_s = (ξ×B)·e_ψ because ξ×B = ξ_s∇ψ + Ξ χ′(∇ζ − q∇θ), so it adds Hermitian bands ⟨μ₀ρ|∇ψ|²/B²⟩, ⟨μ₀ρ u·∇ψ/B²⟩ and ⟨μ₀ρ|u|²/B²⟩ to A, C and H (Glasser 2016 Eqs. 8-9, A6). PerpendicularInertia is a `kinetic` source with kinetic_factor = γ² in s⁻², so the reduction runs on the kinetic (non-Q-factorized) path, where the rational surfaces are regular for γ² > 0 and the integration is one chunk. growth_rate brackets and Brent-solves the γ² at which the least-stable free-boundary energy vanishes, each evaluation one solve. Prototype: constant or user-supplied ρ(ψ), no export yet, no tests. Co-Authored-By: Claude Fable 5.1 --- src/ForceFreeStates/ForceBalance.jl | 7 +- src/ForceFreeStates/ForceFreeStates.jl | 1 + src/ForceFreeStates/Inertia.jl | 168 +++++++++++++++++++++++++ src/ForceFreeStates/Kinetic.jl | 3 + 4 files changed, 177 insertions(+), 2 deletions(-) create mode 100644 src/ForceFreeStates/Inertia.jl diff --git a/src/ForceFreeStates/ForceBalance.jl b/src/ForceFreeStates/ForceBalance.jl index f4db1a240..cf4aba8e5 100644 --- a/src/ForceFreeStates/ForceBalance.jl +++ b/src/ForceFreeStates/ForceBalance.jl @@ -133,9 +133,9 @@ function _prepare_solve(prob::ForceBalanceProblem, alg::AbstractIntegrator; verb error("kinetic_factor scales the kinetic matrices, but no kinetic source is set: pass kinetic=... " * "([ForceFreeStates] kinetic = \"KineticForces\")") # A kinetic model (not supplied or synthetic matrices) is evaluated on the equilibrium's kinetic profiles. - ctrl.kinetic === nothing || ctrl.kinetic isa Union{KineticMatrices,TestKineticMatrices} || equil.kinetic !== nothing || + ctrl.kinetic === nothing || ctrl.kinetic isa Union{KineticMatrices,TestKineticMatrices,PerpendicularInertia} || equil.kinetic !== nothing || error("the KineticForces kinetic source needs kinetic profiles on the equilibrium: set kinetic_file or call attach_kinetic_profiles!") - ctrl.kinetic isa Union{Nothing,KineticMatrices,TestKineticMatrices} || Equilibrium.has_ntv_fields(equil.kinetic) || + ctrl.kinetic isa Union{Nothing,KineticMatrices,TestKineticMatrices,PerpendicularInertia} || Equilibrium.has_ntv_fields(equil.kinetic) || error("the KineticForces kinetic source needs n_e, T_i, T_e and omega_E in the kinetic profiles; the kinetic file lacks at least one") ctrl.kinetic !== nothing && !(alg isa Forward) && error("$(nameof(typeof(alg)))() does not support kinetic runs; use Forward().") @@ -272,6 +272,9 @@ function prepare_force_free_states!( elseif ctrl.kinetic isa TestKineticMatrices intr.verbose && @info "Computing kinetic matrices (synthetic test matrices, scale $(ctrl.kinetic.scale))" km = evaluate_kinetic_matrices(ctrl, equil, mats, intr, metric) + elseif ctrl.kinetic isa PerpendicularInertia + intr.verbose && @info "Computing perpendicular-inertia matrices (γ² = $(ctrl.kinetic_factor) s⁻²)" + km = evaluate_kinetic_matrices(ctrl, equil, mats, intr, metric) elseif ctrl.kinetic !== nothing intr.verbose && @info "Computing kinetic matrices (KineticForces model)" calculated_cb, axis_psi_c = kinetic_model_source(ctrl.kinetic, equil, intr) diff --git a/src/ForceFreeStates/ForceFreeStates.jl b/src/ForceFreeStates/ForceFreeStates.jl index f1374dae3..e2f9bc4ee 100644 --- a/src/ForceFreeStates/ForceFreeStates.jl +++ b/src/ForceFreeStates/ForceFreeStates.jl @@ -72,6 +72,7 @@ include("InputSnapshot.jl") include("Output.jl") include("Run.jl") include("ForceBalance.jl") +include("Inertia.jl") # Small tolerance and root-finder iteration cap used throughout ForceFreeStates global eps = 1e-10 diff --git a/src/ForceFreeStates/Inertia.jl b/src/ForceFreeStates/Inertia.jl new file mode 100644 index 000000000..f8f3b5105 --- /dev/null +++ b/src/ForceFreeStates/Inertia.jl @@ -0,0 +1,168 @@ +""" + PerpendicularInertia(; rho) + +Perpendicular plasma inertia for the force-balance equations, the ELITE inertia model: the kinetic +energy K = ½∫ρ|ξ_⊥|² dV with the parallel inertia dropped (incompressible displacement). With the +normal displacement Ξ = ξ·∇ψ and the in-surface variable ξ_s = (ξ×B)·e_ψ of the energy principle +[Glasser 2016 Eqs. 8-9], ξ×B = ξ_s∇ψ + Ξu with u = χ′(∇ζ − q∇θ), so |ξ_⊥|² = |ξ×B|²/B² is algebraic +in the two retained variables and K adds the Hermitian bands ⟨μ₀ρ|∇ψ|²/B²⟩ to A, ⟨μ₀ρ u·∇ψ/B²⟩ to +C and ⟨μ₀ρ|u|²/B²⟩ to H of [Glasser 2016 Eq. A6], each scaled by γ² for a mode ∝ e^{γt} +(δW + γ²K is stationary). + +Pass it as the `kinetic` source of a [`ForceBalanceProblem`](@ref) with `kinetic_factor = γ²` in +s⁻², or let [`growth_rate`](@ref) find the γ² at which the free-boundary energy vanishes. The F, K, +G reduction then takes the kinetic (non-Q-factorized) path, on which the rational surfaces are +regular points for γ² > 0 and no ideal crossings are needed. + +## Fields + + - `rho` - Mass density in kg/m³: a number, or a function of normalized flux ψ_N. +""" +@kwdef struct PerpendicularInertia + rho::Union{Float64,Function} +end + +_rho_at(src::PerpendicularInertia, psi::Real) = src.rho isa Function ? Float64(src.rho(psi)) : src.rho + +""" + inertia_kinetic_matrices(src::PerpendicularInertia, equil, intr, metric) -> (kw_flat, kt_flat) + +The perpendicular-inertia bands on the metric's (ψ, θ) grid in the (ψ, mode pair, component) layout +of [`KineticMatrices`](@ref): components A, C and H carry ⟨μ₀ρ|∇ψ|²/B²⟩, ⟨μ₀ρ u·∇ψ/B²⟩ and +⟨μ₀ρ|u|²/B²⟩ per unit γ², the others and the torque matrices are zero. The metric stores +J·g_ij with g_ij = e_i·e_j/J² [Glasser 2016 Eq. A1], so the contravariant products are cofactors: +|∇ψ|² = g22g33 − g23², ∇ψ·∇θ = g23g31 − g12g33, ∇ψ·∇ζ = g12g23 − g22g31, |∇θ|² = g11g33 − g31², +|∇ζ|² = g11g22 − g12², ∇θ·∇ζ = g31g12 − g11g23 (all in the stored J·g units), and +B² = χ′²(g22 + 2q g23 + q²g33)/J from B = (χ′/J)(e_θ + q e_ζ) [Glasser 2016 Eq. 5]. +""" +function inertia_kinetic_matrices(src::PerpendicularInertia, equil::Equilibrium.PlasmaEquilibrium, intr::ForceFreeStatesInternal, metric::MetricData) + mpsi, mtheta, mpert = metric.mpsi, metric.mtheta, intr.mpert + chi1 = 2π * equil.psio + integrands = zeros(Float64, mpsi, mtheta, 3) + for ipsi in 1:mpsi + q = equil.profiles.q_spline(metric.xs[ipsi]) + mu0rho = Utilities.MU_0 * _rho_at(src, metric.xs[ipsi]) + for jtheta in 1:mtheta + g11, g22, g33, g23, g31, g12, jac = @view metric.fs[ipsi, jtheta, 1:7] + gradpsi2 = g22 * g33 - g23^2 + gradpsi_theta = g23 * g31 - g12 * g33 + gradpsi_zeta = g12 * g23 - g22 * g31 + gradtheta2 = g11 * g33 - g31^2 + gradzeta2 = g11 * g22 - g12^2 + gradtheta_zeta = g31 * g12 - g11 * g23 + b2 = chi1^2 * (g22 + 2q * g23 + q^2 * g33) / jac + weight = jac * mu0rho / b2 + integrands[ipsi, jtheta, 1] = weight * gradpsi2 + integrands[ipsi, jtheta, 2] = weight * chi1 * (gradpsi_zeta - q * gradpsi_theta) + integrands[ipsi, jtheta, 3] = weight * chi1^2 * (gradzeta2 - 2q * gradtheta_zeta + q^2 * gradtheta2) + end + end + fc = Utilities.FourierCoefficients(metric.xs, metric.ys, integrands, mpert - 1) + + np = intr.numpert_total + kw_flat = zeros(ComplexF64, mpsi, np^2, 6) + kt_flat = zeros(ComplexF64, mpsi, np^2, 6) + band = zeros(ComplexF64, 2mpert - 1) + for ipsi in 1:mpsi, (qty, ic) in ((1, 1), (2, 3), (3, 6)) + # Same band layout as the ideal matrices: coefficient m at mpert−m, conjugates fill the other half + for m in 0:(mpert-1) + band[mpert-m] = Utilities.get_complex_coeff(fc, ipsi, m, qty) + end + for k in 1:(mpert-1) + band[k+mpert] = conj(band[mpert-k]) + end + for ipert_n in 1:intr.npert, ipert_m in 1:mpert, jpert_m in 1:mpert + ipert = ipert_m + (ipert_n - 1) * mpert + jpert = jpert_m + (ipert_n - 1) * mpert + kw_flat[ipsi, ipert+(jpert-1)*np, ic] = band[jpert_m-ipert_m+mpert] + end + end + return kw_flat, kt_flat +end + +""" + GrowthRateResult + +The outcome of [`growth_rate`](@ref). + +## Fields + + - `gamma::Float64` - Growth rate in s⁻¹ of the least-stable free-boundary mode; 0 when the plasma is marginally stable or stable. + - `gamma_sq::Float64` - γ² at which the least-stable free-boundary energy eigenvalue crosses zero. + - `marginal::ForceFreeStatesResult` - The ideal (γ² = 0) solve, whose least-stable energy decides stability. + - `result::Union{Nothing,ForceFreeStatesResult}` - The solve at `gamma_sq`, carrying the eigenmode; `nothing` when stable. + - `history::Vector` - Every (γ², λ_min, nzero) evaluated, in order, with λ_min the least-stable free-boundary energy at that γ². +""" +struct GrowthRateResult + gamma::Float64 + gamma_sq::Float64 + marginal::ForceFreeStatesResult + result::Union{Nothing,ForceFreeStatesResult} + history::Vector{@NamedTuple{gamma_sq::Float64, lambda::Float64, nzero::Int}} +end + +""" + growth_rate(prob::ForceBalanceProblem, alg::Forward; inertia, gamma_sq_init=1e9, gamma_sq_min=1e4, gamma_sq_max=1e16, + rtol=1e-3, verbose=false) -> GrowthRateResult + +Growth rate of the least-stable free-boundary mode of `prob` with the perpendicular inertia `inertia` +(a [`PerpendicularInertia`](@ref)), ELITE's eigenvalue problem solved by shooting: the ideal solve +gives the marginal least-stable energy λ_min(0); when it is negative, the inertia enters as +`kinetic_factor = γ²` and λ_min(γ²), which rises monotonically with γ² once the interior is +positive, is brought to zero. The root is bracketed by decades from `gamma_sq_init` and closed by +Brent's method to a relative tolerance `rtol` in γ²; every evaluation is one force-balance solve. +`gamma_sq_min` is the floor below which the mode counts as marginal, `gamma_sq_max` the ceiling +beyond which bracketing stops with an error. + +The problem's `singfac_min` is set to 0 for the inertial solves unless given: with inertia the +rational surfaces are regular, so the integration runs as a single chunk. A fixed-boundary mode +faster than the free-boundary one shows as `nzero > 0` at the root in `history`; this first version +reports the free-boundary root and leaves the interior count to the caller. +""" +function growth_rate(prob::ForceBalanceProblem, alg::Forward; inertia::PerpendicularInertia, + gamma_sq_init::Real=1e9, gamma_sq_min::Real=1e4, gamma_sq_max::Real=1e16, rtol::Real=1e-3, verbose::Bool=false) + history = @NamedTuple{gamma_sq::Float64, lambda::Float64, nzero::Int}[] + lambda_min(r) = minimum(real, r.free_boundary.et) + + marginal = solve(prob, alg; verbose) + marginal.free_boundary === nothing && error("growth_rate needs the free-boundary energies: set vac_flag=true") + push!(history, (gamma_sq=0.0, lambda=lambda_min(marginal), nzero=marginal.diagnostics.nzero)) + verbose && @info @sprintf("growth_rate: marginal λ_min = %.6e (nzero = %d)", history[end].lambda, history[end].nzero) + history[end].lambda < 0 || return GrowthRateResult(0.0, 0.0, marginal, nothing, history) + + kwargs = copy(prob.ctrl_kwargs) + haskey(kwargs, :singfac_min) || (kwargs[:singfac_min] = 0.0) + kwargs[:kinetic] = inertia + matrices = Ref{Any}(nothing) # the inertia bands, evaluated once and reused at every γ² + last = Ref{Any}(nothing) + function lambda_at(gamma_sq) + kw = copy(kwargs) + matrices[] === nothing || (kw[:kinetic] = matrices[]) + kw[:kinetic_factor] = gamma_sq + inertial = ForceBalanceProblem(prob.equil, prob.wall, kw) + matrices[] === nothing && (matrices[] = kinetic_matrices(inertial, alg; verbose)) + r = solve(inertial, alg; verbose) + last[] = r + push!(history, (gamma_sq=Float64(gamma_sq), lambda=lambda_min(r), nzero=r.diagnostics.nzero)) + verbose && @info @sprintf("growth_rate: γ² = %.4e λ_min = %.6e nzero = %d", gamma_sq, history[end].lambda, history[end].nzero) + return history[end].lambda + end + + # Bracket by decades: down while already positive, up while still negative + hi = Float64(gamma_sq_init) + lam = lambda_at(hi) + while lam > 0 + hi /= 10 + hi < gamma_sq_min && return GrowthRateResult(sqrt(hi), hi, marginal, last[], history) + lam = lambda_at(hi) + end + lo = hi + while lam <= 0 + lo = hi + hi *= 10 + hi > gamma_sq_max && error("growth_rate: λ_min is still negative at γ² = $lo; raise gamma_sq_max or check the inertia density") + lam = lambda_at(hi) + end + gamma_sq = find_zero(lambda_at, (lo, hi), Roots.Brent(); xrtol=rtol) + return GrowthRateResult(sqrt(gamma_sq), gamma_sq, marginal, last[], history) +end diff --git a/src/ForceFreeStates/Kinetic.jl b/src/ForceFreeStates/Kinetic.jl index 51351e556..b1b831a24 100644 --- a/src/ForceFreeStates/Kinetic.jl +++ b/src/ForceFreeStates/Kinetic.jl @@ -142,6 +142,9 @@ function evaluate_kinetic_matrices( if ctrl.kinetic isa TestKineticMatrices kw_flat, kt_flat = fixed_kinetic_matrices(intr.mpert, intr.numpert_total, mpsi, ctrl.kinetic.scale, intr.mlow, mats, xs) return matrices(collect(xs), kw_flat, kt_flat) + elseif ctrl.kinetic isa PerpendicularInertia + kw_flat, kt_flat = inertia_kinetic_matrices(ctrl.kinetic, equil, intr, metric) + return matrices(collect(xs), kw_flat, kt_flat) end isnothing(calculated_source) && error( "The KineticForces kinetic source requires the KineticForces callback: solve the problem, or pass " * From 618cb93cf53720aad50eb94ce8a93fb39a1f11ab Mon Sep 17 00:00:00 2001 From: d-burg Date: Sat, 7 Nov 2026 11:41:14 -0600 Subject: [PATCH 02/12] ForceFreeStates - TEST - Cover the perpendicular-inertia bands, solves and growth rate MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Bands Hermitian and A, H positive definite, linear in ρ; the metric B² identity the bands rest on; the inertial solve regular with no crossings, λ_min approaching the ideal value as √γ² (the inertial-layer term is linear in γ) and rising with γ²; the kinetic crossing machinery and the single chunk agreeing; growth_rate returning γ = 0 on the stable fixture. Memoize growth_rate's evaluations so Brent's endpoint re-evaluations are free. Co-Authored-By: Claude Fable 5.1 --- src/ForceFreeStates/Inertia.jl | 4 +- test/runtests_inertia.jl | 101 +++++++++++++++++++++++++++++++++ test/suites.jl | 1 + 3 files changed, 105 insertions(+), 1 deletion(-) create mode 100644 test/runtests_inertia.jl diff --git a/src/ForceFreeStates/Inertia.jl b/src/ForceFreeStates/Inertia.jl index f8f3b5105..b7fc301cd 100644 --- a/src/ForceFreeStates/Inertia.jl +++ b/src/ForceFreeStates/Inertia.jl @@ -135,7 +135,9 @@ function growth_rate(prob::ForceBalanceProblem, alg::Forward; inertia::Perpendic kwargs[:kinetic] = inertia matrices = Ref{Any}(nothing) # the inertia bands, evaluated once and reused at every γ² last = Ref{Any}(nothing) + seen = Dict{Float64,Float64}() # Brent re-evaluates its bracket ends; each γ² is one solve function lambda_at(gamma_sq) + haskey(seen, gamma_sq) && return seen[gamma_sq] kw = copy(kwargs) matrices[] === nothing || (kw[:kinetic] = matrices[]) kw[:kinetic_factor] = gamma_sq @@ -145,7 +147,7 @@ function growth_rate(prob::ForceBalanceProblem, alg::Forward; inertia::Perpendic last[] = r push!(history, (gamma_sq=Float64(gamma_sq), lambda=lambda_min(r), nzero=r.diagnostics.nzero)) verbose && @info @sprintf("growth_rate: γ² = %.4e λ_min = %.6e nzero = %d", gamma_sq, history[end].lambda, history[end].nzero) - return history[end].lambda + return seen[Float64(gamma_sq)] = history[end].lambda end # Bracket by decades: down while already positive, up while still negative diff --git a/test/runtests_inertia.jl b/test/runtests_inertia.jl new file mode 100644 index 000000000..365392aba --- /dev/null +++ b/test/runtests_inertia.jl @@ -0,0 +1,101 @@ +using TOML +using LinearAlgebra + +# Perpendicular inertia (ELITE's kinetic energy K = ½∫ρ|ξ_⊥|² dV) as a kinetic source of the force-balance +# equations: the A, C, H bands it adds, the B² identity they rely on, the γ²-scaled solves and growth_rate. +@testset "perpendicular inertia" begin + GPEC = GeneralizedPerturbedEquilibrium + FFS = GPEC.ForceFreeStates + template = joinpath(@__DIR__, "test_data", "regression_tj_ideal_example") + deck = TOML.parsefile(joinpath(template, "gpec.toml")) + equil = GPEC.Equilibrium.PlasmaEquilibrium(GPEC.equilibrium_config(deck, template), GPEC.analytic_config(deck)) + wall = GPEC.Vacuum.parse_wall(get(deck, "Wall", nothing)) + ffs_kwargs = Dict{Symbol,Any}(p.first => p.second for p in GPEC.force_balance_keywords(deck, template) if !(p.first in (:nn, :wall))) + ffs_kwargs[:vac_flag] = true + forward = Forward(; eulerlagrange_tolerance=1e-8) + + mlow, mhigh = -3, 4 + mpert = mhigh - mlow + 1 + metric = FFS.make_metric(equil, mpert) + intr = FFS.ForceFreeStatesInternal(; mpert, numpert_total=mpert, mlow, mhigh, nlow=1, nhigh=1, npert=1) + ipsi_mid = cld(metric.mpsi, 2) + bands(kw, ic, ipsi) = reshape(kw[ipsi, :, ic], mpert, mpert) + + @testset "A, C, H bands" begin + kw, kt = FFS.inertia_kinetic_matrices(FFS.PerpendicularInertia(; rho=1e-7), equil, intr, metric) + @test size(kw) == (metric.mpsi, mpert^2, 6) + @test all(iszero, kt) + @test all(iszero, kw[:, :, [2, 4, 5]]) + for ic in (1, 3, 6) + m = bands(kw, ic, ipsi_mid) + @test norm(m - m') <= 1e-12 * norm(m) + @test norm(m) > 0 + end + @test all(>(0), eigvals(Hermitian(bands(kw, 1, ipsi_mid)))) + @test all(>(0), eigvals(Hermitian(bands(kw, 6, ipsi_mid)))) + + kw2, _ = FFS.inertia_kinetic_matrices(FFS.PerpendicularInertia(; rho=2e-7), equil, intr, metric) + @test kw2 ≈ 2kw rtol = 1e-12 + + # A ψ-dependent density scales each surface's bands by ρ(ψ)/ρ₀ and so lowers them toward the edge + kwf, _ = FFS.inertia_kinetic_matrices(FFS.PerpendicularInertia(; rho=psi -> 1e-7 * (1 - psi)), equil, intr, metric) + for ipsi in (2, ipsi_mid, metric.mpsi - 1), ic in (1, 3, 6) + @test bands(kwf, ic, ipsi) ≈ (1 - metric.xs[ipsi]) * bands(kw, ic, ipsi) rtol = 1e-12 + end + diagA = [real(kwf[ipsi, 1, 1]) / real(kw[ipsi, 1, 1]) for ipsi in 1:metric.mpsi] + @test issorted(diagA; rev=true) + end + + @testset "B² from the metric matches the equilibrium field" begin + chi1 = 2π * equil.psio + bnodes = equil.eqfun_B.nodal_derivs.partials + @test size(bnodes)[2:3] == (metric.mpsi, metric.mtheta) + for ipsi in (3, ipsi_mid, metric.mpsi - 2), jtheta in (1, cld(metric.mtheta, 3), cld(metric.mtheta, 2)) + g11, g22, g33, g23, g31, g12, jac = metric.fs[ipsi, jtheta, 1:7] + q = equil.profiles.q_spline(metric.xs[ipsi]) + @test chi1^2 * (g22 + 2q * g23 + q^2 * g33) / jac ≈ bnodes[1, ipsi, jtheta]^2 rtol = 1e-3 + end + end + + lambda_min(r) = minimum(real, r.free_boundary.et) + ideal = solve(equil, forward; nn=1, wall, ffs_kwargs...) + inertial_kwargs = merge(ffs_kwargs, Dict{Symbol,Any}(:kinetic => FFS.PerpendicularInertia(; rho=1e-7), :singfac_min => 0.0)) + + @testset "inertial solves" begin + @test equil.kinetic === nothing + weak = solve(equil, forward; nn=1, wall, inertial_kwargs..., kinetic_factor=1e4) + @test weak.closure === :regular + @test weak.kinetic.kmsing == 0 + + # The inertial layers at the rational surfaces add a term linear in γ, so λ_min approaches the + # ideal value as √γ²: a decade in γ² shrinks the offset by √10, and the offset at γ² = 1e2 is ~1.5e-4 + # of λ (measured 6.5e-5 absolute on this fixture). + weaker = solve(equil, forward; nn=1, wall, inertial_kwargs..., kinetic_factor=1e2) + offset_weak = lambda_min(weak) - lambda_min(ideal) + offset_weaker = lambda_min(weaker) - lambda_min(ideal) + @test 0 < offset_weaker < 3e-4 * lambda_min(ideal) + @test 8 < offset_weak / offset_weaker < 12 + + strong = solve(equil, forward; nn=1, wall, inertial_kwargs..., kinetic_factor=1e10) + @test lambda_min(strong) > lambda_min(weak) + + # The kinetic crossing machinery (singfac_min > 0, det F̄ near-zeros) and the single chunk agree + crossed = solve(equil, forward; nn=1, wall, merge(inertial_kwargs, Dict{Symbol,Any}(:singfac_min => 1e-4))..., kinetic_factor=1e4) + @test crossed.kinetic.kmsing > 0 + @test lambda_min(crossed) ≈ lambda_min(weak) rtol = 1e-5 + end + + @testset "growth_rate" begin + @test_throws ErrorException solve(equil, forward; nn=1, wall, ffs_kwargs..., kinetic_factor=2.0) + + # The fixture is ideal-stable, so growth_rate stops after the marginal solve + @test lambda_min(ideal) > 0 + prob = ForceBalanceProblem(equil; nn=1, wall, ffs_kwargs...) + gr = FFS.growth_rate(prob, forward; inertia=FFS.PerpendicularInertia(; rho=1e-7)) + @test gr.gamma == 0 + @test gr.gamma_sq == 0 + @test gr.result === nothing + @test length(gr.history) == 1 + @test gr.history[1].lambda ≈ lambda_min(ideal) rtol = 1e-12 + end +end diff --git a/test/suites.jl b/test/suites.jl index dfe79f9e4..8374236d6 100644 --- a/test/suites.jl +++ b/test/suites.jl @@ -50,6 +50,7 @@ const TEST_SUITES = [ "runtests_slayer_runner.jl" => "errorfields", "runtests_kinetic.jl" => "core", "runtests_kinetic_validity.jl" => "core", + "runtests_inertia.jl" => "core", "runtests_multiion.jl" => "core", "runtests_fullruns.jl" => "api", "runtests_torque_response.jl" => "api", From f92c98faf31dfc9e98ef8b9c0c12e312ca62a202 Mon Sep 17 00:00:00 2001 From: d-burg Date: Sat, 7 Nov 2026 11:42:34 -0600 Subject: [PATCH 03/12] ForceFreeStates - MINOR - Export the inertia API and include it in the docs Co-Authored-By: Claude Fable 5.1 --- docs/src/stability.md | 2 +- src/GeneralizedPerturbedEquilibrium.jl | 3 ++- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/docs/src/stability.md b/docs/src/stability.md index aa86df2ec..98d61d697 100644 --- a/docs/src/stability.md +++ b/docs/src/stability.md @@ -149,5 +149,5 @@ The Galerkin solver's functions are documented on the [Galerkin Solver](galerkin ```@autodocs Modules = [GeneralizedPerturbedEquilibrium.ForceFreeStates] -Pages = ["ForceFreeStates.jl", "CoreTypes.jl", "Surfaces/Types.jl", "Riccati/Types.jl", "Matching/DeltaPrime.jl", "Result.jl", "Surfaces/ResistEval.jl", "Matching/ResonantMatch.jl", "Matching/MatchProblem.jl", "EulerLagrange.jl", "Surfaces/Finding.jl", "Surfaces/Asymptotics.jl", "Fourfit.jl", "Kinetic.jl", "FixedKineticMatrices.jl", "FixedBoundaryStability.jl", "Utils.jl", "Free.jl", "Riccati/Propagators.jl", "Riccati/Crossings.jl", "Riccati/DeltaPrimeBVP.jl", "Riccati/Driver.jl"] +Pages = ["ForceFreeStates.jl", "CoreTypes.jl", "Surfaces/Types.jl", "Riccati/Types.jl", "Matching/DeltaPrime.jl", "Result.jl", "Surfaces/ResistEval.jl", "Matching/ResonantMatch.jl", "Matching/MatchProblem.jl", "EulerLagrange.jl", "Surfaces/Finding.jl", "Surfaces/Asymptotics.jl", "Fourfit.jl", "Kinetic.jl", "Inertia.jl", "FixedKineticMatrices.jl", "FixedBoundaryStability.jl", "Utils.jl", "Free.jl", "Riccati/Propagators.jl", "Riccati/Crossings.jl", "Riccati/DeltaPrimeBVP.jl", "Riccati/Driver.jl"] ``` diff --git a/src/GeneralizedPerturbedEquilibrium.jl b/src/GeneralizedPerturbedEquilibrium.jl index bda2b7c8d..b753a8e68 100755 --- a/src/GeneralizedPerturbedEquilibrium.jl +++ b/src/GeneralizedPerturbedEquilibrium.jl @@ -71,7 +71,7 @@ include("Rerun.jl") # Scripting-API surface: the integrator selectors, the published result, the equilibrium # constructor and the forcing description, re-exported so a user needs one `using`. -using .ForceFreeStates: AbstractIntegrator, Forward, Riccati, Galerkin, TestKineticMatrices, KineticMatrices, ForceFreeStatesResult +using .ForceFreeStates: AbstractIntegrator, Forward, Riccati, Galerkin, TestKineticMatrices, KineticMatrices, ForceFreeStatesResult, PerpendicularInertia, GrowthRateResult, growth_rate using .ForceFreeStates: MatchProblem using .InnerLayer: GGJModel, SLAYERModel using .Tearing: TearingProblem, TearingConfig, TearingResult, CriticalResonantFieldControl @@ -305,6 +305,7 @@ export main, write_imas export solve, plasma_response, perturbed_equilibrium, kinetic_torque, error_fields export KineticProfileData, IonSpecies, attach_kinetic_profiles! export PlasmaEquilibrium, ForceBalanceProblem, Forward, Riccati, Galerkin, TestKineticMatrices, KineticMatrices, kinetic_matrices, ForceFreeStatesResult, RMPField, PlasmaResponse +export PerpendicularInertia, GrowthRateResult, growth_rate export MatchProblem, TearingProblem, TearingConfig, TearingResult, GGJModel, SLAYERModel export LocalStabilityProblem, Ballooning, LocalStabilityResult export KineticForcesControl, KineticForcesResult, ErrorFieldsResult From d0365caab74d5d9359888036e317438a105ae453 Mon Sep 17 00:00:00 2001 From: d-burg Date: Sat, 7 Nov 2026 11:50:48 -0600 Subject: [PATCH 04/12] ForceFreeStates - PERF - Close the growth-rate root with a bracketed secant MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Brent kept bisecting to its x-tolerance after |λ_min| was already 1e-11, ten extra solves per root on the kink benchmark; an Illinois iteration stopping on |λ_min| < ftol·|λ_min(0)| or a narrow bracket finds the same γ in 5–6 solves in total. Co-Authored-By: Claude Fable 5.1 --- src/ForceFreeStates/Inertia.jl | 34 +++++++++++++++++++++++++++------- 1 file changed, 27 insertions(+), 7 deletions(-) diff --git a/src/ForceFreeStates/Inertia.jl b/src/ForceFreeStates/Inertia.jl index b7fc301cd..ff0a5e738 100644 --- a/src/ForceFreeStates/Inertia.jl +++ b/src/ForceFreeStates/Inertia.jl @@ -103,16 +103,17 @@ end """ growth_rate(prob::ForceBalanceProblem, alg::Forward; inertia, gamma_sq_init=1e9, gamma_sq_min=1e4, gamma_sq_max=1e16, - rtol=1e-3, verbose=false) -> GrowthRateResult + rtol=1e-3, ftol=1e-6, maxiter=30, verbose=false) -> GrowthRateResult Growth rate of the least-stable free-boundary mode of `prob` with the perpendicular inertia `inertia` (a [`PerpendicularInertia`](@ref)), ELITE's eigenvalue problem solved by shooting: the ideal solve gives the marginal least-stable energy λ_min(0); when it is negative, the inertia enters as `kinetic_factor = γ²` and λ_min(γ²), which rises monotonically with γ² once the interior is -positive, is brought to zero. The root is bracketed by decades from `gamma_sq_init` and closed by -Brent's method to a relative tolerance `rtol` in γ²; every evaluation is one force-balance solve. -`gamma_sq_min` is the floor below which the mode counts as marginal, `gamma_sq_max` the ceiling -beyond which bracketing stops with an error. +positive, is brought to zero. The root is bracketed by decades from `gamma_sq_init` and closed by a +bracketed secant (Illinois) iteration that stops when |λ_min| falls below `ftol` times the marginal +|λ_min(0)| or the bracket is narrower than `rtol` in γ², within `maxiter` steps; every evaluation is +one force-balance solve. `gamma_sq_min` is the floor below which the mode counts as marginal, +`gamma_sq_max` the ceiling beyond which bracketing stops with an error. The problem's `singfac_min` is set to 0 for the inertial solves unless given: with inertia the rational surfaces are regular, so the integration runs as a single chunk. A fixed-boundary mode @@ -120,7 +121,8 @@ faster than the free-boundary one shows as `nzero > 0` at the root in `history`; reports the free-boundary root and leaves the interior count to the caller. """ function growth_rate(prob::ForceBalanceProblem, alg::Forward; inertia::PerpendicularInertia, - gamma_sq_init::Real=1e9, gamma_sq_min::Real=1e4, gamma_sq_max::Real=1e16, rtol::Real=1e-3, verbose::Bool=false) + gamma_sq_init::Real=1e9, gamma_sq_min::Real=1e4, gamma_sq_max::Real=1e16, rtol::Real=1e-3, ftol::Real=1e-6, maxiter::Int=30, + verbose::Bool=false) history = @NamedTuple{gamma_sq::Float64, lambda::Float64, nzero::Int}[] lambda_min(r) = minimum(real, r.free_boundary.et) @@ -165,6 +167,24 @@ function growth_rate(prob::ForceBalanceProblem, alg::Forward; inertia::Perpendic hi > gamma_sq_max && error("growth_rate: λ_min is still negative at γ² = $lo; raise gamma_sq_max or check the inertia density") lam = lambda_at(hi) end - gamma_sq = find_zero(lambda_at, (lo, hi), Roots.Brent(); xrtol=rtol) + + # Illinois iteration on the sign-change bracket; the retained end's λ is halved when it is kept twice running + flo, fhi = seen[lo], seen[hi] + gamma_sq, side = hi, 0 + for _ in 1:maxiter + gamma_sq = (lo * fhi - hi * flo) / (fhi - flo) + lam = lambda_at(gamma_sq) + abs(lam) <= ftol * abs(history[1].lambda) && break + if lam < 0 + lo, flo = gamma_sq, lam + side == -1 && (fhi /= 2) + side = -1 + else + hi, fhi = gamma_sq, lam + side == 1 && (flo /= 2) + side = 1 + end + hi - lo <= rtol * hi && break + end return GrowthRateResult(sqrt(gamma_sq), gamma_sq, marginal, last[], history) end From da12c9335503b2312aab98555235c1bab81c6868 Mon Sep 17 00:00:00 2001 From: d-burg Date: Sat, 7 Nov 2026 16:57:14 -0600 Subject: [PATCH 05/12] ForceFreeStates - MINOR - Rename the inertial result type to IdealGrowthRateResult The bare GrowthRateResult is Tearing's dispersion-root type; exporting a second one from ForceFreeStates made the top-level import ambiguous (precompile warning, import ignored). Co-Authored-By: Claude Fable 5.1 --- src/ForceFreeStates/Inertia.jl | 12 ++++++------ src/GeneralizedPerturbedEquilibrium.jl | 5 +++-- 2 files changed, 9 insertions(+), 8 deletions(-) diff --git a/src/ForceFreeStates/Inertia.jl b/src/ForceFreeStates/Inertia.jl index ff0a5e738..20cabfa51 100644 --- a/src/ForceFreeStates/Inertia.jl +++ b/src/ForceFreeStates/Inertia.jl @@ -81,7 +81,7 @@ function inertia_kinetic_matrices(src::PerpendicularInertia, equil::Equilibrium. end """ - GrowthRateResult + IdealGrowthRateResult The outcome of [`growth_rate`](@ref). @@ -93,7 +93,7 @@ The outcome of [`growth_rate`](@ref). - `result::Union{Nothing,ForceFreeStatesResult}` - The solve at `gamma_sq`, carrying the eigenmode; `nothing` when stable. - `history::Vector` - Every (γ², λ_min, nzero) evaluated, in order, with λ_min the least-stable free-boundary energy at that γ². """ -struct GrowthRateResult +struct IdealGrowthRateResult gamma::Float64 gamma_sq::Float64 marginal::ForceFreeStatesResult @@ -103,7 +103,7 @@ end """ growth_rate(prob::ForceBalanceProblem, alg::Forward; inertia, gamma_sq_init=1e9, gamma_sq_min=1e4, gamma_sq_max=1e16, - rtol=1e-3, ftol=1e-6, maxiter=30, verbose=false) -> GrowthRateResult + rtol=1e-3, ftol=1e-6, maxiter=30, verbose=false) -> IdealGrowthRateResult Growth rate of the least-stable free-boundary mode of `prob` with the perpendicular inertia `inertia` (a [`PerpendicularInertia`](@ref)), ELITE's eigenvalue problem solved by shooting: the ideal solve @@ -130,7 +130,7 @@ function growth_rate(prob::ForceBalanceProblem, alg::Forward; inertia::Perpendic marginal.free_boundary === nothing && error("growth_rate needs the free-boundary energies: set vac_flag=true") push!(history, (gamma_sq=0.0, lambda=lambda_min(marginal), nzero=marginal.diagnostics.nzero)) verbose && @info @sprintf("growth_rate: marginal λ_min = %.6e (nzero = %d)", history[end].lambda, history[end].nzero) - history[end].lambda < 0 || return GrowthRateResult(0.0, 0.0, marginal, nothing, history) + history[end].lambda < 0 || return IdealGrowthRateResult(0.0, 0.0, marginal, nothing, history) kwargs = copy(prob.ctrl_kwargs) haskey(kwargs, :singfac_min) || (kwargs[:singfac_min] = 0.0) @@ -157,7 +157,7 @@ function growth_rate(prob::ForceBalanceProblem, alg::Forward; inertia::Perpendic lam = lambda_at(hi) while lam > 0 hi /= 10 - hi < gamma_sq_min && return GrowthRateResult(sqrt(hi), hi, marginal, last[], history) + hi < gamma_sq_min && return IdealGrowthRateResult(sqrt(hi), hi, marginal, last[], history) lam = lambda_at(hi) end lo = hi @@ -186,5 +186,5 @@ function growth_rate(prob::ForceBalanceProblem, alg::Forward; inertia::Perpendic end hi - lo <= rtol * hi && break end - return GrowthRateResult(sqrt(gamma_sq), gamma_sq, marginal, last[], history) + return IdealGrowthRateResult(sqrt(gamma_sq), gamma_sq, marginal, last[], history) end diff --git a/src/GeneralizedPerturbedEquilibrium.jl b/src/GeneralizedPerturbedEquilibrium.jl index b753a8e68..3c53e95f3 100755 --- a/src/GeneralizedPerturbedEquilibrium.jl +++ b/src/GeneralizedPerturbedEquilibrium.jl @@ -71,7 +71,8 @@ include("Rerun.jl") # Scripting-API surface: the integrator selectors, the published result, the equilibrium # constructor and the forcing description, re-exported so a user needs one `using`. -using .ForceFreeStates: AbstractIntegrator, Forward, Riccati, Galerkin, TestKineticMatrices, KineticMatrices, ForceFreeStatesResult, PerpendicularInertia, GrowthRateResult, growth_rate +using .ForceFreeStates: + AbstractIntegrator, Forward, Riccati, Galerkin, TestKineticMatrices, KineticMatrices, ForceFreeStatesResult, PerpendicularInertia, IdealGrowthRateResult, growth_rate using .ForceFreeStates: MatchProblem using .InnerLayer: GGJModel, SLAYERModel using .Tearing: TearingProblem, TearingConfig, TearingResult, CriticalResonantFieldControl @@ -305,7 +306,7 @@ export main, write_imas export solve, plasma_response, perturbed_equilibrium, kinetic_torque, error_fields export KineticProfileData, IonSpecies, attach_kinetic_profiles! export PlasmaEquilibrium, ForceBalanceProblem, Forward, Riccati, Galerkin, TestKineticMatrices, KineticMatrices, kinetic_matrices, ForceFreeStatesResult, RMPField, PlasmaResponse -export PerpendicularInertia, GrowthRateResult, growth_rate +export PerpendicularInertia, IdealGrowthRateResult, growth_rate export MatchProblem, TearingProblem, TearingConfig, TearingResult, GGJModel, SLAYERModel export LocalStabilityProblem, Ballooning, LocalStabilityResult export KineticForcesControl, KineticForcesResult, ErrorFieldsResult From 84374cf108b665084c776dffc7a39037b81e4e36 Mon Sep 17 00:00:00 2001 From: d-burg Date: Sat, 7 Nov 2026 18:36:05 -0600 Subject: [PATCH 06/12] =?UTF-8?q?ForceFreeStates=20-=20FEATURE=20-=20Deriv?= =?UTF-8?q?e=20the=20inertia=20density=20and=20=CF=89=5F*pi=20from=20the?= =?UTF-8?q?=20kinetic=20profiles?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit PerpendicularInertia(equil) builds ρ(ψ_N) = μ m_p n_e from the attached kinetic profiles with the layer's mass per electron; ion_diamagnetic_frequency gives ω_*pi(ψ_N) for a toroidal n and growth_rate reports its peak over the domain for the EPED1 γ > ω_*pi/2 criterion. GrowthRateConfig carries the settings of the coming [GrowthRate] deck section, and growth_rate reuses a marginal solve. Co-Authored-By: Claude Fable 5.1 --- src/ForceFreeStates/Inertia.jl | 93 +++++++++++++++++++++++++++++++--- test/runtests_inertia.jl | 39 ++++++++++++++ 2 files changed, 126 insertions(+), 6 deletions(-) diff --git a/src/ForceFreeStates/Inertia.jl b/src/ForceFreeStates/Inertia.jl index 20cabfa51..f64e21c8f 100644 --- a/src/ForceFreeStates/Inertia.jl +++ b/src/ForceFreeStates/Inertia.jl @@ -24,6 +24,46 @@ end _rho_at(src::PerpendicularInertia, psi::Real) = src.rho isa Function ? Float64(src.rho(psi)) : src.rho +""" + PerpendicularInertia(equil::Equilibrium.PlasmaEquilibrium) + +The inertia of the equilibrium's attached kinetic profiles: ρ(ψ_N) = μ m_p n_e(ψ_N), with μ the mean ion mass per +electron of the main ions and the impurity at the local Z_eff, the density the SLAYER layer uses. Errors when no +kinetic profiles are attached (set `kinetic_file` or call `attach_kinetic_profiles!`). +""" +function PerpendicularInertia(equil::Equilibrium.PlasmaEquilibrium) + data = InnerLayer._layer_kinetic_data(equil, "mass density ρ") + loaded = InnerLayer.SLAYER._slayer_profiles(data, nothing) + mu = InnerLayer.SLAYER._layer_mass_per_electron(data, InnerLayer.SLAYER._layer_zeff(loaded, data)) + n_e = loaded.profiles.n_e + return PerpendicularInertia(; rho=psi -> (mu isa Real ? mu : mu(psi)) * Utilities.M_P * n_e(psi)) +end + +""" + ion_diamagnetic_frequency(equil::Equilibrium.PlasmaEquilibrium, n::Integer; z_i=1) -> ψ_N -> ω_*pi + +The ion diamagnetic frequency in rad/s for toroidal mode number `n`, ω_*pi = n (dp_i/dΨ)/(Z_i e n_i) with Ψ the +poloidal flux per radian, from the attached kinetic profiles with n_i ∝ n_e and p_i = n_i T_i, the layer +convention (positive for the usual outward-decreasing profiles). The EPED1 criterion counts a mode as unstable +when γ > ω_*pi/2 at the pedestal. +""" +function ion_diamagnetic_frequency(equil::Equilibrium.PlasmaEquilibrium, n::Integer; z_i::Real=1) + data = InnerLayer._layer_kinetic_data(equil, "ion diamagnetic frequency ω_*pi") + data.T_i === nothing && throw(ArgumentError("ω_*pi needs the ion temperature, and the kinetic data carries no T_i")) + p = InnerLayer.SLAYER._slayer_profiles(data, nothing).profiles + scale = -n / (z_i * equil.psio) + return psi -> scale * (p.T_i(psi) * p.n_e(psi; deriv=DerivOp(1)) / p.n_e(psi) + p.T_i(psi; deriv=DerivOp(1))) +end + +# Peak |ω_*pi| over the integration domain of `r` at its lowest n; NaN without kinetic profiles that carry T_i. +function _peak_omega_star_pi(r::ForceFreeStatesResult) + data = r.equil.kinetic + (data === nothing || data.T_i === nothing) && return NaN + omega = ion_diamagnetic_frequency(r.equil, r.nlow) + psis = filter(x -> r.psilow <= x <= r.psilim, r.metric.xs) + return isempty(psis) ? NaN : maximum(abs ∘ omega, psis) +end + """ inertia_kinetic_matrices(src::PerpendicularInertia, equil, intr, metric) -> (kw_flat, kt_flat) @@ -89,6 +129,8 @@ The outcome of [`growth_rate`](@ref). - `gamma::Float64` - Growth rate in s⁻¹ of the least-stable free-boundary mode; 0 when the plasma is marginally stable or stable. - `gamma_sq::Float64` - γ² at which the least-stable free-boundary energy eigenvalue crosses zero. + - `omega_star_pi::Float64` - Peak ion diamagnetic frequency |ω_*pi| in rad/s over the integration domain at the problem's lowest n, from the + equilibrium's kinetic profiles (NaN when none carry T_i); the EPED1 criterion counts the mode as unstable when `gamma > omega_star_pi/2`. - `marginal::ForceFreeStatesResult` - The ideal (γ² = 0) solve, whose least-stable energy decides stability. - `result::Union{Nothing,ForceFreeStatesResult}` - The solve at `gamma_sq`, carrying the eigenmode; `nothing` when stable. - `history::Vector` - Every (γ², λ_min, nzero) evaluated, in order, with λ_min the least-stable free-boundary energy at that γ². @@ -96,6 +138,7 @@ The outcome of [`growth_rate`](@ref). struct IdealGrowthRateResult gamma::Float64 gamma_sq::Float64 + omega_star_pi::Float64 marginal::ForceFreeStatesResult result::Union{Nothing,ForceFreeStatesResult} history::Vector{@NamedTuple{gamma_sq::Float64, lambda::Float64, nzero::Int}} @@ -103,7 +146,8 @@ end """ growth_rate(prob::ForceBalanceProblem, alg::Forward; inertia, gamma_sq_init=1e9, gamma_sq_min=1e4, gamma_sq_max=1e16, - rtol=1e-3, ftol=1e-6, maxiter=30, verbose=false) -> IdealGrowthRateResult + rtol=1e-3, ftol=1e-6, maxiter=30, marginal=nothing, verbose=false) -> IdealGrowthRateResult + growth_rate(prob::ForceBalanceProblem, alg::Forward, config::GrowthRateConfig; marginal=nothing, verbose=false) Growth rate of the least-stable free-boundary mode of `prob` with the perpendicular inertia `inertia` (a [`PerpendicularInertia`](@ref)), ELITE's eigenvalue problem solved by shooting: the ideal solve @@ -115,6 +159,7 @@ bracketed secant (Illinois) iteration that stops when |λ_min| falls below `ftol one force-balance solve. `gamma_sq_min` is the floor below which the mode counts as marginal, `gamma_sq_max` the ceiling beyond which bracketing stops with an error. +A `marginal` result of the same problem (an earlier ideal `solve`) is reused instead of solving it again. The problem's `singfac_min` is set to 0 for the inertial solves unless given: with inertia the rational surfaces are regular, so the integration runs as a single chunk. A fixed-boundary mode faster than the free-boundary one shows as `nzero > 0` at the root in `history`; this first version @@ -122,15 +167,16 @@ reports the free-boundary root and leaves the interior count to the caller. """ function growth_rate(prob::ForceBalanceProblem, alg::Forward; inertia::PerpendicularInertia, gamma_sq_init::Real=1e9, gamma_sq_min::Real=1e4, gamma_sq_max::Real=1e16, rtol::Real=1e-3, ftol::Real=1e-6, maxiter::Int=30, - verbose::Bool=false) + marginal::Union{Nothing,ForceFreeStatesResult}=nothing, verbose::Bool=false) history = @NamedTuple{gamma_sq::Float64, lambda::Float64, nzero::Int}[] lambda_min(r) = minimum(real, r.free_boundary.et) - marginal = solve(prob, alg; verbose) + marginal = marginal === nothing ? solve(prob, alg; verbose) : marginal marginal.free_boundary === nothing && error("growth_rate needs the free-boundary energies: set vac_flag=true") + omega_star_pi = _peak_omega_star_pi(marginal) push!(history, (gamma_sq=0.0, lambda=lambda_min(marginal), nzero=marginal.diagnostics.nzero)) verbose && @info @sprintf("growth_rate: marginal λ_min = %.6e (nzero = %d)", history[end].lambda, history[end].nzero) - history[end].lambda < 0 || return IdealGrowthRateResult(0.0, 0.0, marginal, nothing, history) + history[end].lambda < 0 || return IdealGrowthRateResult(0.0, 0.0, omega_star_pi, marginal, nothing, history) kwargs = copy(prob.ctrl_kwargs) haskey(kwargs, :singfac_min) || (kwargs[:singfac_min] = 0.0) @@ -157,7 +203,7 @@ function growth_rate(prob::ForceBalanceProblem, alg::Forward; inertia::Perpendic lam = lambda_at(hi) while lam > 0 hi /= 10 - hi < gamma_sq_min && return IdealGrowthRateResult(sqrt(hi), hi, marginal, last[], history) + hi < gamma_sq_min && return IdealGrowthRateResult(sqrt(hi), hi, omega_star_pi, marginal, last[], history) lam = lambda_at(hi) end lo = hi @@ -186,5 +232,40 @@ function growth_rate(prob::ForceBalanceProblem, alg::Forward; inertia::Perpendic end hi - lo <= rtol * hi && break end - return IdealGrowthRateResult(sqrt(gamma_sq), gamma_sq, marginal, last[], history) + return IdealGrowthRateResult(sqrt(gamma_sq), gamma_sq, omega_star_pi, marginal, last[], history) +end + +""" + GrowthRateConfig(; rho=NaN, gamma_sq_init=1e9, gamma_sq_min=1e4, gamma_sq_max=1e16, rtol=1e-3, ftol=1e-6, maxiter=30) + +The `[GrowthRate]` deck section: its presence runs [`growth_rate`](@ref) after the force-balance stage, and the +same settings drive the three-argument `growth_rate` method. + +## Fields + + - `rho` - Uniform mass density in kg/m³; NaN (the default) derives ρ(ψ_N) from the equilibrium's kinetic profiles. + - `gamma_sq_init` - First γ² in s⁻² tried when bracketing the root. + - `gamma_sq_min` - γ² floor in s⁻² below which the mode counts as marginal. + - `gamma_sq_max` - γ² ceiling in s⁻² beyond which bracketing stops with an error. + - `rtol` - Relative width of the γ² bracket at which the root iteration stops. + - `ftol` - Fraction of the marginal |λ_min| below which |λ_min| counts as zero. + - `maxiter` - Maximum root-iteration steps after bracketing. +""" +@kwdef struct GrowthRateConfig + rho::Float64 = NaN + gamma_sq_init::Float64 = 1e9 + gamma_sq_min::Float64 = 1e4 + gamma_sq_max::Float64 = 1e16 + rtol::Float64 = 1e-3 + ftol::Float64 = 1e-6 + maxiter::Int = 30 +end + +# The inertia a deck section asks for: uniform, or from the equilibrium's kinetic profiles. +inertia_source(config::GrowthRateConfig, equil::Equilibrium.PlasmaEquilibrium) = + isnan(config.rho) ? PerpendicularInertia(equil) : PerpendicularInertia(; rho=config.rho) + +function growth_rate(prob::ForceBalanceProblem, alg::Forward, config::GrowthRateConfig; marginal::Union{Nothing,ForceFreeStatesResult}=nothing, verbose::Bool=false) + return growth_rate(prob, alg; inertia=inertia_source(config, prob.equil), config.gamma_sq_init, config.gamma_sq_min, config.gamma_sq_max, + config.rtol, config.ftol, config.maxiter, marginal, verbose) end diff --git a/test/runtests_inertia.jl b/test/runtests_inertia.jl index 365392aba..33196f82a 100644 --- a/test/runtests_inertia.jl +++ b/test/runtests_inertia.jl @@ -98,4 +98,43 @@ using LinearAlgebra @test length(gr.history) == 1 @test gr.history[1].lambda ≈ lambda_min(ideal) rtol = 1e-12 end + + @testset "rho and omega_star_pi from kinetic profiles" begin + ktemplate = joinpath(@__DIR__, "test_data", "regression_tj_kinetic_calculated") + kdeck = TOML.parsefile(joinpath(ktemplate, "gpec.toml")) + kequil = GPEC.Equilibrium.PlasmaEquilibrium(GPEC.equilibrium_config(kdeck, ktemplate), GPEC.analytic_config(kdeck)) + @test kequil.kinetic !== nothing + src = FFS.PerpendicularInertia(kequil) + # The density the SLAYER layer uses at a rational surface + @test FFS._rho_at(src, 0.5) ≈ only(GPEC.InnerLayer.layer_rho(kequil, [(psifac=0.5,)])) rtol = 1e-12 + @test FFS._rho_at(src, 0.3) > 0 + # ω_*pi = n (dp_i/dΨ)/(Z e n_i) with p_i = n_e T_i, against a finite difference of the profile splines + p = GPEC.InnerLayer.SLAYER._slayer_profiles(kequil.kinetic, nothing).profiles + p_i(ψ) = p.n_e(ψ) * p.T_i(ψ) + h = 1e-4 + for n in (1, 3) + omega = FFS.ion_diamagnetic_frequency(kequil, n) + fd = -n / kequil.psio * (p_i(0.5 + h) - p_i(0.5 - h)) / (2h) / p.n_e(0.5) + @test omega(0.5) ≈ fd rtol = 1e-6 + end + @test FFS.ion_diamagnetic_frequency(kequil, 2)(0.5) ≈ 2 * FFS.ion_diamagnetic_frequency(kequil, 1)(0.5) + @test FFS.ion_diamagnetic_frequency(kequil, 1)(0.5) > 0 # outward-decreasing p_i: the layer sign convention + @test_throws ArgumentError FFS.ion_diamagnetic_frequency(equil, 1) # the ideal fixture carries no kinetic data + @test_throws ArgumentError FFS.PerpendicularInertia(equil) + + # The kinetic TJ fixture is ideal-stable too: γ = 0 with the domain's peak ω_*pi reported + kkw = Dict{Symbol,Any}(q.first => q.second for q in GPEC.force_balance_keywords(kdeck, ktemplate) if !(q.first in (:nn, :wall, :kinetic, :kinetic_factor))) + kkw[:vac_flag] = true + kprob = ForceBalanceProblem(kequil; nn=1, wall=GPEC.Vacuum.parse_wall(get(kdeck, "Wall", nothing)), kkw...) + gr = FFS.growth_rate(kprob, forward; inertia=src) + @test gr.gamma == 0 && isfinite(gr.omega_star_pi) && gr.omega_star_pi > 0 + omega1 = FFS.ion_diamagnetic_frequency(kequil, 1) + @test gr.omega_star_pi ≈ maximum(abs ∘ omega1, filter(x -> gr.marginal.psilow <= x <= gr.marginal.psilim, gr.marginal.metric.xs)) + # The config form reuses a marginal solve and derives the same inertia + cfg = FFS.GrowthRateConfig() + gr2 = FFS.growth_rate(kprob, forward, cfg; marginal=gr.marginal) + @test gr2.gamma == gr.gamma && gr2.omega_star_pi == gr.omega_star_pi && gr2.marginal === gr.marginal + @test FFS._rho_at(FFS.inertia_source(cfg, kequil), 0.4) == FFS._rho_at(src, 0.4) + @test FFS._rho_at(FFS.inertia_source(FFS.GrowthRateConfig(; rho=2.0), kequil), 0.4) == 2.0 + end end From 6e37bc4e7f0f78b68b4c851387e500f3bf050f64 Mon Sep 17 00:00:00 2001 From: d-burg Date: Sat, 7 Nov 2026 18:57:19 -0600 Subject: [PATCH 07/12] ForceFreeStates - FEATURE - Add the [GrowthRate] deck section and its HDF5 output A [GrowthRate] section (keys = GrowthRateConfig fields; rho NaN derives the density from the kinetic profiles) runs growth_rate after the force-free-states stage, reusing that stage's ideal solve as the marginal one; it needs vac_flag and the Forward integrator and cannot share a deck with [Match]. The result is the growth_rate slot of main's tuple and ForceFreeStates/GrowthRate/ in gpec.h5 (gamma, gamma_sq, omega_star_pi, rho on psi, and the root-search History), with the inputs recorded under Input/. Co-Authored-By: Claude Fable 5.1 --- docs/src/api.md | 6 +- docs/src/workflow.md | 1 + src/Deck.jl | 26 ++++++++- src/ForceFreeStates/ForceFreeStates.jl | 1 + src/ForceFreeStates/Inertia.jl | 79 ++++++++++++++------------ src/ForceFreeStates/InertiaOutput.jl | 42 ++++++++++++++ src/GeneralizedPerturbedEquilibrium.jl | 19 +++++-- src/Utilities/HDF5Schema.jl | 1 + 8 files changed, 128 insertions(+), 47 deletions(-) create mode 100644 src/ForceFreeStates/InertiaOutput.jl diff --git a/docs/src/api.md b/docs/src/api.md index 06a117a30..51a401514 100644 --- a/docs/src/api.md +++ b/docs/src/api.md @@ -21,7 +21,7 @@ foreach(r -> save_to_hdf5("gpec.h5", r), (pe, kf, ef, tear)) ``` `solve(eq, alg; kwargs...)` is sugar for the problem form `solve(ForceBalanceProblem(eq; kwargs...), alg)`: the problem names what is solved (the force-balance equation on this equilibrium, with this mode range, wall and closure) and the algorithm names how. -`main` translates each deck section into the problem, algorithm or control it names, calls these functions for every section present, saves each result into `[Run] output_file`, and returns `(; equil, local_stability, ffs, match, resp, pe, kinetic, error_fields, tearing)`, each `nothing` when its section is absent. +`main` translates each deck section into the problem, algorithm or control it names, calls these functions for every section present, saves each result into `[Run] output_file`, and returns `(; equil, local_stability, ffs, match, growth_rate, resp, pe, kinetic, error_fields, tearing)`, each `nothing` when its section is absent. | Stage | Deck section | Call | Result | |---|---|---|---| @@ -29,12 +29,16 @@ foreach(r -> save_to_hdf5("gpec.h5", r), (pe, kf, ef, tear)) | ForceFreeStates | `[ForceFreeStates]` + `[Forward]`, `[Riccati]` or `[Galerkin]` | `solve(ForceBalanceProblem(eq; ...), Forward() \| Riccati() \| Galerkin())` | `ForceFreeStatesResult` | | LocalStability | `[LocalStability]` + `[Ballooning]` | `solve(LocalStabilityProblem(eq; ...), Ballooning(; ...))` | `LocalStabilityResult` | | Match | `[Match]` + `[GGJ]` | `solve(MatchProblem(ffs; ...), GGJModel())` | `ForceFreeStatesResult` with `closure == :matched` | +| GrowthRate | `[GrowthRate]` | `growth_rate(ForceBalanceProblem(eq; ...), Forward(), GrowthRateConfig(; ...); marginal=ffs)` | `IdealGrowthRateResult` | | Drive | `[ForcingTerms]` | `RMPField(ft_ctrl)` (a `ForcingTermsControl`), `RMPField(file)`, `RMPField(modes)` | `RMPField` | | PerturbedEquilibrium | `[PerturbedEquilibrium]` | `plasma_response(ffs; ...)`, then `perturbed_equilibrium(resp, rmp; ...)` | `PlasmaResponse`, `PerturbedEquilibriumState` | | KineticForces | `[KineticForces]` | `kinetic_torque(pe, KineticForcesControl(; ...))` | `KineticForcesResult` | | ErrorFields | `[ErrorFields]` | `error_fields(resp, pe.coil_sets, CoilConfig(...); ...)` | `ErrorFieldsResult` | | Tearing | `[Tearing]` + `[SLAYER]` or `[GGJ]` | `solve(TearingProblem(ffs; ...), SLAYERModel() \| GGJModel(); config=TearingConfig(...))` | `TearingResult` | +`[GrowthRate]` runs the ideal growth-rate solve after the force-free-states stage; its keys are the `GrowthRateConfig` fields, and its `rho`, left at NaN, derives the mass density from the equilibrium's kinetic profiles. +It needs `[ForceFreeStates] vac_flag = true` and `integrator = "Forward"`, and cannot share a deck with `[Match]`, since its marginal solve is the plain ideal one. + The [Workflow](workflow.md) page describes what each stage computes; the module pages document the keywords, which are the deck keys of the matching section. ## Inspecting a result diff --git a/docs/src/workflow.md b/docs/src/workflow.md index 38b4cdd8b..9a280c199 100644 --- a/docs/src/workflow.md +++ b/docs/src/workflow.md @@ -21,6 +21,7 @@ Every stage after `Equilibrium` reads the force-free states, and `KineticForces` |---|---|---| | Is the equilibrium read and mapped correctly; what are q and β_N? | `[Equilibrium]` (+ `[ForceFreeStates]`, which triggers the write) | `Equilibrium/q_95`, `Equilibrium/beta_N`, `Equilibrium/Profiles/` | | Is the plasma ideal-kink stable, with or without a wall? | `[ForceFreeStates]` (+ `[Wall]`) | `ForceFreeStates/FreeBoundaryStability/W_freeboundary`, `…/eigenmode_energies` | +| How fast does the least-stable free-boundary ideal mode grow? | `[GrowthRate]` (+ `[ForceFreeStates] vac_flag = true`, Forward) | `ForceFreeStates/GrowthRate/` | | Is it Mercier, resistive-interchange or ballooning unstable? | `[LocalStability]` (+ `[Ballooning]`) | `LocalStability/D_I`, `LocalStability/D_R`, `LocalStability/ballooning_Delta_prime` | | What is the tearing stability index Δ′ at each rational surface? | `[ForceFreeStates] integrator = "Riccati"` or `"Galerkin"` (+ `[Riccati]` or `[Galerkin]`) | `SingularSurfaces/Delta_prime_matrix` | | How does the plasma respond to an RMP or error field; do islands overlap? | `[ForcingTerms]` + `[PerturbedEquilibrium]` | `PerturbedEquilibrium/SingularCoupling/resonant_area_weighted_field`, `…/island_half_width`, `…/chirikov_parameter` | diff --git a/src/Deck.jl b/src/Deck.jl index ea0f4b17c..c08cf05cb 100644 --- a/src/Deck.jl +++ b/src/Deck.jl @@ -78,7 +78,7 @@ const _RETIRED_SECTIONS = Dict("DEBUG" => "[DEBUG] is retired: output_benchmark_ # Every section a deck may hold, besides the analytic-equilibrium sections. const _DECK_SECTIONS = ("Run", "Equilibrium", "Wall", "LocalStability", "Ballooning", "ForceFreeStates", "Forward", "Riccati", - "Galerkin", "Match", "GGJ", "Tearing", "SLAYER", "ForcingTerms", "PerturbedEquilibrium", "KineticForces", "ErrorFields") + "Galerkin", "Match", "GGJ", "Tearing", "SLAYER", "GrowthRate", "ForcingTerms", "PerturbedEquilibrium", "KineticForces", "ErrorFields") const _INTEGRATORS = ("Forward" => Forward, "Riccati" => Riccati, "Galerkin" => Galerkin) # Stop the run when a parsed gpec.toml section sets keys that no longer exist, so the deck is updated. @@ -377,6 +377,16 @@ function tearing_from_deck(inputs::AbstractDict) return (; problem, model=inner_model == "SLAYER" ? slayer_model(inputs) : ggj_model(inputs), config) end +""" + growth_rate_from_deck(inputs) -> Union{Nothing,GrowthRateConfig} + +The `GrowthRateConfig` of `[GrowthRate]`; `nothing` without that section. +""" +function growth_rate_from_deck(inputs::AbstractDict) + haskey(inputs, "GrowthRate") || return nothing + return GrowthRateConfig(; section_keywords(inputs["GrowthRate"], "GrowthRate", _struct_fields(GrowthRateConfig))...) +end + """ perturbed_equilibrium_from_deck(inputs) -> Union{Nothing,NamedTuple} @@ -447,7 +457,7 @@ function parse_deck(inputs::AbstractDict, path::AbstractString) end has_ffs = haskey(inputs, "ForceFreeStates") if !has_ffs - for section in ("Wall", "Forward", "Riccati", "Galerkin", "Match", "GGJ", "Tearing", "SLAYER", "LocalStability", "Ballooning", + for section in ("Wall", "Forward", "Riccati", "Galerkin", "Match", "GGJ", "Tearing", "SLAYER", "GrowthRate", "LocalStability", "Ballooning", "ForcingTerms", "PerturbedEquilibrium", "KineticForces", "ErrorFields") haskey(inputs, section) && error("[$section] is set but there is no [ForceFreeStates] section, so the run stops after the equilibrium") end @@ -462,6 +472,14 @@ function parse_deck(inputs::AbstractDict, path::AbstractString) error("[KineticForces] is set but neither [PerturbedEquilibrium] (the torque) nor [ForceFreeStates] kinetic = \"KineticForces\" uses it") (haskey(inputs, "KineticForces") || kinetic_source) && isempty(get(get(inputs, "Equilibrium", Dict{String,Any}()), "kinetic_file", "")) && error("this run needs kinetic profiles ([KineticForces], or kinetic = \"KineticForces\"): set kinetic_file in [Equilibrium]") + growth_rate = growth_rate_from_deck(inputs) + # The growth-rate solve needs the plain ideal solve as its marginal one and the free-boundary energy of a Forward solve. + if growth_rate !== nothing + (; keywords, algorithm) = force_balance + haskey(inputs, "Match") && error("[GrowthRate] needs the ideal solve as its marginal one: remove the [Match] section") + get(keywords, :vac_flag, false) || error("[GrowthRate] needs the free-boundary energy: set [ForceFreeStates] vac_flag = true") + algorithm isa Forward || error("[GrowthRate] needs integrator = \"Forward\" in [ForceFreeStates]") + end match = match_from_deck(inputs) # The MatchProblem guards, checked before the force-free-states solve rather than after it. if match !== nothing @@ -483,7 +501,7 @@ function parse_deck(inputs::AbstractDict, path::AbstractString) error("[GGJ] is set but neither [Match] nor [Tearing] inner_model = \"GGJ\" uses it") haskey(inputs, "SLAYER") && tearing_model != "SLAYER" && !(has_ffs && get(inputs["ForceFreeStates"], "psilim_from_layer_overlap", false) === true) && error("[SLAYER] is set but neither [Tearing] inner_model = \"SLAYER\" nor [ForceFreeStates] psilim_from_layer_overlap = true uses it") - return (; run=run_settings(inputs, path), force_balance, local_stability=local_stability_from_deck(inputs), match, tearing, + return (; run=run_settings(inputs, path), force_balance, local_stability=local_stability_from_deck(inputs), match, tearing, growth_rate, kinetic_forces=haskey(inputs, "KineticForces") ? kinetic_forces_control(inputs) : nothing, perturbed_equilibrium, forcing, error_fields=error_fields_from_deck(inputs, path)) end @@ -578,6 +596,7 @@ function deck_inputs(deck::NamedTuple, config::Equilibrium.EquilibriumConfig, an merge!(inputs, _tearing_sections(filter(p -> p.first != :critical_resonant_field, problem), get(problem, :critical_resonant_field, nothing), model, config)) end + deck.growth_rate === nothing || (inputs["GrowthRate"] = _deck_table(deck.growth_rate)) deck.kinetic_forces === nothing || merge!(inputs, deck_sections(deck.kinetic_forces)) pe = deck.perturbed_equilibrium pe === nothing || (inputs["PerturbedEquilibrium"] = _deck_table(merge(pe.response, pe.drive))) @@ -625,6 +644,7 @@ deck_sections(ls::LocalStabilityResult) = Dict{String,Any}("LocalStability" => _ deck_sections(m::ForceFreeStates.MatchResult) = merge(deck_sections(m.model), Dict{String,Any}("Match" => merge(_deck_table(m.problem, (:ffs,)), Dict{String,Any}("inner_model" => "GGJ")))) deck_sections(kf::KineticForcesResult) = deck_sections(kf.control) +deck_sections(g::IdealGrowthRateResult) = g.config === nothing ? Dict{String,Any}() : Dict{String,Any}("GrowthRate" => _deck_table(g.config)) deck_sections(t::TearingResult) = t.model === nothing ? Dict{String,Any}() : _tearing_sections(t.problem, t.problem.critical_resonant_field, t.model, t.config) deck_sections(ef::ErrorFieldsResult) = _error_fields_sections(ef.settings, ef.tolerances === nothing ? "" : _RECORDED_TOLERANCES) diff --git a/src/ForceFreeStates/ForceFreeStates.jl b/src/ForceFreeStates/ForceFreeStates.jl index e2f9bc4ee..7c08fe8ea 100644 --- a/src/ForceFreeStates/ForceFreeStates.jl +++ b/src/ForceFreeStates/ForceFreeStates.jl @@ -73,6 +73,7 @@ include("Output.jl") include("Run.jl") include("ForceBalance.jl") include("Inertia.jl") +include("InertiaOutput.jl") # Small tolerance and root-finder iteration cap used throughout ForceFreeStates global eps = 1e-10 diff --git a/src/ForceFreeStates/Inertia.jl b/src/ForceFreeStates/Inertia.jl index f64e21c8f..11f67cbbf 100644 --- a/src/ForceFreeStates/Inertia.jl +++ b/src/ForceFreeStates/Inertia.jl @@ -120,6 +120,36 @@ function inertia_kinetic_matrices(src::PerpendicularInertia, equil::Equilibrium. return kw_flat, kt_flat end +""" + GrowthRateConfig(; rho=NaN, gamma_sq_init=1e9, gamma_sq_min=1e4, gamma_sq_max=1e16, rtol=1e-3, ftol=1e-6, maxiter=30) + +The `[GrowthRate]` deck section: its presence runs [`growth_rate`](@ref) after the force-balance stage, and the +same settings drive the three-argument `growth_rate` method. + +## Fields + + - `rho` - Uniform mass density in kg/m³; NaN (the default) derives ρ(ψ_N) from the equilibrium's kinetic profiles. + - `gamma_sq_init` - First γ² in s⁻² tried when bracketing the root. + - `gamma_sq_min` - γ² floor in s⁻² below which the mode counts as marginal. + - `gamma_sq_max` - γ² ceiling in s⁻² beyond which bracketing stops with an error. + - `rtol` - Relative width of the γ² bracket at which the root iteration stops. + - `ftol` - Fraction of the marginal |λ_min| below which |λ_min| counts as zero. + - `maxiter` - Maximum root-iteration steps after bracketing. +""" +@kwdef struct GrowthRateConfig + rho::Float64 = NaN + gamma_sq_init::Float64 = 1e9 + gamma_sq_min::Float64 = 1e4 + gamma_sq_max::Float64 = 1e16 + rtol::Float64 = 1e-3 + ftol::Float64 = 1e-6 + maxiter::Int = 30 +end + +# The inertia a deck section asks for: uniform, or from the equilibrium's kinetic profiles. +inertia_source(config::GrowthRateConfig, equil::Equilibrium.PlasmaEquilibrium) = + isnan(config.rho) ? PerpendicularInertia(equil) : PerpendicularInertia(; rho=config.rho) + """ IdealGrowthRateResult @@ -129,6 +159,8 @@ The outcome of [`growth_rate`](@ref). - `gamma::Float64` - Growth rate in s⁻¹ of the least-stable free-boundary mode; 0 when the plasma is marginally stable or stable. - `gamma_sq::Float64` - γ² at which the least-stable free-boundary energy eigenvalue crosses zero. + - `inertia::PerpendicularInertia` - The inertia source the solves used. + - `config::Union{Nothing,GrowthRateConfig}` - The deck settings when the solve came from a `[GrowthRate]` section or the three-argument method; `nothing` otherwise. - `omega_star_pi::Float64` - Peak ion diamagnetic frequency |ω_*pi| in rad/s over the integration domain at the problem's lowest n, from the equilibrium's kinetic profiles (NaN when none carry T_i); the EPED1 criterion counts the mode as unstable when `gamma > omega_star_pi/2`. - `marginal::ForceFreeStatesResult` - The ideal (γ² = 0) solve, whose least-stable energy decides stability. @@ -138,6 +170,8 @@ The outcome of [`growth_rate`](@ref). struct IdealGrowthRateResult gamma::Float64 gamma_sq::Float64 + inertia::PerpendicularInertia + config::Union{Nothing,GrowthRateConfig} omega_star_pi::Float64 marginal::ForceFreeStatesResult result::Union{Nothing,ForceFreeStatesResult} @@ -146,7 +180,7 @@ end """ growth_rate(prob::ForceBalanceProblem, alg::Forward; inertia, gamma_sq_init=1e9, gamma_sq_min=1e4, gamma_sq_max=1e16, - rtol=1e-3, ftol=1e-6, maxiter=30, marginal=nothing, verbose=false) -> IdealGrowthRateResult + rtol=1e-3, ftol=1e-6, maxiter=30, marginal=nothing, config=nothing, verbose=false) -> IdealGrowthRateResult growth_rate(prob::ForceBalanceProblem, alg::Forward, config::GrowthRateConfig; marginal=nothing, verbose=false) Growth rate of the least-stable free-boundary mode of `prob` with the perpendicular inertia `inertia` @@ -159,7 +193,8 @@ bracketed secant (Illinois) iteration that stops when |λ_min| falls below `ftol one force-balance solve. `gamma_sq_min` is the floor below which the mode counts as marginal, `gamma_sq_max` the ceiling beyond which bracketing stops with an error. -A `marginal` result of the same problem (an earlier ideal `solve`) is reused instead of solving it again. +A `marginal` result of the same problem (an earlier ideal `solve`) is reused instead of solving it again; `config` +is recorded on the result and is what the three-argument method passes. The problem's `singfac_min` is set to 0 for the inertial solves unless given: with inertia the rational surfaces are regular, so the integration runs as a single chunk. A fixed-boundary mode faster than the free-boundary one shows as `nzero > 0` at the root in `history`; this first version @@ -167,7 +202,7 @@ reports the free-boundary root and leaves the interior count to the caller. """ function growth_rate(prob::ForceBalanceProblem, alg::Forward; inertia::PerpendicularInertia, gamma_sq_init::Real=1e9, gamma_sq_min::Real=1e4, gamma_sq_max::Real=1e16, rtol::Real=1e-3, ftol::Real=1e-6, maxiter::Int=30, - marginal::Union{Nothing,ForceFreeStatesResult}=nothing, verbose::Bool=false) + marginal::Union{Nothing,ForceFreeStatesResult}=nothing, config::Union{Nothing,GrowthRateConfig}=nothing, verbose::Bool=false) history = @NamedTuple{gamma_sq::Float64, lambda::Float64, nzero::Int}[] lambda_min(r) = minimum(real, r.free_boundary.et) @@ -176,7 +211,7 @@ function growth_rate(prob::ForceBalanceProblem, alg::Forward; inertia::Perpendic omega_star_pi = _peak_omega_star_pi(marginal) push!(history, (gamma_sq=0.0, lambda=lambda_min(marginal), nzero=marginal.diagnostics.nzero)) verbose && @info @sprintf("growth_rate: marginal λ_min = %.6e (nzero = %d)", history[end].lambda, history[end].nzero) - history[end].lambda < 0 || return IdealGrowthRateResult(0.0, 0.0, omega_star_pi, marginal, nothing, history) + history[end].lambda < 0 || return IdealGrowthRateResult(0.0, 0.0, inertia, config, omega_star_pi, marginal, nothing, history) kwargs = copy(prob.ctrl_kwargs) haskey(kwargs, :singfac_min) || (kwargs[:singfac_min] = 0.0) @@ -203,7 +238,7 @@ function growth_rate(prob::ForceBalanceProblem, alg::Forward; inertia::Perpendic lam = lambda_at(hi) while lam > 0 hi /= 10 - hi < gamma_sq_min && return IdealGrowthRateResult(sqrt(hi), hi, omega_star_pi, marginal, last[], history) + hi < gamma_sq_min && return IdealGrowthRateResult(sqrt(hi), hi, inertia, config, omega_star_pi, marginal, last[], history) lam = lambda_at(hi) end lo = hi @@ -232,40 +267,10 @@ function growth_rate(prob::ForceBalanceProblem, alg::Forward; inertia::Perpendic end hi - lo <= rtol * hi && break end - return IdealGrowthRateResult(sqrt(gamma_sq), gamma_sq, omega_star_pi, marginal, last[], history) + return IdealGrowthRateResult(sqrt(gamma_sq), gamma_sq, inertia, config, omega_star_pi, marginal, last[], history) end -""" - GrowthRateConfig(; rho=NaN, gamma_sq_init=1e9, gamma_sq_min=1e4, gamma_sq_max=1e16, rtol=1e-3, ftol=1e-6, maxiter=30) - -The `[GrowthRate]` deck section: its presence runs [`growth_rate`](@ref) after the force-balance stage, and the -same settings drive the three-argument `growth_rate` method. - -## Fields - - - `rho` - Uniform mass density in kg/m³; NaN (the default) derives ρ(ψ_N) from the equilibrium's kinetic profiles. - - `gamma_sq_init` - First γ² in s⁻² tried when bracketing the root. - - `gamma_sq_min` - γ² floor in s⁻² below which the mode counts as marginal. - - `gamma_sq_max` - γ² ceiling in s⁻² beyond which bracketing stops with an error. - - `rtol` - Relative width of the γ² bracket at which the root iteration stops. - - `ftol` - Fraction of the marginal |λ_min| below which |λ_min| counts as zero. - - `maxiter` - Maximum root-iteration steps after bracketing. -""" -@kwdef struct GrowthRateConfig - rho::Float64 = NaN - gamma_sq_init::Float64 = 1e9 - gamma_sq_min::Float64 = 1e4 - gamma_sq_max::Float64 = 1e16 - rtol::Float64 = 1e-3 - ftol::Float64 = 1e-6 - maxiter::Int = 30 -end - -# The inertia a deck section asks for: uniform, or from the equilibrium's kinetic profiles. -inertia_source(config::GrowthRateConfig, equil::Equilibrium.PlasmaEquilibrium) = - isnan(config.rho) ? PerpendicularInertia(equil) : PerpendicularInertia(; rho=config.rho) - function growth_rate(prob::ForceBalanceProblem, alg::Forward, config::GrowthRateConfig; marginal::Union{Nothing,ForceFreeStatesResult}=nothing, verbose::Bool=false) return growth_rate(prob, alg; inertia=inertia_source(config, prob.equil), config.gamma_sq_init, config.gamma_sq_min, config.gamma_sq_max, - config.rtol, config.ftol, config.maxiter, marginal, verbose) + config.rtol, config.ftol, config.maxiter, marginal, config, verbose) end diff --git a/src/ForceFreeStates/InertiaOutput.jl b/src/ForceFreeStates/InertiaOutput.jl new file mode 100644 index 000000000..215082671 --- /dev/null +++ b/src/ForceFreeStates/InertiaOutput.jl @@ -0,0 +1,42 @@ +# Writer for an `IdealGrowthRateResult`. Layout under the caller's parent group (inputs live only under Input/): +# +# ForceFreeStates/GrowthRate/ +# ├── gamma, gamma_sq, omega_star_pi -- growth rate, its square and the EPED1 ω_*pi reference +# ├── psi, rho -- inertia density on the marginal solve's ψ_N grid +# └── History/ -- every (γ², λ_min, nzero) evaluated by the root search + +const GROWTH_RATE_H5_ANNOTATIONS = [ + "gamma" => (; long_name="growth rate γ of the least-stable free-boundary mode with perpendicular inertia (0 when stable)", units="s^-1"), + "gamma_sq" => (; long_name="γ² at which the least-stable free-boundary energy eigenvalue crosses zero (0 when stable)", units="s^-2"), + "omega_star_pi" => (; long_name="peak ion diamagnetic frequency |ω_*pi| over the integration domain at the lowest n (NaN without kinetic T_i)", + units="rad/s"), + "psi" => (; long_name="normalized poloidal flux ψ_N grid of the marginal solve", dims=("psi",), scale="psi"), + "rho" => (; long_name="perpendicular-inertia mass density ρ", units="kg/m^3", dims=("psi",), attach=(1 => "psi",)), + "History/gamma_sq" => (; long_name="γ² evaluated by the growth-rate root search, in order", units="s^-2", dims=("iteration",)), + "History/lambda_min" => (; long_name="power-normalized least-stable free-boundary energy at each evaluated γ² (per unit ⟨|ξ|²⟩)", dims=("iteration",)), + "History/nzero" => (; long_name="number of zero crossings of the fixed-boundary critical eigenvalue at each evaluated γ²", dims=("iteration",)) +] + +""" + save_to_hdf5(parent::Union{HDF5.File,HDF5.Group}, result::IdealGrowthRateResult) + +Write `result` into a `ForceFreeStates/GrowthRate/` subgroup of `parent`, replacing any existing one, and record its +inputs in `Input/`. The growth rate, its square and the peak ω_*pi are scalars; the inertia density is sampled on +the marginal solve's ψ_N grid; `History/` holds every (γ², λ_min, nzero) the root search evaluated. +""" +function save_to_hdf5(parent::Union{HDF5.File,HDF5.Group}, result::IdealGrowthRateResult) + ffs = haskey(parent, "ForceFreeStates") ? parent["ForceFreeStates"] : create_group(parent, "ForceFreeStates") + g = Utilities.HDF5Annotations.replace_group!(ffs, "GrowthRate") + g["gamma"] = result.gamma + g["gamma_sq"] = result.gamma_sq + g["omega_star_pi"] = result.omega_star_pi + psi = collect(Float64, result.marginal.metric.xs) + g["psi"] = psi + g["rho"] = [_rho_at(result.inertia, x) for x in psi] + g["History/gamma_sq"] = [h.gamma_sq for h in result.history] + g["History/lambda_min"] = [h.lambda for h in result.history] + g["History/nzero"] = [h.nzero for h in result.history] + Utilities.HDF5Annotations.annotate!(g, GROWTH_RATE_H5_ANNOTATIONS) + Utilities.record_inputs!(parent, result) + return g +end diff --git a/src/GeneralizedPerturbedEquilibrium.jl b/src/GeneralizedPerturbedEquilibrium.jl index 3c53e95f3..b4caff512 100755 --- a/src/GeneralizedPerturbedEquilibrium.jl +++ b/src/GeneralizedPerturbedEquilibrium.jl @@ -72,7 +72,8 @@ include("Rerun.jl") # Scripting-API surface: the integrator selectors, the published result, the equilibrium # constructor and the forcing description, re-exported so a user needs one `using`. using .ForceFreeStates: - AbstractIntegrator, Forward, Riccati, Galerkin, TestKineticMatrices, KineticMatrices, ForceFreeStatesResult, PerpendicularInertia, IdealGrowthRateResult, growth_rate + AbstractIntegrator, Forward, Riccati, Galerkin, TestKineticMatrices, KineticMatrices, ForceFreeStatesResult, PerpendicularInertia, IdealGrowthRateResult, growth_rate, + GrowthRateConfig, ion_diamagnetic_frequency using .ForceFreeStates: MatchProblem using .InnerLayer: GGJModel, SLAYERModel using .Tearing: TearingProblem, TearingConfig, TearingResult, CriticalResonantFieldControl @@ -141,10 +142,10 @@ equilibrium and writes nothing. The raw forcing file contents are captured in `Input/RawInputs/ForcingTerms/`, and `preloaded_coil_sets` lets the rerun path inject the coil geometry stored in `Input/RawInputs/Coils/`. -Returns `(; equil, local_stability, ffs, match, resp, pe, kinetic, error_fields, tearing)`, each +Returns `(; equil, local_stability, ffs, match, growth_rate, resp, pe, kinetic, error_fields, tearing)`, each `nothing` when its section is absent. `equil` is the equilibrium the stages used (the solve's, on the two-pass grid), `ffs` the result the later stages read (the matched one when the deck has -`[Match]`), and `match` its `MatchResult`. +`[Match]`), `match` its `MatchResult`, and `growth_rate` the `IdealGrowthRateResult` of `[GrowthRate]`. """ function main_from_inputs( inputs::Dict{String,Any}, @@ -159,7 +160,7 @@ function main_from_inputs( runtimes = Vector{Pair{String,Float64}}() deck = parse_deck(inputs, path) (; output_file, verbose) = deck.run - local_stability = ffs = match = resp = pe = kinetic = error_fields = tearing = nothing + local_stability = ffs = match = growth = resp = pe = kinetic = error_fields = tearing = nothing equil = _timed_stage(runtimes, "equilibrium", "Equilibrium") do Equilibrium.PlasmaEquilibrium(eq_config, additional_input) @@ -189,6 +190,12 @@ function main_from_inputs( save_to_hdf5(target, ffs; git_version, local_stability) end @info "Results written to $output_file" + if deck.growth_rate !== nothing + growth = _timed_stage(runtimes, "growth_rate", "GrowthRate") do + growth_rate(ForceBalanceProblem(equil; deck.force_balance.keywords...), deck.force_balance.algorithm, deck.growth_rate; marginal=ffs, verbose) + end + save_to_hdf5(output_file, growth) + end end if deck.perturbed_equilibrium !== nothing @@ -230,7 +237,7 @@ function main_from_inputs( ffs === nothing || _write_runtimes!(output_file, runtimes) @info "\n$_BANNER\n GPEC completed successfully in $(@sprintf("%.3f", total_dt)) s\n$_BANNER" - return (; equil, local_stability, ffs, match, resp, pe, kinetic, error_fields, tearing) + return (; equil, local_stability, ffs, match, growth_rate=growth, resp, pe, kinetic, error_fields, tearing) end """ @@ -306,7 +313,7 @@ export main, write_imas export solve, plasma_response, perturbed_equilibrium, kinetic_torque, error_fields export KineticProfileData, IonSpecies, attach_kinetic_profiles! export PlasmaEquilibrium, ForceBalanceProblem, Forward, Riccati, Galerkin, TestKineticMatrices, KineticMatrices, kinetic_matrices, ForceFreeStatesResult, RMPField, PlasmaResponse -export PerpendicularInertia, IdealGrowthRateResult, growth_rate +export PerpendicularInertia, IdealGrowthRateResult, GrowthRateConfig, growth_rate, ion_diamagnetic_frequency export MatchProblem, TearingProblem, TearingConfig, TearingResult, GGJModel, SLAYERModel export LocalStabilityProblem, Ballooning, LocalStabilityResult export KineticForcesControl, KineticForcesResult, ErrorFieldsResult diff --git a/src/Utilities/HDF5Schema.jl b/src/Utilities/HDF5Schema.jl index f17b53511..f66e0fca2 100644 --- a/src/Utilities/HDF5Schema.jl +++ b/src/Utilities/HDF5Schema.jl @@ -289,6 +289,7 @@ 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/growth_rate" => (; long_name="wall-clock time of the ForceFreeStates growth-rate stage (inertial γ² root search)", units="s"), "Info/Runtimes/match" => (; long_name="wall-clock time of the [Match] inner-layer matching 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"), From 50a92c1fc4ff030eb0fec6fcd005321a962b439f Mon Sep 17 00:00:00 2001 From: d-burg Date: Sat, 7 Nov 2026 18:57:19 -0600 Subject: [PATCH 08/12] ForceFreeStates - TEST - Cover the growth-rate deck path, its HDF5 group and a harness case MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Deck-and-API parity for [GrowthRate] on the stable TJ fixture and on a q_a = 2.9 variant that is free-boundary unstable (γ agrees to 1e-12), the three parse-time guards, the writer's metadata, and a regression case pinning γ on the reversed-shear example with uniform deuterium density. Co-Authored-By: Claude Fable 5.1 --- .../lar_reversed_shear_growth_rate_n1.toml | 79 +++++++++++++++++++ test/runtests_inertia.jl | 16 ++++ test/runtests_solve_api.jl | 57 +++++++++++++ 3 files changed, 152 insertions(+) create mode 100644 regression-harness/cases/lar_reversed_shear_growth_rate_n1.toml diff --git a/regression-harness/cases/lar_reversed_shear_growth_rate_n1.toml b/regression-harness/cases/lar_reversed_shear_growth_rate_n1.toml new file mode 100644 index 000000000..d3351bbb5 --- /dev/null +++ b/regression-harness/cases/lar_reversed_shear_growth_rate_n1.toml @@ -0,0 +1,79 @@ +# Regression case: ELITE-style ideal growth rate of the large-aspect-ratio reversed-shear n=1 example. +# The example is free-boundary unstable at n=1, so [GrowthRate] root-finds a nonzero γ with a uniform +# deuterium mass density (1e20 m⁻³). Tracks γ, γ², the root-find history and the ideal energy it starts from. +# Each [quantities.*] block names an HDF5 path, how to extract it, and its noise floor. +[case] +name = "lar_reversed_shear_growth_rate_n1" +description = "Large-aspect-ratio reversed-shear equilibrium, n=1, ELITE-style ideal growth rate with uniform deuterium density" +example_dir = "examples/LAR_reversed_shear_example" + +[overrides] +"GrowthRate.rho" = 3.3e-7 # Uniform mass density [kg/m³] (deuterium at 1e20 m⁻³); NaN derives ρ(ψ) from the kinetic profiles + +# Growth rate (γ ≈ 3.0e3 1/s); absolute floors match a 1e-10 shift in λ_min near the root +[quantities.gamma] +h5path = "ForceFreeStates/GrowthRate/gamma" +type = "real_scalar" +extract = "value" +label = "growth rate gamma (1/s)" +noise_threshold = 1e-5 +order = 10 + +[quantities.gamma_sq] +h5path = "ForceFreeStates/GrowthRate/gamma_sq" +type = "real_scalar" +extract = "value" +label = "growth rate squared gamma^2 (1/s^2)" +noise_threshold = 1e-1 +order = 11 + +# Root-find history: every γ² evaluated, the least-stable free-boundary energy there, and the fixed-boundary zero count +[quantities.history_gamma_sq] +h5path = "ForceFreeStates/GrowthRate/History/gamma_sq" +type = "real_vector" +extract = "all_real" +label = "root-find gamma^2 history" +noise_threshold = 1e-1 +order = 20 + +[quantities.history_lambda_min] +h5path = "ForceFreeStates/GrowthRate/History/lambda_min" +type = "real_vector" +extract = "all_real" +label = "root-find lambda_min history" +noise_threshold = 1e-10 +order = 21 + +[quantities.history_nzero] +h5path = "ForceFreeStates/GrowthRate/History/nzero" +type = "real_vector" +extract = "all_real" +label = "root-find nzero history (integer counts)" +noise_threshold = 0 +order = 22 + +# Ideal marginal solve the root find starts from +[quantities.et_real] +h5path = "ForceFreeStates/FreeBoundaryStability/eigenmode_energies" +type = "complex_vector" +extract = "real_first" +label = "total energy Re(et[1])" +noise_threshold = 1e-10 +order = 30 + +[quantities.msing] +h5path = "SingularSurfaces/rational_count" +type = "int_scalar" +extract = "value" +label = "# singular surfaces" +noise_threshold = 0 +order = 31 + +# Runtime (special: not from H5) +[quantities.runtime] +h5path = "" +type = "runtime" +extract = "value" +label = "Runtime (s)" +noise_threshold = 0.0 +order = 999 diff --git a/test/runtests_inertia.jl b/test/runtests_inertia.jl index 33196f82a..3f3e99b7c 100644 --- a/test/runtests_inertia.jl +++ b/test/runtests_inertia.jl @@ -1,5 +1,7 @@ using TOML using LinearAlgebra +using HDF5 +include("h5_metadata_check.jl") # Perpendicular inertia (ELITE's kinetic energy K = ½∫ρ|ξ_⊥|² dV) as a kinetic source of the force-balance # equations: the A, C, H bands it adds, the B² identity they rely on, the γ²-scaled solves and growth_rate. @@ -134,7 +136,21 @@ using LinearAlgebra cfg = FFS.GrowthRateConfig() gr2 = FFS.growth_rate(kprob, forward, cfg; marginal=gr.marginal) @test gr2.gamma == gr.gamma && gr2.omega_star_pi == gr.omega_star_pi && gr2.marginal === gr.marginal + @test gr.config === nothing && gr.inertia === src && gr2.config === cfg && isnan(gr2.config.rho) @test FFS._rho_at(FFS.inertia_source(cfg, kequil), 0.4) == FFS._rho_at(src, 0.4) @test FFS._rho_at(FFS.inertia_source(FFS.GrowthRateConfig(; rho=2.0), kequil), 0.4) == 2.0 + # The HDF5 writer: the group, its annotations and the recorded [GrowthRate] inputs + mktempdir() do dir + path = joinpath(dir, "gpec.h5") + FFS.save_to_hdf5(path, gr2) + h5open(path, "r") do f + @test isempty(_collect_metadata_violations(f)) + @test read(f["ForceFreeStates/GrowthRate/gamma"]) == 0.0 + @test read(f["ForceFreeStates/GrowthRate/omega_star_pi"]) == gr2.omega_star_pi + @test length(read(f["ForceFreeStates/GrowthRate/History/lambda_min"])) == 1 + @test read(f["ForceFreeStates/GrowthRate/rho"]) ≈ [FFS._rho_at(src, x) for x in read(f["ForceFreeStates/GrowthRate/psi"])] + @test haskey(TOML.parse(read(f["Input/gpec_toml_raw"])), "GrowthRate") + end + end end end diff --git a/test/runtests_solve_api.jl b/test/runtests_solve_api.jl index 2c17a146b..724de553b 100644 --- a/test/runtests_solve_api.jl +++ b/test/runtests_solve_api.jl @@ -561,6 +561,63 @@ include("display_check.jl") end end + @testset "deck growth rate under [GrowthRate]" begin + config = GrowthRateConfig(; rho=1e-7) + growth_section = Dict{String,Any}("GrowthRate" => Dict{String,Any}("rho" => 1e-7)) + # The fixture is ideal-stable, so both paths stop at the marginal solve with γ = 0. + mktempdir() do dir + ffs = solve(equil, forward; nn=1, wall=wall, ffs_kwargs...) + api = growth_rate(ForceBalanceProblem(equil; nn=1, wall=wall, ffs_kwargs...), forward, config; marginal=ffs) + reference = GPEC.main([_stage_deck(dir; extra_sections=growth_section)]) + + @test reference.growth_rate isa IdealGrowthRateResult + @test reference.growth_rate.gamma == api.gamma == 0 + @test length(reference.growth_rate.history) == length(api.history) + @test [h.gamma_sq for h in reference.growth_rate.history] ≈ [h.gamma_sq for h in api.history] rtol = 1e-10 + @test [h.lambda for h in reference.growth_rate.history] ≈ [h.lambda for h in api.history] rtol = 1e-10 + @test isequal(reference.growth_rate.omega_star_pi, api.omega_star_pi) + @test reference.growth_rate.config == api.config + h5file = joinpath(dir, "gpec.h5") + @test h5open(f -> haskey(f, "Info/Runtimes/growth_rate"), h5file) + @test h5open(f -> haskey(f, "ForceFreeStates/GrowthRate/gamma"), h5file) + # The recorded input carries the [GrowthRate] section. + recorded = TOML.parse(h5open(f -> read(f, "Input/gpec_toml_raw"), h5file)) + @test haskey(recorded, "GrowthRate") + @test recorded["GrowthRate"]["rho"] == 1e-7 + end + + # Lowering q_a to 2.9 makes the fixture free-boundary unstable, so the root find runs and γ > 0 is compared. + kink_deck = merge(deck, Dict{String,Any}("TJAnalytic" => merge(deck["TJAnalytic"], Dict{String,Any}("qa" => 2.9)))) + mktempdir() do dir + kink_equil = GPEC.Equilibrium.PlasmaEquilibrium(GPEC.equilibrium_config(kink_deck, template), GPEC.analytic_config(kink_deck)) + ffs = solve(kink_equil, forward; nn=1, wall=wall, ffs_kwargs...) + @test minimum(real, ffs.free_boundary.et) < 0 + api = growth_rate(ForceBalanceProblem(kink_equil; nn=1, wall=wall, ffs_kwargs...), forward, config; marginal=ffs) + reference = GPEC.main([_stage_deck(dir; extra_sections=merge(growth_section, Dict{String,Any}("TJAnalytic" => kink_deck["TJAnalytic"])))]) + + @test api.gamma > 0 + @test reference.growth_rate.gamma ≈ api.gamma rtol = 1e-12 + @test length(reference.growth_rate.history) == length(api.history) > 1 + @test [h.gamma_sq for h in reference.growth_rate.history] ≈ [h.gamma_sq for h in api.history] rtol = 1e-10 + @test [h.lambda for h in reference.growth_rate.history] ≈ [h.lambda for h in api.history] rtol = 1e-10 + @test [h.nzero for h in reference.growth_rate.history] == [h.nzero for h in api.history] + end + + # The root find needs a Forward free-boundary solve as its plain ideal marginal one, so the deck stops at parse time. + parse_staged(dir; kwargs...) = (deck_dir = _stage_deck(dir; kwargs...); GPEC.parse_deck(TOML.parsefile(joinpath(deck_dir, "gpec.toml")), deck_dir)) + mktempdir() do dir + @test_throws "[GrowthRate] needs the free-boundary energy" parse_staged(dir; ffs_overrides=Dict{String,Any}("vac_flag" => false), extra_sections=growth_section) + end + mktempdir() do dir + @test_throws "[GrowthRate] needs integrator" parse_staged(dir; ffs_overrides=Dict{String,Any}("integrator" => "Riccati"), extra_sections=growth_section) + end + mktempdir() do dir + matched = merge(growth_section, Dict{String,Any}("Match" => Dict{String,Any}("eta" => 1e-7, "rho" => 1e-7, "rotation" => 1.0))) + @test_throws "[GrowthRate] needs the ideal solve as its marginal one" parse_staged(dir; ffs_overrides=Dict{String,Any}("integrator" => "Riccati"), + extra_sections=matched) + end + end + @testset "kinetic_format" begin # A rates-v1 kinetic file read through the deck's [Equilibrium] keys and through the API keyword. mktempdir() do dir From 99bb899e852bfacca19c7edc7bc5a119fd2f0894 Mon Sep 17 00:00:00 2001 From: d-burg Date: Sat, 7 Nov 2026 20:04:22 -0600 Subject: [PATCH 09/12] ForceFreeStates - BUGFIX - Bracket the growth rate on the mode count and keep the inertial solves on the marginal grid MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit From an adversarial review of the branch. The count of modes faster than γ (interior zero crossings plus a negative free-boundary energy) is what the decade bracket now follows, so a fixed-boundary- unstable marginal is no longer reported as γ = 0 and the bracket cannot step over the fastest root; a fixed-boundary fastest mode errors out. main hands growth_rate the problem its marginal solve came from and the inertial solves keep every control of that problem, so the two-pass grid is the same. Also: a converged flag with a warning when maxiter runs out, γ = 0 below gamma_sq_min, errors for a kinetic source and a non-positive density, ω_*pi at the least-stable mode's n and NaN without the profiles, no empty deck recorded from the API form, the inertia bands evaluated once. Co-Authored-By: Claude Fable 5.1 --- src/Deck.jl | 1 + src/ForceFreeStates/Inertia.jl | 157 ++++++++++++++++--------- src/ForceFreeStates/InertiaOutput.jl | 6 +- src/GeneralizedPerturbedEquilibrium.jl | 2 +- test/runtests_inertia.jl | 42 +++++++ test/runtests_solve_api.jl | 23 +++- 6 files changed, 171 insertions(+), 60 deletions(-) diff --git a/src/Deck.jl b/src/Deck.jl index c08cf05cb..1cd61eb50 100644 --- a/src/Deck.jl +++ b/src/Deck.jl @@ -477,6 +477,7 @@ function parse_deck(inputs::AbstractDict, path::AbstractString) if growth_rate !== nothing (; keywords, algorithm) = force_balance haskey(inputs, "Match") && error("[GrowthRate] needs the ideal solve as its marginal one: remove the [Match] section") + get(keywords, :kinetic, nothing) === nothing || error("[GrowthRate] needs an ideal force-balance solve: remove kinetic from [ForceFreeStates]") get(keywords, :vac_flag, false) || error("[GrowthRate] needs the free-boundary energy: set [ForceFreeStates] vac_flag = true") algorithm isa Forward || error("[GrowthRate] needs integrator = \"Forward\" in [ForceFreeStates]") end diff --git a/src/ForceFreeStates/Inertia.jl b/src/ForceFreeStates/Inertia.jl index 11f67cbbf..40ded43aa 100644 --- a/src/ForceFreeStates/Inertia.jl +++ b/src/ForceFreeStates/Inertia.jl @@ -33,6 +33,7 @@ kinetic profiles are attached (set `kinetic_file` or call `attach_kinetic_profil """ function PerpendicularInertia(equil::Equilibrium.PlasmaEquilibrium) data = InnerLayer._layer_kinetic_data(equil, "mass density ρ") + _require_profiles(data, "the inertia density ρ = μ m_p n_e") loaded = InnerLayer.SLAYER._slayer_profiles(data, nothing) mu = InnerLayer.SLAYER._layer_mass_per_electron(data, InnerLayer.SLAYER._layer_zeff(loaded, data)) n_e = loaded.profiles.n_e @@ -49,17 +50,26 @@ when γ > ω_*pi/2 at the pedestal. """ function ion_diamagnetic_frequency(equil::Equilibrium.PlasmaEquilibrium, n::Integer; z_i::Real=1) data = InnerLayer._layer_kinetic_data(equil, "ion diamagnetic frequency ω_*pi") - data.T_i === nothing && throw(ArgumentError("ω_*pi needs the ion temperature, and the kinetic data carries no T_i")) + _require_profiles(data, "ω_*pi") p = InnerLayer.SLAYER._slayer_profiles(data, nothing).profiles scale = -n / (z_i * equil.psio) return psi -> scale * (p.T_i(psi) * p.n_e(psi; deriv=DerivOp(1)) / p.n_e(psi) + p.T_i(psi; deriv=DerivOp(1))) end -# Peak |ω_*pi| over the integration domain of `r` at its lowest n; NaN without kinetic profiles that carry T_i. +# The shared kinetic-profile loader needs n_e, T_e and T_i even where only some of them enter. +_has_profiles(data) = data !== nothing && data.n_e !== nothing && data.T_e !== nothing && data.T_i !== nothing +function _require_profiles(data, what::AbstractString) + _has_profiles(data) || throw(ArgumentError("$what needs kinetic data with n_e, T_e and T_i (the shared profile loader reads all three)")) + return data +end + +# Toroidal mode number of the least-stable free-boundary mode of `r` (its lowest n without free-boundary energies). +_least_stable_n(r::ForceFreeStatesResult) = r.free_boundary === nothing ? r.nlow : r.nlow + r.free_boundary.n_tor_idx[1] + +# Peak |ω_*pi| over the integration domain of `r` at its least-stable mode's n; NaN without the kinetic profiles. function _peak_omega_star_pi(r::ForceFreeStatesResult) - data = r.equil.kinetic - (data === nothing || data.T_i === nothing) && return NaN - omega = ion_diamagnetic_frequency(r.equil, r.nlow) + _has_profiles(r.equil.kinetic) || return NaN + omega = ion_diamagnetic_frequency(r.equil, _least_stable_n(r)) psis = filter(x -> r.psilow <= x <= r.psilim, r.metric.xs) return isempty(psis) ? NaN : maximum(abs ∘ omega, psis) end @@ -81,7 +91,9 @@ function inertia_kinetic_matrices(src::PerpendicularInertia, equil::Equilibrium. integrands = zeros(Float64, mpsi, mtheta, 3) for ipsi in 1:mpsi q = equil.profiles.q_spline(metric.xs[ipsi]) - mu0rho = Utilities.MU_0 * _rho_at(src, metric.xs[ipsi]) + rho = _rho_at(src, metric.xs[ipsi]) + isfinite(rho) && rho > 0 || error("PerpendicularInertia: ρ(ψ_N = $(metric.xs[ipsi])) = $rho; the mass density must be finite and positive") + mu0rho = Utilities.MU_0 * rho for jtheta in 1:mtheta g11, g22, g33, g23, g31, g12, jac = @view metric.fs[ipsi, jtheta, 1:7] gradpsi2 = g22 * g33 - g23^2 @@ -159,10 +171,14 @@ The outcome of [`growth_rate`](@ref). - `gamma::Float64` - Growth rate in s⁻¹ of the least-stable free-boundary mode; 0 when the plasma is marginally stable or stable. - `gamma_sq::Float64` - γ² at which the least-stable free-boundary energy eigenvalue crosses zero. + - `converged::Bool` - Whether the root iteration met its `ftol` or `rtol` stop (always true for a stable or marginal plasma); false when + `maxiter` ran out, in which case `gamma` is the last iterate and `history` brackets the root. - `inertia::PerpendicularInertia` - The inertia source the solves used. - `config::Union{Nothing,GrowthRateConfig}` - The deck settings when the solve came from a `[GrowthRate]` section or the three-argument method; `nothing` otherwise. - - `omega_star_pi::Float64` - Peak ion diamagnetic frequency |ω_*pi| in rad/s over the integration domain at the problem's lowest n, from the - equilibrium's kinetic profiles (NaN when none carry T_i); the EPED1 criterion counts the mode as unstable when `gamma > omega_star_pi/2`. + - `omega_star_pi::Float64` - Peak ion diamagnetic frequency |ω_*pi| in rad/s over the integration domain, at the toroidal n of the + least-stable free-boundary mode of the solve at `gamma_sq` (of the marginal solve when stable), from the equilibrium's kinetic + profiles (NaN without n_e, T_e and T_i); the EPED1 criterion counts the mode as unstable when `gamma > omega_star_pi/2`. The + peak is taken over `[psilow, psilim]`, so a truncated domain can cut off part of a pedestal gradient. - `marginal::ForceFreeStatesResult` - The ideal (γ² = 0) solve, whose least-stable energy decides stability. - `result::Union{Nothing,ForceFreeStatesResult}` - The solve at `gamma_sq`, carrying the eigenmode; `nothing` when stable. - `history::Vector` - Every (γ², λ_min, nzero) evaluated, in order, with λ_min the least-stable free-boundary energy at that γ². @@ -170,6 +186,7 @@ The outcome of [`growth_rate`](@ref). struct IdealGrowthRateResult gamma::Float64 gamma_sq::Float64 + converged::Bool inertia::PerpendicularInertia config::Union{Nothing,GrowthRateConfig} omega_star_pi::Float64 @@ -183,79 +200,103 @@ end rtol=1e-3, ftol=1e-6, maxiter=30, marginal=nothing, config=nothing, verbose=false) -> IdealGrowthRateResult growth_rate(prob::ForceBalanceProblem, alg::Forward, config::GrowthRateConfig; marginal=nothing, verbose=false) -Growth rate of the least-stable free-boundary mode of `prob` with the perpendicular inertia `inertia` -(a [`PerpendicularInertia`](@ref)), ELITE's eigenvalue problem solved by shooting: the ideal solve -gives the marginal least-stable energy λ_min(0); when it is negative, the inertia enters as -`kinetic_factor = γ²` and λ_min(γ²), which rises monotonically with γ² once the interior is -positive, is brought to zero. The root is bracketed by decades from `gamma_sq_init` and closed by a -bracketed secant (Illinois) iteration that stops when |λ_min| falls below `ftol` times the marginal -|λ_min(0)| or the bracket is narrower than `rtol` in γ², within `maxiter` steps; every evaluation is -one force-balance solve. `gamma_sq_min` is the floor below which the mode counts as marginal, -`gamma_sq_max` the ceiling beyond which bracketing stops with an error. - -A `marginal` result of the same problem (an earlier ideal `solve`) is reused instead of solving it again; `config` -is recorded on the result and is what the three-argument method passes. -The problem's `singfac_min` is set to 0 for the inertial solves unless given: with inertia the -rational surfaces are regular, so the integration runs as a single chunk. A fixed-boundary mode -faster than the free-boundary one shows as `nzero > 0` at the root in `history`; this first version -reports the free-boundary root and leaves the interior count to the caller. +Growth rate of the fastest ideal mode of `prob` with the perpendicular inertia `inertia` (a +[`PerpendicularInertia`](@ref)), ELITE's eigenvalue problem solved by shooting. The inertia enters as +`kinetic_factor = γ²`; at each γ² the number of modes growing faster than γ is the interior zero-crossing +count `nzero` plus one when the least-stable free-boundary energy λ_min is negative, and that count cannot +rise with γ² since the kinetic energy is positive. The ideal solve gives the count at γ = 0; when it is +zero the plasma is stable and `gamma` is 0. Otherwise the count is bracketed by decades from +`gamma_sq_init` (down while nothing is faster, up while something still is), and when the lower end has +no interior crossing left, only λ_min changes sign in the bracket and a bracketed secant (Illinois) +iteration brings it to zero, stopping when |λ_min| falls below `ftol` times the reference |λ_min| (the +marginal value, or the bracket's negative end when the marginal energy is positive) or the bracket is +narrower than `rtol` in γ², within `maxiter` steps (a warning and `converged = false` otherwise). Every +evaluation is one force-balance solve. `gamma_sq_min` is the floor below which the mode counts as marginal +(`gamma` = 0), `gamma_sq_max` the ceiling beyond which bracketing stops with an error. A fastest mode that +is fixed-boundary (an interior crossing still present at the bracket's lower end) is reported as an error: +locating a fixed-boundary growth rate is not supported yet. + +A `marginal` result of the same problem (an earlier ideal `solve`) is reused instead of solving it again; +`config` is recorded on the result and is what the three-argument method passes. The problem must carry no +kinetic source. The inertial solves take every other control of the problem, so a two-pass ψ grid is the +same as the marginal one; with `singfac_min > 0` the kinetic surface finder can still place crossings at +the sharp cond(F̂) peaks that appear near rational surfaces when γ²ρ is far below any physical root, where +only the sign of λ_min matters, and `singfac_min = 0` gives a single chunk (the rational surfaces are +regular with inertia). """ function growth_rate(prob::ForceBalanceProblem, alg::Forward; inertia::PerpendicularInertia, gamma_sq_init::Real=1e9, gamma_sq_min::Real=1e4, gamma_sq_max::Real=1e16, rtol::Real=1e-3, ftol::Real=1e-6, maxiter::Int=30, marginal::Union{Nothing,ForceFreeStatesResult}=nothing, config::Union{Nothing,GrowthRateConfig}=nothing, verbose::Bool=false) + get(prob.ctrl_kwargs, :kinetic, nothing) === nothing || + error("growth_rate needs an ideal marginal problem: remove its kinetic source, the inertia is the only kinetic term") history = @NamedTuple{gamma_sq::Float64, lambda::Float64, nzero::Int}[] lambda_min(r) = minimum(real, r.free_boundary.et) + faster(h) = h.nzero + (h.lambda < 0) # modes growing faster than γ; non-increasing in γ² + solves = Dict{Float64,ForceFreeStatesResult}() + function finish(gamma_sq, converged) + r = get(solves, gamma_sq, nothing) + return IdealGrowthRateResult(sqrt(gamma_sq), gamma_sq, converged, inertia, config, _peak_omega_star_pi(something(r, marginal)), marginal, r, history) + end marginal = marginal === nothing ? solve(prob, alg; verbose) : marginal marginal.free_boundary === nothing && error("growth_rate needs the free-boundary energies: set vac_flag=true") - omega_star_pi = _peak_omega_star_pi(marginal) push!(history, (gamma_sq=0.0, lambda=lambda_min(marginal), nzero=marginal.diagnostics.nzero)) verbose && @info @sprintf("growth_rate: marginal λ_min = %.6e (nzero = %d)", history[end].lambda, history[end].nzero) - history[end].lambda < 0 || return IdealGrowthRateResult(0.0, 0.0, inertia, config, omega_star_pi, marginal, nothing, history) + faster(history[1]) == 0 && return finish(0.0, true) kwargs = copy(prob.ctrl_kwargs) - haskey(kwargs, :singfac_min) || (kwargs[:singfac_min] = 0.0) kwargs[:kinetic] = inertia - matrices = Ref{Any}(nothing) # the inertia bands, evaluated once and reused at every γ² - last = Ref{Any}(nothing) - seen = Dict{Float64,Float64}() # Brent re-evaluates its bracket ends; each γ² is one solve - function lambda_at(gamma_sq) - haskey(seen, gamma_sq) && return seen[gamma_sq] + kwargs[:kinetic_factor] = Float64(gamma_sq_init) + kwargs[:kinetic] = kinetic_matrices(ForceBalanceProblem(prob.equil, prob.wall, kwargs), alg; verbose) # the bands once, reused at every γ² + seen = Dict{Float64,Int}() # γ² → history index; each γ² is one solve + function evaluate(gamma_sq) + gamma_sq = Float64(gamma_sq) + haskey(seen, gamma_sq) && return history[seen[gamma_sq]] kw = copy(kwargs) - matrices[] === nothing || (kw[:kinetic] = matrices[]) kw[:kinetic_factor] = gamma_sq - inertial = ForceBalanceProblem(prob.equil, prob.wall, kw) - matrices[] === nothing && (matrices[] = kinetic_matrices(inertial, alg; verbose)) - r = solve(inertial, alg; verbose) - last[] = r - push!(history, (gamma_sq=Float64(gamma_sq), lambda=lambda_min(r), nzero=r.diagnostics.nzero)) + r = solve(ForceBalanceProblem(prob.equil, prob.wall, kw), alg; verbose) + solves[gamma_sq] = r + push!(history, (gamma_sq=gamma_sq, lambda=lambda_min(r), nzero=r.diagnostics.nzero)) + seen[gamma_sq] = length(history) verbose && @info @sprintf("growth_rate: γ² = %.4e λ_min = %.6e nzero = %d", gamma_sq, history[end].lambda, history[end].nzero) - return seen[Float64(gamma_sq)] = history[end].lambda + return history[end] end - # Bracket by decades: down while already positive, up while still negative + # Bracket the count by decades: down while nothing is faster, up while something still is hi = Float64(gamma_sq_init) - lam = lambda_at(hi) - while lam > 0 + h_hi = evaluate(hi) + while faster(h_hi) == 0 hi /= 10 - hi < gamma_sq_min && return IdealGrowthRateResult(sqrt(hi), hi, inertia, config, omega_star_pi, marginal, last[], history) - lam = lambda_at(hi) + if hi < gamma_sq_min + verbose && @info @sprintf("growth_rate: nothing grows faster than γ² = %.2e, below gamma_sq_min: marginal", 10hi) + return finish(0.0, true) + end + h_hi = evaluate(hi) end - lo = hi - while lam <= 0 - lo = hi + lo, h_lo = hi, h_hi + while faster(h_hi) > 0 + lo, h_lo = hi, h_hi hi *= 10 - hi > gamma_sq_max && error("growth_rate: λ_min is still negative at γ² = $lo; raise gamma_sq_max or check the inertia density") - lam = lambda_at(hi) + hi > gamma_sq_max && error("growth_rate: a mode still grows at γ² = $lo; raise gamma_sq_max or check the inertia density") + h_hi = evaluate(hi) end + # The fastest root lies in (lo, hi]; with no interior crossing left at lo only the free-boundary energy changes sign there + h_lo.nzero == 0 || error("growth_rate: the fastest mode is fixed-boundary (nzero = $(h_lo.nzero) at γ² = $lo, 0 at γ² = $hi); " * + "locating a fixed-boundary growth rate is not supported yet") - # Illinois iteration on the sign-change bracket; the retained end's λ is halved when it is kept twice running - flo, fhi = seen[lo], seen[hi] - gamma_sq, side = hi, 0 + # Illinois iteration on the sign change; the retained end's λ is halved when it is kept twice running + scale = history[1].lambda < 0 ? abs(history[1].lambda) : abs(h_lo.lambda) + flo, fhi = h_lo.lambda, h_hi.lambda + gamma_sq, side, converged = hi, 0, false for _ in 1:maxiter gamma_sq = (lo * fhi - hi * flo) / (fhi - flo) - lam = lambda_at(gamma_sq) - abs(lam) <= ftol * abs(history[1].lambda) && break + h = evaluate(gamma_sq) + h.nzero == 0 || error("growth_rate: an interior zero crossing appeared inside the bracket at γ² = $gamma_sq; the free-boundary energy is not monotone there") + lam = h.lambda + if abs(lam) <= ftol * scale + converged = true + break + end if lam < 0 lo, flo = gamma_sq, lam side == -1 && (fhi /= 2) @@ -265,9 +306,13 @@ function growth_rate(prob::ForceBalanceProblem, alg::Forward; inertia::Perpendic side == 1 && (flo /= 2) side = 1 end - hi - lo <= rtol * hi && break + if hi - lo <= rtol * hi + converged = true + break + end end - return IdealGrowthRateResult(sqrt(gamma_sq), gamma_sq, inertia, config, omega_star_pi, marginal, last[], history) + converged || @warn "growth_rate: no convergence in $maxiter steps; γ² lies in [$lo, $hi] and the last iterate $gamma_sq is reported" + return finish(gamma_sq, converged) end function growth_rate(prob::ForceBalanceProblem, alg::Forward, config::GrowthRateConfig; marginal::Union{Nothing,ForceFreeStatesResult}=nothing, verbose::Bool=false) diff --git a/src/ForceFreeStates/InertiaOutput.jl b/src/ForceFreeStates/InertiaOutput.jl index 215082671..7126f5ef6 100644 --- a/src/ForceFreeStates/InertiaOutput.jl +++ b/src/ForceFreeStates/InertiaOutput.jl @@ -1,13 +1,14 @@ # Writer for an `IdealGrowthRateResult`. Layout under the caller's parent group (inputs live only under Input/): # # ForceFreeStates/GrowthRate/ -# ├── gamma, gamma_sq, omega_star_pi -- growth rate, its square and the EPED1 ω_*pi reference +# ├── gamma, gamma_sq, converged, omega_star_pi -- growth rate, its square, the stop flag and the EPED1 ω_*pi reference # ├── psi, rho -- inertia density on the marginal solve's ψ_N grid # └── History/ -- every (γ², λ_min, nzero) evaluated by the root search const GROWTH_RATE_H5_ANNOTATIONS = [ "gamma" => (; long_name="growth rate γ of the least-stable free-boundary mode with perpendicular inertia (0 when stable)", units="s^-1"), "gamma_sq" => (; long_name="γ² at which the least-stable free-boundary energy eigenvalue crosses zero (0 when stable)", units="s^-2"), + "converged" => (; long_name="1 when the root iteration met its stop criterion, 0 when maxiter ran out (gamma is then the last iterate)"), "omega_star_pi" => (; long_name="peak ion diamagnetic frequency |ω_*pi| over the integration domain at the lowest n (NaN without kinetic T_i)", units="rad/s"), "psi" => (; long_name="normalized poloidal flux ψ_N grid of the marginal solve", dims=("psi",), scale="psi"), @@ -29,6 +30,7 @@ function save_to_hdf5(parent::Union{HDF5.File,HDF5.Group}, result::IdealGrowthRa g = Utilities.HDF5Annotations.replace_group!(ffs, "GrowthRate") g["gamma"] = result.gamma g["gamma_sq"] = result.gamma_sq + g["converged"] = Int(result.converged) g["omega_star_pi"] = result.omega_star_pi psi = collect(Float64, result.marginal.metric.xs) g["psi"] = psi @@ -37,6 +39,6 @@ function save_to_hdf5(parent::Union{HDF5.File,HDF5.Group}, result::IdealGrowthRa g["History/lambda_min"] = [h.lambda for h in result.history] g["History/nzero"] = [h.nzero for h in result.history] Utilities.HDF5Annotations.annotate!(g, GROWTH_RATE_H5_ANNOTATIONS) - Utilities.record_inputs!(parent, result) + result.config === nothing || Utilities.record_inputs!(parent, result) # the API form has no deck section to record return g end diff --git a/src/GeneralizedPerturbedEquilibrium.jl b/src/GeneralizedPerturbedEquilibrium.jl index b4caff512..cbb3acd4b 100755 --- a/src/GeneralizedPerturbedEquilibrium.jl +++ b/src/GeneralizedPerturbedEquilibrium.jl @@ -192,7 +192,7 @@ function main_from_inputs( @info "Results written to $output_file" if deck.growth_rate !== nothing growth = _timed_stage(runtimes, "growth_rate", "GrowthRate") do - growth_rate(ForceBalanceProblem(equil; deck.force_balance.keywords...), deck.force_balance.algorithm, deck.growth_rate; marginal=ffs, verbose) + growth_rate(ffs.problem, deck.force_balance.algorithm, deck.growth_rate; marginal=ffs, verbose) # the problem the marginal solve came from end save_to_hdf5(output_file, growth) end diff --git a/test/runtests_inertia.jl b/test/runtests_inertia.jl index 3f3e99b7c..3e329e0fd 100644 --- a/test/runtests_inertia.jl +++ b/test/runtests_inertia.jl @@ -146,6 +146,7 @@ include("h5_metadata_check.jl") h5open(path, "r") do f @test isempty(_collect_metadata_violations(f)) @test read(f["ForceFreeStates/GrowthRate/gamma"]) == 0.0 + @test read(f["ForceFreeStates/GrowthRate/converged"]) == 1 @test read(f["ForceFreeStates/GrowthRate/omega_star_pi"]) == gr2.omega_star_pi @test length(read(f["ForceFreeStates/GrowthRate/History/lambda_min"])) == 1 @test read(f["ForceFreeStates/GrowthRate/rho"]) ≈ [FFS._rho_at(src, x) for x in read(f["ForceFreeStates/GrowthRate/psi"])] @@ -153,4 +154,45 @@ include("h5_metadata_check.jl") end end end + + @testset "root search on a free-boundary-unstable fixture" begin + # q_a = 2.9 makes the TJ fixture external-kink unstable at n = 1 + kink_deck = merge(deck, Dict{String,Any}("TJAnalytic" => merge(deck["TJAnalytic"], Dict{String,Any}("qa" => 2.9)))) + kink_equil = GPEC.Equilibrium.PlasmaEquilibrium(GPEC.equilibrium_config(kink_deck, template), GPEC.analytic_config(kink_deck)) + kink_prob = ForceBalanceProblem(kink_equil; nn=1, wall, ffs_kwargs...) + inertia = FFS.PerpendicularInertia(; rho=1e-7) + gr = FFS.growth_rate(kink_prob, forward; inertia) + @test gr.gamma > 0 && gr.converged && gr.result !== nothing + @test all(h.nzero == 0 for h in gr.history) + @test abs(minimum(real, gr.result.free_boundary.et)) <= 1e-6 * abs(gr.history[1].lambda) + @test gr.result.closure == :regular && isnan(gr.omega_star_pi) + # Running out of iterations is reported, not hidden + short = @test_logs (:warn, r"no convergence in 1 steps") match_mode = :any FFS.growth_rate(kink_prob, forward; inertia, maxiter=1, marginal=gr.marginal) + @test !short.converged && short.gamma > 0 && short.result !== nothing + # A root below gamma_sq_min counts as marginal: γ = 0, no inertial solve reported + floor = FFS.growth_rate(kink_prob, forward; inertia, gamma_sq_min=1e12, marginal=gr.marginal) + @test floor.gamma == 0 && floor.gamma_sq == 0 && floor.result === nothing && floor.converged + @test all(h.lambda > 0 for h in floor.history[2:end]) + # The marginal problem must be ideal, and the density positive + @test_throws "ideal marginal problem" FFS.growth_rate(ForceBalanceProblem(kink_equil; nn=1, wall, ffs_kwargs..., kinetic=GPEC.KineticForcesControl()), forward; inertia) + @test_throws "finite and positive" FFS.inertia_kinetic_matrices(FFS.PerpendicularInertia(; rho=0.0), equil, intr, metric) + end + + @testset "kinetic data without T_e" begin + ktemplate = joinpath(@__DIR__, "test_data", "regression_tj_kinetic_calculated") + kdeck = TOML.parsefile(joinpath(ktemplate, "gpec.toml")) + kequil = GPEC.Equilibrium.PlasmaEquilibrium(GPEC.equilibrium_config(kdeck, ktemplate), GPEC.analytic_config(kdeck)) + data = kequil.kinetic + GPEC.attach_kinetic_profiles!(kequil, GPEC.Equilibrium.KineticProfileData(; psi=data.psi, n_e=data.n_e, T_i=data.T_i)) + @test_throws ArgumentError FFS.PerpendicularInertia(kequil) + @test_throws ArgumentError FFS.ion_diamagnetic_frequency(kequil, 1) + kkw = Dict{Symbol,Any}(q.first => q.second for q in GPEC.force_balance_keywords(kdeck, ktemplate) if !(q.first in (:nn, :wall, :kinetic, :kinetic_factor))) + kkw[:vac_flag] = true + gr = FFS.growth_rate( + ForceBalanceProblem(kequil; nn=1, wall=GPEC.Vacuum.parse_wall(get(kdeck, "Wall", nothing)), kkw...), + forward; + inertia=FFS.PerpendicularInertia(; rho=1e-7) + ) + @test gr.gamma == 0 && isnan(gr.omega_star_pi) + end end diff --git a/test/runtests_solve_api.jl b/test/runtests_solve_api.jl index 724de553b..085ba3bad 100644 --- a/test/runtests_solve_api.jl +++ b/test/runtests_solve_api.jl @@ -603,8 +603,29 @@ include("display_check.jl") @test [h.nzero for h in reference.growth_rate.history] == [h.nzero for h in api.history] end - # The root find needs a Forward free-boundary solve as its plain ideal marginal one, so the deck stops at parse time. parse_staged(dir; kwargs...) = (deck_dir = _stage_deck(dir; kwargs...); GPEC.parse_deck(TOML.parsefile(joinpath(deck_dir, "gpec.toml")), deck_dir)) + # On the two-pass auto grid the deck must hand growth_rate the problem its marginal solve came from, so both paths re-form from pass 1. + auto_overrides = Dict{String,Any}("grid_type" => "auto", "mpsi" => 0, "psi_accuracy" => 3e-2) + mktempdir() do dir + auto_deck = merge(kink_deck, Dict{String,Any}("Equilibrium" => merge(kink_deck["Equilibrium"], auto_overrides))) + pass1 = GPEC.Equilibrium.PlasmaEquilibrium(GPEC.equilibrium_config(auto_deck, template), GPEC.analytic_config(auto_deck)) + api = growth_rate(ForceBalanceProblem(pass1; nn=1, wall=wall, ffs_kwargs...), forward, config) + reference = + GPEC.main([_stage_deck(dir; eq_overrides=auto_overrides, extra_sections=merge(growth_section, Dict{String,Any}("TJAnalytic" => kink_deck["TJAnalytic"])))]) + @test api.gamma > 0 + @test reference.growth_rate.gamma ≈ api.gamma rtol = 1e-10 + @test [h.lambda for h in reference.growth_rate.history] ≈ [h.lambda for h in api.history] rtol = 1e-10 + end + + # A kinetic source in [ForceFreeStates] would make the marginal solve kinetic-MHD while the inertial ones are ideal. + mktempdir() do dir + deck_dir = _stage_deck(dir; eq_overrides=Dict{String,Any}("kinetic_file" => "kinetic.dat"), ffs_overrides=Dict{String,Any}("kinetic" => "KineticForces"), + extra_sections=growth_section) + cp(joinpath(@__DIR__, "test_data", "regression_tj_kinetic_calculated", "kinetic.dat"), joinpath(deck_dir, "kinetic.dat")) + @test_throws "[GrowthRate] needs an ideal force-balance solve" GPEC.parse_deck(TOML.parsefile(joinpath(deck_dir, "gpec.toml")), deck_dir) + end + + # The root find needs a Forward free-boundary solve as its plain ideal marginal one, so the deck stops at parse time. mktempdir() do dir @test_throws "[GrowthRate] needs the free-boundary energy" parse_staged(dir; ffs_overrides=Dict{String,Any}("vac_flag" => false), extra_sections=growth_section) end From 8481a64dd33c9c934e52025cf15443af1aed32bc Mon Sep 17 00:00:00 2001 From: d-burg Date: Sun, 8 Nov 2026 00:45:51 -0600 Subject: [PATCH 10/12] ForceFreeStates - BUGFIX - Bisect the growth-rate bracket on the crossing count and pin the missed-count cases MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit A crossing still present at the bracket's lower end is a fixed-boundary eigenvalue inside the bracket, where the free-boundary energy has a pole. The fastest mode is free-boundary and at least as fast as any fixed-boundary one, so growth_rate now bisects on the count until the lower end clears that eigenvalue instead of erroring; a crossing still present at rtol means the fastest mode has no edge displacement and the upper end is reported. The count comes from the conjugate-point detector, which misses a crossing when a stiff harmonic holds the smallest-magnitude eigenvalue of the inverse energy matrix near zero; on the TJ fixture at p0 = 0.8 that hides a fixed-boundary mode at γ² ≈ 1.2e8, so a search started below it lands on a slower root or reports stability behind a close wall. Both outcomes are pinned as test_broken and the docstring says to start gamma_sq_init above the fastest expected growth rate. Co-Authored-By: Claude Fable 5.1 --- src/ForceFreeStates/Inertia.jl | 45 +++++++++++++++++++++++++--------- test/runtests_inertia.jl | 19 ++++++++++++++ 2 files changed, 52 insertions(+), 12 deletions(-) diff --git a/src/ForceFreeStates/Inertia.jl b/src/ForceFreeStates/Inertia.jl index 40ded43aa..df80d4397 100644 --- a/src/ForceFreeStates/Inertia.jl +++ b/src/ForceFreeStates/Inertia.jl @@ -206,15 +206,23 @@ Growth rate of the fastest ideal mode of `prob` with the perpendicular inertia ` count `nzero` plus one when the least-stable free-boundary energy λ_min is negative, and that count cannot rise with γ² since the kinetic energy is positive. The ideal solve gives the count at γ = 0; when it is zero the plasma is stable and `gamma` is 0. Otherwise the count is bracketed by decades from -`gamma_sq_init` (down while nothing is faster, up while something still is), and when the lower end has -no interior crossing left, only λ_min changes sign in the bracket and a bracketed secant (Illinois) -iteration brings it to zero, stopping when |λ_min| falls below `ftol` times the reference |λ_min| (the -marginal value, or the bracket's negative end when the marginal energy is positive) or the bracket is -narrower than `rtol` in γ², within `maxiter` steps (a warning and `converged = false` otherwise). Every -evaluation is one force-balance solve. `gamma_sq_min` is the floor below which the mode counts as marginal -(`gamma` = 0), `gamma_sq_max` the ceiling beyond which bracketing stops with an error. A fastest mode that -is fixed-boundary (an interior crossing still present at the bracket's lower end) is reported as an error: -locating a fixed-boundary growth rate is not supported yet. +`gamma_sq_init` (down while nothing is faster, up while something still is). An interior crossing still +present at the lower end is a fixed-boundary eigenvalue inside the bracket, where λ_min has a pole; since +the fastest mode is free-boundary and at least as fast as any fixed-boundary one, the lower end is lifted +above it by bisecting on the count. Then only λ_min changes sign in the bracket and a bracketed secant +(Illinois) iteration brings it to zero, stopping when |λ_min| falls below `ftol` times the reference +|λ_min| (the marginal value, or the bracket's negative end when the marginal energy is positive) or the +bracket is narrower than `rtol` in γ², within `maxiter` steps (a warning and `converged = false` +otherwise). Every evaluation is one force-balance solve. `gamma_sq_min` is the floor below which the mode +counts as marginal (`gamma` = 0), `gamma_sq_max` the ceiling beyond which bracketing stops with an error. +When the count bisection reaches `rtol` with a crossing still present, the fastest mode is a fixed-boundary +one (its edge displacement vanishes) and the bracket's upper end is reported as converged. + +The count is only as good as the crossing detector, which follows the smallest-magnitude eigenvalue of W⁻¹ +along ψ and misses a conjugate point whenever a stiff harmonic holds that eigenvalue near zero. A missed +crossing below `gamma_sq_init` makes a positive λ_min look stable, so the search can settle on a slower +mode or report γ = 0. Start `gamma_sq_init` above the fastest growth rate expected: the upward scan needs +only the sign of λ_min, and `history` shows whether λ_min rose monotonically through the bracket. A `marginal` result of the same problem (an earlier ideal `solve`) is reused instead of solving it again; `config` is recorded on the result and is what the three-argument method passes. The problem must carry no @@ -280,9 +288,22 @@ function growth_rate(prob::ForceBalanceProblem, alg::Forward; inertia::Perpendic hi > gamma_sq_max && error("growth_rate: a mode still grows at γ² = $lo; raise gamma_sq_max or check the inertia density") h_hi = evaluate(hi) end - # The fastest root lies in (lo, hi]; with no interior crossing left at lo only the free-boundary energy changes sign there - h_lo.nzero == 0 || error("growth_rate: the fastest mode is fixed-boundary (nzero = $(h_lo.nzero) at γ² = $lo, 0 at γ² = $hi); " * - "locating a fixed-boundary growth rate is not supported yet") + # The fastest root lies in (lo, hi]. A crossing still present at lo is a fixed-boundary eigenvalue inside the bracket, where + # λ_min has a pole; the fastest mode is free-boundary and at least as fast (interlacing), so bisect on the count, which is + # monotone, until lo is above that eigenvalue and only the free-boundary energy changes sign in the bracket + while h_lo.nzero > 0 && hi - lo > rtol * hi + mid = sqrt(lo * hi) + h_mid = evaluate(mid) + if faster(h_mid) == 0 + hi, h_hi = mid, h_mid + else + lo, h_lo = mid, h_mid + end + end + if h_lo.nzero > 0 # the root sits on the pole to within rtol: the fastest mode has no edge displacement + verbose && @info @sprintf("growth_rate: the fastest mode is fixed-boundary to within rtol, γ² in [%.6e, %.6e]", lo, hi) + return finish(hi, true) + end # Illinois iteration on the sign change; the retained end's λ is halved when it is kept twice running scale = history[1].lambda < 0 ? abs(history[1].lambda) : abs(h_lo.lambda) diff --git a/test/runtests_inertia.jl b/test/runtests_inertia.jl index 3e329e0fd..0d3fc88a6 100644 --- a/test/runtests_inertia.jl +++ b/test/runtests_inertia.jl @@ -178,6 +178,25 @@ include("h5_metadata_check.jl") @test_throws "finite and positive" FFS.inertia_kinetic_matrices(FFS.PerpendicularInertia(; rho=0.0), equil, intr, metric) end + @testset "a fixed-boundary mode inside the bracket" begin + # p0 = 0.8 adds a fixed-boundary mode at γ² ≈ 1.2e8 (ρ = 1e-7) below the external kink at ≈ 4.9e9, a pole of λ_min + # between them. The crossing counter watches the smallest-magnitude eigenvalue of W⁻¹, which a stiff harmonic holds + # near zero on this fixture, so it misses the conjugate point; until the count is exact these pin the wrong answers. + hp_deck = merge(deck, Dict{String,Any}("TJAnalytic" => merge(deck["TJAnalytic"], Dict{String,Any}("qa" => 2.9, "p0" => 0.8)))) + hp_equil = GPEC.Equilibrium.PlasmaEquilibrium(GPEC.equilibrium_config(hp_deck, template), GPEC.analytic_config(hp_deck)) + inertia = FFS.PerpendicularInertia(; rho=1e-7) + # Starting below the pole, only the count can carry the search past it to the fastest mode + low_start = FFS.growth_rate(ForceBalanceProblem(hp_equil; nn=1, wall, ffs_kwargs...), forward; inertia, gamma_sq_init=1e8) + @test low_start.converged + @test_broken low_start.gamma_sq > 4e9 + # A close wall makes the marginal free-boundary energy positive, so only the count says the plasma is unstable, and the + # bracket [1e8, 1e9] then holds the pole: the count bisection must lift its lower end before the sign iteration + close_wall = GPEC.Vacuum.ConformalWall(; standoff_fraction=0.05) + walled = FFS.growth_rate(ForceBalanceProblem(hp_equil; nn=1, wall=close_wall, ffs_kwargs...), forward; inertia, gamma_sq_init=1e8) + @test_broken walled.gamma > 0 + @test_broken any(h.nzero > 0 for h in walled.history) + end + @testset "kinetic data without T_e" begin ktemplate = joinpath(@__DIR__, "test_data", "regression_tj_kinetic_calculated") kdeck = TOML.parsefile(joinpath(ktemplate, "gpec.toml")) From c2d6aff8892992cffbc3daff1c38f7fbedc47306 Mon Sep 17 00:00:00 2001 From: d-burg Date: Sun, 8 Nov 2026 09:43:14 -0600 Subject: [PATCH 11/12] ForceFreeStates - BUGFIX - Count conjugate points exactly from the Cayley eigenphases of the solution frame MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The DCON detector behind nzero follows the smallest-magnitude eigenvalue of W⁻¹ between stored steps and misses a conjugate point whenever a stiff harmonic pins a tiny eigenvalue near zero, so growth_rate could bracket a slower mode or report stability. The forward integration now also counts conjugate points exactly, at every accepted step, as the eigenphases of the Cayley transform Ψ = (U₁ − i s U₂)(U₁ + i s U₂)⁻¹ that cross π (its Maslov index): Ψ is invariant under the Gaussian reductions, has no poles, and its eigenvectors identify each phase across a step. Phases are seeded one step into each chunk so the vanishing columns at the axis and after an ideal crossing are not counted, a step with a fast-moving phase is bisected on the integrator's dense output, and a unitarity defect above 1e-3 (a frame that is not Lagrangian to working precision) flags the count. growth_rate switches the count on for its solves through an internal keyword, uses it in place of nzero, and marks a flagged result not converged; the deck path turns it on for the marginal solve when a [GrowthRate] section follows. The three pinned cases on the TJ fixture now pass. Co-Authored-By: Claude Fable 5.1 --- src/ForceFreeStates/CoreTypes.jl | 3 + src/ForceFreeStates/EulerLagrange.jl | 31 +++++ src/ForceFreeStates/FixedBoundaryStability.jl | 123 ++++++++++++++++++ src/ForceFreeStates/ForceBalance.jl | 11 +- src/ForceFreeStates/Inertia.jl | 44 +++++-- src/ForceFreeStates/InertiaOutput.jl | 2 +- src/GeneralizedPerturbedEquilibrium.jl | 4 +- test/runtests_inertia.jl | 6 +- 8 files changed, 201 insertions(+), 23 deletions(-) diff --git a/src/ForceFreeStates/CoreTypes.jl b/src/ForceFreeStates/CoreTypes.jl index 1cb794f65..89eedba73 100644 --- a/src/ForceFreeStates/CoreTypes.jl +++ b/src/ForceFreeStates/CoreTypes.jl @@ -37,6 +37,8 @@ Solve-time scratch state for a force-balance solve: the resolved mode space, sin - `q1lim::Float64` - Safety factor derivative at psilim - `wall_settings::Union{Nothing,Vacuum.WallSettings}` - Wall shape settings for vacuum calculations - `verbose::Bool` - Log solve progress; set from the `verbose` keyword of `solve` + - `count_conjugate_points::Bool` - Run the exact conjugate-point count of the forward integration (`OdeState.nzero_conj`); set by + `growth_rate` for its solves, off otherwise """ @kwdef mutable struct ForceFreeStatesInternal <: ModeSpace mlow::Int = 0 @@ -59,6 +61,7 @@ Solve-time scratch state for a force-balance solve: the resolved mode space, sin q1lim::Float64 = 0.0 wall_settings::Union{Nothing,Vacuum.WallSettings} = nothing verbose::Bool = false + count_conjugate_points::Bool = false """ Inter-surface Δ' matrix, `msing × msing`, PEST3 convention: the odd-parity projection of `delta_prime_raw`, set by `compute_delta_prime_matrix!` (Riccati path only). diff --git a/src/ForceFreeStates/EulerLagrange.jl b/src/ForceFreeStates/EulerLagrange.jl index d40700059..52ad96a6e 100644 --- a/src/ForceFreeStates/EulerLagrange.jl +++ b/src/ForceFreeStates/EulerLagrange.jl @@ -67,6 +67,19 @@ records needed to undo it, and singular-surface asymptotic coefficients. - `ising_start::Int` - Index of the last singular surface inside the starting ψ. - `needs_crossing::Bool` - Whether a rational surface must be crossed after the current region. - `nzero::Int` - Number of detected zero crossings of the critical eigenvalue. + - `count_conj::Bool` - Whether the exact conjugate-point count runs (forward Hamiltonian runs that request it: no kinetic source, or a + `PerpendicularInertia` one); off by default since it costs two LU factorizations and one eigenvalue solve per accepted step. + - `nconj::Int` - Signed running count of Cayley eigenphases that wrapped through π (conjugate points, fixed-boundary modes inside the + ψ interval), reseeded after each ideal crossing so it counts per Newcomb interval. + - `nzero_conj::Int` - Number of conjugate points, `abs(nconj)` at the end of the run; −1 when the count did not run. + - `nconj_flags::Int` - Accepted steps whose count is unreliable (an eigenphase rotation above π/2, or a unitarity defect above 1e-3 that + shows the frame is not Lagrangian to working precision); a nonzero value means `nzero_conj` may be wrong. + - `maslov_theta::Vector{Float64}` - Cayley eigenphases at the last counted step, aligned with the columns of `maslov_vecs`. + - `maslov_vecs::Matrix{ComplexF64}` - Cayley eigenvectors at the last counted step, which identify each eigenphase across a step. + - `maslov_pending::Bool` - Whether the next accepted step seeds the eigenphases and the Cayley scale instead of advancing the count. + - `maslov_s::Float64` - Scale s of the Cayley transform (U₁ − i s U₂)(U₁ + i s U₂)⁻¹, set to ‖U₁‖/‖U₂‖ of the column-normalized frame one step + into each chunk. + - `maslov_defect::Float64` - Largest unitarity defect max| |λ| − 1 | of the Cayley eigenvalues over the run. - `new::Bool` - Whether the next norm call records fresh reference norms `unorm0`. - `unorm::Vector{Float64}` - Current solution-column norms. - `unorm0::Vector{Float64}` - Reference norms since the last reduction. @@ -111,6 +124,17 @@ records needed to undo it, and singular-surface asymptotic coefficients. needs_crossing::Bool = false nzero::Int = 0 + # Exact conjugate-point count + count_conj::Bool = false + nconj::Int = 0 + nzero_conj::Int = -1 + nconj_flags::Int = 0 + maslov_theta::Vector{Float64} = zeros(Float64, numpert_total) + maslov_vecs::Matrix{ComplexF64} = zeros(ComplexF64, numpert_total, numpert_total) + maslov_pending::Bool = false + maslov_s::Float64 = 1.0 + maslov_defect::Float64 = 0.0 + # Used for Gaussian reduction new::Bool = true unorm::Vector{Float64} = zeros(Float64, numpert_total) @@ -285,6 +309,7 @@ function forward_eulerlagrange_integration(ctrl::ForceFreeStatesControl, equil:: odet = OdeState(intr.numpert_total, ctrl.algorithm.numsteps_init, ctrl.algorithm.numunorms_init, intr.msing) initialize_el_at_axis!(odet, ctrl, mats, equil.profiles, intr) + odet.count_conj = intr.count_conjugate_points && kinetic_source(ctrl.kinetic) isa Union{Nothing,PerpendicularInertia} chunks = chunk_el_integration_bounds(odet, ctrl, intr) @@ -293,6 +318,8 @@ function forward_eulerlagrange_integration(ctrl::ForceFreeStatesControl, equil:: end for chunk in chunks + # Reseeds the count one step into the chunk: at the Frobenius start and after an ideal crossing U₁ has vanishing columns + odet.count_conj && maslov_seed!(odet) integrate_el_region!(odet, ctrl, equil, mats, intr, chunk) if verbose @info " ψ = $((@sprintf "%.3f" odet.psifac)), q = $((@sprintf "%.3f" odet.q)), steps = $(odet.total_steps)" @@ -318,6 +345,7 @@ function forward_eulerlagrange_integration(ctrl::ForceFreeStatesControl, equil:: @info "Evaluating fixed-boundary stability criterion" end odet.nzero = evaluate_stability_criterion!(odet, equil.profiles) + odet.count_conj && (odet.nzero_conj = abs(odet.nconj)) # Undo Gaussian reduction to get true solution vectors transform_u!(odet, intr) @@ -663,6 +691,7 @@ function cross_kinetic_singular_surf!( odet.psifac = ksurf.psifac + dpsi sing_der!(du2, odet.u, params, odet.psifac) odet.u .+= (du1 .+ du2) .* dpsi + odet.count_conj && maslov_advance!(odet, odet.u) # compared with the phases of the chunk's last accepted step odet.q = equil.profiles.q_spline(odet.psifac; hint=odet.spline_hint) store_ode_data!(odet, odet.psifac, odet.u) @@ -746,6 +775,8 @@ function integrate_el_region!( compute_solution_norms!(integrator.u, odet, ctrl, intr, false) column_abstol!(integrator.opts.abstol, integrator.u, integrator.opts.reltol) + # The dense output lets a step with a fast-moving eigenphase be bisected for the conjugate-point count + odet.count_conj && maslov_advance!(odet, integrator.u; refine=τ -> integrator(integrator.tprev + τ * (integrator.t - integrator.tprev))) # The first step is always saved, even when q_range ≈ 0 near_start = abs(odet.q - q_start) < near_q_frac * q_range || steps_in_segment[] == 1 diff --git a/src/ForceFreeStates/FixedBoundaryStability.jl b/src/ForceFreeStates/FixedBoundaryStability.jl index 29896ac3d..b34113570 100644 --- a/src/ForceFreeStates/FixedBoundaryStability.jl +++ b/src/ForceFreeStates/FixedBoundaryStability.jl @@ -106,3 +106,126 @@ smallest magnitude. Also flags whether W⁻¹ was non-Hermitian beyond tolerance crit = evals[findmin(abs, evals)[2]] return crit, nonherm end + +""" + cayley_eigen!(θ, V, u1, u2, s) -> defect + +Eigen-decompose the Cayley transform Ψ = (U₁ − i s U₂)(U₁ + i s U₂)⁻¹ of the solution frame `[U₁; U₂]`, +writing the eigenphases in (−π, π] to `θ` and the eigenvectors to the columns of `V`, and return the +unitarity defect max| |λ| − 1 |. For a Lagrangian frame Ψ is unitary; an eigenvalue −1 occurs exactly where +det U₁ = 0, a conjugate point, while an eigenvalue +1 marks a singular U₂, a harmless pole of W⁻¹ (Arnold +1967, Funct. Anal. Appl. 1, 1; Robbin & Salamon 1993, Topology 32, 827). Ψ is invariant under U → U·T for +any invertible T, so Gaussian reduction leaves it unchanged and the frame can be passed with normalized +columns (see [`normalized_frame`](@ref)). When U₁†U₂ drifts from Hermitian, Ψ is not exactly unitary; the +eigenpairs are used regardless and the defect reports the drift. +""" +function cayley_eigen!(θ::AbstractVector{Float64}, V::AbstractMatrix{ComplexF64}, u1::AbstractMatrix{ComplexF64}, u2::AbstractMatrix{ComplexF64}, s::Real) + psi = (u1 .- (im * s) .* u2) / (u1 .+ (im * s) .* u2) + E = eigen!(psi) + defect = 0.0 + for (j, l) in enumerate(E.values) + defect = max(defect, abs(abs(l) - 1)) + a = angle(l) + θ[j] = a <= -π ? Float64(π) : a + end + V .= E.vectors + return defect +end + +""" + normalized_frame(u) -> (u1, u2) + +Copies of the two halves of the solution frame `u = [U₁; U₂]` with each column of the stacked frame scaled +to unit norm, so that the Cayley inverse is conditioned by the frame's directions and not by the column +norms, which span many decades between Gaussian reductions. +""" +function normalized_frame(u::AbstractArray{ComplexF64,3}) + u1 = copy(@view u[:, :, 1]) + u2 = copy(@view u[:, :, 2]) + for j in axes(u1, 2) + c = hypot(norm(@view u1[:, j]), norm(@view u2[:, j])) + c > 0 && (u1[:, j] ./= c; u2[:, j] ./= c) + end + return u1, u2 +end + +""" + maslov_match(θold, Vold, θnew, Vnew) -> (w, maxrot, overlap) + +Follow each eigenphase from one step to the next through the eigenvectors, assigning old to new by greatest +overlap |⟨v_old, v_new⟩| first until every pair is matched (the eigenvectors of a unitary Ψ are orthonormal, +so a clear step gives overlaps near 1; within a cluster of nearly equal phases the assignment is arbitrary +but harmless, since swapped partners rotate by the same small arcs). Each phase moves along the shortest arc +δ = angle(exp(i(θnew − θold))); return the signed number of arcs that cross π (a conjugate point, +1 upward +and −1 downward), the largest |δ|, and the smallest overlap used. +""" +function maslov_match(θold::AbstractVector{Float64}, Vold::AbstractMatrix{ComplexF64}, θnew::AbstractVector{Float64}, Vnew::AbstractMatrix{ComplexF64}) + O = abs.(Vold' * Vnew) + N = length(θnew) + partner = zeros(Int, N) + overlap = 1.0 + for _ in 1:N + o, idx = findmax(O) + j, k = Tuple(idx) + partner[k] = j + overlap = min(overlap, o) + O[j, :] .= -1.0 + O[:, k] .= -1.0 + end + w, maxrot = 0, 0.0 + for k in 1:N + d = angle(cis(θnew[k] - θold[partner[k]])) + maxrot = max(maxrot, abs(d)) + target = θold[partner[k]] + d + w += target > π ? 1 : target <= -π ? -1 : 0 + end + return w, maxrot, overlap +end + +""" + maslov_seed!(odet) + +Mark the conjugate-point count for reseeding at the next accepted step. At a chunk start U₁ has vanishing +columns (every column at the axis, the zeroed resonant column after an ideal crossing), so eigenphases sit +exactly at π there; seeding one step later counts conjugate points in the open interval, as Newcomb's +theorem states it, and sets the Cayley scale s = ‖U₁‖/‖U₂‖ from a nonsingular frame. +""" +function maslov_seed!(odet::OdeState) + odet.maslov_pending = true + return odet +end + +""" + maslov_advance!(odet, u; refine=nothing, depth=0) + +Advance the conjugate-point count from the stored eigenpairs to those of frame `u`: add the signed number +of eigenphases that crossed π to `odet.nconj`. A phase that rotated by more than π/4 in the step leaves its +direction of travel in doubt (an eigenvalue of W⁻¹ passing zero or infinity quickly), so when `refine(τ)` +can supply the frame at fraction τ of the step (the integrator's dense output), the step is bisected, up to +six times, until every arc is short. The step is counted in `odet.nconj_flags` when a rotation above π/2 +remains or the unitarity defect exceeds 1e-3 (the frame is not Lagrangian to working precision, so its +eigenphases are unreliable). A pending seed stores the eigenpairs and the scale instead. +""" +function maslov_advance!(odet::OdeState, u::AbstractArray{ComplexF64,3}; refine=nothing, depth::Int=0) + u1, u2 = normalized_frame(u) + if odet.maslov_pending + n1, n2 = norm(u1), norm(u2) + odet.maslov_s = n1 > 0 && n2 > 0 ? n1 / n2 : 1.0 + odet.maslov_defect = max(odet.maslov_defect, cayley_eigen!(odet.maslov_theta, odet.maslov_vecs, u1, u2, odet.maslov_s)) + odet.maslov_pending = false + return odet + end + θnew, Vnew = similar(odet.maslov_theta), similar(odet.maslov_vecs) + defect = cayley_eigen!(θnew, Vnew, u1, u2, odet.maslov_s) + w, maxrot, _ = maslov_match(odet.maslov_theta, odet.maslov_vecs, θnew, Vnew) + if maxrot > π / 4 && refine !== nothing && depth < 6 + maslov_advance!(odet, refine(0.5); refine=τ -> refine(0.5 * τ), depth=depth + 1) + return maslov_advance!(odet, u; refine=τ -> refine(0.5 + 0.5 * τ), depth=depth + 1) + end + odet.maslov_defect = max(odet.maslov_defect, defect) + (maxrot > π / 2 || defect > 1e-3) && (odet.nconj_flags += 1) + odet.nconj += w + odet.maslov_theta .= θnew + odet.maslov_vecs .= Vnew + return odet +end diff --git a/src/ForceFreeStates/ForceBalance.jl b/src/ForceFreeStates/ForceBalance.jl index cf4aba8e5..57542f7cc 100644 --- a/src/ForceFreeStates/ForceBalance.jl +++ b/src/ForceFreeStates/ForceBalance.jl @@ -92,9 +92,12 @@ ffs = solve(prob, Riccati()) ffs = solve(eq, Riccati(); nn=1, vac_flag=true) # equivalent one-line form ``` """ -function CommonSolve.solve(prob::ForceBalanceProblem, alg::AbstractIntegrator; verbose::Bool=false) +CommonSolve.solve(prob::ForceBalanceProblem, alg::AbstractIntegrator; verbose::Bool=false) = _solve_force_balance(prob, alg; verbose) + +# `solve` with the internal switch for the exact conjugate-point count, which growth_rate turns on for its solves. +function _solve_force_balance(prob::ForceBalanceProblem, alg::AbstractIntegrator; verbose::Bool=false, count_conjugate_points::Bool=false) total_start = time() - ctrl, intr, equil, metric, mats, layer_overlap = _prepare_solve(prob, alg; verbose) + ctrl, intr, equil, metric, mats, layer_overlap = _prepare_solve(prob, alg; verbose, count_conjugate_points) result = run_force_free_states(ctrl, equil, mats, intr, metric; cut_solution=alg isa Galerkin && alg.cut_solution) result = _with(result; problem=prob, algorithm=alg, layer_overlap=layer_overlap) @@ -123,7 +126,7 @@ end # Everything solve does before the integration: the control, the mode space, the (re-formed) equilibrium, the metric, # the matrices, the layer-overlap scan and the kinetic matrices. -function _prepare_solve(prob::ForceBalanceProblem, alg::AbstractIntegrator; verbose::Bool=false) +function _prepare_solve(prob::ForceBalanceProblem, alg::AbstractIntegrator; verbose::Bool=false, count_conjugate_points::Bool=false) equil = prob.equil ctrl_kwargs = copy(prob.ctrl_kwargs) ctrl = ForceFreeStatesControl(; algorithm=alg, ctrl_kwargs...) @@ -148,7 +151,7 @@ function _prepare_solve(prob::ForceBalanceProblem, alg::AbstractIntegrator; verb "nodes; regenerate the file in double precision for converged Δ′." file = eq_config.eq_filename end - intr = ForceFreeStatesInternal(; verbose) + intr = ForceFreeStatesInternal(; verbose, count_conjugate_points) intr.wall_settings = prob.wall resolve_mode_space!(intr, ctrl) diff --git a/src/ForceFreeStates/Inertia.jl b/src/ForceFreeStates/Inertia.jl index df80d4397..91a669c13 100644 --- a/src/ForceFreeStates/Inertia.jl +++ b/src/ForceFreeStates/Inertia.jl @@ -172,7 +172,8 @@ The outcome of [`growth_rate`](@ref). - `gamma::Float64` - Growth rate in s⁻¹ of the least-stable free-boundary mode; 0 when the plasma is marginally stable or stable. - `gamma_sq::Float64` - γ² at which the least-stable free-boundary energy eigenvalue crosses zero. - `converged::Bool` - Whether the root iteration met its `ftol` or `rtol` stop (always true for a stable or marginal plasma); false when - `maxiter` ran out, in which case `gamma` is the last iterate and `history` brackets the root. + `maxiter` ran out, in which case `gamma` is the last iterate and `history` brackets the root, or when any solve flagged an ambiguous + step in its conjugate-point count. - `inertia::PerpendicularInertia` - The inertia source the solves used. - `config::Union{Nothing,GrowthRateConfig}` - The deck settings when the solve came from a `[GrowthRate]` section or the three-argument method; `nothing` otherwise. - `omega_star_pi::Float64` - Peak ion diamagnetic frequency |ω_*pi| in rad/s over the integration domain, at the toroidal n of the @@ -181,7 +182,8 @@ The outcome of [`growth_rate`](@ref). peak is taken over `[psilow, psilim]`, so a truncated domain can cut off part of a pedestal gradient. - `marginal::ForceFreeStatesResult` - The ideal (γ² = 0) solve, whose least-stable energy decides stability. - `result::Union{Nothing,ForceFreeStatesResult}` - The solve at `gamma_sq`, carrying the eigenmode; `nothing` when stable. - - `history::Vector` - Every (γ², λ_min, nzero) evaluated, in order, with λ_min the least-stable free-boundary energy at that γ². + - `history::Vector` - Every (γ², λ_min, nzero) evaluated, in order, with λ_min the least-stable free-boundary energy at that γ² and + `nzero` the exact conjugate-point count of that solve. """ struct IdealGrowthRateResult gamma::Float64 @@ -202,8 +204,8 @@ end Growth rate of the fastest ideal mode of `prob` with the perpendicular inertia `inertia` (a [`PerpendicularInertia`](@ref)), ELITE's eigenvalue problem solved by shooting. The inertia enters as -`kinetic_factor = γ²`; at each γ² the number of modes growing faster than γ is the interior zero-crossing -count `nzero` plus one when the least-stable free-boundary energy λ_min is negative, and that count cannot +`kinetic_factor = γ²`; at each γ² the number of modes growing faster than γ is the number of conjugate +points `nzero` plus one when the least-stable free-boundary energy λ_min is negative, and that count cannot rise with γ² since the kinetic energy is positive. The ideal solve gives the count at γ = 0; when it is zero the plasma is stable and `gamma` is 0. Otherwise the count is bracketed by decades from `gamma_sq_init` (down while nothing is faster, up while something still is). An interior crossing still @@ -218,11 +220,10 @@ counts as marginal (`gamma` = 0), `gamma_sq_max` the ceiling beyond which bracke When the count bisection reaches `rtol` with a crossing still present, the fastest mode is a fixed-boundary one (its edge displacement vanishes) and the bracket's upper end is reported as converged. -The count is only as good as the crossing detector, which follows the smallest-magnitude eigenvalue of W⁻¹ -along ψ and misses a conjugate point whenever a stiff harmonic holds that eigenvalue near zero. A missed -crossing below `gamma_sq_init` makes a positive λ_min look stable, so the search can settle on a slower -mode or report γ = 0. Start `gamma_sq_init` above the fastest growth rate expected: the upward scan needs -only the sign of λ_min, and `history` shows whether λ_min rose monotonically through the bracket. +The conjugate points are counted exactly, as the eigenphases of the Cayley transform of the solution frame +that wrap through π along ψ (its Maslov index), not by the smallest-eigenvalue detector `diagnostics.nzero`; +a step whose phase match is ambiguous is warned about and marks the result not converged. Starting +`gamma_sq_init` above the fastest growth rate expected is cheapest, since the upward scan needs only the sign of λ_min. A `marginal` result of the same problem (an earlier ideal `solve`) is reused instead of solving it again; `config` is recorded on the result and is what the three-argument method passes. The problem must carry no @@ -239,16 +240,30 @@ function growth_rate(prob::ForceBalanceProblem, alg::Forward; inertia::Perpendic error("growth_rate needs an ideal marginal problem: remove its kinetic source, the inertia is the only kinetic term") history = @NamedTuple{gamma_sq::Float64, lambda::Float64, nzero::Int}[] lambda_min(r) = minimum(real, r.free_boundary.et) + counted_solve(p) = _solve_force_balance(p, alg; verbose, count_conjugate_points=true) + counted(r) = r.diagnostics !== nothing && r.diagnostics.nzero_conj >= 0 + flagged = Ref(false) + function check_flags(r, gamma_sq) + nflag = r.diagnostics.nconj_flags + nflag > 0 || return + flagged[] = true + @warn "growth_rate: the conjugate-point count at γ² = $gamma_sq had $nflag ambiguous step(s); the result is marked not converged" + end faster(h) = h.nzero + (h.lambda < 0) # modes growing faster than γ; non-increasing in γ² solves = Dict{Float64,ForceFreeStatesResult}() function finish(gamma_sq, converged) r = get(solves, gamma_sq, nothing) - return IdealGrowthRateResult(sqrt(gamma_sq), gamma_sq, converged, inertia, config, _peak_omega_star_pi(something(r, marginal)), marginal, r, history) + return IdealGrowthRateResult(sqrt(gamma_sq), gamma_sq, converged && !flagged[], inertia, config, _peak_omega_star_pi(something(r, marginal)), marginal, r, history) end - marginal = marginal === nothing ? solve(prob, alg; verbose) : marginal + if marginal !== nothing && !counted(marginal) + verbose && @info "growth_rate: the supplied marginal solve has no conjugate-point count; solving it again with the count" + marginal = nothing + end + marginal = marginal === nothing ? counted_solve(prob) : marginal marginal.free_boundary === nothing && error("growth_rate needs the free-boundary energies: set vac_flag=true") - push!(history, (gamma_sq=0.0, lambda=lambda_min(marginal), nzero=marginal.diagnostics.nzero)) + check_flags(marginal, 0.0) + push!(history, (gamma_sq=0.0, lambda=lambda_min(marginal), nzero=marginal.diagnostics.nzero_conj)) verbose && @info @sprintf("growth_rate: marginal λ_min = %.6e (nzero = %d)", history[end].lambda, history[end].nzero) faster(history[1]) == 0 && return finish(0.0, true) @@ -262,9 +277,10 @@ function growth_rate(prob::ForceBalanceProblem, alg::Forward; inertia::Perpendic haskey(seen, gamma_sq) && return history[seen[gamma_sq]] kw = copy(kwargs) kw[:kinetic_factor] = gamma_sq - r = solve(ForceBalanceProblem(prob.equil, prob.wall, kw), alg; verbose) + r = counted_solve(ForceBalanceProblem(prob.equil, prob.wall, kw)) solves[gamma_sq] = r - push!(history, (gamma_sq=gamma_sq, lambda=lambda_min(r), nzero=r.diagnostics.nzero)) + check_flags(r, gamma_sq) + push!(history, (gamma_sq=gamma_sq, lambda=lambda_min(r), nzero=r.diagnostics.nzero_conj)) seen[gamma_sq] = length(history) verbose && @info @sprintf("growth_rate: γ² = %.4e λ_min = %.6e nzero = %d", gamma_sq, history[end].lambda, history[end].nzero) return history[end] diff --git a/src/ForceFreeStates/InertiaOutput.jl b/src/ForceFreeStates/InertiaOutput.jl index 7126f5ef6..9612e60c8 100644 --- a/src/ForceFreeStates/InertiaOutput.jl +++ b/src/ForceFreeStates/InertiaOutput.jl @@ -15,7 +15,7 @@ const GROWTH_RATE_H5_ANNOTATIONS = [ "rho" => (; long_name="perpendicular-inertia mass density ρ", units="kg/m^3", dims=("psi",), attach=(1 => "psi",)), "History/gamma_sq" => (; long_name="γ² evaluated by the growth-rate root search, in order", units="s^-2", dims=("iteration",)), "History/lambda_min" => (; long_name="power-normalized least-stable free-boundary energy at each evaluated γ² (per unit ⟨|ξ|²⟩)", dims=("iteration",)), - "History/nzero" => (; long_name="number of zero crossings of the fixed-boundary critical eigenvalue at each evaluated γ²", dims=("iteration",)) + "History/nzero" => (; long_name="number of conjugate points (fixed-boundary modes inside the ψ interval) at each evaluated γ²", dims=("iteration",)) ] """ diff --git a/src/GeneralizedPerturbedEquilibrium.jl b/src/GeneralizedPerturbedEquilibrium.jl index cbb3acd4b..504423aac 100755 --- a/src/GeneralizedPerturbedEquilibrium.jl +++ b/src/GeneralizedPerturbedEquilibrium.jl @@ -171,7 +171,9 @@ function main_from_inputs( if deck.force_balance !== nothing ffs = _timed_stage(runtimes, "force_free_states", "Force-Free States") do - solve(ForceBalanceProblem(equil; deck.force_balance.keywords...), deck.force_balance.algorithm; verbose) + # growth_rate reuses this solve as its marginal, so it carries the conjugate-point count when a [GrowthRate] section follows + ForceFreeStates._solve_force_balance(ForceBalanceProblem(equil; deck.force_balance.keywords...), deck.force_balance.algorithm; verbose, + count_conjugate_points=deck.growth_rate !== nothing) end # Every later stage, local stability included, sees the equilibrium the solve formed (re-formed on the two-pass grid). equil = ffs.equil diff --git a/test/runtests_inertia.jl b/test/runtests_inertia.jl index 0d3fc88a6..368dac2f3 100644 --- a/test/runtests_inertia.jl +++ b/test/runtests_inertia.jl @@ -188,13 +188,13 @@ include("h5_metadata_check.jl") # Starting below the pole, only the count can carry the search past it to the fastest mode low_start = FFS.growth_rate(ForceBalanceProblem(hp_equil; nn=1, wall, ffs_kwargs...), forward; inertia, gamma_sq_init=1e8) @test low_start.converged - @test_broken low_start.gamma_sq > 4e9 + @test low_start.gamma_sq > 4e9 # A close wall makes the marginal free-boundary energy positive, so only the count says the plasma is unstable, and the # bracket [1e8, 1e9] then holds the pole: the count bisection must lift its lower end before the sign iteration close_wall = GPEC.Vacuum.ConformalWall(; standoff_fraction=0.05) walled = FFS.growth_rate(ForceBalanceProblem(hp_equil; nn=1, wall=close_wall, ffs_kwargs...), forward; inertia, gamma_sq_init=1e8) - @test_broken walled.gamma > 0 - @test_broken any(h.nzero > 0 for h in walled.history) + @test walled.gamma > 0 + @test any(h.nzero > 0 for h in walled.history) end @testset "kinetic data without T_e" begin From 98fa97ba49948334f53d41c25613b1578e92fe4f Mon Sep 17 00:00:00 2001 From: d-burg Date: Sun, 8 Nov 2026 10:19:46 -0600 Subject: [PATCH 12/12] ForceFreeStates - DOCS - Document the exact conjugate-point count and its references Explain in the stability page how nzero (DCON's smallest-eigenvalue test) and nzero_conj (the Maslov index of the solution frame through its Cayley eigenphases) differ, what keeps the latter exact (seeding one step into each chunk, eigenvector tracking, bisection of fast steps on the dense output, the unitarity-defect flag), its cost and its validation; add Arnold 1967 and Robbin & Salamon 1993 for the Maslov index and Wilson 2002, Snyder 2002 and Lunniss 2016 for ELITE's formulation and solve procedure to the references list. Co-Authored-By: Claude Fable 5.1 --- docs/development/references.md | 19 +++++++++++++++++++ docs/src/stability.md | 9 +++++++++ src/ForceFreeStates/FixedBoundaryStability.jl | 2 +- 3 files changed, 29 insertions(+), 1 deletion(-) diff --git a/docs/development/references.md b/docs/development/references.md index 9dbb3c411..4aa06c67a 100644 --- a/docs/development/references.md +++ b/docs/development/references.md @@ -34,6 +34,25 @@ The ForceFreeStates module (ideal MHD stability analysis) implements methods fro - Published: Physics of Plasmas **25**, 032507 (2018) - Describes: Riccati method for ideal MHD eigenvalue problem +- **Arnold (1967)**: "On a characteristic class entering into conditions of quantization" + - Published: Functional Analysis and Its Applications **1**, 1 (1967) + - Describes: the Maslov index of a path of Lagrangian subspaces, counted through the eigenvalues of a unitary (Cayley) image; basis of the exact conjugate-point count `nzero_conj` in `FixedBoundaryStability.jl` + +- **Robbin & Salamon (1993)**: "The Maslov index for paths" + - Published: Topology **32**, 827 (1993) + - Describes: the Maslov index as a signed count of crossings of a Lagrangian path with a reference subspace, which is what the Cayley eigenphases crossing π compute step by step + +- **Wilson, Snyder, Huysmans & Miller (2002)**: "Numerical studies of edge localized instabilities in tokamaks" + - Published: Physics of Plasmas **9**, 1277 (2002) + - Describes: the ELITE formulation; Eq. 9 is the perpendicular (incompressible) inertia that `PerpendicularInertia` implements; the growth rate is found by shooting with a secant iteration on an edge matching eigenvalue + +- **Snyder, Wilson et al. (2002)**: "Edge localized modes and the pedestal: A model based on coupled peeling–ballooning modes" + - Published: Physics of Plasmas **9**, 2037 (2002) + - Describes: ELITE's peeling–ballooning application and the γ > ω_*pi/2 diamagnetic-stabilization test that `IdealGrowthRateResult.omega_star_pi` serves + +- **Lunniss (2016)**: "Arbitrary-n ELITE" (PhD thesis, University of York; etheses.whiterose.ac.uk/16638) + - Describes: ELITE's radial finite-difference sweep (App. A.1), the edge matching matrix T(γ²) and the fictitious eigenvalue λ(γ²) driven to zero (§2.5.8), and the arbitrary-n inertia options; the comparison point for the count-based bracket of `growth_rate`, which ELITE does not have + ## PerturbedEquilibrium Module The PerturbedEquilibrium module implements GPEC-style perturbed equilibrium calculations from: diff --git a/docs/src/stability.md b/docs/src/stability.md index 98d61d697..74875de15 100644 --- a/docs/src/stability.md +++ b/docs/src/stability.md @@ -141,6 +141,15 @@ The perpendicular diffusivities that set the layer widths are described with the - Both integrators use the adaptive ninth-order `Vern9` solver at `eulerlagrange_tolerance`; the Riccati basis is recovered by renormalizing whenever the solution exceeds `renorm_threshold`. - The two crossing conventions differ slightly, so Riccati and Forward free-boundary energies agree to within a few percent, not to machine precision. +### Counting fixed-boundary modes + +Two counts of the conjugate points (fixed-boundary modes fitting inside the ψ interval) come out of the forward integration. + +- `nzero` is DCON's criterion: a sign change of the smallest-magnitude eigenvalue of W⁻¹ = U₁U₂⁻¹ between stored steps, confirmed at the interpolated zero. It is cheap and is what the `crit` scan in the output records, but it misses a crossing whenever a stiff harmonic holds a tiny eigenvalue of W⁻¹ near zero at every stored step, and it can report crossings where the frame has lost its Hermitian structure. +- `nzero_conj` is the Maslov index of the solution frame (Arnold 1967; Robbin & Salamon 1993), computed at every accepted step and exact for a Lagrangian frame. The Cayley transform Ψ = (U₁ − i s U₂)(U₁ + i s U₂)⁻¹ of the column-normalized frame is unitary when U₁†U₂ is Hermitian, has an eigenvalue −1 exactly where det U₁ = 0 (a conjugate point) and +1 where U₂ is singular (a pole of W⁻¹, harmless), and is invariant under U → U·T, so the Gaussian reductions do not disturb it. Each eigenphase is followed through its eigenvector from one step to the next and counts one conjugate point when its shortest arc crosses π. The phases are seeded one accepted step into each chunk, because at the axis start and after an ideal crossing U₁ has vanishing columns whose phases sit exactly at π and belong to the interval's end, not its interior; a step in which a phase rotates by more than π/4 is bisected on the integrator's dense output until the arcs are short, since an eigenvalue of W⁻¹ can pass zero within one step near a conjugate point. The scale s is set from the frame one step into the chunk and never changed within it: a change of s moves every phase but none through π, yet a large s squashes the phases toward π until a pole of W⁻¹ becomes indistinguishable from a zero crossing, so rescaling is not a substitute for bisection. + +`nconj_flags` counts the steps where the count cannot be trusted: a rotation above π/2 after bisection, or a unitarity defect max||λ|−1| of Ψ above 1e-3, which means the integrated frame is no longer Lagrangian to working precision (the Solovev regression fixture shows this after its q = 2 crossing, where `nzero` also misfires). `maslov_defect` records the largest defect seen. The count costs two matrix inversions and one eigendecomposition per accepted step, roughly 10–15 % of a solve, so it is off by default; `growth_rate` switches it on for all of its solves and uses it, with the sign of the least free-boundary energy, as the number of modes faster than the trial γ, which is what makes its bracket monotone. Checked on the TJ analytic equilibrium at `qa = 2.9, p0 = 0.8`, where a fixed-boundary mode at γ² ≈ 1.2e8 is invisible to `nzero`: `nzero_conj` is 1 below that eigenvalue and 0 above, and the Morse index `nzero_conj` + (negative free-boundary eigenvalues) is non-increasing in γ² and continuous across the pole, as it must be. + Local stability (Mercier, resistive interchange and ballooning) is a separate stage on the same equilibrium; see [Local Stability](ballooning.md). ## API reference diff --git a/src/ForceFreeStates/FixedBoundaryStability.jl b/src/ForceFreeStates/FixedBoundaryStability.jl index b34113570..d557fc04e 100644 --- a/src/ForceFreeStates/FixedBoundaryStability.jl +++ b/src/ForceFreeStates/FixedBoundaryStability.jl @@ -116,7 +116,7 @@ unitarity defect max| |λ| − 1 |. For a Lagrangian frame Ψ is unitary; an eig det U₁ = 0, a conjugate point, while an eigenvalue +1 marks a singular U₂, a harmless pole of W⁻¹ (Arnold 1967, Funct. Anal. Appl. 1, 1; Robbin & Salamon 1993, Topology 32, 827). Ψ is invariant under U → U·T for any invertible T, so Gaussian reduction leaves it unchanged and the frame can be passed with normalized -columns (see [`normalized_frame`](@ref)). When U₁†U₂ drifts from Hermitian, Ψ is not exactly unitary; the +columns (see `normalized_frame`). When U₁†U₂ drifts from Hermitian, Ψ is not exactly unitary; the eigenpairs are used regardless and the defect reports the drift. """ function cayley_eigen!(θ::AbstractVector{Float64}, V::AbstractMatrix{ComplexF64}, u1::AbstractMatrix{ComplexF64}, u2::AbstractMatrix{ComplexF64}, s::Real)