From a6fb05e06445fafb5ec3a38f65cfcf733287f7f8 Mon Sep 17 00:00:00 2001
From: Jake Halpern <jhalpern@purdue.edu>
Date: Wed, 11 Nov 2026 13:28:19 -0500
Subject: [PATCH 1/6] Vacuum - FEATURE! - Grid-scaled singular-correction patch
 and spectral surface tangents for the 3D kernel

The 3D singular correction keeps its algorithm and code; only its parameters change. The patch half-width is
mtheta/3 grid points (was 11), the radial Gauss order max(12, mtheta/2) (was 20, angular twice that) and the Lagrange
order 11 (was 5). A fixed patch shrinks in physical size as the grid is refined and floors the error; growing it with
the grid keeps the high-order convergence (W7-X |m|,|n| <= 8: 8e-4 at 26x34 to 6e-9 at 100x132 per period).

PlasmaGeometry3D and WallGeometry3D take their tangents and normals from the trigonometric interpolant of the grid
samples (spectral differentiation matrix) in place of bicubic splines. Spline tangents converge as h^4 and cap the
kernel there (W7-X 1.9e-6 at 100x132).
---
 src/Vacuum/DataTypes.jl | 29 +++++++----------------------
 src/Vacuum/Utilities.jl | 33 +++++++++++++++++++++++++++++++++
 src/Vacuum/Vacuum.jl    | 10 +++++-----
 3 files changed, 45 insertions(+), 27 deletions(-)

diff --git a/src/Vacuum/DataTypes.jl b/src/Vacuum/DataTypes.jl
index a15f08e330..6b0ca5497a 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
 
@@ -731,7 +723,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 +789,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
 
diff --git a/src/Vacuum/Utilities.jl b/src/Vacuum/Utilities.jl
index 906b9819cf..2cf8ddeb1c 100644
--- a/src/Vacuum/Utilities.jl
+++ b/src/Vacuum/Utilities.jl
@@ -82,6 +82,39 @@ 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 (Trefethen, Spectral Methods in
+MATLAB, ch. 3). 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
+        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
+
 """
     periodic_wrap(x, n) -> Int
 
diff --git a/src/Vacuum/Vacuum.jl b/src/Vacuum/Vacuum.jl
index 8a5b0f405d..bea8dc3f5b 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)

From a7d0c64e891279138ea8454518a6e3e1115a140b Mon Sep 17 00:00:00 2001
From: Jake Halpern <jhalpern@purdue.edu>
Date: Wed, 11 Nov 2026 13:28:30 -0500
Subject: [PATCH 2/6] Vacuum - TEST - Manufactured exterior solution for the 3D
 kernel

Point charges inside a rotating ellipse give a harmonic exterior field; solving K chi = G b with the kernel's own normals measures its error (3.9e-4 at mtheta 16, 2.8e-6 at 32).
---
 test/runtests_vacuum.jl | 27 +++++++++++++++++++++++++++
 1 file changed, 27 insertions(+)

diff --git a/test/runtests_vacuum.jl b/test/runtests_vacuum.jl
index 2459349773..4896890847 100644
--- a/test/runtests_vacuum.jl
+++ b/test/runtests_vacuum.jl
@@ -1203,6 +1203,33 @@ using GeneralizedPerturbedEquilibrium.Vacuum: ConformalWall, EllipticalWall, Dee
             @test all(isfinite, patch_out)
             @test size(patch_out) == (3, 3, 3)
         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
+            surf(θ, φ) = (R=1 + 0.3cos(θ) + 0.06cos(θ - 3φ); [R * cos(φ), R * sin(φ), 0.35sin(θ) - 0.06sin(θ - 3φ)])
+            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
 
     @testset "reset_caches! empties the run-filled caches" begin

From 529990e3d6ef893ed67f2d3626b59d9a29861c1e Mon Sep 17 00:00:00 2001
From: Jake Halpern <jhalpern@purdue.edu>
Date: Thu, 12 Nov 2026 11:30:06 -0500
Subject: [PATCH 3/6] =?UTF-8?q?Vacuum=20-=20IMPROVEMENT=20-=202=C3=97=20tr?=
 =?UTF-8?q?igonometric=20upsampling=20of=20the=20singular-correction=20pat?=
 =?UTF-8?q?ch=20in=20the=203D=20kernel?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

- Interpolate the source geometry onto the polar nodes from a 2× trig-upsampled patch
  (`trig_upsample_matrix`, `trig_upsample`), and fold the scattered weights back to the
  collocation nodes through the restricted upsampling matrices, with field periods folded into Uζ
- Carry (r, ∂r/∂θ, ∂r/∂ζ) as SVector{9} per node; normals are formed from the interpolated
  tangents, replacing `compute_polar_normal!`
- Thread-local `KernelWorkspace` holds all per-observer scratch
- Bind the clamped patch radius to a new name so it is not boxed in the threaded closure
- Map trig-matrix offsets to the symmetric range so sin(πr/n) never rounds near π;
  both matrices now agree with a BigFloat reference to machine precision, odd and even n
- Pin BLAS to one thread around the kernel calls (restored in `finally`): its small
  products run inside the threaded observer loop
- `extract_patch!` is element-generic; test covers periodic wrap

W7-X m16: errors unchanged by the SVector, BLAS and trig changes (2.427e-6 at 56×74,
1.495e-8 at 82×108); 82×108 runs in 6.4 s vs 8.9 s before this change (-t 4).
---
 src/Vacuum/Kernel3D.jl  | 253 ++++++++++++++++++++++------------------
 src/Vacuum/Utilities.jl |  44 ++++++-
 src/Vacuum/Vacuum.jl    |  27 +++--
 test/runtests_vacuum.jl |  17 ++-
 4 files changed, 207 insertions(+), 134 deletions(-)

diff --git a/src/Vacuum/Kernel3D.jl b/src/Vacuum/Kernel3D.jl
index 345671fe3f..c0603dd90f 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
-
-"""
-    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.
+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.
 """
-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.
@@ -334,25 +366,31 @@ 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
+
+    # 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)))
 
-    # 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]
+    # 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)]
+
+    # 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)
 
@@ -375,38 +413,28 @@ function compute_3D_kernel_matrices!(
             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)
@@ -414,22 +442,21 @@ function compute_3D_kernel_matrices!(
             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 = -χ
+            # Far field at this patch node (weighted by 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])
 
-            # 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 2cf8ddeb1c..980116d7a4 100644
--- a/src/Vacuum/Utilities.jl
+++ b/src/Vacuum/Utilities.jl
@@ -85,14 +85,16 @@ end
 """
     trig_derivative_matrix(n) -> Matrix{Float64}
 
