Reduce the Laplace-Beltrami matvec to 6 transforms per species per iteration
Evaluate the implicit solve's geometric correction in flux form
(docs/reduced-transforms.md): the alpha-shift and i*m multiply are exposed
as coefficient-space externals (dthetac/dphic), the flux-metric weights
p1/p2/q2/r are precomputed in f64 alongside the inverse metric quantities,
and all three models' Richardson loops are rewritten to 3 syntheses +
3 analyses per species per iteration, versus Algorithm 4's 12 transforms.
~2x faster steps: 1.13 vs 2.22 ms/step at lmax 127, niter 2, on bumpy.
The 12-transform form ships as schnakenberg-alg4 for live A/B, and
test/fluxChecks.ts holds both forms to the same operator: f64 smoothness
of the analysed fluxes (exactly band-limited on the sphere, ~1e-13 tails
against the non-smooth control's ~1e-2; matching Cartesian-gradient tails
on bumpy), the 6-vs-12 transform count asserted from the compiled op
sequences, and state agreement through a real simulation (3.6e-4 after
20 steps at lmax 63).
Also scripts/sht-accuracy.ts, the fp32 round-trip sweep prescribed by the
doc's Sec 5b. The summation-order changes it was built to judge were
measured and withdrawn: the floor is set by the Legendre recurrence, not
accumulation order, and Kahan compensation folds to a no-op under the
Metal compiler's fast-math. Numbers in the doc's "Measured" section.
23 changed files+1488−246
README.mdmodified+84−64View file
@@ -9,15 +9,12 @@ This is the sibling of
99 solves the same systems on the round sphere. Everything there is here; what is
1010 added is a *surface*.
1111
12-> [!WARNING]
13-> **The geometry is rendered, not yet solved on.** The Laplace–Beltrami
14-> operator in the models is still the round sphere's — the term that carries
15-> the shape is a placeholder that is identically zero. On anything but the
16-> sphere you are looking at the sphere's pattern painted onto that surface, not
17-> the pattern that surface would grow. Everything the correction needs in order
18-> to be dropped in — the embedding, the split of the operator, the iterative
19-> solve, the unrolled loop — is built and tested. See
20-> [The geometry is not in the operator yet](#the-geometry-is-not-in-the-operator-yet).
12+The geometry is in the operator: the models evaluate the surface
13+Laplace–Beltrami operator `lap_g` inside the implicit solve, in a **flux form
14+that costs 6 spherical-harmonic transforms per species per iteration** where
15+the textbook Cartesian-gradient form needs 12. See
16+[The geometry in the operator](#the-geometry-in-the-operator) and
17+[docs/reduced-transforms.md](docs/reduced-transforms.md).
2118
2219 ## What a surface is here
2320
@@ -59,7 +56,7 @@ Four geometries ship: [sphere](geometries/sphere.m) (the reference case),
5956 whose waist is a saddle — and [bumpy](geometries/bumpy.m). Each is editable in
6057 the page, with its own parameters. Changing a shape does not recompile the
6158 solver and does not disturb the run: the geometry is data whose shape in the
62-bindings depends only on the grid, so a swap is six buffer writes and the
59+bindings depends only on the grid, so a swap is sixteen buffer writes and the
6360 pattern carries straight on.
6461
6562 A **morph** slider blends the drawn surface back to the unit sphere. The
@@ -103,57 +100,63 @@ Unew = (B + dt*D*dlap(Unew)) ./ (1 + dt*D*lam)
103100 and the loop iterates it from the round-sphere answer. That is preconditioned
104101 Richardson, with the operator we can invert exactly as the preconditioner; it
105102 converges while `dt*D*dlap` stays small against `(I - dt*D*lap_s)`, which is
106-what would keep the cost to a few transforms per step rather than a full
107-elliptic solve. Written out, the whole of
108-[`models/schnakenberg.m`](models/schnakenberg.m)'s step is:
103+what keeps the cost to a few transforms per step rather than a full elliptic
104+solve (see [docs/richardson-iteration.md](docs/richardson-iteration.md)). One
105+species of [`models/schnakenberg.m`](models/schnakenberg.m)'s solve loop:
109106
110107 ```matlab
111-function [Un, Vn, u, v] = step(U, V, lam, gx, gy, gz, a, b, D1, D2, dt, niter)
112- u = synth(U);
113- v = synth(V);
114- uuv = u .* u .* v;
115-
116- Bu = U + dt * analys(a - u + uuv);
117- Bv = V + dt * analys(b - uuv);
118-
119- Un = Bu ./ (1 + (dt * D1) * lam);
120- Vn = Bv ./ (1 + (dt * D2) * lam);
121-
122- for k = 1:niter
123- dLu = 0 * Un; % <- the placeholder
124- dLv = 0 * Vn;
125- Un = (Bu + (dt * D1) * dLu) ./ (1 + (dt * D1) * lam);
126- Vn = (Bv + (dt * D2) * dLv) ./ (1 + (dt * D2) * lam);
127- end
108+for k = 1:niter
109+ Fu = Un .* filt; % zero the top 2 degrees before differentiating
110+ Ftu = synth(dthetac(Fu)); % sin(theta) * dtheta(u) -- smooth on the sphere
111+ Fpu = synth(dphic(Fu)); % dphi(u) -- smooth on the sphere
112+ Pu = p1 .* Ftu + p2 .* Fpu; % the two fluxes, also smooth: the precomputed
113+ Qu = p2 .* Ftu + q2 .* Fpu; % weights carry every 1/sin(theta) there is
114+ Pcu = analys(Pu) .* filt;
115+ Qcu = analys(Qu) .* filt;
116+ scu = dthetac(Pcu) + dphic(Qcu); % divergence, in coefficient space
117+ lapu = r .* synth(scu); % = lap_g(u) on the grid
118+ dLu = analys(lapu) + lam .* Un; % dlap = lap_g - lap_s
119+ Un = (Bu + (dt * D1) * dLu) ./ (1 + (dt * D1) * lam);
128120 end
129121 ```
130122
131-Written this way rather than as a residual correction on purpose: with `dlap`
132-zero, every iterate is *bit for bit* the first line, with no cancellation to
133-round differently. So the sphere case is not "close to" turing-sphere, it is
134-the same arithmetic, and the tests assert exactly that — the state after 20
135-steps is identical at 0, 1 and 4 iterations.
136-
137-### The geometry is not in the operator yet
138-
139-What belongs where `dLu` is now is `dlap = lap_g - lap_s` applied to the current
140-iterate. Getting it needs two things this repo does not have:
141-
142-1. **The induced metric**, `g_ij = ∂_i X · ∂_j X` for `X = (gx, gy, gz)`. The
143- geometry is static and low-degree, so this is a one-off precomputation, not
144- per-step work — but it needs θ- and φ-derivatives of the embedding.
145-2. **Surface derivatives of the field**, per iteration. In the round frame this
146- is the spheroidal transform pair — SHTNS's `SHsph_to_spat` and
147- `spat_to_SHsph`, i.e. `grad_s` and `div_s` — which lets the operator be
148- written as `div_s(A grad_s f)` with `A` built from the metric, with no
149- explicit `1/sin θ` to go singular at the poles.
150-
151-Both need Legendre *derivative* tables, which the vendored WGSL transforms under
152-[`src/sht/`](src/sht/) do not implement — they are scalar synthesis and analysis
153-only. That is the missing piece, and it is a substantial addition to the
154-transforms rather than a change to the models. Until it lands, the models take
155-`gx, gy, gz` (the surface on the grid) and `Gx, Gy, Gz` (the same surface as
156-coefficients) as arguments and do not use them, and the app says so.
123+On the sphere `dlap` is mathematically zero — `p1 = q2 = 1`, `p2 = 0`,
124+`r = 1/sin²θ`, and the composition collapses to `lap_s` — so the sphere case
125+reproduces turing-sphere to fp32 round-off, and the tests assert the state
126+stays put across 0, 1 and 4 iterations.
127+
128+### The geometry in the operator
129+
130+`dlap = lap_g - lap_s` is applied to the current iterate at every solve
131+iteration, so its transform count is what the whole step's cost scales with.
132+Two formulations ship:
133+
134+1. **The flux form** (above, all three models): `lap_g u` as the weighted
135+ divergence of two weighted fluxes of the sin-scaled derivatives. The
136+ weights `p1, p2, q2, r` are grid arrays precomputed once per surface from
137+ the embedding's θ/φ tangents
138+ ([`src/geom/metric.ts`](src/geom/metric.ts)), chosen so that **every field
139+ that gets analysed is a smooth function on the sphere** — the property
140+ that makes spherical-harmonic analysis meaningful, and the entire
141+ difficulty near the poles. Cost: **6 transforms** per species per
142+ iteration (3 syntheses + 3 analyses; `dthetac`/`dphic` are O(nlm)
143+ coefficient shuffles, not transforms). The derivation, the smoothness
144+ argument and the fp32 error analysis are in
145+ [docs/reduced-transforms.md](docs/reduced-transforms.md).
146+2. **The Cartesian-gradient form** (Algorithm 4 of `docs/algos.pdf`), kept as
147+ a live reference in
148+ [`models/schnakenberg_alg4.m`](models/schnakenberg_alg4.m) and selectable
149+ in the app: the surface gradient carried as three ambient components
150+ through the inverse metric quantities `Vt*/Vp*`. Cost: **12 transforms**
151+ per species per iteration. The tests hold both forms to the same answer on
152+ a curved surface, and both metric formulations are precomputed and
153+ uploaded for every geometry, so either kind of model runs.
154+
155+The θ-derivative machinery both forms need — the α± recurrence
156+(`sin θ ∂θ Y_l^m = α⁺Y_{l+1}^m + α⁻Y_{l-1}^m`) as a coefficient-space shuffle
157+feeding the existing scalar synthesis — lives in
158+[`src/sht/deriv.ts`](src/sht/deriv.ts); no Legendre-derivative tables are
159+required.
157160
158161 ### `for` loops, unrolled
159162
@@ -183,8 +186,9 @@ Two consequences worth stating:
183186 each loop body assigns before the pass and refuses the ones that escape, so
184187 that case is a compile error rather than a stale read.
185188
186-Unrolling is exactly linear in the trip count: 2 GPU ops per species per
187-iteration, asserted in the tests.
189+Unrolling is exactly linear in the trip count: 19 GPU ops per species per
190+iteration (6 transforms, 4 coefficient shuffles, 9 kernels), asserted in the
191+tests.
188192
189193 ## MATLAB, compiled to WebGPU
190194
@@ -198,9 +202,9 @@ operations whose type rules numbl learns from a `.mtoc2.js` workspace file, and
198202 which the backend maps onto the spherical-harmonic pipelines. Anything it cannot
199203 express is refused at compile time with a source position.
200204
201-The Schnakenberg step above compiles to 17 GPU operations at one solve
202-iteration: 4 transforms, 11 generated kernels, and 2 buffer copies feeding the
203-new state back.
205+The Schnakenberg step compiles to 51 GPU operations at one solve iteration:
206+16 transforms, 8 coefficient-space shuffles, 25 generated kernels, and 2
207+buffer copies feeding the new state back.
204208
205209 Two consequences carried over:
206210
@@ -324,7 +328,7 @@ package alone. Its binaries need glibc 2.29+. Other flags: `--steps`,
324328
325329 There is no second implementation of the solver to diff against, so the `.m`
326330 path is checked against **closed-form answers** and against **exact structural
327-properties**. Four modules, run in both environments:
331+properties**. Five modules, run in both environments:
328332
329333 [`test/analyticChecks.ts`](test/analyticChecks.ts) — cases whose evolution is
330334 known exactly, run through the whole real pipeline. All three are statements
@@ -350,11 +354,27 @@ about the round sphere, so all three build on the sphere geometry:
350354 **the same coefficients give the same surface on a 2× grid** — the 2× Gauss
351355 latitudes share no point with the 1× ones, so agreeing there is agreeing
352356 everywhere, which is what "rendered exactly, not subdivided" means;
353-- unrolling is **exactly linear** in the trip count, and the state after 20 steps
354- is **bit-identical** at 0, 1 and 4 iterations;
357+- unrolling is **exactly linear** in the trip count, and on the sphere — where
358+ the geometric correction is mathematically zero — the state after 20 steps
359+ stays within fp32 round-off of the 0-iteration one at 1 and 4 iterations;
355360 - a runtime loop bound is refused at compile time;
356361 - swapping the surface mid-run leaves the spectral state untouched.
357362
363+[`test/fluxChecks.ts`](test/fluxChecks.ts) — the six-transform flux-form
364+Laplace-Beltrami scheme
365+([docs/reduced-transforms.md](docs/reduced-transforms.md)):
366+
367+- on the sphere, the precomputed weights match their closed form and the
368+ analysed fluxes are **exactly band-limited** (beyond-band tails at f64
369+ round-off, ~1e-13), while the deliberately non-smooth control
370+ `Q̃/sin θ` keeps a fat tail (~1e-2) — the discrimination the whole scheme
371+ rests on;
372+- on a non-axisymmetric surface, the flux tails match the Cartesian gradient
373+ component's, the doc's §7.1 criterion;
374+- the compiled op sequences add **6 transforms per species per iteration
375+ against Algorithm 4's 12**, and a real simulation driven by each stays
376+ within fp32 accumulation of the other.
377+
358378 [`test/modelChecks.ts`](test/modelChecks.ts) compiles every model the app offers
359379 and asserts **how many kernels it compiles to**, split into the base step and
360380 what one solve iteration adds. That is a fusion guard: if numbl's inline pass
docs/algos.pdfadded+0−0View file
Binary file not shown.
docs/reduced-transforms.mdadded+388−0View file
@@ -0,0 +1,388 @@
1+# Reducing spherical-harmonic transforms in $\Delta_\Gamma$
2+
3+**Summary.** Algorithm 4 costs 12 transforms per matvec (8 syntheses, 4 analyses). A flux-form
4+reformulation, with weights chosen so that every analyzed field is smooth on $S^2$, evaluates the
5+same operator in **6 transforms** (4 syntheses, 2 analyses). Notation follows `algos.pdf`.
6+
7+---
8+
9+## 1. Where the current cost comes from
10+
11+| Algorithm 4 line | Work | Transforms |
12+|---|---|---|
13+| 1 | $\partial_\theta u,\ \partial_\varphi u$ | 2 $\mathcal{S}$ |
14+| 5 | analysis of 3 Cartesian components of $\nabla_\Gamma u$ | 3 $\mathcal{A}$ |
15+| 6 | $\partial_\theta$ and $\partial_\varphi$ of each of those 3 components | 6 $\mathcal{S}$ |
16+| 8 | final analysis | 1 $\mathcal{A}$ |
17+| | | **12** |
18+
19+Two sources of waste:
20+
21+1. The gradient is carried in **ambient $\mathbb{R}^3$ components** — 3 fields for an intrinsically
22+ 2-dimensional object.
23+2. Line 6 takes **both** derivatives of **each** component, where the divergence needs only one
24+ derivative of each of two fluxes.
25+
26+There is also a possible free win independent of everything below: Algorithm 1 as specified returns
27+all five derivatives. The Laplacian path needs only $\partial_\theta u$ and $\partial_\varphi u$;
28+the second-derivative and mixed-derivative machinery is used exclusively by Algorithm 3 (curvature).
29+If `surface_screened_laplacian` calls Algorithm 1 wholesale, it is doing 5 syntheses where 2 suffice
30+at lines 1 and 6.
31+
32+---
33+
34+## 2. The constraint that shapes the solution
35+
36+The Cartesian design in Algorithm 4 exists to avoid pole singularities, and it is correct to do so.
37+The relevant property is smoothness **as a scalar function on $S^2$**, since that is what controls
38+SH coefficient decay and hence whether $\mathcal{A}$ is meaningful.
39+
40+| Quantity | Smooth on $S^2$? |
41+|---|---|
42+| $\partial_\varphi u$ | yes — exactly band-limited, eq. (2.2) |
43+| $\sin\theta\,\partial_\theta u$ | yes — exactly band-limited, eq. (2.4) |
44+| $\partial_\theta u$ | **no** — bounded, but $\varphi$-dependent limit at the poles |
45+| $g_{\theta\theta},\ g^{\theta\theta},\ V_\theta,\ V_\varphi$ | **no** in general |
46+| $(\nabla_\Gamma u)_x,\ (\nabla_\Gamma u)_y,\ (\nabla_\Gamma u)_z$ | yes |
47+
48+Concretely, for the ellipsoid $X = (a\sin\theta\cos\varphi,\ b\sin\theta\sin\varphi,\ c\cos\theta)$,
49+$|X_\theta|^2 \to a^2\cos^2\varphi + b^2\sin^2\varphi$ as $\theta\to 0$: no limit exists.
50+
51+Algorithm 4 never analyzes anything in the "no" rows — the non-smooth quantities appear only as
52+pointwise grid factors. **Any replacement must preserve this property.** The naive flux form
53+$P = \sqrt{g}\,(g^{\theta\theta}u_\theta + g^{\theta\varphi}u_\varphi)$,
54+$Q = \sqrt{g}\,(g^{\varphi\theta}u_\theta + g^{\varphi\varphi}u_\varphi)$ does not:
55+on the round sphere with $u = x$, $Q = -\sin\varphi$, which is not a function on $S^2$.
56+
57+### The correct weighting
58+
59+$P$ and $Q$ are $\sqrt{g}$ times the contravariant components of $G := \nabla_\Gamma u$. Using
60+$\det[X_\theta, X_\varphi, n] = -\sqrt{g}$:
61+
62+$$P = -\,G\cdot(X_\varphi \times n), \qquad \sin\theta\,Q = -\,G\cdot\big(n \times \sin\theta\,X_\theta\big).$$
63+
64+Every factor on the right is smooth on $S^2$: $G$ smooth, $n$ smooth, $X_\varphi$ smooth, and
65+$\sin\theta\,X_\theta$ smooth because it is exactly band-limited by the recurrence already
66+implemented. Hence $P$ and $\tilde{Q} := \sin\theta\,Q$ are analyzable, on the same footing and for
67+the same structural reason as the Cartesian gradient components.
68+
69+The cross products are the smoothness certificate only — they are not needed in the code.
70+
71+---
72+
73+## 3. Precompute (once per surface, grid space)
74+
75+Replaces `_precompute_metric_quantities()`. From the embedding coefficients $\hat{X}^m_\ell$, obtain
76+$X_\varphi$ and $\sin\theta\,X_\theta$ componentwise via Algorithm 1, where the latter is the
77+**undivided** output of Algorithm 1 line 4, i.e. $\mathcal{S}(v^m_\ell)$ with no $/\sin\theta$.
78+Then, pointwise:
79+
80+$$\tilde{g}_{\theta\theta} := |\sin\theta\,X_\theta|^2, \qquad
81+ \tilde{g}_{\theta\varphi} := (\sin\theta\,X_\theta)\cdot X_\varphi, \qquad
82+ g_{\varphi\varphi} := |X_\varphi|^2$$
83+
84+$$J := \frac{\sqrt{\tilde{g}_{\theta\theta}\,g_{\varphi\varphi} - \tilde{g}_{\theta\varphi}^{\,2}}}{\sin^2\theta}
85+ \qquad\text{so that } \sqrt{\det g} = J\sin\theta$$
86+
87+(The radicand is $\sin^2\theta\det g = J^2\sin^4\theta$, so the square root is $J\sin^2\theta$ — hence
88+$\sin^2\theta$, not $\sin\theta$, in the denominator.)
89+
90+Store four scalar grid arrays:
91+
92+$$p_1 = \frac{g_{\varphi\varphi}}{J\sin^2\theta}, \qquad
93+ p_2 = -\frac{\tilde{g}_{\theta\varphi}}{J\sin^2\theta}, \qquad
94+ q_2 = \frac{\tilde{g}_{\theta\theta}}{J\sin^2\theta}, \qquad
95+ r = \frac{1}{J\sin^2\theta}$$
96+
97+All four are bounded: the $\sin^2\theta$ denominators cancel against vanishing numerators
98+($\tilde{g}_{\theta\varphi} = O(\sin^2\theta)$, $g_{\varphi\varphi} = O(\sin^2\theta)$), the same
99+finite limits the current $V_\theta, V_\varphi$ have. Note $p_1, p_2, q_2$ are bounded where
100+$V_\varphi = O(1/\sin\theta)$ is not.
101+
102+**Three scalar arrays replace the six components of $V_\theta, V_\varphi$.** The surface
103+representation is unchanged: $X_\theta, X_\varphi$ still come componentwise from $\hat{X}^m_\ell$
104+via Algorithm 1.
105+
106+---
107+
108+## 4. The per-matvec algorithm
109+
110+Input $\{u^m_\ell\}$; output $\{(\Delta_\Gamma u)^m_\ell\}$.
111+
112+| # | Step | Transforms |
113+|---|---|---|
114+| 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}$ |
115+| 2 | $B \leftarrow \mathcal{S}(im\,u^m_\ell)$ | $\mathcal{S}$ |
116+| 3 | $P \leftarrow p_1 A + p_2 B$, $\tilde{Q} \leftarrow p_2 A + q_2 B$ — pointwise | — |
117+| 4 | $\hat{P} \leftarrow \mathcal{A}(P)$, $\hat{\tilde{Q}} \leftarrow \mathcal{A}(\tilde{Q})$ | 2 $\mathcal{A}$ |
118+| 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$ | — |
119+| 6 | $\Delta_\Gamma u \leftarrow r \cdot \mathcal{S}(s^m_\ell)$ | $\mathcal{S}$ |
120+| 7 | $\{(\Delta_\Gamma u)^m_\ell\} \leftarrow \mathcal{A}(\Delta_\Gamma u)$; zero $\ell \ge L-2$ | $\mathcal{A}$ |
121+
122+**4 syntheses + 2 analyses = 6**, versus 12.
123+
124+- $A$ and $B$ are exactly $\sin\theta\,\partial_\theta u$ and $\partial_\varphi u$.
125+- Steps 1 and 5 use the **same** precomputed $\alpha^\pm$ table; step 5 is the adjoint-style reuse
126+ of the shift already implemented for step 1. Adding the two flux contributions in coefficient
127+ space before synthesizing is what saves the final pair of transforms.
128+- The **only** division by $\sin\theta$ anywhere is folded into $p_1, p_2, q_2, r$ at precompute
129+ time. The per-matvec path contains none.
130+
131+---
132+
133+## 5. Numerical trade-off
134+
135+Both schemes contain two powers of $\sin\theta$ division in total. What differs is **placement**.
136+
137+- **Algorithm 4** spends them in separate stages, one before line 5 and one after. The intervening
138+ analysis suppresses the polar spike: Gauss–Legendre weights give $w_1 = O(L^{-2})$ at the polar
139+ ring, so a grid error of $\varepsilon L$ there contributes
140+ $\sim L^{-2}\cdot L^{1/2}\cdot \varepsilon L = \varepsilon L^{-1/2}$ to any coefficient. Stage 2
141+ then starts from clean coefficients and incurs a *fresh* $\varepsilon L$. The two amplifications
142+ never multiply. Net grid-space relative error: $\varepsilon L$.
143+- **The new scheme** has no division at all through step 5, then pays for both powers at once in
144+ $r = O(L^2)$ at step 6 — one event, with no intervening analysis to break it in half. Net
145+ grid-space relative error: $\varepsilon L^2$.
146+
147+The mechanism: with $N \approx L+1$ Gauss–Legendre nodes, $1 - x_1 = O(N^{-2})$ so
148+$\sin\theta_1 = O(N^{-1})$. Since $s = J\sin^2\theta\,\Delta_\Gamma u$ is $O(L^{-2})$ at the polar
149+ring but $O(1)$ over the bulk, and synthesis commits roundoff scaled by the field's *global* size at
150+every node alike, multiplying by $r \sim L^2$ recovers the signal and inflates the noise.
151+
152+| | divisions | placement | grid-space relative error |
153+|---|---|---|---|
154+| Algorithm 4 | $\sin\theta$, $\sin\theta$ | separated by $\mathcal{A}$ | $\varepsilon L$ |
155+| Six-transform | $\sin^2\theta$ | all at the end | $\varepsilon L^2$ |
156+
157+**This likely does not reach the returned coefficients.** Step 7's analysis suppresses the spike
158+exactly as line 5 does today: $L^{-2}\cdot L^{1/2}\cdot\varepsilon L^2 = \varepsilon L^{1/2}$,
159+comparable to the ordinary $\varepsilon\sqrt{L}$ accumulation of a transform pair — and the new
160+scheme runs half as many transforms, lowering that baseline. Inside the implicit solve, GMRES sees
161+only coefficients, so the extra power should be invisible.
162+
163+It matters only if grid values of $\Delta_\Gamma u$ are consumed directly: a nonlinear reaction
164+term, max-norm diagnostics, or an adaptive error estimator.
165+
166+Algorithm 1 line 7 already divides by $\sin^2\theta$, so the code is exposed to $\varepsilon L^2$
167+today — just on the second-derivative path, which the Laplacian never touches.
168+
169+---
170+
171+## 5a. float32 / WebGPU
172+
173+Target is WebGPU, which is float32-only: $\varepsilon = 2^{-24} \approx 6\times10^{-8}$ (spacing
174+$2^{-23} \approx 1.2\times10^{-7}$). No float64 fallback exists on device. All estimates in §5 are
175+linear in $\varepsilon$, so they scale directly:
176+
177+| | $\varepsilon\sqrt{L}$ (coeffs) | $\varepsilon L$ (Alg. 4 grid) | $\varepsilon L^2$ (new, grid) |
178+|---|---|---|---|
179+| $L=64$ | $5\times10^{-7}$ | $4\times10^{-6}$ | $2\times10^{-4}$ |
180+| $L=128$ | $7\times10^{-7}$ | $8\times10^{-6}$ | $1\times10^{-3}$ |
181+| $L=256$ | $1\times10^{-6}$ | $1.5\times10^{-5}$ | $4\times10^{-3}$ |
182+
183+**Coefficient space is fine** (~$10^{-6}$), which is the floor a float32 iterative solve sits at
184+anyway. **Grid space is not**: 0.1–0.4% relative on the polar rings at $L\ge128$. For a
185+reaction–diffusion solver this matters only if grid-space $\Delta_\Gamma u$ is consumed outside the
186+matvec. If the IMEX splitting evaluates $f(u)$ from $u$ on the grid (typical), it never is.
187+
188+**Two float32-specific arguments in favour of the new scheme:**
189+
190+- Baseline SHT roundoff accumulates per transform ($\sim\varepsilon\sqrt{L}$ to $\varepsilon L$
191+ each). Running 6 transforms instead of 12 halves that accumulation. On the coefficient-space error
192+ GMRES actually sees, this plausibly outweighs the polar term — the new scheme may be *more*
193+ accurate end-to-end in float32. Not asserted without measurement.
194+- 3 weight arrays instead of 6 halves per-matvec texture/buffer traffic. On GPU that is often the
195+ real bottleneck, independent of arithmetic.
196+
197+**Mitigations available without float64:**
198+
199+- **Pairwise or blocked summation in the Legendre sum over $\ell$.** The single highest-value
200+ float32 change, and it benefits the existing code too. See §5b.
201+- **Double-float (`f32x2`) arithmetic** for the pointwise steps 3 and 6 if needed — cheap, no
202+ transforms involved. Does not help with transform roundoff, which is the dominant term, so try
203+ summation order first.
204+- **CPU precompute in float64.** JS `Number` is float64, so §3 can run on the CPU regardless of
205+ WebGPU's limits, with float32 weights uploaded. Cost is CPU-side SHTs plus upload, paid once per
206+ surface update; viable if the surface evolves slowly or is prescribed analytically, likely too
207+ slow if the metric is rebuilt every timestep. Per the correction below, this is probably
208+ unnecessary.
209+- **Cap $L$.** All the error terms grow with $L$; float32 sets a practical ceiling that float64
210+ would not.
211+
212+---
213+
214+## 5b. Summation order in the Legendre transform
215+
216+This is orthogonal to the 12→6 change, applies equally to the current code, and in float32 is
217+probably worth more than the transform-count reduction. Do it first and independently, so its effect
218+can be measured on its own.
219+
220+**Why.** Every $\varepsilon L$ and $\varepsilon L^2$ in §5 rides on the per-transform roundoff
221+floor, and in float32 that floor is set by *how the sums are accumulated*, not by the mathematics.
222+For each $(m, \theta_i)$ the synthesis evaluates
223+
224+$$u^m(\theta_i) = \sum_{\ell=|m|}^{L} u^m_\ell\,\bar P^m_\ell(\cos\theta_i),$$
225+
226+an $O(L)$-term sum. Error growth by accumulation strategy, for an $N$-term sum:
227+
228+| Strategy | Worst case | Typical (random signs) |
229+|---|---|---|
230+| Sequential | $\varepsilon N$ | $\varepsilon\sqrt{N}$ |
231+| Pairwise / tree | $\varepsilon\log_2 N$ | $\varepsilon\sqrt{\log_2 N}$ |
232+| Kahan compensated | $\varepsilon$ (+ $O(\varepsilon^2 N)$) | $\varepsilon$ |
233+
234+At $L=256$ in float32 that is the difference between $\sim1.5\times10^{-5}$ and $\sim5\times10^{-7}$
235+per transform — more than an order of magnitude, for no change in operation count.
236+
237+**On GPU this may already be partly free.** A workgroup tree reduction over $\ell$ *is* pairwise
238+summation. The failure mode is a serial `for` loop over $\ell$ inside a single thread, which is the
239+natural way to write the shader if each thread owns one $(m,\theta_i)$ pair and is exactly the
240+$\varepsilon N$ row above. Check which shape the kernel has before assuming anything.
241+
242+**Where it applies.**
243+
244+- Synthesis $\mathcal{S}$: the sum over $\ell$, as above. The $\varphi$-direction FFT is already
245+ tree-structured and needs no attention.
246+- Analysis $\mathcal{A}$: the quadrature sum over latitude nodes $\theta_i$ carries the identical
247+ problem and the identical fix. It also matters more here, because this is the step relied on in
248+ §5 to suppress the polar spike — a noisy quadrature sum weakens exactly the mechanism the
249+ six-transform scheme depends on.
250+
251+**Practical notes.**
252+
253+- Blocked summation (accumulate in blocks of 8–32, then combine) captures most of the pairwise
254+ benefit with a simpler kernel and better register behaviour than a full tree.
255+- Kahan costs ~4 flops per term and is usually bandwidth-hidden on GPU; worth benchmarking rather
256+ than assuming it is too expensive.
257+- For $m>0$ near the poles, $\bar P^m_\ell(\cos\theta)$ spans many orders of magnitude across $\ell$.
258+ Summing smallest-magnitude-first helps, and is nearly free here because the terms are already
259+ roughly ordered by $\ell$.
260+- Standard stable recurrences for $\bar P^m_\ell$ (and guarding their under/overflow in float32's
261+ narrower exponent range) are a separate prerequisite — no summation strategy rescues inaccurate
262+ Legendre values.
263+
264+**Measurement.** Transform a band-limited field forward then back and compare to the input, in
265+float32, sweeping $L\in\{64,128,256\}$. Sequential accumulation shows error growing roughly linearly
266+in $L$; pairwise shows near-flat growth. This isolates the transform floor from everything else in
267+§7 and should be run before the validation gate there, since it sets the baseline that gate is
268+measured against.
269+
270+### Measured (2026-08-04, Dawn/Metal, `scripts/sht-accuracy.ts`)
271+
272+The sweep was run and the summation-order changes tried. Outcome: **withdrawn — the floor here is
273+not summation-limited.**
274+
275+| $L$ | grid | rel-$L_2$ roundtrip | worst degree |
276+|---|---|---|---|
277+| 63 | 64×128 | $3.4\times10^{-6}$ | $\ell=62$: $4.5\times10^{-6}$ |
278+| 127 | 128×256 | $4.7\times10^{-6}$ | $\ell=110$: $6.4\times10^{-6}$ |
279+| 255 | 256×512 | $1.1\times10^{-5}$ | $\ell=246$: $1.4\times10^{-5}$ |
280+
281+- **The analysis side already sums pairwise.** The quadrature over latitudes is a workgroup
282+ tree/subgroup reduction (`leg_analys`); only the synthesis has the serial per-thread $\ell$-loop.
283+- **Kahan is unavailable on WebGPU in practice.** Dawn/Metal compiles WGSL with fast-math: a probe
284+ kernel evaluates $((10^8 + 1) - 10^8) - 1$ to $0$, so the compensation folds away and Kahan
285+ compiles to plain summation (bit-identical results, verified).
286+- **Blocked summation (B=16) in the synthesis $\ell$-loop moved nothing**: $1.128\times10^{-5}
287+ \to 1.128\times10^{-5}$ at $L=255$ (low digits shift, confirming the reordering was live), while
288+ costing ~5% per round trip at $L=255$. Reverted.
289+- **Diagnosis:** the worst error concentrates at the top degrees — the signature of the Legendre
290+ *recurrence* error (chains of length $\sim\ell$), not of $\ell$-uniform accumulation noise. This
291+ is the "standard stable recurrences are a separate prerequisite" caveat above: the floor is set
292+ by the accuracy of the $\bar P^m_\ell$ values themselves, and no summation strategy touches it.
293+- The measured floor ($\sim\varepsilon L^{0.85}$, $1.1\times10^{-5}$ at $L=255$) is what the §7
294+ validation gate should be read against.
295+
296+---
297+
298+## 6. Code changes
299+
300+| Location | Change |
301+|---|---|
302+| `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. |
303+| `SurfaceDiffOperator._precompute_metric_quantities()` | Return `p1, p2, q2, r` instead of `V_theta, V_phi`, per §3. |
304+| `src/surface_screened_laplacian::surface_screened_laplacian()` | Replace body with §4. Both `for i in {x,y,z}` loops disappear. |
305+| `SurfaceDiffOperator._precompute_curvature()`, Algorithm 3 | **Unchanged.** Still needs $X_{\theta\theta}, X_{\theta\varphi}, X_{\varphi\varphi}$ and the full Algorithm 1. |
306+| `src/timestepping::make_implicit_op()`, Algorithm 5 | **Unchanged.** Only what line 8 calls changes. |
307+| `src/real_embedding.py` | **Unchanged.** |
308+
309+The deprecated `SurfaceDiffOperator` methods for $\Delta_\Gamma$ and $(I + c\Delta_\Gamma)$ are the
310+natural place to keep the old path as a reference implementation for the validation below.
311+
312+---
313+
314+### Correction: precompute conditioning
315+
316+An earlier draft claimed the polar relative error in $\tilde g_{\theta\theta}$ is $\varepsilon L^2$,
317+making float64 precompute essential. That was wrong by a factor of $L$, in the safe direction.
318+$\tilde g_{\theta\theta}$ is not synthesized directly; it is the square of $\sin\theta\,X_\theta$,
319+which *is* synthesized, is $O(\sin\theta)$ at the poles, and carries absolute error $\varepsilon$ —
320+so relative error $\varepsilon L$, preserved (up to a factor 2) by squaring. Same for
321+$g_{\varphi\varphi}$ and $\tilde g_{\theta\varphi}$. The determinant combination is $O(\sin^4\theta)$
322+and so are both of its terms, so there is no extra cancellation generically; $J$ inherits
323+$\sim\varepsilon L$, i.e. $\sim10^{-5}$ at $L=128$ in float32. Acceptable.
324+
325+Caveat: this assumes the difference is not small compared to its terms, which fails if $X_\theta$
326+becomes nearly parallel to $X_\varphi$ (near-degenerate parametrization). Worth a runtime check on
327+$\det g$ if the surface can deform that far.
328+
329+Precompute error is also a *fixed* perturbation, identical every matvec, so it perturbs which
330+operator is being solved but injects no noise into the Krylov space — GMRES converges normally.
331+
332+---
333+
334+## 7. Validation, in order
335+
336+At $\varepsilon=6\times10^{-8}$ there is no margin for the $\varepsilon\sqrt{L}$ suppression estimate
337+in §5 to be off by an order of magnitude. Step 2 is a **gate**, not a confirmation.
338+
339+1. **Smoothness check (do this first).** On a deformed, non-axisymmetric surface, form $P$ and
340+ $\tilde{Q}$ on the grid and compare their SH coefficient decay against $(\nabla_\Gamma u)_x$ from
341+ the current code. Matching tails confirm both are genuinely smooth on $S^2$. If this fails,
342+ nothing else is worth doing. Run this in float64 on CPU — it is a mathematical check, not a
343+ precision one.
344+2. **Coefficient-space diff in float32 at production $L$**, against a float64 CPU reference
345+ implementation of Algorithm 4. Landing near $10^{-6}$ means the suppression argument holds.
346+ Landing near $10^{-4}$ means the polar spike is surviving the analysis and the polar rings need
347+ separate handling.
348+3. **Sweep $L \in \{64,128,256\}$** and fit the growth exponent of (2). Flat-ish confirms
349+ suppression; growth like $L^2$ means it is not working.
350+4. **Grid-space max-norm diff near the poles.** Expect the extra power of $L$ here. If only this
351+ grows and (2) stays flat, the scheme is fine for use inside the implicit solve.
352+5. **GMRES iteration count and final achieved residual**, float32, versus the current code. The
353+ operator is the same, so iterations should be unchanged; a stall above tolerance that does not
354+ occur in float64 indicates the matvec noise floor is binding.
355+
356+---
357+
358+## 8. Suggestions considered and withdrawn
359+
360+- **Splitting $r$ across steps 3 and 6** to keep $\sin^1$ scaling. Does not work: $P/\sin\theta$ is
361+ not smooth on $S^2$ (round sphere, $u = x$: $P = \sin\theta\cos\theta\cos\varphi$, so
362+ $P/\sin\theta \to \cos\varphi$ at the pole). That $\sin\theta$ must stay in $r$. Algorithm 4 *can*
363+ split because its intermediate — the Cartesian gradient — is smooth; that smoothness is precisely
364+ what the six extra transforms buy.
365+- **Weighting by $J$ to get an SPD operator and use PCG.** Avoiding the $1/\sqrt{g}$ makes
366+ $\mathcal{L}$ self-adjoint, but the mass term becomes multiplication by $J$, costing its own
367+ synthesis/analysis pair. A wash on transform count; worth it only if the CG properties themselves
368+ are wanted. Discrete symmetry would also hold only to quadrature accuracy unless products are
369+ dealiased (3/2 rule).
370+- **Mixed precision (float64 for step 6's synthesis and the $r$ multiply only).** Unavailable:
371+ WebGPU is float32-only. Superseded by the mitigations in §5a.
372+- If the $\varepsilon L^2$ ever does bind, the remaining remedies are a shifted or uniform-in-$\theta$
373+ latitude grid (removing the $O(L^{-2})$ node clustering) or a separate local formula for the polar
374+ rings — both more work than the transform savings justify without a specific reason. In float32
375+ the ceiling on $L$ may bind first and be the cheaper accommodation.
376+
377+## 9. Related: external vector transforms
378+
379+Steps 1–2 and 4–5 together are a vector/spin-weighted spherical harmonic transform. If SHTns
380+(`spat_to_SHsphtor`, `SHsphtor_to_spat`) or SPHEREPACK (`gradgs`, `divgs`) can be linked, the
381+gradient and divergence each become a single library call, the pole divisions are handled internally,
382+and the hand-rolled $\alpha^\pm$ recurrences on this path are no longer needed.
383+
384+## 10. Beyond transform count
385+
386+The other lever is iteration count rather than cost per iteration. Since $M^{-1}A$ approaches
387+multiplication by $1/J$ at high $\ell$, folding a mean or smoothed $J$ into the preconditioner could
388+reduce GMRES iterations by more than any of the above reduces transforms.
models/allencahn.mmodified+12−23View file
@@ -9,36 +9,25 @@ function [U, u] = init(noise)
99 u = synth(U);
1010 end
1111
12-function [Un, u] = step(U, lam, filt, gx, gy, gz, Vtx, Vty, Vtz, Vpx, Vpy, Vpz, eps2, dt, niter)
12+function [Un, u] = step(U, lam, filt, gx, gy, gz, p1, p2, q2, r, eps2, dt, niter)
1313 u = synth(U);
1414
1515 Bu = U + dt * analys(u - u.^3);
1616 Un = Bu ./ (1 + (dt * eps2) * lam);
1717
1818 for k = 1:niter
19- % dlap = lap_g - lap_s, evaluated at the current iterate (see
20- % models/schnakenberg.m and docs/richardson-iteration.md for the
21- % derivation).
19+ % dlap = lap_g - lap_s, evaluated at the current iterate in flux form
20+ % (see models/schnakenberg.m, docs/richardson-iteration.md and
21+ % docs/reduced-transforms.md for the derivation).
2222 Fu = Un .* filt;
23- Ftu = dtheta(Fu);
24- Fpu = dphi(Fu);
25- dux = Ftu .* Vtx + Fpu .* Vpx;
26- duy = Ftu .* Vty + Fpu .* Vpy;
27- duz = Ftu .* Vtz + Fpu .* Vpz;
28- cux = analys(dux) .* filt;
29- cuy = analys(duy) .* filt;
30- cuz = analys(duz) .* filt;
31- Ftcux = dtheta(cux);
32- Fpcux = dphi(cux);
33- Ftcuy = dtheta(cuy);
34- Fpcuy = dphi(cuy);
35- Ftcuz = dtheta(cuz);
36- Fpcuz = dphi(cuz);
37- lapu = Ftcux .* Vtx + Fpcux .* Vpx;
38- lapu = lapu + Ftcuy .* Vty;
39- lapu = lapu + Fpcuy .* Vpy;
40- lapu = lapu + Ftcuz .* Vtz;
41- lapu = lapu + Fpcuz .* Vpz;
23+ Ftu = synth(dthetac(Fu));
24+ Fpu = synth(dphic(Fu));
25+ Pu = p1 .* Ftu + p2 .* Fpu;
26+ Qu = p2 .* Ftu + q2 .* Fpu;
27+ Pcu = analys(Pu) .* filt;
28+ Qcu = analys(Qu) .* filt;
29+ scu = dthetac(Pcu) + dphic(Qcu);
30+ lapu = r .* synth(scu);
4231 dLu = analys(lapu) + lam .* Un;
4332
4433 Un = (Bu + (dt * eps2) * dLu) ./ (1 + (dt * eps2) * lam);
models/brusselator.mmodified+20−42View file
@@ -12,7 +12,7 @@ function [U, V, u, v] = init(noise, A, B)
1212 v = synth(V);
1313 end
1414
15-function [Un, Vn, u, v] = step(U, V, lam, filt, gx, gy, gz, Vtx, Vty, Vtz, Vpx, Vpy, Vpz, A, B, D1, D2, dt, niter)
15+function [Un, Vn, u, v] = step(U, V, lam, filt, gx, gy, gz, p1, p2, q2, r, A, B, D1, D2, dt, niter)
1616 u = synth(U);
1717 v = synth(V);
1818 uuv = u .* u .* v;
@@ -24,51 +24,29 @@ function [Un, Vn, u, v] = step(U, V, lam, filt, gx, gy, gz, Vtx, Vty, Vtz, Vpx,
2424 Vn = Bv ./ (1 + (dt * D2) * lam);
2525
2626 for k = 1:niter
27- % dlap = lap_g - lap_s, evaluated at the current iterate (see
28- % models/schnakenberg.m and docs/richardson-iteration.md for the
29- % derivation).
27+ % dlap = lap_g - lap_s, evaluated at the current iterate in flux form
28+ % (see models/schnakenberg.m, docs/richardson-iteration.md and
29+ % docs/reduced-transforms.md for the derivation).
3030 Fu = Un .* filt;
31- Ftu = dtheta(Fu);
32- Fpu = dphi(Fu);
33- dux = Ftu .* Vtx + Fpu .* Vpx;
34- duy = Ftu .* Vty + Fpu .* Vpy;
35- duz = Ftu .* Vtz + Fpu .* Vpz;
36- cux = analys(dux) .* filt;
37- cuy = analys(duy) .* filt;
38- cuz = analys(duz) .* filt;
39- Ftcux = dtheta(cux);
40- Fpcux = dphi(cux);
41- Ftcuy = dtheta(cuy);
42- Fpcuy = dphi(cuy);
43- Ftcuz = dtheta(cuz);
44- Fpcuz = dphi(cuz);
45- lapu = Ftcux .* Vtx + Fpcux .* Vpx;
46- lapu = lapu + Ftcuy .* Vty;
47- lapu = lapu + Fpcuy .* Vpy;
48- lapu = lapu + Ftcuz .* Vtz;
49- lapu = lapu + Fpcuz .* Vpz;
31+ Ftu = synth(dthetac(Fu));
32+ Fpu = synth(dphic(Fu));
33+ Pu = p1 .* Ftu + p2 .* Fpu;
34+ Qu = p2 .* Ftu + q2 .* Fpu;
35+ Pcu = analys(Pu) .* filt;
36+ Qcu = analys(Qu) .* filt;
37+ scu = dthetac(Pcu) + dphic(Qcu);
38+ lapu = r .* synth(scu);
5039 dLu = analys(lapu) + lam .* Un;
5140
5241 Fv = Vn .* filt;
53- Ftv = dtheta(Fv);
54- Fpv = dphi(Fv);
55- dvx = Ftv .* Vtx + Fpv .* Vpx;
56- dvy = Ftv .* Vty + Fpv .* Vpy;
57- dvz = Ftv .* Vtz + Fpv .* Vpz;
58- cvx = analys(dvx) .* filt;
59- cvy = analys(dvy) .* filt;
60- cvz = analys(dvz) .* filt;
61- Ftcvx = dtheta(cvx);
62- Fpcvx = dphi(cvx);
63- Ftcvy = dtheta(cvy);
64- Fpcvy = dphi(cvy);
65- Ftcvz = dtheta(cvz);
66- Fpcvz = dphi(cvz);
67- lapv = Ftcvx .* Vtx + Fpcvx .* Vpx;
68- lapv = lapv + Ftcvy .* Vty;
69- lapv = lapv + Fpcvy .* Vpy;
70- lapv = lapv + Ftcvz .* Vtz;
71- lapv = lapv + Fpcvz .* Vpz;
42+ Ftv = synth(dthetac(Fv));
43+ Fpv = synth(dphic(Fv));
44+ Pv = p1 .* Ftv + p2 .* Fpv;
45+ Qv = p2 .* Ftv + q2 .* Fpv;
46+ Pcv = analys(Pv) .* filt;
47+ Qcv = analys(Qv) .* filt;
48+ scv = dthetac(Pcv) + dphic(Qcv);
49+ lapv = r .* synth(scv);
7250 dLv = analys(lapv) + lam .* Vn;
7351
7452 Un = (Bu + (dt * D1) * dLu) ./ (1 + (dt * D1) * lam);
models/schnakenberg.mmodified+35−47View file
@@ -8,6 +8,12 @@
88 % spherical-harmonic space (eigenvalues -lam), and the loop iterates the
99 % geometric correction dlap from that exact solve. Grid fields are npts x 1;
1010 % spectral fields are real 2 x nlm. See docs/richardson-iteration.md.
11+%
12+% The correction evaluates lap_g in flux form -- 6 transforms per species
13+% per iteration where the Cartesian-gradient form (Algorithm 4 of
14+% evolving_surface/notes/algos.tex) needs 12. See
15+% docs/reduced-transforms.md, and models/schnakenberg_alg4.m
16+% for the original form kept as a live reference.
1117
1218 function [U, V, u, v] = init(noise, a, b)
1319 us = a + b;
@@ -18,7 +24,7 @@ function [U, V, u, v] = init(noise, a, b)
1824 v = synth(V);
1925 end
2026
21-function [Un, Vn, u, v] = step(U, V, lam, filt, gx, gy, gz, Vtx, Vty, Vtz, Vpx, Vpy, Vpz, a, b, D1, D2, dt, niter)
27+function [Un, Vn, u, v] = step(U, V, lam, filt, gx, gy, gz, p1, p2, q2, r, a, b, D1, D2, dt, niter)
2228 u = synth(U);
2329 v = synth(V);
2430 uuv = u .* u .* v;
@@ -32,56 +38,38 @@ function [Un, Vn, u, v] = step(U, V, lam, filt, gx, gy, gz, Vtx, Vty, Vtz, Vpx,
3238 Vn = Bv ./ (1 + (dt * D2) * lam);
3339
3440 for k = 1:niter
35- % dlap = lap_g - lap_s, evaluated at the current iterate (Algorithm 3 of
36- % evolving_surface/notes/algos.tex): surface gradient of the field,
37- % contracted through the inverse metric quantities Vt*/Vp*; each
38- % Cartesian component re-analysed and differentiated again; recombined
39- % into the surface divergence. lam.*Un adds back -lap_s(Un), since lam
40- % holds +l(l+1). filt zeroes the top two degrees, where the theta/phi
41- % derivative recurrences cannot exactly represent a derivative. See
42- % docs/richardson-iteration.md.
41+ % dlap = lap_g - lap_s at the current iterate, in flux form
42+ % (docs/reduced-transforms.md Sec 4). The sin-weighted
43+ % derivatives sin(theta)*dtheta(u) and dphi(u) -- both smooth on the
44+ % sphere, synthesized straight from the dthetac/dphic coefficient
45+ % shuffles -- are combined pointwise through the precomputed weights
46+ % p1,p2,q2 into two fluxes P,Q, also smooth. Their coefficients are then
47+ % pushed through the *same* shuffles again and summed before the one
48+ % synthesis of the divergence, which r scales into lap_g(u). The only
49+ % division by sin(theta) anywhere is folded into p1,p2,q2,r at precompute
50+ % time. lam.*Un adds back -lap_s(Un), since lam holds +l(l+1). filt
51+ % zeroes the top two degrees, where the derivative recurrences cannot
52+ % exactly represent a derivative.
4353 Fu = Un .* filt;
44- Ftu = dtheta(Fu);
45- Fpu = dphi(Fu);
46- dux = Ftu .* Vtx + Fpu .* Vpx;
47- duy = Ftu .* Vty + Fpu .* Vpy;
48- duz = Ftu .* Vtz + Fpu .* Vpz;
49- cux = analys(dux) .* filt;
50- cuy = analys(duy) .* filt;
51- cuz = analys(duz) .* filt;
52- Ftcux = dtheta(cux);
53- Fpcux = dphi(cux);
54- Ftcuy = dtheta(cuy);
55- Fpcuy = dphi(cuy);
56- Ftcuz = dtheta(cuz);
57- Fpcuz = dphi(cuz);
58- lapu = Ftcux .* Vtx + Fpcux .* Vpx;
59- lapu = lapu + Ftcuy .* Vty;
60- lapu = lapu + Fpcuy .* Vpy;
61- lapu = lapu + Ftcuz .* Vtz;
62- lapu = lapu + Fpcuz .* Vpz;
54+ Ftu = synth(dthetac(Fu));
55+ Fpu = synth(dphic(Fu));
56+ Pu = p1 .* Ftu + p2 .* Fpu;
57+ Qu = p2 .* Ftu + q2 .* Fpu;
58+ Pcu = analys(Pu) .* filt;
59+ Qcu = analys(Qu) .* filt;
60+ scu = dthetac(Pcu) + dphic(Qcu);
61+ lapu = r .* synth(scu);
6362 dLu = analys(lapu) + lam .* Un;
6463
6564 Fv = Vn .* filt;
66- Ftv = dtheta(Fv);
67- Fpv = dphi(Fv);
68- dvx = Ftv .* Vtx + Fpv .* Vpx;
69- dvy = Ftv .* Vty + Fpv .* Vpy;
70- dvz = Ftv .* Vtz + Fpv .* Vpz;
71- cvx = analys(dvx) .* filt;
72- cvy = analys(dvy) .* filt;
73- cvz = analys(dvz) .* filt;
74- Ftcvx = dtheta(cvx);
75- Fpcvx = dphi(cvx);
76- Ftcvy = dtheta(cvy);
77- Fpcvy = dphi(cvy);
78- Ftcvz = dtheta(cvz);
79- Fpcvz = dphi(cvz);
80- lapv = Ftcvx .* Vtx + Fpcvx .* Vpx;
81- lapv = lapv + Ftcvy .* Vty;
82- lapv = lapv + Fpcvy .* Vpy;
83- lapv = lapv + Ftcvz .* Vtz;
84- lapv = lapv + Fpcvz .* Vpz;
65+ Ftv = synth(dthetac(Fv));
66+ Fpv = synth(dphic(Fv));
67+ Pv = p1 .* Ftv + p2 .* Fpv;
68+ Qv = p2 .* Ftv + q2 .* Fpv;
69+ Pcv = analys(Pv) .* filt;
70+ Qcv = analys(Qv) .* filt;
71+ scv = dthetac(Pcv) + dphic(Qcv);
72+ lapv = r .* synth(scv);
8573 dLv = analys(lapv) + lam .* Vn;
8674
8775 Un = (Bu + (dt * D1) * dLu) ./ (1 + (dt * D1) * lam);
models/schnakenberg_alg4.madded+98−0View file
@@ -0,0 +1,98 @@
1+% Schnakenberg reaction-diffusion on a closed surface.
2+%
3+% du/dt = D1*lap_g(u) + a - u + u^2*v
4+% dv/dt = D2*lap_g(v) + b - u^2*v
5+%
6+% Explicit reaction, implicit diffusion (IMEX Euler). The implicit solve
7+% splits lap_g = lap_s + dlap: the round-sphere part lap_s is diagonal in
8+% spherical-harmonic space (eigenvalues -lam), and the loop iterates the
9+% geometric correction dlap from that exact solve. Grid fields are npts x 1;
10+% spectral fields are real 2 x nlm. See docs/richardson-iteration.md.
11+%
12+% This is the 12-transform reference: the correction evaluates lap_g in
13+% Cartesian-gradient form (Algorithm 4 of evolving_surface/notes/algos.tex),
14+% carrying grad_g(u) as three ambient components through the inverse metric
15+% quantities Vt*/Vp*. models/schnakenberg.m computes the same operator in
16+% flux form with 6 transforms per species per iteration
17+% (docs/reduced-transforms.md); this variant is kept live
18+% for A/B comparison, in the app and in the tests.
19+
20+function [U, V, u, v] = init(noise, a, b)
21+ us = a + b;
22+ vs = b / (us * us);
23+ U = analys(us + noise);
24+ V = analys(vs * ones(numel(noise), 1));
25+ u = synth(U);
26+ v = synth(V);
27+end
28+
29+function [Un, Vn, u, v] = step(U, V, lam, filt, gx, gy, gz, Vtx, Vty, Vtz, Vpx, Vpy, Vpz, a, b, D1, D2, dt, niter)
30+ u = synth(U);
31+ v = synth(V);
32+ uuv = u .* u .* v;
33+
34+ % Right-hand side of the implicit solve (I - dt*D*lap_g) Unew = B.
35+ Bu = U + dt * analys(a - u + uuv);
36+ Bv = V + dt * analys(b - uuv);
37+
38+ % Round-sphere solve, then iterate the geometric correction.
39+ Un = Bu ./ (1 + (dt * D1) * lam);
40+ Vn = Bv ./ (1 + (dt * D2) * lam);
41+
42+ for k = 1:niter
43+ % dlap = lap_g - lap_s, evaluated at the current iterate (Algorithm 3 of
44+ % evolving_surface/notes/algos.tex): surface gradient of the field,
45+ % contracted through the inverse metric quantities Vt*/Vp*; each
46+ % Cartesian component re-analysed and differentiated again; recombined
47+ % into the surface divergence. lam.*Un adds back -lap_s(Un), since lam
48+ % holds +l(l+1). filt zeroes the top two degrees, where the theta/phi
49+ % derivative recurrences cannot exactly represent a derivative. See
50+ % docs/richardson-iteration.md.
51+ Fu = Un .* filt;
52+ Ftu = dtheta(Fu);
53+ Fpu = dphi(Fu);
54+ dux = Ftu .* Vtx + Fpu .* Vpx;
55+ duy = Ftu .* Vty + Fpu .* Vpy;
56+ duz = Ftu .* Vtz + Fpu .* Vpz;
57+ cux = analys(dux) .* filt;
58+ cuy = analys(duy) .* filt;
59+ cuz = analys(duz) .* filt;
60+ Ftcux = dtheta(cux);
61+ Fpcux = dphi(cux);
62+ Ftcuy = dtheta(cuy);
63+ Fpcuy = dphi(cuy);
64+ Ftcuz = dtheta(cuz);
65+ Fpcuz = dphi(cuz);
66+ lapu = Ftcux .* Vtx + Fpcux .* Vpx;
67+ lapu = lapu + Ftcuy .* Vty;
68+ lapu = lapu + Fpcuy .* Vpy;
69+ lapu = lapu + Ftcuz .* Vtz;
70+ lapu = lapu + Fpcuz .* Vpz;
71+ dLu = analys(lapu) + lam .* Un;
72+
73+ Fv = Vn .* filt;
74+ Ftv = dtheta(Fv);
75+ Fpv = dphi(Fv);
76+ dvx = Ftv .* Vtx + Fpv .* Vpx;
77+ dvy = Ftv .* Vty + Fpv .* Vpy;
78+ dvz = Ftv .* Vtz + Fpv .* Vpz;
79+ cvx = analys(dvx) .* filt;
80+ cvy = analys(dvy) .* filt;
81+ cvz = analys(dvz) .* filt;
82+ Ftcvx = dtheta(cvx);
83+ Fpcvx = dphi(cvx);
84+ Ftcvy = dtheta(cvy);
85+ Fpcvy = dphi(cvy);
86+ Ftcvz = dtheta(cvz);
87+ Fpcvz = dphi(cvz);
88+ lapv = Ftcvx .* Vtx + Fpcvx .* Vpx;
89+ lapv = lapv + Ftcvy .* Vty;
90+ lapv = lapv + Fpcvy .* Vpy;
91+ lapv = lapv + Ftcvz .* Vtz;
92+ lapv = lapv + Fpcvz .* Vpz;
93+ dLv = analys(lapv) + lam .* Vn;
94+
95+ Un = (Bu + (dt * D1) * dLu) ./ (1 + (dt * D1) * lam);
96+ Vn = (Bv + (dt * D2) * dLv) ./ (1 + (dt * D2) * lam);
97+ end
98+end
scripts/sht-accuracy.tsadded+83−0View file
@@ -0,0 +1,83 @@
1+/**
2+ * The fp32 transform round-trip floor, swept in lmax — the measurement
3+ * docs/reduced-transforms.md Sec 5b prescribes before and
4+ * after any change to summation order in the Legendre kernels.
5+ *
6+ * npx vite-node scripts/sht-accuracy.ts [--lmax 63,127,255] [--seed 42]
7+ *
8+ * For a band-limited spectrum q, analys(synth(q)) = q exactly (Gauss
9+ * quadrature is exact for the band), so the relative round-trip error is the
10+ * transforms' own fp32 round-off with no reference implementation in the
11+ * loop. Sequential accumulation over l in the synthesis shows this floor
12+ * growing roughly linearly in lmax; pairwise/compensated accumulation shows
13+ * it near-flat. The per-degree profile says *where* the error lives (the
14+ * high-l coefficients are the ones the alpha shifts and the l(l+1)
15+ * eigenvalues amplify).
16+ */
17+import { ShtPlan, requestShtDevice, describeAdapter } from '../src/sht/sht.ts';
18+import { gridForLmax, lmIndex, nlmCalc } from '../src/sht/layout.ts';
19+import { randomSpectrum } from '../src/sht/reference.ts';
20+import { installWebGpu, errMsg, NO_ADAPTER_HINT } from './nodeWebGpu.ts';
21+
22+const arg = (name: string): string | undefined => {
23+ const i = process.argv.indexOf(`--${name}`);
24+ return i >= 0 ? process.argv[i + 1] : undefined;
25+};
26+const LMAXES = (arg('lmax') ?? '63,127,255').split(',').map(Number);
27+const SEED = Number(arg('seed') ?? 42);
28+
29+let runtime: string;
30+try {
31+ runtime = await installWebGpu();
32+} catch (e) {
33+ console.error(`sht-accuracy: ${errMsg(e)}\n${NO_ADAPTER_HINT}`);
34+ process.exit(1);
35+}
36+const device = await requestShtDevice();
37+console.log(`sht-accuracy — ${runtime}, ${await describeAdapter(device)}\n`);
38+console.log(' lmax grid rel L2 roundtrip worst degree (rel)');
39+
40+for (const lmax of LMAXES) {
41+ const { nlat, nphi } = gridForLmax(lmax, 1);
42+ const cfg = { lmax, mmax: lmax, nlat, nphi };
43+ const plan = await ShtPlan.create(device, cfg);
44+ const q = randomSpectrum(cfg, SEED);
45+
46+ const grid = await plan.synth(q);
47+ const back = await plan.analys(grid);
48+
49+ // Overall relative L2, and the same per degree — errors concentrate in l.
50+ let num = 0;
51+ let den = 0;
52+ const nlm = nlmCalc(lmax, lmax);
53+ const errL = new Float64Array(lmax + 1);
54+ const magL = new Float64Array(lmax + 1);
55+ for (let m = 0; m <= lmax; m++) {
56+ for (let l = m; l <= lmax; l++) {
57+ const i = lmIndex(lmax, l, m);
58+ const dr = back[2 * i] - q[2 * i];
59+ const di = back[2 * i + 1] - q[2 * i + 1];
60+ const d2 = dr * dr + di * di;
61+ const m2 = q[2 * i] ** 2 + q[2 * i + 1] ** 2;
62+ num += d2;
63+ den += m2;
64+ errL[l] += d2;
65+ magL[l] += m2;
66+ }
67+ }
68+ let worstL = 0;
69+ let worstRel = 0;
70+ for (let l = 0; l <= lmax; l++) {
71+ const rel = Math.sqrt(errL[l] / Math.max(magL[l], 1e-300));
72+ if (rel > worstRel) {
73+ worstRel = rel;
74+ worstL = l;
75+ }
76+ }
77+ console.log(
78+ ` ${String(lmax).padEnd(6)} ${`${nlat}x${nphi}`.padEnd(11)} ` +
79+ `${Math.sqrt(num / den).toExponential(3).padEnd(18)} ` +
80+ `l=${worstL}: ${worstRel.toExponential(3)} (${nlm} coefficients)`,
81+ );
82+ plan.destroy();
83+}
scripts/test-node.tsmodified+2−0View file
@@ -14,6 +14,7 @@ import { transformChecks } from '../test/transformChecks.ts';
1414 import { analyticChecks } from '../test/analyticChecks.ts';
1515 import { modelChecks } from '../test/modelChecks.ts';
1616 import { geometryChecks } from '../test/geometryChecks.ts';
17+import { fluxChecks } from '../test/fluxChecks.ts';
1718
1819 let failures = 0;
1920 const check = (name: string, ok: boolean, detail: string): void => {
@@ -51,6 +52,7 @@ await transformChecks(device, check, log);
5152 await analyticChecks(device, check, log);
5253 await modelChecks(device, check, log);
5354 await geometryChecks(device, check, log);
55+await fluxChecks(device, check, log);
5456
5557 console.log(failures === 0 ? '\nAll tests passed.' : `\n${failures} failed.`);
5658 process.exit(failures === 0 ? 0 : 1);
src/geom/geometry.tsmodified+34−3View file
@@ -31,7 +31,7 @@
3131 import { ShtPlan } from '../sht/sht.ts';
3232 import type { ShtConfig } from '../sht/layout.ts';
3333 import type { DerivPlan } from '../sht/deriv.ts';
34-import { computeMetric } from './metric.ts';
34+import { computeMetric, computeFluxMetric } from './metric.ts';
3535 import { HostBuffers, ModelPlan } from '../mgpu/plan.ts';
3636 import { CompiledModel, type Binding } from '../mgpu/compile.ts';
3737 import { inFunction, inFunctionAsync, inModel } from '../mgpu/errors.ts';
@@ -66,7 +66,8 @@ export class Geometry {
6666 /**
6767 * Inverse metric quantities (src/geom/metric.ts), grid space, npts each.
6868 * Depend only on the geometry, so — like x,y,z,X,Y,Z above — these are a
69- * one-off computed here, not per-solve-step work.
69+ * one-off computed here, not per-solve-step work. Used by the Algorithm-4
70+ * (12-transform) Laplace-Beltrami path.
7071 */
7172 readonly Vtx: Float32Array;
7273 readonly Vty: Float32Array;
@@ -74,12 +75,23 @@ export class Geometry {
7475 readonly Vpx: Float32Array;
7576 readonly Vpy: Float32Array;
7677 readonly Vpz: Float32Array;
78+ /**
79+ * Flux-form metric weights (src/geom/metric.ts computeFluxMetric), grid
80+ * space, npts each — the six-transform Laplace-Beltrami scheme's
81+ * replacement for the six V arrays (docs/reduced-transforms.md
82+ * Sec 3). Both sets are carried so either operator formulation can run.
83+ */
84+ readonly p1: Float32Array;
85+ readonly p2: Float32Array;
86+ readonly q2: Float32Array;
87+ readonly r: Float32Array;
7788
7889 private constructor(init: {
7990 x: Float32Array; y: Float32Array; z: Float32Array;
8091 X: Float32Array; Y: Float32Array; Z: Float32Array;
8192 Vtx: Float32Array; Vty: Float32Array; Vtz: Float32Array;
8293 Vpx: Float32Array; Vpy: Float32Array; Vpz: Float32Array;
94+ p1: Float32Array; p2: Float32Array; q2: Float32Array; r: Float32Array;
8395 }) {
8496 this.x = init.x;
8597 this.y = init.y;
@@ -93,6 +105,10 @@ export class Geometry {
93105 this.Vpx = init.Vpx;
94106 this.Vpy = init.Vpy;
95107 this.Vpz = init.Vpz;
108+ this.p1 = init.p1;
109+ this.p2 = init.p2;
110+ this.q2 = init.q2;
111+ this.r = init.r;
96112 }
97113
98114 /**
@@ -164,7 +180,22 @@ export class Geometry {
164180 const Zp = await deriv.dphi(Z);
165181 const { Vtx, Vty, Vtz, Vpx, Vpy, Vpz } = computeMetric(npts, Xt, Xp, Yt, Yp, Zt, Zp);
166182
167- return new Geometry({ x, y, z, X, Y, Z, Vtx, Vty, Vtz, Vpx, Vpy, Vpz });
183+ // Flux-form metric weights for the six-transform scheme, built from the
184+ // *undivided* theta tangents sin(theta)*X_theta (smooth on the sphere,
185+ // unlike X_theta itself) and the same X_phi as above. Also a one-off;
186+ // the f64 combination happens on the CPU, rounded to f32 for upload.
187+ const sXtx = await deriv.sinDtheta(X);
188+ const sXty = await deriv.sinDtheta(Y);
189+ const sXtz = await deriv.sinDtheta(Z);
190+ const flux = computeFluxMetric(npts, sXtx, sXty, sXtz, Xp, Yp, Zp);
191+
192+ return new Geometry({
193+ x, y, z, X, Y, Z, Vtx, Vty, Vtz, Vpx, Vpy, Vpz,
194+ p1: new Float32Array(flux.p1),
195+ p2: new Float32Array(flux.p2),
196+ q2: new Float32Array(flux.q2),
197+ r: new Float32Array(flux.r),
198+ });
168199 } finally {
169200 plan.destroy();
170201 host.destroy();
src/geom/metric.tsmodified+77−0View file
@@ -67,3 +67,80 @@ export function computeMetric(
6767
6868 return { Vtx, Vty, Vtz, Vpx, Vpy, Vpz };
6969 }
70+
71+/**
72+ * Flux-form metric weights p1, p2, q2, r of the six-transform Laplace-Beltrami
73+ * scheme (docs/reduced-transforms.md Sec 3). Built from the
74+ * *sin-weighted* theta tangent sin(theta)*X_theta — the undivided synthesis of
75+ * the alpha shift, DerivPlan.sinDtheta — and X_phi, both smooth on the sphere:
76+ *
77+ * gtt~ = |sin(theta) X_theta|^2 (= sin^2(theta) g_tt)
78+ * gtp~ = (sin(theta) X_theta).X_phi (= sin(theta) g_tp)
79+ * gpp = |X_phi|^2
80+ * D = sqrt(gtt~ gpp - gtp~^2) (= sin^2(theta) sqrt(det g) / sin(theta)
81+ * = J sin^2(theta), with J = sqrt(det g)/sin(theta))
82+ *
83+ * p1 = gpp / D, p2 = -gtp~ / D, q2 = gtt~ / D, r = 1 / D.
84+ *
85+ * With these, for A = sin(theta) dtheta(u) and B = dphi(u), the two fluxes
86+ *
87+ * P = p1*A + p2*B, Qtilde = p2*A + q2*B
88+ *
89+ * equal sqrt(det g) g^{theta j} u_j and sin(theta) sqrt(det g) g^{phi j} u_j —
90+ * both smooth on the sphere — and Delta_Gamma u = r * (sin(theta) dtheta(P) +
91+ * dphi(Qtilde)). p1, p2, q2 are bounded (the sin^2 in D cancels against the
92+ * vanishing numerators); r ~ 1/sin^2(theta) is finite at the Gauss nodes and
93+ * is the scheme's one concentrated division (Sec 5 of the doc).
94+ *
95+ * All arithmetic is f64 (JS numbers) regardless of the input arrays' storage
96+ * type; results are rounded to f32 only on upload. That is the doc's "CPU
97+ * precompute in float64" mitigation, inherited for free.
98+ */
99+export interface FluxMetricFields {
100+ p1: Float64Array;
101+ p2: Float64Array;
102+ q2: Float64Array;
103+ r: Float64Array;
104+}
105+
106+/**
107+ * sXt* are the Cartesian components of sin(theta)*X_theta, Xp* those of
108+ * X_phi, all grid space, npts each. No sin(theta) input is needed: every
109+ * division the scheme performs is by D, which the sin-weighted inputs build
110+ * directly.
111+ */
112+export function computeFluxMetric(
113+ npts: number,
114+ sXtx: ArrayLike<number>,
115+ sXty: ArrayLike<number>,
116+ sXtz: ArrayLike<number>,
117+ Xpx: ArrayLike<number>,
118+ Xpy: ArrayLike<number>,
119+ Xpz: ArrayLike<number>,
120+): FluxMetricFields {
121+ const p1 = new Float64Array(npts);
122+ const p2 = new Float64Array(npts);
123+ const q2 = new Float64Array(npts);
124+ const r = new Float64Array(npts);
125+
126+ for (let i = 0; i < npts; i++) {
127+ const xt = sXtx[i];
128+ const yt = sXty[i];
129+ const zt = sXtz[i];
130+ const xp = Xpx[i];
131+ const yp = Xpy[i];
132+ const zp = Xpz[i];
133+
134+ const gtt = xt * xt + yt * yt + zt * zt; // sin^2 g_tt
135+ const gtp = xt * xp + yt * yp + zt * zp; // sin g_tp
136+ const gpp = xp * xp + yp * yp + zp * zp; // g_pp
137+ const D = Math.sqrt(gtt * gpp - gtp * gtp); // J sin^2(theta)
138+
139+ p1[i] = gpp / D;
140+ p2[i] = -gtp / D;
141+ q2[i] = gtt / D;
142+ r[i] = 1 / D;
143+ }
144+
145+ return { p1, p2, q2, r };
146+}
src/mgpu/externals.tsmodified+21−6View file
@@ -78,11 +78,16 @@ exports.cBody = function () {
7878 }
7979
8080 /**
81- * Workspace files that make `synth` / `analys` / `dtheta` / `dphi` resolvable
82- * during lowering. `dtheta` and `dphi` (the surface's first partial
83- * derivatives, coefficients -> grid — see src/sht/deriv.ts) have exactly
84- * `synth`'s shape rule: both take spectral coefficients and produce a grid
85- * field.
81+ * Workspace files that make `synth` / `analys` / `dtheta` / `dphi` /
82+ * `dthetac` / `dphic` resolvable during lowering. `dtheta` and `dphi` (the
83+ * surface's first partial derivatives, coefficients -> grid — see
84+ * src/sht/deriv.ts) have exactly `synth`'s shape rule: both take spectral
85+ * coefficients and produce a grid field. `dthetac` and `dphic` are their
86+ * coefficient-space halves alone — the alpha^+/alpha^- shift and the i*m
87+ * multiply, spectral -> spectral — which the six-transform Laplace-Beltrami
88+ * scheme (docs/reduced-transforms.md) applies twice per
89+ * matvec: to the field (gradient side) and to the analysed fluxes
90+ * (divergence side, the same shift, not its transpose).
8691 */
8792 export function externalOpFiles(g: GridSizes): { name: string; source: string }[] {
8893 return [
@@ -102,8 +107,18 @@ export function externalOpFiles(g: GridSizes): { name: string; source: string }[
102107 name: 'dphi.mtoc2.js',
103108 source: transformSource('dphi', 2, g.nlm, g.npts, 1),
104109 },
110+ {
111+ name: 'dthetac.mtoc2.js',
112+ source: transformSource('dthetac', 2, g.nlm, 2, g.nlm),
113+ },
114+ {
115+ name: 'dphic.mtoc2.js',
116+ source: transformSource('dphic', 2, g.nlm, 2, g.nlm),
117+ },
105118 ];
106119 }
107120
108121 /** Names the WGSL backend must implement as GPU encodes rather than kernels. */
109-export const EXTERNAL_OPS = new Set(['synth', 'analys', 'dtheta', 'dphi']);
122+export const EXTERNAL_OPS = new Set([
123+ 'synth', 'analys', 'dtheta', 'dphi', 'dthetac', 'dphic',
124+]);
src/mgpu/model.tsmodified+21−2View file
@@ -68,20 +68,31 @@ export interface GeometryBuffers {
6868 X: Float32Array;
6969 Y: Float32Array;
7070 Z: Float32Array;
71- /** Inverse metric quantities (src/geom/metric.ts), grid space, npts each. */
71+ /** Inverse metric quantities (src/geom/metric.ts), grid space, npts each —
72+ * the Algorithm-4 (12-transform) Laplace-Beltrami path. */
7273 Vtx: Float32Array;
7374 Vty: Float32Array;
7475 Vtz: Float32Array;
7576 Vpx: Float32Array;
7677 Vpy: Float32Array;
7778 Vpz: Float32Array;
79+ /** Flux-form metric weights (src/geom/metric.ts computeFluxMetric), grid
80+ * space, npts each — the six-transform Laplace-Beltrami scheme of
81+ * docs/reduced-transforms.md. */
82+ p1: Float32Array;
83+ p2: Float32Array;
84+ q2: Float32Array;
85+ r: Float32Array;
7886 }
7987
8088 /** Names the .m may take for the grid coordinates and for their coefficients. */
8189 export const GEOMETRY_GRID_NAMES = ['gx', 'gy', 'gz'] as const;
8290 export const GEOMETRY_SPECTRAL_NAMES = ['Gx', 'Gy', 'Gz'] as const;
83-/** Names the .m may take for the inverse metric quantities. */
91+/** Names the .m may take for the inverse metric quantities (Algorithm 4). */
8492 export const METRIC_GRID_NAMES = ['Vtx', 'Vty', 'Vtz', 'Vpx', 'Vpy', 'Vpz'] as const;
93+/** Names the .m may take for the flux-form metric weights (six-transform
94+ * scheme). A model asks for whichever set its loop uses; both are uploaded. */
95+export const FLUX_METRIC_GRID_NAMES = ['p1', 'p2', 'q2', 'r'] as const;
8596
8697 /** Laplace-Beltrami eigenvalues l(l+1), duplicated across re/im so the array
8798 * matches the 2 x nlm spectral layout element for element. */
@@ -184,6 +195,7 @@ export class GpuModel {
184195 for (const g of GEOMETRY_GRID_NAMES) bindings[g] = { kind: 'tensor', shape: [npts, 1] };
185196 for (const g of GEOMETRY_SPECTRAL_NAMES) bindings[g] = { kind: 'tensor', shape: [2, nlm] };
186197 for (const g of METRIC_GRID_NAMES) bindings[g] = { kind: 'tensor', shape: [npts, 1] };
198+ for (const g of FLUX_METRIC_GRID_NAMES) bindings[g] = { kind: 'tensor', shape: [npts, 1] };
187199 }
188200 for (const s of state) bindings[s] = { kind: 'tensor', shape: [2, nlm] };
189201 for (const p of paramNames) bindings[p] = { kind: 'param' };
@@ -212,6 +224,7 @@ export class GpuModel {
212224 for (const g of GEOMETRY_GRID_NAMES) host.ensure(g, npts);
213225 for (const g of GEOMETRY_SPECTRAL_NAMES) host.ensure(g, 2 * nlm);
214226 for (const g of METRIC_GRID_NAMES) host.ensure(g, npts);
227+ for (const g of FLUX_METRIC_GRID_NAMES) host.ensure(g, npts);
215228 }
216229
217230 const initPlan = await inFunctionAsync('init', () =>
@@ -236,6 +249,10 @@ export class GpuModel {
236249 host.upload('Vpx', geometry.Vpx);
237250 host.upload('Vpy', geometry.Vpy);
238251 host.upload('Vpz', geometry.Vpz);
252+ host.upload('p1', geometry.p1);
253+ host.upload('p2', geometry.p2);
254+ host.upload('q2', geometry.q2);
255+ host.upload('r', geometry.r);
239256 }
240257
241258 const readback = device.createBuffer({
@@ -281,6 +298,8 @@ export class GpuModel {
281298 ['Gx', geometry.X], ['Gy', geometry.Y], ['Gz', geometry.Z],
282299 ['Vtx', geometry.Vtx], ['Vty', geometry.Vty], ['Vtz', geometry.Vtz],
283300 ['Vpx', geometry.Vpx], ['Vpy', geometry.Vpy], ['Vpz', geometry.Vpz],
301+ ['p1', geometry.p1], ['p2', geometry.p2],
302+ ['q2', geometry.q2], ['r', geometry.r],
284303 ];
285304 for (const [name, data] of fields) {
286305 if (this.#host.get(name)) this.#host.upload(name, data);
src/mgpu/plan.tsmodified+34−1View file
@@ -121,6 +121,7 @@ type Op =
121121 }
122122 | { kind: 'synth' | 'analys'; binding: ShtBinding; label: string }
123123 | { kind: 'dtheta' | 'dphi'; binding: DerivBinding; label: string }
124+ | { kind: 'dthetac' | 'dphic'; bindGroup: GPUBindGroup; label: string }
124125 | { kind: 'copy'; from: GPUBuffer; to: GPUBuffer; bytes: number; label: string };
125126
126127 export interface PlanSpec {
@@ -368,7 +369,10 @@ export class ModelPlan {
368369 binding: sht.createAnalysBinding(argSlot.buffer, dest.buffer),
369370 label,
370371 });
371- } else if (ext.name === 'dtheta' || ext.name === 'dphi') {
372+ } else if (
373+ ext.name === 'dtheta' || ext.name === 'dphi' ||
374+ ext.name === 'dthetac' || ext.name === 'dphic'
375+ ) {
372376 if (!deriv) {
373377 throw new UnsupportedOnGpu(
374378 `'${ext.name}' needs the surface's derivative transforms, ` +
@@ -376,6 +380,29 @@ export class ModelPlan {
376380 stmt.span,
377381 );
378382 }
383+ if (ext.name === 'dthetac' || ext.name === 'dphic') {
384+ // Coefficient-space shuffles read at l+-1 (dthetac) or in place
385+ // (dphic) and cannot alias their output: WebGPU forbids one buffer
386+ // being readable and writable storage in the same dispatch, and
387+ // there is no scratch-copy fallback here — refuse rather than
388+ // silently reroute.
389+ if (argSlot.buffer === dest.buffer) {
390+ throw new UnsupportedOnGpu(
391+ `'${stmt.name} = ${ext.name}(${ext.argName})' reads and ` +
392+ `writes the same buffer; assign to a new name instead`,
393+ stmt.span,
394+ );
395+ }
396+ ops.push({
397+ kind: ext.name,
398+ bindGroup:
399+ ext.name === 'dthetac'
400+ ? deriv.createDthetacBinding(argSlot.buffer, dest.buffer)
401+ : deriv.createDphicBinding(argSlot.buffer, dest.buffer),
402+ label,
403+ });
404+ return;
405+ }
379406 ops.push(
380407 ext.name === 'dtheta'
381408 ? { kind: 'dtheta', binding: deriv.createDthetaBinding(argSlot.buffer, dest.buffer), label }
@@ -572,6 +599,12 @@ export class ModelPlan {
572599 case 'dphi':
573600 this.#derivInto(inPass(), op);
574601 break;
602+ case 'dthetac':
603+ this.#deriv!.encodeDthetacInto(inPass(), op.bindGroup);
604+ break;
605+ case 'dphic':
606+ this.#deriv!.encodeDphicInto(inPass(), op.bindGroup);
607+ break;
575608 case 'copy':
576609 endPass();
577610 encoder.copyBufferToBuffer(op.from, 0, op.to, 0, op.bytes);
src/mgpu/registry.tsmodified+23−1View file
@@ -12,6 +12,7 @@
1212 * `U`, `V`, ... the corresponding spectral state (uppercase)
1313 */
1414 import schnakenbergSource from '../../models/schnakenberg.m?raw';
15+import schnakenbergAlg4Source from '../../models/schnakenberg_alg4.m?raw';
1516 import brusselatorSource from '../../models/brusselator.m?raw';
1617 import allencahnSource from '../../models/allencahn.m?raw';
1718
@@ -65,6 +66,22 @@ const schnakenberg: MModel = {
6566 source: schnakenbergSource,
6667 };
6768
69+/**
70+ * The same PDE and parameters as `schnakenberg`, with the implicit solve's
71+ * geometric correction in its original Cartesian-gradient form (Algorithm 4,
72+ * 12 transforms per species per iteration) instead of the flux form's 6
73+ * (docs/reduced-transforms.md). Shipped as a live reference:
74+ * the two must agree to fp32 accuracy on any surface, and the tests hold
75+ * them to that.
76+ */
77+const schnakenbergAlg4: MModel = {
78+ ...schnakenberg,
79+ key: 'schnakenberg-alg4',
80+ label: 'Schnakenberg (12-transform reference)',
81+ blurb: 'Same spots, Algorithm-4 Laplace-Beltrami — for A/B against the flux form.',
82+ source: schnakenbergAlg4Source,
83+};
84+
6885 const brusselator: MModel = {
6986 key: 'brusselator',
7087 label: 'Brusselator',
@@ -98,7 +115,7 @@ const allencahn: MModel = {
98115 source: allencahnSource,
99116 };
100117
101-export const mModels: MModel[] = [schnakenberg, brusselator, allencahn];
118+export const mModels: MModel[] = [schnakenberg, brusselator, allencahn, schnakenbergAlg4];
102119
103120 export const mModelByKey = (key: string): MModel | undefined =>
104121 mModels.find((m) => m.key === key);
@@ -133,4 +150,9 @@ export const presets: Preset[] = [
133150 },
134151 { key: 'brussel', label: 'Brusselator — stripes & spots', modelKey: 'brusselator' },
135152 { key: 'allencahn', label: 'Allen–Cahn — coarsening', modelKey: 'allencahn' },
153+ {
154+ key: 'schnak-alg4',
155+ label: 'Schnakenberg — spots (12-transform reference)',
156+ modelKey: 'schnakenberg-alg4',
157+ },
136158 ];
src/mgpu/session.tsmodified+3−2View file
@@ -134,8 +134,9 @@ export class ModelSession {
134134 deriv = await DerivPlan.create(device, sht);
135135 // The surface is built before the model, because the model takes it as
136136 // an argument. It is a one-off: compiled, evaluated, read back, and its
137- // plan discarded — nothing of it survives into the timestep but twelve
138- // buffers of numbers (the embedding and the metric quantities built on it).
137+ // plan discarded — nothing of it survives into the timestep but sixteen
138+ // buffers of numbers (the embedding, and both metric formulations built
139+ // on it: the inverse metric quantities and the flux-form weights).
139140 const geometry = await Geometry.create({
140141 device,
141142 sht,
src/sht/deriv.tsmodified+82−26View file
@@ -10,6 +10,16 @@
1010 * pipeline (ShtPlan.createSynthBinding/encodeSynthInto) unchanged -- neither
1111 * derivative touches the Legendre recurrence stage itself. dtheta
1212 * additionally divides by sin(theta) on the grid afterwards.
13+ *
14+ * The two shuffles are also exposed on their own, coefficients -> coefficients
15+ * (`dthetac`, `dphic`), because the six-transform Laplace-Beltrami operator of
16+ * docs/reduced-transforms.md needs them apart from a
17+ * synthesis, and needs them twice: once on the field (steps 1-2) and once on
18+ * the two fluxes (step 5, which is the *same* alpha^+/alpha^- gather, not its
19+ * transpose). Everything above is then a composition of them:
20+ *
21+ * dphi(U) == synth(dphic(U))
22+ * dtheta(U) == synth(dthetac(U)) / sin(theta)
1323 */
1424 import type { ShtPlan, ShtBinding } from './sht.ts';
1525 import { derivCoeffs } from './derivCoeffs.ts';
@@ -120,17 +130,55 @@ export class DerivPlan {
120130 this.pipeDivide = pDivide;
121131 }
122132
123- /** Bindings for dtheta(qlmIn) -> spatOut, against caller-owned buffers. */
124- createDthetaBinding(qlmIn: GPUBuffer, spatOut: GPUBuffer): DerivBinding {
125- const shuffle = this.device.createBindGroup({
133+ /**
134+ * Bind group for the coefficient-space half of dtheta on its own:
135+ * v_l^m = alpha^+(l-1,m) u_{l-1}^m + alpha^-(l+1,m) u_{l+1}^m, the
136+ * coefficients of sin(theta) * dtheta(u). Input and output must be
137+ * different buffers -- WebGPU forbids binding one buffer as both readable
138+ * and writable storage in a dispatch, and the gather reads l+-1 anyway.
139+ */
140+ createDthetacBinding(qlmIn: GPUBuffer, qlmOut: GPUBuffer): GPUBindGroup {
141+ return this.device.createBindGroup({
126142 layout: this.pipeDtheta.getBindGroupLayout(0),
127143 entries: [
128144 { binding: 0, resource: { buffer: this.bufAPlus } },
129145 { binding: 1, resource: { buffer: this.bufAMinus } },
130146 { binding: 2, resource: { buffer: qlmIn } },
131- { binding: 3, resource: { buffer: this.scratch } },
147+ { binding: 3, resource: { buffer: qlmOut } },
148+ ],
149+ });
150+ }
151+
152+ /** Bind group for the coefficient-space half of dphi on its own:
153+ * (dphi u)_l^m = i*m*u_l^m. Same buffer restriction as dthetac. */
154+ createDphicBinding(qlmIn: GPUBuffer, qlmOut: GPUBuffer): GPUBindGroup {
155+ return this.device.createBindGroup({
156+ layout: this.pipeDphi.getBindGroupLayout(0),
157+ entries: [
158+ { binding: 0, resource: { buffer: this.bufMOf } },
159+ { binding: 1, resource: { buffer: qlmIn } },
160+ { binding: 2, resource: { buffer: qlmOut } },
132161 ],
133162 });
163+ }
164+
165+ /** Record the bare alpha^+/alpha^- shift into an existing compute pass. */
166+ encodeDthetacInto(pass: GPUComputePassEncoder, bindGroup: GPUBindGroup): void {
167+ pass.setPipeline(this.pipeDtheta);
168+ pass.setBindGroup(0, bindGroup);
169+ pass.dispatchWorkgroups(Math.ceil(this.nlm / WG));
170+ }
171+
172+ /** Record the bare i*m multiply into an existing compute pass. */
173+ encodeDphicInto(pass: GPUComputePassEncoder, bindGroup: GPUBindGroup): void {
174+ pass.setPipeline(this.pipeDphi);
175+ pass.setBindGroup(0, bindGroup);
176+ pass.dispatchWorkgroups(Math.ceil(this.nlm / WG));
177+ }
178+
179+ /** Bindings for dtheta(qlmIn) -> spatOut, against caller-owned buffers. */
180+ createDthetaBinding(qlmIn: GPUBuffer, spatOut: GPUBuffer): DerivBinding {
181+ const shuffle = this.createDthetacBinding(qlmIn, this.scratch);
134182 const sht = this.sht.createSynthBinding(this.scratch, spatOut);
135183 const divide = this.device.createBindGroup({
136184 layout: this.pipeDivide.getBindGroupLayout(0),
@@ -144,48 +192,54 @@ export class DerivPlan {
144192
145193 /** Bindings for dphi(qlmIn) -> spatOut, against caller-owned buffers. */
146194 createDphiBinding(qlmIn: GPUBuffer, spatOut: GPUBuffer): DerivBinding {
147- const shuffle = this.device.createBindGroup({
148- layout: this.pipeDphi.getBindGroupLayout(0),
149- entries: [
150- { binding: 0, resource: { buffer: this.bufMOf } },
151- { binding: 1, resource: { buffer: qlmIn } },
152- { binding: 2, resource: { buffer: this.scratch } },
153- ],
154- });
195+ const shuffle = this.createDphicBinding(qlmIn, this.scratch);
155196 const sht = this.sht.createSynthBinding(this.scratch, spatOut);
156197 return { shuffle, sht };
157198 }
158199
159200 /** Record dtheta into an existing compute pass. */
160201 encodeDthetaInto(pass: GPUComputePassEncoder, b: DerivBinding): void {
161- pass.setPipeline(this.pipeDtheta);
162- pass.setBindGroup(0, b.shuffle);
163- pass.dispatchWorkgroups(Math.ceil(this.nlm / WG));
164- this.sht.encodeSynthInto(pass, b.sht);
202+ this.encodeSinDthetaInto(pass, b);
165203 pass.setPipeline(this.pipeDivide);
166204 pass.setBindGroup(0, b.divide!);
167205 pass.dispatchWorkgroups(Math.ceil(this.npts / WG));
168206 }
169207
208+ /** Record dtheta *without* its final division: the grid values of
209+ * sin(theta) * dtheta(u), which unlike dtheta(u) itself is a smooth
210+ * function on the sphere. Takes a dtheta binding and simply stops early. */
211+ encodeSinDthetaInto(pass: GPUComputePassEncoder, b: DerivBinding): void {
212+ this.encodeDthetacInto(pass, b.shuffle);
213+ this.sht.encodeSynthInto(pass, b.sht);
214+ }
215+
170216 /** Record dphi into an existing compute pass. */
171217 encodeDphiInto(pass: GPUComputePassEncoder, b: DerivBinding): void {
172- pass.setPipeline(this.pipeDphi);
173- pass.setBindGroup(0, b.shuffle);
174- pass.dispatchWorkgroups(Math.ceil(this.nlm / WG));
218+ this.encodeDphicInto(pass, b.shuffle);
175219 this.sht.encodeSynthInto(pass, b.sht);
176220 }
177221
178222 /** CPU convenience: qlm (interleaved [re,im], length 2*nlm) -> grid field. */
179223 async dtheta(qlm: Float32Array): Promise<Float32Array> {
180- return this.#runToGrid(qlm, true);
224+ return this.#runToGrid(qlm, 'dtheta');
225+ }
226+
227+ /** CPU convenience: sin(theta) * dtheta(u) on the grid, the undivided
228+ * synthesis of the alpha shift. What the flux-form metric precompute
229+ * (src/geom/metric.ts) is built from. */
230+ async sinDtheta(qlm: Float32Array): Promise<Float32Array> {
231+ return this.#runToGrid(qlm, 'sinDtheta');
181232 }
182233
183234 /** CPU convenience: qlm (interleaved [re,im], length 2*nlm) -> grid field. */
184235 async dphi(qlm: Float32Array): Promise<Float32Array> {
185- return this.#runToGrid(qlm, false);
236+ return this.#runToGrid(qlm, 'dphi');
186237 }
187238
188- async #runToGrid(qlm: Float32Array, withDivide: boolean): Promise<Float32Array> {
239+ async #runToGrid(
240+ qlm: Float32Array,
241+ mode: 'dtheta' | 'sinDtheta' | 'dphi',
242+ ): Promise<Float32Array> {
189243 if (qlm.length !== 2 * this.nlm) throw new Error(`qlm must have length ${2 * this.nlm}`);
190244 const dev = this.device;
191245 const qlmIn = dev.createBuffer({
@@ -205,12 +259,14 @@ export class DerivPlan {
205259 });
206260 try {
207261 dev.queue.writeBuffer(qlmIn, 0, qlm as Float32Array<ArrayBuffer>);
208- const binding = withDivide
209- ? this.createDthetaBinding(qlmIn, spatOut)
210- : this.createDphiBinding(qlmIn, spatOut);
262+ const binding =
263+ mode === 'dphi'
264+ ? this.createDphiBinding(qlmIn, spatOut)
265+ : this.createDthetaBinding(qlmIn, spatOut);
211266 const enc = dev.createCommandEncoder({ label: 'deriv-run' });
212267 const pass = enc.beginComputePass({ label: 'deriv-run' });
213- if (withDivide) this.encodeDthetaInto(pass, binding);
268+ if (mode === 'dtheta') this.encodeDthetaInto(pass, binding);
269+ else if (mode === 'sinDtheta') this.encodeSinDthetaInto(pass, binding);
214270 else this.encodeDphiInto(pass, binding);
215271 pass.end();
216272 enc.copyBufferToBuffer(spatOut, 0, stage, 0, 4 * this.npts);
src/sht/reference.tsmodified+34−13View file
@@ -119,11 +119,15 @@ export class ShtReference {
119119 }
120120
121121 /**
122- * Theta-derivative, f64: v_l^m = alpha^+(l-1,m) u_{l-1}^m + alpha^-(l+1,m)
123- * u_{l+1}^m (algos.tex eq. v_coeffs), then synth(v_l^m) / sin(theta).
122+ * Coefficient-space theta shift, f64: v_l^m = alpha^+(l-1,m) u_{l-1}^m +
123+ * alpha^-(l+1,m) u_{l+1}^m (algos.tex eq. v_coeffs) — the coefficients of
124+ * sin(theta) * dtheta(u). The same shift serves the divergence side of the
125+ * six-transform Laplace-Beltrami scheme (step 5 of
126+ * docs/reduced-transforms.md), which is why it is exposed
127+ * apart from the synthesis.
124128 */
125- dtheta(qlm: ArrayLike<number>): Float64Array {
126- const { lmax, mmax, nlat, nphi } = this.cfg;
129+ dthetac(qlm: ArrayLike<number>): Float64Array {
130+ const { lmax, mmax } = this.cfg;
127131 const v = new Float64Array(2 * this.nlm);
128132 for (let m = 0; m <= mmax; m++) {
129133 for (let l = m; l <= lmax; l++) {
@@ -146,16 +150,11 @@ export class ShtReference {
146150 v[2 * lm + 1] = im;
147151 }
148152 }
149- const grid = this.synth(v);
150- for (let i = 0; i < nlat; i++) {
151- const st = this.st[i];
152- for (let j = 0; j < nphi; j++) grid[i * nphi + j] /= st;
153- }
154- return grid;
153+ return v;
155154 }
156155
157- /** Phi-derivative, f64: (dphi u)_l^m = i*m*u_l^m, then synthesize. */
158- dphi(qlm: ArrayLike<number>): Float64Array {
156+ /** Coefficient-space phi derivative, f64: (dphi u)_l^m = i*m*u_l^m. */
157+ dphic(qlm: ArrayLike<number>): Float64Array {
159158 const { lmax, mmax } = this.cfg;
160159 const v = new Float64Array(2 * this.nlm);
161160 for (let m = 0; m <= mmax; m++) {
@@ -165,7 +164,29 @@ export class ShtReference {
165164 v[2 * lm + 1] = m * qlm[2 * lm];
166165 }
167166 }
168- return this.synth(v);
167+ return v;
168+ }
169+
170+ /** sin(theta) * dtheta(u) on the grid, f64: the undivided synthesis of the
171+ * theta shift. Smooth on the sphere, unlike dtheta(u) itself. */
172+ sinDtheta(qlm: ArrayLike<number>): Float64Array {
173+ return this.synth(this.dthetac(qlm));
174+ }
175+
176+ /** Theta-derivative, f64: synth(dthetac(u)) / sin(theta). */
177+ dtheta(qlm: ArrayLike<number>): Float64Array {
178+ const { nlat, nphi } = this.cfg;
179+ const grid = this.sinDtheta(qlm);
180+ for (let i = 0; i < nlat; i++) {
181+ const st = this.st[i];
182+ for (let j = 0; j < nphi; j++) grid[i * nphi + j] /= st;
183+ }
184+ return grid;
185+ }
186+
187+ /** Phi-derivative, f64: (dphi u)_l^m = i*m*u_l^m, then synthesize. */
188+ dphi(qlm: ArrayLike<number>): Float64Array {
189+ return this.synth(this.dphic(qlm));
169190 }
170191 }
171192
test/fluxChecks.tsadded+401−0View file
@@ -0,0 +1,401 @@
1+/**
2+ * The flux-form (six-transform) Laplace-Beltrami scheme of
3+ * docs/reduced-transforms.md, against the two things that can
4+ * silently go wrong with it:
5+ *
6+ * 1. The smoothness claim (doc Sec 2, validation Sec 7.1). The whole scheme
7+ * rests on the analysed fluxes P and Qtilde being smooth functions on the
8+ * sphere — that is a mathematical property of the p1/p2/q2 weighting, so it
9+ * is checked in f64 on the CPU, where a failure is a wrong formula and not
10+ * round-off. The fields are synthesized and re-analysed on a grid with
11+ * twice the band limit: content beyond the band is exactly the non-smooth
12+ * residue the weighting is supposed to remove.
13+ *
14+ * Two surfaces split the claim's two halves. On the *round sphere* the
15+ * correctly weighted fluxes are exactly band-limited, so their beyond-band
16+ * tail is f64 round-off, while the doc's Sec 8 counterexample
17+ * Qtilde/sin(theta) — bounded but with a phi-dependent polar limit — keeps
18+ * an algebraically decaying tail orders of magnitude above it: the
19+ * decisive smooth-vs-non-smooth discrimination, plus the closed-form check
20+ * p1 = q2 = 1, p2 = 0, r = 1/sin^2(theta). On *bumpy* (non-axisymmetric,
21+ * so the off-diagonal p2 does real work) nothing is band-limited and every
22+ * smooth field's tail is set by the weights' own spectral decay, so the
23+ * check there is the doc's relative one: P and Qtilde must sit on the same
24+ * footing as the Cartesian gradient component Algorithm 4 analyses.
25+ *
26+ * 2. The operator identity (validation Sec 7.2/7.4). The flux form and the
27+ * Cartesian-gradient form (models/schnakenberg.m vs
28+ * models/schnakenberg_alg4.m) are the same operator, so a real simulation
29+ * driven by one must track the other to fp32 accumulation — checked on a
30+ * non-axisymmetric surface, where the off-diagonal weight p2 actually does
31+ * something. The headline transform count (6 vs 12 per species per
32+ * iteration) is asserted from the compiled op sequences, not the doc.
33+ */
34+import { ShtPlan } from '../src/sht/sht.ts';
35+import { DerivPlan } from '../src/sht/deriv.ts';
36+import { ShtReference } from '../src/sht/reference.ts';
37+import { gridForLmax, lmIndex, nlmCalc, type ShtConfig } from '../src/sht/layout.ts';
38+import { ModelSession } from '../src/mgpu/session.ts';
39+import { mModelByKey, defaultParams } from '../src/mgpu/registry.ts';
40+import { Geometry } from '../src/geom/geometry.ts';
41+import { mGeometryByKey, defaultGeometryParams } from '../src/geom/registry.ts';
42+import { computeFluxMetric } from '../src/geom/metric.ts';
43+import type { Check, Log } from './analyticChecks.ts';
44+
45+/** Band limit of the test surface and field. */
46+const LMAX = 24;
47+/** Band limit of the oversampled analysis grid the tails are measured on. */
48+const LMAX_HI = 63;
49+/** Degrees at and above this count as "beyond-band tail": LMAX+1 is the last
50+ * degree with direct content, and the smooth-but-not-band-limited metric
51+ * weights spread it upward with (their own) exponentially decaying spectra,
52+ * so the window starts well above the band edge. */
53+const TAIL_START = 44;
54+
55+/** The part of a transform layout the spectral helpers need. */
56+interface Band {
57+ lmax: number;
58+ mmax: number;
59+}
60+
61+/** Per-degree spectral amplitude: E(l) = sqrt(sum_m |q_l^m|^2). */
62+function degreeEnergy(band: Band, qlm: ArrayLike<number>): Float64Array {
63+ const E = new Float64Array(band.lmax + 1);
64+ for (let m = 0; m <= band.mmax; m++) {
65+ for (let l = m; l <= band.lmax; l++) {
66+ const i = lmIndex(band.lmax, l, m);
67+ E[l] += qlm[2 * i] ** 2 + qlm[2 * i + 1] ** 2;
68+ }
69+ }
70+ for (let l = 0; l <= band.lmax; l++) E[l] = Math.sqrt(E[l]);
71+ return E;
72+}
73+
74+/** max E(l) over l >= TAIL_START, relative to max E(l) overall. */
75+function tailRel(band: Band, qlm: ArrayLike<number>): number {
76+ const E = degreeEnergy(band, qlm);
77+ let bulk = 0;
78+ let tail = 0;
79+ for (let l = 0; l <= band.lmax; l++) {
80+ if (E[l] > bulk) bulk = E[l];
81+ if (l >= TAIL_START && E[l] > tail) tail = E[l];
82+ }
83+ return tail / Math.max(bulk, 1e-300);
84+}
85+
86+/** Re-index coefficients from the lo layout into the hi layout (zero-padded). */
87+function padSpectrum(qlo: ArrayLike<number>, lo: Band, hi: Band): Float64Array {
88+ const out = new Float64Array(2 * nlmCalc(hi.lmax, hi.mmax));
89+ for (let m = 0; m <= lo.mmax; m++) {
90+ for (let l = m; l <= lo.lmax; l++) {
91+ const src = lmIndex(lo.lmax, l, m);
92+ const dst = lmIndex(hi.lmax, l, m);
93+ out[2 * dst] = qlo[2 * src];
94+ out[2 * dst + 1] = qlo[2 * src + 1];
95+ }
96+ }
97+ return out;
98+}
99+
100+/** Deterministic random band-limited spectrum with O(1) coefficients. */
101+function flatSpectrum(band: Band, seed: number): Float64Array {
102+ const nlm = nlmCalc(band.lmax, band.mmax);
103+ const q = new Float64Array(2 * nlm);
104+ let s = seed >>> 0;
105+ const rnd = () => {
106+ s ^= s << 13; s >>>= 0;
107+ s ^= s >> 17;
108+ s ^= s << 5; s >>>= 0;
109+ return (s / 4294967296) * 2 - 1;
110+ };
111+ for (let k = 0; k < 2 * nlm; k++) q[k] = rnd();
112+ for (let l = 0; l <= band.lmax; l++) q[2 * lmIndex(band.lmax, l, 0) + 1] = 0;
113+ return q;
114+}
115+
116+export interface FluxCheckOptions {
117+ /**
118+ * Run the live flux-vs-Algorithm-4 A/B (4 sessions at lmax 63). On by
119+ * default, but — like geometryChecks' sweep, and for the same reason — a
120+ * browser recompiles every session's unrolled step from scratch on software
121+ * WebGPU, so the page leaves it out unless asked (?sweep=1) to keep CI
122+ * short. The f64 smoothness checks always run; they are CPU work.
123+ */
124+ ab?: boolean;
125+}
126+
127+export async function fluxChecks(
128+ device: GPUDevice,
129+ check: Check,
130+ log: Log,
131+ opts: FluxCheckOptions = {},
132+): Promise<void> {
133+ // ---- 1a. round sphere: closed-form weights, decisive discrimination -----
134+ {
135+ const hiGrid = gridForLmax(LMAX_HI, 1);
136+ const hi = { lmax: LMAX_HI, mmax: LMAX_HI, nlat: hiGrid.nlat, nphi: hiGrid.nphi };
137+ const ref = new ShtReference(hi);
138+ const npts = hi.nlat * hi.nphi;
139+
140+ // The unit sphere needs no GPU build: analyse the closed-form embedding
141+ // on the fine grid directly, in f64.
142+ const xg = new Float64Array(npts);
143+ const yg = new Float64Array(npts);
144+ const zg = new Float64Array(npts);
145+ for (let i = 0; i < hi.nlat; i++) {
146+ const ct = ref.ct[i];
147+ const st = ref.st[i];
148+ for (let j = 0; j < hi.nphi; j++) {
149+ const phi = (2 * Math.PI * j) / hi.nphi;
150+ const k = i * hi.nphi + j;
151+ xg[k] = st * Math.cos(phi);
152+ yg[k] = st * Math.sin(phi);
153+ zg[k] = ct;
154+ }
155+ }
156+ const X = ref.analys(xg);
157+ const Y = ref.analys(yg);
158+ const Z = ref.analys(zg);
159+
160+ const sXt = [ref.sinDtheta(X), ref.sinDtheta(Y), ref.sinDtheta(Z)];
161+ const Xp = [ref.dphi(X), ref.dphi(Y), ref.dphi(Z)];
162+ const { p1, p2, q2, r } = computeFluxMetric(
163+ npts, sXt[0], sXt[1], sXt[2], Xp[0], Xp[1], Xp[2],
164+ );
165+
166+ // On the sphere the weights have a closed form: p1 = q2 = 1, p2 = 0,
167+ // r = 1/sin^2(theta) — the flux-form counterpart of geometryChecks'
168+ // closed-form V check, pinning computeFluxMetric before it is buried
169+ // under the operator. f64 throughout, so the tolerance is conditioning
170+ // at the polar rings, not fp32.
171+ let worst = 0;
172+ for (let i = 0; i < hi.nlat; i++) {
173+ const st2 = ref.st[i] * ref.st[i];
174+ for (let j = 0; j < hi.nphi; j++) {
175+ const k = i * hi.nphi + j;
176+ worst = Math.max(
177+ worst,
178+ Math.abs(p1[k] - 1),
179+ Math.abs(p2[k]),
180+ Math.abs(q2[k] - 1),
181+ Math.abs(r[k] * st2 - 1),
182+ );
183+ }
184+ }
185+ check(
186+ 'flux: sphere weights match the closed form (p1 = q2 = 1, p2 = 0, r = 1/sin^2)',
187+ worst < 1e-9,
188+ `max deviation ${worst.toExponential(2)} in f64`,
189+ );
190+
191+ // Flat random u, band-limited at LMAX. The properly weighted fluxes are
192+ // then *exactly* band-limited (P = sin(theta) dtheta u, Qtilde = dphi u),
193+ // so their beyond-band tails are pure round-off; the Sec 8 control
194+ // Qtilde/sin(theta) is not a function on the sphere and keeps a fat tail.
195+ const band = { lmax: LMAX, mmax: LMAX };
196+ const u = padSpectrum(flatSpectrum(band, 777), band, hi);
197+ const A = ref.sinDtheta(u);
198+ const B = ref.dphi(u);
199+ const P = new Float64Array(npts);
200+ const Qt = new Float64Array(npts);
201+ const control = new Float64Array(npts);
202+ for (let i = 0; i < hi.nlat; i++) {
203+ const st = ref.st[i];
204+ for (let j = 0; j < hi.nphi; j++) {
205+ const k = i * hi.nphi + j;
206+ P[k] = p1[k] * A[k] + p2[k] * B[k];
207+ Qt[k] = p2[k] * A[k] + q2[k] * B[k];
208+ control[k] = Qt[k] / st;
209+ }
210+ }
211+ const tails = {
212+ P: tailRel(hi, ref.analys(P)),
213+ Qt: tailRel(hi, ref.analys(Qt)),
214+ control: tailRel(hi, ref.analys(control)),
215+ };
216+ log(
217+ ` flux smoothness on the sphere (f64, band ${LMAX}, analysed to ${LMAX_HI}, ` +
218+ `tail l >= ${TAIL_START}): P ${tails.P.toExponential(2)}, ` +
219+ `Qt ${tails.Qt.toExponential(2)}, control ${tails.control.toExponential(2)}`,
220+ );
221+ check(
222+ 'flux: on the sphere the fluxes are band-limited and the non-smooth control is not',
223+ tails.P < 1e-10 && tails.Qt < 1e-10 &&
224+ tails.control > 1e3 * Math.max(tails.P, tails.Qt, 1e-14),
225+ `P ${tails.P.toExponential(2)}, Qt ${tails.Qt.toExponential(2)}, ` +
226+ `control ${tails.control.toExponential(2)}`,
227+ );
228+ }
229+
230+ // ---- 1b. bumpy: the fluxes sit on the Cartesian gradient's footing ------
231+ {
232+ // The surface: bumpy, the one shipped geometry that is genuinely
233+ // non-axisymmetric (g_thetaphi != 0), so the off-diagonal weight p2 is
234+ // exercised. Built by the real pipeline at LMAX, then everything below is
235+ // CPU f64 from its band-limited coefficients.
236+ const g = mGeometryByKey('bumpy')!;
237+ const { nlat, nphi } = gridForLmax(LMAX, 3);
238+ const cfg = { lmax: LMAX, mmax: LMAX, nlat, nphi };
239+ const sht = await ShtPlan.create(device, cfg);
240+ const deriv = await DerivPlan.create(device, sht);
241+ const geometry = await Geometry.create({
242+ device, sht, cfg,
243+ source: g.source,
244+ paramNames: g.params.map((p) => p.key),
245+ params: defaultGeometryParams(g),
246+ deriv,
247+ });
248+ deriv.destroy();
249+ sht.destroy();
250+
251+ const hiGrid = gridForLmax(LMAX_HI, 1);
252+ const hi = { lmax: LMAX_HI, mmax: LMAX_HI, nlat: hiGrid.nlat, nphi: hiGrid.nphi };
253+ const ref = new ShtReference(hi);
254+ const npts = hi.nlat * hi.nphi;
255+
256+ // Embedding and test field, zero-padded into the fine layout. Both are
257+ // band-limited at LMAX, so on the fine grid every derived field's content
258+ // beyond the band is genuinely the non-band-limited part of the weights —
259+ // the thing being measured — and not aliasing.
260+ const X = padSpectrum(geometry.X, cfg, hi);
261+ const Y = padSpectrum(geometry.Y, cfg, hi);
262+ const Z = padSpectrum(geometry.Z, cfg, hi);
263+ const u = padSpectrum(flatSpectrum(cfg, 777), cfg, hi);
264+
265+ // Tangents, both weightings, all f64.
266+ const sXt = [ref.sinDtheta(X), ref.sinDtheta(Y), ref.sinDtheta(Z)];
267+ const Xp = [ref.dphi(X), ref.dphi(Y), ref.dphi(Z)];
268+ const Xt = [ref.dtheta(X), ref.dtheta(Y), ref.dtheta(Z)];
269+ const { p1, p2, q2 } = computeFluxMetric(
270+ npts, sXt[0], sXt[1], sXt[2], Xp[0], Xp[1], Xp[2],
271+ );
272+
273+ const A = ref.sinDtheta(u); // sin(theta) dtheta u
274+ const B = ref.dphi(u); // dphi u
275+
276+ // The two fluxes. No non-smooth control here: on a deformed surface
277+ // every smooth field's beyond-band tail is set by the weights' own
278+ // (slowly decaying) spectra, which swamps a pole singularity at this
279+ // resolution — the sphere block above is where the discrimination has
280+ // teeth. This block asserts the doc's relative criterion instead.
281+ const P = new Float64Array(npts);
282+ const Qt = new Float64Array(npts);
283+ for (let k = 0; k < npts; k++) {
284+ P[k] = p1[k] * A[k] + p2[k] * B[k];
285+ Qt[k] = p2[k] * A[k] + q2[k] * B[k];
286+ }
287+
288+ // The known-smooth yardstick (doc Sec 2): the x component of the
289+ // Cartesian surface gradient, built the Algorithm-4 way from the inverse
290+ // metric quantities, in f64.
291+ const gradx = new Float64Array(npts);
292+ {
293+ const ut = ref.dtheta(u);
294+ const up = ref.dphi(u);
295+ for (let k = 0; k < npts; k++) {
296+ const gtt = Xt[0][k] ** 2 + Xt[1][k] ** 2 + Xt[2][k] ** 2;
297+ const gtp = Xt[0][k] * Xp[0][k] + Xt[1][k] * Xp[1][k] + Xt[2][k] * Xp[2][k];
298+ const gpp = Xp[0][k] ** 2 + Xp[1][k] ** 2 + Xp[2][k] ** 2;
299+ const det = gtt * gpp - gtp * gtp;
300+ const Vtx = (gpp * Xt[0][k] - gtp * Xp[0][k]) / det;
301+ const Vpx = (gtt * Xp[0][k] - gtp * Xt[0][k]) / det;
302+ gradx[k] = ut[k] * Vtx + up[k] * Vpx;
303+ }
304+ }
305+
306+ const tails = {
307+ P: tailRel(hi, ref.analys(P)),
308+ Qt: tailRel(hi, ref.analys(Qt)),
309+ gradx: tailRel(hi, ref.analys(gradx)),
310+ };
311+ log(
312+ ` flux smoothness on bumpy (f64, band ${LMAX}, analysed to ${LMAX_HI}, ` +
313+ `tail l >= ${TAIL_START}): P ${tails.P.toExponential(2)}, ` +
314+ `Qt ${tails.Qt.toExponential(2)}, gradx ${tails.gradx.toExponential(2)}`,
315+ );
316+ // "Matching tails" (Sec 7.1): same footing as the Cartesian component,
317+ // with an order of magnitude of headroom on top of it. A wrong weighting
318+ // (a missing sin factor, say) puts genuinely non-smooth content into P or
319+ // Qtilde and the tail lands at O(bulk), far above this.
320+ const ceiling = Math.max(30 * tails.gradx, 1e-10);
321+ check(
322+ 'flux: on bumpy, P and Qtilde tails match the Cartesian gradient component',
323+ tails.P < ceiling && tails.Qt < ceiling,
324+ `P ${tails.P.toExponential(2)}, Qt ${tails.Qt.toExponential(2)} vs ` +
325+ `ceiling ${ceiling.toExponential(2)}`,
326+ );
327+ }
328+
329+ // ---- 2. flux form vs Algorithm 4, live, on a curved surface -------------
330+ if (!(opts.ab ?? true)) {
331+ log(
332+ ' flux A/B: skipped — run `npm run test:node` (desktop Dawn) or ' +
333+ '`npm run test:gpu -- --sweep` for the flux-vs-Algorithm-4 comparison.',
334+ );
335+ } else {
336+ const geometry = mGeometryByKey('bumpy')!;
337+ const geometryParams = defaultGeometryParams(geometry);
338+ const LMAX_AB = 63;
339+ const STEPS = 20;
340+ const states: Float32Array[] = [];
341+ const xformsPerIter: number[] = [];
342+
343+ for (const key of ['schnakenberg', 'schnakenberg-alg4']) {
344+ const model = mModelByKey(key)!;
345+ const params = defaultParams(model);
346+ // Real transforms added by one solve iteration: synth/analys ops plus
347+ // dtheta/dphi (each of which contains a synthesis); the coefficient-
348+ // space dthetac/dphic shuffles are O(nlm) index gathers, not transforms.
349+ const counts: number[] = [];
350+ for (const niter of [0, 1]) {
351+ const session = await ModelSession.create({
352+ device, model, params, lmax: LMAX_AB,
353+ geometry, geometryParams, niter,
354+ });
355+ counts.push(
356+ session.describe().step.filter((l) =>
357+ l.startsWith('synth') || l.startsWith('analys') ||
358+ l.startsWith('dtheta ') || l.startsWith('dphi '),
359+ ).length,
360+ );
361+ if (niter === 1) {
362+ session.seed(1);
363+ session.step(STEPS);
364+ states.push(await session.read('U'));
365+ }
366+ session.destroy();
367+ }
368+ xformsPerIter.push(counts[1] - counts[0]);
369+ }
370+
371+ // The headline number, from the compiled op sequences: 6 transforms per
372+ // species per iteration against Algorithm 4's 12 (2 species here).
373+ check(
374+ 'flux: 6 transforms per species per iteration, versus 12',
375+ xformsPerIter[0] === 12 && xformsPerIter[1] === 24,
376+ `flux form adds ${xformsPerIter[0]} transforms/iteration, ` +
377+ `Algorithm 4 adds ${xformsPerIter[1]}`,
378+ );
379+
380+ // Same operator, same discretization, different arithmetic path: after
381+ // STEPS steps the two states may differ only by fp32 accumulation. A
382+ // formulation error (wrong weight, wrong shift, missing sin) would show
383+ // up at O(1), not O(1e-3). Identical states would mean the A/B compared
384+ // one path to itself.
385+ let worst = 0;
386+ let identical = true;
387+ let finite = true;
388+ for (let i = 0; i < states[0].length; i++) {
389+ const d = Math.abs(states[0][i] - states[1][i]);
390+ if (d > worst) worst = d;
391+ if (states[0][i] !== states[1][i]) identical = false;
392+ if (!Number.isFinite(states[0][i]) || !Number.isFinite(states[1][i])) finite = false;
393+ }
394+ check(
395+ 'flux: tracks the Algorithm-4 reference through a real simulation',
396+ finite && !identical && worst < 5e-3,
397+ `max |U_flux - U_alg4| = ${worst.toExponential(2)} after ${STEPS} steps ` +
398+ `on bumpy at lmax ${LMAX_AB}`,
399+ );
400+ }
401+}
test/geometryChecks.tsmodified+7−7View file
@@ -310,14 +310,14 @@ export async function geometryChecks(
310310 `${ops.join(' < ')} ops for ${counts.join(', ')} iterations`,
311311 );
312312 // Unrolling has to be exactly linear in the trip count: the body planned
313- // once per iteration, no more and no less. Per species per iteration: 8
314- // dtheta/dphi + 4 analys transforms (Algorithm 3's cost, applied to the
315- // field and to each of its three Cartesian gradient components) plus 15
316- // generated kernels -- see test/modelChecks.ts's KERNELS_PER_ITERATION,
317- // which counts the kernels alone; this counts every op, transforms
318- // included.
313+ // once per iteration, no more and no less. Per species per iteration: 3
314+ // synths + 3 analyses (the flux-form matvec's six real transforms,
315+ // docs/reduced-transforms.md Sec 4) + 4 coefficient-space
316+ // dthetac/dphic shuffles plus 9 generated kernels -- see
317+ // test/modelChecks.ts's KERNELS_PER_ITERATION, which counts the kernels
318+ // alone; this counts every op, transforms and shuffles included.
319319 const perIteration = ops[1] - ops[0];
320- const want = 54;
320+ const want = 38;
321321 check(
322322 'loop: unrolling is exactly linear in the trip count',
323323 perIteration === want && ops[2] - ops[0] === 4 * perIteration,
test/modelChecks.tsmodified+14−9View file
@@ -28,6 +28,7 @@ const EXPECTED_KERNELS: Record<string, number> = {
2828 schnakenberg: 7,
2929 brusselator: 7,
3030 allencahn: 3,
31+ 'schnakenberg-alg4': 7,
3132 };
3233
3334 /**
@@ -36,16 +37,20 @@ const EXPECTED_KERNELS: Record<string, number> = {
3637 * count is a byproduct of exactly how its expression tree happens to fuse,
3738 * not a clean per-species multiple, so this is measured per model rather
3839 * than derived from `model.species.length`). Each species' correction is
39- * Algorithm 3 of evolving_surface/notes/algos.tex: a surface gradient
40- * (dtheta/dphi contracted through the metric), reanalysed per Cartesian
41- * component and differentiated again, recombined into the divergence, plus
42- * the round-sphere eigenvalue added back — see models/schnakenberg.m and
43- * docs/richardson-iteration.md.
40+ * the flux-form Laplace-Beltrami matvec of
41+ * docs/reduced-transforms.md Sec 4: the two sin-weighted
42+ * derivative synths, the pointwise flux combination through p1/p2/q2, the
43+ * two flux analyses, the re-shifted divergence and its r-scaled synthesis,
44+ * plus the round-sphere eigenvalue added back — see models/schnakenberg.m
45+ * and docs/richardson-iteration.md. `schnakenberg-alg4` keeps the original
46+ * Cartesian-gradient form (Algorithm 3/4 of evolving_surface/notes/algos.tex)
47+ * as a live reference, with its original counts.
4448 */
4549 const KERNELS_PER_ITERATION: Record<string, number> = {
46- schnakenberg: 30,
47- brusselator: 30,
48- allencahn: 14,
50+ schnakenberg: 18,
51+ brusselator: 18,
52+ allencahn: 9,
53+ 'schnakenberg-alg4': 30,
4954 };
5055
5156 const LMAX = 31;
@@ -57,7 +62,7 @@ export async function modelChecks(
5762 check: Check,
5863 log: Log,
5964 ): Promise<void> {
60- check('models: registry populated', mModels.length === 3, `${mModels.length} models`);
65+ check('models: registry populated', mModels.length === 4, `${mModels.length} models`);
6166
6267 // The app formats the run it is showing into a `npm run bench` command and
6368 // the benchmark parses it back. That is only worth anything if the round
test/test-page.tsmodified+2−0View file
@@ -27,6 +27,7 @@ import { transformChecks } from './transformChecks.ts';
2727 import { analyticChecks } from './analyticChecks.ts';
2828 import { modelChecks } from './modelChecks.ts';
2929 import { geometryChecks } from './geometryChecks.ts';
30+import { fluxChecks } from './fluxChecks.ts';
3031
3132 declare global {
3233 interface Window {
@@ -218,6 +219,7 @@ async function main(): Promise<void> {
218219 // The sweep is opt-in here (?sweep=1): it is a few seconds on desktop Dawn
219220 // but minutes in a browser, where each session recompiles its unrolled step.
220221 await geometryChecks(device, check, log, { sweep: q.has('sweep') });
222+ await fluxChecks(device, check, log, { ab: q.has('sweep') });
221223
222224 window.__RESULTS__ = { ok: failures === 0, lines };
223225 log(failures === 0 ? 'ALL PASS' : `${failures} FAILURE(S)`);
test/transformChecks.tsmodified+13−0View file
@@ -110,6 +110,19 @@ export async function transformChecks(
110110 errDtheta < 1e-4 && errDphi < 1e-4,
111111 `dtheta ${errDtheta.toExponential(2)}, dphi ${errDphi.toExponential(2)}`,
112112 );
113+
114+ // The undivided theta derivative sin(theta)*dtheta(u) — the flux-form
115+ // Laplace-Beltrami scheme's step 1 and the flux-metric precompute's
116+ // input — is the same shuffle+synthesis with the divide skipped, so it
117+ // gets the same oracle.
118+ const sinDthetaGpu = await deriv.sinDtheta(new Float32Array(q64));
119+ const sinDthetaCpu = ref.sinDtheta(q64);
120+ const errSinDtheta = relL2(sinDthetaGpu, sinDthetaCpu);
121+ check(
122+ 'deriv: WGSL fp32 sinDtheta (undivided) vs f64 CPU reference',
123+ errSinDtheta < 1e-4,
124+ `sinDtheta ${errSinDtheta.toExponential(2)}`,
125+ );
113126 deriv.destroy();
114127 }
115128