From b9e22e8c3717fb8a19f82806457b62eb87f291c2 Mon Sep 17 00:00:00 2001 From: Jake Halpern Date: Fri, 10 Jul 2026 10:11:21 -0400 Subject: [PATCH 01/10] VACUUM - IMPROVEMENT - moving current matrix calculation into the vacuum module and making it a flag to compute it --- docs/src/vacuum.md | 2 +- docs/src/workflow.md | 5 +- src/ForceFreeStates/ForceFreeStatesStructs.jl | 8 - src/ForceFreeStates/Free.jl | 2 +- src/PerturbedEquilibrium/Response.jl | 10 +- src/PerturbedEquilibrium/SingularCoupling.jl | 90 ++---------- src/Vacuum/Vacuum.jl | 137 ++++++++++++------ test/runtests_vacuum.jl | 70 +++++---- 8 files changed, 150 insertions(+), 174 deletions(-) diff --git a/docs/src/vacuum.md b/docs/src/vacuum.md index a81e6d1e..a24266cf 100644 --- a/docs/src/vacuum.md +++ b/docs/src/vacuum.md @@ -75,7 +75,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 diff --git a/docs/src/workflow.md b/docs/src/workflow.md index d872dcea..556027bf 100644 --- a/docs/src/workflow.md +++ b/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) +- `current_matrix` — Vacuum surface-current matrix Iᵛ when `compute_L=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 current_matrix when computed for PE) | | `perturbed/` | Perturbed equilibrium: ξ, b in mode space, island diagnostics | diff --git a/src/ForceFreeStates/ForceFreeStatesStructs.jl b/src/ForceFreeStates/ForceFreeStatesStructs.jl index f03a7bbb..34b6205c 100644 --- a/src/ForceFreeStates/ForceFreeStatesStructs.jl +++ b/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 @@ -413,8 +409,6 @@ Populated in `Free.jl`. - `et::Vector{ComplexF64}` - Total eigenvalues of plasma + vacuum - `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 vacuum matrix wv, 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) """ @@ -444,8 +438,6 @@ Populated in `Free.jl`. rootA_wv::Array{ComplexF64,2} = fill(complex(NaN), numpert_total, numpert_total) rootA_wt::Array{ComplexF64,2} = fill(complex(NaN), numpert_total, 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 diff --git a/src/ForceFreeStates/Free.jl b/src/ForceFreeStates/Free.jl index 51860c30..163123b8 100644 --- a/src/ForceFreeStates/Free.jl +++ b/src/ForceFreeStates/Free.jl @@ -160,7 +160,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; force_wv_symmetry=ctrl.force_wv_symmetry) - 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 diff --git a/src/PerturbedEquilibrium/Response.jl b/src/PerturbedEquilibrium/Response.jl index ba282724..2e94512f 100644 --- a/src/PerturbedEquilibrium/Response.jl +++ b/src/PerturbedEquilibrium/Response.jl @@ -35,16 +35,12 @@ function compute_plasma_response!( # Plasma inductance Lambda (wt0 formula, Fortran resp_induct_flag=TRUE default) plasma_inductance = calc_plasma_inductance(vac_data, 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, 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) + _, current_matrix, _, _ = Vacuum.compute_vacuum_response(vac_input_2d, wall_nowall; compute_L=true) + surface_inductance = calc_surface_inductance(current_matrix) permeability = calc_permeability(plasma_inductance, surface_inductance) # Reluctance ϱ = L⁻¹·(Λ† − L)·L⁻¹ (Fortran gpresp_reluct: diff_indmats = CONJG(TRANSPOSE(plas_indmats)) − surf_indmats). diff --git a/src/PerturbedEquilibrium/SingularCoupling.jl b/src/PerturbedEquilibrium/SingularCoupling.jl index eb41cd76..51589f42 100644 --- a/src/PerturbedEquilibrium/SingularCoupling.jl +++ b/src/PerturbedEquilibrium/SingularCoupling.jl @@ -40,9 +40,6 @@ function _hermite_cubic_val(u_a, u_b, du_a, du_b, psi_a, psi_b, psi) return @. h00 * u_a + h * h10 * du_a + h01 * u_b + h * h11 * du_b end -# Reflect a periodic theta-space vector θ → -θ (the theta reversal in gpvacuum_flxsurf). -_reverse_theta(v::AbstractVector) = circshift(reverse(v), 1) - """ compute_singular_coupling_metrics!( state::PerturbedEquilibriumState, @@ -160,19 +157,12 @@ function compute_singular_coupling_metrics!( continue end - # Compute Green's functions at this surface for this n (once per pair) + # Surface-current matrix at this surface for this n (once per pair) vac_input = Vacuum.VacuumInput(equil, sing_surf.psifac, mtheta, 1, mlow:mhigh, [nn]) - _, grri_raw, grre_raw, _, _ = Vacuum.compute_vacuum_response(vac_input, wall_settings) - grri = Matrix{ComplexF64}(grri_raw) - grre = Matrix{ComplexF64}(grre_raw) - ffs_intr.sing[s].grri = grri - ffs_intr.sing[s].grre = grre - - # Get ν on the vacuum theta grid (same ν used in the vacuum Fourier basis computation) - ν_vac = Vacuum.PlasmaGeometry(vac_input).ν + _, current_matrix, _, _ = Vacuum.compute_vacuum_response(vac_input, wall_settings; compute_L=true) # 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) + L_surf = calc_surface_inductance(current_matrix) m_idx = m_res - mlow + 1 L_mm = L_surf[m_idx, m_idx] @@ -419,69 +409,17 @@ function compute_current_density( 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). + calc_surface_inductance(current_matrix::AbstractMatrix{ComplexF64})::Matrix{ComplexF64} -## Arguments +Invert the vacuum surface-current matrix into surface inductance `L`. - - `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] +Vacuum returns `current_matrix` from the Green’s→kax/DFT path (`gpvacuum_flxsurf`); +this helper applies regularization, inversion, and Hermitianization. """ -@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) +function calc_surface_inductance(current_matrix::AbstractMatrix{ComplexF64})::Matrix{ComplexF64} - 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) + mpert = size(current_matrix, 1) + μ₀ = 4π * 1e-7 L_surf = zeros(ComplexF64, mpert, mpert) current_mag = maximum(abs.(current_matrix)) @@ -493,11 +431,9 @@ Surface inductance matrix [mpert × mpert] end else try - regularization = 1e-12 * current_mag - current_reg = current_matrix + regularization * I - - L_surf = flux_matrix * inv(current_reg) - + current_reg = current_matrix + 1e-12 * current_mag * I + # L = flux * inv(current) with flux = I + L_surf = inv(current_reg) # Hermitianize (matches Fortran: temp1 = 0.5*(temp1 + CONJG(TRANSPOSE(temp1)))) L_surf = 0.5 * (L_surf + L_surf') catch e diff --git a/src/Vacuum/Vacuum.jl b/src/Vacuum/Vacuum.jl index f485b647..20db1eda 100644 --- a/src/Vacuum/Vacuum.jl +++ b/src/Vacuum/Vacuum.jl @@ -23,22 +23,64 @@ export compute_vacuum_response, compute_vacuum_response!, compute_vacuum_field export extract_plasma_surface_at_psi export PlasmaGeometry +const μ0 = 4π * 1e-7 + +# Vacuum code uses reverse-θ ordering relative to the equilibrium θ grid. +_reverse_theta(v::AbstractVector) = circshift(reverse(v), 1) + """ - _compute_vacuum_response_2d!(vac_data, inputs::VacuumInput, wall_settings::WallShapeSettings) + _surface_current_matrix_from_greens!(current_matrix, grri, grre, mlow, mpert, nn, ν) + +Build the Fourier surface-current matrix from interior/exterior Green's functions +(plasma-surface rows only). Implements the GPEC `gpvacuum_flxsurf` kax/DFT loop; +regularization and inversion to surface inductance `L` are left to PerturbedEquilibrium. +""" +@with_pool pool function _surface_current_matrix_from_greens!( + current_matrix::AbstractMatrix{ComplexF64}, + grri::AbstractMatrix{ComplexF64}, + grre::AbstractMatrix{ComplexF64}, + mlow::Int, + mpert::Int, + nn::Int, + ν::AbstractVector{Float64} +) + mtheta = length(ν) + ft = FourierTransform(mtheta, mpert, mlow) + + kax = zeros!(pool, ComplexF64, mtheta) + grri_surf = @view grri[1:mtheta, :] + grre_surf = @view grre[1:mtheta, :] + phase = cis.(-nn .* ν) + + for i in 1:mpert + # Complex grri/e stores exp(i(mθ-nν)) projection; conjugate for exp(-i(mθ-nν)). + # Eq. 10 of Park 2007 + kax .= conj.(grri_surf[:, i] .+ grre_surf[:, i]) ./ (μ0 * (2π)^2) + g_phased = kax .* phase + # Eq. 21b of Park 2007 + current_matrix[:, i] = ft(_reverse_theta(g_phased)) + end + return current_matrix +end + +""" + _compute_vacuum_response_2d!(vac_data, inputs::VacuumInput, wall_settings::WallShapeSettings; compute_L=false) 2D (axisymmetric) vacuum response calculation. Each toroidal mode `n` decouples in 2D geometry, so the routine loops over `inputs.n_modes`, -building the double-/single-layer operators, solving the exterior and interior systems, and -filling the corresponding diagonal block of the response matrix and the matching column block -of the Green's functions. +building the double-/single-layer operators, solving the exterior system for `wv`, and +optionally the interior system to build `current_matrix` when `compute_L=true`. +Green's functions are internal scratch only. """ -@with_pool pool function _compute_vacuum_response_2d!(vac_data, inputs::VacuumInput, wall_settings::WallShapeSettings) +@with_pool pool function _compute_vacuum_response_2d!(vac_data, inputs::VacuumInput, wall_settings::WallShapeSettings; compute_L::Bool=false) mpert = length(inputs.m_modes) + mlow = inputs.m_modes[1] num_points_surf = inputs.mtheta vac_data.wv .= 0 + compute_L && (vac_data.current_matrix .= 0) # Form the plasma and wall geometries plasma_surf = PlasmaGeometry(inputs) @@ -46,13 +88,11 @@ of the Green's functions. # Loop over all decoupled toroidal modes for (idx_n, n) in enumerate(inputs.n_modes) - ft = FourierTransform(inputs.mtheta, mpert, inputs.m_modes[1]; n=n, ν=plasma_surf.ν) + ft = FourierTransform(inputs.mtheta, mpert, mlow; n=n, ν=plasma_surf.ν) - # Diagonal block of wv and matching column block of the Green's functions + # Diagonal block of wv (and current_matrix when requested) block_idx = ((idx_n-1)*mpert+1):(idx_n*mpert) wv_block = @view vac_data.wv[block_idx, block_idx] - grri_block = @view vac_data.grri[:, block_idx] - grre_block = @view vac_data.grre[:, block_idx] # Active rows for computation (plasma only if no wall, plasma+wall if wall present) num_points_total = wall.nowall ? num_points_surf : 2 * num_points_surf @@ -60,10 +100,7 @@ of the Green's functions. # Local work matrices grad_green = zeros!(pool, num_points_total, num_points_total) green_temp = zeros!(pool, num_points_surf, num_points_surf) - - # Views into output Green's function matrices for the active rows/columns - grre = @view grre_block[1:num_points_total, :] - grri = @view grri_block[1:num_points_total, :] + grre = zeros!(pool, ComplexF64, num_points_total, mpert) # Plasma–Plasma block compute_2D_kernel_matrices!(grad_green, green_temp, plasma_surf, plasma_surf, n) @@ -82,20 +119,29 @@ of the Green's functions. mul!(view(grre, (num_points_surf+1):num_points_total, :), green_temp, ft.basis') end - # Compute both Green's functions: exterior (kernelsign=+1) then interior (kernelsign=-1) - grri .= grre # start from same as exterior - grad_green_interior = similar!(pool, grad_green) - grad_green_interior .= grad_green + if compute_L + # Copy RHS before exterior solve overwrites grre; keep a kernel copy for interior + grri = similar!(pool, grre) + grri .= grre + grad_green_interior = similar!(pool, grad_green) + grad_green_interior .= grad_green + + # Solve exterior first, overwriting grad_green to save memory since we already have the interior kernel + ldiv!(lu!(grad_green), grre) - # Solve exterior first, overwriting grad_green to save memory since we already have the interior kernel - ldiv!(lu!(grad_green), grre) + # Interior flips the normal sign but not the diagonal: -K + 2I + grad_green_interior .*= -1 + for i in 1:num_points_total + grad_green_interior[i, i] += 2.0 + end + ldiv!(lu!(grad_green_interior), grri) - # Interior flips the sign of the normal, but not the diagonal terms, so we multiply by -1 and add 2I to the diagonal - grad_green_interior .*= -1 - for i in 1:num_points_total - grad_green_interior[i, i] += 2.0 + cm_block = @view vac_data.current_matrix[block_idx, block_idx] + _surface_current_matrix_from_greens!(cm_block, grri, grre, mlow, mpert, n, plasma_surf.ν) + else + # Only need exterior system for wv + ldiv!(lu!(grad_green), grre) end - ldiv!(lu!(grad_green_interior), grri) # Project exterior kernel onto observer basis exp(-i*(mθ - nν)) and scale to get the response matrix mul!(wv_block, ft.basis, @view(grre[1:num_points_surf, :])) @@ -116,7 +162,7 @@ of the Green's functions. end """ - _compute_vacuum_response_3d!(vac_data, inputs::VacuumInput, wall_settings::WallShapeSettings) + _compute_vacuum_response_3d!(vac_data, inputs::VacuumInput, wall_settings::WallShapeSettings; compute_L=false) 3D (`inputs.nzeta > 1`) vacuum response via block-circulant field-period reduction. For `nfp == 1` the block-circulant assembly and residue-class loop are skipped in favour of a more efficient @@ -134,14 +180,19 @@ class needs one solve `wv[class k] = (4π²/M)·E_localᴴ·(D̂ₖ \\ Ŝₖ)|_p block-row of the operators is built, so the kernel cost drops by `nfp` and the dense `O(N³)` factorization is replaced by per-class `O(M³)` solves (`M = N/nfp`). -Only `wv` is produced currently; `grri`/`grre` are returned zeroed and are not yet supported in the 3D path. +Only `wv` is produced currently; `current_matrix` is left zeroed when `compute_L=true` +(surface-current / inductance not yet supported in 3D). Extension point: per residue class, apply the per-period basis to `D̂ₖ⁻¹Ŝₖ` for the exterior columns and the interior variant `-D + 2I` for the interior columns, then scatter back into the `[2N × 2·num_modes]` arrays. """ -@with_pool pool function _compute_vacuum_response_3d!(vac_data, inputs::VacuumInput, wall_settings::WallShapeSettings) +@with_pool pool function _compute_vacuum_response_3d!(vac_data, inputs::VacuumInput, wall_settings::WallShapeSettings; compute_L::Bool=false) (; mtheta, nzeta, nfp, m_modes, n_modes) = inputs fill!(vac_data.wv, 0) + fill!(vac_data.current_matrix, 0) + if compute_L + @warn "compute_L=true is not supported for 3D vacuum response; current_matrix left as zeros" maxlog=1 + end # Full-torus geometry for source surface; observers are restricted to one field period full = expand_field_periods(inputs) @@ -216,20 +267,16 @@ interior variant `-D + 2I` for the interior columns, then scatter back into the inputs.force_wv_symmetry && hermitianpart!(vac_data.wv) - # Zero out the Green's function matrices (not tested in 3D yet) - fill!(vac_data.grri, 0) - fill!(vac_data.grre, 0) - # Populate coordinate arrays vac_data.plasma_pts .= plasma_surf.r vac_data.wall_pts .= wall.r end """ - compute_vacuum_response(inputs::VacuumInput, wall_settings::WallShapeSettings) + compute_vacuum_response(inputs::VacuumInput, wall_settings::WallShapeSettings; compute_L=false) -Allocate and return the vacuum response matrix and Green's functions for the given vacuum -inputs. Thin allocating wrapper around the in-place [`compute_vacuum_response!`]: it sizes the +Allocate and return the vacuum response matrix and optional surface-current matrix for the given +vacuum inputs. Thin allocating wrapper around the in-place [`compute_vacuum_response!`]: it sizes the output arrays for the full torus and forwards to the same 2D/3D workers. For performance-critical paths that already own preallocated storage (e.g. `ForceFreeStates.VacuumData`), prefer the in-place method to avoid extra heap allocations. @@ -237,28 +284,27 @@ in-place method to avoid extra heap allocations. # Returns - `wv`: complex vacuum response matrix (`num_modes × num_modes`). - - `grri`, `grre`: interior/exterior Green's functions (zeroed on the 3D nowall path). + - `current_matrix`: Fourier surface-current matrix (`num_modes × num_modes`); zeros unless `compute_L=true`. - `plasma_pts`, `wall_pts`: surface coordinate arrays. """ -function compute_vacuum_response(inputs::VacuumInput, wall_settings::WallShapeSettings) +function compute_vacuum_response(inputs::VacuumInput, wall_settings::WallShapeSettings; compute_L::Bool=false) num_points = inputs.mtheta * inputs.nzeta * inputs.nfp # mtheta for 2D num_modes = length(inputs.m_modes) * length(inputs.n_modes) vac = ( wv=zeros(ComplexF64, num_modes, num_modes), - grri=zeros(ComplexF64, 2 * num_points, num_modes), - grre=zeros(ComplexF64, 2 * num_points, num_modes), + current_matrix=zeros(ComplexF64, num_modes, num_modes), plasma_pts=zeros(num_points, 3), wall_pts=zeros(num_points, 3) ) - compute_vacuum_response!(vac, inputs, wall_settings) + compute_vacuum_response!(vac, inputs, wall_settings; compute_L) - return vac.wv, vac.grri, vac.grre, vac.plasma_pts, vac.wall_pts + return vac.wv, vac.current_matrix, vac.plasma_pts, vac.wall_pts end """ - compute_vacuum_response!(vac_data, inputs::VacuumInput, wall_settings::WallShapeSettings) + compute_vacuum_response!(vac_data, inputs::VacuumInput, wall_settings::WallShapeSettings; compute_L=false) In-place variant that computes the vacuum response and directly populates the arrays stored in `vac_data`. Dispatches on dimensionality only: 2D (`inputs.nzeta == 1`) routes to @@ -268,19 +314,18 @@ The `vac_data` argument is expected to provide the following writable fields wit sizes: - `wv::AbstractMatrix{ComplexF64}` – vacuum response matrix - - `grri::AbstractMatrix{ComplexF64}` – interior Green's functions - - `grre::AbstractMatrix{ComplexF64}` – exterior Green's functions - `plasma_pts::AbstractMatrix{Float64}` – plasma surface coordinates - `wall_pts::AbstractMatrix{Float64}` – wall surface coordinates + - `current_matrix::AbstractMatrix{ComplexF64}` – required when `compute_L=true` This is designed to work with `ForceFreeStates.VacuumData` but does not depend on its concrete type (duck-typed on field names only). """ -function compute_vacuum_response!(vac_data, inputs::VacuumInput, wall_settings::WallShapeSettings) +function compute_vacuum_response!(vac_data, inputs::VacuumInput, wall_settings::WallShapeSettings; compute_L::Bool=false) if inputs.nzeta == 1 - _compute_vacuum_response_2d!(vac_data, inputs, wall_settings) + _compute_vacuum_response_2d!(vac_data, inputs, wall_settings; compute_L) else - _compute_vacuum_response_3d!(vac_data, inputs, wall_settings) + _compute_vacuum_response_3d!(vac_data, inputs, wall_settings; compute_L) end end diff --git a/test/runtests_vacuum.jl b/test/runtests_vacuum.jl index 52a73494..035ed96b 100644 --- a/test/runtests_vacuum.jl +++ b/test/runtests_vacuum.jl @@ -369,7 +369,7 @@ @testset "nowall" begin inputs = _make_inputs() wall_settings = WallShapeSettings(shape="nowall") - wv, grri, grre, plasma_pts, wall_pts = compute_vacuum_response(inputs, wall_settings) + wv, current_matrix, plasma_pts, wall_pts = compute_vacuum_response(inputs, wall_settings) numpoints = inputs.mtheta * inputs.nzeta num_modes = length(inputs.m_modes) * length(inputs.n_modes) @@ -377,10 +377,8 @@ @test size(wv) == (num_modes, num_modes) @test eltype(wv) == ComplexF64 @test all(isfinite, wv) - @test size(grri) == (2 * numpoints, num_modes) - @test size(grre) == (2 * numpoints, num_modes) - @test all(isfinite, grri) - @test all(isfinite, grre) + @test size(current_matrix) == (num_modes, num_modes) + @test all(iszero, current_matrix) # compute_L=false by default @test size(plasma_pts) == (numpoints, 3) @test all(isfinite, plasma_pts) @test size(wall_pts) == (numpoints, 3) @@ -389,15 +387,29 @@ @test isapprox(wv, wv', rtol=1e-12) end + @testset "nowall compute_L=true" begin + inputs = _make_inputs() + wall_settings = WallShapeSettings(shape="nowall") + wv, current_matrix, plasma_pts, wall_pts = compute_vacuum_response(inputs, wall_settings; compute_L=true) + + num_modes = length(inputs.m_modes) * length(inputs.n_modes) + @test size(wv) == (num_modes, num_modes) + @test all(isfinite, wv) + @test size(current_matrix) == (num_modes, num_modes) + @test all(isfinite, current_matrix) + @test !all(iszero, current_matrix) + @test isapprox(wv, wv', rtol=1e-12) + end + @testset "conformal wall" begin inputs = _make_inputs() wall_settings = WallShapeSettings(shape="conformal", a=0.5) - wv, grri, grre, plasma_pts, wall_pts = compute_vacuum_response(inputs, wall_settings) + wv, current_matrix, plasma_pts, wall_pts = compute_vacuum_response(inputs, wall_settings) numpoints = inputs.mtheta * inputs.nzeta num_modes = length(inputs.m_modes) * length(inputs.n_modes) @test size(wv) == (num_modes, num_modes) - @test size(grri) == (2 * numpoints, num_modes) + @test all(iszero, current_matrix) @test all(isfinite, plasma_pts) @test all(isfinite, wall_pts) # plasma_pts layout: col1=R, col2=0, col3=Z @@ -409,19 +421,19 @@ @testset "edge: single poloidal mode mpert=1" begin inputs = _make_inputs(m_modes=[1], n_modes=[1]) wall_settings = WallShapeSettings(shape="nowall") - wv, grri, grre, plasma_pts, wall_pts = compute_vacuum_response(inputs, wall_settings) + wv, current_matrix, plasma_pts, wall_pts = compute_vacuum_response(inputs, wall_settings) @test size(wv) == (1, 1) @test all(isfinite, wv) - @test size(grri, 2) == 1 + @test size(current_matrix) == (1, 1) end @testset "edge: small mtheta" begin # Keep mtheta_eq=17 so boundary has enough points for periodic spline inputs = _make_inputs(mtheta=16, mtheta_eq=17) wall_settings = WallShapeSettings(shape="nowall") - wv, grri, grre, plasma_pts, wall_pts = compute_vacuum_response(inputs, wall_settings) + wv, current_matrix, plasma_pts, wall_pts = compute_vacuum_response(inputs, wall_settings) @test size(wv) == (2, 2) - @test size(grri) == (32, 2) # 2*mtheta, num_modes=2 + @test size(current_matrix) == (2, 2) @test size(plasma_pts) == (16, 3) end @@ -430,17 +442,17 @@ # in-place entry populates caller-owned duck-typed (NamedTuple) storage identically. for wall_settings in (WallShapeSettings(shape="nowall"), WallShapeSettings(shape="conformal", a=0.5)) inputs = _make_inputs() - wv, grri, grre, pp, wp = compute_vacuum_response(inputs, wall_settings) + wv, current_matrix, pp, wp = compute_vacuum_response(inputs, wall_settings; compute_L=true) numpoints = inputs.mtheta * inputs.nzeta num_modes = length(inputs.m_modes) * length(inputs.n_modes) - vac = (wv=zeros(ComplexF64, num_modes, num_modes), grri=zeros(ComplexF64, 2 * numpoints, num_modes), grre=zeros(ComplexF64, 2 * numpoints, num_modes), + vac = (wv=zeros(ComplexF64, num_modes, num_modes), + current_matrix=zeros(ComplexF64, num_modes, num_modes), plasma_pts=zeros(numpoints, 3), wall_pts=zeros(numpoints, 3)) - compute_vacuum_response!(vac, inputs, wall_settings) + compute_vacuum_response!(vac, inputs, wall_settings; compute_L=true) @test vac.wv ≈ wv - @test vac.grri ≈ grri - @test vac.grre ≈ grre + @test vac.current_matrix ≈ current_matrix @test vac.plasma_pts ≈ pp @test vac.wall_pts ≈ wp end @@ -582,17 +594,15 @@ @testset "compute_vacuum_response 3D nowall" begin inputs = _make_3d_inputs(mtheta=32, nzeta=32, mtheta_eq=17) wall_settings = WallShapeSettings(shape="nowall") - wv, grri, grre, plasma_pts, wall_pts = compute_vacuum_response(inputs, wall_settings) + wv, current_matrix, plasma_pts, wall_pts = compute_vacuum_response(inputs, wall_settings) numpoints = inputs.mtheta * inputs.nzeta num_modes = length(inputs.m_modes) * length(inputs.n_modes) @test size(wv) == (num_modes, num_modes) @test eltype(wv) == ComplexF64 @test all(isfinite, wv) - @test size(grri) == (2 * numpoints, num_modes) - @test size(grre) == (2 * numpoints, num_modes) - @test all(isfinite, grri) - @test all(isfinite, grre) + @test size(current_matrix) == (num_modes, num_modes) + @test all(iszero, current_matrix) @test size(plasma_pts) == (numpoints, 3) @test all(isfinite, plasma_pts) @test size(wall_pts) == (numpoints, 3) @@ -605,7 +615,7 @@ @testset "compute_vacuum_response 3D nonaxisymmetric boundary" begin inputs = _make_3d_nonaxis_inputs(mtheta=24, nzeta=24, mtheta_in=12, nzeta_in=12, mpert=2, nlow=0, npert=2) wall_settings = WallShapeSettings(shape="nowall") - wv, grri, grre, plasma_pts, wall_pts = compute_vacuum_response(inputs, wall_settings) + wv, current_matrix, plasma_pts, wall_pts = compute_vacuum_response(inputs, wall_settings) numpoints = inputs.mtheta * inputs.nzeta num_modes = length(inputs.m_modes) * length(inputs.n_modes) @@ -613,10 +623,8 @@ @test size(wv) == (num_modes, num_modes) @test eltype(wv) == ComplexF64 @test all(isfinite, wv) - @test size(grri) == (2 * numpoints, num_modes) - @test size(grre) == (2 * numpoints, num_modes) - @test all(isfinite, grri) - @test all(isfinite, grre) + @test size(current_matrix) == (num_modes, num_modes) + @test all(iszero, current_matrix) @test size(plasma_pts) == (numpoints, 3) @test all(isfinite, plasma_pts) @test size(wall_pts) == (numpoints, 3) @@ -626,12 +634,12 @@ @testset "compute_vacuum_response 3D conformal wall" begin inputs = _make_3d_inputs(mtheta=32, nzeta=32, mtheta_eq=17) wall_settings = WallShapeSettings(shape="conformal", a=0.3) - wv, grri, grre, plasma_pts, wall_pts = compute_vacuum_response(inputs, wall_settings) + wv, current_matrix, plasma_pts, wall_pts = compute_vacuum_response(inputs, wall_settings) numpoints = inputs.mtheta * inputs.nzeta num_modes = length(inputs.m_modes) * length(inputs.n_modes) @test size(wv) == (num_modes, num_modes) - @test size(grri) == (2 * numpoints, num_modes) + @test size(current_matrix) == (num_modes, num_modes) @test all(isfinite, plasma_pts) @test all(isfinite, wall_pts) # Wall and plasma should differ (conformal wall offset from plasma) @@ -677,13 +685,13 @@ wall_settings = WallShapeSettings(shape="nowall") # Reduced (block-circulant) path - wv_red, _, _, plasma_pts_red, _ = compute_vacuum_response(inputs_red, wall_settings) + wv_red, _, plasma_pts_red, _ = compute_vacuum_response(inputs_red, wall_settings) # Full-torus reference: pre-expand so nfp=1 forces the dense path inputs_full = GeneralizedPerturbedEquilibrium.Vacuum.expand_field_periods(inputs_red) @test inputs_full.nfp == 1 @test inputs_full.nzeta == nzeta_p * nfp - wv_full, _, _, _, _ = compute_vacuum_response(inputs_full, wall_settings) + wv_full, _, _, _ = compute_vacuum_response(inputs_full, wall_settings) num_modes = length(m_modes) * length(n_modes) @test size(wv_red) == (num_modes, num_modes) @@ -714,7 +722,7 @@ mtheta=mtheta, nzeta=nzeta_p, force_wv_symmetry=true, nfp=nfp ) - wv_sym, _, _, _, _ = compute_vacuum_response(inputs_sym, wall_settings) + wv_sym, _, _, _ = compute_vacuum_response(inputs_sym, wall_settings) @test isapprox(wv_sym, wv_sym', rtol=1e-12) end From 82923723b0bdc945518e9f28aa35f319825c12b8 Mon Sep 17 00:00:00 2001 From: Jake Halpern Date: Fri, 10 Jul 2026 12:56:34 -0400 Subject: [PATCH 02/10] VACUUM - MINOR - using mul! type basis calculation for surface current matrix --- src/PerturbedEquilibrium/SingularCoupling.jl | 7 +-- src/Vacuum/Vacuum.jl | 47 +++----------------- 2 files changed, 10 insertions(+), 44 deletions(-) diff --git a/src/PerturbedEquilibrium/SingularCoupling.jl b/src/PerturbedEquilibrium/SingularCoupling.jl index 51589f42..9d81cb9e 100644 --- a/src/PerturbedEquilibrium/SingularCoupling.jl +++ b/src/PerturbedEquilibrium/SingularCoupling.jl @@ -414,7 +414,7 @@ end Invert the vacuum surface-current matrix into surface inductance `L`. Vacuum returns `current_matrix` from the Green’s→kax/DFT path (`gpvacuum_flxsurf`); -this helper applies regularization, inversion, and Hermitianization. +this helper applies normalization, regularization, inversion, and Hermitianization. """ function calc_surface_inductance(current_matrix::AbstractMatrix{ComplexF64})::Matrix{ComplexF64} @@ -431,9 +431,10 @@ function calc_surface_inductance(current_matrix::AbstractMatrix{ComplexF64})::Ma end else try - current_reg = current_matrix + 1e-12 * current_mag * I + current_matrix ./= μ₀ * (2π)^2 + current_matrix += 1e-12 * current_mag * I # L = flux * inv(current) with flux = I - L_surf = inv(current_reg) + L_surf = inv(current_matrix) # Hermitianize (matches Fortran: temp1 = 0.5*(temp1 + CONJG(TRANSPOSE(temp1)))) L_surf = 0.5 * (L_surf + L_surf') catch e diff --git a/src/Vacuum/Vacuum.jl b/src/Vacuum/Vacuum.jl index 20db1eda..70f7de39 100644 --- a/src/Vacuum/Vacuum.jl +++ b/src/Vacuum/Vacuum.jl @@ -23,46 +23,6 @@ export compute_vacuum_response, compute_vacuum_response!, compute_vacuum_field export extract_plasma_surface_at_psi export PlasmaGeometry -const μ0 = 4π * 1e-7 - -# Vacuum code uses reverse-θ ordering relative to the equilibrium θ grid. -_reverse_theta(v::AbstractVector) = circshift(reverse(v), 1) - -""" - _surface_current_matrix_from_greens!(current_matrix, grri, grre, mlow, mpert, nn, ν) - -Build the Fourier surface-current matrix from interior/exterior Green's functions -(plasma-surface rows only). Implements the GPEC `gpvacuum_flxsurf` kax/DFT loop; -regularization and inversion to surface inductance `L` are left to PerturbedEquilibrium. -""" -@with_pool pool function _surface_current_matrix_from_greens!( - current_matrix::AbstractMatrix{ComplexF64}, - grri::AbstractMatrix{ComplexF64}, - grre::AbstractMatrix{ComplexF64}, - mlow::Int, - mpert::Int, - nn::Int, - ν::AbstractVector{Float64} -) - mtheta = length(ν) - ft = FourierTransform(mtheta, mpert, mlow) - - kax = zeros!(pool, ComplexF64, mtheta) - grri_surf = @view grri[1:mtheta, :] - grre_surf = @view grre[1:mtheta, :] - phase = cis.(-nn .* ν) - - for i in 1:mpert - # Complex grri/e stores exp(i(mθ-nν)) projection; conjugate for exp(-i(mθ-nν)). - # Eq. 10 of Park 2007 - kax .= conj.(grri_surf[:, i] .+ grre_surf[:, i]) ./ (μ0 * (2π)^2) - g_phased = kax .* phase - # Eq. 21b of Park 2007 - current_matrix[:, i] = ft(_reverse_theta(g_phased)) - end - return current_matrix -end - """ _compute_vacuum_response_2d!(vac_data, inputs::VacuumInput, wall_settings::WallShapeSettings; compute_L=false) @@ -136,8 +96,13 @@ Green's functions are internal scratch only. end ldiv!(lu!(grad_green_interior), grri) + # Build surface-current matrix from interior/exterior Green's functions + # Eq. 10 of Park 2007: J_surface ∝ conj(grri + grre) + # ft.basis = exp(-i(mθ - nν)); conjugate to get exp(-i(mθ + nν)) to match original phase convention cm_block = @view vac_data.current_matrix[block_idx, block_idx] - _surface_current_matrix_from_greens!(cm_block, grri, grre, mlow, mpert, n, plasma_surf.ν) + @views g_sum = conj.(grri[1:num_points_surf, :] .+ grre[1:num_points_surf, :]) + mul!(cm_block, conj.(ft.basis), g_sum) + cm_block ./= num_points_surf else # Only need exterior system for wv ldiv!(lu!(grad_green), grre) From 00e1259b9684311377a051fb0179fb22b5e636dd Mon Sep 17 00:00:00 2001 From: Jake Halpern Date: Fri, 10 Jul 2026 15:43:12 -0400 Subject: [PATCH 03/10] VACUUM - MINOR - cleaning up I_v calculation and adding more detailed comments from Claude --- docs/src/workflow.md | 4 +- src/PerturbedEquilibrium/Response.jl | 4 +- src/PerturbedEquilibrium/SingularCoupling.jl | 24 ++++----- src/Vacuum/DataTypes.jl | 3 +- src/Vacuum/Kernel2D.jl | 4 +- src/Vacuum/Vacuum.jl | 51 +++++++++++--------- test/runtests_vacuum.jl | 48 +++++++++--------- 7 files changed, 74 insertions(+), 64 deletions(-) diff --git a/docs/src/workflow.md b/docs/src/workflow.md index 556027bf..37d79deb 100644 --- a/docs/src/workflow.md +++ b/docs/src/workflow.md @@ -59,7 +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) -- `current_matrix` — Vacuum surface-current matrix Iᵛ when `compute_L=true` (otherwise zeros); PerturbedEquilibrium inverts this to surface inductance `L` +- `I_v` — Vacuum surface-current matrix Iᵛ when `compute_L=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) @@ -181,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 (and current_matrix when computed for PE) | +| `vacuum/` | Vacuum response matrices: wv (and I_v when computed for PE) | | `perturbed/` | Perturbed equilibrium: ξ, b in mode space, island diagnostics | diff --git a/src/PerturbedEquilibrium/Response.jl b/src/PerturbedEquilibrium/Response.jl index 2e94512f..bc179449 100644 --- a/src/PerturbedEquilibrium/Response.jl +++ b/src/PerturbedEquilibrium/Response.jl @@ -39,8 +39,8 @@ function compute_plasma_response!( 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") - _, current_matrix, _, _ = Vacuum.compute_vacuum_response(vac_input_2d, wall_nowall; compute_L=true) - surface_inductance = calc_surface_inductance(current_matrix) + _, I_v, _, _ = Vacuum.compute_vacuum_response(vac_input_2d, wall_nowall; compute_L=true) + surface_inductance = calc_surface_inductance(I_v) permeability = calc_permeability(plasma_inductance, surface_inductance) # Reluctance ϱ = L⁻¹·(Λ† − L)·L⁻¹ (Fortran gpresp_reluct: diff_indmats = CONJG(TRANSPOSE(plas_indmats)) − surf_indmats). diff --git a/src/PerturbedEquilibrium/SingularCoupling.jl b/src/PerturbedEquilibrium/SingularCoupling.jl index 9d81cb9e..483635c7 100644 --- a/src/PerturbedEquilibrium/SingularCoupling.jl +++ b/src/PerturbedEquilibrium/SingularCoupling.jl @@ -159,10 +159,10 @@ function compute_singular_coupling_metrics!( # Surface-current matrix at this surface for this n (once per pair) vac_input = Vacuum.VacuumInput(equil, sing_surf.psifac, mtheta, 1, mlow:mhigh, [nn]) - _, current_matrix, _, _ = Vacuum.compute_vacuum_response(vac_input, wall_settings; compute_L=true) + _, I_v, _, _ = Vacuum.compute_vacuum_response(vac_input, wall_settings; compute_L=true) # Precompute L_surf; only the (m_res, m_res) diagonal element is needed for singflx - L_surf = calc_surface_inductance(current_matrix) + L_surf = calc_surface_inductance(I_v) m_idx = m_res - mlow + 1 L_mm = L_surf[m_idx, m_idx] @@ -409,20 +409,20 @@ function compute_current_density( end """ - calc_surface_inductance(current_matrix::AbstractMatrix{ComplexF64})::Matrix{ComplexF64} + calc_surface_inductance(I_v::AbstractMatrix{ComplexF64})::Matrix{ComplexF64} Invert the vacuum surface-current matrix into surface inductance `L`. -Vacuum returns `current_matrix` from the Green’s→kax/DFT path (`gpvacuum_flxsurf`); -this helper applies normalization, regularization, inversion, and Hermitianization. +Vacuum returns `I_v` from the Green’s functions, this helper applies GPEC normalization, +regularization, inversion, and Hermitianization. """ -function calc_surface_inductance(current_matrix::AbstractMatrix{ComplexF64})::Matrix{ComplexF64} +function calc_surface_inductance(I_v::AbstractMatrix{ComplexF64})::Matrix{ComplexF64} - mpert = size(current_matrix, 1) + mpert = size(I_v, 1) μ₀ = 4π * 1e-7 L_surf = zeros(ComplexF64, mpert, mpert) - current_mag = maximum(abs.(current_matrix)) + current_mag = maximum(abs.(I_v)) if current_mag < 1e-15 @warn "Current matrix is all zeros! Cannot compute surface inductance." maxlog=1 @@ -431,10 +431,10 @@ function calc_surface_inductance(current_matrix::AbstractMatrix{ComplexF64})::Ma end else try - current_matrix ./= μ₀ * (2π)^2 - current_matrix += 1e-12 * current_mag * I - # L = flux * inv(current) with flux = I - L_surf = inv(current_matrix) + I_v ./= μ₀ * (2π)^2 + I_v += 1e-12 * current_mag * I + # L = flux * inv(Iᵛ) with flux = I + L_surf = inv(I_v) # Hermitianize (matches Fortran: temp1 = 0.5*(temp1 + CONJG(TRANSPOSE(temp1)))) L_surf = 0.5 * (L_surf + L_surf') catch e diff --git a/src/Vacuum/DataTypes.jl b/src/Vacuum/DataTypes.jl index 4ba746ba..629e3450 100644 --- a/src/Vacuum/DataTypes.jl +++ b/src/Vacuum/DataTypes.jl @@ -78,7 +78,8 @@ function VacuumInput( r, z, ν = extract_plasma_surface_at_psi(equil, ψ) # Remove the last point to go from the [0, 2π] grid to VACUUM's [0, 2π) grid - # and reverse the arrays for VACUUM's CW θ direction + # and reverse the arrays for VACUUM's CW θ direction (θ_VAC = -θ_GPEC). This handedness is why + # operators returned to GPEC (e.g. the surface-inductance current matrix) are conjugated. return VacuumInput(; x=reverse(r)[1:(end-1)], z=reverse(z)[1:(end-1)], diff --git a/src/Vacuum/Kernel2D.jl b/src/Vacuum/Kernel2D.jl index ea93f224..a5b03397 100644 --- a/src/Vacuum/Kernel2D.jl +++ b/src/Vacuum/Kernel2D.jl @@ -215,7 +215,9 @@ grad_greenfunction is not zeroed since it fills a different block of the end # Normals need to point outward from vacuum region. In VACUUM clockwise θ convention, normal points - # out of vacuum for wall but inward for plasma, so we multiply by -1 for plasma sources + # out of vacuum for wall but inward for plasma, so we multiply by -1 for plasma sources. + # This normal orientation sets the sign of the surface current K = n̂×∇χ (Chance 1997 eq. 184), + # and is the origin of the leading minus in the current-matrix build (Vacuum.jl, compute_L). if source isa PlasmaGeometry grad_greenfunction_block .*= -1 end diff --git a/src/Vacuum/Vacuum.jl b/src/Vacuum/Vacuum.jl index 70f7de39..5dc32de7 100644 --- a/src/Vacuum/Vacuum.jl +++ b/src/Vacuum/Vacuum.jl @@ -30,7 +30,7 @@ export PlasmaGeometry Each toroidal mode `n` decouples in 2D geometry, so the routine loops over `inputs.n_modes`, building the double-/single-layer operators, solving the exterior system for `wv`, and -optionally the interior system to build `current_matrix` when `compute_L=true`. +optionally the interior system to build `I_v` when `compute_L=true`. Green's functions are internal scratch only. """ @with_pool pool function _compute_vacuum_response_2d!(vac_data, inputs::VacuumInput, wall_settings::WallShapeSettings; compute_L::Bool=false) @@ -40,7 +40,7 @@ Green's functions are internal scratch only. num_points_surf = inputs.mtheta vac_data.wv .= 0 - compute_L && (vac_data.current_matrix .= 0) + compute_L && (vac_data.I_v .= 0) # Form the plasma and wall geometries plasma_surf = PlasmaGeometry(inputs) @@ -50,7 +50,7 @@ Green's functions are internal scratch only. for (idx_n, n) in enumerate(inputs.n_modes) ft = FourierTransform(inputs.mtheta, mpert, mlow; n=n, ν=plasma_surf.ν) - # Diagonal block of wv (and current_matrix when requested) + # Diagonal block of wv (and I_v when requested) block_idx = ((idx_n-1)*mpert+1):(idx_n*mpert) wv_block = @view vac_data.wv[block_idx, block_idx] @@ -86,29 +86,36 @@ Green's functions are internal scratch only. grad_green_interior = similar!(pool, grad_green) grad_green_interior .= grad_green - # Solve exterior first, overwriting grad_green to save memory since we already have the interior kernel + # Exterior operator D_ext = 2I + 𝒦 (Chance 1997 eq. 89); solve for the physical + # vacuum-outside potential grre = χ^(vo). Overwrites grad_green to save memory. ldiv!(lu!(grad_green), grre) - # Interior flips the normal sign but not the diagonal: -K + 2I - grad_green_interior .*= -1 + # Interior operator D_int = D_ext - 2I: the double-layer jump between the two one-sided + # boundary limits is 2I here, giving the vacuum-inside potential grri = χ^(vi). for i in 1:num_points_total - grad_green_interior[i, i] += 2.0 + grad_green_interior[i, i] -= 2.0 end ldiv!(lu!(grad_green_interior), grri) - # Build surface-current matrix from interior/exterior Green's functions - # Eq. 10 of Park 2007: J_surface ∝ conj(grri + grre) - # ft.basis = exp(-i(mθ - nν)); conjugate to get exp(-i(mθ + nν)) to match original phase convention - cm_block = @view vac_data.current_matrix[block_idx, block_idx] - @views g_sum = conj.(grri[1:num_points_surf, :] .+ grre[1:num_points_surf, :]) - mul!(cm_block, conj.(ft.basis), g_sum) - cm_block ./= num_points_surf + # Surface-current matrix, Park 2007 eq. 21b: μ₀I^v = χ^(vi) - χ^(vo) = grri - grre + # The minus is because VACUUM builds the operators in its CW-θ frame while GPEC uses CCW-θ, + # flipping the outward-normal sign. + I_v_block = @view vac_data.I_v[block_idx, block_idx] + @views g_sum = -(grri[1:num_points_surf, :] .- grre[1:num_points_surf, :]) + mul!(I_v_block, ft.basis, g_sum) + I_v_block ./= num_points_surf + + # Iᵛ becomes the surface inductance L (Φ = L·Iᵛ), a one-sided operator that + # multiplies GPEC-convention mode vectors. Kernels are real, so conj rotates B·A·B† into + # the GPEC (CCW-θ, conjugate-of-RH) basis. Unlike wv this is not an invariant scalar. + I_v_block .= conj.(I_v_block) else # Only need exterior system for wv ldiv!(lu!(grad_green), grre) end - # Project exterior kernel onto observer basis exp(-i*(mθ - nν)) and scale to get the response matrix + # wv = B·(D_ext⁻¹𝒢)·B† is a Hermitian energy generator (Chance 1997 eqs. 125-126); it is + # contracted on both sides to the real, convention-invariant vacuum energy, so no conjugation. mul!(wv_block, ft.basis, @view(grre[1:num_points_surf, :])) wv_block .*= 4π^2 / num_points_surf end @@ -145,7 +152,7 @@ class needs one solve `wv[class k] = (4π²/M)·E_localᴴ·(D̂ₖ \\ Ŝₖ)|_p block-row of the operators is built, so the kernel cost drops by `nfp` and the dense `O(N³)` factorization is replaced by per-class `O(M³)` solves (`M = N/nfp`). -Only `wv` is produced currently; `current_matrix` is left zeroed when `compute_L=true` +Only `wv` is produced currently; `I_v` is left zeroed when `compute_L=true` (surface-current / inductance not yet supported in 3D). Extension point: per residue class, apply the per-period basis to `D̂ₖ⁻¹Ŝₖ` for the exterior columns and the interior variant `-D + 2I` for the interior columns, then scatter back into the `[2N × 2·num_modes]` arrays. @@ -154,9 +161,9 @@ interior variant `-D + 2I` for the interior columns, then scatter back into the (; mtheta, nzeta, nfp, m_modes, n_modes) = inputs fill!(vac_data.wv, 0) - fill!(vac_data.current_matrix, 0) + fill!(vac_data.I_v, 0) if compute_L - @warn "compute_L=true is not supported for 3D vacuum response; current_matrix left as zeros" maxlog=1 + @warn "compute_L=true is not supported for 3D vacuum response; I_v left as zeros" maxlog=1 end # Full-torus geometry for source surface; observers are restricted to one field period @@ -249,7 +256,7 @@ in-place method to avoid extra heap allocations. # Returns - `wv`: complex vacuum response matrix (`num_modes × num_modes`). - - `current_matrix`: Fourier surface-current matrix (`num_modes × num_modes`); zeros unless `compute_L=true`. + - `I_v`: Vacuum surface-current matrix (`num_modes × num_modes`); zeros unless `compute_L=true`. - `plasma_pts`, `wall_pts`: surface coordinate arrays. """ function compute_vacuum_response(inputs::VacuumInput, wall_settings::WallShapeSettings; compute_L::Bool=false) @@ -259,13 +266,13 @@ function compute_vacuum_response(inputs::VacuumInput, wall_settings::WallShapeSe vac = ( wv=zeros(ComplexF64, num_modes, num_modes), - current_matrix=zeros(ComplexF64, num_modes, num_modes), + I_v=zeros(ComplexF64, num_modes, num_modes), plasma_pts=zeros(num_points, 3), wall_pts=zeros(num_points, 3) ) compute_vacuum_response!(vac, inputs, wall_settings; compute_L) - return vac.wv, vac.current_matrix, vac.plasma_pts, vac.wall_pts + return vac.wv, vac.I_v, vac.plasma_pts, vac.wall_pts end """ @@ -281,7 +288,7 @@ sizes: - `wv::AbstractMatrix{ComplexF64}` – vacuum response matrix - `plasma_pts::AbstractMatrix{Float64}` – plasma surface coordinates - `wall_pts::AbstractMatrix{Float64}` – wall surface coordinates - - `current_matrix::AbstractMatrix{ComplexF64}` – required when `compute_L=true` + - `I_v::AbstractMatrix{ComplexF64}` – required when `compute_L=true` This is designed to work with `ForceFreeStates.VacuumData` but does not depend on its concrete type (duck-typed on field names only). diff --git a/test/runtests_vacuum.jl b/test/runtests_vacuum.jl index 035ed96b..4417e408 100644 --- a/test/runtests_vacuum.jl +++ b/test/runtests_vacuum.jl @@ -369,7 +369,7 @@ @testset "nowall" begin inputs = _make_inputs() wall_settings = WallShapeSettings(shape="nowall") - wv, current_matrix, plasma_pts, wall_pts = compute_vacuum_response(inputs, wall_settings) + wv, I_v, plasma_pts, wall_pts = compute_vacuum_response(inputs, wall_settings) numpoints = inputs.mtheta * inputs.nzeta num_modes = length(inputs.m_modes) * length(inputs.n_modes) @@ -377,8 +377,8 @@ @test size(wv) == (num_modes, num_modes) @test eltype(wv) == ComplexF64 @test all(isfinite, wv) - @test size(current_matrix) == (num_modes, num_modes) - @test all(iszero, current_matrix) # compute_L=false by default + @test size(I_v) == (num_modes, num_modes) + @test all(iszero, I_v) # compute_L=false by default @test size(plasma_pts) == (numpoints, 3) @test all(isfinite, plasma_pts) @test size(wall_pts) == (numpoints, 3) @@ -390,26 +390,26 @@ @testset "nowall compute_L=true" begin inputs = _make_inputs() wall_settings = WallShapeSettings(shape="nowall") - wv, current_matrix, plasma_pts, wall_pts = compute_vacuum_response(inputs, wall_settings; compute_L=true) + wv, I_v, plasma_pts, wall_pts = compute_vacuum_response(inputs, wall_settings; compute_L=true) num_modes = length(inputs.m_modes) * length(inputs.n_modes) @test size(wv) == (num_modes, num_modes) @test all(isfinite, wv) - @test size(current_matrix) == (num_modes, num_modes) - @test all(isfinite, current_matrix) - @test !all(iszero, current_matrix) + @test size(I_v) == (num_modes, num_modes) + @test all(isfinite, I_v) + @test !all(iszero, I_v) @test isapprox(wv, wv', rtol=1e-12) end @testset "conformal wall" begin inputs = _make_inputs() wall_settings = WallShapeSettings(shape="conformal", a=0.5) - wv, current_matrix, plasma_pts, wall_pts = compute_vacuum_response(inputs, wall_settings) + wv, I_v, plasma_pts, wall_pts = compute_vacuum_response(inputs, wall_settings) numpoints = inputs.mtheta * inputs.nzeta num_modes = length(inputs.m_modes) * length(inputs.n_modes) @test size(wv) == (num_modes, num_modes) - @test all(iszero, current_matrix) + @test all(iszero, I_v) @test all(isfinite, plasma_pts) @test all(isfinite, wall_pts) # plasma_pts layout: col1=R, col2=0, col3=Z @@ -421,19 +421,19 @@ @testset "edge: single poloidal mode mpert=1" begin inputs = _make_inputs(m_modes=[1], n_modes=[1]) wall_settings = WallShapeSettings(shape="nowall") - wv, current_matrix, plasma_pts, wall_pts = compute_vacuum_response(inputs, wall_settings) + wv, I_v, plasma_pts, wall_pts = compute_vacuum_response(inputs, wall_settings) @test size(wv) == (1, 1) @test all(isfinite, wv) - @test size(current_matrix) == (1, 1) + @test size(I_v) == (1, 1) end @testset "edge: small mtheta" begin # Keep mtheta_eq=17 so boundary has enough points for periodic spline inputs = _make_inputs(mtheta=16, mtheta_eq=17) wall_settings = WallShapeSettings(shape="nowall") - wv, current_matrix, plasma_pts, wall_pts = compute_vacuum_response(inputs, wall_settings) + wv, I_v, plasma_pts, wall_pts = compute_vacuum_response(inputs, wall_settings) @test size(wv) == (2, 2) - @test size(current_matrix) == (2, 2) + @test size(I_v) == (2, 2) @test size(plasma_pts) == (16, 3) end @@ -442,17 +442,17 @@ # in-place entry populates caller-owned duck-typed (NamedTuple) storage identically. for wall_settings in (WallShapeSettings(shape="nowall"), WallShapeSettings(shape="conformal", a=0.5)) inputs = _make_inputs() - wv, current_matrix, pp, wp = compute_vacuum_response(inputs, wall_settings; compute_L=true) + wv, I_v, pp, wp = compute_vacuum_response(inputs, wall_settings; compute_L=true) numpoints = inputs.mtheta * inputs.nzeta num_modes = length(inputs.m_modes) * length(inputs.n_modes) vac = (wv=zeros(ComplexF64, num_modes, num_modes), - current_matrix=zeros(ComplexF64, num_modes, num_modes), + I_v=zeros(ComplexF64, num_modes, num_modes), plasma_pts=zeros(numpoints, 3), wall_pts=zeros(numpoints, 3)) compute_vacuum_response!(vac, inputs, wall_settings; compute_L=true) @test vac.wv ≈ wv - @test vac.current_matrix ≈ current_matrix + @test vac.I_v ≈ I_v @test vac.plasma_pts ≈ pp @test vac.wall_pts ≈ wp end @@ -594,15 +594,15 @@ @testset "compute_vacuum_response 3D nowall" begin inputs = _make_3d_inputs(mtheta=32, nzeta=32, mtheta_eq=17) wall_settings = WallShapeSettings(shape="nowall") - wv, current_matrix, plasma_pts, wall_pts = compute_vacuum_response(inputs, wall_settings) + wv, I_v, plasma_pts, wall_pts = compute_vacuum_response(inputs, wall_settings) numpoints = inputs.mtheta * inputs.nzeta num_modes = length(inputs.m_modes) * length(inputs.n_modes) @test size(wv) == (num_modes, num_modes) @test eltype(wv) == ComplexF64 @test all(isfinite, wv) - @test size(current_matrix) == (num_modes, num_modes) - @test all(iszero, current_matrix) + @test size(I_v) == (num_modes, num_modes) + @test all(iszero, I_v) @test size(plasma_pts) == (numpoints, 3) @test all(isfinite, plasma_pts) @test size(wall_pts) == (numpoints, 3) @@ -615,7 +615,7 @@ @testset "compute_vacuum_response 3D nonaxisymmetric boundary" begin inputs = _make_3d_nonaxis_inputs(mtheta=24, nzeta=24, mtheta_in=12, nzeta_in=12, mpert=2, nlow=0, npert=2) wall_settings = WallShapeSettings(shape="nowall") - wv, current_matrix, plasma_pts, wall_pts = compute_vacuum_response(inputs, wall_settings) + wv, I_v, plasma_pts, wall_pts = compute_vacuum_response(inputs, wall_settings) numpoints = inputs.mtheta * inputs.nzeta num_modes = length(inputs.m_modes) * length(inputs.n_modes) @@ -623,8 +623,8 @@ @test size(wv) == (num_modes, num_modes) @test eltype(wv) == ComplexF64 @test all(isfinite, wv) - @test size(current_matrix) == (num_modes, num_modes) - @test all(iszero, current_matrix) + @test size(I_v) == (num_modes, num_modes) + @test all(iszero, I_v) @test size(plasma_pts) == (numpoints, 3) @test all(isfinite, plasma_pts) @test size(wall_pts) == (numpoints, 3) @@ -634,12 +634,12 @@ @testset "compute_vacuum_response 3D conformal wall" begin inputs = _make_3d_inputs(mtheta=32, nzeta=32, mtheta_eq=17) wall_settings = WallShapeSettings(shape="conformal", a=0.3) - wv, current_matrix, plasma_pts, wall_pts = compute_vacuum_response(inputs, wall_settings) + wv, I_v, plasma_pts, wall_pts = compute_vacuum_response(inputs, wall_settings) numpoints = inputs.mtheta * inputs.nzeta num_modes = length(inputs.m_modes) * length(inputs.n_modes) @test size(wv) == (num_modes, num_modes) - @test size(current_matrix) == (num_modes, num_modes) + @test size(I_v) == (num_modes, num_modes) @test all(isfinite, plasma_pts) @test all(isfinite, wall_pts) # Wall and plasma should differ (conformal wall offset from plasma) From 5f08ed8530e8a62159abc7013bd8d2abbdc38bee Mon Sep 17 00:00:00 2001 From: Jake Halpern Date: Mon, 13 Jul 2026 10:53:15 -0400 Subject: [PATCH 04/10] GPEC - minor - name change --- src/PerturbedEquilibrium/SingularCoupling.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/PerturbedEquilibrium/SingularCoupling.jl b/src/PerturbedEquilibrium/SingularCoupling.jl index 483635c7..73c091ae 100644 --- a/src/PerturbedEquilibrium/SingularCoupling.jl +++ b/src/PerturbedEquilibrium/SingularCoupling.jl @@ -416,7 +416,7 @@ Invert the vacuum surface-current matrix into surface inductance `L`. Vacuum returns `I_v` from the Green’s functions, this helper applies GPEC normalization, regularization, inversion, and Hermitianization. """ -function calc_surface_inductance(I_v::AbstractMatrix{ComplexF64})::Matrix{ComplexF64} +function calc_surface_inductance(I_v::Matrix{ComplexF64})::Matrix{ComplexF64} mpert = size(I_v, 1) μ₀ = 4π * 1e-7 From 7721276d78fd21f4d5ef7dba70fe46882de47733 Mon Sep 17 00:00:00 2001 From: Jake Halpern Date: Fri, 31 Jul 2026 16:31:08 -0400 Subject: [PATCH 05/10] PE - MINOR - removing non-Fortran consistent try/catch routines and regularization in L calculation and moving to ResponseMatrices.jl --- src/PerturbedEquilibrium/ResponseMatrices.jl | 19 ++++++++++ src/PerturbedEquilibrium/SingularCoupling.jl | 40 -------------------- src/Vacuum/Vacuum.jl | 13 +++---- 3 files changed, 25 insertions(+), 47 deletions(-) diff --git a/src/PerturbedEquilibrium/ResponseMatrices.jl b/src/PerturbedEquilibrium/ResponseMatrices.jl index a4b798a6..0de57868 100644 --- a/src/PerturbedEquilibrium/ResponseMatrices.jl +++ b/src/PerturbedEquilibrium/ResponseMatrices.jl @@ -265,6 +265,25 @@ function calc_plasma_inductance( return inv(temp2) end +""" + calc_surface_inductance(I_v::AbstractMatrix{ComplexF64})::Matrix{ComplexF64} + +Invert the vacuum surface-current matrix into surface inductance `L`. + +Vacuum returns `I_v` from the Green’s functions, this helper applies GPEC normalization, +regularization, inversion, and Hermitianization. +""" +function calc_surface_inductance(I_v::Matrix{ComplexF64})::Matrix{ComplexF64} + + μ₀ = 4π * 1e-7 + I_v ./= μ₀ * (2π)^2 + # L = flux * inv(Iᵛ) with flux = I + L_surf = inv(I_v) + hermitianpart!(L_surf) + + return L_surf +end + """ calc_permeability( plasma_inductance::Matrix{ComplexF64}, diff --git a/src/PerturbedEquilibrium/SingularCoupling.jl b/src/PerturbedEquilibrium/SingularCoupling.jl index 56a46c9a..8c1aadb8 100644 --- a/src/PerturbedEquilibrium/SingularCoupling.jl +++ b/src/PerturbedEquilibrium/SingularCoupling.jl @@ -424,46 +424,6 @@ function compute_current_density( return j_c end -""" - calc_surface_inductance(I_v::AbstractMatrix{ComplexF64})::Matrix{ComplexF64} - -Invert the vacuum surface-current matrix into surface inductance `L`. - -Vacuum returns `I_v` from the Green’s functions, this helper applies GPEC normalization, -regularization, inversion, and Hermitianization. -""" -function calc_surface_inductance(I_v::Matrix{ComplexF64})::Matrix{ComplexF64} - - mpert = size(I_v, 1) - μ₀ = 4π * 1e-7 - L_surf = zeros(ComplexF64, mpert, mpert) - - current_mag = maximum(abs.(I_v)) - - 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 - I_v ./= μ₀ * (2π)^2 - I_v += 1e-12 * current_mag * I - # L = flux * inv(Iᵛ) with flux = I - L_surf = inv(I_v) - # 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, diff --git a/src/Vacuum/Vacuum.jl b/src/Vacuum/Vacuum.jl index 1e651b5b..b98a60ab 100644 --- a/src/Vacuum/Vacuum.jl +++ b/src/Vacuum/Vacuum.jl @@ -127,29 +127,27 @@ Green's functions are internal scratch only. ldiv!(lu!(grad_green_interior), grri) # Surface-current matrix, Park 2007 eq. 21b: μ₀I^v = χ^(vi) - χ^(vo) = grri - grre - # The minus is because VACUUM builds the operators in its CW-θ frame while GPEC uses CCW-θ, - # flipping the outward-normal sign. + # They are flipped because VACUUM builds the operators in its CW-θ frame while GPEC + # uses CCW-θ, flipping the outward-normal sign. I_v_block = @view vac_data.I_v[block_idx, block_idx] - @views g_sum = -(grri[1:num_points_surf, :] .- grre[1:num_points_surf, :]) + @views g_sum = grre[1:num_points_surf, :] .- grri[1:num_points_surf, :] mul!(I_v_block, ft.basis, g_sum) I_v_block ./= num_points_surf # Iᵛ becomes the surface inductance L (Φ = L·Iᵛ), a one-sided operator that # multiplies GPEC-convention mode vectors. Kernels are real, so conj rotates B·A·B† into - # the GPEC (CCW-θ, conjugate-of-RH) basis. Unlike wv this is not an invariant scalar. + # the GPEC (CCW-θ, conjugate-of-RH) basis. Wv is Hermitian so no conjugation is needed I_v_block .= conj.(I_v_block) else # Only need exterior system for wv ldiv!(lu!(grad_green), grre) end - # wv = B·(D_ext⁻¹𝒢)·B† is a Hermitian energy generator (Chance 1997 eqs. 125-126); it is - # contracted on both sides to the real, convention-invariant vacuum energy, so no conjugation. + # Project exterior kernel onto observer basis exp(-i*(mθ - nν)) and scale to get the response matrix mul!(wv_block, ft.basis, @view(grre[1:num_points_surf, :])) wv_block .*= 4π^2 / num_points_surf end - # δW_v = ξ† Wᵛ ξ is real, so Wᵛ must be Hermitian; remove any residual from discretization _symmetrize_vacuum_energy!(vac_data.wv) @@ -268,6 +266,7 @@ interior variant `-D + 2I` for the interior columns, then scatter back into the end end + # δW_v = ξ† Wᵛ ξ is real, so Wᵛ must be Hermitian; remove any residual from discretization _symmetrize_vacuum_energy!(vac_data.wv) # Populate coordinate arrays From 9434b64fc2c24b29fb8812b27c78bdc50cf8aaaa Mon Sep 17 00:00:00 2001 From: Jake Halpern Date: Fri, 31 Jul 2026 16:56:37 -0400 Subject: [PATCH 06/10] VACUUM - MINOR - moving Hermitianization into vacuum module and taking advantage of it to remove some lines --- src/PerturbedEquilibrium/ResponseMatrices.jl | 4 +- src/Vacuum/Vacuum.jl | 49 ++++++++------------ 2 files changed, 20 insertions(+), 33 deletions(-) diff --git a/src/PerturbedEquilibrium/ResponseMatrices.jl b/src/PerturbedEquilibrium/ResponseMatrices.jl index 0de57868..d88f57fa 100644 --- a/src/PerturbedEquilibrium/ResponseMatrices.jl +++ b/src/PerturbedEquilibrium/ResponseMatrices.jl @@ -278,10 +278,8 @@ function calc_surface_inductance(I_v::Matrix{ComplexF64})::Matrix{ComplexF64} μ₀ = 4π * 1e-7 I_v ./= μ₀ * (2π)^2 # L = flux * inv(Iᵛ) with flux = I - L_surf = inv(I_v) - hermitianpart!(L_surf) - return L_surf + return inv(I_v) end """ diff --git a/src/Vacuum/Vacuum.jl b/src/Vacuum/Vacuum.jl index b98a60ab..c8cf0d46 100644 --- a/src/Vacuum/Vacuum.jl +++ b/src/Vacuum/Vacuum.jl @@ -24,32 +24,25 @@ export extract_plasma_surface_at_psi export PlasmaGeometry # Relative anti-Hermitian residual above which we warn that the vacuum grid should be refined. -const _WV_HERMITICITY_WARN_TOL = 1e-4 +const _HERMITICITY_WARN_TOL = 1e-4 """ - _symmetrize_vacuum_energy!(wv) + _warn_and_symmetrize!(mat, name) -Enforce Hermiticity of the vacuum energy matrix Wᵛ in place. - -Wᵛ is the generator of the vacuum magnetic energy, δW_v = ξ† Wᵛ ξ. Because that energy is a real -number for every perturbation ξ, the exact operator is Hermitian. The finite-resolution -boundary-integral quadrature (finite `mtheta`/`nzeta`) breaks exact Hermiticity, leaving a small -anti-Hermitian residual that is a pure discretization artifact and vanishes as the grid is refined. -We replace Wᵛ by its Hermitian part to restore this physical property, warning when the residual -is large enough that the vacuum grid should be refined. +Replace `mat` by its Hermitian part in place, warning first if the anti-Hermitian residual exceeds +`_HERMITICITY_WARN_TOL`. """ -function _symmetrize_vacuum_energy!(wv::AbstractMatrix) - - herm_norm = norm(wv + wv') +function _warn_and_symmetrize!(mat::AbstractMatrix, name::String) + herm_norm = norm(mat + mat') if herm_norm > 0 - # Relative anti-Hermitian residual ‖½(W−W†)‖/‖½(W+W†)‖ - rel_residual = norm(wv - wv') / herm_norm - if rel_residual > _WV_HERMITICITY_WARN_TOL - @warn "Vacuum energy matrix Wᵛ is non-Hermitian above tolerance $(rel_residual) > $(_WV_HERMITICITY_WARN_TOL) before " * + # Relative anti-Hermitian residual ‖½(M−M†)‖/‖½(M+M†)‖ + rel_residual = norm(mat - mat') / herm_norm + if rel_residual > _HERMITICITY_WARN_TOL + @warn "$name is non-Hermitian above tolerance $(rel_residual) > $(_HERMITICITY_WARN_TOL) before " * "symmetrization. Increase vacuum grid resolution to reduce it." end end - hermitianpart!(wv) + hermitianpart!(mat) end """ @@ -133,11 +126,6 @@ Green's functions are internal scratch only. @views g_sum = grre[1:num_points_surf, :] .- grri[1:num_points_surf, :] mul!(I_v_block, ft.basis, g_sum) I_v_block ./= num_points_surf - - # Iᵛ becomes the surface inductance L (Φ = L·Iᵛ), a one-sided operator that - # multiplies GPEC-convention mode vectors. Kernels are real, so conj rotates B·A·B† into - # the GPEC (CCW-θ, conjugate-of-RH) basis. Wv is Hermitian so no conjugation is needed - I_v_block .= conj.(I_v_block) else # Only need exterior system for wv ldiv!(lu!(grad_green), grre) @@ -148,8 +136,9 @@ Green's functions are internal scratch only. wv_block .*= 4π^2 / num_points_surf end - # δW_v = ξ† Wᵛ ξ is real, so Wᵛ must be Hermitian; remove any residual from discretization - _symmetrize_vacuum_energy!(vac_data.wv) + # Remove any non-Hermitian residual from Hermitian matrices due to discretization + _warn_and_symmetrize!(vac_data.wv, "Wᵛ") + compute_L && _warn_and_symmetrize!(vac_data.I_v, "Iᵛ") # Populate coordinate arrays @views begin @@ -191,9 +180,8 @@ interior variant `-D + 2I` for the interior columns, then scatter back into the (; mtheta, nzeta, nfp, m_modes, n_modes) = inputs fill!(vac_data.wv, 0) fill!(vac_data.I_v, 0) - if compute_L - @warn "compute_L=true is not supported for 3D vacuum response; I_v left as zeros" maxlog=1 - end + + compute_L && @warn "compute_L=true is not supported for 3D vacuum response; I_v left as zeros" maxlog=1 # Full-torus geometry for source surface; observers are restricted to one field period full = expand_field_periods(inputs) @@ -266,8 +254,9 @@ interior variant `-D + 2I` for the interior columns, then scatter back into the end end - # δW_v = ξ† Wᵛ ξ is real, so Wᵛ must be Hermitian; remove any residual from discretization - _symmetrize_vacuum_energy!(vac_data.wv) + # Remove any non-Hermitian residual from Hermitian matrices due to discretization + _warn_and_symmetrize!(vac_data.wv, "Wᵛ") + compute_L && _warn_and_symmetrize!(vac_data.I_v, "Iᵛ") # Populate coordinate arrays vac_data.plasma_pts .= plasma_surf.r From a2417ce5b0ff6d3711ca35e1b8467afbbf7e15cc Mon Sep 17 00:00:00 2001 From: Jake Halpern Date: Fri, 31 Jul 2026 17:06:03 -0400 Subject: [PATCH 07/10] PE - MINOR - updating docstring of surface inductance calculation --- src/PerturbedEquilibrium/ResponseMatrices.jl | 16 ++++++---------- 1 file changed, 6 insertions(+), 10 deletions(-) diff --git a/src/PerturbedEquilibrium/ResponseMatrices.jl b/src/PerturbedEquilibrium/ResponseMatrices.jl index d88f57fa..d6970657 100644 --- a/src/PerturbedEquilibrium/ResponseMatrices.jl +++ b/src/PerturbedEquilibrium/ResponseMatrices.jl @@ -266,20 +266,16 @@ function calc_plasma_inductance( end """ - calc_surface_inductance(I_v::AbstractMatrix{ComplexF64})::Matrix{ComplexF64} + calc_surface_inductance(I_v::Matrix{ComplexF64})::Matrix{ComplexF64} -Invert the vacuum surface-current matrix into surface inductance `L`. - -Vacuum returns `I_v` from the Green’s functions, this helper applies GPEC normalization, -regularization, inversion, and Hermitianization. +Surface inductance from the vacuum surface-current matrix, Park 2007 eq. 7: `Φ_x = L·I_v`. +Columns of `I_v` are driven by unit flux harmonics which are consistent since the eq. 3 +weight `1/(J|∇ψ|)` cancels the Jacobian in the vacuum solver's source density — so `Φ_x = 𝕀` +and `L = I_v⁻¹`. `I_v` is normalized in place by `μ₀(2π)²`. """ function calc_surface_inductance(I_v::Matrix{ComplexF64})::Matrix{ComplexF64} - μ₀ = 4π * 1e-7 - I_v ./= μ₀ * (2π)^2 - # L = flux * inv(Iᵛ) with flux = I - - return inv(I_v) + return inv(I_v) .* μ₀ * (2π)^2 end """ From 0a3f5880b2246a80e2cdaa0d72068597dff88c29 Mon Sep 17 00:00:00 2001 From: Jake Halpern Date: Fri, 31 Jul 2026 17:16:33 -0400 Subject: [PATCH 08/10] VACUUM - MINOR - replacing compute_L with the more aptly named compute_Iv --- docs/src/workflow.md | 2 +- src/PerturbedEquilibrium/Response.jl | 2 +- src/PerturbedEquilibrium/SingularCoupling.jl | 6 +-- src/Vacuum/Kernel2D.jl | 2 +- src/Vacuum/Vacuum.jl | 40 ++++++++++---------- test/runtests_vacuum.jl | 10 ++--- 6 files changed, 31 insertions(+), 31 deletions(-) diff --git a/docs/src/workflow.md b/docs/src/workflow.md index 37d79deb..0aa3c843 100644 --- a/docs/src/workflow.md +++ b/docs/src/workflow.md @@ -59,7 +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) -- `I_v` — Vacuum surface-current matrix Iᵛ when `compute_L=true` (otherwise zeros); PerturbedEquilibrium inverts this to surface inductance `L` +- `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) diff --git a/src/PerturbedEquilibrium/Response.jl b/src/PerturbedEquilibrium/Response.jl index bc179449..94318de3 100644 --- a/src/PerturbedEquilibrium/Response.jl +++ b/src/PerturbedEquilibrium/Response.jl @@ -39,7 +39,7 @@ function compute_plasma_response!( 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") - _, I_v, _, _ = Vacuum.compute_vacuum_response(vac_input_2d, wall_nowall; compute_L=true) + _, I_v, _, _ = Vacuum.compute_vacuum_response(vac_input_2d, wall_nowall; compute_Iv=true) surface_inductance = calc_surface_inductance(I_v) permeability = calc_permeability(plasma_inductance, surface_inductance) diff --git a/src/PerturbedEquilibrium/SingularCoupling.jl b/src/PerturbedEquilibrium/SingularCoupling.jl index 8c1aadb8..73431be9 100644 --- a/src/PerturbedEquilibrium/SingularCoupling.jl +++ b/src/PerturbedEquilibrium/SingularCoupling.jl @@ -172,10 +172,10 @@ function compute_singular_coupling_metrics!( # Surface-current matrix at this surface for this n (once per pair) vac_input = Vacuum.VacuumInput(equil, sing_surf.psifac, mtheta, 1, mlow:mhigh, [nn]) - _, I_v, _, _ = Vacuum.compute_vacuum_response(vac_input, wall_settings; compute_L=true) - - # Precompute L_surf; only the (m_res, m_res) diagonal element is needed for singflx + _, I_v, _, _ = Vacuum.compute_vacuum_response(vac_input, wall_settings; compute_Iv=true) L_surf = calc_surface_inductance(I_v) + + # Extract 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] diff --git a/src/Vacuum/Kernel2D.jl b/src/Vacuum/Kernel2D.jl index a5b03397..c290246e 100644 --- a/src/Vacuum/Kernel2D.jl +++ b/src/Vacuum/Kernel2D.jl @@ -217,7 +217,7 @@ grad_greenfunction is not zeroed since it fills a different block of the # Normals need to point outward from vacuum region. In VACUUM clockwise θ convention, normal points # out of vacuum for wall but inward for plasma, so we multiply by -1 for plasma sources. # This normal orientation sets the sign of the surface current K = n̂×∇χ (Chance 1997 eq. 184), - # and is the origin of the leading minus in the current-matrix build (Vacuum.jl, compute_L). + # and is the origin of the leading minus in the current-matrix build. if source isa PlasmaGeometry grad_greenfunction_block .*= -1 end diff --git a/src/Vacuum/Vacuum.jl b/src/Vacuum/Vacuum.jl index c8cf0d46..de653e07 100644 --- a/src/Vacuum/Vacuum.jl +++ b/src/Vacuum/Vacuum.jl @@ -46,23 +46,23 @@ function _warn_and_symmetrize!(mat::AbstractMatrix, name::String) end """ - _compute_vacuum_response_2d!(vac_data, inputs::VacuumInput, wall_settings::WallShapeSettings; compute_L=false) + _compute_vacuum_response_2d!(vac_data, inputs::VacuumInput, wall_settings::WallShapeSettings; compute_Iv=false) 2D (axisymmetric) vacuum response calculation. Each toroidal mode `n` decouples in 2D geometry, so the routine loops over `inputs.n_modes`, building the double-/single-layer operators, solving the exterior system for `wv`, and -optionally the interior system to build `I_v` when `compute_L=true`. +optionally the interior system to build `I_v` when `compute_Iv=true`. Green's functions are internal scratch only. """ -@with_pool pool function _compute_vacuum_response_2d!(vac_data, inputs::VacuumInput, wall_settings::WallShapeSettings; compute_L::Bool=false) +@with_pool pool function _compute_vacuum_response_2d!(vac_data, inputs::VacuumInput, wall_settings::WallShapeSettings; compute_Iv::Bool=false) mpert = length(inputs.m_modes) mlow = inputs.m_modes[1] num_points_surf = inputs.mtheta vac_data.wv .= 0 - compute_L && (vac_data.I_v .= 0) + compute_Iv && (vac_data.I_v .= 0) # Form the plasma and wall geometries plasma_surf = PlasmaGeometry(inputs) @@ -101,7 +101,7 @@ Green's functions are internal scratch only. mul!(view(grre, (num_points_surf+1):num_points_total, :), green_temp, ft.basis') end - if compute_L + if compute_Iv # Copy RHS before exterior solve overwrites grre; keep a kernel copy for interior grri = similar!(pool, grre) grri .= grre @@ -138,7 +138,7 @@ Green's functions are internal scratch only. # Remove any non-Hermitian residual from Hermitian matrices due to discretization _warn_and_symmetrize!(vac_data.wv, "Wᵛ") - compute_L && _warn_and_symmetrize!(vac_data.I_v, "Iᵛ") + compute_Iv && _warn_and_symmetrize!(vac_data.I_v, "Iᵛ") # Populate coordinate arrays @views begin @@ -152,7 +152,7 @@ Green's functions are internal scratch only. end """ - _compute_vacuum_response_3d!(vac_data, inputs::VacuumInput, wall_settings::WallShapeSettings; compute_L=false) + _compute_vacuum_response_3d!(vac_data, inputs::VacuumInput, wall_settings::WallShapeSettings; compute_Iv=false) 3D (`inputs.nzeta > 1`) vacuum response via block-circulant field-period reduction. For `nfp == 1` the block-circulant assembly and residue-class loop are skipped in favour of a more efficient @@ -170,18 +170,18 @@ class needs one solve `wv[class k] = (4π²/M)·E_localᴴ·(D̂ₖ \\ Ŝₖ)|_p block-row of the operators is built, so the kernel cost drops by `nfp` and the dense `O(N³)` factorization is replaced by per-class `O(M³)` solves (`M = N/nfp`). -Only `wv` is produced currently; `I_v` is left zeroed when `compute_L=true` +Only `wv` is produced currently; `I_v` is left zeroed when `compute_Iv=true` (surface-current / inductance not yet supported in 3D). Extension point: per residue class, apply the per-period basis to `D̂ₖ⁻¹Ŝₖ` for the exterior columns and the interior variant `-D + 2I` for the interior columns, then scatter back into the `[2N × 2·num_modes]` arrays. """ -@with_pool pool function _compute_vacuum_response_3d!(vac_data, inputs::VacuumInput, wall_settings::WallShapeSettings; compute_L::Bool=false) +@with_pool pool function _compute_vacuum_response_3d!(vac_data, inputs::VacuumInput, wall_settings::WallShapeSettings; compute_Iv::Bool=false) (; mtheta, nzeta, nfp, m_modes, n_modes) = inputs fill!(vac_data.wv, 0) fill!(vac_data.I_v, 0) - compute_L && @warn "compute_L=true is not supported for 3D vacuum response; I_v left as zeros" maxlog=1 + compute_Iv && @warn "compute_Iv=true is not supported for 3D vacuum response; I_v left as zeros" maxlog=1 # Full-torus geometry for source surface; observers are restricted to one field period full = expand_field_periods(inputs) @@ -256,7 +256,7 @@ interior variant `-D + 2I` for the interior columns, then scatter back into the # Remove any non-Hermitian residual from Hermitian matrices due to discretization _warn_and_symmetrize!(vac_data.wv, "Wᵛ") - compute_L && _warn_and_symmetrize!(vac_data.I_v, "Iᵛ") + compute_Iv && _warn_and_symmetrize!(vac_data.I_v, "Iᵛ") # Populate coordinate arrays vac_data.plasma_pts .= plasma_surf.r @@ -264,7 +264,7 @@ interior variant `-D + 2I` for the interior columns, then scatter back into the end """ - compute_vacuum_response(inputs::VacuumInput, wall_settings::WallShapeSettings; compute_L=false) + compute_vacuum_response(inputs::VacuumInput, wall_settings::WallShapeSettings; compute_Iv=false) Allocate and return the vacuum response matrix and optional surface-current matrix for the given vacuum inputs. Thin allocating wrapper around the in-place [`compute_vacuum_response!`]: it sizes the @@ -275,10 +275,10 @@ in-place method to avoid extra heap allocations. # Returns - `wv`: complex vacuum response matrix (`num_modes × num_modes`). - - `I_v`: Vacuum surface-current matrix (`num_modes × num_modes`); zeros unless `compute_L=true`. + - `I_v`: Vacuum surface-current matrix (`num_modes × num_modes`); zeros unless `compute_Iv=true`. - `plasma_pts`, `wall_pts`: surface coordinate arrays. """ -function compute_vacuum_response(inputs::VacuumInput, wall_settings::WallShapeSettings; compute_L::Bool=false) +function compute_vacuum_response(inputs::VacuumInput, wall_settings::WallShapeSettings; compute_Iv::Bool=false) num_points = inputs.mtheta * inputs.nzeta * inputs.nfp # mtheta for 2D num_modes = length(inputs.m_modes) * length(inputs.n_modes) @@ -289,13 +289,13 @@ function compute_vacuum_response(inputs::VacuumInput, wall_settings::WallShapeSe plasma_pts=zeros(num_points, 3), wall_pts=zeros(num_points, 3) ) - compute_vacuum_response!(vac, inputs, wall_settings; compute_L) + compute_vacuum_response!(vac, inputs, wall_settings; compute_Iv) return vac.wv, vac.I_v, vac.plasma_pts, vac.wall_pts end """ - compute_vacuum_response!(vac_data, inputs::VacuumInput, wall_settings::WallShapeSettings; compute_L=false) + compute_vacuum_response!(vac_data, inputs::VacuumInput, wall_settings::WallShapeSettings; compute_Iv=false) In-place variant that computes the vacuum response and directly populates the arrays stored in `vac_data`. Dispatches on dimensionality only: 2D (`inputs.nzeta == 1`) routes to @@ -307,16 +307,16 @@ sizes: - `wv::AbstractMatrix{ComplexF64}` – vacuum response matrix - `plasma_pts::AbstractMatrix{Float64}` – plasma surface coordinates - `wall_pts::AbstractMatrix{Float64}` – wall surface coordinates - - `I_v::AbstractMatrix{ComplexF64}` – required when `compute_L=true` + - `I_v::AbstractMatrix{ComplexF64}` – required when `compute_Iv=true` This is designed to work with `ForceFreeStates.VacuumData` but does not depend on its concrete type (duck-typed on field names only). """ -function compute_vacuum_response!(vac_data, inputs::VacuumInput, wall_settings::WallShapeSettings; compute_L::Bool=false) +function compute_vacuum_response!(vac_data, inputs::VacuumInput, wall_settings::WallShapeSettings; compute_Iv::Bool=false) if inputs.nzeta == 1 - _compute_vacuum_response_2d!(vac_data, inputs, wall_settings; compute_L) + _compute_vacuum_response_2d!(vac_data, inputs, wall_settings; compute_Iv) else - _compute_vacuum_response_3d!(vac_data, inputs, wall_settings; compute_L) + _compute_vacuum_response_3d!(vac_data, inputs, wall_settings; compute_Iv) end end diff --git a/test/runtests_vacuum.jl b/test/runtests_vacuum.jl index 48fb33be..1e84546b 100644 --- a/test/runtests_vacuum.jl +++ b/test/runtests_vacuum.jl @@ -377,7 +377,7 @@ @test eltype(wv) == ComplexF64 @test all(isfinite, wv) @test size(I_v) == (num_modes, num_modes) - @test all(iszero, I_v) # compute_L=false by default + @test all(iszero, I_v) # compute_Iv=false by default @test size(plasma_pts) == (numpoints, 3) @test all(isfinite, plasma_pts) @test size(wall_pts) == (numpoints, 3) @@ -386,10 +386,10 @@ @test isapprox(wv, wv', rtol=1e-12) end - @testset "nowall compute_L=true" begin + @testset "nowall compute_Iv=true" begin inputs = _make_inputs() wall_settings = WallShapeSettings(shape="nowall") - wv, I_v, plasma_pts, wall_pts = compute_vacuum_response(inputs, wall_settings; compute_L=true) + wv, I_v, plasma_pts, wall_pts = compute_vacuum_response(inputs, wall_settings; compute_Iv=true) num_modes = length(inputs.m_modes) * length(inputs.n_modes) @test size(wv) == (num_modes, num_modes) @@ -441,14 +441,14 @@ # in-place entry populates caller-owned duck-typed (NamedTuple) storage identically. for wall_settings in (WallShapeSettings(shape="nowall"), WallShapeSettings(shape="conformal", a=0.5)) inputs = _make_inputs() - wv, I_v, pp, wp = compute_vacuum_response(inputs, wall_settings; compute_L=true) + wv, I_v, pp, wp = compute_vacuum_response(inputs, wall_settings; compute_Iv=true) numpoints = inputs.mtheta * inputs.nzeta num_modes = length(inputs.m_modes) * length(inputs.n_modes) vac = (wv=zeros(ComplexF64, num_modes, num_modes), I_v=zeros(ComplexF64, num_modes, num_modes), plasma_pts=zeros(numpoints, 3), wall_pts=zeros(numpoints, 3)) - compute_vacuum_response!(vac, inputs, wall_settings; compute_L=true) + compute_vacuum_response!(vac, inputs, wall_settings; compute_Iv=true) @test vac.wv ≈ wv @test vac.I_v ≈ I_v From 7ee5d49b0848408d965d07ef4cd3cca199326e04 Mon Sep 17 00:00:00 2001 From: Jake Halpern Date: Fri, 31 Jul 2026 17:17:12 -0400 Subject: [PATCH 09/10] reducing diff --- src/ForceFreeStates/ForceFreeStatesStructs.jl | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/ForceFreeStates/ForceFreeStatesStructs.jl b/src/ForceFreeStates/ForceFreeStatesStructs.jl index a146e6e0..bb188da8 100644 --- a/src/ForceFreeStates/ForceFreeStatesStructs.jl +++ b/src/ForceFreeStates/ForceFreeStatesStructs.jl @@ -399,8 +399,7 @@ Populated in `Free.jl`. - `ev::Vector{ComplexF64}` - Vacuum energy per eigenmode (power quotient v†·wv·v with v†·N·v = 1) - `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 vacuum matrix wv, clamped to zero - - `plasma_pts::Array{Float64, 3}` - Cartesian coordinates of plasma points, shape (mthvac * nzvac) × 3 for (x, y, z) + - `vacuum_eigenvalue::Float64` - Least stable (minimum) eigenvalue of the pencil (wv, N), clamped to zero - `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) """ @kwdef mutable struct VacuumData From 15c716d66b005ad324b3cc10de10395c7ced6eaf Mon Sep 17 00:00:00 2001 From: Jake Halpern Date: Fri, 31 Jul 2026 17:17:57 -0400 Subject: [PATCH 10/10] reducing diff --- src/ForceFreeStates/ForceFreeStatesStructs.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/ForceFreeStates/ForceFreeStatesStructs.jl b/src/ForceFreeStates/ForceFreeStatesStructs.jl index a146e6e0..31e3c6f6 100644 --- a/src/ForceFreeStates/ForceFreeStatesStructs.jl +++ b/src/ForceFreeStates/ForceFreeStatesStructs.jl @@ -399,7 +399,7 @@ Populated in `Free.jl`. - `ev::Vector{ComplexF64}` - Vacuum energy per eigenmode (power quotient v†·wv·v with v†·N·v = 1) - `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 vacuum matrix wv, clamped to zero + - `vacuum_eigenvalue::Float64` - Least stable (minimum) eigenvalue of the pencil (wv, N), clamped to zero - `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) """