diff --git a/.claude/agent-memory/fortran-physics-reviewer/fortran_correspondence_map.md b/.claude/agent-memory/fortran-physics-reviewer/fortran_correspondence_map.md index 4ba8494e8..452de38d1 100644 --- a/.claude/agent-memory/fortran-physics-reviewer/fortran_correspondence_map.md +++ b/.claude/agent-memory/fortran-physics-reviewer/fortran_correspondence_map.md @@ -15,7 +15,7 @@ relative to this repository. - `gpec/gpeq.f` (gpeq_sol, gpeq_contra, gpeq_surface, gpeq_normal) -> `src/PerturbedEquilibrium/FieldReconstruction.jl` + `src/PerturbedEquilibrium/ResponseMatrices.jl` - `gpec/gpresp.f` (gpresp_pinduct, gpresp_sinduct, gpresp_permeab) -> `src/PerturbedEquilibrium/ResponseMatrices.jl` - `gpec/gpout.f` (gpout_singcoup, gpout_xbnormal) -> `src/PerturbedEquilibrium/SingularCoupling.jl` + `src/PerturbedEquilibrium/FieldReconstruction.jl` -- `gpec/gpvacuum.f` (gpvacuum_flxsurf) -> `src/PerturbedEquilibrium/SingularCoupling.jl` (compute_surface_inductance_from_greens) +- `gpec/gpvacuum.f` (gpvacuum_flxsurf) -> `src/Vacuum/Vacuum.jl` (compute_Iv branch of _compute_vacuum_response_2d!) + `src/PerturbedEquilibrium/ResponseMatrices.jl` (calc_surface_inductance) ## KineticForces (NTV, Fortran `pentrc/`) — see kinetic_ntv_map.md for the audit checklist - `pentrc/torque.F90` -> `src/KineticForces/Torque.jl` (tpsi! single-surface torque) diff --git a/docs/development/architecture.md b/docs/development/architecture.md index 2eb4a2422..ef756c9c0 100644 --- a/docs/development/architecture.md +++ b/docs/development/architecture.md @@ -46,13 +46,15 @@ GPEC consists of **seven main modules** organized in `src/`: 4. **Vacuum** (`src/Vacuum/`) - Vacuum field calculations and Green's functions - Computes vacuum response matrices for ideal MHD analysis - - Calculates both **interior** (grri) and **exterior** (grre) Green's functions + - Solves the exterior boundary-integral system for the vacuum energy matrix `wv`, and optionally + (`compute_Iv=true`) the interior system as well to build the surface-current matrix `I_v` + (Park 2007 eq. 21b). The interior/exterior Green's functions themselves are internal scratch. - Main functions: - - `compute_vacuum_response()` - Pure Julia implementation + - `compute_vacuum_response()` / `compute_vacuum_response!()` - allocating and in-place entry points - Key files: - - `VacuumStructs.jl` - Data structures - - `VacuumInternals.jl` - Core algorithms - - `VacuumFromEquilibrium.jl` - Integration with equilibrium data + - `DataTypes.jl` - Data structures (`VacuumInput`, `PlasmaGeometry`, `WallGeometry`) + - `Kernel2D.jl` / `Kernel3D.jl` - Single-/double-layer kernel assembly + - `Field.jl` - Vacuum field and potential evaluation off the surface - Status: **Pure Julia implementation complete and available** 5. **ForceFreeStates** (`src/ForceFreeStates/`) - Ideal MHD stability analysis (DCON-style) @@ -130,8 +132,7 @@ The complete GPEC analysis pipeline: 2. **Vacuum Response**: - Initialize plasma and wall surfaces from equilibrium - - Compute vacuum response matrices (wv, grri, grre) - - Calculate both interior and exterior Green's functions + - Compute the vacuum energy matrix `wv` (and `I_v` when `compute_Iv=true`) - Pure Julia implementation 3. **Stability Analysis** (ForceFreeStates): @@ -169,11 +170,9 @@ The complete GPEC analysis pipeline: ### Stability - `SingType` - Singular surface data including: - - Rational surface location (ψ, q = m/n) - - Δ' (tearing stability parameter) - - Eigenmode structure at singular surface - - Green's functions (grri, grre) at interior singular surfaces - - Surface inductance + - Rational surface location (ψ, ρ, q = m/n, dq/dψ) + - Δ' (tearing stability parameter) — **stub**; the valid Δ' is `ForceFreeStatesInternal.delta_prime_matrix` + - Asymptotic solution bases at the inner-layer boundaries ### Perturbed Equilibrium - `PerturbedEquilibriumControl` - User-facing TOML configuration parameters diff --git a/docs/src/vacuum.md b/docs/src/vacuum.md index 40b7c08fb..ddbbdcdf4 100644 --- a/docs/src/vacuum.md +++ b/docs/src/vacuum.md @@ -73,7 +73,7 @@ wall_settings = GeneralizedPerturbedEquilibrium.Vacuum.WallShapeSettings( equal_arc_wall = true # Use equal arc length spacing ) -# Compute vacuum response; returns a VacuumResponse with wv, grri, grre, plasma_pts, wall_pts +# Compute vacuum response; returns a VacuumResponse with wv, I_v, plasma_pts, wall_pts vac = GeneralizedPerturbedEquilibrium.Vacuum.compute_vacuum_response(inputs, wall_settings) ``` diff --git a/docs/src/workflow.md b/docs/src/workflow.md index 6929ddc3c..6ae9eef2a 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) diff --git a/src/ForceFreeStates/ForceFreeStatesStructs.jl b/src/ForceFreeStates/ForceFreeStatesStructs.jl index c495619f1..00eaf0ba7 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 @@ -540,26 +536,26 @@ and a small set of temporary matrices and factors used to compute singular-layer - `ising_start::Int` - Index of the starting singular surface to be crossed during integration. - # Initialization parameters - - `psimax::Float64` - Maximum psi value for which the integrator is allowed to run in next integration region. - `needs_crossing::Bool` - Flag indicating whether a rational surface needs to be crossed after the current integration region. - `nzero::Int` - Count of detected zero crossings (used for diagnostics). - # Saved data throughout integration - - `new::Bool` - Flag indicating whether a new `unorm0` should be computed after a fixup. -# Total ODE solver steps taken (all steps, not just saved ones) + # Initialization parameters - `unorm::Vector{Float64}` - Current norms of the solution vectors (length `numpert_total`). - `unorm0::Vector{Float64}` - Reference/initial norms of the solution vectors (length `numpert_total`). + # Saved data throughout integration + - `ifix::Int` - Number of normalization operations performed (index into normalization arrays). +# Total ODE solver steps taken (all steps, not just saved ones) + - `index::Array{Int,2}` - Index matrix used for sorting solution norms with shape `(numpert_total, numunorms_init)`. - `sing_flag::Vector{Bool}` - Boolean flags indicating which stored normalizations correspond to singular solutions # Edge dW scan state and results (disabled sentinel when psiedge >= psilim, i.e. no edge scan) @@ -567,8 +563,7 @@ and a small set of temporary matrices and factors used to compute singular-layer - `zeroed_idx::Vector{Vector{Int}}` - For each ideal rational surface jump, a vector of indices of solutions that were zeroed. # Data for integrator - - `fixfac::Array{ComplexF64,3}` - Fix-up factors for Gaussian reduction with shape # Initialization parameters - `(numpert_total, numpert_total, numunorms_init)`. + - `fixfac::Array{ComplexF64,3}` - Fix-up factors for Gaussian reduction with shape `(numpert_total, numpert_total, numunorms_init)`. - `fixstep::Vector{Int64}` - Step indices (psi step positions) at which normalization/fixups were performed (length `numunorms_init`). """ diff --git a/src/PerturbedEquilibrium/Response.jl b/src/PerturbedEquilibrium/Response.jl index e501bad29..d9a985204 100644 --- a/src/PerturbedEquilibrium/Response.jl +++ b/src/PerturbedEquilibrium/Response.jl @@ -36,16 +36,9 @@ function compute_plasma_response!( # Plasma inductance Lambda (wt0 formula, Fortran resp_induct_flag=TRUE default) plasma_inductance = calc_plasma_inductance(wt0, 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, mthvac, 1, ffs_intr.mlow:ffs_intr.mhigh, [nn]) - wall_nowall = Vacuum.WallShapeSettings(; shape="nowall") - vac_2d = Vacuum.compute_vacuum_response(vac_input_2d, wall_nowall) - grri_2d = Matrix{ComplexF64}(vac_2d.grri) - grre_2d = Matrix{ComplexF64}(vac_2d.grre) - ν_vac = Vacuum.PlasmaGeometry(vac_input_2d).ν - surface_inductance = compute_surface_inductance_from_greens(grri_2d, grre_2d, ffs_intr, nn, ν_vac) + surface_inductance = calc_surface_inductance(equil, ffs_intr.psilim, mthvac, ffs_intr.mlow:ffs_intr.mhigh, nn) 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 791ed9c57..4a8ec46b8 100644 --- a/src/PerturbedEquilibrium/ResponseMatrices.jl +++ b/src/PerturbedEquilibrium/ResponseMatrices.jl @@ -111,10 +111,12 @@ where ξ_ψ is the radial displacement eigenfunction. ## Arguments - `boundary_data`: Output from extract_boundary_displacements() + + ξ_psi_boundary: Boundary displacement [numpert_total, numpert_total] + dPsi_drho: Flux surface spacing at boundary (scalar) + q_boundary: Safety factor at boundary (scalar) + psi_boundary: Normalized flux at boundary (scalar) + - `intr`: ForceFreeStates internal state with mode arrays (mlow, mhigh, nlow, etc.) ## Returns @@ -262,6 +264,48 @@ function calc_plasma_inductance( return inv(temp2) end +""" + calc_surface_inductance( + equil::Equilibrium.PlasmaEquilibrium, + ψ::Float64, + mtheta::Int, + m_modes::AbstractUnitRange{Int}, + nn::Int + )::Matrix{ComplexF64} + +Surface inductance L at the flux surface ψ, from the vacuum surface-current matrix, solved +against a no-wall vacuum. + +Solves the 2D vacuum problem with `compute_Iv=true` and inverts. The driving flux harmonics are +unit columns (`Φ = 𝕀` by construction, matching Fortran `gpvacuum_flxsurf`'s unit driving +harmonics), so `Φ = L·I^v` (Park 2007, eq. 7 and following text) gives `L = μ₀(2π)²·I_v⁻¹`. + +## Arguments + + - `equil`: Plasma equilibrium + - `ψ`: Normalized poloidal flux of the surface + - `mtheta`: Number of vacuum poloidal grid points + - `m_modes`: Poloidal mode range mlow:mhigh + - `nn`: Toroidal mode number + +## Returns + + - Surface inductance matrix in henries [mpert × mpert] +""" +function calc_surface_inductance( + equil::Equilibrium.PlasmaEquilibrium, + ψ::Float64, + mtheta::Int, + m_modes::AbstractUnitRange{Int}, + nn::Int +)::Matrix{ComplexF64} + vac_input = Vacuum.VacuumInput(equil, ψ, mtheta, 1, m_modes, [nn]) + wall_settings = Vacuum.WallShapeSettings(; shape="nowall") + I_v = Vacuum.compute_vacuum_response(vac_input, wall_settings; compute_Iv=true).I_v + μ₀ = 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 417474dc7..8a2c12bd5 100644 --- a/src/PerturbedEquilibrium/SingularCoupling.jl +++ b/src/PerturbedEquilibrium/SingularCoupling.jl @@ -51,9 +51,6 @@ function _hermite_cubic_deriv(u_a, u_b, du_a, du_b, psi_a, psi_b, psi) return @. d00 * u_a + d10 * du_a + d01 * u_b + d11 * du_b end -# Reflect a periodic theta-space vector θ → -θ (the theta reversal in gpvacuum_flxsurf). -_reverse_theta(v::AbstractVector) = circshift(reverse(v), 1) - """ _chord_solution_at(psi, resnum, odet, nstep) -> (u, du) @@ -137,7 +134,7 @@ function _solution_at( # Same-side candidate nodes around the bracket, trimmed to the 4 nearest psi. side = sign(psi - psi_surf) - idxs = [j for j in max(1, il - 3):min(nstep, ir + 3) if sign(odet.psi_store[j] - psi_surf) == side] + idxs = [j for j in max(1, il-3):min(nstep, ir+3) if sign(odet.psi_store[j] - psi_surf) == side] while length(idxs) > 4 abs(odet.psi_store[idxs[1]] - psi) > abs(odet.psi_store[idxs[end]] - psi) ? popfirst!(idxs) : pop!(idxs) end @@ -280,7 +277,6 @@ function compute_singular_coupling_metrics!( chi1 = 2π * equil.psio twopi = 2π mtheta = mthvac - wall_settings = Vacuum.WallShapeSettings(; shape="nowall") # Phase 1: Collect all resonant (surface, n) pairs in psi order resonant_pairs = Tuple{Int,Int}[] @@ -362,17 +358,10 @@ function compute_singular_coupling_metrics!( continue end - # Compute Green's functions at this surface for this n (once per pair) - vac_input = Vacuum.VacuumInput(equil, sing_surf.psifac, mtheta, 1, mlow:mhigh, [nn]) - vac = Vacuum.compute_vacuum_response(vac_input, wall_settings) - grri = Matrix{ComplexF64}(vac.grri) - grre = Matrix{ComplexF64}(vac.grre) - - # Get ν on the vacuum theta grid (same ν used in the vacuum Fourier basis computation) - ν_vac = Vacuum.PlasmaGeometry(vac_input).ν + # Surface inductance at this surface for this n (once per pair) + L_surf = calc_surface_inductance(equil, sing_surf.psifac, mtheta, mlow:mhigh, nn) - # 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) + # Only 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] @@ -614,99 +603,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 63dec0b9e..da0b8167f 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)], @@ -152,13 +153,15 @@ boundary-integral solve produces along the way. ## Fields - `wv::Matrix{ComplexF64}`: Vacuum energy matrix Wᵛ (`num_modes × num_modes`), block-diagonal in n for 2D - - `grri`, `grre::Matrix{ComplexF64}`: Interior/exterior Green's functions (`2·num_points × num_modes`), zeroed on every 3D path + - `I_v::Matrix{ComplexF64}`: Vacuum surface-current matrix Iᵛ (`num_modes × num_modes`), left zeroed + unless `compute_vacuum_response` is called with `compute_Iv=true` (2D only). Stored without the + `μ₀`/`4π²` normalization: the physical surface inductance is `μ₀(2π)²·I_v⁻¹` + (see `PerturbedEquilibrium.calc_surface_inductance`). - `plasma_pts`, `wall_pts::Matrix{Float64}`: Cartesian surface coordinates (`num_points × 3`) """ struct VacuumResponse wv::Matrix{ComplexF64} - grri::Matrix{ComplexF64} - grre::Matrix{ComplexF64} + I_v::Matrix{ComplexF64} plasma_pts::Matrix{Float64} wall_pts::Matrix{Float64} end @@ -173,8 +176,7 @@ function VacuumResponse(inputs::VacuumInput) num_modes = length(inputs.m_modes) * length(inputs.n_modes) return VacuumResponse( zeros(ComplexF64, num_modes, num_modes), - zeros(ComplexF64, 2 * num_points, num_modes), - zeros(ComplexF64, 2 * num_points, num_modes), + zeros(ComplexF64, num_modes, num_modes), zeros(num_points, 3), zeros(num_points, 3) ) diff --git a/src/Vacuum/Vacuum.jl b/src/Vacuum/Vacuum.jl index 30af03b9a..3fd13eb22 100644 --- a/src/Vacuum/Vacuum.jl +++ b/src/Vacuum/Vacuum.jl @@ -24,52 +24,49 @@ 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. +Replace `mat` by its Hermitian part in place, warning first if the anti-Hermitian residual exceeds +`_HERMITICITY_WARN_TOL`. -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. +The matrices passed here are Hermitian in exact arithmetic — δW_v = ξ†Wᵛξ is a real energy, and Iᵛ +is the inverse of the Hermitian surface inductance up to a real scalar — so any anti-Hermitian part +is a discretization artifact that vanishes as the vacuum grid is refined. """ -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::VacuumResponse, inputs::VacuumInput, wall_settings::WallShapeSettings) + _compute_vacuum_response_2d!(vac_data::VacuumResponse, 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::VacuumResponse, inputs::VacuumInput, wall_settings::WallShapeSettings) +@with_pool pool function _compute_vacuum_response_2d!(vac_data::VacuumResponse, inputs::VacuumInput, wall_settings::WallShapeSettings; compute_Iv::Bool=false) mpert = length(inputs.m_modes) + mlow = inputs.m_modes[1] num_points_surf = inputs.mtheta fill!(vac_data.wv, 0) - fill!(vac_data.grri, 0) - fill!(vac_data.grre, 0) + fill!(vac_data.I_v, 0) # Form the plasma and wall geometries plasma_surf = PlasmaGeometry(inputs) @@ -77,13 +74,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 @@ -91,10 +86,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) @@ -113,29 +105,46 @@ 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); the solve gives + # grre = -(2π)²χ^(vo), the vacuum-outside potential. 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] + g_diff = @view grri[1:num_points_surf, :] + g_diff .= @view(grre[1:num_points_surf, :]) .- g_diff + mul!(I_v_block, ft.basis, g_diff) + conj!(I_v_block) # Flip θ_VAC → -θ_VAC to get I^v in GPEC's CCW-θ frame. + 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 @@ -149,7 +158,7 @@ of the Green's functions. end """ - _compute_vacuum_response_3d!(vac_data::VacuumResponse, inputs::VacuumInput, wall_settings::WallShapeSettings) + _compute_vacuum_response_3d!(vac_data::VacuumResponse, 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 @@ -167,14 +176,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::VacuumResponse, inputs::VacuumInput, wall_settings::WallShapeSettings) +@with_pool pool function _compute_vacuum_response_3d!(vac_data::VacuumResponse, 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) @@ -247,11 +260,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 @@ -259,30 +270,31 @@ interior variant `-D + 2I` for the interior columns, then scatter back into the end """ - compute_vacuum_response(inputs::VacuumInput, wall_settings::WallShapeSettings) -> VacuumResponse + compute_vacuum_response(inputs::VacuumInput, wall_settings::WallShapeSettings; compute_Iv=false) -> VacuumResponse Compute the vacuum response for the given inputs. Allocating wrapper around [`compute_vacuum_response!`](@ref); pass a preallocated [`VacuumResponse`](@ref) to that method -instead when reusing storage across calls. +instead when reusing storage across calls. Pass `compute_Iv=true` to additionally populate the +surface-current matrix `I_v` (2D only). """ -function compute_vacuum_response(inputs::VacuumInput, wall_settings::WallShapeSettings) +function compute_vacuum_response(inputs::VacuumInput, wall_settings::WallShapeSettings; compute_Iv::Bool=false) vac = VacuumResponse(inputs) - compute_vacuum_response!(vac, inputs, wall_settings) + compute_vacuum_response!(vac, inputs, wall_settings; compute_Iv) return vac end """ - compute_vacuum_response!(vac_data::VacuumResponse, inputs::VacuumInput, wall_settings::WallShapeSettings) + compute_vacuum_response!(vac_data::VacuumResponse, inputs::VacuumInput, wall_settings::WallShapeSettings; compute_Iv=false) In-place variant that populates the arrays of an existing [`VacuumResponse`](@ref). Dispatches on dimensionality only: 2D (`inputs.nzeta == 1`) routes to [`_compute_vacuum_response_2d!`], 3D to [`_compute_vacuum_response_3d!`]. """ -function compute_vacuum_response!(vac_data::VacuumResponse, inputs::VacuumInput, wall_settings::WallShapeSettings) +function compute_vacuum_response!(vac_data::VacuumResponse, 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_resist_eval.jl b/test/runtests_resist_eval.jl index 315645494..58aee3d9b 100644 --- a/test/runtests_resist_eval.jl +++ b/test/runtests_resist_eval.jl @@ -70,7 +70,6 @@ dq = deriv_view(equil.profiles.q_spline, 1) s1 = SingType(; psifac=0.3, rho=sqrt(0.3), m=[2], n=[1], q=2.0, q1=dq(0.3), - grri=zeros(Float64, 0, 0), grre=zeros(Float64, 0, 0), delta_prime=ComplexF64[], delta_prime_col=zeros(ComplexF64, 0, 0), ua_left=zeros(ComplexF64, 0, 0, 0), @@ -78,7 +77,6 @@ psi_ua_left=0.0, psi_ua_right=0.0) s2 = SingType(; psifac=0.7, rho=sqrt(0.7), m=[3], n=[1], q=3.0, q1=dq(0.7), - grri=zeros(Float64, 0, 0), grre=zeros(Float64, 0, 0), delta_prime=ComplexF64[], delta_prime_col=zeros(ComplexF64, 0, 0), ua_left=zeros(ComplexF64, 0, 0, 0), @@ -113,7 +111,6 @@ dq = deriv_view(equil.profiles.q_spline, 1) s1 = SingType(; psifac=0.3, rho=sqrt(0.3), m=[2], n=[1], q=2.0, q1=dq(0.3), - grri=zeros(Float64, 0, 0), grre=zeros(Float64, 0, 0), delta_prime=ComplexF64[], delta_prime_col=zeros(ComplexF64, 0, 0), ua_left=zeros(ComplexF64, 0, 0, 0), @@ -155,7 +152,6 @@ s_unpop = SingType(; psifac=0.5, rho=sqrt(0.5), m=[2], n=[1], q=2.0, q1=1.0, - grri=zeros(Float64, 0, 0), grre=zeros(Float64, 0, 0), delta_prime=ComplexF64[], delta_prime_col=zeros(ComplexF64, 0, 0), ua_left=zeros(ComplexF64, 0, 0, 0), @@ -178,7 +174,6 @@ dq = deriv_view(equil.profiles.q_spline, 1) s1 = SingType(; psifac=0.3, rho=sqrt(0.3), m=[2], n=[1], q=2.0, q1=dq(0.3), - grri=zeros(Float64, 0, 0), grre=zeros(Float64, 0, 0), delta_prime=ComplexF64[], delta_prime_col=zeros(ComplexF64, 0, 0), ua_left=zeros(ComplexF64, 0, 0, 0), diff --git a/test/runtests_slayer_inputs.jl b/test/runtests_slayer_inputs.jl index 5e2c5a34f..c0b7d2d7d 100644 --- a/test/runtests_slayer_inputs.jl +++ b/test/runtests_slayer_inputs.jl @@ -9,25 +9,24 @@ # Load the Solovev analytic equilibrium shipped with the examples. # This exercise gets run once for all LayerInputs tests. dir_path = joinpath(dirname(@__DIR__), "examples", "Solovev_ideal_example") - inputs = TOML.parsefile(joinpath(dir_path, "gpec.toml")) - eq_cfg = Equilibrium.EquilibriumConfig(inputs["Equilibrium"], dir_path) - sol_cfg = Equilibrium.SolovevConfig(inputs["SOL_INPUT"]) - equil = Equilibrium.setup_equilibrium(eq_cfg, sol_cfg) + inputs = TOML.parsefile(joinpath(dir_path, "gpec.toml")) + eq_cfg = Equilibrium.EquilibriumConfig(inputs["Equilibrium"], dir_path) + sol_cfg = Equilibrium.SolovevConfig(inputs["SOL_INPUT"]) + equil = Equilibrium.setup_equilibrium(eq_cfg, sol_cfg) # Synthetic profiles (simple linear-in-ψ temperature decrease) - psi_pts = collect(0.0:0.1:1.0) + psi_pts = collect(0.0:0.1:1.0) profiles = KineticProfiles(; psi=psi_pts, - n_e=fill(5.0e19, length(psi_pts)), - T_e=1000.0 .* (1.0 .- 0.7 .* psi_pts), - T_i=1000.0 .* (1.0 .- 0.6 .* psi_pts), - omega=fill(0.0, length(psi_pts)), - omega_e=fill(1.0e4, length(psi_pts)), - omega_i=fill(5.0e3, length(psi_pts))) + n_e=fill(5.0e19, length(psi_pts)), + T_e=1000.0 .* (1.0 .- 0.7 .* psi_pts), + T_i=1000.0 .* (1.0 .- 0.6 .* psi_pts), + omega=fill(0.0, length(psi_pts)), + omega_e=fill(1.0e4, length(psi_pts)), + omega_i=fill(5.0e3, length(psi_pts))) # Helper to build a minimal SingType without touching unused fields _mk_sing(; psi, q, q1, m, n, delta_prime=-10.0+0im) = SingType( psifac=psi, rho=sqrt(psi), m=[m], n=[n], q=q, q1=q1, - grri=zeros(Float64, 0, 0), grre=zeros(Float64, 0, 0), delta_prime=ComplexF64[delta_prime], delta_prime_col=zeros(ComplexF64, 0, 0), ua_left=zeros(ComplexF64, 0, 0, 0), @@ -58,15 +57,15 @@ # Near ψ=0 and ψ=1, the function falls back to one-sided FD and # should still produce a finite positive number (minor radius is # still increasing). - d_near_axis = surface_da_dpsi(equil, 1e-6) - d_near_edge = surface_da_dpsi(equil, 1.0 - 1e-6) + d_near_axis = surface_da_dpsi(equil, 1e-6) + d_near_edge = surface_da_dpsi(equil, 1.0 - 1e-6) @test isfinite(d_near_axis) && d_near_axis > 0 @test isfinite(d_near_edge) && d_near_edge > 0 end @testset "build_slayer_inputs: returns correct per-surface data" begin sings = [_mk_sing(psi=0.3, q=2.0, q1=1.5, m=2, n=1), - _mk_sing(psi=0.6, q=3.0, q1=2.5, m=3, n=1)] + _mk_sing(psi=0.6, q=3.0, q1=2.5, m=3, n=1)] # dr_val=0.0 bypasses the build_slayer_inputs requirement that sing.restype be # pre-populated by ForceFreeStates.resist_eval_all! — the test sings here are # minimal stubs without restype, so we supply dr_val explicitly. @@ -74,7 +73,7 @@ # rather than being recomputed from n_e/T_e/T_i gradients — required for the Q_e == # -tauk·omega_e(ψ) identity check below. sl = build_slayer_inputs(equil, sings, profiles; bt=2.0, dr_val=0.0, - compute_omega_star=false) + compute_omega_star=false) @test length(sl) == 2 @test sl[1] isa SLAYERParameters @@ -112,19 +111,19 @@ # Scalar (dr_val=0.0 bypasses the sing.restype requirement; see comment above) sl_s = build_slayer_inputs(equil, sings, profiles; - bt=2.0, chi_perp=2.0, chi_tor=1.5, dr_val=0.0) + bt=2.0, chi_perp=2.0, chi_tor=1.5, dr_val=0.0) # Callable with matching value chi_p(psi) = 2.0 + 0.0*psi chi_t(psi) = 1.5 + 0.0*psi sl_c = build_slayer_inputs(equil, sings, profiles; - bt=2.0, chi_perp=chi_p, chi_tor=chi_t, dr_val=0.0) + bt=2.0, chi_perp=chi_p, chi_tor=chi_t, dr_val=0.0) @test sl_s[1].P_perp ≈ sl_c[1].P_perp - @test sl_s[1].P_tor ≈ sl_c[1].P_tor + @test sl_s[1].P_tor ≈ sl_c[1].P_tor # Callable with ψ-dependence changes the result chi_p_var(psi) = 1.0 + 10.0 * psi # χ⊥(0.5) = 6.0 > 2.0 sl_var = build_slayer_inputs(equil, sings, profiles; - bt=2.0, chi_perp=chi_p_var, chi_tor=1.5, dr_val=0.0) + bt=2.0, chi_perp=chi_p_var, chi_tor=1.5, dr_val=0.0) # P_perp = τ_r · χ⊥ / r² grows with χ⊥, so the varying-χ case at # ψ=0.5 (χ⊥=6) gives a *larger* P_perp than the scalar χ⊥=2. @test sl_var[1].P_perp > sl_s[1].P_perp @@ -136,17 +135,17 @@ # dc_type=:none and dr_val=0.0 → dc_tmp = 0 regardless of dr_val sl_none = build_slayer_inputs(equil, sings, profiles; - bt=2.0, dc_type=:none, dr_val=0.0) + bt=2.0, dc_type=:none, dr_val=0.0) @test sl_none[1].dc_tmp == 0.0 # dc_type=:rfitzp with dr_val = 0 still gives zero sl_rf0 = build_slayer_inputs(equil, sings, profiles; - bt=2.0, dc_type=:rfitzp, dr_val=0.0) + bt=2.0, dc_type=:rfitzp, dr_val=0.0) @test sl_rf0[1].dc_tmp == 0.0 # dc_type=:rfitzp with dr_val > 0 → nonzero negative offset sl_rf = build_slayer_inputs(equil, sings, profiles; - bt=2.0, dc_type=:rfitzp, dr_val=0.01) + bt=2.0, dc_type=:rfitzp, dr_val=0.01) @test sl_rf[1].dc_tmp < 0 @test isfinite(sl_rf[1].dc_tmp) end diff --git a/test/runtests_vacuum.jl b/test/runtests_vacuum.jl index 630c60436..19f80b94d 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,28 @@ @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 +419,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,37 +440,32 @@ # in-place entry populates caller-owned storage identically. for wall_settings in (WallShapeSettings(shape="nowall"), WallShapeSettings(shape="conformal", a=0.5)) inputs = _make_inputs() - ref = compute_vacuum_response(inputs, wall_settings) + ref = compute_vacuum_response(inputs, wall_settings; compute_Iv=true) vac = VacuumResponse(inputs) - compute_vacuum_response!(vac, inputs, wall_settings) + compute_vacuum_response!(vac, inputs, wall_settings; compute_Iv=true) @test vac.wv ≈ ref.wv - @test vac.grri ≈ ref.grri - @test vac.grre ≈ ref.grre + @test vac.I_v ≈ ref.I_v @test vac.plasma_pts ≈ ref.plasma_pts @test vac.wall_pts ≈ ref.wall_pts end end @testset "in-place compute_vacuum_response! clears a reused buffer" begin - # The nowall path writes only the plasma rows of grri/grre, so a buffer left - # over from a wall run must not leak its wall rows into the next result. + # A buffer left over from an earlier run must not leak into the next result: + # I_v in particular is only written when compute_Iv=true. inputs = _make_inputs() - wall_rows = (inputs.mtheta+1):(2*inputs.mtheta) vac = VacuumResponse(inputs) - compute_vacuum_response!(vac, inputs, WallShapeSettings(; shape="conformal", a=0.5)) - @test any(!iszero, view(vac.grre, wall_rows, :)) + compute_vacuum_response!(vac, inputs, WallShapeSettings(; shape="conformal", a=0.5); compute_Iv=true) + @test !all(iszero, vac.I_v) compute_vacuum_response!(vac, inputs, WallShapeSettings(; shape="nowall")) fresh = compute_vacuum_response(inputs, WallShapeSettings(; shape="nowall")) - @test all(iszero, view(vac.grri, wall_rows, :)) - @test all(iszero, view(vac.grre, wall_rows, :)) + @test all(iszero, vac.I_v) @test vac.wv ≈ fresh.wv - @test vac.grri ≈ fresh.grri - @test vac.grre ≈ fresh.grre end end @@ -481,6 +487,19 @@ r_out, z_out, _ = extract_plasma_surface_at_psi(pe, 0.8) @test extent(r_out, z_out) > extent(r_in, z_in) # minor radius grows outward in ψ end + + @testset "calc_surface_inductance" begin + # Same Solovev recipe as above; exercises the PerturbedEquilibrium helper end-to-end + eq_config = Equilibrium.EquilibriumConfig(; eq_type="sol", eq_filename="unused", jac_type="pest", grid_type="ldp", psilow=1e-4, psihigh=0.99999, mpsi=64, mtheta=128) + sol_config = Equilibrium.SolovevConfig(64, 64, 64, 1.6, 0.33, 1.0, 1.9, 1.0, 1.0, 1.0) + pe = Equilibrium.equilibrium_solver(Equilibrium.sol_run(eq_config, sol_config)) + + m_modes = -2:5 + L = PerturbedEquilibrium.calc_surface_inductance(pe, 0.9, 64, m_modes, 1) + @test size(L) == (length(m_modes), length(m_modes)) + @test all(isfinite, L) + @test isapprox(L, L', rtol=1e-8) # Hermitian inductance + end end # 3D vacuum: nzeta > 1, full (m,n) coupling, PlasmaGeometry3D, WallGeometry3D @@ -598,17 +617,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) @@ -621,7 +638,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) @@ -629,10 +646,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) @@ -642,12 +657,11 @@ @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)