Skip to content

ForceFreeStates - BUGFIX - Count every fixed-boundary zero crossing: inside Riccati chunks, and after a switch of the smallest eigenvalue - #668

Draft
StuartBenjamin wants to merge 2 commits into
OpenFUSIONToolkit:developfrom
StuartBenjamin:fix/riccati-zero-crossings-v2
Draft

StuartBenjamin wants to merge 2 commits into
OpenFUSIONToolkit:developfrom
StuartBenjamin:fix/riccati-zero-crossings-v2

Conversation

@StuartBenjamin

Copy link
Copy Markdown
Collaborator

Problem. With integrator = "Riccati", nzero tests crit only at chunk endpoints, about 8 per singular surface. A crossing inside a chunk is missed, and two crossings in one chunk cancel. On a DIII-D equilibrium that RDCON and Fortran STRIDE find fixed-boundary unstable at n = 1 and 2, the Riccati path reports 0 at both. The forward path, which checks every step, has a separate gap. It confirms a crit sign change with an interpolated midpoint, and that check rejects a real crossing when the smallest eigenvalue switched at the step before.

Fix.

  • Riccati, inside each chunk: each chunk counts the zero crossings of its principal solution (U₁ = 0 at its start) at every ODE step, in the propagator callback. This uses the inertia of U₂ᴴU₁ from an LDLᴴ factorization; crit is evaluated only on the steps where the inertia changes.
  • Riccati, at assembly: the crossings of the axis solution beyond those come from the chunk's end data alone. They are the negative eigenvalues of dir·X_sᴴP₁₂⁻¹X_e, a Morse-index split that holds for F̄ > 0.
  • Riccati, outer plasma: it is still counted from its stored steps.
  • Forward: a sign change of crit counts when the number of negative eigenvalues of W⁻¹ drops. With F̄ > 0, a crossing raises an eigenvalue through 0; a switch of the smallest eigenvalue or a passage through ∞ does not.
Test Result
DIII-D equilibrium 174956 at 1715 ms (draw 12), n = 1 / 2; RDCON and Fortran STRIDE give 2 / 1 Riccati 0 / 0 before, 1 / 1 after; forward 0 / 0 before, 1 / 0 after. At n = 2 the forward solution's crit never changes sign (closest |crit| = 0.2 at ψ = 0.13, the same on develop), so there is nothing for it to count. On the previous develop (728a7e7) both paths gave 1 / 1, and Riccati matched an exact count of the code's own axis solution (eigenphases of (X + iU)(X − iU)⁻¹, 5000 steps per chunk). The remaining n = 1 gap to Fortran comes from small equilibrium differences, not the counting.
test/runtests_riccati.jl, test/runtests_eulerlagrange.jl pass; new: oscillator (per-step count, Morse split, basis and scale independence, backward chunks), a crossing right after an eigenvalue switch (fails on develop), Solovev verdict (Riccati said stable on develop)
Full test suite all files pass (9,662 tests after runtests_equil.jl); its one failure, "Load TokaMaker i-file" (runtests_equil.jl:312), fails on develop too
Regression harness (27 cases, develop vs this branch) all 27 bitwise

The Solovev example is Mercier unstable at q = 2: its crossings pile up log-periodically into that surface, so their number depends on resolution, and the test compares only the verdict.

🤖 Generated with Claude Code

StuartBenjamin and others added 2 commits November 2, 2026 15:48
…ach chunk, not only at chunk endpoints

The Riccati integrator evaluated the fixed-boundary stability criterion only on the stored chunk-endpoint states in the core, so crossings inside a chunk were missed or cancelled (DIII-D 174956 at 1715 ms: 0 at n=1 and n=2 where RDCON finds 2 and 1). Each chunk now counts the crossings of its principal solution at every step of its propagator solve (step_zero_crossings: inertia of U₂ᴴU₁ by Bunch–Kaufman, crit checked only where it changes), and the assembly adds the crossings of the actual solution beyond those from the chunk's end data (crossings_beyond_principal, a Morse index split exact for F̄ > 0). The re-integrated outer plasma keeps its dense check. Same method as the Fortran STRIDE fix.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
… by the inertia of W⁻¹, not an interpolated midpoint

The midpoint test rejects a real crossing when crit changed eigenvalue at the step before. A crossing lowers the number of negative eigenvalues of W⁻¹ (F̄ > 0); a switch of the smallest eigenvalue or a passage through ∞ does not.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

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

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant