Repository navigation
ForceFreeStates - Repair or retire the Frobenius axis start #464
Description
Activity
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).- The Re(a) comparison ties (−1e-19 vs 0.0).
- The |U₁|-component fallback also ties: 0.0018029472842171165 vs …156.
- Roundoff therefore picks the branch, and multi-n vs single-n arithmetic lands on opposite branches.
- Either branch gives U₁ = 1.7e-4 ∓ 1.8e-3i. With U₂ = I, a non-real diagonal U₁ violates U₁†U₂ = U₂†U₁, which the ODE conserves, so W_p carries an anti-Hermitian part to the edge.
n=1 W_p (DIII-D, ψ_low = 1e-4) ‖anti-Herm‖/‖Herm‖ multi-n vs single-n Frobenius start 4.6e-2 9.1e-2 Fixed start ( frobenius_psi_max = 0)1.5e-6 7.4e-9 Implications:
- This doesn't depend on the equilibrium. Any axis-started run with m=0 in range is exposed, single-n included. Multi-n only flips the branch.
- Eigenvalue-based checks don't see it. The lowest W_t still matches Fortran at 1e-4, as ForceFreeStates - BUGFIX! - Use the fixed start for integrations that begin away from the axis #460's table shows.
- A repaired start needs a U₁ that keeps U₁†U₂ Hermitian (real, for a diagonal U₁ with U₂ = I), and a deterministic rule for degenerate m=0. The retire option clears this immediately.
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.
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:m = 0 start W_plasma ‖W − W†‖/2‖W‖ toroidal torque PE plasma energy PE surface energy et[1] Frobenius, branch Im a < 0 (develop today) 2.8e-2 +0.050875 3.422677 5.826684 0.8012318 + 4.1e-5 i Frobenius, branch Im a > 0 2.8e-2 −0.050643 3.422678 5.842136 0.8012319 + 1.2e-4 i fixed start ( frobenius_psi_max = 0), develop0e68a05531.4e-6 +0.000116 3.420570 5.828844 0.8012355 + 8.2e-5 i So 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.
Summary
compute_axis_initimplements 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):
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:
et[1]reports a spurious −1e5 "instability" where Fortran DCON is stable;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.01default:frobenius_psi_max = 0.01)frobenius_psi_max = 0(fixed start)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:
compute_axis_initbuilds one 2×2 Frobenius block per mode from the diagonal entries of F̄⁻¹, K and G, discarding all mode coupling;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_brokenlines 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
grid_type = "ldp"withmpsi = 128, not the currentautodefault; the ψ_low = 0.01 numbers should be re-confirmed onautobefore anyone leans on them.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
frobenius_psi_max); the mitigation, not the repair.crit(ψ)trace comparison and the scan figures.