PerturbedEquilibrium - BUGFIX - change singular coupling calculations to match Fortran, using du_store instead of splined xi'

[matt-pharr]

Re-implemented du_store saving in ForceFreeStates so b1 can be calculated from this at rational surfaces rather than the spline. This closes #347. Implementation matches the fortran formulation.

Regression table:

Fixed grid (grid_type = "ldp", mpsi = 256)

quantity	0.4.17	0.4.15	diff	(issue, pre-fix)
resonant area-weighted field b^r	[5 elem]	[5 elem]	1.5e-5	24.62%
‖resonant area-weighted field‖	5.220818e-04	5.220814e-04	7e-7	10.59%
Chirikov parameter	[5 elem]	[5 elem]	3.8e-5	8.81%
island half-widths	[5 elem]	[5 elem]	1.2e-5	3.94%
delta prime	[5 elem]	[5 elem]	4.7e-6	~18%
ODE steps (saved)	1481	1466	1.01%	1.01%
PE toroidal torque	−5.060620660e-02	−5.060620658e-02	2.6e-10	1.3e-11
singular psi locations	[5 elem]	[5 elem]	0.0	—

Auto grid (grid_type = "auto", mpsi = 0)

quantity	0.4.17	0.4.15	diff	(issue, pre-fix)
PE toroidal torque	−5.043144946e-02	−5.043144996e-02	9.9e-9	sign flip
resonant area-weighted field b^r	[5 elem]	[5 elem]	2.3e-5	4.45%
Mercier D_I profile (checksum)	3d6d04ae…	577cbb7d…	different	different
ODE steps (total)	1987	1941	2.32%	3.45%
singular psi locations	[5 elem]	[5 elem]	5.5e-15	2e-14

[matt-pharr]

Review — PR #348: singular coupling from du_store instead of splined ξ′

Overview

This PR fixes the two-point b¹ evaluation at rational surfaces (#347) by (1) storing the true ODE right-hand side du/dψ (du_store) at every saved step, (2) reconstructing ξ′ at ψ_s ± spot from it via two new evaluators, and (3) fixing a conjugation bug in the C-coefficient contraction. Regression tables show the previously 4–25% discrepancies vs. Fortran collapsing to ≤4e-5, with PE toroidal torque and singular psi locations essentially unmoved — consistent with a targeted fix.

Correctness — verified against the codebase

• dot → transpose(...)*ck: this is the headline bugfix. dot conjugates its first argument; the Fortran contraction is unconjugated. Correct, and likely the dominant contribution to the error collapse.
• _el_solution_at reproduces the ideal sing_der! sequence exactly — Q⁻¹F̄⁻¹(Q⁻¹u₂ − K̄u₁) with the same splined Cholesky-L convention (ldiv! L then Lᴴ, matching compute_axis_init), and the same vec([... for m, n]) m-fastest singfac ordering used by sing_der!/compute_axis_init. Using the very operators that generated the solution is the right reconstruction, and Hermite-interpolating u₁/u₂ with the stored du nodes is self-consistent.
• _solution_at (kinetic fallback): interpolating singfac·ξ′ and dividing the pole back out is a sound change of variables; row-wise broadcasting orientation (mode index = dim 1) checks out against the u_store[resnum, :, ...] layout.
• Store plumbing is complete: both crossing functions store du1 recomputed after the post-crossing solution surgery (consistent with u_store); resize_storage!/trim_storage! handle the new array; and both basis-transform sites (transform_u!, free_run!) apply identical solution-mixing transforms to du_store — right-multiplication commutes with d/dψ, so transformed du remains the derivative of transformed u.
• Path gating: use_du_store (zeros sentinel) correctly routes the Riccati/gal-matched odet — which never populates du_store — to the old chord-slope path, and kinetic_populated is set only under ctrl.kinetic_factor > 0 (make_kinetic_matrix is gated in main), so ideal runs always get the EL reconstruction. Replay (Rerun.jl) re-runs the pipeline rather than restoring OdeState from HDF5, so no snapshot-staleness issue.
• Save-time sing_der!: recomputing at the accepted point before every save is a quiet accuracy improvement to ud_store itself (previously it could hold the last internal RK-stage value). No feedback into the integrator trajectory — callback saves don't touch integrator state.

Comments (see line notes)

1. kinetic_populated inline comment misstates which matrices the kinetic path overwrites (minor, wording only).
2. du_store slot 1 now exactly duplicates ud_store slot 1 — the added memory (~size of u_store) could be halved; non-blocking, plus the zero-init sentinel deserves a protective comment.
3. The bare elseif odet.new branch lost its explanatory comment and now reads as dead code.
4. _solution_at stencil: same-side filter doesn't guard against a neighboring singular surface inside the window; 1-node degeneration is silent. Both unlikely, flagged for awareness.

Style / conventions

• New struct fields documented in the OdeState docstring (and the ud_store slot-2 clarification is a genuinely useful correction) — follows project conventions.
• No PR/issue references added to source; comments are concise; formatting (kwarg spacing, line width) conforms.

Tests / regression

No unit tests added; validation rests on the regression harness, and the PR description includes both fixed-grid and auto-grid tables per project policy — good. One gap worth closing before merge: the tables appear to exercise the ideal (EL) path only. The other two branches — _solution_at (kinetic runs) and the chord-slope fallback (gal-matched runs) — are new/changed routing that the shown tables don't cover. A kinetic-case regression run (or a note that one was done) would complete the picture.

Verdict

Physics-faithful, well-plumbed, and the regression evidence is strong. Only minor comment-accuracy and robustness notes; nothing blocking from my side.

⚠️ THIS PR MUST NOT BE MERGED WITHOUT THIRD-PARTY HUMAN REVIEW ⚠️

Per project policy (CLAUDE.md), no pull request is ever merged into develop without a human reviewer's approval. This review is advisory only and does not satisfy that requirement — it is non-negotiable. Human reviewers are already requested (@logan-nc, @jhalpern30); please wait for their approval.

---

Generated by Claude Code

[matt-pharr]

Minor: this comment states the reason slightly wrong. _compute_fkg_matrices! overwrites amats/bmats/cmats — fmats_lower and kmats are left untouched and remain the ideal EL operators. The actual reason the gate is needed is that the kinetic solution obeys the FKG kinetic ODE (F̄/K̄ built from f0mats/pmats/... with singfac absorbed, see the kinetic_factor > 0 branch of sing_der!), so the ideal EL relation in _el_solution_at does not reconstruct its Ξ′. The gating logic itself is correct — suggest rewording to something like: # set by make_kinetic_matrix; the solution then obeys the FKG ODE, not the ideal EL relation.

---

Generated by Claude Code

[matt-pharr]

Observation (non-blocking): slot 1 of du_store now duplicates slot 1 of ud_store at every saved step. sing_der! ends with odet.ud[:, :, 1] .= du1, and every save site now calls sing_der! immediately before storing both arrays from the same state, so du_store[:,:,1,:] ≡ ud_store[:,:,1,:] (both transform sites also apply identical transforms to both). The genuinely new information is only slot 2 (du₂/dψ vs. the Clebsch Ξ_s in ud_store).

This array is numpert_total² × 2 × nsteps ComplexF64 — for a large multi-n run that's an extra allocation on the order of the existing u_store (+50% of solution-storage memory), half of it redundant. If memory ever becomes a concern, du_store could hold only du₂/dψ (with _solution_at/_el_solution_at reading du₁/dψ from ud_store), halving the cost. Keeping them separate is defensible for clarity/robustness, so fine to leave as-is — but worth a note in the docstring that slot 1 mirrors ud_store slot 1 on the serial path.

Also: zeros(...) here (vs. undef for the siblings) is load-bearing — it's the sentinel that use_du_store = any(!iszero, ...) relies on, and resize_storage! correctly preserves that with zeros too. A one-line comment saying the zero-init is the "unpopulated" sentinel would protect it from a future "optimize to undef" cleanup.

---

Generated by Claude Code

[matt-pharr]

The old block kept a comment explaining why this call exists; the elseif branch now looks like dead code (its du_buffer result is discarded). Suggest restoring a one-liner, e.g. # Gaussian reduction modified u: recompute for the odet.ud side effect (Fortran ode_output.f sing_der before euler.bin write) — otherwise a future cleanup pass will delete it.

---

Generated by Claude Code

[matt-pharr]

Two minor robustness notes on the stencil selection:

1. The same-side filter only checks the current surface. The ±3 window can, in principle, reach across a neighboring singular surface (where the stored solution is also discontinuous from the crossing fixup) and include its far-side nodes in the Lagrange stencil. With dense near-boundary saving this needs two rational surfaces within ~4 saved steps of each other, so it's unlikely in practice — but a closely-spaced-surface case (high-n, low shear) could hit it silently.

2. If the window yields only 1 same-side node the "Lagrange cubic" degenerates to a constant with no warning. Guaranteed ≥1 by construction, but a length(idxs) >= 2 || @warn ... (maxlog=1) would make degradation visible.

Neither is blocking given the current save density; flagging for awareness.

---

Generated by Claude Code

[logan-nc]

@jhalpern30 this seems important but I can't get to a thorough review for over a week... if you could help give it a human read through, it would be much appreciated. Maybe include asking the AI why it was done how it was in the first place (based on reading past commit messages) too so we know the history well

[jhalpern30]

@matt-pharr before I dig into reviewing this, can you confirm:

1. is the regression table confirming that you fixed this? Is the main takeaway that ODE steps is changing between .17 and .15 but now the resonant metrics are similar?
2. are you done editing this branch? It looks like Claude made some comments, I didn't know if you planned on addressing them
3. does this exactly match the Fortran logic now, or are there any key differences?