VACUUM - Migrate surface current calculation into the VACUUM module

.claude/agent-memory/fortran-physics-reviewer/fortran_correspondence_map.md

@@ -15,7 +15,7 @@ relative to this repository.
 - `gpec/gpeq.f` (gpeq_sol, gpeq_contra, gpeq_surface, gpeq_normal) -> `src/PerturbedEquilibrium/FieldReconstruction.jl` + `src/PerturbedEquilibrium/ResponseMatrices.jl`
 - `gpec/gpresp.f` (gpresp_pinduct, gpresp_sinduct, gpresp_permeab) -> `src/PerturbedEquilibrium/ResponseMatrices.jl`
 - `gpec/gpout.f` (gpout_singcoup, gpout_xbnormal) -> `src/PerturbedEquilibrium/SingularCoupling.jl` + `src/PerturbedEquilibrium/FieldReconstruction.jl`
-- `gpec/gpvacuum.f` (gpvacuum_flxsurf) -> `src/PerturbedEquilibrium/SingularCoupling.jl` (compute_surface_inductance_from_greens)
+- `gpec/gpvacuum.f` (gpvacuum_flxsurf) -> `src/Vacuum/Vacuum.jl` (compute_Iv branch of _compute_vacuum_response_2d!) + `src/PerturbedEquilibrium/ResponseMatrices.jl` (calc_surface_inductance)
 
 ## KineticForces (NTV, Fortran `pentrc/`) — see kinetic_ntv_map.md for the audit checklist
 - `pentrc/torque.F90` -> `src/KineticForces/Torque.jl` (tpsi! single-surface torque)

docs/development/architecture.md

@@ -46,13 +46,15 @@ GPEC consists of **seven main modules** organized in `src/`:
 
 4. **Vacuum** (`src/Vacuum/`) - Vacuum field calculations and Green's functions
    - Computes vacuum response matrices for ideal MHD analysis
-   - Calculates both **interior** (grri) and **exterior** (grre) Green's functions
+   - Solves the exterior boundary-integral system for the vacuum energy matrix `wv`, and optionally
+     (`compute_Iv=true`) the interior system as well to build the surface-current matrix `I_v`
+     (Park 2007 eq. 21b). The interior/exterior Green's functions themselves are internal scratch.
    - Main functions:
-     - `compute_vacuum_response()` - Pure Julia implementation
+     - `compute_vacuum_response()` / `compute_vacuum_response!()` - allocating and in-place entry points
    - Key files:
-     - `VacuumStructs.jl` - Data structures
-     - `VacuumInternals.jl` - Core algorithms
-     - `VacuumFromEquilibrium.jl` - Integration with equilibrium data
+     - `DataTypes.jl` - Data structures (`VacuumInput`, `PlasmaGeometry`, `WallGeometry`)
+     - `Kernel2D.jl` / `Kernel3D.jl` - Single-/double-layer kernel assembly
+     - `Field.jl` - Vacuum field and potential evaluation off the surface
    - Status: **Pure Julia implementation complete and available**
 
 5. **ForceFreeStates** (`src/ForceFreeStates/`) - Ideal MHD stability analysis (DCON-style)
@@ -130,8 +132,7 @@ The complete GPEC analysis pipeline:
 
 2. **Vacuum Response**:
    - Initialize plasma and wall surfaces from equilibrium
-   - Compute vacuum response matrices (wv, grri, grre)
-   - Calculate both interior and exterior Green's functions
+   - Compute the vacuum energy matrix `wv` (and `I_v` when `compute_Iv=true`)
    - Pure Julia implementation
 
 3. **Stability Analysis** (ForceFreeStates):
@@ -169,11 +170,9 @@ The complete GPEC analysis pipeline:
 
 ### Stability
 - `SingType` - Singular surface data including:
-  - Rational surface location (ψ, q = m/n)
-  - Δ' (tearing stability parameter)
-  - Eigenmode structure at singular surface
-  - Green's functions (grri, grre) at interior singular surfaces
-  - Surface inductance
+  - Rational surface location (ψ, ρ, q = m/n, dq/dψ)
+  - Δ' (tearing stability parameter) — **stub**; the valid Δ' is `ForceFreeStatesInternal.delta_prime_matrix`
+  - Asymptotic solution bases at the inner-layer boundaries
 
 ### Perturbed Equilibrium
 - `PerturbedEquilibriumControl` - User-facing TOML configuration parameters

