PE - Removing AI slop from PerturbedEquilibrium module

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

@@ -74,7 +74,7 @@ wall_settings = GeneralizedPerturbedEquilibrium.Vacuum.WallShapeSettings(
 )
 
 # Compute vacuum response matrix
-wv, grri, xzpts = GeneralizedPerturbedEquilibrium.Vacuum.compute_vacuum_response(inputs, wall_settings)
+wv, _, _, _ = GeneralizedPerturbedEquilibrium.Vacuum.compute_vacuum_response(inputs, wall_settings)
 ```
 
 ### Vacuum Field Calculation at Observation Points

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)
 
@@ -182,5 +181,5 @@ All results are written to a single HDF5 file (default: `gpec.h5`). The file is
 | `locstab/` | Local stability: Mercier criterion, shear |
 | `integration/` | ODE integration results: energy matrices, eigenvalues |
 | `singular/` | Per-surface data: ψ_s, m/n, Δ', small solution coefficients |
-| `vacuum/` | Vacuum response matrices: wv, grri, grre |
+| `vacuum/` | Vacuum response matrices: wv (and I_v when computed for PE) |
 | `perturbed/` | Perturbed equilibrium: ξ, b in mode space, island diagnostics |

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
@@ -422,8 +418,6 @@ Populated in `Free.jl`.
   - `et::Vector{ComplexF64}` - Total energy eigenvalues of the pencil (W, N): power-normalized and invariant to the working (Jacobian) coordinate; et = ep + ev per mode
   - `n_tor_idx::Vector{Int}` -  0-based toroidal mode number index of each sorted eigenvalue (numpert_total). Needed in `write_imas`
   - `vacuum_eigenvalue::Float64` - Least stable (minimum) eigenvalue of the pencil (wv, N), clamped to zero
-  - `grri::Array{ComplexF64, 2}` - Interior Green's function matrices (2 * mthvac * nzvac × numpert_total)
-  - `grre::Array{ComplexF64, 2}` - Exterior Green's function matrices (2 * mthvac * nzvac × numpert_total)
   - `plasma_pts::Array{Float64, 3}` - Cartesian coordinates of plasma points, shape (mthvac * nzvac) × 3 for (x, y, z)
   - `wall_pts::Array{Float64, 3}` - Cartesian coordinates of wall points, shape (mthvac * nzvac) × 3 for (x, y, z)
 """
@@ -441,8 +435,6 @@ Populated in `Free.jl`.
     et::Vector{ComplexF64} = Vector{ComplexF64}(undef, numpert_total)
     n_tor_idx::Vector{Int} = zeros(Int, numpert_total)
     vacuum_eigenvalue::Float64 = NaN
-    grri::Array{ComplexF64,2} = Array{ComplexF64}(undef, 2 * numpoints, numpert_total)
-    grre::Array{ComplexF64,2} = Array{ComplexF64}(undef, 2 * numpoints, numpert_total)
     plasma_pts::Array{Float64,2} = Array{Float64}(undef, numpoints, 3)
     wall_pts::Array{Float64,2} = Array{Float64}(undef, numpoints, 3)
 end

src/ForceFreeStates/Free.jl

@@ -175,7 +175,7 @@ q-window minimum.
 
         # Compute raw vacuum matrix at the actual scan psi (singfac NOT applied; free_compute_total applies it analytically)
         vac_inputs = Vacuum.VacuumInput(equil, psi_array[i], ctrl.mthvac, ctrl.nzvac, intr.mlow:intr.mhigh, intr.nlow:intr.nhigh)
-        wv, _, _, _, _ = Vacuum.compute_vacuum_response(vac_inputs, intr.wall_settings)
+        wv, _, _, _ = Vacuum.compute_vacuum_response(vac_inputs, intr.wall_settings)
         @views wv_array[i, :, :] .= wv
     end
 

src/GeneralizedPerturbedEquilibrium.jl

@@ -461,12 +461,10 @@ function main_from_inputs(
             @info "Computing free boundary energies ($wall_desc)"
         end
         vac_data = free_run!(odet, ctrl, equil, ffit, intr)
-        if real(vac_data.et[1]) < 0
-            if ctrl.verbose
+        if ctrl.verbose
+            if real(vac_data.et[1]) < 0
                 @warn "Free-boundary mode unstable for n = $nstring"
-            end
-        else
-            if ctrl.verbose
+            else
                 @info "All free-boundary modes stable for n = $nstring"
             end
         end
@@ -554,6 +552,30 @@ function main_from_inputs(
             slayer=slayer_result)
     end
 
+    # No perturbed equilibrium calculations if vacuum data is not available
+    if !ctrl.vac_flag
+        @warn "Vacuum data not available. Skipping perturbed equilibrium calculations. Set vac_flag=true in [ForceFreeStates] section."
+        return
+    end
+
+    # No perturbed equilibrium calculations if free-boundary mode is unstable
+    if real(vac_data.et[1]) < 0
+        @warn "Since a free-boundary mode is unstable, perturbed equilibrium calculations will not run."
+        return
+    end
+
+    # No perturbed equilibrium calculations if wall shape is not nowall
+    if intr.wall_settings.shape !== "nowall"
+        @warn "Perturbed equilibrium calculations are only supported for nowall cases - skipping perturbed equilibrium calculations."
+        return
+    end
+
+    # Perturbed equilibrium calculations still need multi-n support added
+    if intr.npert > 1
+        @warn "Perturbed equilibrium calculations are only supported for single-n cases - skipping perturbed equilibrium calculations."
+        return
+    end
+
     # ----------------------------------------------------------------
     # Perturbed Equilibrium
     # ----------------------------------------------------------------
@@ -672,17 +694,11 @@ function main_from_inputs(
     end
     slayer_result = _run_slayer_stage(pe_file)
 
-    # ----------------------------------------------------------------
-    # Done
-    # ----------------------------------------------------------------
     @info "\n$_BANNER\n  GPEC completed successfully in $(@sprintf("%.3f", time() - total_start)) s\n$_BANNER"
 
-    # TODO: Do not allow perturbed equilibrium calculations if zero crossings are found
-
     return (ctrl=ctrl, equil=equil, intr=intr, ffit=ffit, odet=odet,
         vac_data=ctrl.vac_flag ? vac_data : nothing,
         slayer=slayer_result)
-
 end
 
 """
@@ -693,10 +709,6 @@ This combines the functionality of several pieces of the Fortran code in `ode_ou
 primarily `ode_output_open` and the various `bin_euler` writes that occur throughout the
 integration. Some parameters are only dumped in their respective flags are true, e.g.
 vacuum data if `vac_flag` is true.
