3Spherical harmonic transforms on **WebGPU** (browser, fp32), modeled on
4[SHTNS](https://nschaeff.bitbucket.io/shtns/). This is a from-scratch
5TypeScript + WGSL implementation of the scalar transforms, structured after
6the SHTNS CUDA backend (`sht_gpu.cu` / `SHT/cuda_legendre.gen.cu`).
8**Live demo:** https://concept-collection.github.io/shtns-webgpu/
9(validation suite: [test.html](https://concept-collection.github.io/shtns-webgpu/test.html))
11## Scope (v0.1)
13- **Scalar transforms of real fields**, both directions:
14 - `synth()` — spectral → spatial (SHTNS `SH_to_spat`)
15 - `analys()` — spatial → spectral (SHTNS `spat_to_SH`)
16- **Gauss–Legendre grid** in latitude, uniform longitude grid.
17- **fp32 on the GPU** end to end; all precomputation (Gauss nodes/weights,
18 recurrence coefficients, twiddle factors) is done host-side in f64.
19- Not (yet) implemented: vector (spheroidal/toroidal) transforms, complex
20 fields, regular grids, `mres > 1`, Schmidt/4π normalizations, on-the-fly
21 truncation (`llim < lmax`).
23## Conventions (= SHTNS defaults)
25- **Orthonormal** spherical harmonics **with Condon–Shortley phase**.
26- Spectral coefficients `Q_lm` are complex, stored for `m >= 0` with the
27 SHTNS LM ordering: `for m = 0..mmax: for l = m..lmax`, interleaved
28 `[re, im]` (`Float32Array` of length `2*nlm`). Real fields imply
29 `Q_{l,-m} = (-1)^m conj(Q_lm)`; `m = 0` coefficients must be real.
30- Spatial fields are `Float32Array[nlat * nphi]`, **phi-contiguous**,
31 latitudes ordered north → south (colatitude increasing).
33## Usage
35```ts
36import { ShtPlan, requestShtDevice, lmIndex } from 'shtns-webgpu';
38const device = await requestShtDevice(); // or your own GPUDevice
39const plan = await ShtPlan.create(device, { lmax: 127, mmax: 127, nlat: 128, nphi: 256 });
41const qlm = new Float32Array(2 * plan.nlm);
42qlm[2 * lmIndex(127, 8, 5)] = 1.0; // Y_8^5
44const spat = await plan.synth(qlm); // nlat*nphi field
45const qBack = await plan.analys(spat); // back to spectral
46plan.destroy();
47```
49For GPU-resident pipelines (no readback), use `plan.encodeSynth(encoder)` /
50`plan.encodeAnalys(encoder)` with the exposed `qlmIn` / `qlmOut` / `spatBuf`
51buffers.
53Constraints checked at plan creation: `nlat > lmax` (Gauss quadrature
54exactness), `nphi >= 2*mmax + 1` (no aliasing).
56## How it works
58Same two-stage split as SHTNS:
601. **Legendre stage** (`src/wgsl/leg.ts`): associated Legendre functions are
61 generated *on the fly* inside the shader by the standard 3-term
62 recurrence over `l` (coefficients from `legendre_precomp()`-equivalent
63 host code, `src/coeffs.ts`). Underflow of `sin^m(theta)` — fatal in fp32
64 beyond `m ≈ 75` — is handled with the SHTNS extended-range scheme
65 (`SHT_SCALE_FACTOR = 2^56`, `SHT_ACCURACY = 1e-15`, per-thread integer
66 exponent), ported from the `HI_LLIM` path of `SHT/cuda_legendre.gen.cu`.
67 Synthesis runs one thread per latitude and one workgroup row per `m`;
68 analysis runs one workgroup per `m` with a shared-memory tree reduction
69 over latitudes (the portable equivalent of SHTNS's warp shuffles).
702. **Fourier stage** (`src/wgsl/fourier.ts`): batched radix-2 Stockham FFT
71 in workgroup memory (one workgroup per latitude row) when `nphi` is a
72 power of two that fits (`16*nphi <= maxComputeWorkgroupStorageSize`);
73 otherwise a direct band-limited DFT. Twiddles come from a host-computed
74 f64 table — device `sin`/`cos` is only guaranteed to ~2^-11 under
75 Vulkan, which would otherwise dominate the error budget.
77All problem sizes are baked into the WGSL at plan creation (the WGSL
78equivalent of SHTNS's NVRTC runtime compilation).
80## Accuracy (fp32)
82Relative L2 errors vs the double-precision reference (`src/reference.ts`),
83random spectra, measured on SwiftShader (results on hardware GPUs are the
84same to within noise since the arithmetic is IEEE fp32):
86| lmax | synthesis | analysis | round trip |
87|-----:|----------:|---------:|-----------:|
88| 15 | 6e-7 | 3e-7 | 6e-7 |
89| 63 | 5e-6 | 2e-6 | 3e-6 |
90| 127 | 7e-6 | 3e-6 | 5e-6 |
91| 255 | 2e-5 | 6e-6 | 1e-5 |
92| 399 | 9e-5 | 1e-5 | 2e-5 |
94SHTNS itself switches its fp32 GPU recurrence to f64 above `lmax = 128`
95(`SHT_L_RESCALE_FLY_FLOAT`); WGSL has no f64, so past that point accuracy
96degrades gracefully as above. Fine for visualization; for scientific use
97keep `lmax ≲ 128` or wait for the float-float recurrence (planned).
99## Develop / test
101```sh
102npm install
103npm run dev # demo at http://localhost:5173
104npm run test:node # f64 math tests (no GPU needed)
105npm run test:gpu # builds, then runs the browser suite in headless Chrome
106 # (falls back to SwiftShader software WebGPU; CHROME_PATH to override)
107```
109## Roadmap
111- Vector transforms (spheroidal/toroidal), gradients — port of `leg_m_kernel<1>`.
112- Float-float (double-single) recurrence option for full accuracy at high lmax.
113- Latitude parity folding (2× Legendre work reduction, as in SHTNS).
114- Subgroup (warp) reductions where available, replacing the shared-memory tree.
115- `mres > 1`, truncated transforms, regular grids.