diff --git a/docs/src/vacuum.md b/docs/src/vacuum.md index 46ca5f44..633df587 100644 --- a/docs/src/vacuum.md +++ b/docs/src/vacuum.md @@ -74,7 +74,7 @@ wall_settings = GeneralizedPerturbedEquilibrium.Vacuum.WallShapeSettings( ) # Compute vacuum response matrix -wv, grri, xzpts = GeneralizedPerturbedEquilibrium.Vacuum.compute_vacuum_response(inputs, wall_settings) +wv, _, _, _ = GeneralizedPerturbedEquilibrium.Vacuum.compute_vacuum_response(inputs, wall_settings) ``` ### Vacuum Field Calculation at Observation Points diff --git a/docs/src/workflow.md b/docs/src/workflow.md index d872dcea..0aa3c843 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) +- `I_v` — Vacuum surface-current matrix Iᵛ when `compute_Iv=true` (otherwise zeros); PerturbedEquilibrium inverts this to surface inductance `L` **Key references**: [Chance et al. (1997)](citations.md#Vacuum-Module), [Chance et al. (2007)](citations.md#Vacuum-Module) @@ -182,5 +181,5 @@ All results are written to a single HDF5 file (default: `gpec.h5`). The file is | `locstab/` | Local stability: Mercier criterion, shear | | `integration/` | ODE integration results: energy matrices, eigenvalues | | `singular/` | Per-surface data: ψ_s, m/n, Δ', small solution coefficients | -| `vacuum/` | Vacuum response matrices: wv, grri, grre | +| `vacuum/` | Vacuum response matrices: wv (and I_v when computed for PE) | | `perturbed/` | Perturbed equilibrium: ξ, b in mode space, island diagnostics | diff --git a/src/ForceFreeStates/ForceFreeStatesStructs.jl b/src/ForceFreeStates/ForceFreeStatesStructs.jl index 7f727ea5..31e3c6f6 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 @@ -404,8 +400,6 @@ Populated in `Free.jl`. - `et::Vector{ComplexF64}` - Total energy eigenvalues of the pencil (W, N): power-normalized and invariant to the working (Jacobian) coordinate; et = ep + ev per mode - `n_tor_idx::Vector{Int}` - 0-based toroidal mode number index of each sorted eigenvalue (numpert_total). Needed in `write_imas` - `vacuum_eigenvalue::Float64` - Least stable (minimum) eigenvalue of the pencil (wv, N), clamped to zero - - `grri::Array{ComplexF64, 2}` - Interior Green's function matrices (2 * mthvac * nzvac × numpert_total) - - `grre::Array{ComplexF64, 2}` - Exterior Green's function matrices (2 * mthvac * nzvac × numpert_total) - `plasma_pts::Array{Float64, 3}` - Cartesian coordinates of plasma points, shape (mthvac * nzvac) × 3 for (x, y, z) - `wall_pts::Array{Float64, 3}` - Cartesian coordinates of wall points, shape (mthvac * nzvac) × 3 for (x, y, z) """ @@ -423,8 +417,6 @@ Populated in `Free.jl`. et::Vector{ComplexF64} = Vector{ComplexF64}(undef, numpert_total) n_tor_idx::Vector{Int} = zeros(Int, numpert_total) vacuum_eigenvalue::Float64 = NaN - grri::Array{ComplexF64,2} = Array{ComplexF64}(undef, 2 * numpoints, numpert_total) - grre::Array{ComplexF64,2} = Array{ComplexF64}(undef, 2 * numpoints, numpert_total) plasma_pts::Array{Float64,2} = Array{Float64}(undef, numpoints, 3) wall_pts::Array{Float64,2} = Array{Float64}(undef, numpoints, 3) end diff --git a/src/ForceFreeStates/Free.jl b/src/ForceFreeStates/Free.jl index 52acb9ba..0ce2af36 100644 --- a/src/ForceFreeStates/Free.jl +++ b/src/ForceFreeStates/Free.jl @@ -175,7 +175,7 @@ q-window minimum. # Compute raw vacuum matrix at the actual scan psi (singfac NOT applied; free_compute_total applies it analytically) vac_inputs = Vacuum.VacuumInput(equil, psi_array[i], ctrl.mthvac, ctrl.nzvac, intr.mlow:intr.mhigh, intr.nlow:intr.nhigh) - wv, _, _, _, _ = Vacuum.compute_vacuum_response(vac_inputs, intr.wall_settings) + wv, _, _, _ = Vacuum.compute_vacuum_response(vac_inputs, intr.wall_settings) @views wv_array[i, :, :] .= wv end diff --git a/src/PerturbedEquilibrium/Response.jl b/src/PerturbedEquilibrium/Response.jl index ba282724..94318de3 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) + _, 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) # Reluctance ϱ = L⁻¹·(Λ† − L)·L⁻¹ (Fortran gpresp_reluct: diff_indmats = CONJG(TRANSPOSE(plas_indmats)) − surf_indmats). diff --git a/src/PerturbedEquilibrium/ResponseMatrices.jl b/src/PerturbedEquilibrium/ResponseMatrices.jl index a4b798a6..d6970657 100644 --- a/src/PerturbedEquilibrium/ResponseMatrices.jl +++ b/src/PerturbedEquilibrium/ResponseMatrices.jl @@ -265,6 +265,19 @@ function calc_plasma_inductance( return inv(temp2) end +""" + calc_surface_inductance(I_v::Matrix{ComplexF64})::Matrix{ComplexF64} + +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 + return inv(I_v) .* μ₀ * (2π)^2 +end + """ calc_permeability( plasma_inductance::Matrix{ComplexF64}, diff --git a/src/PerturbedEquilibrium/SingularCoupling.jl b/src/PerturbedEquilibrium/SingularCoupling.jl index d9cd1ebc..73431be9 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, @@ -173,17 +170,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) - - # Get ν on the vacuum theta grid (same ν used in the vacuum Fourier basis computation) - ν_vac = Vacuum.PlasmaGeometry(vac_input).ν + _, I_v, _, _ = Vacuum.compute_vacuum_response(vac_input, wall_settings; compute_Iv=true) + L_surf = calc_surface_inductance(I_v) - # 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) + # 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] @@ -432,99 +424,6 @@ function compute_current_density( return j_c end -""" - compute_surface_inductance_from_greens( - grri::Matrix{ComplexF64}, - grre::Matrix{ComplexF64}, - ffs_intr::ForceFreeStatesInternal, - nn::Int, - ν::Vector{Float64} - )::Matrix{ComplexF64} - -Compute surface inductance matrix from Green's functions at flux surface. - -Implements the GPEC `gpvacuum_flxsurf` algorithm. - -The Julia vacuum code uses SFL Fourier basis `cos(m*θ - n*ν)` in the column transform, -so the row DFT must apply the matching toroidal phase correction `exp(-i*n*ν)` before -the DFT (matching Fortran `gpvacuum_flxsurf`'s `EXP(-ifac*nn*dphi)` phase correction). - -## Arguments - - - `grri`: Interior Green's function [mtheta, mpert] - - `grre`: Exterior Green's function [mtheta, mpert] - - `ffs_intr`: ForceFreeStates internal state - - `nn`: Toroidal mode number - - `ν`: Toroidal angle offset on the vacuum theta grid [mtheta] - -## Returns - -Surface inductance matrix [mpert × mpert] -""" -@with_pool pool function compute_surface_inductance_from_greens( - grri::Matrix{ComplexF64}, - grre::Matrix{ComplexF64}, - ffs_intr::ForceFreeStatesInternal, - nn::Int, - ν::Vector{Float64} -)::Matrix{ComplexF64} - mpert = ffs_intr.mpert - mtheta = length(ν) - μ₀ = 4π * 1e-7 - - ft = FourierTransforms.FourierTransform(mtheta, mpert, ffs_intr.mlow) - - flux_matrix = zeros!(pool, ComplexF64, mpert, mpert) - current_matrix = zeros!(pool, ComplexF64, mpert, mpert) - - kax = zeros!(pool, ComplexF64, mtheta) - grri_surf = @view grri[1:mtheta, :] - grre_surf = @view grre[1:mtheta, :] - - # Toroidal phase correction: exp(-i*n*ν) - phase = cis.(-nn .* ν) - - for i in 1:mpert - flux_matrix[i, i] = 1.0 - - # Complex grri/e stores exp(i(mθ-nν)) projection, need conjugate for exp(-i(mθ-nν)) - kax .= conj.(grri_surf[:, i] .+ grre_surf[:, i]) ./ (μ₀ * (2π)^2) - - # Port of Fortran gpvacuum_flxsurf: apply toroidal phase, reverse theta, forward-DFT. - g_phased = kax .* phase - current_matrix[:, i] = ft(_reverse_theta(g_phased)) - end - - # Compute surface inductance: L_surf = flux * inv(current) = inv(current) - L_surf = zeros(ComplexF64, mpert, mpert) - - current_mag = maximum(abs.(current_matrix)) - - if current_mag < 1e-15 - @warn "Current matrix is all zeros! Cannot compute surface inductance." maxlog=1 - for i in 1:mpert - L_surf[i, i] = μ₀ * 1e-6 - end - else - try - regularization = 1e-12 * current_mag - current_reg = current_matrix + regularization * I - - L_surf = flux_matrix * inv(current_reg) - - # Hermitianize (matches Fortran: temp1 = 0.5*(temp1 + CONJG(TRANSPOSE(temp1)))) - L_surf = 0.5 * (L_surf + L_surf') - catch e - @warn "Surface inductance inversion failed: $e" maxlog=1 - for i in 1:mpert - L_surf[i, i] = μ₀ * 1e-6 - end - end - end - - return L_surf -end - """ compute_surface_area( equil::Equilibrium.PlasmaEquilibrium, diff --git a/src/Vacuum/DataTypes.jl b/src/Vacuum/DataTypes.jl index 06ee4a6a..fc8bd0ed 100644 --- a/src/Vacuum/DataTypes.jl +++ b/src/Vacuum/DataTypes.jl @@ -73,7 +73,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..c290246e 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. if source isa PlasmaGeometry grad_greenfunction_block .*= -1 end diff --git a/src/Vacuum/Vacuum.jl b/src/Vacuum/Vacuum.jl index b060f59f..de653e07 100644 --- a/src/Vacuum/Vacuum.jl +++ b/src/Vacuum/Vacuum.jl @@ -24,50 +24,45 @@ 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 """ - _compute_vacuum_response_2d!(vac_data, inputs::VacuumInput, wall_settings::WallShapeSettings) + _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 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 `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) +@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_Iv && (vac_data.I_v .= 0) # Form the plasma and wall geometries plasma_surf = PlasmaGeometry(inputs) @@ -75,13 +70,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 I_v 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 @@ -89,10 +82,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) @@ -111,29 +101,44 @@ 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 - - # Solve exterior first, overwriting grad_green to save memory since we already have the interior kernel - ldiv!(lu!(grad_green), grre) - - # 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 + if compute_Iv + # 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 + + # 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 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 + end + ldiv!(lu!(grad_green_interior), grri) + + # Surface-current matrix, Park 2007 eq. 21b: μ₀I^v = χ^(vi) - χ^(vo) = grri - grre + # 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 = grre[1:num_points_surf, :] .- grri[1:num_points_surf, :] + mul!(I_v_block, ft.basis, g_sum) + I_v_block ./= num_points_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, :])) 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_Iv && _warn_and_symmetrize!(vac_data.I_v, "Iᵛ") # Populate coordinate arrays @views begin @@ -147,7 +152,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_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 @@ -165,14 +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; `grri`/`grre` are returned zeroed and are not yet supported in the 3D path. +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) +@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_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) @@ -245,11 +254,9 @@ interior variant `-D + 2I` for the interior columns, then scatter back into the end end - _symmetrize_vacuum_energy!(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) + # Remove any non-Hermitian residual from Hermitian matrices due to discretization + _warn_and_symmetrize!(vac_data.wv, "Wᵛ") + compute_Iv && _warn_and_symmetrize!(vac_data.I_v, "Iᵛ") # Populate coordinate arrays vac_data.plasma_pts .= plasma_surf.r @@ -257,10 +264,10 @@ interior variant `-D + 2I` for the interior columns, then scatter back into the end """ - compute_vacuum_response(inputs::VacuumInput, wall_settings::WallShapeSettings) + compute_vacuum_response(inputs::VacuumInput, wall_settings::WallShapeSettings; compute_Iv=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. @@ -268,28 +275,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). + - `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) +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) vac = ( wv=zeros(ComplexF64, num_modes, num_modes), - grri=zeros(ComplexF64, 2 * num_points, num_modes), - grre=zeros(ComplexF64, 2 * num_points, 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_vacuum_response!(vac, inputs, wall_settings; compute_Iv) - return vac.wv, vac.grri, vac.grre, vac.plasma_pts, vac.wall_pts + return vac.wv, vac.I_v, 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_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 @@ -299,19 +305,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 + - `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) +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_vacuum_response_2d!(vac_data, inputs, wall_settings; compute_Iv) else - _compute_vacuum_response_3d!(vac_data, inputs, wall_settings) + _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 7e8b51b4..1e84546b 100644 --- a/test/runtests_vacuum.jl +++ b/test/runtests_vacuum.jl @@ -368,7 +368,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, 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) @@ -376,10 +376,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(I_v) == (num_modes, num_modes) + @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) @@ -388,15 +386,29 @@ @test isapprox(wv, wv', rtol=1e-12) end + @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_Iv=true) + + num_modes = length(inputs.m_modes) * length(inputs.n_modes) + @test size(wv) == (num_modes, num_modes) + @test all(isfinite, wv) + @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, grri, grre, 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(grri) == (2 * numpoints, num_modes) + @test all(iszero, I_v) @test all(isfinite, plasma_pts) @test all(isfinite, wall_pts) # plasma_pts layout: col1=R, col2=0, col3=Z @@ -408,19 +420,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, I_v, 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(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, grri, grre, 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(grri) == (32, 2) # 2*mtheta, num_modes=2 + @test size(I_v) == (2, 2) @test size(plasma_pts) == (16, 3) end @@ -429,17 +441,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, 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), grri=zeros(ComplexF64, 2 * numpoints, num_modes), grre=zeros(ComplexF64, 2 * numpoints, num_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_vacuum_response!(vac, inputs, wall_settings; compute_Iv=true) @test vac.wv ≈ wv - @test vac.grri ≈ grri - @test vac.grre ≈ grre + @test vac.I_v ≈ I_v @test vac.plasma_pts ≈ pp @test vac.wall_pts ≈ wp end @@ -581,17 +593,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, 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(grri) == (2 * numpoints, num_modes) - @test size(grre) == (2 * numpoints, num_modes) - @test all(isfinite, grri) - @test all(isfinite, grre) + @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) @@ -604,7 +614,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, 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) @@ -612,10 +622,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(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) @@ -625,12 +633,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, 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(grri) == (2 * numpoints, 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) @@ -676,13 +684,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)