ForceFreeStates - BUGFIX! - Keep the sign of crit so fixed-boundary crossings are detected

[jhalpern30]

Release note

• Audience: users
• Numerical impact: no tracked quantity moves; nzero, the fixed-boundary instability count, changes from always-zero to reporting genuine conjugate-point crossings (harness @ 39bfde6)
• Migration: re-run any fixed-boundary stability analysis whose conclusion rested on nzero; equilibria previously reported stable may now be reported unstable

ForceFreeStates (FFS) reported every equilibrium as fixed-boundary stable. The critical eigenvalue crit was stored as a magnitude rather than as a signed value, so the sign-change test that detects conjugate points could never fire and nzero was always zero. Fixed-boundary instabilities are now detected.

Regression report

Regression Report: diiid_n1
=================================================================================================
Ref 1: 2cb4f643  @ 2cb4f643 (2026-08-21)
       env: julia 1.12.6, arm64-apple-darwin24.0.0, manifest b75b6327 (pinned), 6 threads/6 BLAS
Ref 2: local  @ local (2026-08-21)
       env: julia 1.12.6, arm64-apple-darwin24.0.0, manifest b75b6327 (pinned), 6 threads/6 BLAS
-------------------------------------------------------------------------------------------------
Quantity                                      2cb4f643         local            Diff       Status
-------------------------------------------------------------------------------------------------
total energy Re(et[1])                        8.013943e-01     8.013943e-01     0.0e+00    OK    
total energy Im(et[1])                        1.233500e-04     1.233500e-04     0.0e+00    OK    
plasma energy Re(ep[1])                       -1.348114e+00    -1.348114e+00    0.0e+00    OK    
vacuum energy Re(ev[1])                       2.149508e+00     2.149508e+00     0.0e+00    OK    
vacuum matrix min eigenvalue                  1.873975e-01     1.873975e-01     0.0e+00    OK    
plasma energy (all)                           [35 elem]        [35 elem]        0.0e+00    OK    
vacuum energy (all)                           [35 elem]        [35 elem]        0.0e+00    OK    
total energy (all)                            [35 elem]        [35 elem]        0.0e+00    OK    
ODE steps (saved)                             2655             2655             0.0e+00    OK    
ODE steps (total)                             4738             4738             0.0e+00    OK    
q0                                            1.204212e+00     1.204212e+00     0.0e+00    OK    
q95                                           4.781723e+00     4.781723e+00     0.0e+00    OK    
beta_t                                        1.327024e-02     1.327024e-02     0.0e+00    OK    
beta_n                                        1.372511e+00     1.372511e+00     0.0e+00    OK    
internal inductance li1                       8.842230e-01     8.842230e-01     0.0e+00    OK    
internal inductance li2                       7.080727e-01     7.080727e-01     0.0e+00    OK    
internal inductance li3                       7.304309e-01     7.304309e-01     0.0e+00    OK    
poloidal beta betap1                          6.680738e-01     6.680738e-01     0.0e+00    OK    
poloidal beta betap2                          5.349836e-01     5.349836e-01     0.0e+00    OK    
poloidal beta betap3                          5.518763e-01     5.518763e-01     0.0e+00    OK    
# singular surfaces                           5                5                0.0e+00    OK    
singular psi locations                        [5 elem]         [5 elem]         0.0e+00    OK    
singular q values                             [5 elem]         [5 elem]         0.0e+00    OK    
current beta betaj                            4.236479e-01     4.236479e-01     0.0e+00    OK    
plasma volume                                 1.829472e+01     1.829472e+01     0.0e+00    OK    
plasma current                                1.152130e+00     1.152130e+00     0.0e+00    OK    
mpert                                         35               35               0.0e+00    OK    
npert                                         1                1                0.0e+00    OK    
toroidal field bt0                            2.006573e+00     2.006573e+00     0.0e+00    OK    
wall field bwall                              3.880145e-01     3.880145e-01     0.0e+00    OK    
aspect ratio                                  2.845746e+00     2.845746e+00     0.0e+00    OK    
elongation kappa                              1.708322e+00     1.708322e+00     0.0e+00    OK    
q profile (checksum)                          ed7c21fd61df...  ed7c21fd61df...  identical  OK    
pressure profile (checksum)                   e15550827bf1...  e15550827bf1...  identical  OK    
Mercier D_I profile (checksum)                5a6fcb1c3a97...  5a6fcb1c3a97...  identical  OK    
resistive interchange D_R profile (checksum)  6284a4c9a75a...  6284a4c9a75a...  identical  OK    
ballooning Delta' profile (checksum)          44bf968c25d5...  44bf968c25d5...  identical  OK    
island half-widths                            [5 elem]         [5 elem]         0.0e+00    OK    
Chirikov parameter                            [5 elem]         [5 elem]         0.0e+00    OK    
||resonant area-weighted field||              5.207739e-04     5.207739e-04     0.0e+00    OK    
PE plasma energy                              3.422586e+00     3.422586e+00     0.0e+00    OK    
PE vacuum energy                              3.174509e+00     3.174509e+00     0.0e+00    OK    
PE surface energy                             5.841099e+00     5.841099e+00     0.0e+00    OK    
PE toroidal torque                            -5.062793e-02    -5.062793e-02    0.0e+00    OK    
NTV torque FGAR [N·m]                         5.587060e-01     5.587060e-01     0.0e+00    OK    
NTV kinetic energy dW FGAR [J]                6.903492e-02     6.903492e-02     0.0e+00    OK    
Runtime (s)                                   146.7s           146.4s                      --    
resonant area-weighted field b^r              [5 elem]         [5 elem]         0.0e+00    OK    
=================================================================================================
Summary: 47 unchanged

