PE - BUGFIX! - Regularize stored xi_s in kinetic Clebsch displacements

[jhalpern30]

Release note

• Audience: users
• Numerical impact: all regularized PerturbedEquilibrium quantities from a kinetic run move — xi_clebsch_alpha by 18%, Jb_theta_reg/Jb_zeta_reg and Jxi_theta_reg/Jxi_zeta_reg by order their own magnitude. Ideal runs and reg_spot = 0 are bit-identical. (harness @ f1cd7b7)
• Migration: none

In a kinetic run the Clebsch displacements were computed by Cholesky-factorizing the kinetic A matrix, which is not Hermitian — the upper triangle was silently discarded and the run completed without error, returning wrong regularized displacements and fields. Kinetic PerturbedEquilibrium results produced before this fix should be treated as suspect.

Regression report

Run at f1cd7b7b against the stack parent refactor/freeze-fourfitvars (bc15b720):

Regression Report: diiid_n1
=============================================================================================================
Ref 1: refactor/freeze-fourfitvars  @ bc15b720 (2026-08-18)
       env: julia 1.12.6, arm64-apple-darwin24.0.0, manifest ece8f5d1 (pinned), 6 threads/6 BLAS
Ref 2: local  @ local (2026-08-18)
       env: julia 1.12.6, arm64-apple-darwin24.0.0, manifest ece8f5d1 (pinned), 6 threads/6 BLAS
-------------------------------------------------------------------------------------------------------------
Quantity                                      refactor/freeze-fourfitvars  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)                                   159.7s                       171.3s                      --    
resonant area-weighted field b^r              [5 elem]                     [5 elem]         0.0e+00    OK    
=============================================================================================================
Summary: 47 unchanged

Regression Report: solovev_n1
=============================================================================================
Ref 1: refactor/freeze-fourfitvars  @ bc15b720 (2026-08-18)
       env: julia 1.12.6, arm64-apple-darwin24.0.0, manifest ece8f5d1 (pinned), 6 threads/6 BLAS
Ref 2: local  @ local (2026-08-18)
       env: julia 1.12.6, arm64-apple-darwin24.0.0, manifest ece8f5d1 (pinned), 6 threads/6 BLAS
---------------------------------------------------------------------------------------------
Quantity                      refactor/freeze-fourfitvars  local            Diff       Status
---------------------------------------------------------------------------------------------
total energy Re(et[1])        7.924296e-01                 7.924296e-01     0.0e+00    OK    
total energy Im(et[1])        1.013520e-03                 1.013520e-03     0.0e+00    OK    
plasma energy Re(ep[1])       -9.568089e+00                -9.568089e+00    0.0e+00    OK    
vacuum energy Re(ev[1])       1.036052e+01                 1.036052e+01     0.0e+00    OK    
vacuum matrix min eigenvalue  2.171547e+00                 2.171547e+00     0.0e+00    OK    
plasma energy (all)           [32 elem]                    [32 elem]        0.0e+00    OK    
vacuum energy (all)           [32 elem]                    [32 elem]        0.0e+00    OK    
total energy (all)            [32 elem]                    [32 elem]        0.0e+00    OK    
ODE steps (saved)             384                          384              0.0e+00    OK    
ODE steps (total)             605                          605              0.0e+00    OK    
q0                            1.900006e+00                 1.900006e+00     0.0e+00    OK    
q95                           3.147422e+00                 3.147422e+00     0.0e+00    OK    
beta_t                        4.620628e-02                 4.620628e-02     0.0e+00    OK    
beta_n                        3.215363e+00                 3.215363e+00     0.0e+00    OK    
# singular surfaces           2                            2                0.0e+00    OK    
singular psi locations        [2 elem]                     [2 elem]         0.0e+00    OK    
singular q values             [2 elem]                     [2 elem]         0.0e+00    OK    
mpert                         32                           32               0.0e+00    OK    
npert                         1                            1                0.0e+00    OK    
q profile (checksum)          01d32a2c9bfb...              01d32a2c9bfb...  identical  OK    
pressure profile (checksum)   8c3cdcd12941...              8c3cdcd12941...  identical  OK    
Runtime (s)                   116.7s                       117.9s                      --    
=============================================================================================
Summary: 21 unchanged

Regression Report: solovev_kinetic_calculated
=====================================================================================
Ref 1: refactor/freeze-fourfitvars  @ bc15b720 (2026-08-18)
       env: julia 1.12.6, arm64-apple-darwin24.0.0, manifest ece8f5d1 (pinned), 6 threads/6 BLAS
Ref 2: local  @ local (2026-08-18)
       env: julia 1.12.6, arm64-apple-darwin24.0.0, manifest ece8f5d1 (pinned), 6 threads/6 BLAS
-------------------------------------------------------------------------------------
Quantity                  refactor/freeze-fourfitvars  local          Diff     Status
-------------------------------------------------------------------------------------
vacuum energy Re(ev[1])   1.037963e+01                 1.037963e+01   0.0e+00  OK    
q0                        1.900003e+00                 1.900003e+00   0.0e+00  OK    
||xi_clebsch_alpha||      N/A                          6.666572e-04            N/A   
total energy Re(et[1])    1.861536e+00                 1.861536e+00   0.0e+00  OK    
total energy Im(et[1])    -1.418013e+00                -1.418013e+00  0.0e+00  OK    
plasma energy Re(ep[1])   -8.518098e+00                -8.518098e+00  0.0e+00  OK    
ODE steps (saved)         600                          600            0.0e+00  OK    
# singular surfaces       2                            2              0.0e+00  OK    
singular psi locations    [2 elem]                     [2 elem]       0.0e+00  OK    
total energy (all)        [32 elem]                    [32 elem]      0.0e+00  OK    
ODE steps (total)         735                          735            0.0e+00  OK    
mpert                     32                           32             0.0e+00  OK    
q95                       3.147422e+00                 3.147422e+00   0.0e+00  OK    
||Jb_theta_reg||          N/A                          2.744765e-04            N/A   
Runtime (s)               112.6s                       120.9s                  --    
||dxi_clebsch_psi/dpsi||  N/A                          4.285790e-03            N/A   
singular q values         [2 elem]                     [2 elem]       0.0e+00  OK    
npert                     1                            1              0.0e+00  OK    
=====================================================================================
Summary: 14 unchanged, 3 missing/N/A