-Derivative at the nodes of the trigonometric interpolant on `n` equispaced periodic nodes (Trefethen, Spectral Methods in
-MATLAB, ch. 3). For even `n` the Nyquist mode is a cosine, whose derivative vanishes at the nodes.
+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
@@ -115,6 +117,44 @@ function trig_tangents(r::Matrix{Float64}, mtheta::Int, nzeta::Int)
     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 bea8dc3f5b..86c8e8afac 100644
--- a/src/Vacuum/Vacuum.jl
+++ b/src/Vacuum/Vacuum.jl
@@ -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 4896890847..6969adee53 100644
--- a/test/runtests_vacuum.jl
+++ b/test/runtests_vacuum.jl
@@ -1191,17 +1191,16 @@ 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 e0d1ee998636841c651be3e8dd0ce13e7d2b4bf4 Mon Sep 17 00:00:00 2001
From: Jake Halpern <jhalpern@purdue.edu>
Date: Fri, 13 Nov 2026 08:13:34 -0500
Subject: [PATCH 4/6] VACUUM - MINOR - SVector optimizations and
 guards/warnings

---
 src/Vacuum/DataTypes.jl | 12 ++++++++++--
 src/Vacuum/Kernel3D.jl  | 21 ++++++++++++---------
 2 files changed, 22 insertions(+), 11 deletions(-)

diff --git a/src/Vacuum/DataTypes.jl b/src/Vacuum/DataTypes.jl
index 6b0ca5497a..c6cbc18c65 100644
--- a/src/Vacuum/DataTypes.jl
+++ b/src/Vacuum/DataTypes.jl
@@ -477,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,
@@ -820,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 c0603dd90f..15958da35d 100644
--- a/src/Vacuum/Kernel3D.jl
+++ b/src/Vacuum/Kernel3D.jl
@@ -342,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
@@ -382,6 +383,10 @@ function compute_3D_kernel_matrices!(
     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
@@ -398,17 +403,16 @@ 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
@@ -440,12 +444,11 @@ 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)
+            d, idx_col = period_col[idx_src]
 
             # Far field at this patch node (weighted by 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])
+            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])
 
             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ζ)

From 3e4ba53c5f93cf3a0385a8fe8fde74148cf62381 Mon Sep 17 00:00:00 2001
From: Jake Halpern <jhalpern@purdue.edu>
Date: Fri, 13 Nov 2026 08:35:42 -0500
Subject: [PATCH 5/6] Vacuum - TEST - Build the plasma outside the no-warning
 wall check

The 24 x 24 periodic test grid has toroidal spacing 5.7x the poloidal one, so PlasmaGeometry3D
now warns about it inside @test_logs. The check is about the wall's sub-cell and near-fold
warnings only. Vacuum unit tests 431/431 (Feynman).
---
 test/runtests_vacuum.jl | 10 ++++------
 1 file changed, 4 insertions(+), 6 deletions(-)

diff --git a/test/runtests_vacuum.jl b/test/runtests_vacuum.jl
index 6969adee53..aa2159a375 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

From 59483ecfd6e8b0496be2ca9a745f9b01f908d0f4 Mon Sep 17 00:00:00 2001
From: Jake Halpern <jhalpern@purdue.edu>
Date: Fri, 13 Nov 2026 12:19:32 -0500
Subject: [PATCH 6/6] Vacuum - MINOR - Format the manufactured-solution test

---
 test/runtests_vacuum.jl | 5 ++++-
 1 file changed, 4 insertions(+), 1 deletion(-)

diff --git a/test/runtests_vacuum.jl b/test/runtests_vacuum.jl
index aa2159a375..d2d9295473 100644
--- a/test/runtests_vacuum.jl
+++ b/test/runtests_vacuum.jl
@@ -1205,7 +1205,10 @@ using GeneralizedPerturbedEquilibrium.Vacuum: ConformalWall, EllipticalWall, Dee
             # χ 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
-            surf(θ, φ) = (R=1 + 0.3cos(θ) + 0.06cos(θ - 3φ); [R * cos(φ), R * sin(φ), 0.35sin(θ) - 0.06sin(θ - 3φ)])
+            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)