Repository navigation
Conversation
…te solve 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 <noreply@anthropic.com>
…s and growth rate 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 <noreply@anthropic.com>
…e docs Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
…secant 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 <noreply@anthropic.com>
…wthRateResult 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 <noreply@anthropic.com>
… the kinetic profiles 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 <noreply@anthropic.com>
… 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 <noreply@anthropic.com>
…oup and a harness case 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 <noreply@anthropic.com>
…and keep the inertial solves on the marginal grid 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 <noreply@anthropic.com>
…sing count and pin the missed-count cases 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 <noreply@anthropic.com>
…yley eigenphases of the solution frame 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 <noreply@anthropic.com>
… 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 <noreply@anthropic.com>
…owth-rate search 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 <noreply@anthropic.com>
5 tasks
…holds 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 <noreply@anthropic.com>
Collaborator
|
@d-burg this needs plots or an artifact containing plots please |
…ith 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 <noreply@anthropic.com>
logan-nc
reviewed
Nov 12, 2026
logan-nc
left a comment
logan-nc
left a comment
Collaborator
There was a problem hiding this comment.
@d-burg this is a good start. The algebra holds up, and the count-based bracket is the right design. I went through the diff at 98fa97b and re-derived the pieces by hand against Glasser, Phys. Plasmas 23, 072505 (2016), Wilson et al., Phys. Plasmas 9, 1277 (2002) and Snyder et al., Phys. Plasmas 16, 056118 (2009). Nothing below was run. These pieces check out:
- the identity ξ×B = ξs∇ψ + Ξu and its sign;
- the metric cofactors and B²;
- the normalization, with no stray 2π;
- the A/C/H bands acting on the true ξs through the kinetic change of variable;
- the ξs elimination redone with A + γ²AK;
- the cylinder formula.
The physics comes first below, because it decides how γ can be described and used. These are also the points Glasser would most likely raise against a growth rate from DCON.
Physics
-
The number this computes is the Γ = 0 growth rate, and that is an upper bound on the ideal-MHD rate, not the rate itself.
- DCON has already minimized δW over ξ∥. Dropping ξ∥ inertia on top of that is exactly ideal MHD with Γ = 0: δW then contains no ξ∥, and the Rayleigh quotient is maximized at ξ∥ = 0. So γ is an exact eigenvalue of a well-posed self-adjoint problem, which is good.
- For Γ = 5/3 it is an upper bound, because δWfull ≥ δW⊥ and Kfull ≥ K⊥. The two share only the marginal point.
- For γ ≫ cs/qR the gap is O(β). For γ ≪ cs/qR, incompressibility needs parallel return flow, and its inertia raises the effective mass by the Pfirsch–Schlüter factor. Then γ⊥/γideal → √(1+2q²) at large aspect ratio, which is √3 for a q = 1 internal kink.
- Near marginal stability every mode is slow, so this is where the gap is largest. Wilson 2002 says this just before Eq. 9: the model is "valid at marginal stability".
- The zero-β cylinder benchmark cannot see any of this.
- Requested: (a) docstrings, the release note and the HDF5
long_nameshould say something like "growth rate with perpendicular inertia only (Γ = 0); an upper bound on the ideal-MHD growth rate". (b) Add a toroidal internal-kink known-answer test on a TJ case with q0 < 1; the code should return √3 × the Bussac rate. - Caveat on the factor: I have √(1+2q²) from our own derivation and from memory of Glasser–Greene–Johnson and Bussac. It should be pinned to a source before it goes in a docstring.
-
The count in
faster(h) = h.nzero + (h.lambda < 0)is the right verdict.- It is the Morse index of δW + γ²K with a free end point: the number of modes faster than γ equals the number of conjugate points plus the number of negative eigenvalues of WP(γ²) + WV.
- Your own qa = 2.9, p0 = 0.8 TJ table shows why the sign of λmin alone is not enough. At γ² = 7.5e7–1e8, λmin > 0 while a mode at γ² = 4.9155e9 exists.
- In a pedestal this is the normal case. High-n ballooning modes are nearly fixed-boundary, so their free-boundary root sits just above a pole.
- Requested:
- Count the negative eigenvalues instead of using
(λ_min < 0), so the integer is the full index. This does not change which root the solver finds, but it matters when the count itself is the verdict. That is what an ELITE/EPED replacement would use: "is anything growing faster than ω*/2?" with no root solve. - State the theorem and its condition (K > 0 on (ξs, Ξ)) in
docs/src/stability.md, with a citation. This is Morse / σ-stability (Goedbloed & Sakanaka, Phys. Fluids 17, 908 (1974)), not Glasser's γ = 0 theorem. - Add a regression test that the fixture above returns "unstable" at γ² = 9e7.
- Count the negative eigenvalues instead of using
-
The small-γ limit is non-analytic, and the √γ² test only holds at zero pressure.
- The inertial layer at a rational surface sets cL/cS ∝ γ2μ, so WP(γ²) − WP(0) ∝ γ2μ with μ = √(−DI).
- μ = ½ at zero β gives your linear-in-γ result. In a pedestal the exponent moves. It goes to zero as a surface approaches Mercier marginality, and for DI > 0 roots pile up at γ = 0.
- The
8 < ratio < 12window runs on a low-β fixture, so it confirms the exponent only where it is trivially 1. - Requested: fit d ln[λmin(γ²) − λmin(0)]/d ln γ on a finite-β case and compare it with 2√(−DI) from
resist_eval_all!at the dominant surface. - The same physics means that near marginality γ ∝ (−δW)1/2μ. The layer half-width, about γτAq in |m − nq|, quickly drops below ρi.
- So the small-γτA numbers are the most model-dependent output of the feature, including the 0.023 in the n = 8 example. It would help to write the layer widths, in metres, next to γ.
-
ω*pi and the EPED1 threshold are off by a factor of two.
- Snyder et al. 2009, Sec. II, gives the threshold as γ > ω*pi/2, with ω*pi taken at its half-maximum value across the edge barrier.
_peak_omega_star_pireturns the full maximum over [psilow, psilim], and the docstring then statesgamma > omega_star_pi/2. The threshold as written is therefore twice EPED1's.- Because the maximum is over the whole domain, a core gradient can set it, and a psilim inside the pedestal foot can cut it off.
- Requested: either return the profile and let the caller choose, or take the maximum over a stated pedestal window and apply the half-maximum rule explicitly. Add a test on a tanh pedestal with a steeper core gradient.
- Separately, γ ∝ ρ−½, so uniform versus profile density can move γ/γ0 by up to about 2 in a pedestal. That deserves a switch and a note on which convention ELITE uses.
-
For the record, though not a problem in this PR (the guard against combining with a kinetic source is right): any later kinetic growth rate needs the KineticForces resonances evaluated at ω = iγ inside the root search.
- At γ ~ 1e4–1e5 s⁻¹, the static ω = 0 matrices are the wrong functional.
- Once the matrices are non-Hermitian, the monotone count and the bracket no longer hold.
Code
-
growth_ratedoes not setsingfac_min, so the inertial solves inherit the default 1e-4.- The kinetic surface finder then makes a single trapezoid crossing wherever cond(F̄) > 1e8. That happens once γτA falls below roughly 1e-3 to 1e-4, near marginality.
- The "Solve path" bullet in the description says the opposite; the adversarial-review section agrees with the code.
- The only evidence that the crossing is accurate is one point at γ² = 1e4.
- Requested: run the inertial solves as a single chunk, as described; or record
kmsingper solve and mark the result not converged when the root used a crossing. Also comparesingfac_min = 0against the default at γτA ≈ 1e-4.
-
nzero_conj = abs(odet.nconj)can cancel silently.- With a Lagrangian frame every crossing goes the same way, so a +1 and a −1 can only come from a miscounted step. Taking the absolute value hides that.
- Requested: accumulate the up and down wraps separately and flag any solve where both are non-zero.
- The ideal frame is the one that loses the Lagrangian property on Solovev. Consider deciding stability from an inertial single-chunk solve at
gamma_sq_minand keeping the ideal count as a diagnostic. - It would also help to report the flag counts and
maslov_defectfor the DIII-D-like marginal solve.
-
The
ftol * scalestop can returnconverged = truewith a wide bracket.- When the marginal λ is positive,
scaleis |λmin(lo)|. λmin → −∞ just above a pole, so if the bisection leaveslonear one, the first secant iterate can pass the test. - Requested: require the bracket to meet
rtolas well, and store the final [lo, hi] on the result.
- When the marginal λ is positive,
-
CK and the non-orthogonal part of HK are untested.
- In the cylinder they enter γ² only at O(ε²). For up-down symmetric equilibria CK has no diagonal and is purely imaginary, so a sign or conjugation error would pass every current fixture. My derivation agrees with your signs, but no test would notice if they were wrong.
- Cheap checks:
- the pointwise identity |∇ψ×u|² = B², i.e.
integrands[1]*integrands[3] - integrands[2]^2 ≈ (jac*mu0rho)^2/b2; - joint positive-definiteness of [AK CK; CK† HK];
- the same γ across
jac_type= hamada / pest / boozer / equal_arc on the DIII-D-like case.
- the pointwise identity |∇ψ×u|² = B², i.e.
- Truncation in m is a Ritz restriction, so γ converges from below. An under-resolved run looks more stable, which is the wrong direction for a stability tool. Please show m-convergence on that case too.
-
Smaller items:
- A root below
gamma_sq_mincomes back as γ = 0 withconverged = true, so it cannot be told apart from stable. An "unstable below floor" flag would fix that. - Reject
kinetic_factor <= 0forPerpendicularInertia; the LU path has no definiteness check. - The γ² controls (1e4 / 1e9 / 1e16 s⁻²) are dimensional. In units of 1/τA² they would mean the same thing on every device.
- The HDF5 annotation still says "at the lowest n", but the code uses the least-stable mode's n.
docs/src/stability.mddoes not yet describe the model or the solver.
- A root below
Review drafted with Claude from a read of the diff and the cited papers; nothing here was run.
Generated by Claude Code
…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 <noreply@anthropic.com>
…uple The stage result tuple gained a growth_rate entry with the [GrowthRate] deck section; the forward-run test still listed the old fields. Co-Authored-By: Claude <noreply@anthropic.com>
…ussian reduction The reductions only combine columns, so the column norms grew without bound (1e238 at n=25). Each reduction now ends by scaling every column to unit norm, with the scale on the diagonal of fixfac so transform_u! undoes it. W = U2/U1 is unchanged; ca_left/ca_right on the forward path change with the column basis. Co-Authored-By: Claude <noreply@anthropic.com>
…tion records consistent with the column rescaling The fixed-boundary zero-crossing check blends two stored steps before the reductions are undone. A reduction between them now rescales the columns, so the earlier step is brought to the same scales first; without this a genuine crossing could be rejected as noise. The rescaling errors on a non-finite solution column, skips a column whose inverse norm is not finite, and transform_u! warns when the compounded scales fall below the floating-point range. Co-Authored-By: Claude <noreply@anthropic.com>
Collaborator
|
This is the missing piece for RMP studies of the pedestal. Ideal δW alone does not settle pedestal stability: the EPED test compares a growth rate with an ad hoc damping (γ > ω_*pi/2), which is why ELITE has been needed. With this PR, GPEC can do that step itself. The chain it completes:
Together this is a full prediction of the RMP's impact on pedestal stability, from the coil currents to the stability verdict. As far as we know, no one has done that so far. 🤖 Generated with Claude Code |
d-burg
changed the title
github-actions
Bot
added
the
changed-results
This branch has not been deployed
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Release note
solovev_n1leading total energy 7.398220e-3 → 7.362675e-3 (0.48 %, a ~1400× cancellation of plasma against vacuum energy, each of which moves 6e-6),tj_n1toroidal torque andtj_kinetic_calculatedeigenfunction norms 0.03 %, NTV torques ~1e-5, all other energies ≤ 3e-6; the forward-integratorca_left/ca_rightcoefficients change basis; growth rate oflar_reversed_shear_growth_rate_n1unchanged (2.957092e3 s⁻¹) (harness @ 9b43744)GPEC can now compute the ideal-MHD growth rate of the least-stable free-boundary mode, the way ELITE does:
PerpendicularInertia(; rho)is akineticsource carrying the perpendicular (incompressible) plasma inertia,kinetic_factor = γ²scales it, andgrowth_rate(prob, Forward(); inertia)finds the γ² at which the free-boundary energy eigenvalue crosses zero. A[GrowthRate]deck section runs the same solve after the force-free-states stage, with the density taken from the equilibrium's kinetic profiles unlessrhois given, and reports the peak ion diamagnetic frequency ω_*pi for the EPED1 γ > ω_*pi/2 test. Existing ideal, kinetic and resistive paths are untouched.What this adds
src/ForceFreeStates/Inertia.jl). They are evaluated once per problem and reused at every γ².kinetic_factor = γ². With inertia the rational surfaces are regular points for γ² > 0, sosingfac_mindefaults to 0 for the inertial solves and the integration is a single chunk with no Frobenius crossings.growth_rateruns the ideal (γ² = 0) solve and counts the modes growing faster than γ as the interior zero crossingsnzeroplus one when the least-stable free-boundary energy λ_min is negative; that count cannot rise with γ² since K > 0. It brackets the count by decades (down while nothing is faster, up while something still is) and, once no interior crossing is left at the lower end, closes λ_min = 0 with an Illinois (bracketed secant) iteration. Each evaluation is one force-balance solve; the root typically needs 5–6 solves after the marginal one. Returns anIdealGrowthRateResult(γ, γ²,converged, the marginal and the inertial solves, and the full (γ², λ_min, nzero) history). A crossing still present at the lower bracket end is bisected away on the count (see the count section below), a root belowgamma_sq_minis reported as γ = 0, and running out ofmaxiteras a warning withconverged = false.PerpendicularInertia(equil)builds ρ(ψ_N) = μ m_p n_e from the attached kinetic profiles with the SLAYER layer's mass per electron (main ions plus the impurity at the local Z_eff), i.e. the same densityInnerLayer.layer_rhouses.ion_diamagnetic_frequency(equil, n)gives ω_*pi(ψ_N) = n (dp_i/dΨ)/(Z_i e n_i) in the layer convention (n_i ∝ n_e, p_i = n_i T_i, Ψ the poloidal flux per radian, positive for outward-decreasing profiles);IdealGrowthRateResult.omega_star_piis its peak over the integration domain at the lowest n (NaN without kinetic T_i), so the EPED1 criterion isgamma > omega_star_pi/2.[GrowthRate](keys =GrowthRateConfigfields:rhoin kg/m³ with NaN meaning from the kinetic profiles,gamma_sq_init/min/max,rtol,ftol,maxiter) runs the stage after the ForceFreeStates write on the very problem that stage solved, reusing its ideal solve as the marginal one; it needsvac_flag = true, the Forward integrator and no kinetic source, and cannot share a deck with[Match](parse-time errors). The result ismain'sgrowth_rateslot andForceFreeStates/GrowthRate/in gpec.h5:gamma,gamma_sq,omega_star_pi,rhoonpsi, andHistory/{gamma_sq, lambda_min, nzero}, with the section recorded underInput/.docs/src/api.mdandworkflow.mdcarry the section.diagnostics.nzero_conj, withnconj_flagsandmaslov_defect), the Maslov index of the solution frame; it is off by default andgrowth_rateswitches it on for its solves through an internal keyword (no deck key, publicsolveunchanged). Details in the count section below.test/runtests_inertia.jl(core suite, 86 tests): the bands are Hermitian with A and H positive definite and linear in ρ; the metric's B² matches the equilibrium's to 1e-3; inertial solves close as:regularwith no crossings; λ_min(γ²) converges to the ideal value as √γ² (the inertial-layer contribution is linear in γ, checked by a factor-of-ten ratio between γ² = 1e4 and 1e2); the crossing path (singfac_min > 0) and the single chunk agree to 1e-5;growth_rateon a stable fixture returns γ = 0 after the marginal solve alone; on a q_a = 2.9 TJ external kink it converges with nzero = 0 throughout; on the q_a = 2.9, p0 = 0.8 variant with a fixed-boundary mode inside the bracket it reaches the fastest root fromgamma_sq_init = 1e8and finds γ > 0 behind a standoff-0.05 wall withnzero > 0in the history, reportsconverged = falsewith a warning atmaxiter = 1, returns γ = 0 when the root lies belowgamma_sq_min, and rejects a kinetic source and a non-positive density; ρ(ψ) from the kinetic TJ fixture matcheslayer_rho; ω_*pi matches a finite difference of n_e T_i and scales with n; the HDF5 group passes the repo's metadata check.test/runtests_solve_api.jlgains a"deck and API parity"case: deck and API agree on the stable TJ fixture (γ = 0) and on a q_a = 2.9 variant that is free-boundary unstable (γ = 1.35e4 s⁻¹, agreement to 1e-12, identical 8-entry histories), plus the four parse-time guards, and the same kink on the two-pass auto grid (grid_type = "auto",mpsi = 0), where deck and API agree only if the deck handsgrowth_ratethe problem the marginal solve came from.Verification outside the test suite
Uniform-current large-aspect-ratio cylinder, external kink m = 2, n = 1, against the analytic
(γτ_A)² = 2(m − nq_a)(nq_a − m + 1)/q_a² with τ_A = R₀√(μ₀ρ)/B₀ (large-aspect-ratio equilibrium, zero pressure, no wall):
The residual is O(ε²), as expected for a cylinder limit of a toroidal code. Other checks: the metric's B² reproduces the equilibrium B² to 2e-5; λ_min(γ²) is monotone over six decades of γ² on the DIII-D n = 2 example; on the TJ fixture the offset λ_min(γ²) − λ_min(0) falls 6.8e-3 → 2.1e-3 → 6.6e-4 → 2.1e-4 → 6.5e-5 for γ² = 1e6 … 1e2 (√10 per decade). On the DIII-D-like example at n = 8 with ρ = 5e19 m⁻³ deuterium: γ = 5.6e4 s⁻¹, γτ_A = 0.023, 9 solves, 21 min on a laptop.
Benchmark against ELITE 5.8 (in progress, 2026-11-13)
A coarse VARYPED pressure/current scan of a DIII-D pedestal equilibrium, run through ELITE 5.8 and through
growth_rateon the same g-files and profile files, n = 5, 10, 15, 20, 25. Both are normalized the same way (γ over the peak ion diamagnetic frequency, the OMFITgamws16convention). GPEC is q-matched to ELITE's edge q per n and runs on ψ_N from 1e-4 to 0.995; ELITE runs on its usual edge window.Caveats, stated plainly:
U₂/U₁solve for the plasma response loses all accuracy when the frame's column magnitudes span more than a hundred decades and the last Gaussian reduction falls on the final step. Scaling the columns first recovers a Hermitian matrix. The fix is commitbce191e30on this branch; see the next section. The affected scan points have not been rerun with it yet.No refinement of the scan has been done yet. The figure is temporary and diagnostic.
Column-scaled plasma response solve (
bce191e30)plasma_response_matrix(u)scales every frame column to unitU₁norm before formingW = U₂ U₁⁻¹. This leavesWunchanged mathematically and is used by the free-boundary solve, the edge scan and the result assembly.Verification on the two failing high-n inertial solves:
Regression harness, all 28 cases,
eb37c4503(parent) againstbce191e30: 20 cases identical, 8 with rows above their noise floor. No noise floor was changed.tj_n1leading-mode plasma energy 3.771920e-1 → 3.772271e-1 (0.01 %), vacuum energy by the opposite amount, largest shift over all modes 3.03 absolute;lar_reversed_shear_n11.5e-7 on the leading mode and 6e-4 over all modes;diiid_n1,diiid_multi_n,diiid_kinetic_multi_n,solovev_n1andsolovev_kinetic_nuzeroat most 6e-9.diiid_n1toroidal torque moves by 1.6e-12 absolute (1.7e-6 relative).The split moves because the unscaled solve was already inaccurate on the standard ideal examples. Evaluating the same final frame both ways, the non-Hermiticity of
Wis 1.2e-4 unscaled against 5.6e-9 scaled on the TJ ideal example (cond(U₁) = 4e19), 9.0e-6 against 3.2e-7 on the reversed-shear example, and 8.2e-7 either way on Solovev. The difference between the twoWmatrices is off-diagonal in the eigenmode basis (TJ: diagonal 1.8e-7, off-diagonal 0.22), so eigenvalues are unaffected to first order and only eigenvector-dependent quantities move.The Gaussian reductions themselves are rescaled in the next section. The full report below predates this commit (
c2d6aff88againstdevelop); the three commits between it andeb37c4503touch only the growth-rate path (root search, conjugate-point count, output) and its documentation; they were not compared againstdevelopseparately.Column rescaling after each Gaussian reduction (
a0354b860,9b437442b)After every Gaussian reduction the solution columns are scaled to unit norm over both blocks, and the scale is stored on the diagonal of that reduction's record so
transform_u!undoes it together with the triangular part.W = U₂ U₁⁻¹is unchanged mathematically. Without it the column norms grow without bound between reductions (5e238 at n = 25) and the high-n solves loseW.Review. Two independent read-only reviews (numerics; downstream consumers) confirmed the reconstruction algebra, the growth-ratio bookkeeping, the truncation order,
W, δW, the normalized eigenfunctions, the conjugate-point count and the Riccati path. Both found one defect, fixed in9b437442b: the fixed-boundary zero-crossing check blends two stored steps before the reductions are undone, so a reduction between them left the two steps on different column scales and a real crossing could be rejected. The earlier step is now brought to the later step's scales, which restores the pre-rescaling behavior of that check; a unit test fails without the fix and passes with it. The same commit errors on a non-finite column norm and warns when the compounded scales fall below the floating-point range (total growth beyond ~1e308 would underflow the stored solutions nearest the axis on reconstruction).Does the rescaling change the physics or only the rounding? The resonant lead column at an ideal crossing is picked by largest resonant amplitude, which is not scale-invariant, so the rescaling changes the lead column at every crossing. A test that forces the original lead columns back while keeping the rescaling shows that this is not what moves the numbers:
The solve is covariant under column scaling to integration tolerance, not to rounding: any rounding-level difference changes the adaptive step sequence. On the Solovev deck the plasma energy spreads 3e-5 across the three variants and the total energy is a ~1400× cancellation of plasma against vacuum energy, which turns that into −0.5 % / +2.3 %. The 0.48 % move of
solovev_n1is therefore the fixture's own tolerance-level sensitivity. The fixed-boundary count is 0 in all three variants on both decks.Regression harness, all 28 cases,
8279da47a(parent) /a0354b860/9b437442b, exit code 0, no noise floor changed:a0354b860against the parent:solovev_n1Re et[1] 7.398220e-3 → 7.362676e-3 (0.48 %), ep[1] and ev[1] by 6e-6 and 3e-6;tj_n1toroidal torque 5.903723e-15 → 5.905384e-15 (0.03 %);tj_kinetic_calculated‖ξ_α‖ and ‖dξ_ψ/dψ‖ by 0.013 %; NTV torques and energies in the kinetic cases ~1e-5; ODE step counts +3 onsolovev_n1andtj_n1;ca_leftchecksum changes (basis-dependent). Growth rate unchanged.9b437442bagainsta0354b860: 21 cases unchanged, 7 with rounding-level rows (energies ≤ 2e-8 relative, one ODE step fewer onsolovev_n1andtj_n1, the twotj_kinetic_calculatedeigenfunction norms by 0.04 %). The fix commit multiplies each column by the stored inverse norm instead of dividing by the norm, so the applied and recorded scales agree exactly; that last-bit change is amplified by the adaptive stepper as above.Left as they are on
developand noted for follow-up: the zero-crossing check still blends two steps that differ by the unit-triangular part of a reduction;plot_delta_primebuilds its diagonal from forward-integrator coefficients, which are basis-dependent; the harness tracks no fixed-boundary count on the ideal cases.Regression report
Adversarial review (two Opus agents, read-only, after the first push)
The physics reviewer re-derived the inertia bands and found them correct: the metric cofactors, B², the J weighting, the Fourier normalization and band layout, the (ξ_s, ξ_s)/(ξ_s, Ξ)/(Ξ, Ξ) placement of A/C/H, the sign of u (hence of C), the kinetic Schur reduction, and that the pencil's normalization is γ²-independent; the cylinder benchmark excludes any constant factor. Defects found and fixed in the last commit:
mainrebuilt the problem from the re-formed equilibrium, which re-forms again on a third grid: deck and API disagreed there.mainnow passesffs.problem; a parity case on the auto grid covers it.singfac_minto 0, which changes the two-pass grid relative to the marginal solve. They now take every control of the marginal problem (grid identical); crossings from the kinetic surface finder can still appear far below any physical root, where only the sign of λ_min matters, and are documented.maxiterexhaustion was silent (now a warning andconverged = false, also in HDF5); thegamma_sq_minexit returned a never-evaluated nonzero γ with the wrong solve (now γ = 0, marginal); a kinetic source on the problem was silently replaced (now an error, API and deck); ω_*pi used the problem's lowest n (now the least-stable mode's n); missing T_e crashed the ω_*pi evaluation (now NaN, with clear errors from the constructors); the API form wrote an empty deck intoInput/gpec_toml_raw(skipped without a config); the first γ² prepared the problem twice; ρ ≤ 0 or non-finite is now rejected.Not addressed, carried as limitations: per-n ω_*pi for multi-n runs, and ω_*pi's n_i ∝ n_e assumption with Z_eff-aware μ in ρ. The fixed-boundary branch and the count are handled below.
The exact conjugate-point count
Why. The bracket of
growth_rateis provably on the fastest mode only if the number of modes growing faster than the trial γ,nconj + (λ_min < 0), is exact.nconjis the number of conjugate points (fixed-boundary modes fitting inside the ψ interval). DCON's detectornzero(sign change of the smallest-magnitude eigenvalue of W⁻¹ = U₁U₂⁻¹ between stored steps) misses one whenever a stiff harmonic pins a tiny eigenvalue near zero. On the TJ fixture atqa = 2.9, p0 = 0.8a fixed-boundary mode at γ² ≈ 1.2e8 was invisible to it, so a search fromgamma_sq_init = 1e8converged to the slower free-boundary root (6.5e7 instead of 4.9e9) and the plasma behind a standoff-0.05 wall was reported stable.nzerois kept as a diagnostic; the new count sits beside it.What. The forward integration counts conjugate points as the Maslov index of the solution frame (Arnold 1967; Robbin & Salamon 1993). With U = [U₁; U₂] the frame of regular solutions, Ψ = (U₁ − i s U₂)(U₁ + i s U₂)⁻¹ is unitary whenever U is Lagrangian (U₁†U₂ 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 and a column normalization leave it unchanged. At every accepted step the eigenphases of Ψ are computed, each is followed by its eigenvector (the eigenvectors are orthonormal), and a phase whose shortest arc crosses π counts one conjugate point. Result fields:
diagnostics.nzero_conj,nconj_flags,maslov_defect.How it stays exact, each item the fix for a failure seen on the fixtures:
growth_ratewarn and reportconverged = false. A rotation still above π/2 after bisection flags as well.Cost: 10–15 % of an inertial solve. Off by default;
growth_rateswitches it on for its solves through the internal_solve_force_balance(…; count_conjugate_points=true)(no deck key; publicsolveunchanged), andmainturns it on for the marginal solve when a[GrowthRate]section follows.Validation (TJ
qa = 2.9, p0 = 0.8, ρ = 1e-7, n = 1;neg= negative free-boundary eigenvalues; the Morse indexnconj + negmust be non-increasing and continuous across the pole):nzero(DCON)nzero_conjnegSame table with the standoff-0.05 wall (its
negdiffers, the index is again monotone). Flags 0 and defect ≤ 5e-10 at every γ².growth_ratefromgamma_sq_init = 1e8: γ² = 4.9155e9 without the wall and 7.17449e8 with it, both converged in 11 solves; the three cases pinned as@test_brokenat the previous commit now pass as@test.A fixture no W⁻¹-based count can handle, and what it says about
nzero. On the Solovev ideal fixture the frame loses the Lagrangian property after the q = 2 crossing (defect 0.2–4; the raw frame's relative non-Hermiticity stays at 1e-8 because its large columns dominate), with non-monotone phase motion even in the first chunk at tolerances 1e-7 and 1e-9. The new count flags it on 132 steps, as designed. DCON'snzeroreports 6 conjugate points there (7 at 1e-9) on a fixture whose free-boundary energy has eight positive eigenvalues and none negative, hence Morse index 0: a tolerance-dependent false positive of the existing detector, not tracked by the harness. An issue is being drafted separately.Relation to ELITE (Snyder et al. user guide, Wilson 2002, Lunniss 2016). ELITE shoots too: a finite-difference sweep in ψ, an edge matching matrix T(γ²), and a secant iteration on one of its eigenvalues from the user's
gamsqguess, with no documented mode counting, γ² scan or treatment of poles, so it converges to whichever root the guess leads to. The count-based bracket here is therefore a step beyond ELITE. ELITE's inertia (Wilson 2002 Eq. 9) is the perpendicular incompressible form used here. The method and its references are written up indocs/src/stability.mdanddocs/development/references.md.Notes for reviewers
IdealGrowthRateResult; the bareGrowthRateResultis Tearing's dispersion-root type, and the first export clashed with it (precompile warning), hence the name.nzerofalse positive above).lar_reversed_shear_growth_rate_n1(reversed-shear LAR example, uniform deuterium density 3.3e-7 kg/m³): γ = 2957 s⁻¹ after 9 solves, nonzero > 0in the history. Its noise floors are absolute (the harness has no relative option): 1e-5 on γ and 1e-1 on γ² and the γ² history, which is the shift a 1e-10 change in λ_min near the root produces (dλ/dγ² ≈ 1.8e-9 per s⁻² there); 1e-10 on the λ_min history as elsewhere. These are new floors on a new case, not relaxed existing ones; tighten if you prefer.🤖 Generated with Claude Code