PerturbedEquilibrium - 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 @ 0fc6b4a)
• Migration: Recompute kinetic PerturbedEquilibrium outputs produced before this fix; prior regularized quantities are unreliable.

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

Ran at 0fc6b4a2, 8 cases, origin/develop (d7e6e6fe) vs working tree:

diiid_n1                    Summary: 47 unchanged
diiid_n1_riccati            Summary: 17 unchanged
solovev_n1                  Summary: 22 unchanged
solovev_kinetic_calculated  Summary: 14 unchanged, 3 missing/N/A
solovev_kinetic_nuzero      Summary: 14 unchanged
solovev_kinetic_ntv         Summary:  6 unchanged
solovev_kinetic_multiion    Summary:  6 unchanged
gal_resistive_diiid         Summary: 10 unchanged

136 tracked quantities, 0 changed.

The 3 N/A are the PE quantities this PR adds. They have no baseline because develop's copy of
the deck has no [PerturbedEquilibrium] section, so nothing is produced to compare against.
Values on this branch:

||xi_clebsch_alpha||      6.506410e-04
||dxi_clebsch_psi/dpsi||  4.170592e-03
||Jb_theta_reg||          2.732935e-04

solovev_kinetic_nuzero returning 14/14 unchanged is the direct check that adding
[PerturbedEquilibrium] to the shared deck does not perturb that case's tracked quantities — it
only costs it the extra PE runtime.

