1# Reducing spherical-harmonic transforms in $\Delta_\Gamma$
3**Summary.** Algorithm 4 costs 12 transforms per matvec (8 syntheses, 4 analyses). A flux-form
4reformulation, with weights chosen so that every analyzed field is smooth on $S^2$, evaluates the
5same operator in **6 transforms** (4 syntheses, 2 analyses), or **7** with the divergence split
6against the round sphere that §5 turned out to require. Notation follows `algos.pdf`.
8---
10## 1. Where the current cost comes from
12| Algorithm 4 line | Work | Transforms |
13|---|---|---|
14| 1 | $\partial_\theta u,\ \partial_\varphi u$ | 2 $\mathcal{S}$ |
15| 5 | analysis of 3 Cartesian components of $\nabla_\Gamma u$ | 3 $\mathcal{A}$ |
16| 6 | $\partial_\theta$ and $\partial_\varphi$ of each of those 3 components | 6 $\mathcal{S}$ |
17| 8 | final analysis | 1 $\mathcal{A}$ |
18| | | **12** |
20Two sources of waste:
221. The gradient is carried in **ambient $\mathbb{R}^3$ components** — 3 fields for an intrinsically
23 2-dimensional object.
242. Line 6 takes **both** derivatives of **each** component, where the divergence needs only one
25 derivative of each of two fluxes.
27There is also a possible free win independent of everything below: Algorithm 1 as specified returns
28all five derivatives. The Laplacian path needs only $\partial_\theta u$ and $\partial_\varphi u$;
29the second-derivative and mixed-derivative machinery is used exclusively by Algorithm 3 (curvature).
30If `surface_screened_laplacian` calls Algorithm 1 wholesale, it is doing 5 syntheses where 2 suffice
31at lines 1 and 6.
33---
35## 2. The constraint that shapes the solution
37The Cartesian design in Algorithm 4 exists to avoid pole singularities, and it is correct to do so.
38The relevant property is smoothness **as a scalar function on $S^2$**, since that is what controls
39SH coefficient decay and hence whether $\mathcal{A}$ is meaningful.
41| Quantity | Smooth on $S^2$? |
42|---|---|
43| $\partial_\varphi u$ | yes — exactly band-limited, eq. (2.2) |
44| $\sin\theta\,\partial_\theta u$ | yes — exactly band-limited, eq. (2.4) |
45| $\partial_\theta u$ | **no** — bounded, but $\varphi$-dependent limit at the poles |
46| $g_{\theta\theta},\ g^{\theta\theta},\ V_\theta,\ V_\varphi$ | **no** in general |
47| $(\nabla_\Gamma u)_x,\ (\nabla_\Gamma u)_y,\ (\nabla_\Gamma u)_z$ | yes |
49Concretely, for the ellipsoid $X = (a\sin\theta\cos\varphi,\ b\sin\theta\sin\varphi,\ c\cos\theta)$,
50$|X_\theta|^2 \to a^2\cos^2\varphi + b^2\sin^2\varphi$ as $\theta\to 0$: no limit exists.
52Algorithm 4 never analyzes anything in the "no" rows — the non-smooth quantities appear only as
53pointwise grid factors. **Any replacement must preserve this property.** The naive flux form
54$P = \sqrt{g}\,(g^{\theta\theta}u_\theta + g^{\theta\varphi}u_\varphi)$,
55$Q = \sqrt{g}\,(g^{\varphi\theta}u_\theta + g^{\varphi\varphi}u_\varphi)$ does not:
56on the round sphere with $u = x$, $Q = -\sin\varphi$, which is not a function on $S^2$.
58### The correct weighting
60$P$ and $Q$ are $\sqrt{g}$ times the contravariant components of $G := \nabla_\Gamma u$. Using
61$\det[X_\theta, X_\varphi, n] = -\sqrt{g}$:
63$$P = -\,G\cdot(X_\varphi \times n), \qquad \sin\theta\,Q = -\,G\cdot\big(n \times \sin\theta\,X_\theta\big).$$
65Every factor on the right is smooth on $S^2$: $G$ smooth, $n$ smooth, $X_\varphi$ smooth, and
66$\sin\theta\,X_\theta$ smooth because it is exactly band-limited by the recurrence already
67implemented. Hence $P$ and $\tilde{Q} := \sin\theta\,Q$ are analyzable, on the same footing and for
68the same structural reason as the Cartesian gradient components.
70The cross products are the smoothness certificate only — they are not needed in the code.
72---
74## 3. Precompute (once per surface, grid space)
76Replaces `_precompute_metric_quantities()`. From the embedding coefficients $\hat{X}^m_\ell$, obtain
77$X_\varphi$ and $\sin\theta\,X_\theta$ componentwise via Algorithm 1, where the latter is the
78**undivided** output of Algorithm 1 line 4, i.e. $\mathcal{S}(v^m_\ell)$ with no $/\sin\theta$.
79Then, pointwise:
81$$\tilde{g}_{\theta\theta} := |\sin\theta\,X_\theta|^2, \qquad
82 \tilde{g}_{\theta\varphi} := (\sin\theta\,X_\theta)\cdot X_\varphi, \qquad
83 g_{\varphi\varphi} := |X_\varphi|^2$$
85$$J := \frac{\sqrt{\tilde{g}_{\theta\theta}\,g_{\varphi\varphi} - \tilde{g}_{\theta\varphi}^{\,2}}}{\sin^2\theta}
86 \qquad\text{so that } \sqrt{\det g} = J\sin\theta$$
88(The radicand is $\sin^2\theta\det g = J^2\sin^4\theta$, so the square root is $J\sin^2\theta$ — hence
89$\sin^2\theta$, not $\sin\theta$, in the denominator.)
91Store four scalar grid arrays:
93$$p_1 = \frac{g_{\varphi\varphi}}{J\sin^2\theta}, \qquad
94 p_2 = -\frac{\tilde{g}_{\theta\varphi}}{J\sin^2\theta}, \qquad
95 q_2 = \frac{\tilde{g}_{\theta\theta}}{J\sin^2\theta}, \qquad
96 r = \frac{1}{J\sin^2\theta}$$
98All four are bounded: the $\sin^2\theta$ denominators cancel against vanishing numerators
99($\tilde{g}_{\theta\varphi} = O(\sin^2\theta)$, $g_{\varphi\varphi} = O(\sin^2\theta)$), the same
100finite limits the current $V_\theta, V_\varphi$ have. Note $p_1, p_2, q_2$ are bounded where
101$V_\varphi = O(1/\sin\theta)$ is not.
103**Three scalar arrays replace the six components of $V_\theta, V_\varphi$.** The surface
104representation is unchanged: $X_\theta, X_\varphi$ still come componentwise from $\hat{X}^m_\ell$
105via Algorithm 1.
107---
109## 4. The per-matvec algorithm
111Input $\{u^m_\ell\}$; output $\{(\Delta_\Gamma u)^m_\ell\}$.
113| # | Step | Transforms |
114|---|---|---|
115| 1 | $v^m_\ell \leftarrow \alpha^+(\ell-1,m)u^m_{\ell-1} + \alpha^-(\ell+1,m)u^m_{\ell+1}$; $A \leftarrow \mathcal{S}(v^m_\ell)$ | $\mathcal{S}$ |
116| 2 | $B \leftarrow \mathcal{S}(im\,u^m_\ell)$ | $\mathcal{S}$ |
117| 3 | $P \leftarrow p_1 A + p_2 B$, $\tilde{Q} \leftarrow p_2 A + q_2 B$ — pointwise | — |
118| 4 | $\hat{P} \leftarrow \mathcal{A}(P)$, $\hat{\tilde{Q}} \leftarrow \mathcal{A}(\tilde{Q})$ | 2 $\mathcal{A}$ |
119| 5 | $s^m_\ell \leftarrow \alpha^+(\ell-1,m)\hat{P}^m_{\ell-1} + \alpha^-(\ell+1,m)\hat{P}^m_{\ell+1} + im\,\hat{\tilde{Q}}^m_\ell$ | — |
120| 6 | $\Delta_\Gamma u \leftarrow r \cdot \mathcal{S}(s^m_\ell)$ | $\mathcal{S}$ |
121| 7 | $\{(\Delta_\Gamma u)^m_\ell\} \leftarrow \mathcal{A}(\Delta_\Gamma u)$; zero $\ell \ge L-2$ | $\mathcal{A}$ |
123**4 syntheses + 2 analyses = 6**, versus 12.
125- $A$ and $B$ are exactly $\sin\theta\,\partial_\theta u$ and $\partial_\varphi u$.
126- Steps 1 and 5 use the **same** precomputed $\alpha^\pm$ table; step 5 is the adjoint-style reuse
127 of the shift already implemented for step 1. Adding the two flux contributions in coefficient
128 space before synthesizing is what saves the final pair of transforms.
129- The **only** division by $\sin\theta$ anywhere is folded into $p_1, p_2, q_2, r$ at precompute
130 time. The per-matvec path contains none.
132### Implemented variation (2026-08-05): the $\varphi$-flux never needs the Legendre basis
134Step 4's analysis of $\tilde{Q}$ exists only so step 5 can apply $\partial_\varphi$ — but
135$\partial_\varphi$ is diagonal in the Fourier index, so the implementation differentiates
136$\tilde{Q}$ on the grid instead: FFT each latitude row, multiply mode $m$ by $im$ (zeroing
137$m \ge L-2$ to mirror the top-degree filter; the Fourier analysis stage truncates $m > m_{\max}$
138for free), inverse FFT. **5 Legendre transforms + one Legendre-free FFT derivative**, versus 6.
139The caveat is that the grid route skips $\tilde Q$'s band projection in $\ell$; measured, this
140does not bite — the band-edge spectra are identical to the 6-transform route's (the $m$ mask and
141the final analysis's projection contain it), the Algorithm-4 A/B agreement is unchanged
142($3.6\times10^{-4}$ after 20 steps at $L=63$), and the step gets ~8% faster at $L=255$
143(~2% at $L=127$, where transform batching had already amortized most of what this removes).
145---
147## 5. Numerical trade-off
149Both schemes contain two powers of $\sin\theta$ division in total. What differs is **placement**.
151- **Algorithm 4** spends them in separate stages, one before line 5 and one after. The intervening
152 analysis suppresses the polar spike: Gauss–Legendre weights give $w_1 = O(L^{-2})$ at the polar
153 ring, so a grid error of $\varepsilon L$ there contributes
154 $\sim L^{-2}\cdot L^{1/2}\cdot \varepsilon L = \varepsilon L^{-1/2}$ to any coefficient. Stage 2
155 then starts from clean coefficients and incurs a *fresh* $\varepsilon L$. The two amplifications
156 never multiply. Net grid-space relative error: $\varepsilon L$.
157- **The new scheme** has no division at all through step 5, then pays for both powers at once in
158 $r = O(L^2)$ at step 6 — one event, with no intervening analysis to break it in half. Net
159 grid-space relative error: $\varepsilon L^2$.
161The mechanism: with $N \approx L+1$ Gauss–Legendre nodes, $1 - x_1 = O(N^{-2})$ so
162$\sin\theta_1 = O(N^{-1})$. Since $s = J\sin^2\theta\,\Delta_\Gamma u$ is $O(L^{-2})$ at the polar
163ring but $O(1)$ over the bulk, and synthesis commits roundoff scaled by the field's *global* size at
164every node alike, multiplying by $r \sim L^2$ recovers the signal and inflates the noise.
166| | divisions | placement | grid-space relative error |
167|---|---|---|---|
168| Algorithm 4 | $\sin\theta$, $\sin\theta$ | separated by $\mathcal{A}$ | $\varepsilon L$ |
169| Six-transform | $\sin^2\theta$ | all at the end | $\varepsilon L^2$ |
171**It does reach the returned coefficients, and it matters.** The suppression argument above is
172right as far as it goes — step 7's analysis knocks the spike down to $\varepsilon L^{1/2}$, and this
173document originally concluded from that the extra power would be invisible inside the solve. It is
174not, and the reason is not about accuracy. Measured on the default ellipsoid at $L=63$: starting
175from the exact uniform steady state, three Richardson iterations per step leave a standing
176coefficient-space perturbation $50\times$ the no-correction floor ($1.4\times10^{-5}$ vs
177$2.8\times10^{-7}$), against $\sim\!1\times$ for Algorithm 4. That perturbation is static, polar,
178and re-injected every step. In a Turing problem the pattern is seeded by whatever is largest in the
179unstable band, so a forcing four orders below the field selects the nucleation site: the run grows a
180spot at the pole, on every seed, regardless of the initial condition.
182**The fix is to keep $r$ off the round sphere.** Write $p_1 = 1 + \delta p_1$, $q_2 = 1 + \delta q_2$
183($p_2$ is already zero on the sphere). The sphere's share of the divergence is the cancelling part,
184and it is known in closed form: $\sin\theta\,\partial_\theta A + \partial_\varphi B =
185-\sin^2\theta\,\Delta_{S^2}u$, and $\Delta_{S^2}$ is diagonal. So
187$$\Delta_\Gamma u = -\frac{1}{J}\,\Delta_{S^2}u \;+\; r\,(\sin\theta\,\partial_\theta P'
188 + \partial_\varphi \tilde{Q}'), \qquad P' = \delta p_1 A + p_2 B, \quad
189 \tilde{Q}' = p_2 A + \delta q_2 B$$
191with $1/J = r\sin^2\theta$ bounded. Only the geometry *deviation* now meets the concentrated
192division. Cost: one extra synthesis per species per iteration for $-\lambda u$ — 7 transforms, not
1936 — which batches into the gradient's existing grouped call and measures at ~10% of a step, against
194$3\times$ for reverting to Algorithm 4. $\delta p_1$ and $\delta q_2$ must be formed in float64 at
195precompute time (`src/geom/geometry.ts`): on a near-sphere they *are* the small quantity, and
196subtracting 1 in float32 on device would lose them.
198Measured against Algorithm 4 through a real run (relative $L^2$ of $u$ at $t=8$, $\texttt{niter}=6$):
200| | plain flux | sphere-split |
201|---|---|---|
202| ellipsoid, $L=63$ | $3.5\times10^{-4}$ | $8.9\times10^{-6}$ |
203| blob, $L=63$ | $4.4\times10^{-4}$ | $8.5\times10^{-6}$ |
204| ellipsoid, $L=127$ | $7.2\times10^{-3}$ | $6.0\times10^{-6}$ |
206and the polar noise gain is asserted in `test/fluxChecks.ts`, which fails at $50\times$ on the
207unsplit form.
209The residual $\varepsilon L^2$ still applies to grid values of $\Delta_\Gamma u$ consumed directly
210— a nonlinear reaction term, max-norm diagnostics, an adaptive error estimator — for the deviation
211part alone.
213Algorithm 1 line 7 already divides by $\sin^2\theta$, so the code is exposed to $\varepsilon L^2$
214today — just on the second-derivative path, which the Laplacian never touches.
216---
218## 5a. float32 / WebGPU
220Target is WebGPU, which is float32-only: $\varepsilon = 2^{-24} \approx 6\times10^{-8}$ (spacing
221$2^{-23} \approx 1.2\times10^{-7}$). No float64 fallback exists on device. All estimates in §5 are
222linear in $\varepsilon$, so they scale directly:
224| | $\varepsilon\sqrt{L}$ (coeffs) | $\varepsilon L$ (Alg. 4 grid) | $\varepsilon L^2$ (new, grid) |
225|---|---|---|---|
226| $L=64$ | $5\times10^{-7}$ | $4\times10^{-6}$ | $2\times10^{-4}$ |
227| $L=128$ | $7\times10^{-7}$ | $8\times10^{-6}$ | $1\times10^{-3}$ |
228| $L=256$ | $1\times10^{-6}$ | $1.5\times10^{-5}$ | $4\times10^{-3}$ |
230**Coefficient space is fine** (~$10^{-6}$), which is the floor a float32 iterative solve sits at
231anyway. **Grid space is not**: 0.1–0.4% relative on the polar rings at $L\ge128$. For a
232reaction–diffusion solver this matters only if grid-space $\Delta_\Gamma u$ is consumed outside the
233matvec. If the IMEX splitting evaluates $f(u)$ from $u$ on the grid (typical), it never is.
235**Two float32-specific arguments in favour of the new scheme:**
237- Baseline SHT roundoff accumulates per transform ($\sim\varepsilon\sqrt{L}$ to $\varepsilon L$
238 each). Running 6 transforms instead of 12 halves that accumulation. On the coefficient-space error
239 GMRES actually sees, this plausibly outweighs the polar term — the new scheme may be *more*
240 accurate end-to-end in float32. Not asserted without measurement.
241- 3 weight arrays instead of 6 halves per-matvec texture/buffer traffic. On GPU that is often the
242 real bottleneck, independent of arithmetic.
244**Mitigations available without float64:**
246- **Pairwise or blocked summation in the Legendre sum over $\ell$.** The single highest-value
247 float32 change, and it benefits the existing code too. See §5b.
248- **Double-float (`f32x2`) arithmetic** for the pointwise steps 3 and 6 if needed — cheap, no
249 transforms involved. Does not help with transform roundoff, which is the dominant term, so try
250 summation order first.
251- **CPU precompute in float64.** JS `Number` is float64, so §3 can run on the CPU regardless of
252 WebGPU's limits, with float32 weights uploaded. Cost is CPU-side SHTs plus upload, paid once per
253 surface update; viable if the surface evolves slowly or is prescribed analytically, likely too
254 slow if the metric is rebuilt every timestep. Per the correction below, this is probably
255 unnecessary.
256- **Cap $L$.** All the error terms grow with $L$; float32 sets a practical ceiling that float64
257 would not.
259---
261## 5b. Summation order in the Legendre transform
263This is orthogonal to the 12→6 change, applies equally to the current code, and in float32 is
264probably worth more than the transform-count reduction. Do it first and independently, so its effect
265can be measured on its own.
267**Why.** Every $\varepsilon L$ and $\varepsilon L^2$ in §5 rides on the per-transform roundoff
268floor, and in float32 that floor is set by *how the sums are accumulated*, not by the mathematics.
269For each $(m, \theta_i)$ the synthesis evaluates
271$$u^m(\theta_i) = \sum_{\ell=|m|}^{L} u^m_\ell\,\bar P^m_\ell(\cos\theta_i),$$
273an $O(L)$-term sum. Error growth by accumulation strategy, for an $N$-term sum:
275| Strategy | Worst case | Typical (random signs) |
276|---|---|---|
277| Sequential | $\varepsilon N$ | $\varepsilon\sqrt{N}$ |
278| Pairwise / tree | $\varepsilon\log_2 N$ | $\varepsilon\sqrt{\log_2 N}$ |
279| Kahan compensated | $\varepsilon$ (+ $O(\varepsilon^2 N)$) | $\varepsilon$ |
281At $L=256$ in float32 that is the difference between $\sim1.5\times10^{-5}$ and $\sim5\times10^{-7}$
282per transform — more than an order of magnitude, for no change in operation count.
284**On GPU this may already be partly free.** A workgroup tree reduction over $\ell$ *is* pairwise
285summation. The failure mode is a serial `for` loop over $\ell$ inside a single thread, which is the
286natural way to write the shader if each thread owns one $(m,\theta_i)$ pair and is exactly the
287$\varepsilon N$ row above. Check which shape the kernel has before assuming anything.
289**Where it applies.**
291- Synthesis $\mathcal{S}$: the sum over $\ell$, as above. The $\varphi$-direction FFT is already
292 tree-structured and needs no attention.
293- Analysis $\mathcal{A}$: the quadrature sum over latitude nodes $\theta_i$ carries the identical
294 problem and the identical fix. It also matters more here, because this is the step relied on in
295 §5 to suppress the polar spike — a noisy quadrature sum weakens exactly the mechanism the
296 six-transform scheme depends on.
298**Practical notes.**
300- Blocked summation (accumulate in blocks of 8–32, then combine) captures most of the pairwise
301 benefit with a simpler kernel and better register behaviour than a full tree.
302- Kahan costs ~4 flops per term and is usually bandwidth-hidden on GPU; worth benchmarking rather
303 than assuming it is too expensive.
304- For $m>0$ near the poles, $\bar P^m_\ell(\cos\theta)$ spans many orders of magnitude across $\ell$.
305 Summing smallest-magnitude-first helps, and is nearly free here because the terms are already
306 roughly ordered by $\ell$.
307- Standard stable recurrences for $\bar P^m_\ell$ (and guarding their under/overflow in float32's
308 narrower exponent range) are a separate prerequisite — no summation strategy rescues inaccurate
309 Legendre values.
311**Measurement.** Transform a band-limited field forward then back and compare to the input, in
312float32, sweeping $L\in\{64,128,256\}$. Sequential accumulation shows error growing roughly linearly
313in $L$; pairwise shows near-flat growth. This isolates the transform floor from everything else in
314§7 and should be run before the validation gate there, since it sets the baseline that gate is
315measured against.
317### Measured (2026-08-04, Dawn/Metal, `scripts/sht-accuracy.ts`)
319The sweep was run and the summation-order changes tried. Outcome: **withdrawn — the floor here is
320not summation-limited.**
322| $L$ | grid | rel-$L_2$ roundtrip | worst degree |
323|---|---|---|---|
324| 63 | 64×128 | $3.4\times10^{-6}$ | $\ell=62$: $4.5\times10^{-6}$ |
325| 127 | 128×256 | $4.7\times10^{-6}$ | $\ell=110$: $6.4\times10^{-6}$ |
326| 255 | 256×512 | $1.1\times10^{-5}$ | $\ell=246$: $1.4\times10^{-5}$ |
328- **The analysis side already sums pairwise.** The quadrature over latitudes is a workgroup
329 tree/subgroup reduction (`leg_analys`); only the synthesis has the serial per-thread $\ell$-loop.
330- **Kahan is unavailable on WebGPU in practice.** Dawn/Metal compiles WGSL with fast-math: a probe
331 kernel evaluates $((10^8 + 1) - 10^8) - 1$ to $0$, so the compensation folds away and Kahan
332 compiles to plain summation (bit-identical results, verified).
333- **Blocked summation (B=16) in the synthesis $\ell$-loop moved nothing**: $1.128\times10^{-5}
334 \to 1.128\times10^{-5}$ at $L=255$ (low digits shift, confirming the reordering was live), while
335 costing ~5% per round trip at $L=255$. Reverted.
336- **Diagnosis:** the worst error concentrates at the top degrees — the signature of the Legendre
337 *recurrence* error (chains of length $\sim\ell$), not of $\ell$-uniform accumulation noise. This
338 is the "standard stable recurrences are a separate prerequisite" caveat above: the floor is set
339 by the accuracy of the $\bar P^m_\ell$ values themselves, and no summation strategy touches it.
340- The measured floor ($\sim\varepsilon L^{0.85}$, $1.1\times10^{-5}$ at $L=255$) is what the §7
341 validation gate should be read against.
343---
345## 6. Code changes
347| Location | Change |
348|---|---|
349| `src/surface_gradient/partial_derivatives` | Expose $\mathcal{S}(v^m_\ell)$ **pre-division** (flag or separate entry point). Needed by both the precompute and step 1. |
350| `SurfaceDiffOperator._precompute_metric_quantities()` | Return `p1, p2, q2, r` instead of `V_theta, V_phi`, per §3. |
351| `src/surface_screened_laplacian::surface_screened_laplacian()` | Replace body with §4. Both `for i in {x,y,z}` loops disappear. |
352| `SurfaceDiffOperator._precompute_curvature()`, Algorithm 3 | **Unchanged.** Still needs $X_{\theta\theta}, X_{\theta\varphi}, X_{\varphi\varphi}$ and the full Algorithm 1. |
353| `src/timestepping::make_implicit_op()`, Algorithm 5 | **Unchanged.** Only what line 8 calls changes. |
354| `src/real_embedding.py` | **Unchanged.** |
356The deprecated `SurfaceDiffOperator` methods for $\Delta_\Gamma$ and $(I + c\Delta_\Gamma)$ are the
357natural place to keep the old path as a reference implementation for the validation below.
359---
361### Correction: precompute conditioning
363An earlier draft claimed the polar relative error in $\tilde g_{\theta\theta}$ is $\varepsilon L^2$,
364making float64 precompute essential. That was wrong by a factor of $L$, in the safe direction.
365$\tilde g_{\theta\theta}$ is not synthesized directly; it is the square of $\sin\theta\,X_\theta$,
366which *is* synthesized, is $O(\sin\theta)$ at the poles, and carries absolute error $\varepsilon$ —
367so relative error $\varepsilon L$, preserved (up to a factor 2) by squaring. Same for
368$g_{\varphi\varphi}$ and $\tilde g_{\theta\varphi}$. The determinant combination is $O(\sin^4\theta)$
369and so are both of its terms, so there is no extra cancellation generically; $J$ inherits
370$\sim\varepsilon L$, i.e. $\sim10^{-5}$ at $L=128$ in float32. Acceptable.
372Caveat: this assumes the difference is not small compared to its terms, which fails if $X_\theta$
373becomes nearly parallel to $X_\varphi$ (near-degenerate parametrization). Worth a runtime check on
374$\det g$ if the surface can deform that far.
376Precompute error is also a *fixed* perturbation, identical every matvec, so it perturbs which
377operator is being solved but injects no noise into the Krylov space — GMRES converges normally.
379---
381## 7. Validation, in order
383At $\varepsilon=6\times10^{-8}$ there is no margin for the $\varepsilon\sqrt{L}$ suppression estimate
384in §5 to be off by an order of magnitude. Step 2 is a **gate**, not a confirmation.
3861. **Smoothness check (do this first).** On a deformed, non-axisymmetric surface, form $P$ and
387 $\tilde{Q}$ on the grid and compare their SH coefficient decay against $(\nabla_\Gamma u)_x$ from
388 the current code. Matching tails confirm both are genuinely smooth on $S^2$. If this fails,
389 nothing else is worth doing. Run this in float64 on CPU — it is a mathematical check, not a
390 precision one.
3912. **Coefficient-space diff in float32 at production $L$**, against a float64 CPU reference
392 implementation of Algorithm 4. Landing near $10^{-6}$ means the suppression argument holds.
393 Landing near $10^{-4}$ means the polar spike is surviving the analysis and the polar rings need
394 separate handling.
3953. **Sweep $L \in \{64,128,256\}$** and fit the growth exponent of (2). Flat-ish confirms
396 suppression; growth like $L^2$ means it is not working.
3974. **Grid-space max-norm diff near the poles.** Expect the extra power of $L$ here. If only this
398 grows and (2) stays flat, the scheme is fine for use inside the implicit solve.
3995. **GMRES iteration count and final achieved residual**, float32, versus the current code. The
400 operator is the same, so iterations should be unchanged; a stall above tolerance that does not
401 occur in float64 indicates the matvec noise floor is binding.
403---
405## 8. Suggestions considered and withdrawn
407- **Splitting $r$ across steps 3 and 6** to keep $\sin^1$ scaling. Does not work: $P/\sin\theta$ is
408 not smooth on $S^2$ (round sphere, $u = x$: $P = \sin\theta\cos\theta\cos\varphi$, so
409 $P/\sin\theta \to \cos\varphi$ at the pole). That $\sin\theta$ must stay in $r$. Algorithm 4 *can*
410 split because its intermediate — the Cartesian gradient — is smooth; that smoothness is precisely
411 what the six extra transforms buy.
412- **Weighting by $J$ to get an SPD operator and use PCG.** Avoiding the $1/\sqrt{g}$ makes
413 $\mathcal{L}$ self-adjoint, but the mass term becomes multiplication by $J$, costing its own
414 synthesis/analysis pair. A wash on transform count; worth it only if the CG properties themselves
415 are wanted. Discrete symmetry would also hold only to quadrature accuracy unless products are
416 dealiased (3/2 rule).
417- **Mixed precision (float64 for step 6's synthesis and the $r$ multiply only).** Unavailable:
418 WebGPU is float32-only. Superseded by the mitigations in §5a.
419- If the $\varepsilon L^2$ ever does bind, the remaining remedies are a shifted or uniform-in-$\theta$
420 latitude grid (removing the $O(L^{-2})$ node clustering) or a separate local formula for the polar
421 rings — both more work than the transform savings justify without a specific reason. In float32
422 the ceiling on $L$ may bind first and be the cheaper accommodation.
424## 9. Related: external vector transforms
426Steps 1–2 and 4–5 together are a vector/spin-weighted spherical harmonic transform. If SHTns
427(`spat_to_SHsphtor`, `SHsphtor_to_spat`) or SPHEREPACK (`gradgs`, `divgs`) can be linked, the
428gradient and divergence each become a single library call, the pole divisions are handled internally,
429and the hand-rolled $\alpha^\pm$ recurrences on this path are no longer needed.
431## 10. Beyond transform count
433The other lever is iteration count rather than cost per iteration. Since $M^{-1}A$ approaches
434multiplication by $1/J$ at high $\ell$, folding a mean or smoothed $J$ into the preconditioner could
435reduce GMRES iterations by more than any of the above reduces transforms.
437### Measured (2026-08-05, Richardson iteration, fp32)
439Implemented, with two corrections the measurements forced.
441**The right constant is the minimax over the symbol, not over $J$.** The high-$\ell$ per-mode
442factor is governed by the full principal symbol: in the orthonormal frame the symbol matrix is
443$S = (1/J)\begin{pmatrix} p_1 & p_2 \\ p_2 & q_2\end{pmatrix}$, whose eigenvalues $\mu(x)$ are the
444inverse squared principal stretches — direction matters. Preconditioning with $\lambda/\hat J$
445contracts every mode and direction iff $\hat J\mu \in (0,2)$, so
446$$\hat J = 2/(\mu_{\min} + \mu_{\max}), \qquad \text{rate} = (\mu_{\max}-\mu_{\min})/(\mu_{\max}+\mu_{\min}) < 1.$$
447The det-based mean of $J$ (this section's original suggestion; $\mu$'s geometric mean, exact only
448for conformal surfaces) is insufficient: on the shipped ellipsoid it leaves directional
449high-degree bands with amplification $> 1$ — patterns went qualitatively high-frequency at
450moderate settings and diverged as niter or $L$ grew. With the symbol-based constant every
451niter/geometry combination in the test sweep converges (peanut: $\mu \in [0.44, 6.2]$, plain rate
4525.2, preconditioned rate 0.87).
454**The correction must be band-projected.** Algorithm 5's "zero $\ell \ge L-2$" is load-bearing:
455without applying the same mask to the correction $d\Delta u$, the top two degrees iterate toward
456the *undiffused* right-hand side — each Richardson iteration strips more of their implicit
457diffusion, at species-dependent rates, manufacturing a spurious Turing band at the band edge
458(observed on the round sphere: top-degree energy growing $\sim 3\%$/step at $L=127$, 8 iterations).
460**Payoff shape:** on mildly deformed surfaces one iteration already reaches the $\sim10^{-4}$ fp32
461accumulation floor, so iteration counts do not drop — the speculation above does not hold at fp32.
462The gain is reach and correctness: stiff geometries and high niter/$L$ combinations that
463previously diverged (or silently shifted the pattern's wavelength) now converge with a
464resolution-independent spectrum.