VACUUM - Migrate surface current calculation into the VACUUM module

[jhalpern30]

This PR moves the surface inductance calculation (or at least the majority of it) into the VACUUM module. The surface inductance calculation is fundamentally a vacuum quantity. It was hacked into GPEC to avoid cross-code changes, but can more simply be integrated into the VACUUM module and passed into the PerturbedEquilibrium module.

This also makes it a boolean if it is computed or not. For the standard stability calculation where only Wv is needed, computing the surface current calculation costs an extra matrix inversion + storage since it requires forming and solving the interior kernel, which can add up especially when making the Wv spline in the edge. I didn't benchmark the speedup in 2D since its not massive, but this will be very helpful for 3D.

I have looked at the Fortran and checked that no not-yet-ported code will require the old way of doing this. And if anything new does, its more suitable to be coded up in the VACUUM module instead anyway.

Some notes on specific changes:

• the grri/grre matrices are now internal to the VACUUM module only. The only quantity the rest of the code sees is the surface current matrix I_v
• The machinery for compute_surface_current_from_greens is now drastically simplified
• I removed the regularization logic surrounding the formation of the surface inductance matrix in SingularCoupling.jl. I looked at the Fortran and this was unfounded, and screamed "Claude did this because it had one failing test case while porting over the code and added try/catch logic". This same type of code appeared when we were doing the edge matrix splining for psiedge < 1 logic.
• I've added more documentation on where code that handles the difference in theta orientation between GPEC and VACUUM appears and what it does
• I've moved the Hermitization of the surface inductance into the VACUUM module

---

Update (2026-08-14): Addressed review feedback and merged latest develop:

