diff --git a/src/Vacuum/DataTypes.jl b/src/Vacuum/DataTypes.jl index de1acfe02..1fe4bc609 100644 --- a/src/Vacuum/DataTypes.jl +++ b/src/Vacuum/DataTypes.jl @@ -154,7 +154,7 @@ boundary-integral solve produces along the way. - `wv::Matrix{ComplexF64}`: Vacuum energy matrix Wᵛ (`num_modes × num_modes`), block-diagonal in n for 2D - `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 + unless `compute_vacuum_response` is called with `compute_Iv=true`. 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`) diff --git a/src/Vacuum/Kernel3D.jl b/src/Vacuum/Kernel3D.jl index 3d135295b..b2e5f3a51 100644 --- a/src/Vacuum/Kernel3D.jl +++ b/src/Vacuum/Kernel3D.jl @@ -121,9 +121,11 @@ function SingularQuadratureData(PATCH_RAD::Int, RAD_DIM::Int, INTERP_ORDER::Int) x0 = 0.5 + 0.5 * qx[ir] * cos(dθ * (ia - 1)) x1 = 0.5 + 0.5 * qx[ir] * sin(dθ * (ia - 1)) - # Lower-left corner indices of INTERP_ORDER × INTERP_ORDER stencil centered on (x0,x1) - y0 = clamp(trunc(Int, x0 * (PATCH_DIM - 1) - (INTERP_ORDER - 1) ÷ 2), 0, PATCH_DIM - INTERP_ORDER) - y1 = clamp(trunc(Int, x1 * (PATCH_DIM - 1) - (INTERP_ORDER - 1) ÷ 2), 0, PATCH_DIM - INTERP_ORDER) + # Lower-left corner indices of INTERP_ORDER × INTERP_ORDER stencil centered on (x0,x1). + # Round to the nearest node before offsetting so the selection is equivariant under the patch's + # π-rotation, which is what lets a stellarator-symmetric surface produce a symmetric operator. + y0 = clamp(round(Int, x0 * (PATCH_DIM - 1)) - (INTERP_ORDER - 1) ÷ 2, 0, PATCH_DIM - INTERP_ORDER) + y1 = clamp(round(Int, x1 * (PATCH_DIM - 1)) - (INTERP_ORDER - 1) ÷ 2, 0, PATCH_DIM - INTERP_ORDER) # Local coordinates within INTERP_ORDER×INTERP_ORDER stencil, normalized to [0,1] z0 = (x0 * (PATCH_DIM - 1) - y0) * h @@ -321,7 +323,7 @@ function KernelWorkspace(PATCH_DIM::Int, RAD_DIM::Int, ANG_DIM::Int) end """ - compute_3D_kernel_matrices!(grad_greenfunction, greenfunction, observer, source, PATCH_RAD, RAD_DIM, INTERP_ORDER) + compute_3D_kernel_matrices!(grad_blocks, green_blocks, observer, source, PATCH_RAD, RAD_DIM, INTERP_ORDER, phases, sym=nothing) Compute boundary integral kernel matrices for 3D geometries with the singular correction algorithm from [Malhotra Plasma Phys. and Cont. Fusion 2019 024004]. @@ -334,43 +336,53 @@ grad_greenfunction is the double-layer kernel matrix, where each entry is ∇_{x_src} φ(x_obs, x_src) · n_src, and greenfunction is the single-layer kernel matrix, where each entry is φ(x_obs, x_src). -Takes advantage of field periodicity to evaluate the kernel only over a single field period. +Takes advantage of field periodicity to evaluate the kernel only over a single field period, in both +the observer and the source index: a source in field period `d` is accumulated onto the period-0 column +with weight `phases[d+1]`, so the block-circulant reduction is applied as the kernel is written and the +full-torus source blocks are never stored. # Arguments - - `grad_greenfunction`: Double-layer kernel matrix (Nobs × Nsrc) filled in place - - `greenfunction`: Single-layer kernel matrix (Nobs × Nsrc); filled only when `source` is plasma + - `grad_blocks`: Double-layer operator as a list of blocks, filled in place. Without `sym` this is + a one-element list holding the whole `(Nobs × Nsrc)` matrix. + - `green_blocks`: Single-layer operator, same layout; filled only when `source` is plasma - `observer`: Observer geometry (PlasmaGeometry3D) - `source`: Source geometry (PlasmaGeometry3D) - `PATCH_RAD`: Number of points adjacent to source point to treat as singular - `RAD_DIM`: Polar radial quadrature order. Angular order = 2 * RAD_DIM - `INTERP_ORDER`: Lagrange interpolation order, must be ≤ (2 * PATCH_RAD + 1) + - `phases`: Field-period phases `ω^{k d}` for `d = 0 … nfp-1`, with `nfp = length(phases)`. Pass the + real `[1.0]` for a single period, which keeps the output matrices real. + - `sym`: [`StellaratorBasis`](@ref) for this residue class, or `nothing` to build the untransformed + operator. When given, only the involution orbit representatives are evaluated — half the observer + points — and each row is emitted in the basis that makes the operator real and, for a + self-conjugate class, block diagonal. """ function compute_3D_kernel_matrices!( - grad_greenfunction::AbstractMatrix{Float64}, - greenfunction::AbstractMatrix{Float64}, + grad_blocks::AbstractVector{<:AbstractMatrix{<:Number}}, + green_blocks::AbstractVector{<:AbstractMatrix{<:Number}}, observer::Union{PlasmaGeometry3D,WallGeometry3D}, source::Union{PlasmaGeometry3D,WallGeometry3D}, PATCH_RAD::Int, RAD_DIM::Int, - INTERP_ORDER::Int + INTERP_ORDER::Int, + phases::AbstractVector{<:Number}, + sym::Union{Nothing,StellaratorBasis}=nothing ) num_points = observer.mtheta * observer.nzeta - n_obs = size(greenfunction, 1) # num_points ÷ nfp + num_points_per_fp = num_points ÷ length(phases) # observer/source points in one field period dθdζ = 4π^2 / num_points - # Get block of grad green function matrix + # Surface sub-block of the operator this call fills col_index = (source isa PlasmaGeometry3D ? 1 : 2) row_index = (observer isa PlasmaGeometry3D ? 1 : 2) - grad_greenfunction_block = view( - grad_greenfunction, - ((row_index-1)*n_obs+1):(row_index*n_obs), - ((col_index-1)*num_points+1):(col_index*num_points) - ) # 𝒢ⁿ only needed for plasma as source term (RHS of eqs. 26/27 in Chance 1997) populate_greenfunction = source isa PlasmaGeometry3D - populate_greenfunction && fill!(greenfunction, 0.0) + + # With the symmetry the partner of each involution orbit follows from its representative, so only + # the representatives are evaluated + observers = sym === nothing ? (1:num_points_per_fp) : sym.orbit_rep # This allows the code to run at lower resolution without erroring out, but will warn the user. if PATCH_RAD > (min(source.mtheta, source.nzeta) - 1) ÷ 2 @@ -382,18 +394,29 @@ function compute_3D_kernel_matrices!( quad_data = get_singular_quadrature(PATCH_RAD, RAD_DIM, INTERP_ORDER) (; PATCH_DIM, PATCH_RAD, ANG_DIM, RAD_DIM, Ppou, Gpou, P2G) = quad_data - # Allocate thread-local workspaces (one per thread) + # Allocate thread-local workspaces (one per thread). One operator row is accumulated at a time so + # the field-period fold and the basis change both happen before anything is stored. max_threadid = Threads.maxthreadid() workspaces = [KernelWorkspace(PATCH_DIM, RAD_DIM, ANG_DIM) for _ in 1:max_threadid] + Trow = eltype(phases) + rows_double = [zeros(Trow, num_points_per_fp) for _ in 1:max_threadid] + rows_single = [zeros(Trow, num_points_per_fp) for _ in 1:max_threadid] + partner_rows = [zeros(ComplexF64, sym === nothing ? 0 : num_points_per_fp) for _ in 1:max_threadid] # Parallel loop through observer points - Threads.@threads for idx_obs in 1:n_obs + Threads.@threads for i_orbit in eachindex(observers) # Get thread-local workspace - ws = workspaces[Threads.threadid()] + tid = Threads.threadid() + ws = workspaces[tid] (; r_patch, dr_dθ_patch, dr_dζ_patch, r_polar, dr_dθ_polar, dr_dζ_polar, n_polar, M_polar_single, M_polar_double, M_grid_single_flat, M_grid_double_flat) = ws + row_double = rows_double[tid] + row_single = rows_single[tid] + fill!(row_double, 0) + populate_greenfunction && fill!(row_single, 0) # Convert linear index to 2D indices + idx_obs = observers[i_orbit] i_obs = mod1(idx_obs, observer.mtheta) j_obs = (idx_obs - 1) ÷ observer.mtheta + 1 r_obs = @view observer.r[idx_obs, :] @@ -408,11 +431,12 @@ function compute_3D_kernel_matrices!( n_src = @view source.normal[idx_src, :] far_single, far_double = laplace_kernel(r_obs[1], r_obs[2], r_obs[3], r_src[1], r_src[2], r_src[3], n_src[1], n_src[2], n_src[3]) - # Apply weights (periodic trapezoidal rule = constant weights) + # Periodic trapezoidal rule (constant weights); fold this source's field period onto period 0 + d, idx_col = fldmod1(idx_src, num_points_per_fp) if populate_greenfunction - greenfunction[idx_obs, idx_src] = far_single * dθdζ + row_single[idx_col] += phases[d] * (far_single * dθdζ) end - grad_greenfunction_block[idx_obs, idx_src] = far_double * dθdζ + row_double[idx_col] += phases[d] * (far_double * dθdζ) end # ============================================================ @@ -458,6 +482,7 @@ function compute_3D_kernel_matrices!( idx_pol = periodic_wrap(i_obs - PATCH_RAD + i - 1, source.mtheta) idx_tor = periodic_wrap(j_obs - PATCH_RAD + j - 1, source.nzeta) idx_src = idx_pol + source.mtheta * (idx_tor - 1) + d, idx_col = fldmod1(idx_src, num_points_per_fp) # Remainder of far-field contribution on the singular grid: Gpou = -χ r_src = @view source.r[idx_src, :] @@ -466,21 +491,34 @@ function compute_3D_kernel_matrices!( # Apply near + far contributions if populate_greenfunction - greenfunction[idx_obs, idx_src] += M_grid_single[i, j] + far_single * Gpou[i, j] * dθdζ + row_single[idx_col] += phases[d] * (M_grid_single[i, j] + far_single * Gpou[i, j] * dθdζ) end - grad_greenfunction_block[idx_obs, idx_src] += M_grid_double[i, j] + far_double * Gpou[i, j] * dθdζ + row_double[idx_col] += phases[d] * (M_grid_double[i, j] + far_double * Gpou[i, j] * dθdζ) end - end - # Use the same normalization as in the 2D kernel so we can just add I to the diagonal - # This makes the grri logic identical to the 2D kernel. - grad_greenfunction_block ./= 2π - populate_greenfunction && (greenfunction ./= 2π) + # Normalize so the Green's-identity term below is a unit shift. The exterior/interior jump is + # then 2I, the same scalar shift the 2D kernel carries, so the grri logic is identical. + row_double ./= 2π + populate_greenfunction && (row_single ./= 2π) + + if sym === nothing + emit_plain_row!(grad_blocks, row_double, idx_obs, row_index, col_index) + populate_greenfunction && emit_plain_row!(green_blocks, row_single, idx_obs, row_index, 1) + else + partner = partner_rows[tid] + emit_symmetric_row!(grad_blocks, sym, row_double, partner, i_orbit, row_index, col_index) + populate_greenfunction && emit_symmetric_row!(green_blocks, sym, row_single, partner, i_orbit, row_index, 1) + end + end - # Add the term that comes from the volume integral of Green's identity + # Add the term that comes from the volume integral of Green's identity. The identity is invariant + # under the basis change, so it is the same unit shift on each transformed block. if typeof(source) == typeof(observer) - for i in 1:n_obs - grad_greenfunction_block[i, i] += 1.0 + for (b, H) in enumerate(grad_blocks) + nsize = sym === nothing ? num_points_per_fp : sym.block_size[b] + for i in 1:nsize + H[(row_index-1)*nsize+i, (col_index-1)*nsize+i] += 1.0 + end end end end diff --git a/src/Vacuum/Symmetry3D.jl b/src/Vacuum/Symmetry3D.jl new file mode 100644 index 000000000..d455d5def --- /dev/null +++ b/src/Vacuum/Symmetry3D.jl @@ -0,0 +1,260 @@ +""" + StellaratorBasis + +Change of basis for one toroidal residue class `k` that makes the reduced double-layer operator +`D̂ₖ` real, and block-diagonal when the class is self-conjugate. + +A stellarator-symmetric boundary is invariant under `(θ, ζ) → (−θ, −ζ)` with `Z → −Z`, i.e. under +the rotation `diag(1, −1, −1)`. Both Laplace kernels depend only on `|r_obs − r_src|` and +`n_src·(r_obs − r_src)`, so the operator inherits the symmetry. Under the field-period fold the +involution acts *antilinearly* on the reduced operator, + + D̂ₖ[σp, σq] = ω^{k(aₚ − a_q)} · conj(D̂ₖ[p, q]), ω = exp(-2πi/nfp), + +with `σ` the within-period involution and `a = 0` on the `ζ = 0` symmetry plane, `1` elsewhere. +Two cases follow, and the operator memory halves in both: + + - **Self-conjugate class** (`mod(2k, nfp) == 0`, which is every class when `nfp == 1`): `ω^k = ±1` + and `D̂ₖ` is real, commuting with the real signed involution `J = ΛΠ`, `Λ[p] = (ω^k)^{aₚ}`. The + operator splits into the two `J`-eigenspaces — two real blocks of roughly half size each. + - **Otherwise**: the involution is antiunitary with `Θ² = +1`, so the half-twist `λ = ω^{k a / 2}`, + the parity basis, and a factor `i` on the odd half make the operator real at full size. + +The partner row of an orbit is never evaluated — it follows from the relation above — so the kernel +visits only half the observer points. + +## Fields + + - `σ`: within-period involution as a grid-index permutation of `1:num_points_per_fp` + - `self_conjugate`: whether the class splits into two blocks + - `λ`: half-twist `ω^{k a_p / 2}` per grid point (non-self-conjugate case) + - `Λ`: real sign `(ω^k)^{a_p}` per grid point (self-conjugate case) + - `orbit_rep`, `orbit_partner`: the two grid points of each involution orbit, equal for a fixed point + - `orbit_cols`: the one or two basis columns each orbit produces, as indices into `block`/`slot` + - `col_p`, `col_q`, `col_cp`, `col_cq`: basis column `c` is `col_cp[c]·e_{col_p[c]} + col_cq[c]·e_{col_q[c]}`, + with `col_q == col_p` and `col_cq == 0` for a fixed point of `σ` + - `block`, `slot`: destination block of each basis column and its position within that block + - `block_size`: number of basis columns per surface in each block +""" +struct StellaratorBasis + σ::Vector{Int} + self_conjugate::Bool + λ::Vector{ComplexF64} + Λ::Vector{Int} + orbit_rep::Vector{Int} + orbit_partner::Vector{Int} + orbit_cols::Vector{Vector{Int}} + col_p::Vector{Int} + col_q::Vector{Int} + col_cp::Vector{ComplexF64} + col_cq::Vector{ComplexF64} + block::Vector{Int} + slot::Vector{Int} + block_size::Vector{Int} +end + +""" + stellarator_involution(plasma, wall, nfp) -> Union{Vector{Int},Nothing} + +Within-period grid-index involution `σ` if both surfaces are stellarator symmetric, `nothing` +otherwise. The test compares `r(−θ, −ζ)` against `diag(1,−1,−1)·r(θ, ζ)` on the full torus; a +symmetric boundary matches to round-off while an asymmetric one is off by a fraction of the minor +radius, so the threshold is not delicate. +""" +function stellarator_involution(plasma::PlasmaGeometry3D, wall::WallGeometry3D, nfp::Int) + mtheta, nzeta_full = plasma.mtheta, plasma.nzeta + nzeta_full % nfp == 0 || return nothing + num_points = mtheta * nzeta_full + σ_full = [mod1(2 - mod1(p, mtheta), mtheta) + mtheta * (mod1(2 - ((p - 1) ÷ mtheta + 1), nzeta_full) - 1) for p in 1:num_points] + + function symmetric(r) + tol = 1e-10 * maximum(abs, r) + return all(p -> abs(r[σ_full[p], 1] - r[p, 1]) ≤ tol && abs(r[σ_full[p], 2] + r[p, 2]) ≤ tol && abs(r[σ_full[p], 3] + r[p, 3]) ≤ tol, 1:num_points) + end + symmetric(plasma.r) || return nothing + (wall.nowall || symmetric(wall.r)) || return nothing + + # Within one field period the involution is composed with a field-period rotation, so the ζ = 0 + # line maps to itself and every other line reflects about it. + nzeta = nzeta_full ÷ nfp + return [mod1(2 - mod1(p, mtheta), mtheta) + mtheta * (mod1(nzeta + 2 - ((p - 1) ÷ mtheta + 1), nzeta) - 1) for p in 1:(mtheta*nzeta)] +end + +""" + StellaratorBasis(σ, mtheta, k, nfp) + +Build the basis for toroidal residue class `k` from the within-period involution `σ`. +""" +function StellaratorBasis(σ::Vector{Int}, mtheta::Int, k::Int, nfp::Int) + npts = length(σ) + a(p) = ((p - 1) ÷ mtheta + 1) == 1 ? 0 : 1 # ζ = 0 plane is the symmetry plane and carries no twist + self_conjugate = mod(2k, nfp) == 0 + ωk = cis(-2π * k / nfp) + λ = ComplexF64[cis(-π * k * a(p) / nfp) for p in 1:npts] + Λ = Int[round(Int, real(ωk))^a(p) for p in 1:npts] + + orbit_rep = Int[] + orbit_partner = Int[] + seen = falses(npts) + for p in 1:npts + seen[p] && continue + seen[p] = true + seen[σ[p]] = true + push!(orbit_rep, p) + push!(orbit_partner, σ[p]) + end + + # Each orbit contributes one column (fixed point) or a symmetric/antisymmetric pair. In the + # self-conjugate case the pair splits by J-eigenvalue ±Λ[p]; otherwise everything is one block. + block = Int[] + col_p = Int[] + col_q = Int[] + col_cp = ComplexF64[] + col_cq = ComplexF64[] + orbit_cols = [Int[] for _ in orbit_rep] + rt = 1 / √2 + function add_column!(blk, p, q, cp, cq) + push!(block, blk) + push!(col_p, p) + push!(col_q, q) + push!(col_cp, cp) + push!(col_cq, cq) + return length(block) + end + + for (o, (p, q)) in enumerate(zip(orbit_rep, orbit_partner)) + # Symmetric combination first, antisymmetric second; the half-twist is absorbed into the + # coefficients when the class is not self-conjugate. + c1 = self_conjugate ? ComplexF64(1) : λ[p] + blk1 = self_conjugate && Λ[p] == -1 ? 2 : 1 + if q == p + push!(orbit_cols[o], add_column!(blk1, p, p, c1, 0)) + else + push!(orbit_cols[o], add_column!(blk1, p, q, c1 * rt, c1 * rt)) + c2 = self_conjugate ? ComplexF64(rt) : im * λ[p] * rt + blk2 = self_conjugate ? (Λ[p] == -1 ? 1 : 2) : 1 + push!(orbit_cols[o], add_column!(blk2, p, q, c2, -c2)) + end + end + + nb = self_conjugate ? 2 : 1 + block_size = [count(==(b), block) for b in 1:nb] + slot = zeros(Int, length(block)) + filled = zeros(Int, nb) + for c in eachindex(block) + filled[block[c]] += 1 + slot[c] = filled[block[c]] + end + return StellaratorBasis(σ, self_conjugate, λ, Λ, orbit_rep, orbit_partner, orbit_cols, col_p, col_q, col_cp, col_cq, block, slot, block_size) +end + +""" + emit_symmetric_row!(dest, basis, row, work, orbit, row_index, col_index) + +Scatter the raw operator row of one involution-orbit representative into the transformed blocks. + +`row` holds the untransformed operator row over one source surface and `work` is scratch of the same +length. `dest[b]` receives block `b`, with `row_index`/`col_index` (1 = plasma, 2 = wall) selecting +the surface sub-block. The partner row is reconstructed from `row`, so the caller evaluates only the +orbit representatives. +""" +function emit_symmetric_row!( + dest::AbstractVector{<:AbstractMatrix{Float64}}, + basis::StellaratorBasis, + row::AbstractVector{<:Number}, + work::AbstractVector{ComplexF64}, + orbit::Int, + row_index::Int, + col_index::Int +) + (; σ, self_conjugate, λ, Λ, orbit_rep, orbit_partner, orbit_cols, block, slot) = basis + p, q = orbit_rep[orbit], orbit_partner[orbit] + fixed_row = p == q + rt = √2 + + # Reconstruct the partner row from the representative: real and signed when the class is + # self-conjugate, conjugated in the half-twisted frame otherwise. + if self_conjugate + @inbounds for x in eachindex(row) + work[x] = Λ[p] * Λ[x] * real(row[σ[x]]) + end + else + cλp = conj(λ[p]) + @inbounds for x in eachindex(row) + work[x] = cλp * row[x] * λ[x] + end + end + + for (n, c) in enumerate(orbit_cols[orbit]) + blk = block[c] + H = dest[blk] + nsize = basis.block_size[blk] + r = (row_index - 1) * nsize + slot[c] + col_offset = (col_index - 1) * nsize + even_row = n == 1 + sgn = even_row ? 1.0 : -1.0 + @inbounds for oc in eachindex(orbit_rep) + x, y = orbit_rep[oc], orbit_partner[oc] + fixed_col = x == y + if self_conjugate + vx = fixed_row ? real(row[x]) : (real(row[x]) + sgn * real(work[x])) / rt + vy = fixed_col ? 0.0 : (fixed_row ? real(row[y]) : (real(row[y]) + sgn * real(work[y])) / rt) + for (m, c2) in enumerate(orbit_cols[oc]) + # Only the diagonal block survives; the rest vanishes by the symmetry + block[c2] == blk || continue + H[r, col_offset+slot[c2]] = fixed_col ? vx : (m == 1 ? (vx + vy) / rt : (vx - vy) / rt) + end + else + cx = work[x] + cy = fixed_col ? cx : work[y] + for (m, c2) in enumerate(orbit_cols[oc]) + val = if fixed_col + fixed_row ? real(cx) : rt * (even_row ? real(cx) : imag(cx)) + elseif fixed_row + m == 1 ? rt * real(cx) : -rt * imag(cx) + elseif even_row + m == 1 ? real(cx) + real(cy) : imag(cy) - imag(cx) + else + m == 1 ? imag(cx) + imag(cy) : real(cx) - real(cy) + end + H[r, col_offset+slot[c2]] = val + end + end + end + end + return nothing +end + +""" + emit_plain_row!(dest, row, idx_obs, row_index, col_index) + +Write one accumulated operator row into the untransformed operator, the `sym === nothing` fall +through. `dest` is a one-element block list holding the full matrix. +""" +function emit_plain_row!(dest::AbstractVector{<:AbstractMatrix}, row::AbstractVector{<:Number}, idx_obs::Int, row_index::Int, col_index::Int) + H = dest[1] + npts = length(row) + r = (row_index - 1) * npts + idx_obs + col_offset = (col_index - 1) * npts + @inbounds for x in eachindex(row) + H[r, col_offset+x] = row[x] + end + return nothing +end + +""" + transform_mode_basis!(dest, basis_matrix, sym) + +Express the Fourier mode basis in the symmetry-adapted basis, `Ẽ = E·U`, writing block `b` into +`dest[b]`. Substituting `Ẽ` for `E` carries the change of basis through the right-hand side, the +`wv` projection and `I_v` unchanged, since `U` is unitary. +""" +function transform_mode_basis!(dest::AbstractVector{<:AbstractMatrix{ComplexF64}}, basis_matrix::AbstractMatrix, sym::StellaratorBasis) + for c in eachindex(sym.block) + out = @view dest[sym.block[c]][:, sym.slot[c]] + p, q = sym.col_p[c], sym.col_q[c] + @views out .= sym.col_cp[c] .* basis_matrix[:, p] + q != p && (@views out .+= sym.col_cq[c] .* basis_matrix[:, q]) + end + return nothing +end diff --git a/src/Vacuum/Vacuum.jl b/src/Vacuum/Vacuum.jl index 9c4921618..bf9d348c7 100644 --- a/src/Vacuum/Vacuum.jl +++ b/src/Vacuum/Vacuum.jl @@ -13,6 +13,7 @@ using ..Utilities.FourierTransforms: FourierTransform, compute_fourier_coefficie include("Utilities.jl") include("DataTypes.jl") +include("Symmetry3D.jl") include("PnQuadCache.jl") include("Kernel2D.jl") include("Kernel3D.jl") @@ -49,6 +50,24 @@ function _warn_and_symmetrize!(mat::AbstractMatrix, name::String) hermitianpart!(mat) end +""" + _fill_Iv_block!(I_v_block, basis, grre, grri, num_points) + +Surface-current matrix Iᵛ from the exterior and interior potentials, Park 2007 eq. 21b: +`μ₀I^v = χ^(vi) - χ^(vo)`. Overwrites the plasma-observer rows of `grri`. + +The difference is taken as `grre - grri` and the result conjugated because VACUUM builds the +operators in its CW-θ frame while GPEC uses CCW-θ, flipping the outward-normal sign. +""" +function _fill_Iv_block!(I_v_block::AbstractMatrix, basis::AbstractMatrix, grre::AbstractMatrix, grri::AbstractMatrix, num_points::Int) + g_diff = @view grri[1:num_points, :] + g_diff .= @view(grre[1:num_points, :]) .- g_diff + mul!(I_v_block, basis, g_diff) + conj!(I_v_block) # Flip θ_VAC → -θ_VAC to get I^v in GPEC's CCW-θ frame. + I_v_block ./= num_points + return I_v_block +end + """ _compute_vacuum_response_2d!(vac_data::VacuumResponse, inputs::VacuumInput, wall_settings::WallShapeSettings; compute_Iv=false) @@ -113,7 +132,7 @@ Green's functions are internal scratch only. 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. + # grre = -(2π)²χ^(vo), the vacuum-outside potential. Overwrites the block to save memory. ldiv!(lu!(grad_green), grre) # Interior operator D_int = D_ext - 2I: the double-layer jump between the two one-sided @@ -123,15 +142,7 @@ Green's functions are internal scratch only. 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 + _fill_Iv_block!(@view(vac_data.I_v[block_idx, block_idx]), ft.basis, grre, grri, num_points_surf) else # Only need exterior system for wv ldiv!(lu!(grad_green), grre) @@ -158,105 +169,180 @@ Green's functions are internal scratch only. end """ - _compute_vacuum_response_3d!(vac_data::VacuumResponse, inputs::VacuumInput, wall_settings::WallShapeSettings; compute_Iv=false) + _compute_vacuum_response_3d!(vac_data::VacuumResponse, inputs::VacuumInput, wall_settings::WallShapeSettings; compute_Iv=false, use_symmetry=true) -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 -direct real solve on the built kernel matrices (equivalent to 2D). +3D (`inputs.nzeta > 1`) vacuum response via block-circulant field-period reduction. -The `nfp`-periodic boundary makes the single-/double-layer operators `S`, `D` block-circulant in -the field-period index. Writing the response as `wv = (4π²/N)·Eᴴ·D⁻¹S·E` (`E` the complex Fourier -basis, `N = mtheta·nzeta_full`), the structure block-diagonalizes the problem by toroidal residue -class `k = mod(n, nfp)`: modes with different `k` do not couple, and within a class +The `nfp`-periodic boundary makes the single-/double-layer operators `S`, `D` block-circulant in the +field-period index, so the problem block-diagonalizes by toroidal residue class `k = mod(n, nfp)`: +modes with different `k` do not couple, and within a class D̂ₖ = Σ_d D_d ω^{k d}, Ŝₖ = Σ_d S_d ω^{k d}, ω = exp(-2πi/nfp), -with `D_d`, `S_d` the blocks coupling observers in field period 0 to sources in period `d`. Each -class needs one solve `wv[class k] = (4π²/M)·E_localᴴ·(D̂ₖ \\ Ŝₖ)|_plasma·E_local`. Only the first -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_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 `D_d`, `S_d` the blocks coupling observers in field period 0 to sources in period `d`. Each class +needs one solve `wv[class k] = (4π²/M)·E_localᴴ·(D̂ₖ \\ Ŝₖ)|_plasma·E_local` (`E` the complex Fourier +basis, `M = mtheta·nzeta`), so the routine loops over classes exactly as the 2D routine loops over +decoupled `n`. The phase sum is folded into the kernel write, so only the reduced `[nb·M × nb·M]` +operator is ever stored; for `nfp == 1` every phase is unity and the operators stay real. + +When both surfaces are stellarator symmetric the class operator is additionally transformed by the +[`StellaratorBasis`](@ref) for that class, which makes it real and — for a self-conjugate class — +splits it into two blocks of roughly half the size. That halves the operator memory and the kernel +work, and cuts the factorization ~2.6-2.8×. Pass `use_symmetry=false` to force the untransformed +solve; a boundary that fails the symmetry test falls through to it automatically. + +With `compute_Iv=true` each class also solves the interior operator `D_int = D_ext - 2I`, as in 2D. +That shift is block-diagonal in the field-period index and invariant under the basis change, so both +reductions stay exact. """ -@with_pool pool function _compute_vacuum_response_3d!(vac_data::VacuumResponse, inputs::VacuumInput, wall_settings::WallShapeSettings; compute_Iv::Bool=false) +@with_pool pool function _compute_vacuum_response_3d!( + vac_data::VacuumResponse, + inputs::VacuumInput, + wall_settings::WallShapeSettings; + compute_Iv::Bool=false, + use_symmetry::Bool=true +) (; 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) plasma_surf = PlasmaGeometry3D(full) wall = WallGeometry3D(full, plasma_surf, wall_settings) - num_points_per_fp = mtheta * nzeta # points per field period - num_points = num_points_per_fp * nfp # full-torus point count + num_points_per_fp = mtheta * nzeta # points per field period mpert = length(m_modes) nb = wall.nowall ? 1 : 2 # surface blocks: plasma, or [plasma; wall] - n_obs = nb * num_points_per_fp # number of observer points per field period + n_obs = nb * num_points_per_fp # observer rows: plasma (and wall) points of one field period # Complex Fourier basis exp(-i(mθ-nζ)) on a single field period exp_mn_basis = compute_fourier_coefficients(mtheta, m_modes, nzeta * nfp, n_modes; nfp=nfp) - # Local work matrices - grad_green = zeros!(pool, n_obs, nb * num_points) # double-layer D (with +I on the period-0 diagonals) - green_temp = zeros!(pool, n_obs, num_points) # single-layer S (plasma sources only) - # Kernel parameters, hardcoded for now PATCH_RAD = 11 RAD_DIM = 20 INTERP_ORDER = 5 - # Plasma–plasma block - compute_3D_kernel_matrices!(grad_green, view(green_temp, 1:num_points_per_fp, :), plasma_surf, plasma_surf, PATCH_RAD, RAD_DIM, INTERP_ORDER) - - if !wall.nowall - # Plasma–Wall block - compute_3D_kernel_matrices!(grad_green, view(green_temp, 1:num_points_per_fp, :), plasma_surf, wall, PATCH_RAD, RAD_DIM, INTERP_ORDER) - # Wall–Plasma block - compute_3D_kernel_matrices!(grad_green, view(green_temp, (num_points_per_fp+1):(2*num_points_per_fp), :), wall, plasma_surf, PATCH_RAD, RAD_DIM, INTERP_ORDER) - # Wall–Wall block - compute_3D_kernel_matrices!(grad_green, view(green_temp, (num_points_per_fp+1):(2*num_points_per_fp), :), wall, wall, PATCH_RAD, RAD_DIM, INTERP_ORDER) + # Stellarator symmetry halves both the operator and the kernel work when the surfaces allow it + σ = use_symmetry ? stellarator_involution(plasma_surf, wall, nfp) : nothing + classes = unique(mod.(n_modes, nfp)) + bases = [σ === nothing ? nothing : StellaratorBasis(σ, mtheta, k, nfp) for k in classes] + block_sizes = [b === nothing ? [num_points_per_fp] : b.block_size for b in bases] + + # One flat buffer per operator, carved into this class's blocks each pass. Sized for the widest + # class so the pool does not grow with the number of classes. + T = σ === nothing ? (nfp == 1 ? Float64 : ComplexF64) : Float64 + grad_len = maximum(sum((nb * sz)^2 for sz in szs) for szs in block_sizes) + green_len = maximum(sum(nb * sz * sz for sz in szs) for szs in block_sizes) + grad_buffer = zeros!(pool, T, grad_len) + green_buffer = zeros!(pool, T, green_len) + interior_buffer = compute_Iv ? zeros!(pool, T, grad_len) : zeros!(pool, T, 0) + grre = zeros!(pool, ComplexF64, n_obs, mpert * length(n_modes)) + grri = compute_Iv ? similar!(pool, grre) : zeros!(pool, ComplexF64, 0, 0) + basis_buffer = zeros!(pool, ComplexF64, mpert * length(n_modes), σ === nothing ? 0 : num_points_per_fp) + + # Carve a flat buffer into this class's blocks; a reshaped contiguous view stays strided, so the + # factorizations and matrix products below stay on the BLAS path. + function carve(buffer, szs, ncols) + offsets = cumsum([0; [nb * sz * ncols(sz) for sz in szs]]) + return [reshape(view(buffer, (offsets[i]+1):offsets[i+1]), nb * szs[i], ncols(szs[i])) for i in eachindex(szs)] end - if nfp == 1 - # nfp == 1 case: direct real solve on the built kernel matrices (equivalent to 2D) - ldiv!(lu!(grad_green), green_temp) - G_plasma = @view green_temp[1:num_points_per_fp, :] - vac_data.wv .= (4π^2 / num_points_per_fp) .* (exp_mn_basis * (G_plasma * exp_mn_basis')) - else - # nfp > 1 case: block-circulant solve on the built kernel matrices - # Solve one reduced (nb·M)×(nb·M) system per toroidal residue class k = mod(n, nfp) - D̂ = zeros!(pool, ComplexF64, n_obs, n_obs) - Ŝ = zeros!(pool, ComplexF64, n_obs, num_points_per_fp) - for k in unique(mod.(n_modes, nfp)) - # Reduced operators D̂ₖ = Σ_d D_d ω^{k d}, Ŝₖ = Σ_d S_d ω^{k d} (fused, allocation-free). - # D̂ keeps the nb-block source structure (plasma cols 1:M, wall cols M+1:2M); Ŝ is plasma- - # source-only over all nb·M observer rows. - fill!(D̂, 0) - fill!(Ŝ, 0) - for d in 0:(nfp-1) - phase = cis(-2π * (k * d) / nfp) - for s in 0:(nb-1) - src_cols = (s*num_points+d*num_points_per_fp+1):(s*num_points+(d+1)*num_points_per_fp) - dst_cols = (s*num_points_per_fp+1):((s+1)*num_points_per_fp) - @views @. D̂[:, dst_cols] += phase * grad_green[:, src_cols] + # Loop over all decoupled toroidal residue classes + for (idx_k, k) in enumerate(classes) + cols = [(idx_m + (idx_n-1)*mpert) for (idx_n, n) in enumerate(n_modes) if mod(n, nfp) == k for idx_m in 1:mpert] + # A contiguous class (always so for nfp == 1) keeps the basis and output views strided, and so on the BLAS path + mode_cols = length(cols) == cols[end] - cols[1] + 1 ? (cols[1]:cols[end]) : cols + # Diagonal block of wv (and I_v when requested) + wv_block = @view vac_data.wv[mode_cols, mode_cols] + E = @view exp_mn_basis[mode_cols, :] + sym = bases[idx_k] + szs = block_sizes[idx_k] + + # Phases are real whenever ω^k = ±1, which keeps the whole class on the real BLAS path + phases = if (σ === nothing ? nfp == 1 : mod(2k, nfp) == 0) + sgn = mod(2k, nfp) == 0 && isodd(2k ÷ nfp) ? -1.0 : 1.0 + Float64[sgn^d for d in 0:(nfp-1)] + else + ComplexF64[cis(-2π * (k * d) / nfp) for d in 0:(nfp-1)] + end + + grad_blocks = carve(grad_buffer, szs, sz -> nb * sz) + green_blocks = carve(green_buffer, szs, sz -> sz) + + # Plasma–Plasma block + compute_3D_kernel_matrices!(grad_blocks, green_blocks, plasma_surf, plasma_surf, PATCH_RAD, RAD_DIM, INTERP_ORDER, phases, sym) + + if !wall.nowall + # Plasma–Wall block + compute_3D_kernel_matrices!(grad_blocks, green_blocks, plasma_surf, wall, PATCH_RAD, RAD_DIM, INTERP_ORDER, phases, sym) + # Wall–Plasma block + compute_3D_kernel_matrices!(grad_blocks, green_blocks, wall, plasma_surf, PATCH_RAD, RAD_DIM, INTERP_ORDER, phases, sym) + # Wall–Wall block + compute_3D_kernel_matrices!(grad_blocks, green_blocks, wall, wall, PATCH_RAD, RAD_DIM, INTERP_ORDER, phases, sym) + end + + # Mode basis in the symmetry-adapted basis; Ẽ = E·U carries the transform through the solve + mode_basis = if sym === nothing + [E] + else + mb = [view(basis_buffer, 1:length(mode_cols), (sum(szs[1:(i-1)])+1):sum(szs[1:i])) for i in eachindex(szs)] + transform_mode_basis!(mb, E, sym) + mb + end + + interior_blocks = compute_Iv ? carve(interior_buffer, szs, sz -> nb * sz) : grad_blocks + + row_offset = 0 + for (b, sz) in enumerate(szs) + nrow = nb * sz + rows = (row_offset+1):(row_offset+nrow) + Ẽ = mode_basis[b] + grre_k = @view grre[rows, 1:length(mode_cols)] + + # The mode basis acts on columns and D⁻¹ on rows, so (D⁻¹S)Eᴴ == D⁻¹(SEᴴ): projecting the + # RHS before the solve is exact and carries this class's modes instead of one column per point. + mul!(grre_k, green_blocks[b], Ẽ') + + if compute_Iv + # Copy RHS before exterior solve overwrites grre; keep a kernel copy for interior + grri_k = @view grri[rows, 1:length(mode_cols)] + grri_k .= grre_k + interior_blocks[b] .= grad_blocks[b] + + # Exterior operator D_ext = 2I + 𝒦 (Chance 1997 eq. 89); the solve gives + # grre = -(2π)²χ^(vo), the vacuum-outside potential. Overwrites the block to save memory. + ldiv!(lu!(grad_blocks[b]), grre_k) + + # 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:nrow + interior_blocks[b][i, i] -= 2.0 end - scols = (d*num_points_per_fp+1):((d+1)*num_points_per_fp) - @views @. Ŝ += phase * green_temp[:, scols] + ldiv!(lu!(interior_blocks[b]), grri_k) + + # μ₀Iᵛ = χ^(vi) - χ^(vo) (Park 2007 eq. 21b), accumulated over the parity blocks + g_diff = @view grri_k[1:sz, :] + g_diff .= @view(grre_k[1:sz, :]) .- g_diff + mul!(@view(vac_data.I_v[mode_cols, mode_cols]), Ẽ, g_diff, 1, 1) + else + # Only need exterior system for wv + ldiv!(lu!(grad_blocks[b]), grre_k) end - # Exterior response operator for this sector: solve D̂ₖ G = Ŝₖ in place (Ŝ ← G = D̂ₖ⁻¹Ŝₖ) - ldiv!(lu!(D̂), Ŝ) - mode_cols = [(idx_m + (idx_n-1)*mpert) for (idx_n, n) in enumerate(n_modes) if mod(n, nfp) == k for idx_m in 1:mpert] - E = @view exp_mn_basis[mode_cols, :] - G_plasma = @view Ŝ[1:num_points_per_fp, :] # plasma-observer rows of the combined exterior operator - vac_data.wv[mode_cols, mode_cols] .= (4π^2 / num_points_per_fp) .* (E * (G_plasma * E')) + # Project exterior kernel onto observer basis exp(-i(mθ-nζ)), summed over the blocks + mul!(wv_block, Ẽ, @view(grre_k[1:sz, :]), 1, 1) + row_offset += nrow + end + wv_block .*= 4π^2 / num_points_per_fp + + if compute_Iv + # Flip θ_VAC → -θ_VAC to get I^v in GPEC's CCW-θ frame, and normalize + Iv_block = @view vac_data.I_v[mode_cols, mode_cols] + conj!(Iv_block) + Iv_block ./= num_points_per_fp end end @@ -270,31 +356,31 @@ interior variant `-D + 2I` for the interior columns, then scatter back into the end """ - compute_vacuum_response(inputs::VacuumInput, wall_settings::WallShapeSettings; compute_Iv=false) -> VacuumResponse + compute_vacuum_response(inputs::VacuumInput, wall_settings::WallShapeSettings; compute_Iv=false, use_symmetry=true) -> 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. Pass `compute_Iv=true` to additionally populate the -surface-current matrix `I_v` (2D only). +surface-current matrix `I_v`, and `use_symmetry=false` to skip the 3D stellarator-symmetry solve. """ -function compute_vacuum_response(inputs::VacuumInput, wall_settings::WallShapeSettings; compute_Iv::Bool=false) +function compute_vacuum_response(inputs::VacuumInput, wall_settings::WallShapeSettings; compute_Iv::Bool=false, use_symmetry::Bool=true) vac = VacuumResponse(inputs) - compute_vacuum_response!(vac, inputs, wall_settings; compute_Iv) + compute_vacuum_response!(vac, inputs, wall_settings; compute_Iv, use_symmetry) return vac end """ - compute_vacuum_response!(vac_data::VacuumResponse, inputs::VacuumInput, wall_settings::WallShapeSettings; compute_Iv=false) + compute_vacuum_response!(vac_data::VacuumResponse, inputs::VacuumInput, wall_settings::WallShapeSettings; compute_Iv=false, use_symmetry=true) 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!`]. +[`_compute_vacuum_response_3d!`]. `use_symmetry` applies to the 3D path only. """ -function compute_vacuum_response!(vac_data::VacuumResponse, inputs::VacuumInput, wall_settings::WallShapeSettings; compute_Iv::Bool=false) +function compute_vacuum_response!(vac_data::VacuumResponse, inputs::VacuumInput, wall_settings::WallShapeSettings; compute_Iv::Bool=false, use_symmetry::Bool=true) if inputs.nzeta == 1 _compute_vacuum_response_2d!(vac_data, inputs, wall_settings; compute_Iv) else - _compute_vacuum_response_3d!(vac_data, inputs, wall_settings; compute_Iv) + _compute_vacuum_response_3d!(vac_data, inputs, wall_settings; compute_Iv, use_symmetry) end end diff --git a/test/runtests_slayer_params.jl b/test/runtests_slayer_params.jl index 330ba7297..bab15e393 100644 --- a/test/runtests_slayer_params.jl +++ b/test/runtests_slayer_params.jl @@ -165,4 +165,29 @@ @test_throws ArgumentError r_based_shear(0.5, 2.0, 1.0, 0.0) @test_throws ArgumentError r_based_shear(0.5, 0.0, 1.0, 0.5) end + + @testset "Test 3: reverse-shear invariance" begin + # The layer timescales and widths depend on |dq/dr|, not its sign: a negative-shear + # surface must reduce to its positive-shear mirror, with only the recorded sval_r + # diagnostic keeping the sign. dc_type=:lar with nonzero dr_val exercises the Wd + # iteration and the critical-Δ square roots as well as tau_h. + base = _ref_kwargs(; dr_val=-0.1, dc_type=:lar) + pos = slayer_parameters(; base...) + neg = slayer_parameters(; merge(base, (; sval_r=-1.0))...) + + # The sign survives where it is a diagnostic, not a magnitude. + @test pos.sval_r == 1.0 + @test neg.sval_r == -1.0 + + # Every normalized layer quantity is bit-identical: abs(-1.0) === 1.0, + # so the whole downstream chain reproduces exactly. + for f in (:tau, :lu, :c_beta, :D_norm, :P_perp, :P_tor, :Q_e, :Q_i, + :iota_e, :tauk, :tau_r, :delta_n, :eta, :d_beta, :dc_tmp) + @test getfield(neg, f) == getfield(pos, f) + end + + # tau_h > 0 keeps the Lundquist number positive, so S^(1/3) is defined on reverse shear. + @test neg.lu > 0 + @test neg.tauk > 0 + end end diff --git a/test/runtests_vacuum.jl b/test/runtests_vacuum.jl index 2db363573..d93bc6da7 100644 --- a/test/runtests_vacuum.jl +++ b/test/runtests_vacuum.jl @@ -778,6 +778,50 @@ @test isapprox(wv, wv', rtol=1e-12) end + @testset "compute_vacuum_response 3D compute_Iv=true" begin + num_modes(inp) = length(inp.m_modes) * length(inp.n_modes) + for wall_settings in (WallShapeSettings(shape="nowall"), WallShapeSettings(shape="conformal", a=0.3)) + inputs = _make_3d_inputs(mtheta=32, nzeta=32, mtheta_eq=17) + (; I_v) = compute_vacuum_response(inputs, wall_settings; compute_Iv=true) + @test size(I_v) == (num_modes(inputs), num_modes(inputs)) + @test all(isfinite, I_v) + @test !all(iszero, I_v) + @test isapprox(I_v, I_v', rtol=1e-8) + end + + # A reused buffer must not keep a stale I_v from an earlier compute_Iv=true solve + inputs = _make_3d_inputs(mtheta=32, nzeta=32, mtheta_eq=17) + vac = GeneralizedPerturbedEquilibrium.Vacuum.VacuumResponse(inputs) + compute_vacuum_response!(vac, inputs, WallShapeSettings(shape="nowall"); compute_Iv=true) + @test !all(iszero, vac.I_v) + compute_vacuum_response!(vac, inputs, WallShapeSettings(shape="nowall")) + @test all(iszero, vac.I_v) + end + + # The 3D interior operator is the 2D one shifted by the same scalar, D_int = D_ext - 2I, so an + # axisymmetric boundary driven through both paths must give the same Iᵛ. Tolerances are loose + # because Iᵛ is a difference of two solves, which amplifies the 3D toroidal discretization error + # (~8e-3 on wv here) by an order of magnitude; a wrong shift or sign gives O(1) instead. + @testset "compute_vacuum_response 3D I_v matches the 2D path" begin + mtheta = 48 + θ = range(; start=0, length=mtheta, step=2π/mtheta) + # Up-down asymmetric so Iᵛ is genuinely complex and the θ_VAC → -θ_VAC conjugation is observable + R = 1.7 .+ 0.3 .* cos.(θ) + Z = 0.3 .* sin.(θ) .+ 0.08 .* sin.(2θ) .+ 0.08 .* cos.(θ) + # Arrays are reversed for VACUUM's CW θ, as the equilibrium-based constructor does + make(nzeta) = VacuumInput(x=collect(reverse(R)), z=collect(reverse(Z)), ν=zeros(mtheta), + mtheta_in=mtheta, nzeta_in=1, m_modes=[-2, -1, 0, 1, 2], n_modes=[1], mtheta=mtheta, nzeta=nzeta) + nowall = WallShapeSettings(shape="nowall") + + r2d = compute_vacuum_response(make(1), nowall; compute_Iv=true) + r3d = compute_vacuum_response(make(mtheta), nowall; compute_Iv=true) + + @test norm(r3d.wv - r2d.wv) / norm(r2d.wv) < 2e-2 + @test norm(r3d.I_v - r2d.I_v) / norm(r2d.I_v) < 0.2 + # The imaginary parts must agree in sign, not be opposed — this is what pins the conjugation + @test norm(imag.(r3d.I_v) - imag.(r2d.I_v)) < norm(imag.(r3d.I_v) + imag.(r2d.I_v)) + end + # Field-periodic (layer-2) reduction: an nfp-periodic boundary makes the boundary-integral # operators block-circulant, so the reduced per-residue-class solve must reproduce the full # torus result (the n_stride=1 bridge case spans several residue classes mod nfp). @@ -848,6 +892,105 @@ # Hermitian part is enforced after assembly on both paths @test isapprox(wv_red, wv_red', rtol=1e-12) + + # The interior solve is a scalar shift of the exterior one, so it decomposes by the same + # residue class and Iᵛ must reduce exactly like wv does. + Iv_red = compute_vacuum_response(inputs_red, wall_settings; compute_Iv=true).I_v + Iv_full = compute_vacuum_response(inputs_full, wall_settings; compute_Iv=true).I_v + @test all(isfinite, Iv_red) + @test !all(iszero, Iv_red) + @test isapprox(Iv_red, Iv_full; rtol=1e-6, atol=1e-7) + for in1 in eachindex(n_modes), in2 in eachindex(n_modes) + if classes[in1] != classes[in2] + @test all(iszero, Iv_red[((in1-1)*mpert+1):(in1*mpert), ((in2-1)*mpert+1):(in2*mpert)]) + end + end + + # A wall adds a second source block to the operator, so repeat the check with one present: + # the field-period fold has to land the wall columns in the right block for both to agree. + walled = WallShapeSettings(shape="conformal", a=0.2, equal_arc_wall=false) + vac_red_wall = compute_vacuum_response(inputs_red, walled; compute_Iv=true) + vac_full_wall = compute_vacuum_response(inputs_full, walled; compute_Iv=true) + @test isapprox(vac_red_wall.wv, vac_full_wall.wv; rtol=1e-6, atol=1e-7) + @test isapprox(vac_red_wall.I_v, vac_full_wall.I_v; rtol=1e-6, atol=1e-7) + end + + @testset "compute_vacuum_response 3D stellarator symmetry" begin + # Rotating ellipse: R(-θ,-ζ) = R(θ,ζ) and Z(-θ,-ζ) = -Z(θ,ζ), so the surface is + # stellarator symmetric. Adding `odd` breaks that symmetry without changing anything else. + _stell_boundary(; mtheta, nzeta_p, nfp, R0=1.7, a=0.3, b=0.09, odd=0.0) = begin + X = Float64[] + Y = Float64[] + Z = Float64[] + for j in 1:nzeta_p + ζ = (j - 1) * 2π / (nzeta_p * nfp) + for i in 1:mtheta + θi = (i - 1) * 2π / mtheta + R = R0 + a * cos(θi) + b * cos(θi - nfp * ζ) + odd * sin(θi - nfp * ζ) + push!(X, R * cos(ζ)) + push!(Y, R * sin(ζ)) + push!(Z, -a * sin(θi) + b * sin(θi - nfp * ζ) + odd * cos(2θi - nfp * ζ)) + end + end + return X, Y, Z + end + _stell_inputs(; mtheta, nzeta_p, nfp, n_modes, odd=0.0) = begin + X, Y, Z = _stell_boundary(; mtheta=mtheta, nzeta_p=nzeta_p, nfp=nfp, odd=odd) + return VacuumInput( + x=X, y=Y, z=Z, + mtheta_in=mtheta, nzeta_in=nzeta_p, + m_modes=collect(-1:1), n_modes=n_modes, + mtheta=mtheta, nzeta=nzeta_p, + nfp=nfp + ) + end + + mtheta, nzeta_p = 24, 8 + nowall = WallShapeSettings(shape="nowall") + walled = WallShapeSettings(shape="conformal", a=0.2, equal_arc_wall=false) + + # The whole item rests on the operator inheriting the involution, so assert it directly. + inputs = _stell_inputs(; mtheta=mtheta, nzeta_p=nzeta_p, nfp=3, n_modes=[1]) + full = GeneralizedPerturbedEquilibrium.Vacuum.expand_field_periods(inputs) + plasma = GeneralizedPerturbedEquilibrium.Vacuum.PlasmaGeometry3D(full) + wall = GeneralizedPerturbedEquilibrium.Vacuum.WallGeometry3D(full, plasma, walled) + npts = plasma.mtheta * plasma.nzeta + σ_full = [mod1(2 - mod1(p, plasma.mtheta), plasma.mtheta) + + plasma.mtheta * (mod1(2 - ((p - 1) ÷ plasma.mtheta + 1), plasma.nzeta) - 1) for p in 1:npts] + D = zeros(npts, npts) + S = zeros(npts, npts) + GeneralizedPerturbedEquilibrium.Vacuum.compute_3D_kernel_matrices!([D], [S], plasma, plasma, 11, 20, 5, [1.0]) + @test isapprox(D[σ_full, σ_full], D; rtol=1e-9, atol=1e-9 * maximum(abs, D)) + @test isapprox(S[σ_full, σ_full], S; rtol=1e-9, atol=1e-9 * maximum(abs, S)) + + # Detection: symmetric surfaces are recognised, an odd-parity perturbation is not + @test GeneralizedPerturbedEquilibrium.Vacuum.stellarator_involution(plasma, wall, 3) !== nothing + asym = _stell_inputs(; mtheta=mtheta, nzeta_p=nzeta_p, nfp=3, n_modes=[1], odd=0.07) + asym_full = GeneralizedPerturbedEquilibrium.Vacuum.expand_field_periods(asym) + asym_plasma = GeneralizedPerturbedEquilibrium.Vacuum.PlasmaGeometry3D(asym_full) + asym_wall = GeneralizedPerturbedEquilibrium.Vacuum.WallGeometry3D(asym_full, asym_plasma, walled) + @test GeneralizedPerturbedEquilibrium.Vacuum.stellarator_involution(asym_plasma, asym_wall, 3) === nothing + + # The symmetry-adapted solve must reproduce the untransformed one. nfp = 1 and k = 0 split + # into two real half-size blocks; k ≠ 0 becomes real at full size; nfp = 4, k = 2 is the + # self-conjugate class that needs the signed involution rather than the half-twist. + for (nfp, n_modes, wall_settings) in [ + (1, [1], nowall), (1, [1], walled), + (3, [0], walled), (3, [1], nowall), (3, [1], walled), + (3, [1, 2, 3], walled), (4, [2], walled), (2, [1], walled) + ] + inp = _stell_inputs(; mtheta=mtheta, nzeta_p=nzeta_p, nfp=nfp, n_modes=n_modes) + sym = compute_vacuum_response(inp, wall_settings; compute_Iv=true, use_symmetry=true) + ref = compute_vacuum_response(inp, wall_settings; compute_Iv=true, use_symmetry=false) + @test isapprox(sym.wv, ref.wv; rtol=1e-9, atol=1e-9 * maximum(abs, ref.wv)) + @test isapprox(sym.I_v, ref.I_v; rtol=1e-9, atol=1e-9 * maximum(abs, ref.I_v)) + end + + # An asymmetric boundary must fall through to exactly the untransformed solve + sym = compute_vacuum_response(asym, walled; compute_Iv=true, use_symmetry=true) + ref = compute_vacuum_response(asym, walled; compute_Iv=true, use_symmetry=false) + @test sym.wv == ref.wv + @test sym.I_v == ref.I_v end @testset "Kernel3D laplace_kernel" begin