-
-### TODOs
-
-Combine spline unpacking if possible, too many extra lines
 """
 function write_outputs_to_HDF5(
     ctrl::ForceFreeStatesControl,

src/PerturbedEquilibrium/PerturbedEquilibrium.jl

@@ -1,5 +1,7 @@
 module PerturbedEquilibrium
 
+const μ0 = 4π * 1e-7
+
 # Imports
 using HDF5
 using Printf
@@ -78,7 +80,7 @@ function compute_perturbed_equilibrium(
 
     state = PerturbedEquilibriumState()
 
-    # Step 0: Initialize mode arrays for convenient indexing
+    # Initialize mode arrays for convenient indexing
     initialize_mode_arrays!(intr, ffs_intr)
 
     # Load forcing data. On the gpec.h5 replay path the caller preloads
@@ -113,27 +115,16 @@ function compute_perturbed_equilibrium(
         end
     end
 
-    # Step 2: Compute plasma response
+    # Compute plasma response
     if ctrl.compute_response
-        if vac_data === nothing
-            @warn "Vacuum data not available. Skipping plasma response calculation. Set vac_flag=true in [ForceFreeStates] section."
-        else
-            compute_plasma_response!(state, equil, ForceFreeStates_results, vac_data, ffs_intr, intr, ctrl, metric, ffit)
-        end
+        compute_plasma_response!(state, equil, ForceFreeStates_results, vac_data, ffs_intr, intr, ctrl, metric, ffit)
     end
 
-    # Step 3: Compute singular coupling metrics
+    # Compute singular coupling metrics
     if ctrl.compute_singular_coupling
-        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)
-        end
+        compute_singular_coupling_metrics!(state, equil, ForceFreeStates_results, vac_data, ffs_intr, intr, ctrl)
     end
 
-    # Step 4: Output eigenmode fields (integrated into HDF5 output)
-    # This is handled by write_outputs_to_HDF5 in main()
-
     return state
 end
 

src/PerturbedEquilibrium/Response.jl

@@ -25,40 +25,31 @@ function compute_plasma_response!(
     metric::MetricData,
     ffit::FourFitVars
 )
-    if ctrl.verbose
-        @info "Computing plasma response (wt0-based inductance)"
-    end
+    ctrl.verbose && @info "Computing plasma response"
 
     # Build flux matrix from ForceFreeStates eigenmodes [mode × eigenmode]
-    flux_matrix = build_flux_matrix(equil, ForceFreeStates_results, vac_data, ffs_intr)
+    flux_matrix = build_flux_matrix(equil, ForceFreeStates_results, ffs_intr)
 
-    # Plasma inductance Lambda (wt0 formula, Fortran resp_induct_flag=TRUE default)
-    plasma_inductance = calc_plasma_inductance(vac_data, ffs_intr, equil.psio)
+    # Compute plasma inductance
+    plasma_inductance = calc_plasma_inductance(ffs_intr, vac_data.wt0, 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, vac_data.mthvac, 1, ffs_intr.mlow:ffs_intr.mhigh, [nn])
-    wall_nowall = Vacuum.WallShapeSettings(; shape="nowall")
-    _, grri_2d_raw, grre_2d_raw, _, _ = Vacuum.compute_vacuum_response(vac_input_2d, wall_nowall)
-    grri_2d = Matrix{ComplexF64}(grri_2d_raw)
-    grre_2d = Matrix{ComplexF64}(grre_2d_raw)
-    ν_vac = Vacuum.PlasmaGeometry(vac_input_2d).ν
-    surface_inductance = compute_surface_inductance_from_greens(grri_2d, grre_2d, ffs_intr, nn, ν_vac)
-    permeability = calc_permeability(plasma_inductance, surface_inductance)
-
-    # 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
+    _, I_v, _, _ = Vacuum.compute_vacuum_response(vac_input_2d, wall_nowall; compute_Iv=true)
+    surface_inductance = calc_surface_inductance(I_v)
 
-    # Store permeability in internal state for singular coupling / field reconstruction.
-    # These consumers operate on the physical control-surface flux Φ_x, so the internal
-    # copy stays in flux space; only the stored/output quantities are conformed to fields below.
+    # Compute permeability P = Λ·L⁻¹ and store in internal state for singular coupling / field reconstruction.
+    permeability = plasma_inductance / surface_inductance
     intr.plasma_response = permeability
 
+    # Compute reluctance ϱ = L⁻¹·(Λ† − L)·L⁻¹
+    # Λ 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
+
     # Conform the control-surface matrices to the coordinate-invariant root-area-weighted
-    # field (b̃) space for output (issue #233 / Pharr 2026). Store the b̃→b̄ operator S = Σ/√A
+    # field (b̃) space for output (Pharr 2026). Store the b̃→b̄ operator S = Σ/√A
     # and the scalar surface area A so users can recover the area-weighted field (b̄ = S·b̃) or
     # flux (Φ = A·b̄) — see Utils.jl output docs.
     rootarea_to_area_weight, surface_area = build_control_surface_rootarea_to_area_weight(equil, ffs_intr)
@@ -70,11 +61,9 @@ function compute_plasma_response!(
     state.rootarea_to_area_weight = rootarea_to_area_weight
     state.surface_area = surface_area
 
-    # Forcing and response on the control surface. Flux Φ appears only as a brief internal bridge:
-    # forcing arrives as Φ_x, the field reconstruction below consumes Φ_tot, and the b̃ spectra are
-    # formed via the conform operator R = S·A (Φ = R·b̃).
+    # Compute actual flux from external flux and permiability, Φ = P Φ^x
     forcing_flux = map_forcing_to_eigenmodes(intr.forcing_modes, ffs_intr)
-    response_flux = compute_plasma_response_vector(permeability, forcing_flux)
+    response_flux = permeability * forcing_flux
 
     # Output forcing/response in the three Pharr field representations (all tesla):
     #   b̃ (root-area-weighted) = R⁻¹·Φ,  b (bare) = Σ⁻¹·b̃,  b̄ (area-weighted) = S·b̃.
@@ -119,7 +108,5 @@ function compute_plasma_response!(
     state.b_n_modes = b_n_modes
     state.xi_n_modes = xi_n_modes
 
-    if ctrl.verbose
-        @info "Response complete: $(length(intr.forcing_modes)) forcing modes, max amplitude = $(@sprintf("%.3e", maximum(abs.(response_flux))))"
-    end
+    ctrl.verbose && @info "Response complete: $(length(intr.forcing_modes)) forcing modes, max amplitude = $(@sprintf("%.3e", maximum(abs.(response_flux))))"
 end