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/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/Inertia.jl b/src/ForceFreeStates/Inertia.jl index 91a669c131..9f26585ae9 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 @@ -197,6 +198,126 @@ struct IdealGrowthRateResult history::Vector{@NamedTuple{gamma_sq::Float64, lambda::Float64, nzero::Int}} end +""" + inertia_band_splines(km::KineticMatrices) -> (A, C, H) + +ψ-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). +""" +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 + 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, lam_lo = x, lam, lam + side == -1 && (fhi /= 2) + else + 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 +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 @@ -211,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 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 @@ -250,7 +374,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) @@ -271,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) @@ -278,10 +405,14 @@ 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) + 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 + 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 @@ -321,33 +452,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, 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/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 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 diff --git a/test/runtests_inertia.jl b/test/runtests_inertia.jl index 368dac2f3a..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 @@ -197,6 +213,42 @@ 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 + # 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 + end + @testset "kinetic data without T_e" begin ktemplate = joinpath(@__DIR__, "test_data", "regression_tj_kinetic_calculated") kdeck = TOML.parsefile(joinpath(ktemplate, "gpec.toml"))