From ef9f2036c59fd1b5f58d927235edcceeef77287c Mon Sep 17 00:00:00 2001 From: Jake Halpern Date: Mon, 31 Aug 2026 17:32:32 -0400 Subject: [PATCH 1/2] VACUUM - IMPROVEMENT - adding 3D surface current matrix calculation + speedups with RHS multiplication by Fourier operator before inversion + incorporating the mode family phase into the kernel construction and looping through each mode family separately. In 3D, the lu! dominates the cost, so adding the extra kernel calls is outweighed by reducing peak memory and the lu size --- src/Vacuum/DataTypes.jl | 2 +- src/Vacuum/Kernel3D.jl | 46 ++++++----- src/Vacuum/Vacuum.jl | 164 ++++++++++++++++++++++------------------ test/runtests_vacuum.jl | 65 ++++++++++++++++ 4 files changed, 183 insertions(+), 94 deletions(-) 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..bcbf7bc6a 100644 --- a/src/Vacuum/Kernel3D.jl +++ b/src/Vacuum/Kernel3D.jl @@ -321,7 +321,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_greenfunction, greenfunction, observer, source, PATCH_RAD, RAD_DIM, INTERP_ORDER, phases) Compute boundary integral kernel matrices for 3D geometries with the singular correction algorithm from [Malhotra Plasma Phys. and Cont. Fusion 2019 024004]. @@ -334,7 +334,10 @@ 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 @@ -345,18 +348,21 @@ Takes advantage of field periodicity to evaluate the kernel only over a single f - `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. """ function compute_3D_kernel_matrices!( - grad_greenfunction::AbstractMatrix{Float64}, - greenfunction::AbstractMatrix{Float64}, + grad_greenfunction::AbstractMatrix{<:Number}, + greenfunction::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} ) 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 @@ -364,13 +370,15 @@ function compute_3D_kernel_matrices!( 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) + ((row_index-1)*num_points_per_fp+1):(row_index*num_points_per_fp), + ((col_index-1)*num_points_per_fp+1):(col_index*num_points_per_fp) ) # 𝒢ⁿ 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) + # All field periods fold onto the same columns, so every write below accumulates + fill!(grad_greenfunction_block, 0) + populate_greenfunction && fill!(greenfunction, 0) # 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 @@ -387,7 +395,7 @@ function compute_3D_kernel_matrices!( workspaces = [KernelWorkspace(PATCH_DIM, RAD_DIM, ANG_DIM) for _ in 1:max_threadid] # Parallel loop through observer points - Threads.@threads for idx_obs in 1:n_obs + Threads.@threads for idx_obs in 1:num_points_per_fp # Get thread-local workspace ws = workspaces[Threads.threadid()] (; r_patch, dr_dθ_patch, dr_dζ_patch, r_polar, dr_dθ_polar, dr_dζ_polar, @@ -408,11 +416,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ζ + greenfunction[idx_obs, idx_col] += phases[d] * (far_single * dθdζ) end - grad_greenfunction_block[idx_obs, idx_src] = far_double * dθdζ + grad_greenfunction_block[idx_obs, idx_col] += phases[d] * (far_double * dθdζ) end # ============================================================ @@ -458,6 +467,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,20 +476,20 @@ 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ζ + greenfunction[idx_obs, 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ζ + grad_greenfunction_block[idx_obs, 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. + # 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. grad_greenfunction_block ./= 2π populate_greenfunction && (greenfunction ./= 2π) # Add the term that comes from the volume integral of Green's identity if typeof(source) == typeof(observer) - for i in 1:n_obs + for i in 1:num_points_per_fp grad_greenfunction_block[i, i] += 1.0 end end diff --git a/src/Vacuum/Vacuum.jl b/src/Vacuum/Vacuum.jl index 9c4921618..fe4fdeced 100644 --- a/src/Vacuum/Vacuum.jl +++ b/src/Vacuum/Vacuum.jl @@ -49,6 +49,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) @@ -123,15 +141,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) @@ -160,26 +170,23 @@ end """ _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 -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`). +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. -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 `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, so it enters only the `d = 0` block and the +reduction stays exact. """ @with_pool pool function _compute_vacuum_response_3d!(vac_data::VacuumResponse, inputs::VacuumInput, wall_settings::WallShapeSettings; compute_Iv::Bool=false) @@ -187,77 +194,84 @@ interior variant `-D + 2I` for the interior columns, then scatter back into the 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) + # Work matrices, real for nfp == 1 where every field-period phase is unity + T = nfp == 1 ? Float64 : ComplexF64 + grad_green = zeros!(pool, T, n_obs, nb * num_points_per_fp) # reduced double-layer D̂ₖ (with +I on the diagonals) + green_temp = zeros!(pool, T, n_obs, num_points_per_fp) # reduced single-layer Ŝₖ (plasma sources only) + grre = zeros!(pool, ComplexF64, n_obs, mpert * length(n_modes)) + grri = compute_Iv ? similar!(pool, grre) : zeros!(pool, ComplexF64, 0, 0) + grad_green_interior = compute_Iv ? similar!(pool, grad_green) : zeros!(pool, T, 0, 0) + + # Loop over all decoupled toroidal residue classes + for k in unique(mod.(n_modes, nfp)) + 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, :] + phases = T[cis(-2π * (k * d) / nfp) for d in 0:(nfp-1)] + grre_k = @view grre[:, 1:length(mode_cols)] - 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) - end + # 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, phases) - 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] - end - scols = (d*num_points_per_fp+1):((d+1)*num_points_per_fp) - @views @. Ŝ += phase * green_temp[:, scols] + 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, phases) + # Wall–Plasma block + compute_3D_kernel_matrices!(grad_green, view(green_temp, (num_points_per_fp+1):n_obs, :), wall, plasma_surf, PATCH_RAD, RAD_DIM, INTERP_ORDER, phases) + # Wall–Wall block + compute_3D_kernel_matrices!(grad_green, view(green_temp, (num_points_per_fp+1):n_obs, :), wall, wall, PATCH_RAD, RAD_DIM, INTERP_ORDER, phases) + end + + # Project observers onto source basis exp(i(mθ-nζ)); D⁻¹ acts on rows, so this commutes with the solve + mul!(grre_k, green_temp, E') + + if compute_Iv + # Copy RHS before exterior solve overwrites grre; keep a kernel copy for interior + grri_k = @view grri[:, 1:length(mode_cols)] + grri_k .= grre_k + 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_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:n_obs + grad_green_interior[i, i] -= 2.0 end + ldiv!(lu!(grad_green_interior), grri_k) - # 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')) + _fill_Iv_block!(@view(vac_data.I_v[mode_cols, mode_cols]), E, grre_k, grri_k, num_points_per_fp) + else + # Only need exterior system for wv + ldiv!(lu!(grad_green), grre_k) end + + # Project exterior kernel onto observer basis exp(-i(mθ-nζ)) and scale to get the response matrix + mul!(wv_block, E, @view(grre_k[1:num_points_per_fp, :])) + wv_block .*= 4π^2 / num_points_per_fp end # Remove any non-Hermitian residual from Hermitian matrices due to discretization @@ -275,7 +289,7 @@ end 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`. """ function compute_vacuum_response(inputs::VacuumInput, wall_settings::WallShapeSettings; compute_Iv::Bool=false) vac = VacuumResponse(inputs) diff --git a/test/runtests_vacuum.jl b/test/runtests_vacuum.jl index 2db363573..05b94c0f3 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,27 @@ # 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 "Kernel3D laplace_kernel" begin From a61d7dde64e27631bef42fcbcab387ca31cb13dd Mon Sep 17 00:00:00 2001 From: Jake Halpern Date: Tue, 1 Sep 2026 08:14:10 -0400 Subject: [PATCH 2/2] VACUUM - IMPROVEMENT - taking advtange of stellarator symmetry. Needs cleaning up since I have no idea what this does, but Claude is reporting a large speedup and memory reduction --- src/Vacuum/Kernel3D.jl | 94 +++++++++----- src/Vacuum/Symmetry3D.jl | 260 +++++++++++++++++++++++++++++++++++++++ src/Vacuum/Vacuum.jl | 170 +++++++++++++++++-------- test/runtests_vacuum.jl | 78 ++++++++++++ 4 files changed, 520 insertions(+), 82 deletions(-) create mode 100644 src/Vacuum/Symmetry3D.jl diff --git a/src/Vacuum/Kernel3D.jl b/src/Vacuum/Kernel3D.jl index bcbf7bc6a..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, phases) + 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]. @@ -341,8 +343,9 @@ 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 @@ -350,35 +353,36 @@ full-torus source blocks are never stored. - `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{<:Number}, - greenfunction::AbstractMatrix{<:Number}, + 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, - phases::AbstractVector{<:Number} + phases::AbstractVector{<:Number}, + sym::Union{Nothing,StellaratorBasis}=nothing ) num_points = observer.mtheta * observer.nzeta 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)*num_points_per_fp+1):(row_index*num_points_per_fp), - ((col_index-1)*num_points_per_fp+1):(col_index*num_points_per_fp) - ) # 𝒢ⁿ only needed for plasma as source term (RHS of eqs. 26/27 in Chance 1997) populate_greenfunction = source isa PlasmaGeometry3D - # All field periods fold onto the same columns, so every write below accumulates - fill!(grad_greenfunction_block, 0) - populate_greenfunction && fill!(greenfunction, 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 @@ -390,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:num_points_per_fp + 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, :] @@ -419,9 +434,9 @@ function compute_3D_kernel_matrices!( # 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_col] += phases[d] * (far_single * dθdζ) + row_single[idx_col] += phases[d] * (far_single * dθdζ) end - grad_greenfunction_block[idx_obs, idx_col] += phases[d] * (far_double * dθdζ) + row_double[idx_col] += phases[d] * (far_double * dθdζ) end # ============================================================ @@ -476,21 +491,34 @@ function compute_3D_kernel_matrices!( # Apply near + far contributions if populate_greenfunction - greenfunction[idx_obs, idx_col] += phases[d] * (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_col] += phases[d] * (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 - # 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. - 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:num_points_per_fp - 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 fe4fdeced..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") @@ -131,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 @@ -168,7 +169,7 @@ 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. @@ -184,11 +185,23 @@ basis, `M = mtheta·nzeta`), so the routine loops over classes exactly as the 2D 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, so it enters only the `d = 0` block and the -reduction stays exact. +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) @@ -212,66 +225,125 @@ reduction stays exact. RAD_DIM = 20 INTERP_ORDER = 5 - # Work matrices, real for nfp == 1 where every field-period phase is unity - T = nfp == 1 ? Float64 : ComplexF64 - grad_green = zeros!(pool, T, n_obs, nb * num_points_per_fp) # reduced double-layer D̂ₖ (with +I on the diagonals) - green_temp = zeros!(pool, T, n_obs, num_points_per_fp) # reduced single-layer Ŝₖ (plasma sources only) + # 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) - grad_green_interior = compute_Iv ? similar!(pool, grad_green) : zeros!(pool, T, 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 # Loop over all decoupled toroidal residue classes - for k in unique(mod.(n_modes, nfp)) + 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, :] - phases = T[cis(-2π * (k * d) / nfp) for d in 0:(nfp-1)] - grre_k = @view grre[:, 1:length(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_green, view(green_temp, 1:num_points_per_fp, :), plasma_surf, plasma_surf, PATCH_RAD, RAD_DIM, INTERP_ORDER, phases) + 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_green, view(green_temp, 1:num_points_per_fp, :), plasma_surf, wall, PATCH_RAD, RAD_DIM, INTERP_ORDER, phases) + 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_green, view(green_temp, (num_points_per_fp+1):n_obs, :), wall, plasma_surf, PATCH_RAD, RAD_DIM, INTERP_ORDER, phases) + 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_green, view(green_temp, (num_points_per_fp+1):n_obs, :), wall, wall, PATCH_RAD, RAD_DIM, INTERP_ORDER, phases) + compute_3D_kernel_matrices!(grad_blocks, green_blocks, wall, wall, PATCH_RAD, RAD_DIM, INTERP_ORDER, phases, sym) end - # Project observers onto source basis exp(i(mθ-nζ)); D⁻¹ acts on rows, so this commutes with the solve - mul!(grre_k, green_temp, E') - - if compute_Iv - # Copy RHS before exterior solve overwrites grre; keep a kernel copy for interior - grri_k = @view grri[:, 1:length(mode_cols)] - grri_k .= grre_k - 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_k) + # 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 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:n_obs - grad_green_interior[i, i] -= 2.0 + 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 + 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 - ldiv!(lu!(grad_green_interior), grri_k) - _fill_Iv_block!(@view(vac_data.I_v[mode_cols, mode_cols]), E, grre_k, grri_k, num_points_per_fp) - else - # Only need exterior system for wv - ldiv!(lu!(grad_green), grre_k) + # 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 - - # Project exterior kernel onto observer basis exp(-i(mθ-nζ)) and scale to get the response matrix - mul!(wv_block, E, @view(grre_k[1:num_points_per_fp, :])) 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 # Remove any non-Hermitian residual from Hermitian matrices due to discretization @@ -284,31 +356,31 @@ reduction stays exact. 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`. +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_vacuum.jl b/test/runtests_vacuum.jl index 05b94c0f3..d93bc6da7 100644 --- a/test/runtests_vacuum.jl +++ b/test/runtests_vacuum.jl @@ -915,6 +915,84 @@ @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 G, K = GeneralizedPerturbedEquilibrium.Vacuum.laplace_kernel(1.0, 0.0, 0.0, 2.0, 0.0, 0.0, 1.0, 0.0, 0.0) # Kernel returns 1/|r_obs - r_src| (4π factor applied elsewhere in BIE)