• calc_surface_inductance now takes ψ and does the full VacuumInput → compute_vacuum_response(compute_Iv=true) → inversion chain internally; both call sites collapse to one line (per review)
• Removed the one non-pooled allocation in the compute_Iv branch (in-place difference into grri scratch)
• Corrected the interior-operator sign-convention comments (the operator is negated relative to Fortran's kernelsign=-1 convention; the θ handedness is handled entirely by conj!), restored the Hermitization physics rationale, and documented I_v's stored normalization and the Hermitize-before-inversion deviation from Fortran
• Added a unit test for calc_surface_inductance (Solovev, nowall); fixed a stale test comment
• Merged origin/develop (tearing tests now build SingType without the removed grri/grre fields; kept develop's new _hermite_cubic_deriv, dropped the now-unused _reverse_theta)

[jhalpern30]

Regression tests - this only really affects PerturbedEquilibrium so I ran DIIID

Regression Report: diiid_n1
================================================================================================================
Ref 1: local @ local (2026-08-03)
Ref 2: develop @ 167ada76 (2026-07-30)
----------------------------------------------------------------------------------------------------------------
Quantity local develop Diff Status ----------------------------------------------------------------------------------------------------------------
total energy Re(et[1]) 7.923042e-01 7.923042e-01 0.0e+00 OK total energy Im(et[1]) 1.261466e-04 1.261466e-04 0.0e+00 OK plasma energy Re(ep[1]) -1.375573e+00 -1.375573e+00 0.0e+00 OK vacuum energy Re(ev[1]) 2.167877e+00 2.167877e+00 0.0e+00 OK vacuum matrix min eigenvalue 1.873975e-01 1.873975e-01 0.0e+00 OK plasma energy (all) [35 elem] [35 elem] 0.0e+00 OK vacuum energy (all) [35 elem] [35 elem] 0.0e+00 OK total energy (all) [35 elem] [35 elem] 0.0e+00 OK ODE steps (saved) 1353 1353 0.0e+00 OK ODE steps (total) 1990 1990 0.0e+00 OK q0 1.204202e+00 1.204202e+00 0.0e+00 OK q95 4.781724e+00 4.781724e+00 0.0e+00 OK beta_t 1.327082e-02 1.327082e-02 0.0e+00 OK beta_n 1.372570e+00 1.372570e+00 0.0e+00 OK internal inductance li1 8.842261e-01 8.842261e-01 0.0e+00 OK internal inductance li2 7.080767e-01 7.080767e-01 0.0e+00 OK internal inductance li3 7.304351e-01 7.304351e-01 0.0e+00 OK poloidal beta betap1 6.681028e-01 6.681028e-01 0.0e+00 OK poloidal beta betap2 5.350080e-01 5.350080e-01 0.0e+00 OK poloidal beta betap3 5.519015e-01 5.519015e-01 0.0e+00 OK # singular surfaces 5 5 0.0e+00 OK singular psi locations [5 elem] [5 elem] 0.0e+00 OK singular q values [5 elem] [5 elem] 0.0e+00 OK current beta betaj 4.236772e-01 4.236772e-01 0.0e+00 OK plasma volume 1.829476e+01 1.829476e+01 0.0e+00 OK plasma current 1.152130e+00 1.152130e+00 0.0e+00 OK mpert 35 35 0.0e+00 OK npert 1 1 0.0e+00 OK toroidal field bt0 2.006573e+00 2.006573e+00 0.0e+00 OK wall field bwall 3.880145e-01 3.880145e-01 0.0e+00 OK aspect ratio 2.845746e+00 2.845746e+00 0.0e+00 OK elongation kappa 1.708322e+00 1.708322e+00 0.0e+00 OK q profile (checksum) 6e5ed3ba20f0... 6e5ed3ba20f0... identical OK pressure profile (checksum) ba7ac449bf30... ba7ac449bf30... identical OK Mercier D_I profile (checksum) 414eb69dd7de... 414eb69dd7de... identical OK resistive interchange D_R profile (checksum) 1978c3923c51... 1978c3923c51... identical OK ballooning Delta' profile (checksum) 908bad88b3c5... 908bad88b3c5... identical OK delta prime (BVP diagonal) [5 elem] [5 elem] 0.0e+00 OK island half-widths [5 elem] [5 elem] 1.5e-13 OK Chirikov parameter [5 elem] [5 elem] 1.3e-12 OK ||resonant area-weighted field|| 3.953877e-04 3.953877e-04 6.4e-16 OK PE plasma energy 3.444911e+00 3.444911e+00 3.3e-11 OK PE vacuum energy 3.174509e+00 3.174509e+00 1.3e-11 OK PE surface energy 5.962637e+00 5.962637e+00 1.196e-10 (0.00%) ** CHANGED **
PE toroidal torque -5.031207e-02 -5.031207e-02 2.9e-13 OK NTV torque FGAR [N·m] 5.772805e-01 5.772805e-01 1.3e-11 OK NTV kinetic energy dW FGAR [J] 7.235890e-02 7.235890e-02 1.6e-12 OK Runtime (s) 180.2s 170.1s -- resonant area-weighted field b^r [5 elem] [5 elem] 2.1e-15 OK ================================================================================================================
Summary: 1 changed, 47 unchanged

The small changes are due to slight variations in the order that Hermitization is performed. The 1 changed value is not physical, since it is just slightly above the arbitrary limit of 1e-10 that we set

[jhalpern30]

@claude review

[claude]

Review: #345 — Migrate surface current calculation into VACUUM module

Verified the core refactor algebraically: the interior-kernel solve now uses A - 2I (subtracting 2 from the un-negated grad_green diagonal) instead of the old -A + 2I, and the downstream combination correspondingly switched from grri + grre to grre - grri. Working through both, grre_new - grri_new == grre_old + grri_old, so the new I_v-based path is algebraically equivalent to the removed compute_surface_inductance_from_greens — this is a real refactor, not a silent behavior change. The compute_Iv boolean is a sensible perf win: it avoids the interior solve + Green's-function storage on the common Wv-only path.

Issues

1. VacuumData struct is missing the I_v field it's documented to need (src/ForceFreeStates/ForceFreeStatesStructs.jl). The PR removes grri/grre from VacuumData but doesn't add I_v, even though compute_vacuum_response!'s docstring (src/Vacuum/Vacuum.jl) lists I_v::AbstractMatrix{ComplexF64} as part of the required duck-typed contract "when compute_Iv=true". Today nothing breaks because the only call site against VacuumData (Free.jl:56) never passes compute_Iv=true, but the moment the in-place path is wired up for 3D (which the PR description flags as the actual motivation) it will throw has no field I_v. Worth adding the field now, or narrowing the docstring to say VacuumData doesn't support the I_v path yet.

2. Test gap that hides Port over splines from fortran to Julia #1: the "in-place compute_vacuum_response! matches wrapper" test in test/runtests_vacuum.jl uses a hand-built NamedTuple (with I_v already present) rather than an actual VacuumData, so it wouldn't catch the struct/docstring mismatch above. Also, calc_surface_inductance (new public function in ResponseMatrices.jl) has no direct unit test — only the DIIID regression case (reported in PR comments) exercises it end-to-end. A small targeted test (e.g. checking the μ₀(2π)² normalization / inv() relation on a synthetic I_v) would pin this down independent of the full pipeline.

3. Stale docs: docs/development/architecture.md (~lines 49, 133-134, 175) still describes grri/grre as the Vacuum module's outputs and as per-singular-surface data. This file wasn't touched by the PR and is now inconsistent with the I_v/compute_Iv design — docs/src/workflow.md and docs/src/vacuum.md were correctly updated, this one was missed.

4. Minor: broken comment in SingularCoupling.jl: "Extract the (m_res, m_res) diagonal element is needed for singflx" reads as a leftover half-edit of "only the (m_res,m_res) diagonal element is needed" — should be reworded (e.g. "Only the (m_res, m_res) diagonal element is needed for singflx").

Regression harness

Good to see this was already run and posted (diiid_n1, local vs develop) per the CLAUDE.md mandate — 47/48 quantities identical, one (PE surface energy) differs by 1.2e-10 absolute, attributed to Hermitization-ordering changes and explained as noise. That's a reasonable result. Since the SingularCoupling diagonal-coupling path (island half-widths, Chirikov parameter) is one of the more direct consumers of the new I_v → calc_surface_inductance route, it'd be worth a second look at whether any existing/new case stresses that path more than diiid_n1 does, but nothing in the report looks alarming.

Removed try/catch regularization in surface inductance

Agree with removing the ad hoc regularization/try-catch fallback in the old compute_surface_inductance_from_greens. Per the PR description this had no Fortran basis and was a one-off patch for a specific failing case. Fine to let inv(I_v) fail loudly if it's ever singular rather than silently returning a placeholder.

---

No pull request should be merged without human review. Per this repo's CLAUDE.md, this is non-negotiable — please make sure a third-party human reviewer signs off before merging into develop, regardless of any automated review feedback (including this one).

[logan-nc]

Will review weds during cleanup session

[logan-nc]

Great cleanup. One suggestion regarding the helper function, which left repeated 3-line chunks that could all be in a one-line helper.

[logan-nc]

This is a nice helper, but it seems unnecessarily small. All calls are proceeded by the same set of lines to set the vac_input at a psi, compute_vacuum_response, and then call this to invert Iv. the vac_input and Iv are never needed... It seems cleaner to have the helper function just take the psi of interest, do those two lines and this inversion all internally, and return the L matrix.

[logan-nc]

Regression harness after the review fixes + develop merge (diiid_n1, local working tree vs develop): 47 of 48 quantities unchanged (zero diff, all profile checksums identical). The single mover is the same one reported before:

Quantity develop local Diff Status
PE surface energy 5.826847e+00 5.826847e+00 1.161e-10 (0.00%) ** CHANGED **

This is the known Hermitization-ordering roundoff (Iᵛ is now Hermitized before inversion rather than L after, equal to O(ε²) in the anti-Hermitian discretization residual — now annotated at the call site in Vacuum.jl). Tests pass in runtests_vacuum.jl (incl. a new calc_surface_inductance unit test), and in runtests_resist_eval.jl/runtests_slayer_inputs.jl after the merge resolution.

Two things flagged during review for human eyes, no code change made:

1. Wall settings at singular surfaces (pre-existing):compute_singular_coupling_metrics! passes the caller's general wall_settings into the surface-inductance solve, whereas Fortran gpvacuum_flxsurf always uses no-wall/far-wall (Response.jl does force nowall). If the Fortran behavior is intended everywhere, the helper call in SingularCoupling.jl should hard-code nowall — reviewer's call.
2. Note for @jhalpern30: the 3D extension-point docstring in Vacuum.jl still describes the interior variant as -D + 2I, but the 2D compute_Iv path now uses D − 2I (overall sign flip). Implementing 3D from the note as written would give an opposite-sign I_v. Left untouched since the 3D path is your placeholder — flagging so the stub note can be updated whenever you bring it online.

[logan-nc]

Updated to address original review. Approved now