Their noise_threshold was re-pinned 1e-12 -> 1e-8: every other extract = "norm" quantity in
the harness uses 1e-6 or 1e-8, and 1e-12 on norms of order 1e-4 is ~1e-9 relative, tight
enough to trip spuriously.

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 develop. Earlier revisions targeted refactor/freeze-fourfitvars (the head branch of
#383), because this fix touches code that #383 renames. #383 has since merged into develop as
d7e6e6fe, so develop is now the correct base and the stack note no longer applies — including
the caveat about placing a numerics-moving change inside a branch claimed to be numerically inert.

[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.

[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

[matt-pharr]

@jhalpern30 @logan-nc how did the harness for kinetic return everything unchanged if it was a bug that changes what the A matrix is for kinetic? Shouldn't that affect delta W?

[jhalpern30]

@jhalpern30 @logan-nc how did the harness for kinetic return everything unchanged if it was a bug that changes what the A matrix is for kinetic? Shouldn't that affect delta W?

The bug in this PR only affected the FieldReconstruction files, so its only the compute_clebsch_displacements outputs that mattered

[matt-pharr]

Ohhh, I understand now. OK, is it ready for merge? It is reviewed so we should get it in if so. Looks like just a toml conflicts.

[jhalpern30]

Yeah it was just waiting on the rest of the stack and then I didn't have a chance to fix the merge conflicts yesterday.

I'll just note that I couldn't get an exact benchmark between the Fortran and Julia on this - I had Claude look into it but it came back saying it couldn't get the ideal case to agree so I left it there. I basically stopped there so it doesn't necessarily mean its broken though

[logan-nc]

Hm. I've done a couple "check that julia gets the same result" runs for an ITPEA benchmarking exercise recently and it was coming up ok... I'd be interested to a figure of how this is changing the clebsch xi profiles fed to KF. You could over plot ideal, kinetic pre-fix and kinetic post-fix even if there is no fortran comparison... but I don't see WHY there wouldn't be a fortran comparison if you toggle the xclebsch_flag in fortran it outputs the exact xi_s that gets fed to PENTRC.

As a reminder, I recommend this in your claude.md which would have probably automated the suggested figure above for you.

## Major Pull Requests

For substantial implementation efforts (roughly anything beyond a small bug fix), generate a self-contained review package in addition to the normal PR description.

The review package should consist of a single self-contained static HTML file that helps both the author and reviewers quickly understand the intent, scope, implementation, and scientific or engineering consequences of the change. Embed figures directly into the HTML whenever practical.

Emphasize visual communication over prose. Include only material that substantially improves understanding beyond what GitHub already provides through the PR description and diff.

Include only the elements that materially improve understanding:
- workflow or architecture diagrams (prefer D2)
- before/after figures illustrating the behavior or results of the change
- plots, metrics, diagnostics, or benchmarks that support the technical claims made by the PR
- concise explanations of the motivation, design decisions, implementation approach, and consequences
- links to relevant source files, commits, generated outputs, or supporting documentation where useful

Select figures, diagrams, and diagnostics that are specific to the work performed in this PR. Avoid generic templates. Ask: *"What visual evidence would help a skeptical reviewer understand and evaluate these changes most quickly?"*

The review package should complement the GitHub review experience, not duplicate it. Do not repeat file lists, diffs, or other information that GitHub already presents effectively. Brief summaries and key result or performance tables are appropriate when they help orient the reader or support the accompanying figures.

The review package is intended both as a development aid for the author and as supporting material for code reviewers. Produce it automatically for major implementation efforts without requiring the user to request it explicitly.

[jhalpern30]

@logan-nc — both of your asks, answered. Short version: you were right that the Fortran comparison was available, and the xclebsch_flag route settled it. Details below.

1. Fortran parity, verified at the source

You don't actually need the run comparison to answer this one — the Fortran source is unambiguous.

GPEC_main/gpec/gpeq.f:117-123:

IF (reg_flag) THEN
   xmp1_mn=xsp1_mn*(singfac**2/(singfac**2+reg_spot**2))
   IF (kin_flag) THEN
      xms_mn=xss_mn*(singfac**2/(singfac**2+reg_spot**2))
   ELSE
      xms_mn=-MATMUL(bmat,xmp1_mn)-MATMUL(cmat,xsp_mn)
   ENDIF

Under kin_flag it is a pure scalar rescale of the stored xss_mn — no factorization, no solve. It can do that because at gpeq.f:74-79 the kinetic branch reads xss_mn=u4%f straight off the DCON spline, so ξ_s already exists and never needs rebuilding from A. The zhetrf/zhetrs inversion at gpeq.f:93-95 lives only inside the non-kinetic ELSE.

And idcon_matrix really is ideal-only — idcon.f:772-774 builds amat purely from metric tensors:

amat(ipert,jpert)=twopi**2*(nn*nn*g22(dm)
     $     +nn*(m1+m2)*g23(dm)+m1*m2*g33(dm))

with no kinetic term and no kin_flag dependence anywhere in the assembly. Corroborating: idcon.f:797 factorizes that same amat with zhetrf and hard-stops on failure — which would break immediately on a non-Hermitian kinetic A. Fortran's A at this site is structurally guaranteed Hermitian. That is exactly why transplanting cholesky!(Hermitian(...)) onto the Julia kinetic A was invalid.

This also confirms the asymmetry is faithful, not an artifact: Fortran's ideal branch feeds the regularized xmp1_mn into the B matmul while leaving xsp_mn unregularized in the C term — matching the Julia ideal branch line for line. So both branches mirror Fortran, not just the one being fixed.

One nuance worth recording: Fortran gates all of this on reg_flag, and with reg_flag=f (gpeq.f:128-129) sets xmp1_mn=xsp1_mn; xms_mn=xss_mn unscaled. Julia has no reg_flag and keys off reg_spot > 0; the two are equivalent (reg_spot=0 ≡ reg_flag=f), but that equivalence is currently implicit.

2. The figure — and the answer is more specific than "18%"

Overplotted ideal / kinetic pre-fix / kinetic post-fix, |ξ^α/χ₁| vs ψ_N, on Solovev_kinetic_calculated_example at reg_spot=0.05. The resonant harmonics barely move; the non-resonant sidebands drop by about an order of magnitude.

m	pre-fix peak	post-fix peak	change
2 (resonant, ψ_s=0.122)	1.7295e-04	1.7293e-04	−0.01%
3 (resonant, ψ_s=0.880)	7.1387e-05	7.0290e-05	−1.54%
−1	9.45e-06	5.44e-07	−94%
−6	7.37e-07	5.26e-08	−93%
−8	3.01e-07	3.08e-08	−90%
−11	1.63e-07	1.50e-08	−91%

All-m L2 change: 13.6%. So the buggy Cholesky path was inflating the non-resonant sideband structure of the ξ^α spectrum fed to PENTRC, while leaving the resonant harmonics essentially intact. Since NTV consumes the whole spectrum, that is where the downstream consequence lives — worth knowing if you are looking at torque rather than at the resonant surfaces.

Also notable: the pre-fix run does not throw PosDefException. Reflecting the lower triangle still leaves a positive-definite matrix, so the corruption is completely silent.

Caveat, stated plainly: the ideal curve sits ~71% away from the kinetic one regardless of the fix, because kinetic_factor=0 changes the entire ForceFreeStates solve (different ξ, 640- vs 600-node ψ grid requiring interpolation), not just the reconstruction step. It is a shape reference, not a target the kinetic curve should converge toward — please don't read it as validation in either direction. Item 1 is the real evidence.

3. Why the earlier Fortran benchmark didn't agree

Worth recording, since it was left as "couldn't get the ideal case to agree." Diffing the two decks turned up several mismatches, of which two would each independently wreck the comparison:

• reg_flag=f in the shipped docs/examples/solovev_kinetic_example — regularization switched off, so Fortran never exercises the reg_spot path being compared at all.
• The kinetic profiles differ ~30x in both n and T. Fortran kinetic.txt is n≈1.52e21 m⁻³ / T≈88 eV core; Julia kinetic.dat is n≈4.9e19 / T≈2700 eV. Both give the same total pressure (2nT·e ≈ 4.27e4 Pa = the Solovev P₀), so an ideal run agrees fine — but collisionality and drift frequencies key off the n/T split, not the product.
• Different drives: Fortran mode_flag=t, mode=1 (single DCON eigenfunction) vs Julia's forcing.dat harmonic. Reconcilable via harmonic_flag=t, cosmn(2)=1e-4, which matches forcing.dat's 1 2 1e-4 0.0 exactly.

Plus psilow 1e-4 vs 1e-1, psihigh 0.9995 vs 0.99999, mpsi 16 vs 64, delta_mlow/mhigh 8 vs 12, mthvac 64 vs 512, ucrit, and D-T mix vs single ion with electron_flag=f. Encouragingly, the entire Solovev block (mr=mz=ma=128, e=1.6, a=0.33, r0=1.0, q0=1.9) is already identical between the two, so the equilibrium was never the problem.

An aligned deck is prepared, but the run comparison is not finished — kinetic DCON was still integrating when I stopped. Given that item 1 answers the question at the source, I'd rather not block on it. Flagging one limit for whenever someone does run it: Fortran's kinetic matrices come from PENTRC (nl=6, nutype=harmonic, nlmda=128, xmax=128) with no exposed Julia counterpart, so some residual disagreement in xi_s is expected for reasons unrelated to this PR.

---

Regression harness at 0fc6b4a2 vs develop (d7e6e6fe): 8 cases, 136 tracked quantities, 0 changed. solovev_kinetic_nuzero returns 14/14 unchanged, confirming the added [PerturbedEquilibrium] section costs that case only runtime.

[jhalpern30]

@logan-nc I am getting a bit out of my depth here physics-wise - what do you think of the plots? My interpretation is that this fix is an actual bugfix, but doesn't always produce drastically different results. And that everything I can find indicates that this is faithful to the Fortran. It also pointed out differences between the Fortran and Julia examples, which I am not sure were intentional or not

[logan-nc]

Isn't the solovev kinetic example pretty garbage for the sake of speeding the harness? I'd suggest using DIII-D for any benchmark. I also don't see fortran overlayed so it's pretty hard to determine which is "right". I am getting a bit confused by "So both branches mirror Fortran, not just the one being fixed." - Does this mean claude is overly proactive and it was all matching fortran fine to start with?

Note that this is a bit a goose chase if we automatically turn reg_flag off for the (self regularizing) kinetic runs anyways... but I am surprised to see the profiles changing far from the rational (reg_spot should mean the "regularization" is only impacting a small region near the rational) - so that does indeed need to be chased down.