docs/src/vacuum.md

@@ -73,7 +73,7 @@ wall_settings = GeneralizedPerturbedEquilibrium.Vacuum.WallShapeSettings(
     equal_arc_wall = true     # Use equal arc length spacing
 )
 
-# Compute vacuum response; returns a VacuumResponse with wv, grri, grre, plasma_pts, wall_pts
+# Compute vacuum response; returns a VacuumResponse with wv, I_v, plasma_pts, wall_pts
 vac = GeneralizedPerturbedEquilibrium.Vacuum.compute_vacuum_response(inputs, wall_settings)
 ```
 

docs/src/workflow.md

@@ -59,8 +59,7 @@ The single `gpec.toml` file supplies user-selected options to every module. The
 
 **Outputs**:
 - `wv` — Vacuum response matrix (scaled by the singular factor (m - nq)(m' - nq), see Chance 1997)
-- `grri` — Interior Green's function matrix (plasma boundary → plasma boundary)
-- `grre` — Exterior Green's function matrix (plasma boundary → wall)
+- `I_v` — Vacuum surface-current matrix Iᵛ when `compute_Iv=true` (otherwise zeros); PerturbedEquilibrium inverts this to surface inductance `L`
 
 **Key references**: [Chance et al. (1997)](citations.md#Vacuum-Module), [Chance et al. (2007)](citations.md#Vacuum-Module)
 

src/ForceFreeStates/ForceFreeStatesStructs.jl

@@ -11,8 +11,6 @@ A mutable struct holding data related to the singular surfaces in the equilibriu
   - `n::Vector{Int}` - Toroidal mode number(s)
   - `q::Float64` - Safety factor (= m/n)
   - `q1::Float64` - Derivative of safety factor with respect to ψ
-  - `grri::Array{ComplexF64,2}` - Interior Green's function at this surface [mthvac, mpert]
-  - `grre::Array{ComplexF64,2}` - Exterior Green's function at this surface [mthvac, mpert]
   - `delta_prime::Vector{ComplexF64}` - **STUB (not physically valid)**. Per-surface ca-based Δ' estimate retained for future work / debugging only. The physically valid Δ' is `ForceFreeStatesInternal.delta_prime_matrix`, computed via the STRIDE global BVP (Glasser 2018 PoP 25, 032501). Do not use this field for tearing-stability analysis; do not expect agreement with `delta_prime_matrix`.
   - `delta_prime_col::Matrix{ComplexF64}` - **STUB (not physically valid)**. Per-surface ca-based Δ' column retained for future work / debugging only. Shape (numpert_total × n_res_modes); `delta_prime_col[j, i] = (ca_r[j,ipert_res_i,2] - ca_l[j,ipert_res_i,2]) / (4π²·psio)`. The diagonal element matches the (also stubbed) `delta_prime[i]`. Only populated for the Riccati/parallel FM paths. The physically valid Δ' is `ForceFreeStatesInternal.delta_prime_matrix`; this field exists for future development on intra-surface coupling diagnostics, not for production use.
 """
@@ -23,8 +21,6 @@ A mutable struct holding data related to the singular surfaces in the equilibriu
     n::Vector{Int} = Int[]
     q::Float64 = 0.0
     q1::Float64 = 0.0
-    grri::Array{ComplexF64,2} = Array{ComplexF64}(undef, 0, 0)
-    grre::Array{ComplexF64,2} = Array{ComplexF64}(undef, 0, 0)
     delta_prime::Vector{ComplexF64} = ComplexF64[]
     delta_prime_col::Matrix{ComplexF64} = Matrix{ComplexF64}(undef, 0, 0)
     ua_left::Array{ComplexF64,3} = Array{ComplexF64}(undef, 0, 0, 0)   # asymptotic basis at left inner-layer boundary
@@ -540,35 +536,34 @@ and a small set of temporary matrices and factors used to compute singular-layer
 
   - `ising_start::Int` - Index of the starting singular surface to be crossed during integration.
 
-    # Initialization parameters
-
   - `psimax::Float64` - Maximum psi value for which the integrator is allowed to run in next integration region.
 
   - `needs_crossing::Bool` - Flag indicating whether a rational surface needs to be crossed after the current integration region.
 
   - `nzero::Int` - Count of detected zero crossings (used for diagnostics).
 
-    # Saved data throughout integration
-
   - `new::Bool` - Flag indicating whether a new `unorm0` should be computed after a fixup.
 
-# Total ODE solver steps taken (all steps, not just saved ones)
+    # Initialization parameters
 
   - `unorm::Vector{Float64}` - Current norms of the solution vectors (length `numpert_total`).
 
   - `unorm0::Vector{Float64}` - Reference/initial norms of the solution vectors (length `numpert_total`).
 
+    # Saved data throughout integration
+
   - `ifix::Int` - Number of normalization operations performed (index into normalization arrays).
 
+# Total ODE solver steps taken (all steps, not just saved ones)
+
   - `index::Array{Int,2}` - Index matrix used for sorting solution norms with shape `(numpert_total, numunorms_init)`.
 
   - `sing_flag::Vector{Bool}` - Boolean flags indicating which stored normalizations correspond to singular solutions    # Edge dW scan state and results (disabled sentinel when psiedge >= psilim, i.e. no edge scan)
     (length `numunorms_init`).
 
   - `zeroed_idx::Vector{Vector{Int}}` - For each ideal rational surface jump, a vector of indices of solutions that were zeroed.    # Data for integrator
 
-  - `fixfac::Array{ComplexF64,3}` - Fix-up factors for Gaussian reduction with shape    # Initialization parameters
-    `(numpert_total, numpert_total, numunorms_init)`.
+  - `fixfac::Array{ComplexF64,3}` - Fix-up factors for Gaussian reduction with shape `(numpert_total, numpert_total, numunorms_init)`.
 
   - `fixstep::Vector{Int64}` - Step indices (psi step positions) at which normalization/fixups were performed (length `numunorms_init`).
 """

src/PerturbedEquilibrium/Response.jl

@@ -36,16 +36,9 @@ function compute_plasma_response!(
     # Plasma inductance Lambda (wt0 formula, Fortran resp_induct_flag=TRUE default)
     plasma_inductance = calc_plasma_inductance(wt0, ffs_intr, equil.psio)
 
-    # Surface inductance L from Green's functions at psilim.
-    # Requires a 2D (nzvac=1) vacuum response so rows are theta points only,
+    # Surface inductance L from vacuum surface-current matrix at psilim.
     nn = ffs_intr.nlow
-    vac_input_2d = Vacuum.VacuumInput(equil, ffs_intr.psilim, mthvac, 1, ffs_intr.mlow:ffs_intr.mhigh, [nn])
-    wall_nowall = Vacuum.WallShapeSettings(; shape="nowall")
-    vac_2d = Vacuum.compute_vacuum_response(vac_input_2d, wall_nowall)
-    grri_2d = Matrix{ComplexF64}(vac_2d.grri)
-    grre_2d = Matrix{ComplexF64}(vac_2d.grre)
-    ν_vac = Vacuum.PlasmaGeometry(vac_input_2d).ν
-    surface_inductance = compute_surface_inductance_from_greens(grri_2d, grre_2d, ffs_intr, nn, ν_vac)
+    surface_inductance = calc_surface_inductance(equil, ffs_intr.psilim, mthvac, ffs_intr.mlow:ffs_intr.mhigh, nn)
     permeability = calc_permeability(plasma_inductance, surface_inductance)
 
     # Reluctance ϱ = L⁻¹·(Λ† − L)·L⁻¹ (Fortran gpresp_reluct: diff_indmats = CONJG(TRANSPOSE(plas_indmats)) − surf_indmats).

src/PerturbedEquilibrium/ResponseMatrices.jl

@@ -111,10 +111,12 @@ where ξ_ψ is the radial displacement eigenfunction.
 ## Arguments
 
   - `boundary_data`: Output from extract_boundary_displacements()
+
       + ξ_psi_boundary: Boundary displacement [numpert_total, numpert_total]
       + dPsi_drho: Flux surface spacing at boundary (scalar)
       + q_boundary: Safety factor at boundary (scalar)
       + psi_boundary: Normalized flux at boundary (scalar)
+
   - `intr`: ForceFreeStates internal state with mode arrays (mlow, mhigh, nlow, etc.)
 
 ## Returns
@@ -262,6 +264,48 @@ function calc_plasma_inductance(
     return inv(temp2)
 end
 
+"""
+    calc_surface_inductance(
+        equil::Equilibrium.PlasmaEquilibrium,
+        ψ::Float64,
+        mtheta::Int,
+        m_modes::AbstractUnitRange{Int},
+        nn::Int
+    )::Matrix{ComplexF64}
+
+Surface inductance L at the flux surface ψ, from the vacuum surface-current matrix, solved
+against a no-wall vacuum.
+
+Solves the 2D vacuum problem with `compute_Iv=true` and inverts. The driving flux harmonics are
+unit columns (`Φ = 𝕀` by construction, matching Fortran `gpvacuum_flxsurf`'s unit driving
+harmonics), so `Φ = L·I^v` (Park 2007, eq. 7 and following text) gives `L = μ₀(2π)²·I_v⁻¹`.
+
+## Arguments
+
+  - `equil`: Plasma equilibrium
+  - `ψ`: Normalized poloidal flux of the surface
+  - `mtheta`: Number of vacuum poloidal grid points
+  - `m_modes`: Poloidal mode range mlow:mhigh
+  - `nn`: Toroidal mode number
+
+## Returns
+
+  - Surface inductance matrix in henries [mpert × mpert]
+"""
+function calc_surface_inductance(
+    equil::Equilibrium.PlasmaEquilibrium,
+    ψ::Float64,
+    mtheta::Int,
+    m_modes::AbstractUnitRange{Int},
+    nn::Int
+)::Matrix{ComplexF64}
+    vac_input = Vacuum.VacuumInput(equil, ψ, mtheta, 1, m_modes, [nn])
+    wall_settings = Vacuum.WallShapeSettings(; shape="nowall")
+    I_v = Vacuum.compute_vacuum_response(vac_input, wall_settings; compute_Iv=true).I_v
+    μ₀ = 4π * 1e-7
+    return inv(I_v) .* (μ₀ * (2π)^2)
+end
+
 """
     calc_permeability(
         plasma_inductance::Matrix{ComplexF64},

src/PerturbedEquilibrium/SingularCoupling.jl

@@ -51,9 +51,6 @@ function _hermite_cubic_deriv(u_a, u_b, du_a, du_b, psi_a, psi_b, psi)
     return @. d00 * u_a + d10 * du_a + d01 * u_b + d11 * du_b
 end
 
-# Reflect a periodic theta-space vector θ → -θ (the theta reversal in gpvacuum_flxsurf).
-_reverse_theta(v::AbstractVector) = circshift(reverse(v), 1)
-
 """
     _chord_solution_at(psi, resnum, odet, nstep) -> (u, du)
 
@@ -137,7 +134,7 @@ function _solution_at(
 
     # 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]
+    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
@@ -280,7 +277,6 @@ function compute_singular_coupling_metrics!(
     chi1 = 2π * equil.psio
     twopi = 2π
     mtheta = mthvac
-    wall_settings = Vacuum.WallShapeSettings(; shape="nowall")
 
     # Phase 1: Collect all resonant (surface, n) pairs in psi order
     resonant_pairs = Tuple{Int,Int}[]
@@ -362,17 +358,10 @@ function compute_singular_coupling_metrics!(
                 continue
             end
 
-            # Compute Green's functions at this surface for this n (once per pair)
-            vac_input = Vacuum.VacuumInput(equil, sing_surf.psifac, mtheta, 1, mlow:mhigh, [nn])
-            vac = Vacuum.compute_vacuum_response(vac_input, wall_settings)
-            grri = Matrix{ComplexF64}(vac.grri)
-            grre = Matrix{ComplexF64}(vac.grre)
-
-            # Get ν on the vacuum theta grid (same ν used in the vacuum Fourier basis computation)
-            ν_vac = Vacuum.PlasmaGeometry(vac_input).ν
+            # Surface inductance at this surface for this n (once per pair)
+            L_surf = calc_surface_inductance(equil, sing_surf.psifac, mtheta, mlow:mhigh, nn)
 
-            # Precompute L_surf; only the (m_res, m_res) diagonal element is needed for singflx
-            L_surf = compute_surface_inductance_from_greens(grri, grre, ffs_intr, nn, ν_vac)
+            # Only the (m_res, m_res) diagonal element is needed for singflx
             m_idx = m_res - mlow + 1
             L_mm = L_surf[m_idx, m_idx]
 
@@ -614,99 +603,6 @@ function compute_current_density(
     return j_c
 end
 
-"""
-    compute_surface_inductance_from_greens(
-        grri::Matrix{ComplexF64},
-        grre::Matrix{ComplexF64},
-        ffs_intr::ForceFreeStatesInternal,
-        nn::Int,
-        ν::Vector{Float64}
-    )::Matrix{ComplexF64}
-
-Compute surface inductance matrix from Green's functions at flux surface.
-
-Implements the GPEC `gpvacuum_flxsurf` algorithm.
-
-The Julia vacuum code uses SFL Fourier basis `cos(m*θ - n*ν)` in the column transform,
-so the row DFT must apply the matching toroidal phase correction `exp(-i*n*ν)` before
-the DFT (matching Fortran `gpvacuum_flxsurf`'s `EXP(-ifac*nn*dphi)` phase correction).
-
-## Arguments
-
-  - `grri`: Interior Green's function [mtheta, mpert]
-  - `grre`: Exterior Green's function [mtheta, mpert]
-  - `ffs_intr`: ForceFreeStates internal state
-  - `nn`: Toroidal mode number
-  - `ν`: Toroidal angle offset on the vacuum theta grid [mtheta]
-
-## Returns
-
-Surface inductance matrix [mpert × mpert]
-"""
-@with_pool pool function compute_surface_inductance_from_greens(
-    grri::Matrix{ComplexF64},
-    grre::Matrix{ComplexF64},
-    ffs_intr::ForceFreeStatesInternal,
-    nn::Int,
-    ν::Vector{Float64}
-)::Matrix{ComplexF64}
-    mpert = ffs_intr.mpert
-    mtheta = length(ν)
-    μ₀ = 4π * 1e-7
-
-    ft = FourierTransforms.FourierTransform(mtheta, mpert, ffs_intr.mlow)
-
-    flux_matrix = zeros!(pool, ComplexF64, mpert, mpert)
-    current_matrix = zeros!(pool, ComplexF64, mpert, mpert)
-
-    kax = zeros!(pool, ComplexF64, mtheta)
-    grri_surf = @view grri[1:mtheta, :]
-    grre_surf = @view grre[1:mtheta, :]
-
-    # Toroidal phase correction: exp(-i*n*ν)
-    phase = cis.(-nn .* ν)
-
-    for i in 1:mpert
-        flux_matrix[i, i] = 1.0
-
-        # Complex grri/e stores exp(i(mθ-nν)) projection, need conjugate for exp(-i(mθ-nν))
-        kax .= conj.(grri_surf[:, i] .+ grre_surf[:, i]) ./ (μ₀ * (2π)^2)
-
-        # Port of Fortran gpvacuum_flxsurf: apply toroidal phase, reverse theta, forward-DFT.
-        g_phased = kax .* phase
-        current_matrix[:, i] = ft(_reverse_theta(g_phased))
-    end
-
-    # Compute surface inductance: L_surf = flux * inv(current) = inv(current)
-    L_surf = zeros(ComplexF64, mpert, mpert)
-
-    current_mag = maximum(abs.(current_matrix))
-
-    if current_mag < 1e-15
-        @warn "Current matrix is all zeros! Cannot compute surface inductance." maxlog=1
-        for i in 1:mpert
-            L_surf[i, i] = μ₀ * 1e-6
-        end
-    else
-        try
-            regularization = 1e-12 * current_mag
-            current_reg = current_matrix + regularization * I
-
-            L_surf = flux_matrix * inv(current_reg)
-
-            # Hermitianize (matches Fortran: temp1 = 0.5*(temp1 + CONJG(TRANSPOSE(temp1))))
-            L_surf = 0.5 * (L_surf + L_surf')
-        catch e
-            @warn "Surface inductance inversion failed: $e" maxlog=1
-            for i in 1:mpert
-                L_surf[i, i] = μ₀ * 1e-6
-            end
-        end
-    end
-
-    return L_surf
-end
-
 """
     compute_surface_area(
         equil::Equilibrium.PlasmaEquilibrium,