Merge pull request #306 from OpenFUSIONToolkit/feature/3Dvac_nfp

VACUUM (3D) - taking advantage of field period symmetry for nonaxisymmetric vacuum calculations

Author: jhalpern30

benchmarks/benchmark_fourier_transforms.jl

@@ -54,7 +54,7 @@ w(θ) = √(J·|∇ψ|).
 
 Operationally, sqrtamat is the mode-space √weight operator: for a field b with
 Fourier coefficients b_fft, it satisfies the identity
-  `‖sqrtamat·b_fft‖² = N² · ∫ |b|² · J|∇ψ| dθ`
+`‖sqrtamat·b_fft‖² = N² · ∫ |b|² · J|∇ψ| dθ`
 which is Jacobian-invariant on a given flux surface (see
 `scripts/test_power_norm_invariance.jl`).
 
@@ -78,13 +78,8 @@ function compute_sqrtamat(
         e_k .= 0.0
         e_k[k] = 1.0 + 0.0im
 
-        # Standard backward FT: f(θ_j) = (1/N) Σ_m c_m exp(-imθ_j)
-        # exp(-imθ) = cos(mθ) - i·sin(mθ), so:
-        #   Re(f) = (1/N)(cslth·Re(c) + snlth·Im(c))
-        #   Im(f) = (1/N)(cslth·Im(c) - snlth·Re(c))
-        real_part = (ft.cslth * real.(e_k) .+ ft.snlth * imag.(e_k)) ./ mtheta
-        imag_part = (ft.cslth * imag.(e_k) .- ft.snlth * real.(e_k)) ./ mtheta
-        theta_vec = complex.(real_part, imag_part)
+        # Standard backward FT: f(θ_j) = (1/N) Σ_m c_m exp(-imθ_j) = (transpose(basis) * c) / N
+        theta_vec = (transpose(ft.basis) * e_k) ./ mtheta
 
         # Multiply pointwise by √(J·|∇ψ|) in theta-space
         theta_vec .*= sqrt_jdp

src/Equilibrium/CoordinateInvariant.jl

@@ -11,8 +11,8 @@ 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{Float64,2}` - Interior Green's function at this surface [2*mthvac, 2*mpert]
-  - `grre::Array{Float64,2}` - Exterior Green's function at this surface [2*mthvac, 2*mpert]
+  - `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 +23,8 @@ 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{Float64,2} = Array{Float64}(undef, 0, 0)
-    grre::Array{Float64,2} = Array{Float64}(undef, 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
@@ -410,8 +410,8 @@ 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{Float64, 2}` - Interior Green's function matrices (2 * mthvac * nzvac × 2 * numpert_total)
-  - `grre::Array{Float64, 2}` - Exterior Green's function matrices (2 * mthvac * nzvac × 2 * numpert_total)
+  - `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)
 """
@@ -429,9 +429,8 @@ 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{Float64,2} = Array{Float64}(undef, 2 * numpoints, 2 * numpert_total)
-    grre::Array{Float64,2} = Array{Float64}(undef, 2 * numpoints, 2 * numpert_total)
+    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/ForceFreeStatesStructs.jl

@@ -54,7 +54,7 @@ and data dumping.
     end
 
     # Compute vacuum response matrix in-place (handles 2D single-n, 2D multi-n block-diagonal, and 3D)
-    vac_inputs = Vacuum.VacuumInput(equil, psilim, ctrl.mthvac, ctrl.nzvac, mpert, mlow, npert, nlow)
+    vac_inputs = Vacuum.VacuumInput(equil, psilim, ctrl.mthvac, ctrl.nzvac, mlow:mhigh, nlow:nhigh)
     Vacuum.compute_vacuum_response!(vac_data, vac_inputs, wall_settings)
 
     # Scale by (m - n*q)(m' - n'*q) [Chance Phys. Plasmas 1997 2161 eq. 126]
@@ -176,7 +176,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.mpert, intr.mlow, intr.npert, intr.nlow)
+        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)
         @views wv_array[i, :, :] .= wv
     end