3D Vacuum Part 2: Passing structs instead of arrays into kernel, cleaning up spline interface, and matrix renaming

src/ForceFreeStates/ForceFreeStatesStructs.jl

@@ -470,60 +470,3 @@ end
 
 # Initialize function for OdeState with relevant parameters for array initialization
 OdeState(numpert_total::Int, numsteps_init::Int, numunorms_init::Int, msing::Int) = OdeState(; numpert_total, numsteps_init, numunorms_init, msing)
-
-# Below here are debug output structs used for benchmarking and unit testing
-
-
-
-
-"""
-    VacuumBenchmarkInputs
-
-A struct to hold all inputs required for vacuum benchmarking between Fortran and Julia implementations.
-
-## Fields
-
-  - `wv_block::Matrix{ComplexF64}` - Vacuum response matrix block
-  - `mpert::Int` - Number of poloidal modes
-  - `mtheta_eq::Int` - Number of poloidal grid points in input equilibrium (corresponds to `mtheta_eq` in VacuumInput)
-  - `mthvac::Int` - Number of poloidal grid points in vacuum calculations (corresponds to `mtheta` in VacuumInput)
-  - `complex_flag::Bool` - Flag indicating if complex arithmetic is used
-  - `kernelsign::Float64` - Sign of the kernel for vacuum calculation
-  - `wall_flag::Bool` - Flag indicating presence of wall
-  - `farwall_flag::Bool` - Flag indicating presence of far wall
-  - `grri::Matrix{Float64}` - Green's function response matrix
-  - `xzpts::Matrix{Float64}` - Coordinate points on plasma boundary [R, Z]
-  - `ahg_file::String` - Filename for AHG data
-  - `dir_path::String` - Directory path for input/output files
-  - `vac_inputs::Vacuum.VacuumInput` - VacuumInput struct for Julia vacuum code
-  - `wall_settings::Vacuum.WallShapeSettings` - Wall shape settings
-  - `n::Int` - Toroidal mode number
-  - `ipert_n::Int` - Index of perturbed toroidal mode
-  - `psifac::Float64` - Normalized flux coordinate
-"""
-@kwdef struct VacuumBenchmarkInputs
-    # Vacuum computation parameters
-    wv_block::Matrix{ComplexF64}
-    mpert::Int
-    mtheta_eq::Int
-    mthvac::Int
-    complex_flag::Bool
-    kernelsign::Float64
-    wall_flag::Bool
-    farwall_flag::Bool
-    grri::Matrix{Float64}
-    xzpts::Matrix{Float64}
-    ahg_file::String
-    dir_path::String
-
-    # VacuumInput struct for Julia code
-    vac_inputs::Vacuum.VacuumInput
-
-    # Wall settings
-    wall_settings::Vacuum.WallShapeSettings
-
-    # Additional context
-    n::Int
-    ipert_n::Int
-    psifac::Float64
-end

src/ForceFreeStates/Free.jl