gal_resistive_pe was also run: it extracts N/A for all 8 quantities on both refs, i.e. it is
already broken on the parent branch, independent of this change. Flagged below.

Notes for reviewers

What was broken

compute_clebsch_displacements factorized the active A with cholesky!(Hermitian(amat, :L)).
In a kinetic run that A is the kinetic A — amat + kwmat[:,:,1] + ktmat[:,:,1]
(src/ForceFreeStates/Kinetic.jl:145) — which is non-Hermitian; Kinetic.jl itself factorizes it
with lu, and Fortran fourfit.F:1153 comments it ! invert non-hermitian a matrix.
Hermitian(amat, :L) discards the upper triangle, and cholesky! then presumes definiteness.

It fails quietly. On the Solovev kinetic deck the unfixed code exits 0 — no PosDefException —
so this was silent corruption, not a latent crash. Pre- vs post-fix on that deck:

Quantity	rel. change
xi_clebsch_alpha	18%
Jb_theta_reg / Jb_zeta_reg	103% / 101%
Jxi_theta_reg / Jxi_zeta_reg	101% / 102%
xi_clebsch_psi, dxi_clebsch_psidpsi, xi_psi	0 (bit-identical)

Pre-existing, not introduced by the stack: on develop the same line read ffit.amats, which in a
kinetic run was already the kinetic-overwritten A.

The fix, and why not just swap in an LU

gpeq.f:117-123: under kin_flag GPEC applies the scalar regularizer directly to the stored
xss_mn and inverts nothing. It never holds a kinetic A at this site at all — idcon_matrix
(idcon.f:773-777) rebuilds amat analytically from the metric tensors, ideal-only. Our stored
ξ_s is the right analogue: compute_node_xi_s! already builds it with the correct LU on the
kinetic A.

Swapping cholesky! for lu! and keeping the -A⁻¹(B·xmp1 + C·xsp) form would be a Julia-only
deviation rather than a restoration, and the two are not equivalent at finite reg_spot: the
LU form regularizes xsp1 before B multiplies it and leaves C·xsp unregularized, whereas Fortran
regularizes the assembled ξ_s afterward. These do not commute; they agree only as reg_spot → 0.

@logan-nc — this is the main thing I want your read on. Is "regularize the assembled ξ_s" what
PENTRC should be consuming, or does the ordering matter for how you use these downstream?

Coverage

There was none. No example ran kinetic + PE, so this branch had never executed in any example or
regression case. Added [ForcingTerms]/[PerturbedEquilibrium] + forcing.dat to
Solovev_kinetic_calculated_example and three quantities to its regression case. Costs ~8 s there,
and the same on solovev_kinetic_nuzero, which reuses the deck.

Second question for @logan-nc: happy with that, or would you rather have a separate deck to keep
those cases lean?

Out of scope — worth separate issues

1. The ideal path uses cholesky! where Fortran uses zhetrf (Bunch-Kaufman, Hermitian
   indefinite) while deliberately using zpbtrf for fmat in the same routine — and has an
   operator-facing error at idcon.f:800: "zhetrf: amat singular at psi = ..., reduce
   delta_mband". Same answer when A is positive definite; a bare PosDefException instead of that
   actionable message when it is not.
2. gal_resistive_pe extracts N/A for all quantities on both refs — already broken on the parent
   branch.

Base branch

Targets refactor/freeze-fourfitvars (the head branch of #383) rather than develop, because the
fix touches code that #383 renames. Not managed as a stack entry — merge this into
refactor/freeze-fourfitvars directly; #383's own path onward is handled separately.

Note for whoever merges: #383's stated claim is that it is numerically inert. Merging this into that
branch puts a numerics-moving change (kinetic PE only — see the table above) inside it, so that
claim will no longer hold for the combined branch. Ideal-path results stay bit-identical either way.

[jhalpern30]

@logan-nc this was something that came up during #383 that Claude flagged - I had it look into the Fortran source for comparisons and this is what it came up with, which seems legit, but I don't know much of the physics here so I need a pair of kinetically-inclined eyes

[logan-nc]

huh. The more we dig the more skeletons we find, eh?

This diff looks fine and clean to me. Assuming the description is correct in that this matches the fortran (I don't have time to confirm), then yeah we should just merge and match fortran for now. I suppose an issue should be raised asking someone to do the task of benchmarking the impact of correctly implementing the full factorization approach for kinetic cases. Issues are piling up, but that is no reason not to record things like this there.

[github-actions]

This pull request is missing a reviewer.

If you are not ready to name them, mark this pull request as a draft.
docs/development/contributors.md suggests lead developers to ask.
Merging is not blocked here, but no pull request may be merged without human review.

[jhalpern30]

Great, thanks for looking into it. I have Claude looking into benchmarking the kinetic (and also potential ideal issue it flagged) as a side task and will confirm that before I merge anything