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

[jhalpern30]

Idk anything about kinetic physics, you can judge if this is needed or not

[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

[jhalpern30]

I think this is a good suggestion - I think what I would propose is:

1. Only have one ud / du parameter. This should store exactly the derivative of u, such that we store both u and u' in eq. 24 of Glasser 2016 at every step. I vote either named du, u_der, or maybe u_prime. Something that makes it super clear what the variable is
2. Make a new storage parameter xi_s_store of size numpert_total² × nsteps, storing eq. 18 of Glasser 2016 on its own. I always though the u, ud, du naming was confusing, and wasn't sure why we lumped dXi/dPsi and Xi_s together like that

The new struct then looks like

u_store::Array{ComplexF64,4} = zeros(ComplexF64, numpert_total, numpert_total, 2, numsteps_init)
du_store::Array{ComplexF64,4} = zeros(ComplexF64, numpert_total, numpert_total, 2, numsteps_init)
xi_s_store::Array{ComplexF64,3} = zeros(ComplexF64, numpert_total, numpert_total, numsteps_init)

[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

[jhalpern30]

Agreed - I think is dead code. The way this is set up, sing_der! is called again whenever we're saving anyway. A sing_der! call on its own doesn't do anything, it just filled in the buffer which becomes unused. It only does anything if paired with 1) an integrator which uses it to take a step or 2) to recompute the derivatives when u changed so we can write them to memory

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

[matt-pharr]

@jhalpern30 I asked claude to review it to give you a head start lol. If you agree with any of its comments just say so and I can change stuff. The regression table is confirming the fix, yeah, and it is ready for your review. This solution is essentially equivalent to the fortran version now, though I reused @logan-nc's fancy bisecting/cubic splining code.

[jhalpern30]

This definitely looks like a good fix. From my quick catch-up on this code, it seemed like the Julia had a relatively hacky way of computing the field jump that didn't match the more sophisticated Fortran methods, and this PR adds in those two methods via el_solution_at and solution_at using new stored ODE derivatives. The jump methods now differ based on ideal, kinetic, or other branches.

Because it seems like the regression cases are now working, I'm trusting all of the physics in this PR. I also had Claude take a look at the Fortran and this branch and it didn't flag any essential missing functionality, so my review is mainly from the code cleanup perspective/making things clear for unfamiliar users.

I'm leaving this as a comment - @matt-pharr you can choose how much of my comments to integrate and can merge whenever you'd like. If you want me to take another look, just ping me again. I'd say the main essential ones are restructuring the redundant ud/du store and adding some helpful comments.

[jhalpern30]

I think this is a good suggestion - I think what I would propose is:

1. Only have one ud / du parameter. This should store exactly the derivative of u, such that we store both u and u' in eq. 24 of Glasser 2016 at every step. I vote either named du, u_der, or maybe u_prime. Something that makes it super clear what the variable is
2. Make a new storage parameter xi_s_store of size numpert_total² × nsteps, storing eq. 18 of Glasser 2016 on its own. I always though the u, ud, du naming was confusing, and wasn't sure why we lumped dXi/dPsi and Xi_s together like that

The new struct then looks like

u_store::Array{ComplexF64,4} = zeros(ComplexF64, numpert_total, numpert_total, 2, numsteps_init)
du_store::Array{ComplexF64,4} = zeros(ComplexF64, numpert_total, numpert_total, 2, numsteps_init)
xi_s_store::Array{ComplexF64,3} = zeros(ComplexF64, numpert_total, numpert_total, numsteps_init)

[jhalpern30]

Agreed - I think is dead code. The way this is set up, sing_der! is called again whenever we're saving anyway. A sing_der! call on its own doesn't do anything, it just filled in the buffer which becomes unused. It only does anything if paired with 1) an integrator which uses it to take a step or 2) to recompute the derivatives when u changed so we can write them to memory

[jhalpern30]

This feels super inefficient - more or less every integrator step (if save_interval = 1), we compute the derivative twice for no reason, not just when we renormalize the solutions.

If you redo OdeState based on my comment below, you could just make it such that it stores all the necessary information at the end of sing_der! in odet.u and new variables odet.du and odet.xi_s. You could then keep the same logic of only recomputing if u has changed.

I think there are even cleaner ways to do this but I think this would be relatively simple but helpful

[jhalpern30]

If you're feeling adventurous, I'm not realizing this could just be a Utility function since we call it 4 separate times

function store_ode_data!(odet, step, psi, u)
  odet.psi_store[odet.step] = psi
  @views odet.u_store[:, :, :, odet.step] .= u
  odet.q_store[odet.step] = odet.q
  @views odet.ud_store[:, :, :, odet.step] .= odet.ud
  @views odet.xi_s_store[:, :, :, odet.step] .= odet.xi_s
  odet.step += 1
end

[jhalpern30]

Idk anything about kinetic physics, you can judge if this is needed or not

[jhalpern30]

I think this logic is kinda confusing. My quick understanding from this is that:

• use_el is the ideal branch for when we have both u_store and can interpolate it and its derivative using the ideal ODE
• use_du_store is the kinetic path, where you also interpolate the derivative
• the non use_du_store path is for resistive/Ricatti runs where du_store does not exist so we must use the old method

If this is correct, I have a few comments

1. I don't like conditioning it on iszero - this seems finicky and requires a different initialization of the other store variables (Undef vs. zeroed). Can the logic instead be made to be based on the input flags? It seems like kin_flag and gal_flag should be able to do this?
2. Add some one-line comments/info prints to terminal in the conditionals around line 310 below to make it clear what method is being used:

# Evaluate u and dξ/dψ at lpsi and rpsi from the stored ODE solution
if !use_du_store
  # Interpolate u and finite-difference for dxi/dpsi
...
elseif use_el
  # Interpolate u and compute du from EL equation
etc.

[jhalpern30]

In the long term, I think 3 different methods here seems overkill, but from what I can find this is a faithful reproduction of the current state of the Fortran so its a good starting point, just wanted to note

[jhalpern30]

I also feel like a large amount of this file is unpacking indices and its ugly. Again, a future fix tho

[jhalpern30]

Bad Claude comment - it always tries to reference code it got rid of. Its ok to just fix the bug here and not explain why code that is no longer there was removed

[jhalpern30]

Can you pass resnum to this function too to make this more readable? Such that _solution_at and _el_solution_at just return u_l, u_r, ud_l, and ud_r directly? This section of code is already super difficult to read as a human because of the amount of lines dedicated to naming left and right variables