diff --git a/src/Vacuum/DataTypes.jl b/src/Vacuum/DataTypes.jl index a15f08e330..c6cbc18c65 100644 --- a/src/Vacuum/DataTypes.jl +++ b/src/Vacuum/DataTypes.jl @@ -417,7 +417,7 @@ end Build the 3D plasma surface on the `mtheta × nzeta` grid in (θ, ζ). An axisymmetric boundary (`nzeta_in == 1`) is revolved from its 2D contour, placing each point at ϕ = ζ + ν(θ); a 3D boundary is interpolated from its input grid. -Tangents come from periodic bicubic splines and normals point into the plasma. +Tangents are those of the trigonometric interpolant of the grid samples, and normals point into the plasma. # Arguments @@ -439,8 +439,6 @@ function PlasmaGeometry3D(inputs::VacuumInput) # Allocate output arrays r = zeros(num_points, 3) - dr_dθ = zeros(num_points, 3) - dr_dζ = zeros(num_points, 3) normal = zeros(num_points, 3) if inputs.nzeta_in == 1 @@ -468,15 +466,9 @@ function PlasmaGeometry3D(inputs::VacuumInput) assert_counterclockwise(hypot.(r[1:mtheta, 1], r[1:mtheta, 2]), r[1:mtheta, 3], "Plasma boundary ζ = 0 cross-section") end - # Compute tangent vectors and normal vectors via periodic bicubic splines - itps = [cubic_interp((θ_grid, ζ_grid), reshape(r[:, k], mtheta, nzeta); bc=(PeriodicBC(; endpoint=:exclusive), PeriodicBC(; endpoint=:exclusive))) for k in 1:3] - for i in 1:mtheta, j in 1:nzeta - idx = i + (j - 1) * mtheta - for k in 1:3 - partials = itps[k].nodal_derivs.partials - dr_dθ[idx, k] = partials[2, i, j] - dr_dζ[idx, k] = partials[3, i, j] - end + # Tangent and normal vectors of the trigonometric interpolant, which converge spectrally with the grid + dr_dθ, dr_dζ = trig_tangents(r, mtheta, nzeta) + for idx in 1:num_points cross3!(normal, dr_dθ, dr_dζ, idx) end @@ -485,12 +477,14 @@ function PlasmaGeometry3D(inputs::VacuumInput) normal_orient = normal[idx, 1] < 0 ? 1 : -1 normal .*= normal_orient - # Warn if grid spacing is highly anisotropic + # Warn for unresolved grids spacing_θ = sqrt(sum(abs2, dr_dθ) / size(dr_dθ, 1)) * dθ spacing_ζ = sqrt(sum(abs2, dr_dζ) / size(dr_dζ, 1)) * dζ aspect_ratio = spacing_ζ / spacing_θ @info "Average grid spacing [m]: dθ=$(round(spacing_θ, digits=4)), dζ=$(round(spacing_ζ, digits=4)), aspect ratio=$(round(aspect_ratio, digits=2))" - aspect_ratio > 10.0 && @warn "Grid aspect ratio is highly anisotropic, which may degrade quadrature accuracy" + inputs.nzeta_in > 1 && aspect_ratio > 1.75 && + @warn "Toroidal grid spacing is $(round(aspect_ratio, digits=2))× the poloidal one, which degrades the quadrature accuracy. " * + "Raise nzeta by about $(round(aspect_ratio, digits=1))× for equal spacing." return PlasmaGeometry3D( mtheta, @@ -731,7 +725,6 @@ function WallGeometry3D(inputs::VacuumInput, plasma_surf::PlasmaGeometry3D, wall dθ = 2π / mtheta dζ = 2π / nzeta - θ_grid = range(; start=0, length=mtheta, step=dθ) ϕ_grid = range(; start=0, length=nzeta, step=dζ) num_points = mtheta * nzeta @@ -798,15 +791,9 @@ function WallGeometry3D(inputs::VacuumInput, plasma_surf::PlasmaGeometry3D, wall error("Wall shape $(nameof(typeof(wall_settings))) is not available for a non-axisymmetric boundary (nzeta_in > 1).") end - # Compute tangent vectors and normal vectors via periodic bicubic splines - itps = [cubic_interp((θ_grid, ϕ_grid), reshape(r[:, k], mtheta, nzeta); bc=(PeriodicBC(; endpoint=:exclusive), PeriodicBC(; endpoint=:exclusive))) for k in 1:3] - for i in 1:mtheta, j in 1:nzeta - idx = i + (j - 1) * mtheta - for k in 1:3 - partials = itps[k].nodal_derivs.partials - dr_dθ[idx, k] = partials[2, i, j] - dr_dζ[idx, k] = partials[3, i, j] - end + # Tangent and normal vectors of the trigonometric interpolant, as on the plasma + dr_dθ, dr_dζ = trig_tangents(r, mtheta, nzeta) + for idx in 1:num_points cross3!(normal, dr_dθ, dr_dζ, idx) end @@ -835,6 +822,12 @@ function WallGeometry3D(inputs::VacuumInput, plasma_surf::PlasmaGeometry3D, wall ) (min_align < 0.5 || min_area_ratio < 0.5) && @warn "Conformal wall is close to folding: min(n̂_wall·n̂_plasma)=$((@sprintf "%.3f" min_align)), min area-element ratio=$((@sprintf "%.3f" min_area_ratio)). Reduce the wall distance a." + # An exact normal offset has n̂_wall ∥ n̂_plasma, so their angle is the interpolation error of the wall surface: + # a sharp convex plasma edge fans out over the gap faster than the grid follows it. 2° ≈ 1e-4 relative error in Wᵛ. + misalign = acosd(clamp(min_align, -1.0, 1.0)) + min_align >= 0.5 && misalign > 2.0 && + @warn "Conformal wall is under-resolved by the $(mtheta)×$(nzeta) grid: wall and plasma normals differ by up to $((@sprintf "%.1f" misalign))° " * + "(0° when resolved), where a sharp convex edge of the plasma fans out across the gap. Refine mtheta and nzeta until this is below 2°, or reduce the wall distance a." end # The 3D singular quadrature resolves the double-layer kernel well only while the plasma-wall gap spans a diff --git a/src/Vacuum/Kernel3D.jl b/src/Vacuum/Kernel3D.jl index 345671fe3f..15958da35d 100644 --- a/src/Vacuum/Kernel3D.jl +++ b/src/Vacuum/Kernel3D.jl @@ -8,18 +8,52 @@ # Original implementation: Dhairya Malhotra # # Modifications Copyright (C) 2026 Columbia Fusion Research Center and Seoul National University: -# - Julia translation; runtime patch, quadrature and interpolation orders (default Lagrange -# order 5), with the patch clamped on coarse grids +# - Julia translation; runtime patch, quadrature and interpolation orders, with the patch clamped on coarse grids # - stencil corner by rounding, so odd orders are equivariant under the patch's pi-rotation # - one sparse interpolation matrix, applied forward to the geometry and transposed to scatter weights # - VACUUM kernel convention: 1/r and n.grad'(1/r) over 2pi, unit Green's-identity diagonal # - far field by the trapezoidal rule on the collocation grid, without upsampling +# - patch geometry upsampled 2x by trigonometric interpolation to refine the polar interpolation stencil, +# correction weights folded back to the collocation nodes by its transpose (BIEST instead upsamples the +# whole surface and density and shrinks the patch) # - singular correction applied to plasma-wall cross-surface blocks # - dense operator row assembly, field-period folding, stellarator-symmetric basis, threading # # SPDX-License-Identifier: Apache-2.0 # Full text in LICENSES/Apache-2.0.txt; attribution in NOTICE. +""" +Thread-local scratch for `compute_3D_kernel_matrices!`, one per thread so the parallel observer loop does not race. +`patch` and `polar` hold (r, ∂r/∂θ, ∂r/∂ζ) per node; `T` is the operator row's element type. +""" +struct KernelWorkspace{T<:Number} + patch::Matrix{SVector{9,Float64}} + polar::Matrix{SVector{9,Float64}} + M_polar_single::Matrix{Float64} + M_polar_double::Matrix{Float64} + M_grid_flat::Vector{Float64} + fold_tmp::Matrix{Float64} + folded_patch::Matrix{T} + row_single::Vector{T} + row_double::Vector{T} + sym_work::Vector{ComplexF64} +end + +function KernelWorkspace{T}(UP_DIM::Int, RAD_DIM::Int, ANG_DIM::Int, mtheta::Int, nzeta_fp::Int, nrow::Int, nsym::Int) where {T} + return KernelWorkspace{T}( + zeros(SVector{9,Float64}, UP_DIM, UP_DIM), + zeros(SVector{9,Float64}, RAD_DIM, ANG_DIM), + zeros(RAD_DIM, ANG_DIM), + zeros(RAD_DIM, ANG_DIM), + zeros(UP_DIM^2), + zeros(UP_DIM, mtheta), + zeros(T, nzeta_fp, mtheta), + zeros(T, nrow), + zeros(T, nrow), + zeros(ComplexF64, nsym) + ) +end + """ SingularQuadratureData @@ -29,8 +63,9 @@ Precomputed polar singular-correction quadrature (Malhotra JCP 397 (2019) 108791 - `Gpou::Matrix{Float64}`: Partition of unity on the Cartesian patch, `Gpou = -χ` (PATCH_DIM × PATCH_DIM) - `Ppou::Matrix{Float64}`: Partition of unity times polar quadrature weights (RAD_DIM × ANG_DIM) - - `P2G::SparseMatrixCSC{Float64,Int}`: Lagrange interpolation matrix (Ngrid × Npolar): `polar = P2G' * patch`, `grid = P2G * polar` + - `P2G::SparseMatrixCSC{Float64,Int}`: Lagrange interpolation matrix on the 2× upsampled patch: `polar = P2G' * patch`, `grid = P2G * polar` - `PATCH_DIM::Int`: Patch width `2 * PATCH_RAD + 1` + - `UP_DIM::Int`: Upsampled patch width `2 * PATCH_DIM - 1` - `PATCH_RAD::Int`: Number of points on each side of the source point treated as singular - `ANG_DIM::Int`: Number of angular quadrature points - `RAD_DIM::Int`: Number of radial quadrature points @@ -41,6 +76,7 @@ struct SingularQuadratureData Ppou::Matrix{Float64} P2G::SparseMatrixCSC{Float64,Int} PATCH_DIM::Int + UP_DIM::Int PATCH_RAD::Int ANG_DIM::Int RAD_DIM::Int @@ -51,7 +87,7 @@ end SingularQuadratureData(PATCH_RAD, RAD_DIM, INTERP_ORDER) Build the polar quadrature, partitions of unity, and Lagrange interpolant for the singular patch. -`ANG_DIM = 2 * RAD_DIM`. `INTERP_ORDER` must be `≤ 2 * PATCH_RAD + 1`. +`ANG_DIM = 2 * RAD_DIM`. `INTERP_ORDER` must be `≤ 4 * PATCH_RAD + 1`, the upsampled patch width. # Arguments @@ -67,7 +103,9 @@ function SingularQuadratureData(PATCH_RAD::Int, RAD_DIM::Int, INTERP_ORDER::Int) # Total size of square patch extracted around singular point (odd number: 2*PATCH_RAD+1) PATCH_DIM = 2 * PATCH_RAD + 1 - @assert INTERP_ORDER <= PATCH_DIM "Must have INTERP_ORDER <= PATCH_DIM, got INTERP_ORDER=$INTERP_ORDER, PATCH_DIM=$PATCH_DIM" + # We upsample the patch by 2x in each direction for better accuracy + UP_DIM = 2 * PATCH_DIM - 1 + @assert INTERP_ORDER <= UP_DIM "Must have INTERP_ORDER <= UP_DIM, got INTERP_ORDER=$INTERP_ORDER, UP_DIM=$UP_DIM" # Number of angular quadrature nodes in polar coordinates (uniformly distributed around circle) ANG_DIM = 2 * RAD_DIM @@ -76,7 +114,7 @@ function SingularQuadratureData(PATCH_RAD::Int, RAD_DIM::Int, INTERP_ORDER::Int) qx = (qx_raw .+ 1) ./ 2 # [-1, 1] → [0, 1] qw = qw_raw ./ 2 # Adjust weights for interval change - # χ(r) = exp(-36 r^p), p from PATCH_DIM + # χ(r) = exp(-36 r^p), p from PATCH_DIM: χ blends on the collocation grid, so it keys on the coarse patch, not UP_DIM pou_power = PATCH_DIM > 45 ? 10 : (PATCH_DIM > 20 ? 8 : 6) pou(r) = r ≥ 1.0 ? 0.0 : exp(-36.0 * r^pou_power) @@ -115,7 +153,7 @@ function SingularQuadratureData(PATCH_RAD::Int, RAD_DIM::Int, INTERP_ORDER::Int) end # grid = P2G * polar, polar = P2G' * grid; each column is the INTERP_ORDER² Lagrange stencil - Ngrid = PATCH_DIM * PATCH_DIM + Ngrid = UP_DIM * UP_DIM Npolar = RAD_DIM * ANG_DIM # Preallocate COO storage: @@ -134,12 +172,12 @@ function SingularQuadratureData(PATCH_RAD::Int, RAD_DIM::Int, INTERP_ORDER::Int) # Lower-left corner indices of INTERP_ORDER × INTERP_ORDER stencil centered on (x0,x1). # Round, don't truncate, and keep INTERP_ORDER odd: the stencil must be equivariant under the patch's π-rotation or stellarator symmetry is lost. - 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) + y0 = clamp(round(Int, x0 * (UP_DIM - 1)) - (INTERP_ORDER - 1) ÷ 2, 0, UP_DIM - INTERP_ORDER) + y1 = clamp(round(Int, x1 * (UP_DIM - 1)) - (INTERP_ORDER - 1) ÷ 2, 0, UP_DIM - INTERP_ORDER) # Local coordinates within INTERP_ORDER×INTERP_ORDER stencil, normalized to [0,1] - z0 = (x0 * (PATCH_DIM - 1) - y0) * h - z1 = (x1 * (PATCH_DIM - 1) - y1) * h + z0 = (x0 * (UP_DIM - 1) - y0) * h + z1 = (x1 * (UP_DIM - 1) - y1) * h # Polar point index (column in P2G) j_polar = ir + RAD_DIM * (ia - 1) @@ -147,7 +185,7 @@ function SingularQuadratureData(PATCH_RAD::Int, RAD_DIM::Int, INTERP_ORDER::Int) # Populate stencil contributions for this polar node for i0 in 1:INTERP_ORDER, i1 in 1:INTERP_ORDER # Grid point index (row in P2G), using column-major layout - i_grid = (y0 + i0) + PATCH_DIM * (y1 + i1 - 1) + i_grid = (y0 + i0) + UP_DIM * (y1 + i1 - 1) I_coo[idx] = i_grid J_coo[idx] = j_polar V_coo[idx] = lagrange_interp(z0, z1, i0 - 1, i1 - 1) @@ -158,7 +196,7 @@ function SingularQuadratureData(PATCH_RAD::Int, RAD_DIM::Int, INTERP_ORDER::Int) # Assemble sparse interpolation matrix P2G = sparse(I_coo, J_coo, V_coo, Ngrid, Npolar) - return SingularQuadratureData(Gpou, Ppou, P2G, PATCH_DIM, PATCH_RAD, ANG_DIM, RAD_DIM, INTERP_ORDER) + return SingularQuadratureData(Gpou, Ppou, P2G, PATCH_DIM, UP_DIM, PATCH_RAD, ANG_DIM, RAD_DIM, INTERP_ORDER) end # Cached singular quadrature; rebuilt if PATCH_RAD / RAD_DIM / INTERP_ORDER change @@ -211,77 +249,70 @@ end """ Extract a periodically wrapped `PATCH_DIM × PATCH_DIM` patch centered at `(idx_pol_center, idx_tor_center)`. """ -function extract_patch!(patch::Array{Float64,3}, data::Matrix{Float64}, idx_pol_center::Int, idx_tor_center::Int, npol::Int, ntor::Int, PATCH_DIM::Int) +function extract_patch!(patch::AbstractMatrix, data::AbstractVector, idx_pol_center::Int, idx_tor_center::Int, npol::Int, ntor::Int, PATCH_DIM::Int) PATCH_RAD = (PATCH_DIM - 1) ÷ 2 @inbounds for j in 1:PATCH_DIM, i in 1:PATCH_DIM idx_pol = periodic_wrap(idx_pol_center - PATCH_RAD + i - 1, npol) idx_tor = periodic_wrap(idx_tor_center - PATCH_RAD + j - 1, ntor) - idx_src = idx_pol + npol * (idx_tor - 1) - patch[i, j, 1] = data[idx_src, 1] - patch[i, j, 2] = data[idx_src, 2] - patch[i, j, 3] = data[idx_src, 3] + patch[i, j] = data[idx_pol+npol*(idx_tor-1)] end end """ -Interpolate a Cartesian patch onto polar quadrature nodes: `polar = P2G' * patch`. -""" -function interpolate_to_polar!(polar_data::Array{Float64,3}, patch::Array{Float64,3}, P2G::SparseMatrixCSC{Float64,Int}) - patch_flat = reshape(patch, :, size(patch, 3)) - mul!(reshape(polar_data, :, size(patch, 3)), P2G', patch_flat) -end - +Interpolate the patch onto the polar nodes, `polar = P2G' * patch`. Each node carries all nine fields as one `SVector`, +so `P2G` is swept once rather than once per field as `mul!` would. """ - compute_polar_normal!(n_polar, dr_dθ, dr_dζ, normal_orient) - -`n = ∂r/∂θ × ∂r/∂ζ` at the polar nodes. Re-apply `normal_orient`: these normals are rebuilt from -interpolated tangents, so they do not inherit the stored surface orientation. -""" -function compute_polar_normal!(n_polar::Array{Float64,3}, dr_dθ::Array{Float64,3}, dr_dζ::Array{Float64,3}, normal_orient::Int) - @inbounds for ia in axes(dr_dθ, 2), ir in axes(dr_dθ, 1) - n_polar[ir, ia, 1] = dr_dθ[ir, ia, 2] * dr_dζ[ir, ia, 3] - dr_dθ[ir, ia, 3] * dr_dζ[ir, ia, 2] - n_polar[ir, ia, 2] = dr_dθ[ir, ia, 3] * dr_dζ[ir, ia, 1] - dr_dθ[ir, ia, 1] * dr_dζ[ir, ia, 3] - n_polar[ir, ia, 3] = dr_dθ[ir, ia, 1] * dr_dζ[ir, ia, 2] - dr_dθ[ir, ia, 2] * dr_dζ[ir, ia, 1] +function interpolate_to_polar!(polar::Matrix{SVector{9,Float64}}, patch::Matrix{SVector{9,Float64}}, P2G::SparseMatrixCSC{Float64,Int}) + rv, nz = rowvals(P2G), nonzeros(P2G) + @inbounds for c in axes(P2G, 2) + # Accumulates (r, ∂r/∂θ, ∂r/∂ζ) at polar node c, no allocations + acc = zero(SVector{9,Float64}) + for p in nzrange(P2G, c) + # Lagrange weight × all nine fields at that patch node + acc += nz[p] * patch[rv[p]] + end + polar[c] = acc end - n_polar .*= normal_orient end """ -Thread-local scratch for `compute_3D_kernel_matrices!`. One instance per thread so the parallel observer loop does not race. +Add the polar correction weights `M_polar` to the operator row: `P2G` scatters them onto the upsampled patch, and +`Uζ_patch * M_up' * Uθ_patch` (the transposed trigonometric upsampling about the first node, field periods folded into +`Uζ_patch`) onto the collocation nodes, shifted to the observer `(i_obs, j_obs)` by translation invariance. """ -struct KernelWorkspace - r_patch::Array{Float64,3} - dr_dθ_patch::Array{Float64,3} - dr_dζ_patch::Array{Float64,3} - r_polar::Array{Float64,3} - dr_dθ_polar::Array{Float64,3} - dr_dζ_polar::Array{Float64,3} - n_polar::Array{Float64,3} - M_polar_single::Matrix{Float64} - M_polar_double::Matrix{Float64} - M_grid_single_flat::Vector{Float64} - M_grid_double_flat::Vector{Float64} -end - -""" - KernelWorkspace(PATCH_DIM, RAD_DIM, ANG_DIM) - -Preallocated patch, polar, and grid buffers for one observer row. -""" -function KernelWorkspace(PATCH_DIM::Int, RAD_DIM::Int, ANG_DIM::Int) - return KernelWorkspace( - zeros(PATCH_DIM, PATCH_DIM, 3), - zeros(PATCH_DIM, PATCH_DIM, 3), - zeros(PATCH_DIM, PATCH_DIM, 3), - zeros(RAD_DIM, ANG_DIM, 3), - zeros(RAD_DIM, ANG_DIM, 3), - zeros(RAD_DIM, ANG_DIM, 3), - zeros(RAD_DIM, ANG_DIM, 3), - zeros(RAD_DIM, ANG_DIM), - zeros(RAD_DIM, ANG_DIM), - zeros(PATCH_DIM^2), - zeros(PATCH_DIM^2) - ) +function scatter_polar_to_row!( + row::Vector{T}, + M_polar::Matrix{Float64}, + P2G::SparseMatrixCSC{Float64,Int}, + Uθ_patch::Matrix{Float64}, + Uζ_patch::Matrix{T}, + ws::KernelWorkspace{T}, + i_obs::Int, + j_obs::Int, + wrap_phase::T +) where {T<:Number} + + (; M_grid_flat, fold_tmp, folded_patch) = ws + UP_DIM = size(Uθ_patch, 1) + + # Scatter the polar corrections back to the upsampled patch + mul!(M_grid_flat, P2G, vec(M_polar)) + + # Fold the upsampled patch back to the collocation nodes + mul!(fold_tmp, transpose(reshape(M_grid_flat, UP_DIM, UP_DIM)), Uθ_patch) + mul!(folded_patch, Uζ_patch, fold_tmp) + + nzeta_fp, mtheta = size(folded_patch) + @inbounds for j in 1:nzeta_fp + # Columns the shift pushes into the next field period fold back with one more phase step; + # this assumes geometric phases, phases[d] = wrap_phase^(d-1) + j_fold = j - j_obs + 1 + phase = j_fold < 1 ? wrap_phase : one(T) + j_fold < 1 && (j_fold += nzeta_fp) + for i in 1:mtheta + row[i+mtheta*(j-1)] += phase * folded_patch[j_fold, periodic_wrap(i - i_obs + 1, mtheta)] + end + end end """ @@ -289,7 +320,8 @@ end 3D single- and double-layer kernels with Malhotra's polar singular correction (PPCF 62 (2020) 024004). -Far field: trapezoidal rule. Near field: polar quadrature blended by a partition of unity. +Far field: trapezoidal rule. Near field: polar quadrature blended by a partition of unity, on source data +upsampled 2× by trigonometric interpolation. A source in field period `d` is accumulated onto the period-0 column with weight `phases[d+1]`, so the block-circulant reduction is written in and the full-torus source blocks are never stored. @@ -310,6 +342,7 @@ function compute_3D_kernel_matrices!( phases::AbstractVector{<:Number}, sym_basis::Union{Nothing,StellSymBasis}=nothing ) + @assert (observer.mtheta, observer.nzeta) == (source.mtheta, source.nzeta) "observer and source must share one grid" num_points = observer.mtheta * observer.nzeta @assert num_points % length(phases) == 0 "full-torus grid ($num_points points) is not divisible by nfp ($(length(phases)))" num_points_per_fp = num_points ÷ length(phases) # observer/source points in one field period @@ -334,25 +367,35 @@ function compute_3D_kernel_matrices!( end quad_data = get_singular_quadrature(patch_rad, RAD_DIM, INTERP_ORDER) - (; PATCH_DIM, ANG_DIM, Ppou, Gpou, P2G) = quad_data + (; PATCH_DIM, UP_DIM, ANG_DIM, Ppou, Gpou, P2G) = quad_data - # One operator row is accumulated per observer, 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] + # Source r, ∂r/∂θ, ∂r/∂ζ on the 2× upsampled grid, where node 2i-1 is collocation node i + Uθ, Uζ = trig_upsample_matrix(source.mtheta), trig_upsample_matrix(source.nzeta) + src_up = SVector{9}.(eachrow(trig_upsample(hcat(source.r, source.dr_dθ, source.dr_dζ), Uθ, Uζ, source.mtheta, source.nzeta))) + + # Upsampling restricted to the patch about the first node, with the field periods folded onto period 0 in ζ (see `scatter_polar_to_row!`) 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] - # `store_kernel_row!` needs two rows of scratch: the reconstructed partner and their combination - sym_work = [zeros(ComplexF64, sym_basis === nothing ? 0 : 2 * num_points_per_fp) for _ in 1:max_threadid] + nzeta_fp = source.nzeta ÷ length(phases) + Uθ_patch = Uθ[periodic_wrap.((1:UP_DIM) .- 2 * patch_rad, 2 * source.mtheta), :] + Uζ_patch = zeros(Trow, nzeta_fp, UP_DIM) + for i in 1:UP_DIM, d in eachindex(phases), j in 1:nzeta_fp + Uζ_patch[j, i] += phases[d] * Uζ[periodic_wrap(i - 2 * patch_rad, 2 * source.nzeta), j+(d-1)*nzeta_fp] + end + wrap_phase = phases[min(2, end)] + + # Source point and normal packed per node, and each node's (field period, period-0 column), read on every observer row + src = [SVector{6}(source.r[k, 1], source.r[k, 2], source.r[k, 3], source.normal[k, 1], source.normal[k, 2], source.normal[k, 3]) for k in 1:num_points] + period_col = [fldmod1(k, num_points_per_fp) for k in 1:num_points] + + # One operator row is accumulated per observer, so the field-period fold and the basis change both happen before anything is stored. + # `store_kernel_row!` needs two rows of scratch: the reconstructed partner and their combination. + nsym = sym_basis === nothing ? 0 : 2 * num_points_per_fp + workspaces = [KernelWorkspace{Trow}(UP_DIM, RAD_DIM, ANG_DIM, source.mtheta, nzeta_fp, num_points_per_fp, nsym) for _ in 1:Threads.maxthreadid()] # :static pins each task to its thread, so the threadid()-indexed scratch below stays private Threads.@threads :static for i_pair in eachindex(observers) - 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] + ws = workspaces[Threads.threadid()] + (; patch, polar, M_polar_single, M_polar_double, row_single, row_double, sym_work) = ws fill!(row_double, 0) populate_greenfunction && fill!(row_single, 0) @@ -360,76 +403,63 @@ function compute_3D_kernel_matrices!( idx_obs = observers[i_pair] i_obs = mod1(idx_obs, observer.mtheta) j_obs = (idx_obs - 1) ÷ observer.mtheta + 1 - r_obs = @view observer.r[idx_obs, :] + r_obs = SVector(observer.r[idx_obs, 1], observer.r[idx_obs, 2], observer.r[idx_obs, 3]) # Far field: trapezoidal rule over all sources (the kernels return zero at r_src = r_obs) @inbounds for idx_src in 1:num_points # Evaluate kernels at grid points - r_src = @view source.r[idx_src, :] - 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]) + s = src[idx_src] + far_single, far_double = laplace_kernel(r_obs[1], r_obs[2], r_obs[3], s[1], s[2], s[3], s[4], s[5], s[6]) # Periodic trapezoidal rule (constant weights); fold this source's field period onto period 0 - d, idx_col = fldmod1(idx_src, num_points_per_fp) + d, idx_col = period_col[idx_src] populate_greenfunction && (row_single[idx_col] += phases[d] * (far_single * dθdζ)) row_double[idx_col] += phases[d] * (far_double * dθdζ) end - # Near field: polar quadrature with singular correction. Extract patches of source data around the singular point (size = PATCH_DIM x PATCH_DIM x dof) - extract_patch!(r_patch, source.r, i_obs, j_obs, source.mtheta, source.nzeta, PATCH_DIM) - extract_patch!(dr_dθ_patch, source.dr_dθ, i_obs, j_obs, source.mtheta, source.nzeta, PATCH_DIM) - extract_patch!(dr_dζ_patch, source.dr_dζ, i_obs, j_obs, source.mtheta, source.nzeta, PATCH_DIM) + # Near field: polar quadrature with singular correction, on the upsampled patch about the observer + extract_patch!(patch, src_up, 2 * i_obs - 1, 2 * j_obs - 1, 2 * source.mtheta, 2 * source.nzeta, UP_DIM) - # Interpolate coordinates and tangent vectors to polar quadrature points - interpolate_to_polar!(r_polar, r_patch, P2G) - interpolate_to_polar!(dr_dθ_polar, dr_dθ_patch, P2G) - interpolate_to_polar!(dr_dζ_polar, dr_dζ_patch, P2G) - - # Compute normal vectors at polar points from interpolated tangent vectors - compute_polar_normal!(n_polar, dr_dθ_polar, dr_dζ_polar, source.normal_orient) + # Interpolate the geometry onto the polar nodes + interpolate_to_polar!(polar, patch, P2G) # Evaluate kernels and apply quadrature weights: area element × POU, where POU contains rdrdθ already @inbounds for ia in 1:ANG_DIM, ir in 1:RAD_DIM - r_src = @view r_polar[ir, ia, :] - n_src = @view n_polar[ir, ia, :] - near_single, near_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]) + q = polar[ir, ia] # = (r, ∂r/∂θ, ∂r/∂ζ) at this polar node + # n = ∂r/∂θ × ∂r/∂ζ from the interpolated tangents, need to reapply orientation + n = source.normal_orient * cross(q[SVector(4, 5, 6)], q[SVector(7, 8, 9)]) + near_single, near_double = laplace_kernel(r_obs[1], r_obs[2], r_obs[3], q[1], q[2], q[3], n[1], n[2], n[3]) M_polar_single[ir, ia] = near_single * Ppou[ir, ia] * dθdζ M_polar_double[ir, ia] = near_double * Ppou[ir, ia] * dθdζ end - # Scatter the polar corrections back to the Cartesian patch: grid = P2G * polar - mul!(M_grid_double_flat, P2G, vec(M_polar_double)) - M_grid_double = reshape(M_grid_double_flat, PATCH_DIM, PATCH_DIM) - if populate_greenfunction - mul!(M_grid_single_flat, P2G, vec(M_polar_single)) - M_grid_single = reshape(M_grid_single_flat, PATCH_DIM, PATCH_DIM) - end + # Scatter the singular correction on the polar nodes back to the collocation nodes + scatter_polar_to_row!(row_double, M_polar_double, P2G, Uθ_patch, Uζ_patch, ws, i_obs, j_obs, wrap_phase) + populate_greenfunction && scatter_polar_to_row!(row_single, M_polar_single, P2G, Uθ_patch, Uζ_patch, ws, i_obs, j_obs, wrap_phase) - # POU correction: singular correction + (1 + Gpou) * far-field terms + # Remove the trapezoidal far field where the polar correction applies: Gpou = -χ, net weight is (1 - χ) @inbounds for j in 1:PATCH_DIM, i in 1:PATCH_DIM # Map back to global indices 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) + d, idx_col = period_col[idx_src] - # Remainder of far-field contribution on the singular grid: Gpou = -χ - r_src = @view source.r[idx_src, :] - 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]) + # Far field at this patch node (weighted by Gpou = −χ) + s = src[idx_src] + far_single, far_double = laplace_kernel(r_obs[1], r_obs[2], r_obs[3], s[1], s[2], s[3], s[4], s[5], s[6]) - # Apply near + far contributions - populate_greenfunction && (row_single[idx_col] += phases[d] * (M_grid_single[i, j] + far_single * Gpou[i, j] * dθdζ)) - row_double[idx_col] += phases[d] * (M_grid_double[i, j] + far_double * Gpou[i, j] * dθdζ) + populate_greenfunction && (row_single[idx_col] += phases[d] * (far_single * Gpou[i, j] * dθdζ)) + row_double[idx_col] += phases[d] * (far_double * Gpou[i, j] * dθdζ) end # Normalize so the Green's-identity term below is a unit shift and the exterior/interior jump is 2I, as in the 2D kernel row_double ./= 2π populate_greenfunction && (row_single ./= 2π) - store_kernel_row!(grad_blocks, sym_basis, row_double, sym_work[tid], i_pair, row_index, col_index) - populate_greenfunction && store_kernel_row!(green_blocks, sym_basis, row_single, sym_work[tid], i_pair, row_index, 1) + store_kernel_row!(grad_blocks, sym_basis, row_double, sym_work, i_pair, row_index, col_index) + populate_greenfunction && store_kernel_row!(green_blocks, sym_basis, row_single, sym_work, i_pair, row_index, 1) end # Volume-integral term of Green's identity; invariant under the basis change, so it is the same unit shift on every block diff --git a/src/Vacuum/Utilities.jl b/src/Vacuum/Utilities.jl index 906b9819cf..980116d7a4 100644 --- a/src/Vacuum/Utilities.jl +++ b/src/Vacuum/Utilities.jl @@ -82,6 +82,79 @@ function distribute_to_equal_arc_grid(xin::Vector{Float64}, zin::Vector{Float64} return xout, zout end +""" + trig_derivative_matrix(n) -> Matrix{Float64} + +Derivative at the nodes of the trigonometric interpolant on `n` equispaced periodic nodes. +For even `n` the Nyquist mode is a cosine, whose derivative vanishes at the nodes. +""" +function trig_derivative_matrix(n::Int) + D = zeros(n, n) + for j in 1:n, i in 1:n + r = mod(i - j, n) + r == 0 && continue + # Symmetric range keeps sin(πr/n) away from π, where it loses relative precision; the sign factor absorbs the shift + r > n ÷ 2 && (r -= n) + s, c = sincospi(r / n) + D[i, j] = (isodd(r) ? -0.5 : 0.5) * (iseven(n) ? c / s : 1 / s) + end + return D +end + +""" + trig_tangents(r, mtheta, nzeta) -> (dr_dθ, dr_dζ) + +Tangents of the trigonometric interpolant of the `(mtheta * nzeta, 3)` grid samples `r`. +""" +function trig_tangents(r::Matrix{Float64}, mtheta::Int, nzeta::Int) + Dθ, Dζ = trig_derivative_matrix(mtheta), trig_derivative_matrix(nzeta) + dr_dθ, dr_dζ = similar(r), similar(r) + for k in 1:3 + R = reshape(r[:, k], mtheta, nzeta) + dr_dθ[:, k] = vec(Dθ * R) + dr_dζ[:, k] = vec(R * transpose(Dζ)) + end + return dr_dθ, dr_dζ +end + +""" + trig_upsample_matrix(n) -> Matrix{Float64} + +Values at the `2n` nodes of the twice finer grid of the trigonometric interpolant on `n` equispaced periodic nodes. +Odd rows are the identity. +""" +function trig_upsample_matrix(n::Int) + U = zeros(2n, n) + for j in 1:n, i in 1:2n + # Offset of fine node i from coarse node j, in fine spacings + r = mod(i - 2j + 1, 2n) + if iseven(r) + U[i, j] = r == 0 + else + # Symmetric range, as in `trig_derivative_matrix`; mod(r, 4) flips with the shift exactly when csc does (odd n) + r > n && (r -= 2n) + s, c = sincospi(r / 2n) + # Periodic sinc at a midpoint: sign alternates every other midpoint; cot for even n (Nyquist cosine), csc for odd n + U[i, j] = (mod(r, 4) == 1 ? 1 : -1) * (iseven(n) ? c / s : 1 / s) / n + end + end + return U +end + +""" + trig_upsample(a, Uθ, Uζ, mtheta, nzeta) -> Matrix{Float64} + +Columns of the `(mtheta * nzeta, k)` grid samples `a` on the 2× finer grid; `Uθ`, `Uζ` from [`trig_upsample_matrix`](@ref). +""" +function trig_upsample(a::AbstractMatrix{Float64}, Uθ::Matrix{Float64}, Uζ::Matrix{Float64}, mtheta::Int, nzeta::Int) + # 2× in θ and 2× in ζ → 4× the nodes, still θ-fastest + a_up = Matrix{Float64}(undef, 4 * mtheta * nzeta, size(a, 2)) + for k in axes(a, 2) + a_up[:, k] = vec(Uθ * reshape(a[:, k], mtheta, nzeta) * transpose(Uζ)) + end + return a_up +end + """ periodic_wrap(x, n) -> Int diff --git a/src/Vacuum/Vacuum.jl b/src/Vacuum/Vacuum.jl index 8a5b0f405d..86c8e8afac 100644 --- a/src/Vacuum/Vacuum.jl +++ b/src/Vacuum/Vacuum.jl @@ -241,11 +241,11 @@ class operator real and splitting a self-conjugate class into two half-size bloc # exp(-i(mθ − nζ)) on one field period exp_mn_basis = compute_fourier_coefficients(mtheta, m_modes, nzeta * nfp, n_modes; nfp=nfp) - # 23×23 patch, 20 radial / 40 angular polar nodes, 5-point Lagrange. Malhotra JCP 397 (2019) 108791 sec 3.2.2 - # grows these as N^(1/4); over our grids that spans M = 19..35 with no measured accuracy gain, so they are fixed. - PATCH_RAD = 11 - RAD_DIM = 20 - INTERP_ORDER = 5 + # Patch half-width mtheta/3 and max(12, mtheta/2) radial / twice as many angular polar nodes, 11-point Lagrange. The patch + # and quadrature grow with the grid; a fixed patch floors the error (Malhotra JCP 397 (2019) 108791 sec 3.2.2). + PATCH_RAD = mtheta ÷ 3 + RAD_DIM = max(12, round(Int, mtheta / 2)) + INTERP_ORDER = 11 # Shared within-period reflection when both surfaces allow it; `nothing` leaves the operator untransformed σ_map = stell_sym_map(plasma_surf, wall, nfp) @@ -290,16 +290,23 @@ class operator real and splitting a self-conjugate class into two half-size bloc grad_blocks = [zeros!(pool, T, nb * sz, nb * sz) for sz in sizes] green_blocks = [zeros!(pool, T, nb * sz, sz) for sz in sizes] - # Plasma–Plasma block - compute_3D_kernel_matrices!(grad_blocks, green_blocks, plasma_surf, plasma_surf, PATCH_RAD, RAD_DIM, INTERP_ORDER, phases, sym_basis) - - 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_basis) - # Wall–Wall block - compute_3D_kernel_matrices!(grad_blocks, green_blocks, wall, wall, PATCH_RAD, RAD_DIM, INTERP_ORDER, phases, sym_basis) - # Wall–Plasma block - compute_3D_kernel_matrices!(grad_blocks, green_blocks, wall, plasma_surf, PATCH_RAD, RAD_DIM, INTERP_ORDER, phases, sym_basis) + # The kernel's small per-row BLAS products run inside its threaded loop; a threaded BLAS would oversubscribe the cores + blas_nthreads = BLAS.get_num_threads() + BLAS.set_num_threads(1) + try + # Plasma–Plasma block + compute_3D_kernel_matrices!(grad_blocks, green_blocks, plasma_surf, plasma_surf, PATCH_RAD, RAD_DIM, INTERP_ORDER, phases, sym_basis) + + 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_basis) + # Wall–Wall block + compute_3D_kernel_matrices!(grad_blocks, green_blocks, wall, wall, PATCH_RAD, RAD_DIM, INTERP_ORDER, phases, sym_basis) + # Wall–Plasma block + compute_3D_kernel_matrices!(grad_blocks, green_blocks, wall, plasma_surf, PATCH_RAD, RAD_DIM, INTERP_ORDER, phases, sym_basis) + end + finally + BLAS.set_num_threads(blas_nthreads) end # Factor once per group (conjugate partner reuses it) diff --git a/test/runtests_vacuum.jl b/test/runtests_vacuum.jl index 2459349773..d2d9295473 100644 --- a/test/runtests_vacuum.jl +++ b/test/runtests_vacuum.jl @@ -958,12 +958,10 @@ using GeneralizedPerturbedEquilibrium.Vacuum: ConformalWall, EllipticalWall, Dee # or an offset built on the wrong side, which finiteness and Hermiticity both pass. @test real(wv[1, 1]) > real(compute_vacuum_response(inputs, nothing).wv[1, 1]) - # This configuration clears the one-cell floor, so it draws no sub-cell or near-fold warning - @test_logs min_level = Base.CoreLogging.Warn GeneralizedPerturbedEquilibrium.Vacuum.WallGeometry3D( - inputs, - GeneralizedPerturbedEquilibrium.Vacuum.PlasmaGeometry3D(inputs), - settings - ) + # This configuration clears the one-cell floor, so the wall draws no sub-cell or near-fold warning. + # The plasma is built outside the check: its 24 x 24 grid is anisotropic and warns about that. + plasma = GeneralizedPerturbedEquilibrium.Vacuum.PlasmaGeometry3D(inputs) + @test_logs min_level = Base.CoreLogging.Warn GeneralizedPerturbedEquilibrium.Vacuum.WallGeometry3D(inputs, plasma, settings) end # Field-periodic (layer-2) reduction: an nfp-periodic boundary makes the boundary-integral @@ -1191,17 +1189,46 @@ using GeneralizedPerturbedEquilibrium.Vacuum: ConformalWall, EllipticalWall, Dee end @testset "Kernel3D extract_patch!" begin - # 4×4 grid (npol=4, ntor=4), PATCH_DIM=3, center at (2,2); data is (16, 3) - data = zeros(16, 3) - data[:, 1] .= 1.0:16.0 - data[:, 2] .= 1.0 - data[:, 3] .= 0.0 - patch_out = zeros(3, 3, 3) + # 4×4 grid (npol=4, ntor=4), PATCH_DIM=3, center at (2,2) + data = collect(1.0:16.0) + patch_out = zeros(3, 3) GeneralizedPerturbedEquilibrium.Vacuum.extract_patch!(patch_out, data, 2, 2, 4, 4, 3) # Center of patch should be data at (2,2) = index 2+4*(2-1)=6 - @test isapprox(patch_out[2, 2, 1], 6.0) - @test all(isfinite, patch_out) - @test size(patch_out) == (3, 3, 3) + @test patch_out[2, 2] == 6.0 + @test patch_out[1, 1] == 1.0 + # Periodic wrap: centered at (1,1), corner (0,0) wraps to (4,4) = 16 + GeneralizedPerturbedEquilibrium.Vacuum.extract_patch!(patch_out, data, 1, 1, 4, 4, 3) + @test patch_out[1, 1] == 16.0 + end + + @testset "Kernel3D manufactured exterior solution converges at high order" begin + # χ from point charges inside a rotating ellipse is harmonic outside it, so K χ = G b with b = ∇χ · n on the + # surface's own (area-weighted) normals, which the kernel uses; solving for χ measures the kernel's error. + VAC = GeneralizedPerturbedEquilibrium.Vacuum + function surf(θ, φ) + R = 1 + 0.3cos(θ) + 0.06cos(θ - 3φ) + return [R * cos(φ), R * sin(φ), 0.35sin(θ) - 0.06sin(θ - 3φ)] + end + charges = [([cos(a), sin(a), z], q) for (a, z, q) in ((0.3, 0.05, 1.0), (2.1, -0.04, -0.7), (4.0, 0.0, 0.5))] + χ(x) = sum(q / norm(x - y) for (y, q) in charges) + ∇χ(x) = -sum(q * (x - y) / norm(x - y)^3 for (y, q) in charges) + function manufactured_error(mt) + nz = 3mt + pts = [surf(2π * (i - 1) / mt, 2π * (j - 1) / nz) for j in 1:nz for i in 1:mt] + inp = VacuumInput(; x=getindex.(pts, 1), y=getindex.(pts, 2), z=getindex.(pts, 3), mtheta_in=mt, nzeta_in=nz, + m_modes=[0], n_modes=[0], mtheta=mt, nzeta=nz, nfp=1) + plasma = VAC.PlasmaGeometry3D(VAC.expand_field_periods(inp)) + K, G = [zeros(mt * nz, mt * nz)], [zeros(mt * nz, mt * nz)] + # The patch, radial order and Lagrange order compute_vacuum_response uses + VAC.compute_3D_kernel_matrices!(K, G, plasma, plasma, mt ÷ 3, max(12, round(Int, mt / 2)), 11, [1.0]) + x = [χ(plasma.r[k, :]) for k in axes(plasma.r, 1)] + b = [dot(∇χ(plasma.r[k, :]), plasma.normal[k, :]) for k in axes(plasma.r, 1)] + return norm(K[1] \ (G[1] * b) - x) / norm(x) + end + e16, e32 = manufactured_error(16), manufactured_error(32) + @test e16 < 1e-3 + @test e32 < 5e-6 + @test e32 < 1e-2 * e16 # faster than h⁶ end end