Repository navigation
ForceFreeStates - Repair or retire the Frobenius axis start #464
Description
Summary
compute_axis_init implements the Frobenius axis start [Glasser 2016 Eq. 51], whose docstring promises a solution behaving as ~ψ_low^(|m|/2). It does not deliver that limit at any starting surface anyone actually uses, and the resulting initial condition is wrong by enough to move — and on one configuration, to destroy — the free-boundary energies.
Two independent investigations have now measured this from different directions. #460 bounds where the start is used; neither it nor the superseded #440 repairs the start itself. This issue tracks the repair.
Evidence
Away from the axis (#460, DIII-D ideal example, against Fortran DCON):
| ψ_low | innermost coupling row, Julia/Fortran norm | lowest W_t, Julia vs Fortran |
|---|---|---|
| 1e-4 | 1.000 | 1e-4 |
| 0.01 | 1.010 | 0.03 % |
| 0.03 | 1.042 | 0.1 % |
| 0.10 | 1.243 | 2.3 % |
| 0.30 | 2.140 | 7 % |
At ψ_low = 0.01 (#440, synthetic diverted DIII-D-like ramp-up, q0 = 3.41, n = 1, no wall, ldp grid mpsi 128, riccati). The returned U₁ diagonal is 0.11…0.23 for m = −14…−3 and +6.1/−9.8/−13.7 for m = 3/4/5 — nowhere near ψ_low^(|m|/2). Consequences on that equilibrium:
- one eigenvalue of W_p passes through zero, and
et[1]reports a spurious −1e5 "instability" where Fortran DCON is stable; - on a clean neighbouring slice the same stiff eigenvalue is +2.76e5 against DCON's +1.61e4, a factor 17.
The second point is the important one: the eigenvalue is wrong everywhere in that scan, and the apparent "pole" is merely where it crosses zero. A threshold on the starting ψ does not address this.
Measured on #460's head (e93df8f) with #440's fixture, which starts at exactly ψ_N = 0.01 and so stays on the Frobenius path under the frobenius_psi_max = 0.01 default:
| setting | et[1] | max eig W_p | negative W_p eigs |
|---|---|---|---|
default (frobenius_psi_max = 0.01) |
−1.29e5 | 4.70e3 | 1 |
frobenius_psi_max = 0 (fixed start) |
+1.504 | 1.6090e4 | 0 |
Fortran DCON on the same deck, same mlow/mhigh/psilim/dmlim: stiff W_p +1.6092e4, no negative eigenvalue.
Suspects
Both named in #440, neither investigated:
- the diagonal 2×2 truncation of A₀ —
compute_axis_initbuilds one 2×2 Frobenius block per mode from the diagonal entries of F̄⁻¹, K and G, discarding all mode coupling; - the regular-branch selection rule — the eigenvector of the 2×2 block is picked by dominant |U₁| component, with a fallback to U₁ = 1, U₂ = 0 when the second component underflows.
Resolution
Either the start delivers its documented asymptotic limit, verified against that limit on an analytic equilibrium, or it is demoted to an opt-in with the docstring's promise removed. The present situation — a start that is documented as asymptotically correct, used by default near the axis, and measurably not — is the part to end.
Once it is repaired, the @test_broken lines in the fixture-based test (to be carried forward from #440 after #460 merges) report an unexpected pass and should be promoted.
Caveats on the numbers above
- ForceFreeStates - BUGFIX! - Default to the DCON fixed-axis initial condition #440's fixture uses
grid_type = "ldp"withmpsi = 128, not the currentautodefault; the ψ_low = 0.01 numbers should be re-confirmed onautobefore anyone leans on them. - ForceFreeStates - BUGFIX! - Default to the DCON fixed-axis initial condition #440's regression table dates from 2026-09-02 and is stale against current
develop. The two-row table in this issue was measured 2026-09-17 on ForceFreeStates - BUGFIX! - Use the fixed start for integrations that begin away from the axis #460's head.
References
- ForceFreeStates - BUGFIX! - Use the fixed start for integrations that begin away from the axis #460 — bounds the Frobenius start by starting ψ (
frobenius_psi_max); the mitigation, not the repair. - ForceFreeStates - BUGFIX! - Default to the DCON fixed-axis initial condition #440 — superseded and closed; carries the full investigation, the DCON
crit(ψ)trace comparison and the scan figures.
Activity
jhalpern30 commentedon Sep 25, 2026
Adding onto this: another failure of the regular-branch selection, at ψ_low = 1e-4 (inside #460's Frobenius range)
This was found while comparing a multi-n DIII-D run (n = 1–3,
psilow = 1e-4) with single-n runs. The multi-n n=1 block of W_p was 9% off the single-n result. All of that difference is in the anti-Hermitian part. The Hermitian parts agree to 2e-7 and the eigenvalues are identical.Cause: for m=0, the 2×2 Frobenius block in
compute_axis_inithas a purely imaginary eigenpair (≈ ±0.0037i).frobenius_psi_max = 0)Implications:
Opinion: as ψ → 0 the m=0 block tends to a Jordan block whose only eigenvector is U₁ = 0, so the returned |U₁| ≈ √(b/c) = 1.8e-3 is a finite-ψ artifact; a full Frobenius repair is research-grade for little gain over the fixed start.
I think making the fixed start the default and keep Frobenius opt-in, at the very least until it has been properly tested and verified.
d-burg commentedon Oct 5, 2026
More evidence on the m = 0 tie, from single-n runs and the regression harness
Measured on
examples/DIIID-like_ideal_example(n = 1, ψ_low = 1e-4, ideal), julia 1.11.9, x86_64, pinned manifest.1. Equilibrium rounding alone flips the branch. Tightening the equilibrium ODE absolute tolerance (#491, a 1e-8-level change) makes the eigen-solver return the m = 0 pair in the other order, and the fallback then picks the other solution: W₀₀ = 57.0 − 552 i becomes 57.0 + 552 i. The resulting W_plasma is the adjoint of develop's (‖W_new − W_old†‖/‖W‖ = 1.9e-4, against ‖W_new − W_old‖/‖W‖ = 5.6 %). Forcing the other branch on unchanged develop code reproduces it.
2. The harness's
PerturbedEquilibrium/Energies/toroidal_torqueis this defect and nothing else. Recomputing −2n·Im(b̃†Λ̃⁻¹b̃)/4 from the saved matrices reproduces the stored value to all digits, and its sign is the branch:frobenius_psi_max = 0), develop0e68a0553So a branch flip shows up in a regression report as a torque sign change, a 0.27 % move in the PE surface energy, and Im(et[1]) switching between 4.1e-5 and 1.2e-4, with no physics behind any of it. The auto-grid torque sign flip in #347 looks like the same thing.
3. What the flip does not touch. On
examples/DIIID-like_riccati_deltaprime_examplethe Δ′ matrix, its raw form and Δ_coil are bit-identical under either branch. The ideal-crossing elimination is not involved (#489 leaves W_plasma unchanged to 1e-12).4. Cost of the fixed start on this case. Against today's default it moves the PE plasma energy by −0.06 %, the surface energy by +0.04 %, et[1] by 5e-6 and the resonant fields by 0.15–0.2 %, and it removes the anti-Hermitian part and the spurious torque. That supports making it the default, as proposed above.
Suggested follow-ups: until the start is repaired or retired, treat
toroidal_torqueondiiid_n1as a defect measure rather than a converged quantity in the regression harness; and if Frobenius stays available, give the degenerate m = 0 case a deterministic rule.