Skip to content

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

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
matt-pharr wants to merge 1 commit into
base: develop
Choose a base branch
from
Open

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

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
18 changes: 12 additions & 6 deletions src/ForceFreeStates/EulerLagrange.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -645,6 +645,7 @@ function cross_ideal_singular_surf!(
odet.q_store[odet.step] = odet.q
odet.u_store[:, :, :, odet.step] = odet.u
odet.ud_store[:, :, :, odet.step] = odet.ud
odet.du_store[:, :, :, odet.step] = du1
odet.step += 1
end

Expand Down Expand Up @@ -696,6 +697,7 @@ function cross_kinetic_singular_surf!(
odet.q_store[odet.step] = odet.q
odet.u_store[:, :, :, odet.step] = odet.u
odet.ud_store[:, :, :, odet.step] = odet.ud
odet.du_store[:, :, :, odet.step] = du1
odet.step += 1
end

Expand Down Expand Up @@ -757,12 +759,6 @@ function integrate_el_region!(

compute_solution_norms!(integrator.u, odet, ctrl, intr, false)

# If Gaussian reduction modified u, recompute ud to keep it consistent.
# Matches Fortran ode_output.f which calls sing_der before writing euler.bin.
if odet.new
sing_der!(du_buffer, integrator.u, integrator.p, integrator.t)
end

# Save near segment boundaries (symmetric, in q not psi) and every Nth step.
# The step-count fallback (== 1) guarantees the first step is always saved
# even for near-degenerate segments where q_range ≈ 0.
Expand All @@ -775,11 +771,15 @@ function integrate_el_region!(
if odet.step >= size(odet.u_store, 4)
resize_storage!(odet)
end
sing_der!(du_buffer, integrator.u, integrator.p, integrator.t)

@jhalpern30 jhalpern30 Aug 3, 2026

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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 jhalpern30 Aug 3, 2026

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

odet.psi_store[odet.step] = integrator.t
@views odet.u_store[:, :, :, odet.step] .= integrator.u
odet.q_store[odet.step] = odet.q
@views odet.ud_store[:, :, :, odet.step] .= odet.ud
@views odet.du_store[:, :, :, odet.step] .= du_buffer
odet.step += 1
elseif odet.new
sing_der!(du_buffer, integrator.u, integrator.p, integrator.t)

@matt-pharr matt-pharr Jul 31, 2026

Copy link
Copy Markdown
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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 jhalpern30 Aug 3, 2026

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

end
end

Expand All @@ -794,10 +794,12 @@ function integrate_el_region!(
if odet.step >= size(odet.u_store, 4)
resize_storage!(odet)
end
sing_der!(du_buffer, sol.u[end], (ctrl, equil, ffit, intr, odet, chunk), sol.t[end])
odet.psi_store[odet.step] = sol.t[end]
@views odet.u_store[:, :, :, odet.step] .= sol.u[end]
odet.q_store[odet.step] = odet.q
@views odet.ud_store[:, :, :, odet.step] .= odet.ud
@views odet.du_store[:, :, :, odet.step] .= du_buffer
odet.step += 1
end

Expand Down Expand Up @@ -1035,6 +1037,10 @@ function transform_u!(odet::OdeState, intr::ForceFreeStatesInternal)
odet.ud_store[:, :, 1, istep] .= gauss_buffer
mul!(gauss_buffer, odet.ud_store[:, :, 2, istep], transforms[:, :, ifix])
odet.ud_store[:, :, 2, istep] .= gauss_buffer
mul!(gauss_buffer, odet.du_store[:, :, 1, istep], transforms[:, :, ifix])
odet.du_store[:, :, 1, istep] .= gauss_buffer
mul!(gauss_buffer, odet.du_store[:, :, 2, istep], transforms[:, :, ifix])
odet.du_store[:, :, 2, istep] .= gauss_buffer
end
jfix = kfix + 1
end
Expand Down
7 changes: 6 additions & 1 deletion src/ForceFreeStates/ForceFreeStatesStructs.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -346,6 +346,8 @@ end
# Kinetic torque matrix splines: 6 components
ktmats::Vector{S} = [_empty_series_interp_complex(numpert_total^2, itp_opts) for _ in 1:6]

kinetic_populated::Bool = false # set by make_kinetic_matrix; fmats_lower/kmats are then not the EL operators

@matt-pharr matt-pharr Jul 31, 2026

Copy link
Copy Markdown
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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 jhalpern30 Aug 3, 2026

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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


# Pre-computed FKG kinetic matrices (populated by make_kinetic_matrix)
f0mats::S = _empty_series_interp_complex(numpert_total^2, itp_opts)
pmats::S = _empty_series_interp_complex(numpert_total^2, itp_opts)
Expand Down Expand Up @@ -492,7 +494,9 @@ and a small set of temporary matrices and factors used to compute singular-layer
- `u_store::Array{ComplexF64,4}` - Stored solution arrays at each saved step with shape
`(numpert_total, numpert_total, 2, numsteps_init)` (complex solution state used by the solver).

- `ud_store::Array{ComplexF64,4}` - Stored derivatives of the solution at each saved step with same shape as `u_store`.
- `ud_store::Array{ComplexF64,4}` - Stored derivatives of the solution at each saved step with same shape as `u_store`; slot 2 holds the Clebsch Ξ_s, not du₂/dψ.

- `du_store::Array{ComplexF64,4}` - du₁/dψ and du₂/dψ at the accepted point of each saved step, populated only by the serial EL path.

- `crit_store::Vector{Float64}` - Stored crit parameter values (smallest eigenvalue of W⁻ꜝ) (length `numsteps_init`).

Expand Down Expand Up @@ -560,6 +564,7 @@ and a small set of temporary matrices and factors used to compute singular-layer
q_store::Vector{Float64} = Vector{Float64}(undef, numsteps_init)
u_store::Array{ComplexF64,4} = Array{ComplexF64}(undef, numpert_total, numpert_total, 2, numsteps_init)
ud_store::Array{ComplexF64,4} = Array{ComplexF64}(undef, numpert_total, numpert_total, 2, numsteps_init)
du_store::Array{ComplexF64,4} = zeros(ComplexF64, numpert_total, numpert_total, 2, numsteps_init)

@matt-pharr matt-pharr Jul 31, 2026

Copy link
Copy Markdown
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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 jhalpern30 Aug 3, 2026

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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)

crit_store::Vector{Float64} = Vector{Float64}(undef, numsteps_init)
ca_r::Array{ComplexF64,4} = Array{ComplexF64}(undef, numpert_total, numpert_total, 2, msing)
ca_l::Array{ComplexF64,4} = Array{ComplexF64}(undef, numpert_total, numpert_total, 2, msing)
Expand Down
4 changes: 4 additions & 0 deletions src/ForceFreeStates/Free.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -130,6 +130,10 @@ and data dumping.
odet.ud_store[:, :, 1, istep] .= tmp_mat
mul!(tmp_mat, odet.ud_store[:, :, 2, istep], coeffs)
odet.ud_store[:, :, 2, istep] .= tmp_mat
mul!(tmp_mat, odet.du_store[:, :, 1, istep], coeffs)
odet.du_store[:, :, 1, istep] .= tmp_mat
mul!(tmp_mat, odet.du_store[:, :, 2, istep], coeffs)
odet.du_store[:, :, 2, istep] .= tmp_mat
end

# Write energies to screen
Expand Down
1 change: 1 addition & 0 deletions src/ForceFreeStates/Kinetic.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -257,6 +257,7 @@ function _compute_fkg_matrices!(
ffit.amats = cubic_interp(xs, Series(ak_flat); ffit.itp_opts...)
ffit.bmats = cubic_interp(xs, Series(bk_flat); ffit.itp_opts...)
ffit.cmats = cubic_interp(xs, Series(ck_flat); ffit.itp_opts...)
ffit.kinetic_populated = true

return nothing
end
4 changes: 4 additions & 0 deletions src/ForceFreeStates/Utils.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -12,18 +12,21 @@ function resize_storage!(odet::OdeState)
# Allocate new arrays
u_new = Array{ComplexF64,4}(undef, odet.numpert_total, odet.numpert_total, 2, newlen)
ud_new = Array{ComplexF64,4}(undef, odet.numpert_total, odet.numpert_total, 2, newlen)
du_new = zeros(ComplexF64, odet.numpert_total, odet.numpert_total, 2, newlen)
psi_new = Vector{Float64}(undef, newlen)
q_new = Vector{Float64}(undef, newlen)

# Copy old data
u_new[:, :, :, 1:odet.step] = odet.u_store[:, :, :, 1:odet.step]
ud_new[:, :, :, 1:odet.step] = odet.ud_store[:, :, :, 1:odet.step]
du_new[:, :, :, 1:odet.step] = odet.du_store[:, :, :, 1:odet.step]
psi_new[1:odet.step] = odet.psi_store[1:odet.step]
q_new[1:odet.step] = odet.q_store[1:odet.step]

# Replace old arrays
odet.u_store = u_new
odet.ud_store = ud_new
odet.du_store = du_new
odet.psi_store = psi_new
odet.q_store = q_new
end
Expand All @@ -40,4 +43,5 @@ function trim_storage!(odet::OdeState)
resize!(odet.q_store, odet.step)
odet.u_store = odet.u_store[:, :, :, 1:odet.step]
odet.ud_store = odet.ud_store[:, :, :, 1:odet.step]
odet.du_store = odet.du_store[:, :, :, 1:odet.step]
end
2 changes: 1 addition & 1 deletion src/PerturbedEquilibrium/PerturbedEquilibrium.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -127,7 +127,7 @@ function compute_perturbed_equilibrium(
if vac_data === nothing
@warn "Vacuum data not available. Skipping singular coupling calculation. Set vac_flag=true in [ForceFreeStates] section."
else
compute_singular_coupling_metrics!(state, equil, ForceFreeStates_results, vac_data, ffs_intr, intr, ctrl)
compute_singular_coupling_metrics!(state, equil, ForceFreeStates_results, vac_data, ffs_intr, intr, ctrl, ffit)
end
end

Expand Down
191 changes: 154 additions & 37 deletions src/PerturbedEquilibrium/SingularCoupling.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -43,6 +43,106 @@ end
# Reflect a periodic theta-space vector θ → -θ (the theta reversal in gpvacuum_flxsurf).
_reverse_theta(v::AbstractVector) = circshift(reverse(v), 1)

"""
_solution_at(psi, psi_surf, odet, equil, ffs_intr, nstep) -> (u1, du1)

Evaluate Ξ_ψ and Ξ′_ψ at `psi` from the stored ODE solution: cubic Hermite for the value,
Lagrange cubic through the singfac-weighted `du_store` nodes for the derivative. Works for
any formulation; ideal runs use the more accurate `_el_solution_at`. The Ξ′ stencil stays
on `psi`'s side of the singular surface `psi_surf`, where the stored solution is
discontinuous. Requires a populated `odet.du_store`.
"""
function _solution_at(
psi::Float64,
psi_surf::Float64,
odet::OdeState,
equil::Equilibrium.PlasmaEquilibrium,
ffs_intr::ForceFreeStatesInternal,
nstep::Int
)
il, ir, _ = _psi_bracket(odet.psi_store, psi, nstep)
psi_a, psi_b = odet.psi_store[il], odet.psi_store[ir]

u1_a = @view odet.u_store[:, :, 1, il]
u1_b = @view odet.u_store[:, :, 1, ir]
du1_a = @view odet.du_store[:, :, 1, il]
du1_b = @view odet.du_store[:, :, 1, ir]

u1_e = _hermite_cubic_val(u1_a, u1_b, du1_a, du1_b, psi_a, psi_b, psi)

# Same-side candidate nodes around the bracket, trimmed to the 4 nearest psi.
side = sign(psi - psi_surf)
idxs = [j for j in max(1, il - 3):min(nstep, ir + 3) if sign(odet.psi_store[j] - psi_surf) == side]
while length(idxs) > 4
abs(odet.psi_store[idxs[1]] - psi) > abs(odet.psi_store[idxs[end]] - psi) ? popfirst!(idxs) : pop!(idxs)
end
Comment on lines +73 to +78

@matt-pharr matt-pharr Jul 31, 2026

Copy link
Copy Markdown
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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


# Interpolate singfac·Ξ′ and divide the ~1/singfac pole back out at psi.
singfac_at(p) = vec([m - equil.profiles.q_spline(p) * n for m in ffs_intr.mlow:ffs_intr.mhigh, n in ffs_intr.nlow:ffs_intr.nhigh])
du1_e = zeros(ComplexF64, size(odet.du_store, 1), size(odet.du_store, 2))
for k in eachindex(idxs)
w = 1.0
for j in eachindex(idxs)
j == k && continue
w *= (psi - odet.psi_store[idxs[j]]) / (odet.psi_store[idxs[k]] - odet.psi_store[idxs[j]])
end
du1_e .+= (w .* singfac_at(odet.psi_store[idxs[k]])) .* @view(odet.du_store[:, :, 1, idxs[k]])
end
du1_e ./= singfac_at(psi)

return u1_e, du1_e
end

"""
_el_solution_at(psi, odet, ffit, equil, ffs_intr, nstep) -> (u1, du1)

Evaluate Ξ_ψ and Ξ′_ψ at `psi` from the stored ODE solution via the ideal Euler-Lagrange
relation Ξ′ = Q⁻¹·F̄⁻¹·(Q⁻¹·u₂ − K̄·u₁) [Glasser 2016 eqs. 22-24], with u₁, u₂ Hermite-
interpolated to `psi`. Only valid for ideal runs where `ffit.fmats_lower` and `kmats`
generated the solution. Requires a populated `odet.du_store`.
"""
function _el_solution_at(
psi::Float64,
odet::OdeState,
ffit::FourFitVars,
equil::Equilibrium.PlasmaEquilibrium,
ffs_intr::ForceFreeStatesInternal,
nstep::Int
)
npert = ffs_intr.numpert_total
il, ir, _ = _psi_bracket(odet.psi_store, psi, nstep)
psi_a, psi_b = odet.psi_store[il], odet.psi_store[ir]

u1_a = @view odet.u_store[:, :, 1, il]
u1_b = @view odet.u_store[:, :, 1, ir]
u2_a = @view odet.u_store[:, :, 2, il]
u2_b = @view odet.u_store[:, :, 2, ir]
du1_a = @view odet.du_store[:, :, 1, il]
du1_b = @view odet.du_store[:, :, 1, ir]
du2_a = @view odet.du_store[:, :, 2, il]
du2_b = @view odet.du_store[:, :, 2, ir]

hint = Ref(1)
kmat = Matrix{ComplexF64}(undef, npert, npert)

u1_e = _hermite_cubic_val(u1_a, u1_b, du1_a, du1_b, psi_a, psi_b, psi)
u2_e = _hermite_cubic_val(u2_a, u2_b, du2_a, du2_b, psi_a, psi_b, psi)

# Ξ′ = Q⁻¹·F̄⁻¹·(Q⁻¹·u₂ − K̄·u₁) with Q⁻¹ = diag(1/(m − n·q))
q_e = equil.profiles.q_spline(psi)
singfac_inv = vec([1.0 / (m - q_e * n) for m in ffs_intr.mlow:ffs_intr.mhigh, n in ffs_intr.nlow:ffs_intr.nhigh])
fmat_lower = Matrix{ComplexF64}(undef, npert, npert)
ffit.fmats_lower(vec(fmat_lower), psi; hint=hint)
ffit.kmats(vec(kmat), psi; hint=hint)
du1_e = u2_e .* singfac_inv
du1_e .-= kmat * u1_e
ldiv!(LowerTriangular(fmat_lower), du1_e)
ldiv!(UpperTriangular(fmat_lower'), du1_e)
du1_e .*= singfac_inv

return u1_e, du1_e
end

"""
compute_singular_coupling_metrics!(
state::PerturbedEquilibriumState,
Expand All @@ -51,7 +151,8 @@ _reverse_theta(v::AbstractVector) = circshift(reverse(v), 1)
vac_data::VacuumData,
ffs_intr::ForceFreeStatesInternal,
intr::PerturbedEquilibriumInternal,
ctrl::PerturbedEquilibriumControl
ctrl::PerturbedEquilibriumControl,
ffit::FourFitVars
)

Compute singular layer coupling matrices and applied resonant vectors.
Expand Down Expand Up @@ -82,7 +183,8 @@ function compute_singular_coupling_metrics!(
vac_data::VacuumData,
ffs_intr::ForceFreeStatesInternal,
intr::PerturbedEquilibriumInternal,
ctrl::PerturbedEquilibriumControl
ctrl::PerturbedEquilibriumControl,
ffit::FourFitVars
)
ctrl.verbose && @info "Computing singular coupling metrics (GPEC method)"

Expand Down Expand Up @@ -159,6 +261,10 @@ function compute_singular_coupling_metrics!(
# @threads region.
psi_store_all = ForceFreeStates_results.psi_store
nstep = ForceFreeStates_results.step
# ξ′ evaluation preference: ideal EL relation, then interpolated stored RHS for kinetic,
# then chord slope for paths that never populate du_store like the gal-matched odet.
use_du_store = any(!iszero, ForceFreeStates_results.du_store)

@jhalpern30 jhalpern30 Aug 3, 2026

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

use_el = use_du_store && !ffit.kinetic_populated
_blas_nthreads = BLAS.get_num_threads()
BLAS.set_num_threads(1)
try
Expand All @@ -169,7 +275,7 @@ function compute_singular_coupling_metrics!(

resnum = findfirst(j -> intr.m_modes[j] == m_res && intr.n_modes[j] == nn, 1:numpert_total)
if resnum === nothing
@warn "Could not find index for resonant mode (m=$m_res, n=$nn)" maxlog=1
@warn "Could not find index for resonant mode (m=$m_res, n=$nn)" maxlog = 1
continue
end

Expand Down Expand Up @@ -201,36 +307,46 @@ function compute_singular_coupling_metrics!(
lpsi = sing_surf.psifac - spot_psi
rpsi = sing_surf.psifac + spot_psi

# Interpolate u and du/dψ at lpsi and rpsi using stored ODE solution.
# Value (u): Hermite cubic using ud_store for smoother interpolation than linear.
# Derivative (ud): chord slope from u_store only — see comment below.
il_l, ir_l, _ = _psi_bracket(psi_store_all, lpsi, nstep)
il_r, ir_r, _ = _psi_bracket(psi_store_all, rpsi, nstep)

u_node = ForceFreeStates_results.u_store
ud_node = ForceFreeStates_results.ud_store
ua_l = u_node[resnum, :, 1, il_l]
ub_l = u_node[resnum, :, 1, ir_l]
ua_r = u_node[resnum, :, 1, il_r]
ub_r = u_node[resnum, :, 1, ir_r]
dua_l = ud_node[resnum, :, 1, il_l]
dub_l = ud_node[resnum, :, 1, ir_l]
dua_r = ud_node[resnum, :, 1, il_r]
dub_r = ud_node[resnum, :, 1, ir_r]

psi_il_l = psi_store_all[il_l]
psi_ir_l = psi_store_all[ir_l]
psi_il_r = psi_store_all[il_r]
psi_ir_r = psi_store_all[ir_r]

u_l = _hermite_cubic_val(ua_l, ub_l, dua_l, dub_l, psi_il_l, psi_ir_l, lpsi)
u_r = _hermite_cubic_val(ua_r, ub_r, dua_r, dub_r, psi_il_r, psi_ir_r, rpsi)
# Derivative (ud): chord slope from u_store only — ud_store can be systematically off near
# outer surfaces (q=4/5) where the ODE solution varies rapidly; near-cancellation in bwp1
# then amplifies even a ~10% ud_store error into a large Delta' error. Chord slope avoids
# this by using only u values, which are accurately stored by the ODE integrator.
ud_l = (ub_l .- ua_l) ./ (psi_ir_l - psi_il_l)
ud_r = (ub_r .- ua_r) ./ (psi_ir_r - psi_il_r)
# Evaluate u and dξ/dψ at lpsi and rpsi from the stored ODE solution
if !use_du_store

@jhalpern30 jhalpern30 Aug 3, 2026

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

il_l, ir_l, _ = _psi_bracket(psi_store_all, lpsi, nstep)
il_r, ir_r, _ = _psi_bracket(psi_store_all, rpsi, nstep)

u_node = ForceFreeStates_results.u_store
ud_node = ForceFreeStates_results.ud_store
ua_l = u_node[resnum, :, 1, il_l]

@jhalpern30 jhalpern30 Aug 3, 2026

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

ub_l = u_node[resnum, :, 1, ir_l]
ua_r = u_node[resnum, :, 1, il_r]
ub_r = u_node[resnum, :, 1, ir_r]
dua_l = ud_node[resnum, :, 1, il_l]
dub_l = ud_node[resnum, :, 1, ir_l]
dua_r = ud_node[resnum, :, 1, il_r]
dub_r = ud_node[resnum, :, 1, ir_r]

psi_il_l = psi_store_all[il_l]
psi_ir_l = psi_store_all[ir_l]
psi_il_r = psi_store_all[il_r]
psi_ir_r = psi_store_all[ir_r]

u_l = _hermite_cubic_val(ua_l, ub_l, dua_l, dub_l, psi_il_l, psi_ir_l, lpsi)
u_r = _hermite_cubic_val(ua_r, ub_r, dua_r, dub_r, psi_il_r, psi_ir_r, rpsi)
ud_l = (ub_l .- ua_l) ./ (psi_ir_l - psi_il_l)
ud_r = (ub_r .- ua_r) ./ (psi_ir_r - psi_il_r)
elseif use_el
u1m_l, du1m_l = _el_solution_at(lpsi, ForceFreeStates_results, ffit, equil, ffs_intr, nstep)
u1m_r, du1m_r = _el_solution_at(rpsi, ForceFreeStates_results, ffit, equil, ffs_intr, nstep)
u_l = u1m_l[resnum, :]

@jhalpern30 jhalpern30 Aug 3, 2026

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

u_r = u1m_r[resnum, :]
ud_l = du1m_l[resnum, :]
ud_r = du1m_r[resnum, :]
else
u1m_l, du1m_l = _solution_at(lpsi, sing_surf.psifac, ForceFreeStates_results, equil, ffs_intr, nstep)
u1m_r, du1m_r = _solution_at(rpsi, sing_surf.psifac, ForceFreeStates_results, equil, ffs_intr, nstep)
u_l = u1m_l[resnum, :]
u_r = u1m_r[resnum, :]
ud_l = du1m_l[resnum, :]
ud_r = du1m_r[resnum, :]
end

q_l = equil.profiles.q_spline(lpsi)
q1_l = equil.profiles.q_deriv(lpsi)
Expand All @@ -242,10 +358,11 @@ function compute_singular_coupling_metrics!(
jump_vec = Vector{ComplexF64}(undef, numpert_total)
for k in 1:numpert_total
ck = @view C_coeffs[:, k]
xsp_l = dot(u_l, ck)
xsp1_l = dot(ud_l, ck)
xsp_r = dot(u_r, ck)
xsp1_r = dot(ud_r, ck)
# unconjugated contraction, dot would conjugate u

@jhalpern30 jhalpern30 Aug 3, 2026

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

xsp_l = transpose(u_l) * ck
xsp1_l = transpose(ud_l) * ck
xsp_r = transpose(u_r) * ck
xsp1_r = transpose(ud_r) * ck
bwp1_l = 2π * im * chi1 * (singfac_l * xsp1_l - nn * q1_l * xsp_l)
bwp1_r = 2π * im * chi1 * (singfac_r * xsp1_r - nn * q1_r * xsp_r)
jump_vec[k] = bwp1_r - bwp1_l
Expand Down
Loading