PERTURBED EQUILIBRIUM - BUG FIX - Use Λ† in reluctance to match Fortran gpresp_reluct

[logan-nc]

Summary

The flux-space reluctance matrix in src/PerturbedEquilibrium/Response.jl was computed with the plasma inductance Λ where the Fortran GPEC reference (gpresp_reluct) uses its adjoint Λ†:

	formula
Julia (before)	ϱ = L⁻¹·(Λ − L)·L⁻¹
Fortran gpresp_reluct	ϱ = L⁻¹·(Λ† − L)·L⁻¹
Julia (after)	ϱ = L⁻¹·(Λ† − L)·L⁻¹ ✅

Λ (plasma inductance) is not Hermitian — its anti-Hermitian part carries the dissipative/torque response — so dropping the adjoint introduces a small but real error in the reported reluctance. For a perfectly Hermitian Λ the two forms coincide, which is why it went unnoticed.

This was surfaced by a fortran-physics-reviewer audit of the coordinate-invariant field-space work (merged in #277). It is a pre-existing inconsistency — it predates #277 and was only flagged during that review — hence this separate bugfix branched from develop.

The fix

# Reluctance ϱ = L⁻¹·(Λ† − L)·L⁻¹ (Fortran gpresp_reluct: diff_indmats = CONJG(TRANSPOSE(plas_indmats)) − surf_indmats).
# Λ (plasma inductance) is not Hermitian — its anti-Hermitian part is the dissipative/torque response — so the adjoint matters.
L_inv = inv(surface_inductance)
reluctance = L_inv * (plasma_inductance' - surface_inductance) * L_inv

(' is the Julia adjoint = conjugate transpose, i.e. CONJG(TRANSPOSE(...)).)

Scope / blast radius

reluctance is an output-only quantity: it is stored in PerturbedEquilibriumState, conformed to field space by field_space_response_matrices, and written to HDF5 — nothing downstream consumes it (permeability is built directly from Λ and L, not from ϱ). The fix therefore only changes the reluctance dataset itself.

Verification

• using GeneralizedPerturbedEquilibrium loads cleanly.
• test/runtests_coordinate_invariant.jl: 15/15 + energy net-zero still green (those tests exercise the congruence conform on synthetic inputs, not this production formula).
• Regression harness diiid_n1, develop vs this branch: 46 quantities unchanged, 0 moved — confirms the change is isolated (no pinned quantity moves; reluctance is not pinned).

Fortran reference (for convenience — exact lines, no need to pull up the source)

~/Code/gpec/gpec/gpresp.f, SUBROUTINE gpresp_reluct (lines 423–461). temp1 is set to the identity then zhetrf/zhetrs solve surf_indmats·X = I, so temp1 = L⁻¹:

         temp1=0
         DO i=1,mpert
            temp1(i,i)=1
         ENDDO
         temp2=surf_indmats
         CALL zhetrf('L',mpert,temp2,mpert,ipiv,work2,mpert*mpert,info)
         CALL zhetrs('L',mpert,mpert,temp2,mpert,ipiv,temp1,mpert,info)
         diff_indmats(j,:,:)=CONJG(TRANSPOSE(plas_indmats(j,:,:)))
     $        -surf_indmats
         reluctmats(j,:,:)=MATMUL(temp1,
     $        MATMUL(diff_indmats(j,:,:),temp1))

Key lines:

• L449–450: diff_indmats = CONJG(TRANSPOSE(plas_indmats)) − surf_indmats → Λ† − L
• L451–452: reluctmats = MATMUL(temp1, MATMUL(diff_indmats, temp1)) → L⁻¹·(Λ† − L)·L⁻¹

🤖 Generated with Claude Code