Skip to content

ForceFreeStates - PERF! - Close the growth-rate root with a Newton step from the eigenmode - #674

Draft
d-burg wants to merge 5 commits into
feature/growth-rate-inertiafrom
performance/growth-rate-newton-step
Draft

d-burg wants to merge 5 commits into
feature/growth-rate-inertiafrom
performance/growth-rate-newton-step

Conversation

@d-burg

@d-burg d-burg commented Nov 13, 2026 •

Copy link
Copy Markdown
Collaborator

Release note

  • Audience: users
  • Numerical impact: on lar_reversed_shear_growth_rate_n1 only: γ moves by 1.9e-3 s⁻¹ (6e-7 relative, inside the root search's own rtol/ftol) and the root-find history shortens from 9 to 7 entries; the other 27 cases are unchanged (harness @ ddc1ada)
  • Migration: re-baseline any pinned [GrowthRate] output: γ can differ within the root search's rtol/ftol, and the History/ arrays hold fewer entries

growth_rate now needs fewer force-balance solves to find γ. Each inertial solve also yields the slope of its least-stable energy in γ², taken from the eigenmode it already saved, so the root is closed by Newton steps instead of a bracketed secant iteration. The API, the stopping rules (ftol, rtol, maxiter) and the decade bracketing are unchanged.

Stacked on #663 (base feature/growth-rate-inertia). That branch has four local commits not yet on the remote (88a770e76, ffc9d47e6, eb37c4503, bce191e30); this branch sits on top of them, so they appear in this diff until #663's branch is pushed. Only ddc1adab3 belongs to this PR.

What changed

  • Slope from the eigenmode. The edge energy Ξ†WΞ is the stationary value of δW + γ²K over the interior displacement at fixed edge Ξ, so by the Hellmann–Feynman theorem, for the Hermitian pencil W·v = λ·N·v of the free-boundary step,

    dλ_min/dγ² = ∫ (Ξ_s†A_ρΞ_s + 2Re Ξ_s†C_ρΞ + Ξ†H_ρΞ) dψ / (ψ₀² Ξ_a†NΞ_a).

    Ξ is the saved least-stable eigenmode (column 1 of the stored solution, which normalize_eigenfunctions! has already put in the eigenvector basis), Ξ_s = −A⁻¹(BΞ′ + CΞ) the eliminated in-surface displacement (Glasser 2016 Eq. 18) with the inertial A and C, A_ρ, C_ρ, H_ρ the inertia bands per unit γ², and N the power norm at the edge. The kinetic reduction the inertial solves take is exactly the Glasser functional with A, C, H shifted by γ² times the bands (its B-like term bkmat is the same shift written against the factorized B), so the stored xi_s_store is the Ξ_s the formula needs.

  • Quadrature. Simpson's rule over each saved ODE step, with the midpoint Ξ and Ξ′ from the cubic Hermite interpolant of the stored Ξ, Ξ′ and one LU of A per midpoint. The trapezoid rule on the saved nodes alone is not good enough: it is 1.3–7 % high on the TJ fixture (58 nodes), where Simpson agrees with a finite difference to 2e-5–3e-4. The cost is a small fraction of one solve.

  • Root iteration. After the decade bracket (and the count bisection, unchanged), each step is a Newton step on γ² from the bracket end with the smaller |λ_min|. A step that would leave the bracket is an Illinois step instead (with its bisection safeguard, as before), and after a Newton step that fails to halve |λ_min| the slope is no longer used. A missing (NaN) or non-positive slope always lands outside the bracket, so it falls back automatically.

Verification

Derivative against a central difference of λ_min(γ²), TJ analytic fixture with q_a = 2.9 (external kink, ρ = 1e-7 kg/m³), relative step 1e-3 (1e-2 gives the same to 4e-6):

γ² (s⁻²) finite difference eigenmode slope ratio
1.8343e8 (the root) 3.26126e-10 3.26168e-10 1.00013
1e9 2.01450e-10 2.01454e-10 1.00002
1e8 3.71104e-10 3.71219e-10 1.00031

Solve counts, including the marginal solve, default rtol = 1e-3, ftol = 1e-6. "Before" is this branch with the slope switched off, which the tests show reproduces the base branch's iteration exactly:

case before after γ before (s⁻¹) γ after (s⁻¹)
TJ fixture, q_a = 2.9 8 6 13543.6306 13543.6305
TJ fixture, q_a = 2.9, p0 = 0.8 12 6 70110.8000 70110.8003
LAR reversed shear (harness case) 9 7 2957.09182 2957.08995

Uniform-current cylinder, external kink m = 2, n = 1, zero pressure, no wall, against (γτ_A)² = 2(m − nq_a)(nq_a − m + 1)/q_a², τ_A = R₀√(μ₀ρ)/B₀:

ε = a/R₀ q_a solves before solves after γτ_A before γτ_A after analytic residual |λ_min|/|λ_min(0)| before → after
0.10 1.600 6 5 0.431202 0.431206 0.43303 1.7e-5 → 5.7e-8
0.10 1.300 6 5 0.491194 0.491198 0.49852 1.8e-5 → 2.0e-8
0.10 1.800 6 5 0.314443 0.314444 0.31429 9.7e-6 → 4.0e-8
0.05 1.600 5 5 0.432539 0.432539 0.43302 6.2e-7 → 3.1e-10

The cylinder decks are my reconstruction of the #663 benchmark (large-aspect-ratio equilibrium with sigma_type = "wesson", p_sig = 0, beta0 = 0, on the reversed-shear example's numerics), so q_a and the ratios differ slightly from that PR's table. On these cases the old iteration stopped on the bracket width with |λ_min| still ~1e-5 of the marginal value; the Newton iteration stops on ftol, so it returns a root about 300 times closer in one fewer solve. The last row saves nothing: the first Newton step lands within 1e-6 of the root but one more solve is needed to certify it.

Tests

test/runtests_inertia.jl: 103 pass (92 on the base branch plus 11 new). New: the slope against a central difference on the TJ kink fixture to 1 % (measured 1.3e-4; the reference's own error at h = 1e-3 is ~1e-5); the fixture's root in fewer than the previous 8 solves with γ within rtol; at the function level, Newton steps reach the same root as the Illinois iteration in fewer evaluations, a NaN or wrong-signed slope reproduces the Illinois iteration exactly, and a slope 100 times too steep still converges. test/runtests_solve_api.jl: 209 pass.

Regression report

regress --refs origin/develop,feature/growth-rate-inertia,performance/growth-rate-newton-step   (all 28 cases)
"last vs prev" compares this branch (ddc1adab3) with its base (bce191e30); origin/develop (f11d75a0a) fails the growth-rate case because it has no [GrowthRate] section.

tj_n1 : 25 unchanged
diiid_slayer_n1_toroidal : 5 unchanged
gal_resistive_diiid_galinner : 11 unchanged
ggj_ray_q500i : 4 unchanged
solovev_kinetic_nuzero : 14 unchanged
diiid_error_field : 10 unchanged
diiid_slayer_n1 : 22 unchanged
diiid_truncation_n1 : 12 unchanged
tj_kinetic_calculated : 20 unchanged
diiid_slayer_n1_coupled : 7 unchanged
gal_resistive_diiid : 11 unchanged
solovev_n1 : 22 unchanged
ggj_reference : 4 unchanged
tj_kinetic_ntv : 7 unchanged
solovev_multi_n : 15 unchanged
diiid_n1 : 53 unchanged
tj_kinetic_ntv_rates : 2 unchanged
diiid_slayer_n1_chi_perp_e : 22 unchanged
diiid_n1_riccati : 17 unchanged
lar_reversed_shear_growth_rate_n1 : 5 changed, 2 unchanged
lar_ggj_tearing_n1 : 3 unchanged
diiid_multi_n : 19 unchanged
diiid_slayer_n1_winding : 22 unchanged
efit_fixedbdy_separatrix : 5 unchanged
tj_kinetic_multiion : 6 unchanged
diiid_kinetic_multi_n : 8 unchanged
diiid_n1_riccati_precompiled : 17 unchanged
lar_reversed_shear_n1 : 22 unchanged

================================================================
Case: lar_reversed_shear_growth_rate_n1 — Large-aspect-ratio reversed-shear equilibrium, n=1, ELITE-style ideal growth rate with uniform deuterium density
================================================================

Regression Tracking: lar_reversed_shear_growth_rate_n1
=====================================================================================================================================================================
Ref 1: origin/develop  @ f11d75a0a (2026-11-12) (FAILED)
       env: environment unknown (cached before fingerprinting)
  Ref 1 error: in expression starting at /var/folders/wz/x_qd5n8d62s0qj0_kcyyqf780000gn/T/jl_kxaZ29GDgH.jl:5
Ref 2: feature/growth-rate-inertia  @ bce191e30 (2026-11-13)
       env: julia 1.11.6, arm64-apple-darwin24.0.0, manifest f1e0eb38 (pinned), 8 threads/8 BLAS, jit build, warm
Ref 3: performance/growth-rate-newton-step  @ ddc1adab3 (2026-11-13)
       env: julia 1.11.6, arm64-apple-darwin24.0.0, manifest f1e0eb38 (pinned), 8 threads/8 BLAS, jit build, warm
---------------------------------------------------------------------------------------------------------------------------------------------------------------------
Quantity                                  origin/develop  feature/growth-rate-inertia  performance/growth-rate-newton-step  Δ prev             Status
---------------------------------------------------------------------------------------------------------------------------------------------------------------------
growth rate gamma (1/s)                   FAILED          2.957092e+03                 2.957090e+03                         1.867e-03 (0.00%)  ** CHANGED **
growth rate squared gamma^2 (1/s^2)       FAILED          8.744392e+06                 8.744381e+06                         1.104e+01 (0.00%)  ** CHANGED **
root-find gamma^2 history                 FAILED          [9 elem]                     [7 elem]                             N/A                ** CHANGED (length) **
root-find lambda_min history              FAILED          [9 elem]                     [7 elem]                             N/A                ** CHANGED (length) **
root-find nzero history (integer counts)  FAILED          [9 elem]                     [7 elem]                             N/A                ** CHANGED (length) **
total energy Re(et[1])                    FAILED          -2.263762e-02                -2.263762e-02                        0.0e+00            OK
# singular surfaces                       FAILED          2                            2                                    0.0e+00            OK
Runtime (s)                               FAILED          46.6s                        11.9s                                                   --
=====================================================================================================================================================================

Notes for reviewers

  • The slope assumes the start of the integration at psilow does not depend on γ². The finite-difference agreement above bounds that term at the 1e-4 level on the fixture; it is not checked on shaped equilibria with a Frobenius axis start.
  • The formula holds for a Hermitian pencil, which the inertial solves are up to roundoff. A wrong slope cannot give a wrong root, only extra solves, because every step stays inside the sign-change bracket.
  • Open decision, not taken here: stopping when the Newton estimate of the distance to the root, |λ_min|/slope, falls below rtol·γ² would save one more solve per root (5 → 4 on the cylinder), but it is a new stopping rule, so it needs approval.
  • CI: the api jobs fail on test/runtests_result_struct.jl:48, which still expects main's return keys without growth_rate. That failure is inherited from ForceFreeStates - FEATURE! - Add perpendicular inertia and an ELITE-style growth-rate solve #663 (its own CI fails the same way) and is not touched here.
  • The PR title carries ! because the tracked γ of the growth-rate harness case moves above its noise floor; the commit subject does not.

NO MERGE WITHOUT HUMAN REVIEW. Draft, assignee @d-burg, no reviewer requested yet.

🤖 Generated with Claude Code

d-burg and others added 5 commits November 9, 2026 17:18
…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>
…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>
…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>
…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>
…p from the eigenmode

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 <noreply@anthropic.com>
@d-burg d-burg self-assigned this Nov 13, 2026
@github-actions github-actions Bot added changed-results Results move or an interface breaks - read before upgrading perf Same answers, less time or memory labels Nov 13, 2026
@github-actions

Copy link
Copy Markdown
Contributor

This pull request is missing a reviewer.

If you are not ready to name them, mark this pull request as a draft.
docs/development/contributors.md suggests lead developers to ask.
Merging is not blocked here, but no pull request may be merged without human review.

This branch has not been deployed

No deployments
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

changed-results Results move or an interface breaks - read before upgrading perf Same answers, less time or memory

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant