ForceFreeStates.Galerkin - FEATURE! - Solve the Galerkin system by banded Cholesky, with LU fallback

[StuartBenjamin]

Release note

• Audience: users
• Numerical impact: (harness @ c827d4d)
  1. The Galerkin Δ′ is solved by Cholesky wherever the matrix is positive definite (gal_resistive_diiid: PEST3 Δ diagonal 0.40 %, Δ′ matrix norm 2e-5, at cd290c6).
  2. Resonant Frobenius powers are the eigenvalues of each side's M₀ block (ForceFreeStates - BUGFIX! - Use trace-aware Frobenius exponents at resonant surfaces #517). Imaginary Δ′ at resonant surfaces falls to near round-off; real Δ′ and the energies move ≤ 3e-5, Galerkin real Δ′ ≤ 6e-4. The Δ′ matrix is Hermitian with a real diagonal.
  3. A warning flags a resonant block whose trace is not small against α.
• Migration: To reproduce the previous Galerkin results, set gal_solver = "LU". extended_precision_bvp now defaults to false (it was true); set it true for the previous Double64 Δ′ BVP solve. Anything pinned to imaginary Δ′ or to Galerkin alpha must be re-baselined.

The Galerkin Δ′ solve now factors the matrix by banded Cholesky and falls back to LU when it is not positive definite (gal_solver = "auto", the new default, as @matt-pharr suggested). Before, gal_solver = "cholesky" crashed (UndefVarError: pbtrf!).

This branch also carries #517 (resonant Frobenius powers from each side's M₀ block), merged into it, and the review follow-ups: the same α on both sides of a Mercier-unstable surface, extended_precision_bvp back as an opt-in (default false), and a warning when a resonant M₀ block's trace is not small against α.

Regression report

regress --refs 278057d96,c827d4d98 (develop vs this head, all 18 cases; julia 1.11.7, pinned manifest, 2 threads):

Case	What moved
gal_resistive_diiid	Im Δ′ at q = 5: −15182 → 24.9, so the Δ′ matrix norm goes 12442 → 2183. Re diagonal ≤ 6e-4 relative (−2700.25 → −2698.58 at q = 5)
diiid_n1_riccati (and _precompiled)	Im Δ′ at q = 5: 286 → 5e-6. Re diagonal ≤ 1e-7 relative
diiid_n1, diiid_multi_n, diiid_truncation_n1, diiid_error_field	energies, PE and torques ≤ 1e-6 relative; island widths and Chirikov ≤ 0.02 %; ODE steps −6 of 1431
solovev_n1, solovev_multi_n, solovev_kinetic_*	Re energies ≤ 2e-5 relative (Im parts ~1e-5 absolute); NTV ≤ 0.06 %
diiid_slayer_n1	Q, ω, γ ≤ 0.17 %
the rest	unchanged; gal_resistive_pe is N/A on both, as on develop

With regress --refs 278057d96,4b44dcd02,c827d4d98, the three review commits (a968945, 2c02244, c827d4d) change none of these numbers: every case reports c827d4d unchanged against 4b44dcd.

Galerkin Cholesky alone, before #517 was merged (6fc89d0 vs cd290c6)

regress --cases gal_resistive_diiid,gal_resistive_pe --refs 6fc89d021,cd290c659 (develop merge-base vs the Cholesky commit; julia 1.11.7, pinned manifest, 2 threads). gal_resistive_pe is N/A on both refs, as on develop.

Case: gal_resistive_diiid — DIII-D-like, n=1, RDCON outer-region Galerkin Δ′ with rpec coil columns (delta_coil)
================================================================

Regression Report: gal_resistive_diiid
=====================================================================================
Ref 1: 6fc89d021  @ 6fc89d021 (2026-10-09)
       env: julia 1.11.7, x86_64-linux-gnu, manifest fdfdedfa (pinned), 2 threads/2 BLAS, jit build
Ref 2: cd290c659  @ cd290c659 (2026-10-13)
       env: julia 1.11.7, x86_64-linux-gnu, manifest fdfdedfa (pinned), 2 threads/2 BLAS, jit build
-------------------------------------------------------------------------------------
Quantity                 6fc89d021     cd290c659     Diff               Status       
-------------------------------------------------------------------------------------
gal # singular surfaces  4             4             0.0e+00            OK           
gal singular q values    [4 elem]      [4 elem]      0.0e+00            OK           
gal PEST3 Δ diagonal     [4 elem]      [4 elem]      6.140e+01 (0.40%)  ** CHANGED **
||gal Δ′ matrix||        1.243825e+04  1.243849e+04  2.367e-01 (0.00%)  ** CHANGED **
||gal Δ_coil block||     2.498141e+02  2.498168e+02  2.686e-03 (0.00%)  ** CHANGED **
gal D_I per surface      [4 elem]      [4 elem]      0.0e+00            OK           
gal α per surface        [4 elem]      [4 elem]      0.0e+00            OK           
||gal inner-layer Δ||    1.474630e+06  1.474630e+06  0.0e+00            OK           
||gal match cout||       3.061595e-02  3.061584e-02  1.2e-07            OK           
gal match residual       1.901643e-16  2.612202e-16  7.1e-17            OK           
Runtime (s)              226.0s        185.4s                           --           
=====================================================================================
Summary: 3 changed, 7 unchanged

Notes for reviewers