@@ -28,28 +28,14 @@ function free_run!(odet::OdeState, ctrl::ForceFreeStatesControl, equil::Equilibr
     for ipert_n in 1:intr.npert
         # Set VACUUM run parameters and boundary shape
         n = ipert_n - 1 + intr.nlow
-        vac_inputs = compute_vacuum_inputs(intr.psilim, n, equil, intr, ctrl)
+        vac_inputs = Vacuum.VacuumInput(equil, intr.psilim, ctrl.mthvac, intr.mpert, intr.mlow, n; force_wv_symmetry=ctrl.force_wv_symmetry)
         fill!(vac.grri, 0.0)
         fill!(vac.grre, 0.0)
         fill!(vac.xzpts, 0.0)
 
         # Compute vacuum energy matrix and both Green's functions
         wv_block, vac.grri, vac.grre, vac.xzpts = Vacuum.compute_vacuum_response(vac_inputs, intr.wall_settings)
 
-        # Output data for unit testing and benchmarking
-        if intr.debug_settings.output_benchmark_data
-            @info "Outputting top level vacuum debug data for n = $n"
-            farwall_flag = intr.wall_settings.shape == "nowall" ? true : false
-            benchmark_inputs = VacuumBenchmarkInputs(
-                wv_block, intr.mpert, equil.control.mtheta, ctrl.mthvac,
-                true, vac_inputs.kernelsign, false,
-                farwall_flag, vac.grri, vac.xzpts, "ahg2msc_dcon.out", intr.dir_path,
-                vac_inputs, intr.wall_settings,
-                n, ipert_n, intr.psilim
-            )
-            @save "vacuum_response_inputs.jld2" benchmark_inputs
-        end
-
         # Equation 126 in Chance 1997 - scale by (m - n*q)(m' - n*q)
         singfac = collect(intr.mlow:intr.mhigh) .- (n * intr.qlim)
         @inbounds for ipert in 1:intr.mpert
@@ -129,66 +115,6 @@ function free_run!(odet::OdeState, ctrl::ForceFreeStatesControl, equil::Equilibr
     return vac
 end
 
-"""
-    compute_vacuum_inputs(psifac::Float64, n::Int, equil::Equilibrium.PlasmaEquilibrium, intr::ForceFreeStatesInternal, ctrl::ForceFreeStatesControl)
-
-Compute the necessary input paramters to the VACUUM code for computing the vacuum response matrix.
-Performs the same function as `free_write_msc` in the Fortran code, except we create a data struct
-to pass in memory to the Julia VACUUM code instead of writing to a file. Computed quantities include
-the r, z, and ν values at the plasma boundary, as well as mode numbers and number of poloidal points.
-
-### Arguments
-
-  - `psifac`: Flux surface value at the plasma boundary (Float64)
-  - `n`: Toroidal mode number (Int)
-  - `equil`: Plasma equilibrium data (Equilibrium.PlasmaEquilibrium)
-  - `intr`: Internal ForceFreeStates parameters
-  - `ctrl`: ForceFreeStates control parameters
-"""
-function compute_vacuum_inputs(psifac::Float64, n::Int, equil::Equilibrium.PlasmaEquilibrium, intr::ForceFreeStatesInternal, ctrl::ForceFreeStatesControl)
-
-    # Allocations
-    theta_norm = equil.rzphi_ys
-    mtheta = equil.config.mtheta + 1
-    angle = zeros(Float64, mtheta)
-    r = zeros(Float64, mtheta)
-    z = zeros(Float64, mtheta)
-    ν = zeros(Float64, mtheta)
-    rfac = zeros(Float64, mtheta)
-
-    # Compute r, z, and ν at the plasma boundary
-    qa = equil.profiles.q_spline(psifac)
-    for itheta in 1:mtheta
-        # Evaluate geometric quantities at (ψ, θ)
-        r2 = equil.rzphi_rsquared((psifac, theta_norm[itheta]))
-        offset = equil.rzphi_offset((psifac, theta_norm[itheta]))
-        nu_val = equil.rzphi_nu((psifac, theta_norm[itheta]))
-
-        rfac[itheta] = sqrt(r2)
-        angle[itheta] = 2π * (theta_norm[itheta] + offset)
-        ν[itheta] = nu_val
-    end
-    r .= equil.ro .+ rfac .* cos.(angle)
-    z .= equil.zo .+ rfac .* sin.(angle)
-
-    # Invert values for n < 0
-    if n < 0
-        ν .= -ν
-        n = -n
-    end
-
-    return Vacuum.VacuumInput(;
-        r=reverse(r),
-        z=reverse(z),
-        ν=reverse(ν),
-        mlow=intr.mlow,
-        mpert=intr.mpert,
-        n=n,
-        mtheta=ctrl.mthvac,
-        force_wv_symmetry=ctrl.force_wv_symmetry
-    )
-end
-
 """
     free_compute_wv_spline(ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium, intr::ForceFreeStatesInternal)
 
@@ -231,7 +157,7 @@ function free_compute_wv_spline(ctrl::ForceFreeStatesControl, equil::Equilibrium
         for ipert_n in 1:intr.npert
             # Compute vacuum matrix
             n = ipert_n - 1 + intr.nlow
-            vac_inputs = compute_vacuum_inputs(intr.psilim, n, ctrl, equil, intr)
+            vac_inputs = Vacuum.VacuumInput(equil, intr.psilim, intr.mtheta, intr.mpert, intr.mlow, n; force_wv_symmetry=ctrl.force_wv_symmetry)
             wv_block, _, _ = Vacuum.compute_vacuum_response(vac_inputs, intr.wall_settings)
 
             # Apply singular factor scaling

src/PerturbedEquilibrium/SingularCoupling.jl

@@ -28,28 +28,32 @@ Calculates the coupling between forcing modes and resonant surfaces, which deter
 the effectiveness of external perturbations at rational surfaces where q = m/n.
 
 Implements GPEC algorithm from gpout.f:
-1. For each forcing mode, apply unit field and compute plasma response
-2. For each singular surface, evaluate field jump and compute coupling quantities
-3. Calculate diagnostic island widths and overlap parameters
+
+ 1. For each forcing mode, apply unit field and compute plasma response
+ 2. For each singular surface, evaluate field jump and compute coupling quantities
+ 3. Calculate diagnostic island widths and overlap parameters
 
 ## Arguments
-- `state`: Output state structure to store results
-- `equil`: Equilibrium solution with q-profile and flux surfaces
-- `ForceFreeStates_results`: ForceFreeStates stability calculation results
-- `vac_data`: Vacuum response data including Green's functions
-- `ffs_intr`: ForceFreeStates internal state with singular surface data
-- `intr`: PerturbedEquilibrium internal state with permeability matrix
-- `ctrl`: Control parameters
+
+  - `state`: Output state structure to store results
+  - `equil`: Equilibrium solution with q-profile and flux surfaces
+  - `ForceFreeStates_results`: ForceFreeStates stability calculation results
+  - `vac_data`: Vacuum response data including Green's functions
+  - `ffs_intr`: ForceFreeStates internal state with singular surface data
+  - `intr`: PerturbedEquilibrium internal state with permeability matrix
+  - `ctrl`: Control parameters
 
 ## Modifies
+
 Populates the following fields in `state`:
-- `resonant_flux[msing, numpert_total]`: Normalized resonant flux Φ_r/A
-- `resonant_current[msing, numpert_total]`: Resonant current density
-- `island_width_sq[msing, numpert_total]`: Square of island half-width (w/2)²
-- `penetrated_field[msing, numpert_total]`: Normal field at resonant surface
-- `delta_prime[msing, numpert_total]`: Tearing stability parameter Δ'
-- `island_half_width[msing]`: Dimensional island half-width
-- `chirikov_parameter[msing]`: Island overlap metric
+
+  - `resonant_flux[msing, numpert_total]`: Normalized resonant flux Φ_r/A
+  - `resonant_current[msing, numpert_total]`: Resonant current density
+  - `island_width_sq[msing, numpert_total]`: Square of island half-width (w/2)²
+  - `penetrated_field[msing, numpert_total]`: Normal field at resonant surface
+  - `delta_prime[msing, numpert_total]`: Tearing stability parameter Δ'
+  - `island_half_width[msing]`: Dimensional island half-width
+  - `chirikov_parameter[msing]`: Island overlap metric
 
 Note: numpert_total = mpert × npert handles all (m,n) mode combinations
 """
@@ -113,14 +117,13 @@ function compute_singular_coupling_metrics!(
     end
 
     # Get vacuum calculation parameters
-    mtheta_eq = length(equil.rzphi_ys)  # Equilibrium poloidal grid size
     mtheta = vac_data.mthvac  # Vacuum poloidal grid size
     mlow = ffs_intr.mlow
     nlow = ffs_intr.nlow
     nhigh = ffs_intr.nhigh
 
     # Perturbed equilibrium calculations always use nowall
-    wall_settings = Vacuum.WallShapeSettings(shape="nowall")
+    wall_settings = Vacuum.WallShapeSettings(; shape="nowall")
 
     if ctrl.verbose
         println("  Processing toroidal modes n = $nlow:$nhigh")
@@ -156,16 +159,8 @@ function compute_singular_coupling_metrics!(
             end
 
             # Compute Green's functions at this surface for this n
-            vac_input = Vacuum.create_vacuum_input_at_psi(
-                equil,
-                sing_surf.psifac,
-                mtheta_eq,
-                mtheta,
-                mpert,
-                mlow,
-                nn
-            )
-            grri, grre = Vacuum.compute_greens_functions_only(vac_input, wall_settings)
+            vac_input = Vacuum.VacuumInput(equil, sing_surf.psifac, mtheta, mpert, mlow, nn)
+            _, grri, grre, _ = Vacuum.compute_vacuum_response(vac_input, wall_settings; green_only=true)
 
             # Store in singular surface struct (overwrites for each n)
             ffs_intr.sing[s].grri = grri
@@ -311,7 +306,7 @@ Interpolate magnetic field at arbitrary flux surface.
 
 Interpolates displacement from ForceFreeStates solution and converts to field using
 ideal MHD relation from flux coordinates:
-    b^ψ = i * χ₁ * (m - n*q) * ξ_ψ
+b^ψ = i * χ₁ * (m - n*q) * ξ_ψ
 
 This matches the field reconstruction in FieldReconstruction.jl.
 """
@@ -367,10 +362,10 @@ end
 Compute effective current density coefficient at given flux surface.
 
 Implements GPEC's j_c calculation (gpout.f line 560-578):
-    j_c = χ₁² * q / (μ₀ * integral)
+j_c = χ₁² * q / (μ₀ * integral)
 
 where the integral is computed via flux surface integration:
-    integral = ∫ (jac * |∇ψ| * sqreqb / |∇ψ|³) dθ
+integral = ∫ (jac * |∇ψ| * sqreqb / |∇ψ|³) dθ
 
 ## GPEC Formula
 
@@ -394,9 +389,10 @@ j_c(ising)=1.0/j_c(ising)*chi1**2*sq%f(4)/mu0
 
 Uses trapezoidal rule integration around the flux surface with metric quantities
 from the equilibrium bicubic spline (rzphi). The integrand includes:
-- jac: Jacobian of flux coordinates
-- |∇ψ|: Flux gradient magnitude (delpsi)
-- sqreqb: Magnetic field quantity (F² + χ₁²|∇ψ|²)/(2πR)²
+
+  - jac: Jacobian of flux coordinates
+  - |∇ψ|: Flux gradient magnitude (delpsi)
+  - sqreqb: Magnetic field quantity (F² + χ₁²|∇ψ|²)/(2πR)²
 """
 function compute_current_density(
     equil::Equilibrium.PlasmaEquilibrium,
@@ -436,8 +432,8 @@ function compute_current_density(
         r2 = equil.rzphi_rsquared((psi, theta))           # rfac²
         deta = equil.rzphi_offset((psi, theta))           # angle offset
         jac = equil.rzphi_jac((psi, theta))               # Jacobian
-        r2_y = equil.rzphi_rsquared((psi, theta); deriv=(0,1))  # ∂(rfac²)/∂theta
-        deta_y = equil.rzphi_offset((psi, theta); deriv=(0,1))  # ∂(deta)/∂theta
+        r2_y = equil.rzphi_rsquared((psi, theta); deriv=(0, 1))  # ∂(rfac²)/∂theta
+        deta_y = equil.rzphi_offset((psi, theta); deriv=(0, 1))  # ∂(deta)/∂theta
 
         rfac = sqrt(abs(r2))
         fy_rfac2 = r2_y
@@ -493,11 +489,13 @@ Apply Green's function to Fourier mode coefficients.
 Simplified version that extracts only plasma surface rows (first mtheta rows).
 
 ## Arguments
-- `green`: Green's function matrix [2*mtheta, 2*mpert]
-- `mode_coeffs`: Complex Fourier coefficients [mpert]
+
+  - `green`: Green's function matrix [2*mtheta, 2*mpert]
+  - `mode_coeffs`: Complex Fourier coefficients [mpert]
 
 ## Returns
-- Potential at plasma surface theta points [mtheta]
+
+  - Potential at plasma surface theta points [mtheta]
 """
 function apply_green_function_simple(
     green::Matrix{Float64},
@@ -509,7 +507,7 @@ function apply_green_function_simple(
     # Pack complex to real/imag format
     packed = zeros(Float64, 2 * mpert)
     for i in 1:mpert
-        packed[2*i - 1] = real(mode_coeffs[i])
+        packed[2*i-1] = real(mode_coeffs[i])
         packed[2*i] = imag(mode_coeffs[i])
     end
 
@@ -530,21 +528,25 @@ end
 Compute surface inductance matrix from Green's functions at flux surface.
 
 This implements the full GPEC gpvacuum_flxsurf algorithm (gpvacuum.f line 299-331):
-1. For each mode, apply Green's functions to unit flux
-2. Compute surface current from potential jump: kax = (chi - che) / μ₀
-3. Fourier transform to mode space
-4. Return inductance: L_surf = flux * inv(current)
+
+ 1. For each mode, apply Green's functions to unit flux
+ 2. Compute surface current from potential jump: kax = (chi - che) / μ₀
+ 3. Fourier transform to mode space
+ 4. Return inductance: L_surf = flux * inv(current)
 
 ## Arguments
-- `grri`: Interior Green's function [2*mtheta, 2*mpert] (stored in sing_surf)
-- `grre`: Exterior Green's function [2*mtheta, 2*mpert] (stored in sing_surf)
-- `ffs_intr`: ForceFreeStates internal state
-- `intr`: PerturbedEquilibrium internal state with mode arrays
+
+  - `grri`: Interior Green's function [2*mtheta, 2*mpert] (stored in sing_surf)
+  - `grre`: Exterior Green's function [2*mtheta, 2*mpert] (stored in sing_surf)
+  - `ffs_intr`: ForceFreeStates internal state
+  - `intr`: PerturbedEquilibrium internal state with mode arrays
 
 ## Returns
+
 Surface inductance matrix [numpert_total, numpert_total]
 
 ## GPEC Reference
+
 ```fortran
 DO i=1,mpert
    vbwp_mn=0; vbwp_mn(i)=1.0
@@ -650,7 +652,7 @@ end
 Compute singular flux from resonant current using surface inductance.
 
 Uses surface inductance matrix at the singular surface:
-    Φ = L_surf * I
+Φ = L_surf * I
 
 where L_surf is computed from Green's functions stored in sing_surf.
 
@@ -718,7 +720,7 @@ end
 Compute flux surface area at given ψ.
 
 Implements GPEC's area calculation (gpout.f line 568):
-    area = ∫ jac * |∇ψ| dθ
+area = ∫ jac * |∇ψ| dθ
 
 where the integral is computed around the flux surface.
 
@@ -740,8 +742,9 @@ area(ising)=area(ising)-jac*delpsi(mthsurf)/mthsurf  ! trapezoidal rule
 ## Implementation
 
 Uses trapezoidal rule integration around the flux surface with:
-- jac: Jacobian of flux coordinates from rzphi
-- |∇ψ|: Flux gradient magnitude (delpsi) from metric tensor
+
+  - jac: Jacobian of flux coordinates from rzphi
+  - |∇ψ|: Flux gradient magnitude (delpsi) from metric tensor
 """
 function compute_surface_area(
     equil::Equilibrium.PlasmaEquilibrium,
@@ -772,8 +775,8 @@ function compute_surface_area(
         r2 = equil.rzphi_rsquared((psi, theta))
         jac = equil.rzphi_jac((psi, theta))
         deta = equil.rzphi_offset((psi, theta))
-        r2_y = equil.rzphi_rsquared((psi, theta); deriv=(0,1))
-        deta_y = equil.rzphi_offset((psi, theta); deriv=(0,1))
+        r2_y = equil.rzphi_rsquared((psi, theta); deriv=(0, 1))
+        deta_y = equil.rzphi_offset((psi, theta); deriv=(0, 1))
 
         # Compute rfac
         rfac = sqrt(abs(r2))