From 88a770e7659768ef52d02733313e284f4f093891 Mon Sep 17 00:00:00 2001 From: d-burg Date: Mon, 9 Nov 2026 16:41:43 -0500 Subject: [PATCH 1/5] ForceFreeStates - BUGFIX - Keep only the bracket-end solves in the growth-rate search MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit growth_rate stored the full solve of every γ² it evaluated, so memory grew by one solve per iterate (3 GB for one n = 5 solve, 21 GB after nine; 75 GB at n = 15). Only a bracket end can be returned, so keep the latest solve on each side of the root and drop the rest. The iterates and the returned result are unchanged. Co-Authored-By: Claude --- src/ForceFreeStates/Inertia.jl | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/src/ForceFreeStates/Inertia.jl b/src/ForceFreeStates/Inertia.jl index 91a669c131..54a58f05fa 100644 --- a/src/ForceFreeStates/Inertia.jl +++ b/src/ForceFreeStates/Inertia.jl @@ -250,7 +250,8 @@ function growth_rate(prob::ForceBalanceProblem, alg::Forward; inertia::Perpendic @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}() + solves = Dict{Float64,ForceFreeStatesResult}() # the latest solve on each side of the root: only a bracket end can be returned + bracket_end = Dict{Bool,Float64}() # faster(h) > 0 → the γ² of that side's retained solve function finish(gamma_sq, converged) r = get(solves, gamma_sq, nothing) return IdealGrowthRateResult(sqrt(gamma_sq), gamma_sq, converged && !flagged[], inertia, config, _peak_omega_star_pi(something(r, marginal)), marginal, r, history) @@ -278,10 +279,13 @@ function growth_rate(prob::ForceBalanceProblem, alg::Forward; inertia::Perpendic kw = copy(kwargs) kw[:kinetic_factor] = gamma_sq r = counted_solve(ForceBalanceProblem(prob.equil, prob.wall, kw)) - solves[gamma_sq] = r 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) + below = faster(history[end]) > 0 # this γ² is below the fastest root + haskey(bracket_end, below) && delete!(solves, bracket_end[below]) + bracket_end[below] = gamma_sq + solves[gamma_sq] = r verbose && @info @sprintf("growth_rate: γ² = %.4e λ_min = %.6e nzero = %d", gamma_sq, history[end].lambda, history[end].nzero) return history[end] end From ffc9d47e6e2a9ea37fd3e4e9671feebe358e8310 Mon Sep 17 00:00:00 2001 From: d-burg Date: Tue, 10 Nov 2026 23:37:47 -0500 Subject: [PATCH 2/5] ForceFreeStates - DOCS - Document the diamagnetic stabilization thresholds and correct the EPED1 factor growth_rate applies no threshold to the ideal growth rate. State the EPED1 and EPED1.6 conventions in terms of gamma and omega_star_pi, their sources and caveats, and replace the docstring statement gamma > omega_star_pi/2 with omega_star_pi/4 (the effective frequency is half the peak). Co-Authored-By: Claude --- docs/src/citations.md | 16 ++++++++++++++++ docs/src/stability.md | 23 +++++++++++++++++++++++ src/ForceFreeStates/Inertia.jl | 9 +++++---- src/ForceFreeStates/InertiaOutput.jl | 2 +- 4 files changed, 45 insertions(+), 5 deletions(-) diff --git a/docs/src/citations.md b/docs/src/citations.md index 014a258393..a17da09914 100644 --- a/docs/src/citations.md +++ b/docs/src/citations.md @@ -68,6 +68,22 @@ Reformulates the DCON eigenvalue problem as a Riccati matrix ODE, enabling paral --- +> P. B. Snyder, R. J. Groebner, A. W. Leonard, T. H. Osborne, and H. R. Wilson, "Development and validation of a predictive model for the pedestal height," +> *Physics of Plasmas* **16**, 056118 (2009). +> DOI: [10.1063/1.3122146](https://doi.org/10.1063/1.3122146) + +The EPED1 pedestal model. Source of the diamagnetic threshold γ > ω_*pi/2 with an effective ω_*pi linear in n, against which the ideal growth rate from `growth_rate` is compared (see the thresholds section of the stability page). + +--- + +> P. B. Snyder, R. J. Groebner, J. W. Hughes, T. H. Osborne, M. Beurskens, A. W. Leonard, H. R. Wilson, and X. Q. Xu, "A first-principles predictive model of the pedestal height and width: development, testing and ITER optimization with the EPED model," +> *Nuclear Fusion* **51**, 103016 (2011). +> DOI: [10.1088/0029-5515/51/10/103016](https://doi.org/10.1088/0029-5515/51/10/103016) + +The EPED1.6 model. Replaces the linear-in-n diamagnetic threshold with a two-segment fit to two-fluid calculations that rolls over at high n. + +--- + ## PerturbedEquilibrium > J.-K. Park, A. H. Boozer, and A. H. Glasser, "Computation of three-dimensional tokamak and spherical torus equilibria," diff --git a/docs/src/stability.md b/docs/src/stability.md index 74875de157..8f9a35a0e4 100644 --- a/docs/src/stability.md +++ b/docs/src/stability.md @@ -150,6 +150,29 @@ Two counts of the conjugate points (fixed-boundary modes fitting inside the ψ i `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. +### Diamagnetic stabilization thresholds + +`growth_rate` returns the ideal-MHD γ and the peak ion diamagnetic frequency `omega_star_pi` (the maximum of |ω_*pi| over the integration domain, at the mode's n). It applies no stability threshold. Pedestal models call a peeling–ballooning mode unstable only when γ exceeds half of an effective diamagnetic frequency, and the choice of that frequency moves the stability boundary, so every boundary quoted from these numbers must name the convention it used. + +| Convention | Unstable when | Ratio in terms of the outputs | +|:--|:--|:--| +| Ideal | γ > 0 | none | +| EPED1 ([Snyder (2009)](citations.md#ForceFreeStates)) | γ > ω_eff/2 with ω_eff half the pedestal maximum of ω_*pi, so ω_eff ∝ n | R₁ = `gamma / (omega_star_pi/4)` > 1 | +| EPED1.6 ([Snyder (2011)](citations.md#ForceFreeStates)) | γ > ω_eff/2 with ω_eff from a two-segment fit to two-fluid calculations: ∝ n up to n_crit, a weaker slope above | R₁.₆ = R₁/√2 for n ≤ n_crit, and R₁/(√2 (n_crit/n + 0.168 (n − n_crit)/n)) above; unstable when > 1 | + +R₁ is the quantity ELITE prints as γ/(ω_*pi,max/4). The EPED1.6 row is written in the form used by the stability-map post-processing of the OMFIT ELITE module, which defaults to n_crit = 6 and is what the conventional pedestal stability diagrams (pressure gradient α against edge current, contour at 1, maximum over n) show. + +The motivation for the second fit is that ω_*pi varies strongly across a pedestal while a high-n mode is radially narrow, so the frequency the mode feels grows more slowly than n; a threshold linear in n over-stabilizes high n. + +Points to keep in view when using either ratio: + +- The two ratios are equal at n ≈ 1.54 n_crit (n ≈ 9 for n_crit = 6). Below that the EPED1.6 form is the more stable one, by √2 at n ≤ n_crit; above it the less stable, by a factor 1.92 at n = 25 and 2.64 at n = 50. The boundary of a map taken as the maximum over n therefore moves, and so does the most-unstable n. +- With the EPED1.6 form the maximum over n is often set by the highest n computed, so the n range is part of the result. +- n_crit is a model parameter, not a constant: the value 6 is a default, and a published DIII-D analysis scales it with the edge safety factor (n_crit ≈ 21.7/q at the peak pressure gradient). +- The factor √2 and the slope 0.168 are quoted from that post-processing; check them against the fit in [Snyder (2011)](citations.md#ForceFreeStates) before relying on them outside a like-for-like comparison. +- `omega_star_pi` comes from the attached kinetic profiles with n_i ∝ n_e and Z_i = 1, and its peak is taken over `[psilow, psilim]`; a truncated domain can cut off part of the pedestal gradient. +- A comparison of two codes is independent of the convention as long as both growth rates go through the same one; the boundary location is not. + 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/Inertia.jl b/src/ForceFreeStates/Inertia.jl index 54a58f05fa..81a4732fe3 100644 --- a/src/ForceFreeStates/Inertia.jl +++ b/src/ForceFreeStates/Inertia.jl @@ -45,8 +45,8 @@ end 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. +convention (positive for the usual outward-decreasing profiles). Pedestal models compare γ with half of an +effective ω_*pi built from this profile; the conventions in use are in the stability documentation. """ function ion_diamagnetic_frequency(equil::Equilibrium.PlasmaEquilibrium, n::Integer; z_i::Real=1) data = InnerLayer._layer_kinetic_data(equil, "ion diamagnetic frequency ω_*pi") @@ -178,8 +178,9 @@ The outcome of [`growth_rate`](@ref). - `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 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. + profiles (NaN without n_e, T_e and T_i). No diamagnetic threshold is applied to `gamma`; with the EPED1 convention (effective + ω_*pi = half the peak) a mode counts as unstable when `gamma > omega_star_pi/4`, and the stability documentation lists the + other conventions. 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 γ² and diff --git a/src/ForceFreeStates/InertiaOutput.jl b/src/ForceFreeStates/InertiaOutput.jl index 9612e60c8c..e409e201b1 100644 --- a/src/ForceFreeStates/InertiaOutput.jl +++ b/src/ForceFreeStates/InertiaOutput.jl @@ -1,7 +1,7 @@ # Writer for an `IdealGrowthRateResult`. Layout under the caller's parent group (inputs live only under Input/): # # ForceFreeStates/GrowthRate/ -# ├── gamma, gamma_sq, converged, omega_star_pi -- growth rate, its square, the stop flag and the EPED1 ω_*pi reference +# ├── gamma, gamma_sq, converged, omega_star_pi -- growth rate, its square, the stop flag and the peak ω_*pi (no threshold applied) # ├── psi, rho -- inertia density on the marginal solve's ψ_N grid # └── History/ -- every (γ², λ_min, nzero) evaluated by the root search From eb37c4503be1d4ab9df925476dff37d9d7f5d5a6 Mon Sep 17 00:00:00 2001 From: d-burg Date: Thu, 12 Nov 2026 11:29:38 -0500 Subject: [PATCH 3/5] ForceFreeStates - BUGFIX - Safeguard the growth-rate root iteration with bisection When the free-boundary energy at one bracket end dwarfs the other, the interpolated point equals the opposite end in floating point and the search stalled there, then stopped on the bracket width. The iteration now lives in _sign_change_root and takes a geometric bisection step when the interpolated point is not strictly inside the bracket or the same end has been kept three times running. Co-Authored-By: Claude --- src/ForceFreeStates/Inertia.jl | 59 ++++++++++++++++++++-------------- test/runtests_inertia.jl | 18 +++++++++++ 2 files changed, 53 insertions(+), 24 deletions(-) diff --git a/src/ForceFreeStates/Inertia.jl b/src/ForceFreeStates/Inertia.jl index 81a4732fe3..0d38bde160 100644 --- a/src/ForceFreeStates/Inertia.jl +++ b/src/ForceFreeStates/Inertia.jl @@ -198,6 +198,37 @@ struct IdealGrowthRateResult history::Vector{@NamedTuple{gamma_sq::Float64, lambda::Float64, nzero::Int}} end +""" + _sign_change_root(lambda_at, lo, hi, flo, fhi; rtol, ftol, scale, maxiter) -> (x, converged, lo, hi) + +Root of `lambda_at` in the bracket `(lo, hi)` with `flo < 0 < fhi`, by the Illinois iteration (the retained +end's value is halved when it is kept twice running) safeguarded by geometric bisection. The interpolated +point is replaced by `sqrt(lo * hi)` when it is not strictly inside the bracket, which happens in floating +point when one end value dwarfs the other, or when the same end has been kept three times running. +Stops when `|lambda_at(x)| <= ftol * scale` or the bracket is narrower than `rtol * hi`. +""" +function _sign_change_root(lambda_at, lo::Float64, hi::Float64, flo::Float64, fhi::Float64; rtol::Real, ftol::Real, scale::Real, maxiter::Int) + x, side, kept = hi, 0, 0 + for _ in 1:maxiter + x = (lo * fhi - hi * flo) / (fhi - flo) + (kept >= 3 || !(lo < x < hi)) && (x = sqrt(lo * hi)) + lam = lambda_at(x) + abs(lam) <= ftol * scale && return x, true, lo, hi + newside = lam < 0 ? -1 : 1 + kept = newside == side ? kept + 1 : 1 + if lam < 0 + lo, flo = x, lam + side == -1 && (fhi /= 2) + else + hi, fhi = x, lam + side == 1 && (flo /= 2) + end + side = newside + hi - lo <= rtol * hi && return x, true, lo, hi + end + return x, false, lo, hi +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, config=nothing, verbose=false) -> IdealGrowthRateResult @@ -213,7 +244,7 @@ zero the plasma is stable and `gamma` is 0. Otherwise the count is bracketed by 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 +(Illinois) iteration, safeguarded by bisection, 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 @@ -326,33 +357,13 @@ function growth_rate(prob::ForceBalanceProblem, alg::Forward; inertia::Perpendic 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) - 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) + function lambda_at(gamma_sq) 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) - side = -1 - else - hi, fhi = gamma_sq, lam - side == 1 && (flo /= 2) - side = 1 - end - if hi - lo <= rtol * hi - converged = true - break - end + return h.lambda end + gamma_sq, converged, lo, hi = _sign_change_root(lambda_at, lo, hi, h_lo.lambda, h_hi.lambda; rtol, ftol, scale, maxiter) 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 diff --git a/test/runtests_inertia.jl b/test/runtests_inertia.jl index 368dac2f3a..2c751a9f3a 100644 --- a/test/runtests_inertia.jl +++ b/test/runtests_inertia.jl @@ -197,6 +197,24 @@ include("h5_metadata_check.jl") @test any(h.nzero > 0 for h in walled.history) end + @testset "root iteration with a dominant end value" begin + # An end value 21 orders above the other puts the interpolated point on the opposite end in floating point + f(x) = x < 1.5e11 ? -1e21 : log(x / 3e11) + x, converged, lo, hi = FFS._sign_change_root(f, 1e11, 1e12, f(1e11), f(1e12); rtol=1e-3, ftol=1e-6, scale=1.0, maxiter=30) + @test converged + @test abs(x / 3e11 - 1) < 2e-3 + # A well-conditioned sign change keeps the plain interpolation steps + nev = Ref(0) + g(x) = (nev[] += 1; x / 4e10 - 1) + x, converged, _, _ = FFS._sign_change_root(g, 1e10, 1e11, g(1e10), g(1e11); rtol=1e-3, ftol=1e-6, scale=1.0, maxiter=30) + @test converged + @test x ≈ 4e10 rtol = 1e-6 + @test nev[] == 3 + # maxiter exhausted is reported, never an end of the bracket as converged + x, converged, lo, hi = FFS._sign_change_root(f, 1e11, 1e12, f(1e11), f(1e12); rtol=1e-12, ftol=0.0, scale=1.0, maxiter=3) + @test !converged && lo < x <= hi && x < 1e12 + 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 bce191e30ab3651a1920075b4002f38091eeb890 Mon Sep 17 00:00:00 2001 From: d-burg Date: Fri, 13 Nov 2026 10:45:13 -0500 Subject: [PATCH 4/5] ForceFreeStates - BUGFIX - Scale the frame columns before the plasma response solve The Euler-Lagrange frame's column magnitudes can differ by well over a hundred decades at the edge of a high-n integration. The unscaled U2/U1 solve then loses all accuracy when taken directly after a Gaussian reduction, giving a complex, non-Hermitian plasma energy. Scaling each column to unit norm leaves W unchanged and restores the solve. Co-Authored-By: Claude --- src/ForceFreeStates/Free.jl | 16 ++++++++++++++-- src/ForceFreeStates/Result.jl | 2 +- 2 files changed, 15 insertions(+), 3 deletions(-) diff --git a/src/ForceFreeStates/Free.jl b/src/ForceFreeStates/Free.jl index 0d971fa07f..08925c196d 100644 --- a/src/ForceFreeStates/Free.jl +++ b/src/ForceFreeStates/Free.jl @@ -57,6 +57,18 @@ function power_norm_matrix!(Nmat::AbstractMatrix{ComplexF64}, jmat::AbstractVect return Nmat end +""" + plasma_response_matrix(u::AbstractArray{ComplexF64,3}) -> Matrix{ComplexF64} + +Return the plasma response matrix `W = U₂ U₁⁻¹` of the Euler-Lagrange frame `u`. Each column is scaled to unit `U₁` norm before the +solve, which leaves `W` unchanged and keeps the solve accurate when the column magnitudes differ by many decades. +""" +function plasma_response_matrix(u::AbstractArray{ComplexF64,3}) + @views u1, u2 = u[:, :, 1], u[:, :, 2] + column_scale = Diagonal([inv(norm(column)) for column in eachcol(u1)]) + return (u2 * column_scale) / (u1 * column_scale) +end + """ normalize_eigenfunctions!(odet::OdeState, wt::AbstractMatrix{ComplexF64}, psio::Float64) -> OdeState @@ -128,7 +140,7 @@ plasma/vacuum/total energies as a [`FreeBoundaryResult`](@ref). # Compute plasma response matrix W = U₂ * U₁⁻¹ wp = zeros(ComplexF64, numpert_total, numpert_total) - @views wp .= (odet.u[:, :, 2] / odet.u[:, :, 1]) ./ equil.psio^2 + wp .= plasma_response_matrix(odet.u) ./ equil.psio^2 wv, vac = compute_scaled_wv(ctrl, equil, intr) @@ -267,7 +279,7 @@ NamedTuple `(total_eigenvalue, plasma_energy, vacuum_energy, vacuum_eigenvalue)` dV_dpsi = equil.profiles.dVdpsi_spline(odet.psifac) # Compute plasma response matrix - @views wp .= (odet.u[:, :, 2] / odet.u[:, :, 1]) ./ equil.psio^2 + wp .= plasma_response_matrix(odet.u) ./ equil.psio^2 # Retrieve raw vacuum matrix from spline and apply singfac analytically at the local q. # Singfac is not pre-applied in the spline (see free_compute_wv_spline) to avoid interpolating diff --git a/src/ForceFreeStates/Result.jl b/src/ForceFreeStates/Result.jl index 2d2affa457..ab58e56ff7 100644 --- a/src/ForceFreeStates/Result.jl +++ b/src/ForceFreeStates/Result.jl @@ -227,7 +227,7 @@ function build_result( wp = if free_energies !== nothing free_energies.wp elseif odet !== nothing - (odet.u[:, :, 2] / odet.u[:, :, 1]) ./ equil.psio^2 + plasma_response_matrix(odet.u) ./ equil.psio^2 else nothing end From ddc1adab3441e087285b6f0871d2d46df4358ff2 Mon Sep 17 00:00:00 2001 From: d-burg Date: Fri, 13 Nov 2026 11:07:14 -0500 Subject: [PATCH 5/5] ForceFreeStates - PERF - Close the growth-rate root with a Newton step from the eigenmode MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The slope of the least-stable free-boundary energy in γ² is the inertial energy of that solve's own eigenmode (Hellmann–Feynman), integrated from the saved solution with the eliminated in-surface displacement, so each solve gives a Newton step at no extra solve. Steps that leave the bracket fall back to the Illinois iteration; the stopping rules are unchanged. Co-Authored-By: Claude --- src/ForceFreeStates/Inertia.jl | 123 +++++++++++++++++++++++++++++---- test/runtests_inertia.jl | 34 +++++++++ 2 files changed, 143 insertions(+), 14 deletions(-) diff --git a/src/ForceFreeStates/Inertia.jl b/src/ForceFreeStates/Inertia.jl index 0d38bde160..9f26585ae9 100644 --- a/src/ForceFreeStates/Inertia.jl +++ b/src/ForceFreeStates/Inertia.jl @@ -199,31 +199,120 @@ struct IdealGrowthRateResult end """ - _sign_change_root(lambda_at, lo, hi, flo, fhi; rtol, ftol, scale, maxiter) -> (x, converged, lo, hi) + inertia_band_splines(km::KineticMatrices) -> (A, C, H) -Root of `lambda_at` in the bracket `(lo, hi)` with `flo < 0 < fhi`, by the Illinois iteration (the retained -end's value is halved when it is kept twice running) safeguarded by geometric bisection. The interpolated -point is replaced by `sqrt(lo * hi)` when it is not strictly inside the bracket, which happens in floating -point when one end value dwarfs the other, or when the same end has been kept three times running. -Stops when `|lambda_at(x)| <= ftol * scale` or the bracket is narrower than `rtol * hi`. +ψ-interpolants of the perpendicular-inertia bands per unit γ² (components A, C and H of `km.kw`), the +kernel of the inertial energy in [`lambda_min_slope`](@ref). """ -function _sign_change_root(lambda_at, lo::Float64, hi::Float64, flo::Float64, fhi::Float64; rtol::Real, ftol::Real, scale::Real, maxiter::Int) +inertia_band_splines(km::KineticMatrices) = ntuple(i -> cubic_interp(km.psi, Series(km.kw[:, :, (1, 3, 6)[i]]); extrap=ExtendExtrap()), 3) + +""" + lambda_min_slope(r::ForceFreeStatesResult, bands) -> Float64 + +Derivative dλ_min/dγ² of the least-stable free-boundary energy of the inertial solve `r`, with `bands` +the inertia bands of [`inertia_band_splines`](@ref). The edge energy Ξ†WΞ is the stationary value of +δW + γ²K over the interior displacement at fixed edge Ξ, so its γ² derivative is the inertial energy of +the solution itself, and for the Hermitian pencil W·v = λ·N·v of the free-boundary step +(Hellmann–Feynman) + + dλ_min/dγ² = ∫ (Ξ_s†A_ρΞ_s + 2Re Ξ_s†C_ρΞ + Ξ†H_ρΞ) dψ / (ψ₀² Ξ_a†NΞ_a), + +with Ξ the saved least-stable eigenmode, Ξ_s = −A⁻¹(BΞ′ + CΞ) its eliminated in-surface displacement +[Glasser 2016 Eq. 18] with the inertial A and C, A_ρ, C_ρ, H_ρ the bands per unit γ², Ξ_a the edge value +and N the power norm of the pencil. The integral is Simpson's rule over each saved step, with the midpoint +Ξ and Ξ′ from the cubic Hermite interpolant of the stored Ξ, Ξ′ (both fourth-order there). The start +of the integration at `psilow` is taken as γ²-independent. NaN when `r` carries no saved solution, no +free-boundary energies or no kinetic matrices. +""" +function lambda_min_slope(r::ForceFreeStatesResult, bands) + sol, kin = r.solution, r.mats.kinetic + (sol === nothing || r.free_boundary === nothing || kin === nothing) && return NaN + nstep = size(sol.du_store, 3) + (nstep >= 2 && nstep == size(sol.u_store, 4)) || return NaN + np = r.numpert_total + amat, bmat, cmat = (Matrix{ComplexF64}(undef, np, np) for _ in 1:3) + xi_mid, dxi_mid, xi_s_mid = (Vector{ComplexF64}(undef, np) for _ in 1:3) + hint = Ref(1) + function density(psi, xi, xi_s) # inertial energy density per unit γ² + bands[1](vec(amat), psi; hint) + bands[2](vec(cmat), psi; hint) + bands[3](vec(bmat), psi; hint) + return real(dot(xi_s, amat, xi_s) + 2 * dot(xi_s, cmat, xi) + dot(xi, bmat, xi)) + end + mode(i) = (@view(sol.u_store[:, 1, 1, i]), @view(sol.du_store[:, 1, i]), @view(sol.xi_s_store[:, 1, i])) + + energy = 0.0 + xi0, dxi0, xi_s0 = mode(1) + d0 = density(sol.psi_store[1], xi0, xi_s0) + for i in 2:nstep + xi1, dxi1, xi_s1 = mode(i) + d1 = density(sol.psi_store[i], xi1, xi_s1) + h = sol.psi_store[i] - sol.psi_store[i-1] + if h > 0 + psi = sol.psi_store[i-1] + h / 2 + @. xi_mid = (xi0 + xi1) / 2 + h * (dxi0 - dxi1) / 8 + @. dxi_mid = 3 * (xi1 - xi0) / (2h) - (dxi0 + dxi1) / 4 + kin.A_spline(vec(amat), psi; hint) + kin.B_spline(vec(bmat), psi; hint) + kin.C_spline(vec(cmat), psi; hint) + mul!(xi_s_mid, bmat, dxi_mid) + mul!(xi_s_mid, cmat, xi_mid, -1, -1) + ldiv!(lu!(amat), xi_s_mid) + energy += h * (d0 + 4 * density(psi, xi_mid, xi_s_mid) + d1) / 6 + end + xi0, dxi0, d0 = xi1, dxi1, d1 + end + + nmat = Matrix{ComplexF64}(undef, np, np) + jmat = Vector{ComplexF64}(undef, 2 * r.mpert - 1) + r.mats.ideal.J_spline(jmat, r.psilim; hint=Ref(1)) + power_norm_matrix!(nmat, jmat, r.mpert, r.npert, r.equil.profiles.dVdpsi_spline(r.psilim)) + xi_edge = @view sol.u_store[:, 1, 1, nstep] + return energy / (r.equil.psio^2 * real(dot(xi_edge, nmat, xi_edge))) +end + +""" + _sign_change_root(lambda_at, lo, hi, flo, fhi; rtol, ftol, scale, maxiter, slope=nothing) -> (x, converged, lo, hi) + +Root of `lambda_at` in the bracket `(lo, hi)` with `flo < 0 < fhi`. With `slope` (a function giving +dλ/dx at a point `lambda_at` has evaluated, the bracket ends included), each step is a Newton step from +the bracket end with the smaller |λ| whenever that lands strictly inside the bracket; after a Newton step +that does not at least halve |λ| the slope is no longer used. The fallback, and the whole iteration +without `slope`, is the Illinois iteration (the retained end's value is halved when it is kept twice +running) safeguarded by geometric bisection: the interpolated point is replaced by `sqrt(lo * hi)` when it +is not strictly inside the bracket, which happens in floating point when one end value dwarfs the other, +or when the same end has been kept three times running. Stops when `|lambda_at(x)| <= ftol * scale` or the +bracket is narrower than `rtol * hi`. +""" +function _sign_change_root(lambda_at, lo::Float64, hi::Float64, flo::Float64, fhi::Float64; rtol::Real, ftol::Real, scale::Real, maxiter::Int, slope=nothing) x, side, kept = hi, 0, 0 + lam_lo, lam_hi = flo, fhi # the end values themselves; flo and fhi carry the Illinois halvings + newton_ok = slope !== nothing for _ in 1:maxiter - x = (lo * fhi - hi * flo) / (fhi - flo) - (kept >= 3 || !(lo < x < hi)) && (x = sqrt(lo * hi)) + base, lam_base = abs(lam_lo) <= abs(lam_hi) ? (lo, lam_lo) : (hi, lam_hi) + x = newton_ok ? base - lam_base / slope(base) : NaN + newton = lo < x < hi # false for a NaN or non-positive slope + if !newton + x = (lo * fhi - hi * flo) / (fhi - flo) + (kept >= 3 || !(lo < x < hi)) && (x = sqrt(lo * hi)) + end lam = lambda_at(x) abs(lam) <= ftol * scale && return x, true, lo, hi + halved = newton && 2 * abs(lam) <= abs(lam_base) + newton_ok = newton_ok && (!newton || halved) newside = lam < 0 ? -1 : 1 kept = newside == side ? kept + 1 : 1 if lam < 0 - lo, flo = x, lam + lo, flo, lam_lo = x, lam, lam side == -1 && (fhi /= 2) else - hi, fhi = x, lam + hi, fhi, lam_hi = x, lam, lam side == 1 && (flo /= 2) end side = newside + if halved # a successful Newton step restarts the Illinois bookkeeping from the end values + flo, fhi, side, kept = lam_lo, lam_hi, 0, 0 + end hi - lo <= rtol * hi && return x, true, lo, hi end return x, false, lo, hi @@ -243,8 +332,11 @@ zero the plasma is stable and `gamma` is 0. Otherwise the count is bracketed by `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, safeguarded by bisection, brings it to zero, stopping when |λ_min| falls below `ftol` times the reference +above it by bisecting on the count. Then only λ_min changes sign in the bracket and Newton steps on γ² bring it to +zero, each from the bracket end with the smaller |λ_min| and with the slope dλ_min/dγ² taken from that solve's own +eigenmode as its inertial energy ([`lambda_min_slope`](@ref)), at no further solve; a step that would leave the +bracket, and every step after a Newton step that did not halve |λ_min|, is a bracketed secant (Illinois) step +safeguarded by bisection instead. The iteration stops 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 @@ -304,6 +396,8 @@ function growth_rate(prob::ForceBalanceProblem, alg::Forward; inertia::Perpendic kwargs[:kinetic] = inertia 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 γ² + bands = inertia_band_splines(kwargs[:kinetic]) + slopes = Dict{Float64,Float64}() # γ² → dλ_min/dγ² of that solve's eigenmode seen = Dict{Float64,Int}() # γ² → history index; each γ² is one solve function evaluate(gamma_sq) gamma_sq = Float64(gamma_sq) @@ -314,6 +408,7 @@ function growth_rate(prob::ForceBalanceProblem, alg::Forward; inertia::Perpendic 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) + slopes[gamma_sq] = lambda_min_slope(r, bands) below = faster(history[end]) > 0 # this γ² is below the fastest root haskey(bracket_end, below) && delete!(solves, bracket_end[below]) bracket_end[below] = gamma_sq @@ -363,7 +458,7 @@ function growth_rate(prob::ForceBalanceProblem, alg::Forward; inertia::Perpendic 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") return h.lambda end - gamma_sq, converged, lo, hi = _sign_change_root(lambda_at, lo, hi, h_lo.lambda, h_hi.lambda; rtol, ftol, scale, maxiter) + gamma_sq, converged, lo, hi = _sign_change_root(lambda_at, lo, hi, h_lo.lambda, h_hi.lambda; rtol, ftol, scale, maxiter, slope=x -> slopes[x]) 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 diff --git a/test/runtests_inertia.jl b/test/runtests_inertia.jl index 2c751a9f3a..86ade13e6d 100644 --- a/test/runtests_inertia.jl +++ b/test/runtests_inertia.jl @@ -166,6 +166,22 @@ include("h5_metadata_check.jl") @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) + # The root is closed by Newton steps: with the Illinois iteration alone this search took 8 solves (7 inertial after the + # marginal one) to γ = 13543.63 s⁻¹, and rtol = 1e-3 on γ² is 5e-4 on γ + @test length(gr.history) < 8 + @test gr.gamma ≈ 13543.63 rtol = 5e-4 + # dλ_min/dγ² from the eigenmode's inertial energy against a central difference of λ_min. At the relative step + # h = 1e-3 the difference's truncation error is ~h²/6 (λ_min is smooth in γ² at this scale, λ″γ²/λ′ ≈ −0.25) and + # the 1e-8 integration tolerance adds ~1e-5, so the 1 % bound is two decades above the reference's own error; + # the measured disagreement is 1.3e-4. + gamma_sq, h = gr.gamma_sq, 1e-3 + inertial_prob(factor, kinetic) = ForceBalanceProblem(kink_equil; nn=1, wall, ffs_kwargs..., singfac_min=0.0, kinetic, kinetic_factor=factor) + km = FFS.kinetic_matrices(inertial_prob(gamma_sq, inertia), forward) + below, at_root, above = (solve(inertial_prob(gamma_sq * f, km), forward) for f in (1 - h, 1, 1 + h)) + slope = FFS.lambda_min_slope(at_root, FFS.inertia_band_splines(km)) + @test slope > 0 + @test slope ≈ (lambda_min(above) - lambda_min(below)) / (2h * gamma_sq) rtol = 1e-2 + @test isnan(FFS.lambda_min_slope(gr.marginal, FFS.inertia_band_splines(km))) # the ideal solve carries no inertia # 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 @@ -210,6 +226,24 @@ include("h5_metadata_check.jl") @test converged @test x ≈ 4e10 rtol = 1e-6 @test nev[] == 3 + # Newton steps from the end values' slopes reach the same root in fewer evaluations + function count_evaluations(; kwargs...) + n = Ref(0) + l(x) = (n[] += 1; log(x / 3e11)) + x, converged, _, _ = FFS._sign_change_root(l, 1e11, 1e12, log(1 / 3), log(10 / 3); rtol=1e-3, ftol=1e-6, scale=1.0, maxiter=30, kwargs...) + return x, converged, n[] + end + x_secant, converged_secant, n_secant = count_evaluations() + x_newton, converged_newton, n_newton = count_evaluations(; slope=x -> 1 / x) + @test converged_secant && converged_newton + @test abs(x_newton / 3e11 - 1) < 1e-3 && abs(x_secant / 3e11 - 1) < 1e-3 + @test n_newton < n_secant + # A slope that is missing or points out of the bracket leaves the bracketed iteration unchanged + @test count_evaluations(; slope=x -> NaN) == (x_secant, converged_secant, n_secant) + @test count_evaluations(; slope=x -> -1 / x) == (x_secant, converged_secant, n_secant) + # A Newton step that stalls hands over to the bracketed iteration: a slope 100 times too steep still converges + x_steep, converged_steep, _ = count_evaluations(; slope=x -> 100 / x) + @test converged_steep && abs(x_steep / 3e11 - 1) < 1e-3 # maxiter exhausted is reported, never an end of the bracket as converged x, converged, lo, hi = FFS._sign_change_root(f, 1e11, 1e12, f(1e11), f(1e12); rtol=1e-12, ftol=0.0, scale=1.0, maxiter=3) @test !converged && lo < x <= hi && x < 1e12