• gal_solver: "auto" tries Cholesky, then LU; "LU" skips Cholesky; "cholesky" errors instead of falling back. The key stays because the TPS wrapper passes it.
• The matrix is always assembled with the full band, so the LU fallback sees the same matrix as before; the Cholesky-only assembly branches are gone. The Dirichlet rows (axis, fixed edge, coil edge) are identity rows with their columns kept; the Cholesky path moves those columns to the right-hand side, which leaves a Hermitian system.
• Cholesky uses only the lower band and the real part of the diagonal. The resonant diagonal entries carry a small imaginary error (up to 5.6e-5 relative on gal_resistive_diiid, 1.6e-7 with ForceFreeStates - BUGFIX! - Use trace-aware Frobenius exponents at resonant surfaces #517), which LU keeps and Cholesky drops: that is the 0.40 % change above.
• zpbtrf's info is checked. Fortran RDCON's gal.f ignores it and solves with a partial factor (fixed in RDCON - BUGFIX - Galerkin factorization failures: Cholesky falls back to LU, a failed LU stops PrincetonUniversity/GPEC#301).
• Unit test runtests_galerkin_cholesky.jl: a banded system with Dirichlet rows and a nonzero prescribed value is solved by Cholesky and matches a dense solve; made indefinite, it falls back to LU with the same answer, and "cholesky" errors.

---

🤖 Generated with Claude Code

[matt-pharr]

Couple cleanup items and one main question, being whether we should do away with the flag altogether and just make every galerkin solve try to use cholesky first and then use LU as a fallback. See the relevant comment.

[d-burg]

#517 was merged into this branch, so I'm reviewing its exponent change here, where it will reach develop.

The exponent change is correct, and it is the paper's definition rather than a departure from it, so that's good. Glasser (2016) defines the powers by det[M₀ − αI] = 0 (Eq. 46). Of α = ±√(−D_I) (Eqs. 48–49) the paper says it could not be proved analytically, only that the roots converge to it as Fourier components are added, and it proposes that agreement as a Fourier-convergence check. So sqrt(−det) was Eq. 46 under an unstated tr = 0 assumption.

From claude: regarding @matt-pharr's question (more right, or only more self-consistent): as far as I can derive it, the resonant block has the form [−κ+σ/2 φ; γ conj(κ)−σ/2] with φ, γ real, so tr = −2i·Im κ. That is a genuine term of the truncated operator, not an assembly slip. Matching to a basis with any other powers makes Δ′ depend on the matching distance, so this is more right for the truncated system. It does hide one thing: the old imaginary Δ′ was the only visible symptom of a non-converged resonant block, so I'd like the trace itself reported. Log Im tr / 2 and α per surface (ideally also next to alpha in the HDF5 output), and warn when |tr| is not small against α. The debug print in Asymptotics.jl (lines 125–126, outside this diff) still prints di next to an α that is no longer √(−di); something like:

        block_det = m0mat[1, 1]*m0mat[2, 2] - m0mat[2, 1]*m0mat[1, 2]
        block_tr = m0mat[1, 1] + m0mat[2, 2]
        msg *= @sprintf("  det= %+.12e, tr= %+.12e %+.12ei\n", real(block_det), real(block_tr), imag(block_tr))
        msg *= @sprintf("  alpha= %+.12e %+.12ei\n", real(alpha[1]), imag(alpha[1]))

Requests before this goes to develop:

1. Release note. This PR now carries #517's results-moving change and the removal of extended_precision_bvp, but the release note covers only gal_solver. Suggested Migration: "To reproduce the previous Galerkin results, set gal_solver = "LU". Remove extended_precision_bvp from [ForceFreeStates] (the TOML key is ignored with a warning; passing it as a keyword to solve is now an error); the Δ′ BVP always solves in Float64. Anything pinned to imaginary Δ′ or to Galerkin alpha must be re-baselined." The Numerical-impact line needs the exponent change too.
2. CI is red at test/runtests_eulerlagrange.jl:706: the opposed-powers testset still checks bitwise a branch that was removed (suggestion inline).
3. Harness at the current head. The stamp here (cd290c659) predates the #517 merge, and #517's own report was against a merge-base older than the axis-start change (#514). The near-axis surface is the one most exposed to that.
4. Wording: "Hermitian with a real diagonal" rather than "Δ′ must be real", as Matt suggested on #517.
5. Solovev. Solovev is up-down symmetric, so its blocks should be traceless and the exponent change should not move it, yet #517 reported Im Δ′ going from −1.8e9 to −27 there. What is Im tr at its q = 2 surface? If it is round-off level, that change needs another explanation; if not, the matrices have an imaginary part on a symmetric equilibrium.
6. Near-axis surface at n = 4. A residual Im Δ′ of −25 (about 1 % of the diagonal) means the series and the integrated solution disagree by about 1 % there; I'd treat that as the error bar. A cheap diagnostic: run with JULIA_DEBUG=GeneralizedPerturbedEquilibrium, compare Im tr, α and dpsi^(−2α) at q = 1.25 with the other surfaces, then a three-point singfac_min scan (3e-4, 1e-4, 3e-5) to see whether Δ′₁₁ plateaus.
7. Small: using DoubleFloats in src/ForceFreeStates/ForceFreeStates.jl (line 22) is now dead; GGJ has its own using.

[d-burg]

The opposed-powers branch this testset pinned is gone, so its bitwise check now fails CI (line 706). The second half of the testset already checks the general path exactly, so I'd drop the first half rather than relax it to ≈.

Suggested change

	@testset "column scaling stays /dpsi^α when the powers are opposed" begin
	α = ComplexF64(1.5)
	dpsi = 1.3e-4
	big = ones(ComplexF64, 2, 2)
	small = ones(ComplexF64, 2, 2)
	FFS._apply_column_powers!(big, small, -α, α, dpsi)
	pfac = dpsi^α
	big_ref = ones(ComplexF64, 2, 2)
	small_ref = ones(ComplexF64, 2, 2)
	big_ref ./= pfac
	small_ref .*= pfac
	@test big == big_ref
	@test small == small_ref
	@testset "column scaling multiplies each column by dpsi^power" begin
	dpsi = 1.3e-4

[StuartBenjamin]

Done in a968945: the opposed-powers half is dropped, and the general half stays.

[d-burg]

With the traceless branch removed this is the only path, so a real block with complex roots throws DomainError and one with real roots returns Float64. Converting keeps it total and type-stable without changing any complex result.

Suggested change

	α = sqrt(tr * tr / 4 - det)
	α = sqrt(ComplexF64(tr * tr / 4 - det))

[StuartBenjamin]

Done: α = sqrt(ComplexF64(tr * tr / 4 - det)).

[d-burg]

When tr²/4 − det is negative (a Mercier-unstable surface), it is real by construction, so the sign of the root comes from the sign of a zero or round-off imaginary part. Left and right can then label opposite solutions "big"; a standalone copy of this helper does that for a real block and its mirror. The old shared α made the pairing consistent by construction, and it is now gone entirely. Could α be canonicalised the same way on both sides (Re α > 0, or Im α > 0 when Re α is zero to round-off), or the pair computed once per surface and passed to both sides?

[StuartBenjamin]

Done in a968945 and 2c02244. When Re α is round-off, α is taken with Im α > 0, so a Mercier-unstable surface gets the same α, and the same big root, on both sides. A test puts ±1e-17i in the radicand.

[d-burg]

These tests pin the 2×2 algebra but not the claim. The assembled block always has a purely imaginary trace, which none of the inputs here has, and the unit-test fixtures are up-down symmetric, so CI never exercises the new path. I'd add a case with the real block structure (before the column-scaling testset):

    @testset "Hermitian-structured block: shared imaginary shift, real eigenvectors, same powers on both sides" begin
        # Shape the assembly produces: [−κ+σ/2 φ; γ conj(κ)−σ/2] with φ, γ real and κ → −κ on the left side.
        κ, φ, γ = 0.3 + 0.02im, 0.7, 2.5
        m_right = ComplexF64[-κ+0.5 φ; γ conj(κ)-0.5]
        m_left = ComplexF64[κ-0.5 φ; γ -conj(κ)+0.5]
        α_right, big_right, small_right = FFS.resonant_block_exponents(m_right, 1.0)
        α_left, big_left, small_left = FFS.resonant_block_exponents(m_left, -1.0)
        @test α_right ≈ sqrt((real(κ) - 0.5)^2 + φ * γ)
        @test α_left ≈ α_right
        @test imag(big_right) ≈ -imag(κ)
        @test imag(small_right) ≈ -imag(κ)
        @test big_left ≈ big_right
        @test small_left ≈ small_right
        # ξ′ entries of both leading eigenvectors are real, so the big/small cross-Wronskian 2α/φ is real.
        @test abs(imag(-(m_right[1, 1] - big_right) / φ)) < 1e-14
        @test abs(imag(-(m_right[1, 1] - small_right) / φ)) < 1e-14
        @test -(m_right[1, 1] - small_right) / φ + (m_right[1, 1] - big_right) / φ ≈ 2 * α_right / φ
    end

Ideally also one integration-level assertion on an asymmetric fixture: left and right powers agree, and the Riccati Δ′ diagonal is real to a stated tolerance.

[StuartBenjamin]

Added as you wrote it, before the column-scaling testset (a968945).

[d-burg]

I don't think the evidence supports deleting extended_precision_bvp before this reaches develop, only changing its default. It covers only singfac_min = 1e-4; #517's own 2 % miss at singfac_min = 1e-6 (Solovev q₀ = 1.85) is a case where Float64 is measurably wrong, and the code can no longer detect or repair it. A singfac_min scan is also exactly how one would verify the new exponents. The exponent fix and the precision of this LU solve are separate error sources. Could we keep the option as opt-in (default false, which keeps the whole speedup), or replace it with one step of iterative refinement (Float64 LU, residual in Double64) plus a warning when the correction to the read-off entries is large? What do you think?

[StuartBenjamin]

Agreed. It is back as an opt-in, default false (2c02244), so the speedup stays and a singfac_min scan can still check Float64 (point 6 is that case). A test checks that the Double64 solve matches the Float64 Δ′ diagonal on the DIII-D BVP case.

[StuartBenjamin]

Thanks for the careful review. Most of it concerns #517, which is merged into this branch, so I've answered those points here too. Pushed c827d4d (three commits on top of 4b44dcd). Point by point:

1. Release note: updated in the PR body. It now covers ForceFreeStates - BUGFIX! - Use trace-aware Frobenius exponents at resonant surfaces #517's exponent change and extended_precision_bvp.
2. CI: the opposed-powers half of the column-scaling testset is gone, as you suggested. The general half stays.
3. Harness at the head against develop (regress --refs 278057d96,c827d4d98, all 18 cases, julia 1.11.7, pinned manifest, 2 threads):

Case	What moved
gal_resistive_diiid	Im Δ′ at q = 5: −15182 → 24.9, so the Δ′ matrix norm goes 12442 → 2183. Re diagonal ≤ 6e-4 relative (−2700.25 → −2698.58 at q = 5)
diiid_n1_riccati (and _precompiled)	Im Δ′ at q = 5: 286 → 5e-6. Re diagonal ≤ 1e-7 relative
diiid_n1, diiid_multi_n, diiid_truncation_n1, diiid_error_field	energies, PE and torques ≤ 1e-6 relative; island widths and Chirikov ≤ 0.02 %; ODE steps −6 of 1431
solovev_n1, solovev_multi_n, solovev_kinetic_*	Re energies ≤ 2e-5 relative (Im parts ~1e-5 absolute); NTV ≤ 0.06 %
diiid_slayer_n1	Q, ω, γ ≤ 0.17 %
the rest	unchanged; gal_resistive_pe is N/A on both, as on develop

My three commits on top of 4b44dcd change none of these numbers. With regress --refs 278057d96,4b44dcd02,c827d4d98, every case reports c827d4d unchanged against 4b44dcd: 17 cases run, and gal_resistive_pe is N/A on all three refs.

4. Wording: the release note now says "Hermitian with a real diagonal".
5. Solov'ev q₀ = 1.85, n = 1. Im tr/2 is 3.8e-6 at q = 2 (α = 1.89) and −6.2e-6 at q = 3 (α = 0.96), so it is not round-off: the M₀ blocks carry an imaginary part of about 1e-6 of their entries on this up-down-symmetric equilibrium.
   • It does not shrink with mtheta: 3.8e-6, 3.2e-6, 3.3e-6 at q = 2 for mtheta 256, 512, 1024.
   • It changes sign with mpsi (128 → 256: −7e-7 at q = 2, −4.2e-5 at q = 3), so it looks like radial-spline noise in the M₀ entries.
   • An imaginary power of that size, amplified by dpsi^(−2α) with α = 1.89, accounts for the old −1.8e9.
6. DIII-D-like n = 4 (riccati Δ′ example, nn = 4):
   • At q = 1.25, Im tr/2 ÷ α is 2.6e-6, inside the range of the other 19 surfaces (5e-8 to 3e-5), so the trace is not what is special there.
   • Its block has |M₀₁₁| = 27.6 against α = 0.83 (1 to 4 elsewhere), so α comes out of a large cancellation in det.
   • singfac_min scan (3e-4, 1e-4, 3e-5): Δ′₁₁ = −1196.3, −1279.3, −1272.1, a plateau to 0.6 % below 1e-4. Δ′₁₂ is stable to 0.2 %. Δ′₂₁ = −2.92, −0.20, −4.55, not converged.
   • The other diagonals move ≤ 1.5e-3. Their Im grows from 1e-7 to 5e-4 as singfac_min falls, which is Float64 losing digits at small matching distance.
   • This deck gives |Δ′₁₁| ≈ 1270 and Im ≈ −1.5 at the head, not the 2300 and −25 quoted in ForceFreeStates - BUGFIX! - Use trace-aware Frobenius exponents at resonant surfaces #517.
7. using DoubleFloats in ForceFreeStates.jl stays, because the restored opt-in path uses it.

You also asked for the trace to be reported. The debug print now shows det, tr and α. A warning fires when a resonant block's |tr|/2 exceeds 1e-3·|α|. On the n = 4 and Solov'ev decks above the ratio is at most 3e-5, so it stays quiet there.

[StuartBenjamin]

Just human stuart adding here: @jhalpern30 most of these comments refer to ur original PR, I tried my best to respond to d-burg's review, but feel free to check out the responses if you want

[matt-pharr]

@d-burg @StuartBenjamin Following up on point 5 and on what τ = Im tr(M₀)/2 actually is. I ran a few more tests on feyn. In short: τ is numerical, it has two identifiable sources, and both go away with resolution. The eigenvalue exponents are still the right choice.

What τ is. In every run, τ equals Im[(F̄⁻¹K̄)ᵣᵣ]/(n q′) at the rational surface (r = the resonant mode) to within 1e-9. That combination is large away from the surface and should pass exactly through zero at ψₛ, because two terms cancel there: the resonant term (F̄⁻¹)ᵣᵣ·Im K̄ᵣᵣ and the sum over the other modes. τ is whatever is left over from that cancellation.

It isn't geometry. Mirroring the DIII-D-like g-file in Z does not flip the sign of τ (the mirror/original ratios are +0.2 to +3.2), and the up-down symmetric Solovev has the largest τ of the cases I ran. The off-diagonal Δ′ phases, on the other hand, conjugate exactly under the mirror (e.g. Δ′₁₂ −16.57° → +16.57°). So the complex off-diagonals are genuine geometry: Δ′ is Hermitian, and it is real only for up-down symmetric equilibria.

Where the leftover comes from (DIII-D-like deck):

1. Radial interpolation. L̄ and K̄ are splined separately in ψ, and between knots their combination isn't exactly zero at ψₛ. Putting ψₛ exactly on a grid knot cuts τ by 100–10,000× (e.g. q=3: 1.4e-7 → 1.6e-9 at mpsi = 512). This part shrinks with mpsi, but unevenly, depending on where ψₛ falls within its cell.
2. Poloidal mode truncation. Even with ψₛ on a knot, a small residue remains, largest near the edge: at q=6 it is 2.7e-7 relative to Re(F̄⁻¹K̄)ᵣᵣ with delta_m = 8. Widening the mode band removes it. At q=6 it drops about 1000× at delta_m = 16 and reaches about 1e-11 at 32; at q=2 it reaches round-off by 16.

Neither mtheta (256 to 1024) nor etol (1e-10 to 1e-14) changes either part.

Solovev (point 5). Its τ (1e-6 to 1e-4) wanders and flips sign as mpsi and mtheta change, and the Mercier D_I itself drifts with the Solovev grid. So the floor comes from the analytic Solovev equilibrium rather than from the ψ/θ grids. Amplified near the pole, that gives the old −1.8e9, which fits Stuart's numbers.

What it means for this PR. At any finite resolution the discretized operator really does have this small trace. The eigenvalue exponents keep the asymptotic series consistent with it, so Δ′ no longer depends on the matching distance. Forcing tr = 0 would undo that. Nothing suggests Glasser's ±√(−D_I) is a symmetric-only result; it is the converged limit.

On the new warning. τ is only a loose proxy for an under-resolved surface. At mpsi = 128 the q=4 surface has det(M₀) 34 % off the Mercier D_I, yet τ/α ≈ 4e-5, below the 1e-3 threshold. Comparing det(M₀) directly with the Mercier D_I would catch that case. I'd suggest two follow-ups for an issue rather than this PR: that check, and inserting each rational surface into the radial grid as a knot, which removes the interpolation part of τ by construction.