Repository navigation
Vacuum - FEATURE! - High-order singular quadrature for the 3D vacuum kernel - #673
jhalpern30 wants to merge 6 commits into
Conversation
…l 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).
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).
…rrection patch in the 3D kernel
- 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).
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).
github-actions
Bot
added
changed-results
|
This pull request is missing a reviewer. If you are not ready to name them, mark this pull request as a draft. |
|
@logan-nc going to wait on merging this in until I have finalized the reviewer responses, but this is the code update motivated by the reviewer comments + APS discussions I was talking about. The TL/DR is the before and after for the paper convergence figures - we actually achieve spectral convergence now instead of ~h^4 and error absolutely TANKS. Can get down to almost machine precision on some cases on a laptop. Basically required me going error-hunting, since even though I had a theoretically spectrally convergent scheme there were specific implementation details (both in the singular quadrature and also relying on bicubic derivatives) limiting my convergence. The funny thing is for Fig. 2 below, the 3D code actually beats the 2D runtime for some low error rates because it converges so much faster 
|
Release note
wvmoves; 2D path unchanged (harness @ e0d1ee9)wv) used for results and re-check its grid convergence, since 3D values move; 3D grids now needmtheta≥ 9 and at least 7 toroidal points over the full torus (smaller grids raise an error), and the new grid-aspect and wall-resolution warnings flag grids that cannot reach full accuracyThe 3D vacuum response now converges spectrally with
mthetaandnzeta. On W7-X it self-converges to ~1e-11, where it used to stall at 1e-5 to 1e-4. A 3D run of an axisymmetric case now reaches the 2D code's answer to its floor (4e-7 without a wall, 6e-8 with one). New warnings flag grids that cannot reach that accuracy.Summary
Review package with figures: https://claude.ai/artifact/4HaFFsohHYdSqRvuKTKyK6
The 3D kernel's singular correction had a floor that refining the grid did not remove. Three changes remove it, and all three are needed (the ablation figure in the package removes them one at a time):
PlasmaGeometry3DandWallGeometry3Dtake ∂r/∂θ and ∂r/∂ζ from the trigonometric interpolant (trig_tangents) instead of periodic bicubic splines, which were O(h⁴) and capped the whole scheme.PATCH_RAD = mtheta ÷ 3,RAD_DIM = max(12, mtheta/2)(ANG_DIM = 2·RAD_DIM),INTERP_ORDER = 11. The patch is now a fixed physical region, so the smooth-remainder error decays spectrally. Malhotra's N^¼ growth is a cost choice for FMM codes; the dense build here has no FMM to protect.trig_upsample_matrix,trig_upsample). The 11-point Lagrange stencil runs on the fine patch, and the polar weights are folded back to the collocation nodes by the transposed upsampling, with field periods folded into Uζ (scatter_polar_to_row!). The far field stays the trapezoidal rule on the collocation grid.Guards (new warnings and assert)
PlasmaGeometry3D: warns when the RMS spacing ratio hζ/hθ > 1.75 on a 3D boundary and suggests annzeta. This replaces the old warning at > 10.WallGeometry3D(conformal wall): warns when the wall and plasma normals differ by more than 2°. This is about 1e-4 in Wᵛ: a sharp convex plasma edge fans out over the gap faster than the grid can follow.compute_3D_kernel_matrices!: asserts that observer and source share one grid.Performance
SVector{9}, soP2Gis swept once. Per-thread scratch now lives in a parametricKernelWorkspace{T}.finally, because the per-row products run inside the threaded observer loop.Tests
Kernel3D manufactured exterior solution converges at high order. Point charges inside a rotating ellipse give a harmonic exterior field, and solving K χ = G b on the kernel's own normals measures its error. The test asserts e16 < 1e-3, e32 < 5e-6 and a rate faster than h⁶.extract_patch!test rewritten for the element-generic signature, including the periodic wrap.Regression report
Baseline is the branch point a1149c5. The harness has no 3D vacuum case, so this checks only that 2D results are unchanged. Vacuum unit tests:
runtests_vacuum.jl431/431 at 3e4ba53.Notes for reviewers
scatter_polar_to_row!. The folded patch is built once about node 1 and shifted to each observer. Columns pushed into the next field period pick up one more phase step, which assumes geometric phases (phases[d] = wrap_phase^(d-1), ascompute_vacuum_responsebuilds them).