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>
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
lar_reversed_shear_growth_rate_n1(harness @ c2d6aff)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.
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