Notes for reviewers

The mechanism. findmin(abs, evals) returns the pair (min|λ|, index). compute_smallest_eigenvalue took element [1] — the magnitude — so crit was non-negative by construction. check_for_zero_crossings! detects a conjugate point via crit_store[istep] * crit_store[istep-1] < 0, and a product of two non-negative numbers is never negative, so nzero stayed at 0 for every run. The fix takes element [2] and indexes back into evals, which is what the Fortran ode_output_get_crit does: it sorts on -ABS(evalsi) to locate the index, then returns evalsi at that index, keeping the sign.

Why the harness shows no movement. diiid_n1 is fixed-boundary stable, so crit never crosses zero and nzero is 0 both before and after. No tracked case currently exercises a fixed-boundary-unstable equilibrium, which is precisely the regime this bug hid. Suggested follow-up: add a regression case with a fixed-boundary-unstable equilibrium so nzero is actually pinned by the harness. I did not add one here to keep the change minimal, but it would be worth doing.

Why the !. Results move for any user running a fixed-boundary-unstable equilibrium: the reported stability verdict flips from stable to unstable. Nothing moves for stable cases.

Deliberately out of scope. Switching eigvals! to eigen here would additionally expose the null vector at the crossing, which is what one needs to reconstruct fixed-boundary marginal eigenfunctions. That is a capability change, not a bug fix, and belongs in its own PR.

---

Opened as a draft: no human reviewer is named yet. Per docs/development/contributors.md the natural reviewers for ForceFreeStates are @logan-nc, @matt-pharr, or @d-burg — please pick one and mark ready for review.

🤖 Generated with Claude Code

https://claude.ai/code/session_018YkwX4YMDVFaDi1NM9WcJD

[jhalpern30]

@logan-nc randomly came across this doing something completely different trying to understand benchmarking between 3D GPEC and TERPSICHORE. I had Claude look into something for me and it was like "by the way, here's a bug I found in your code"... idk how we missed this one

On an unrelated note, Claude seems to think we can construct the fixed boundary marginal eigenfunctions with a few lines of code change (which is what its final note is about). Not sure if we would ever need this

[logan-nc]

oopsies! Good catch.

IDK why we'd ever need the eigenvector, but good to have the note on record I suppose.