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