Split the flux-form divergence against the round sphere - concept-collection/turing-surface

concept-collection / turing-surface

Split the flux-form divergence against the round sphere

The six-transform Laplace-Beltrami correction pays both powers of 1/sin(theta) at once, in r = 1/(J sin^2 theta), multiplying a bracket that must cancel to O(sin^2 theta) at the poles. Absolute fp32 round-off in the bracket therefore comes back amplified by 1/sin^2(theta) ~ L^2 -- the epsilon*L^2 that docs/reduced-transforms.md Sec 5 predicted, against Algorithm 4's epsilon*L. The doc concluded that this would not matter, because the final analysis suppresses the polar spike back down to ~epsilon*sqrt(L) in the coefficients the solve actually sees. The suppression is real; the conclusion is not. The residue is a *static* polar forcing, re-injected every step, and a Turing pattern is seeded by whatever is largest in its unstable band -- so a forcing four orders below the field selects the nucleation site. Every run grew a spot at the pole, on every seed, on the sphere and on curved surfaces alike. Measured on the default ellipsoid at lmax 63: from the exact uniform steady state with the Turing band switched off, three Richardson iterations leave a standing perturbation 50x the no-correction floor (1.4e-5 vs 2.8e-7), where Algorithm 4 leaves ~1x. Keep r off the round sphere. With p1 = 1 + dp1 and q2 = 1 + dq2 (p2 is already zero on the sphere), the sphere's share of the divergence is the cancelling part and is known in closed form -- lap_s is diagonal -- so lap_g u = -jinv * lap_s(u) + r * (sin(theta) dtheta(P') + dphi(Q')) with P' = dp1*A + p2*B, Q' = p2*A + dq2*B and jinv = 1/J = r sin^2(theta) bounded. Only the geometry deviation meets the concentrated division. dp1 and dq2 are formed in f64 at precompute time: on a near-sphere they are the small quantity, and subtracting 1 in f32 on device would lose them. Cost is one extra synthesis per species per iteration for lap_s -- 7 transforms, not 6 -- and it rides in the gradient's existing grouped call, so it measures at ~10% of a step against 3x for reverting to Algorithm 4. Results on the ellipsoid, against the 12-transform reference: flux split alg4 noise floor, niter=3 1.4e-5 2.4e-7 4.0e-7 nucleation, lmax 127 t=14 9.4 t=23 0.4 t=23 0.4 rel L2 vs alg4, L=63 3.5e-4 8.9e-6 -- rel L2 vs alg4, L=127 7.2e-3 6.0e-6 -- ms/step, niter=3 0.68 0.74 2.26 The split does not merely stop the pole nucleating: it lands on Algorithm 4's answer to two decimals everywhere tested (ellipsoid, blob, sphere; lmax 63 and 127), and is 40-1200x closer to the reference than the unsplit form. The polar term was dominating accuracy, not only the seeding. test/fluxChecks.ts gains the check that catches this: from the exact uniform state with D1 = D2, the niter-3 round-off must stay under 5x the niter-0 floor. It fails at 50.4x on the previous models and passes at 0.9x now -- the operator A/B tests cannot see it, since they compare operators rather than watch what a run nucleates from. Two adjacent fixes to counts the change moves: the op-count assertions across the test suite, and two fixtures in modelChecks.ts whose string surgery had silently stopped matching allencahn.m and were passing vacuously -- they now throw if the fixture drifts again. README's op counts were already stale before this change (51 ops / 19 per species where the compiler emitted 50 / 16), so they are set to measured values rather than adjusted.

Dan Fortunato <dan.fortunato@gmail.com> committed 13 days ago commit 3d078cf1e121 parent beac00a Browse files

10 changed files+261−74

README.mdmodified+18−12View file

@@ -273,15 +273,21 @@ Two formulations ship:

273273

274274 1. **The flux form** (above, all three models): `lap_g u` as the weighted

275275 divergence of two weighted fluxes of the sin-scaled derivatives. The

276- weights `p1, p2, q2, r` are grid arrays precomputed once per surface from

277- the embedding's θ/φ tangents

276+ weights `p2, r, dp1, dq2, jinv` are grid arrays precomputed once per

277+ surface from the embedding's θ/φ tangents

278278 ([`src/geom/metric.ts`](src/geom/metric.ts)), chosen so that **every field

279279 that gets analysed is a smooth function on the sphere** — the property

280280 that makes spherical-harmonic analysis meaningful, and the entire

281- difficulty near the poles. Cost: **5 Legendre transforms** per species

282- per iteration (3 syntheses + 2 analyses; the phi flux never needs the

283- Legendre basis — `dphig` differentiates it on the grid with two FFT

284- stages, masking m past the top-degree filter — and `dthetac`/`dphic`

281+ difficulty near the poles. The divergence is split against the round

282+ sphere: the sphere's share of it is `-jinv .* lap_s(u)`, exact in

283+ spectral space, so `r ~ 1/sin²θ` multiplies only the geometry deviation.

284+ Without that split `r` amplifies the polar round-off of the whole flux

285+ into a static forcing that nucleates a spot at the pole on every seed.

286+ Cost: **6 Legendre transforms** per species per iteration (4 syntheses —

287+ two gradient, one divergence, one for the sphere's `-lam .* u`, which

288+ rides in the gradient's batch — plus 2 analyses; the phi flux never needs

289+ the Legendre basis, `dphig` differentiates it on the grid with two FFT

290+ stages, masking m past the top-degree filter, and `dthetac`/`dphic`

285291 are O(nlm) coefficient shuffles). The derivation, the smoothness

286292 argument and the fp32 error analysis are in

287293 [docs/reduced-transforms.md](docs/reduced-transforms.md).

@@ -328,8 +334,8 @@ Two consequences worth stating:

328334 each loop body assigns before the pass and refuses the ones that escape, so

329335 that case is a compile error rather than a stale read.

330336

331-Unrolling is exactly linear in the trip count: 19 GPU ops per species per

332-iteration (6 transforms, 4 coefficient shuffles, 9 kernels), asserted in the

337+Unrolling is exactly linear in the trip count: 18 GPU ops per species per

338+iteration (7 transforms, 3 coefficient shuffles, 8 kernels), asserted in the

333339 tests.

334340

335341 ## MATLAB, compiled to WebGPU

@@ -345,8 +351,8 @@ operations whose type rules numbl learns from a `.mtoc2.js` workspace file, and

345351 which the backend maps onto the spherical-harmonic pipelines. Anything it cannot

346352 express is refused at compile time with a source position.

347353

348-The Schnakenberg step compiles to 51 GPU operations at one solve iteration:

349-16 transforms, 8 coefficient-space shuffles, 25 generated kernels, and 2

354+The Schnakenberg step compiles to 50 GPU operations at one solve iteration:

355+18 transforms, 6 coefficient-space shuffles, 24 generated kernels, and 2

350356 buffer copies feeding the new state back.

351357

352358 **Transforms batch.** The expensive part of every Legendre stage is

@@ -356,7 +362,7 @@ take multiple fields, and a grouped call runs as one batched dispatch: one

356362 walk of the recurrence, one accumulator lane per field —

357363

358364 ```matlab

359-[Ftu, Fpu, Ftv, Fpv] = synth(vtu, vpu, vtv, vpv); % one Legendre dispatch

365+[Ftu, Fpu, Ftv, Fpv, Su, Sv] = synth(vtu, vpu, vtv, vpv, lam .* Fu, lam .* Fv);

360366 ```

361367

362368 The grouping is a promise of independence, never of a lane width: the

@@ -367,7 +373,7 @@ the same source runs anywhere. Ungrouped transforms that happen to sit on

367373 consecutive independent lines are batched the same way. Per-lane arithmetic

368374 is identical to the scalar kernels', so batched and scalar plans produce

369375 bit-identical states, asserted in the tests along with compile-time refusal

370-of a group that drops one of its outputs. All 16 transforms of the step

376+of a group that drops one of its outputs. All 16 Legendre transforms of the step

371377 above land in batches, worth ~25% of the whole step (0.88 vs 1.14 ms/step at

372378 lmax 127, 2 iterations, on bumpy).

373379

docs/reduced-transforms.mdmodified+42−8View file

@@ -2,7 +2,8 @@

22

33 **Summary.** Algorithm 4 costs 12 transforms per matvec (8 syntheses, 4 analyses). A flux-form

44 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`.

5+same operator in **6 transforms** (4 syntheses, 2 analyses), or **7** with the divergence split

6+against the round sphere that §5 turned out to require. Notation follows `algos.pdf`.

67

78 ---

89

@@ -167,14 +168,47 @@ every node alike, multiplying by $r \sim L^2$ recovers the signal and inflates t

167168 | Algorithm 4 | $\sin\theta$, $\sin\theta$ | separated by $\mathcal{A}$ | $\varepsilon L$ |

168169 | Six-transform | $\sin^2\theta$ | all at the end | $\varepsilon L^2$ |

169170

170-**This likely does not reach the returned coefficients.** Step 7's analysis suppresses the spike

171-exactly as line 5 does today: $L^{-2}\cdot L^{1/2}\cdot\varepsilon L^2 = \varepsilon L^{1/2}$,

172-comparable to the ordinary $\varepsilon\sqrt{L}$ accumulation of a transform pair — and the new

173-scheme runs half as many transforms, lowering that baseline. Inside the implicit solve, GMRES sees

174-only coefficients, so the extra power should be invisible.

171+**It does reach the returned coefficients, and it matters.** The suppression argument above is

172+right as far as it goes — step 7's analysis knocks the spike down to $\varepsilon L^{1/2}$, and this

173+document originally concluded from that the extra power would be invisible inside the solve. It is

174+not, and the reason is not about accuracy. Measured on the default ellipsoid at $L=63$: starting

175+from the exact uniform steady state, three Richardson iterations per step leave a standing

176+coefficient-space perturbation $50\times$ the no-correction floor ($1.4\times10^{-5}$ vs

177+$2.8\times10^{-7}$), against $\sim\!1\times$ for Algorithm 4. That perturbation is static, polar,

178+and re-injected every step. In a Turing problem the pattern is seeded by whatever is largest in the

179+unstable band, so a forcing four orders below the field selects the nucleation site: the run grows a

180+spot at the pole, on every seed, regardless of the initial condition.

181+

182+**The fix is to keep $r$ off the round sphere.** Write $p_1 = 1 + \delta p_1$, $q_2 = 1 + \delta q_2$

183+($p_2$ is already zero on the sphere). The sphere's share of the divergence is the cancelling part,

184+and it is known in closed form: $\sin\theta\,\partial_\theta A + \partial_\varphi B =

185+-\sin^2\theta\,\Delta_{S^2}u$, and $\Delta_{S^2}$ is diagonal. So

186+

187+$$\Delta_\Gamma u = -\frac{1}{J}\,\Delta_{S^2}u \;+\; r\,(\sin\theta\,\partial_\theta P'

188+ + \partial_\varphi \tilde{Q}'), \qquad P' = \delta p_1 A + p_2 B, \quad

189+ \tilde{Q}' = p_2 A + \delta q_2 B$$

190+

191+with $1/J = r\sin^2\theta$ bounded. Only the geometry *deviation* now meets the concentrated

192+division. Cost: one extra synthesis per species per iteration for $-\lambda u$ — 7 transforms, not

193+6 — which batches into the gradient's existing grouped call and measures at ~10% of a step, against

194+$3\times$ for reverting to Algorithm 4. $\delta p_1$ and $\delta q_2$ must be formed in float64 at

195+precompute time (`src/geom/geometry.ts`): on a near-sphere they *are* the small quantity, and

196+subtracting 1 in float32 on device would lose them.

197+

198+Measured against Algorithm 4 through a real run (relative $L^2$ of $u$ at $t=8$, $\texttt{niter}=6$):

199+

200+| | plain flux | sphere-split |

201+|---|---|---|

202+| ellipsoid, $L=63$ | $3.5\times10^{-4}$ | $8.9\times10^{-6}$ |

203+| blob, $L=63$ | $4.4\times10^{-4}$ | $8.5\times10^{-6}$ |

204+| ellipsoid, $L=127$ | $7.2\times10^{-3}$ | $6.0\times10^{-6}$ |

205+

206+and the polar noise gain is asserted in `test/fluxChecks.ts`, which fails at $50\times$ on the

207+unsplit form.

175208

176-It matters only if grid values of $\Delta_\Gamma u$ are consumed directly: a nonlinear reaction

177-term, max-norm diagnostics, or an adaptive error estimator.

209+The residual $\varepsilon L^2$ still applies to grid values of $\Delta_\Gamma u$ consumed directly

210+— a nonlinear reaction term, max-norm diagnostics, an adaptive error estimator — for the deviation

211+part alone.

178212

179213 Algorithm 1 line 7 already divides by $\sin^2\theta$, so the code is exposed to $\varepsilon L^2$

180214 today — just on the second-derivative path, which the Laplacian never touches.

models/allencahn.mmodified+6−6View file

@@ -2,7 +2,7 @@

22 %

33 % du/dt = eps2*lap_g(u) + u - u^3

44 %

5-% Same scheme as models/schnakenberg.m.

5+% Same scheme as models/schnakenberg.m, sphere-split flux divergence included.

66

77 % Seeded from a smooth random field -- see models/schnakenberg.m.

88 function [U, u] = init(lam3, gx, gy, gz)

@@ -10,7 +10,7 @@ function [U, u] = init(lam3, gx, gy, gz)

1010 u = synth(U);

1111 end

1212

13-function [Un, u] = step(U, lam, filt, gx, gy, gz, p1, p2, q2, r, jhat, eps2, dt, niter)

13+function [Un, u] = step(U, lam, filt, gx, gy, gz, p2, r, dp1, dq2, jinv, jhat, eps2, dt, niter)

1414 u = synth(U);

1515

1616 Bu = U + dt * analys(u - u.^3);

@@ -27,15 +27,15 @@ function [Un, u] = step(U, lam, filt, gx, gy, gz, p1, p2, q2, r, jhat, eps2, dt,

2727 Fu = Un .* filt;

2828 vtu = dthetac(Fu);

2929 vpu = dphic(Fu);

30- [Ftu, Fpu] = synth(vtu, vpu);

31- Pu = p1 .* Ftu + p2 .* Fpu;

32- Qu = p2 .* Ftu + q2 .* Fpu;

30+ [Ftu, Fpu, Su] = synth(vtu, vpu, lam .* Fu);

31+ Pu = dp1 .* Ftu + p2 .* Fpu;

32+ Qu = p2 .* Ftu + dq2 .* Fpu;

3333 PAu = analys(Pu);

3434 Pcu = PAu .* filt;

3535 scu = dthetac(Pcu);

3636 Lu = synth(scu);

3737 dQu = dphig(Qu);

38- lapu = r .* (Lu + dQu);

38+ lapu = r .* (Lu + dQu) - jinv .* Su;

3939 dLu = (analys(lapu) + lamJ .* Un) .* filt;

4040

4141 Un = (Bu + (dt * eps2) * dLu) ./ (1 + (dt * eps2) * lamJ);

models/brusselator.mmodified+11−9View file

@@ -4,7 +4,9 @@

44 % dv/dt = D2*lap_g(v) + B*u - u^2*v

55 %

66 % Same scheme as models/schnakenberg.m, including the grouped transforms:

7-% [a, b] = synth(x, y) runs the group as batched Legendre dispatches.

7+% [a, b] = synth(x, y) runs the group as batched Legendre dispatches, and the

8+% sphere-split flux divergence that keeps r ~ 1/sin^2(theta) off the round

9+% sphere's share of the operator.

810

911 % Seeded from a smooth random field -- see models/schnakenberg.m.

1012 function [U, V, u, v] = init(lam3, gx, gy, gz, A, B)

@@ -13,7 +15,7 @@ function [U, V, u, v] = init(lam3, gx, gy, gz, A, B)

1315 [u, v] = synth(U, V);

1416 end

1517

16-function [Un, Vn, u, v] = step(U, V, lam, filt, gx, gy, gz, p1, p2, q2, r, jhat, A, B, D1, D2, dt, niter)

18+function [Un, Vn, u, v] = step(U, V, lam, filt, gx, gy, gz, p2, r, dp1, dq2, jinv, jhat, A, B, D1, D2, dt, niter)

1719 [u, v] = synth(U, V);

1820 uuv = u .* u .* v;

1921

@@ -38,11 +40,11 @@ function [Un, Vn, u, v] = step(U, V, lam, filt, gx, gy, gz, p1, p2, q2, r, jhat,

3840 vpu = dphic(Fu);

3941 vtv = dthetac(Fv);

4042 vpv = dphic(Fv);

41- [Ftu, Fpu, Ftv, Fpv] = synth(vtu, vpu, vtv, vpv);

42- Pu = p1 .* Ftu + p2 .* Fpu;

43- Qu = p2 .* Ftu + q2 .* Fpu;

44- Pv = p1 .* Ftv + p2 .* Fpv;

45- Qv = p2 .* Ftv + q2 .* Fpv;

43+ [Ftu, Fpu, Ftv, Fpv, Su, Sv] = synth(vtu, vpu, vtv, vpv, lam .* Fu, lam .* Fv);

44+ Pu = dp1 .* Ftu + p2 .* Fpu;

45+ Qu = p2 .* Ftu + dq2 .* Fpu;

46+ Pv = dp1 .* Ftv + p2 .* Fpv;

47+ Qv = p2 .* Ftv + dq2 .* Fpv;

4648 [PAu, PAv] = analys(Pu, Pv);

4749 Pcu = PAu .* filt;

4850 Pcv = PAv .* filt;

@@ -51,8 +53,8 @@ function [Un, Vn, u, v] = step(U, V, lam, filt, gx, gy, gz, p1, p2, q2, r, jhat,

5153 [Lu, Lv] = synth(scu, scv);

5254 dQu = dphig(Qu);

5355 dQv = dphig(Qv);

54- lapu = r .* (Lu + dQu);

55- lapv = r .* (Lv + dQv);

56+ lapu = r .* (Lu + dQu) - jinv .* Su;

57+ lapv = r .* (Lv + dQv) - jinv .* Sv;

5658 [LAu, LAv] = analys(lapu, lapv);

5759 dLu = (LAu + lamJ .* Un) .* filt;

5860 dLv = (LAv + lamJ .* Vn) .* filt;

models/schnakenberg.mmodified+35−14View file

@@ -9,11 +9,20 @@

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.

1111 %

12-% The correction evaluates lap_g in flux form -- 6 transforms per species

12+% The correction evaluates lap_g in flux form -- 7 transforms per species

1313 % per iteration where the Cartesian-gradient form (Algorithm 4 of

1414 % evolving_surface/notes/algos.tex) needs 12. See

1515 % docs/reduced-transforms.md, and models/schnakenberg_alg4.m

1616 % for the original form kept as a live reference.

17+%

18+% The flux divergence is split against the round sphere: the sphere's share

19+% of it is -jinv*lap_s(u), exact in spectral space, and only the geometry

20+% *deviation* meets r ~ 1/sin^2(theta). Without that split the concentrated

21+% division amplifies the polar roundoff of the whole flux, and since a

22+% Turing pattern is seeded by whatever is largest in its unstable band, the

23+% amplified polar noise -- static, and re-injected every step -- picks the

24+% nucleation site and grows a spot at the pole. See docs/reduced-transforms.md

25+% Sec 5.

1726

1827 % The uniform steady state, perturbed by a smooth random field: chebfun's

1928 % randnfun3 on the surface's bounding box, restricted to the surface by

@@ -28,7 +37,7 @@ function [U, V, u, v] = init(lam3, gx, gy, gz, a, b)

2837 [u, v] = synth(U, V);

2938 end

3039

31-function [Un, Vn, u, v] = step(U, V, lam, filt, gx, gy, gz, p1, p2, q2, r, jhat, a, b, D1, D2, dt, niter)

40+function [Un, Vn, u, v] = step(U, V, lam, filt, gx, gy, gz, p2, r, dp1, dq2, jinv, jhat, a, b, D1, D2, dt, niter)

3241 % Grouped transforms -- [a, b] = synth(x, y) -- are explicit batching:

3342 % output k is the transform of input k, and the whole group runs as one

3443 % batched Legendre dispatch, or as many as the device's lane width allows

@@ -61,16 +70,28 @@ function [Un, Vn, u, v] = step(U, V, lam, filt, gx, gy, gz, p1, p2, q2, r, jhat,

6170 for k = 1:niter

6271 % dlap = lap_g - lap_s at the current iterate, in flux form

6372 % (docs/reduced-transforms.md Sec 4). The sin-weighted derivatives

64- % sin(theta)*dtheta(u) and dphi(u) -- both smooth on the sphere,

73+ % A = sin(theta)*dtheta(u) and B = dphi(u) -- both smooth on the sphere,

6574 % synthesized straight from the dthetac/dphic coefficient shuffles --

66- % are combined pointwise through the precomputed weights p1,p2,q2 into

67- % two fluxes P,Q, also smooth. The theta flux P goes back to

75+ % are combined pointwise through the precomputed weights into two

76+ % fluxes P,Q, also smooth. The theta flux P goes back to

6877 % coefficients, through the same shuffle again, and is synthesized as

6978 % sin(theta)*dtheta(P); the phi flux Q never leaves the grid -- d/dphi

7079 % is diagonal in the Fourier index, so dphig differentiates it with two

7180 % FFT stages and no Legendre work (masking m past filt's reach). Their

7281 % sum, scaled by r, is lap_g(u). The only division by sin(theta)

73- % anywhere is folded into p1,p2,q2,r at precompute time.

82+ % anywhere is folded into the weights at precompute time.

83+ %

84+ % The weights here are the *sphere-subtracted* ones: p1 = 1 + dp1 and

85+ % q2 = 1 + dq2 (p2 is zero on the sphere already), so P,Q below are the

86+ % deviation fluxes P' = P - A, Q' = Q - B. What that leaves out is the

87+ % round sphere's own divergence, sin(theta)*dtheta(A) + dphi(B) =

88+ % -sin^2(theta)*lap_s(u), which needs no flux machinery at all: lap_s is

89+ % diagonal, so it is -lam.*Fu synthesized once (S below, riding along in

90+ % the gradient's batched synthesis) and scaled by the bounded

91+ % jinv = 1/J = r*sin^2(theta). r therefore multiplies only the deviation

92+ % -- the difference between this and multiplying the whole flux is two

93+ % orders of magnitude of polar roundoff, and it is what keeps a pattern

94+ % from nucleating at the pole (src/geom/geometry.ts, dp1/dq2/jinv).

7495 % lamJ.*Un adds back the preconditioner's -lap_s(Un)/jhat, since lam

7596 % holds +l(l+1). filt zeroes the top two degrees, where the derivative

7697 % recurrences cannot exactly represent a derivative -- and the correction

@@ -80,7 +101,7 @@ function [Un, Vn, u, v] = step(U, V, lam, filt, gx, gy, gz, p1, p2, q2, r, jhat,

80101 % point is the undiffused Bu), and the two species un-diffuse at

81102 % different rates -- a spurious Turing band at the band edge.

82103 %

83- % The two species share each grouped call: the four gradient

104+ % The two species share each grouped call: the six gradient-and-sphere

84105 % syntheses, the two theta-flux analyses, the two divergence syntheses

85106 % and the two final analyses each run as one batched dispatch.

86107 Fu = Un .* filt;

@@ -89,11 +110,11 @@ function [Un, Vn, u, v] = step(U, V, lam, filt, gx, gy, gz, p1, p2, q2, r, jhat,

89110 vpu = dphic(Fu);

90111 vtv = dthetac(Fv);

91112 vpv = dphic(Fv);

92- [Ftu, Fpu, Ftv, Fpv] = synth(vtu, vpu, vtv, vpv);

93- Pu = p1 .* Ftu + p2 .* Fpu;

94- Qu = p2 .* Ftu + q2 .* Fpu;

95- Pv = p1 .* Ftv + p2 .* Fpv;

96- Qv = p2 .* Ftv + q2 .* Fpv;

113+ [Ftu, Fpu, Ftv, Fpv, Su, Sv] = synth(vtu, vpu, vtv, vpv, lam .* Fu, lam .* Fv);

114+ Pu = dp1 .* Ftu + p2 .* Fpu;

115+ Qu = p2 .* Ftu + dq2 .* Fpu;

116+ Pv = dp1 .* Ftv + p2 .* Fpv;

117+ Qv = p2 .* Ftv + dq2 .* Fpv;

97118 [PAu, PAv] = analys(Pu, Pv);

98119 Pcu = PAu .* filt;

99120 Pcv = PAv .* filt;

@@ -102,8 +123,8 @@ function [Un, Vn, u, v] = step(U, V, lam, filt, gx, gy, gz, p1, p2, q2, r, jhat,

102123 [Lu, Lv] = synth(scu, scv);

103124 dQu = dphig(Qu);

104125 dQv = dphig(Qv);

105- lapu = r .* (Lu + dQu);

106- lapv = r .* (Lv + dQv);

126+ lapu = r .* (Lu + dQu) - jinv .* Su;

127+ lapv = r .* (Lv + dQv) - jinv .* Sv;

107128 [LAu, LAv] = analys(lapu, lapv);

108129 dLu = (LAu + lamJ .* Un) .* filt;

109130 dLv = (LAv + lamJ .* Vn) .* filt;

src/geom/geometry.tsmodified+34−0View file

@@ -94,6 +94,27 @@ export class Geometry {

9494 readonly p2: Float32Array;

9595 readonly q2: Float32Array;

9696 readonly r: Float32Array;

97+ /**

98+ * The same flux weights with the round sphere subtracted off, plus the

99+ * bounded 1/J — what lets a model evaluate lap_g without ever multiplying

100+ * the *whole* flux divergence by r ~ 1/sin^2(theta). Writing p1 = 1 + dp1,

101+ * q2 = 1 + dq2 (p2 is already a pure deviation, zero on the sphere) splits

102+ * the divergence into a round-sphere part, whose cancelling bracket

103+ * sin(theta) dtheta(A) + dphi(B) = -sin^2(theta) lap_s u is known exactly in

104+ * spectral space, and a remainder:

105+ *

106+ * lap_g u = -jinv * lap_s u + r * (sin(theta) dtheta(P') + dphi(Q'))

107+ *

108+ * with P' = dp1*A + p2*B, Q' = p2*A + dq2*B. Only the remainder meets the

109+ * concentrated division, so the polar roundoff gain drops by |P'|/|P|

110+ * instead of applying to the full flux. Subtracting 1 in f64 here is the

111+ * point: on a near-sphere dp1 is the small quantity, and forming it as an

112+ * f32 difference in the .m would lose it. See docs/reduced-transforms.md

113+ * Sec 5 and models/schnakenberg.m.

114+ */

115+ readonly dp1: Float32Array;

116+ readonly dq2: Float32Array;

117+ readonly jinv: Float32Array;

97118 /**

98119 * Preconditioner scale for the implicit solve (docs/reduced-transforms.md

99120 * Sec 10). At high degree the Richardson iteration's per-mode factor is

@@ -129,6 +150,7 @@ export class Geometry {

129150 Vtx: Float32Array; Vty: Float32Array; Vtz: Float32Array;

130151 Vpx: Float32Array; Vpy: Float32Array; Vpz: Float32Array;

131152 p1: Float32Array; p2: Float32Array; q2: Float32Array; r: Float32Array;

153+ dp1: Float32Array; dq2: Float32Array; jinv: Float32Array;

132154 Jhat: number; muMin: number; muMax: number; Jmin: number; Jmax: number;

133155 }) {

134156 this.x = init.x;

@@ -147,6 +169,9 @@ export class Geometry {

147169 this.p2 = init.p2;

148170 this.q2 = init.q2;

149171 this.r = init.r;

172+ this.dp1 = init.dp1;

173+ this.dq2 = init.dq2;

174+ this.jinv = init.jinv;

150175 this.Jhat = init.Jhat;

151176 this.muMin = init.muMin;

152177 this.muMax = init.muMax;

@@ -211,12 +236,20 @@ export class Geometry {

211236 let muMax = 0;

212237 let Jmin = Infinity;

213238 let Jmax = 0;

239+ // The sphere-subtracted weights ride along on this loop: 1/J is already

240+ // being formed here, and dp1/dq2 want the same f64 arithmetic.

241+ const dp1 = new Float32Array(npts);

242+ const dq2 = new Float32Array(npts);

243+ const jinv = new Float32Array(npts);

214244 for (let i = 0; i < cfg.nlat; i++) {

215245 const ct = sht.cosTheta[i];

216246 const st2 = Math.max(0, 1 - ct * ct);

217247 for (let j = 0; j < cfg.nphi; j++) {

218248 const k = i * cfg.nphi + j;

219249 const invJ = flux.r[k] * st2;

250+ dp1[k] = flux.p1[k] - 1;

251+ dq2[k] = flux.q2[k] - 1;

252+ jinv[k] = invJ;

220253 const s11 = flux.p1[k] * invJ;

221254 const s12 = flux.p2[k] * invJ;

222255 const s22 = flux.q2[k] * invJ;

@@ -237,6 +270,7 @@ export class Geometry {

237270 p2: new Float32Array(flux.p2),

238271 q2: new Float32Array(flux.q2),

239272 r: new Float32Array(flux.r),

273+ dp1, dq2, jinv,

240274 Jhat, muMin, muMax, Jmin, Jmax,

241275 });

242276 }

src/mgpu/model.tsmodified+11−1View file

@@ -84,6 +84,12 @@ export interface GeometryBuffers {

8484 p2: Float32Array;

8585 q2: Float32Array;

8686 r: Float32Array;

87+ /** The same weights with the round sphere subtracted (Geometry.dp1/dq2/jinv)

88+ * — the sphere-split form of the flux divergence, which keeps r off the

89+ * round-sphere part of the operator. */

90+ dp1: Float32Array;

91+ dq2: Float32Array;

92+ jinv: Float32Array;

8793 /** Mean-J preconditioner scale (Geometry.Jhat): folded into every

8894 * setParams upload as the 'jhat' uniform, so a .m that takes jhat is

8995 * never left with the zero a missing parameter would default to. An

@@ -99,7 +105,7 @@ export const GEOMETRY_SPECTRAL_NAMES = ['Gx', 'Gy', 'Gz'] as const;

99105 export const METRIC_GRID_NAMES = ['Vtx', 'Vty', 'Vtz', 'Vpx', 'Vpy', 'Vpz'] as const;

100106 /** Names the .m may take for the flux-form metric weights (six-transform

101107 * scheme). A model asks for whichever set its loop uses; both are uploaded. */

102-export const FLUX_METRIC_GRID_NAMES = ['p1', 'p2', 'q2', 'r'] as const;

108+export const FLUX_METRIC_GRID_NAMES = ['p1', 'p2', 'q2', 'r', 'dp1', 'dq2', 'jinv'] as const;

103109

104110 /** Laplace-Beltrami eigenvalues l(l+1), duplicated across re/im so the array

105111 * matches the 2 x nlm spectral layout element for element. */

@@ -270,6 +276,9 @@ export class GpuModel {

270276 host.upload('p2', geometry.p2);

271277 host.upload('q2', geometry.q2);

272278 host.upload('r', geometry.r);

279+ host.upload('dp1', geometry.dp1);

280+ host.upload('dq2', geometry.dq2);

281+ host.upload('jinv', geometry.jinv);

273282 }

274283

275284 const readback = device.createBuffer({

@@ -320,6 +329,7 @@ export class GpuModel {

320329 ['Vpx', geometry.Vpx], ['Vpy', geometry.Vpy], ['Vpz', geometry.Vpz],

321330 ['p1', geometry.p1], ['p2', geometry.p2],

322331 ['q2', geometry.q2], ['r', geometry.r],

332+ ['dp1', geometry.dp1], ['dq2', geometry.dq2], ['jinv', geometry.jinv],

323333 ];

324334 for (const [name, data] of fields) {

325335 if (this.#host.get(name)) this.#host.upload(name, data);

test/fluxChecks.tsmodified+77−4View file

@@ -30,6 +30,17 @@

3030 * non-axisymmetric surface, where the off-diagonal weight p2 actually does

3131 * something. The headline transform count (6 vs 12 per species per

3232 * iteration) is asserted from the compiled op sequences, not the doc.

33+ *

34+ * 3. The polar conditioning of the divergence (doc Sec 5). r ~ 1/sin^2(theta)

35+ * multiplies a bracket that must cancel to O(sin^2(theta)) at the poles,

36+ * so it amplifies the polar round-off of whatever it is handed. Splitting

37+ * the round sphere out of the divergence (models/schnakenberg.m, and

38+ * dp1/dq2/jinv in src/geom/geometry.ts) keeps r off all but the geometry

39+ * deviation; without the split, the amplified round-off is a static polar

40+ * forcing that a Turing instability grows into a spot at the pole,

41+ * regardless of the seed. That is the failure this checks for: it is

42+ * invisible to 1 and 2, which compare operators rather than watch what a

43+ * run nucleates from.

3344 */

3445 import { ShtPlan } from '../src/sht/sht.ts';

3546 import { DerivPlan } from '../src/sht/deriv.ts';

@@ -368,13 +379,15 @@ export async function fluxChecks(

368379 xformsPerIter.push(counts[1] - counts[0]);

369380 }

370381

371- // The headline number, from the compiled op sequences: 5 Legendre

382+ // The headline number, from the compiled op sequences: 6 Legendre

372383 // transforms per species per iteration against Algorithm 4's 12

373- // (2 species here). The phi flux's derivative runs as dphig -- two

384+ // (2 species here). Five of the six are the flux matvec; the sixth is

385+ // the round-sphere synthesis the divergence split buys its polar

386+ // conditioning with. The phi flux's derivative runs as dphig -- two

374387 // Fourier stages, no Legendre work -- and is deliberately not counted.

375388 check(

376- 'flux: 5 Legendre transforms per species per iteration, versus 12',

377- xformsPerIter[0] === 10 && xformsPerIter[1] === 24,

389+ 'flux: 6 Legendre transforms per species per iteration, versus 12',

390+ xformsPerIter[0] === 12 && xformsPerIter[1] === 24,

378391 `flux form adds ${xformsPerIter[0]} transforms/iteration, ` +

379392 `Algorithm 4 adds ${xformsPerIter[1]}`,

380393 );

@@ -399,5 +412,65 @@ export async function fluxChecks(

399412 `max |U_flux - U_alg4| = ${worst.toExponential(2)} after ${STEPS} steps ` +

400413 `on bumpy at lmax ${LMAX_AB}`,

401414 );

415+

416+ // ---- 3. the correction must not manufacture its own perturbation -----

417+ //

418+ // From the exact uniform steady state, with the Turing band switched off

419+ // (D1 = D2) so nothing can grow on its own, the only thing driving the

420+ // state away from uniform is round-off. niter = 0 never touches the flux

421+ // machinery and sets the floor; niter = 3 runs it three times per step.

422+ // The ratio is the correction's noise gain. Sphere-split it is O(1); with

423+ // r multiplying the whole divergence it was ~50 at lmax 63, and that

424+ // margin is what decides where a pattern nucleates. The ellipsoid is the

425+ // case to run it on: axisymmetric grid, strongly non-spherical geometry.

426+ const model = mModelByKey('schnakenberg')!;

427+ const quiet = model.source.replace(

428+ /function \[U, V, u, v\] = init\([\s\S]*?\nend/,

429+ `function [U, V, u, v] = init(lam3, gx, gy, gz, a, b)

430+ us = a + b;

431+ vs = b / (us * us);

432+ [U, V] = analys(us * ones(numel(gx), 1), vs * ones(numel(gx), 1));

433+ [u, v] = synth(U, V);

434+end`,

435+ );

436+ if (quiet === model.source) throw new Error('quiet-start fixture no longer matches schnakenberg.m');

437+ const ell = mGeometryByKey('ellipsoid')!;

438+ const noise: number[] = [];

439+ for (const niter of [0, 3]) {

440+ const session = await ModelSession.create({

441+ device,

442+ model,

443+ params: { ...defaultParams(model), D2: defaultParams(model).D1 },

444+ lmax: LMAX_AB,

445+ source: quiet,

446+ niter,

447+ geometry: ell,

448+ geometryParams: defaultGeometryParams(ell),

449+ });

450+ await session.seed(1);

451+ session.step(400);

452+ const U = await session.read('U');

453+ // Everything above the mean: l = 0, m = 0 is the uniform state itself.

454+ let sum = 0;

455+ for (let m = 0; m <= LMAX_AB; m++) {

456+ for (let l = Math.max(m, 1); l <= LMAX_AB; l++) {

457+ const i = lmIndex(LMAX_AB, l, m);

458+ sum += (U[2 * i] ** 2 + U[2 * i + 1] ** 2) * (m === 0 ? 1 : 2);

459+ }

460+ }

461+ noise.push(Math.sqrt(sum));

462+ session.destroy();

463+ }

464+ const gain = noise[1] / noise[0];

465+ log(

466+ ` flux polar noise gain on ellipsoid: ||U'|| ${noise[0].toExponential(2)} ` +

467+ `at niter 0, ${noise[1].toExponential(2)} at niter 3`,

468+ );

469+ check(

470+ 'flux: the geometric correction does not amplify polar round-off',

471+ Number.isFinite(gain) && gain < 5,

472+ `niter-3 round-off is ${gain.toFixed(1)}x the niter-0 floor ` +

473+ `(sphere-split: ~1; r on the whole divergence: ~50)`,

474+ );

402475 }

403476 }

test/geometryChecks.tsmodified+7−6View file

@@ -310,14 +310,15 @@ 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: 3

313+ // once per iteration, no more and no less. Per species per iteration: 4

314314 // synths + 2 analyses (the flux-form matvec's five Legendre transforms,

315- // docs/reduced-transforms.md Sec 4 with the dphig variation) + the

316- // grid-space phi-derivative + 3 coefficient-space shuffles plus 7

317- // generated kernels -- see test/modelChecks.ts's KERNELS_PER_ITERATION,

318- // which counts the kernels alone; this counts every op.

315+ // docs/reduced-transforms.md Sec 4 with the dphig variation, plus the

316+ // round-sphere synthesis of the divergence split) + the grid-space

317+ // phi-derivative + 3 coefficient-space shuffles plus 8 generated kernels

318+ // -- see test/modelChecks.ts's KERNELS_PER_ITERATION, which counts the

319+ // kernels alone; this counts every op.

319320 const perIteration = ops[1] - ops[0];

320- const want = 32;

321+ const want = 36;

321322 check(

322323 'loop: unrolling is exactly linear in the trip count',

323324 perIteration === want && ops[2] - ops[0] === 4 * perIteration,

test/modelChecks.tsmodified+20−14View file

@@ -41,15 +41,19 @@ const EXPECTED_KERNELS: Record<string, number> = {

4141 * docs/reduced-transforms.md Sec 4: the two sin-weighted

4242 * derivative synths, the pointwise flux combination through p1/p2/q2, the

4343 * two flux analyses, the re-shifted divergence and its r-scaled synthesis,

44- * plus the round-sphere eigenvalue added back — see models/schnakenberg.m

44+ * the round-sphere share of the divergence subtracted off through jinv, plus

45+ * the round-sphere eigenvalue added back — see models/schnakenberg.m

4546 * and docs/richardson-iteration.md. `schnakenberg-alg4` keeps the original

4647 * Cartesian-gradient form (Algorithm 3/4 of evolving_surface/notes/algos.tex)

4748 * as a live reference, with its original counts.

4849 */

4950 const KERNELS_PER_ITERATION: Record<string, number> = {

50- schnakenberg: 14,

51- brusselator: 14,

52- allencahn: 7,

51+ // 14 / 14 / 7 before the divergence was split against the round sphere:

52+ // forming lam .* F for the sphere term, and subtracting jinv .* S from the

53+ // deviation's r-scaled divergence, is one extra kernel per species.

54+ schnakenberg: 16,

55+ brusselator: 16,

56+ allencahn: 8,

5357 // 30 before the correction gained its band projection (.* filt on dLu):

5458 // that line fused into the state update in this model's expression shape,

5559 // and no longer does — one extra 2 x nlm kernel per species per iteration.

@@ -178,15 +182,16 @@ export async function modelChecks(

178182 }

179183 // Every batchable run at one solve iteration: the u/v syntheses and the

180184 // reaction analyses outside the loop (2 + 2), the four gradient

181- // syntheses, two theta-flux analyses, two divergence syntheses and two

182- // final analyses inside it (4 + 2 + 2 + 2; the phi flux goes through

183- // dphig, which has no Legendre stage to batch). Lane counts are

184- // batch-width invariant: a x4 run is one batch at K = 4 and two at

185- // K = 2, but the lanes annotated are the same 14 either way.

185+ // syntheses and the two round-sphere syntheses riding in the same group,

186+ // two theta-flux analyses, two divergence syntheses and two final

187+ // analyses inside it (6 + 2 + 2 + 2; the phi flux goes through dphig,

188+ // which has no Legendre stage to batch). Lane counts are batch-width

189+ // invariant: a x4 run is one batch at K = 4 and two at K = 2, but the

190+ // lanes annotated are the same 16 either way.

186191 check(

187192 'batch: the compiled step batches every adjacent transform pair',

188- batchedLanes === 14,

189- `${batchedLanes} batched transform lanes (expected 14)`,

193+ batchedLanes === 16,

194+ `${batchedLanes} batched transform lanes (expected 16)`,

190195 );

191196 let worst = 0;

192197 for (let i = 0; i < states[0].length; i++) {

@@ -208,7 +213,7 @@ export async function modelChecks(

208213 const cases: [string, string, string][] = [

209214 [

210215 'a single output bound to a grouped call',

211- 'Ftu = synth(vtu, vpu);',

216+ 'Ftu = synth(vtu, vpu, lam .* Fu);',

212217 'bind each one',

213218 ],

214219 [

@@ -216,12 +221,13 @@ export async function modelChecks(

216221 // Fpu is reassigned so the only error left is the dropped slot

217222 // itself, which the planner refuses (numbl would otherwise catch

218223 // the undefined 'Fpu' first, masking the check under test).

219- '[Ftu, ~] = synth(vtu, vpu);\n Fpu = Ftu;',

224+ '[Ftu, ~, Su] = synth(vtu, vpu, lam .* Fu);\n Fpu = Ftu;',

220225 'must be bound',

221226 ],

222227 ];

223228 for (const [what, bad, expect] of cases) {

224- const source = model.source.replace('[Ftu, Fpu] = synth(vtu, vpu);', bad);

229+ const source = model.source.replace('[Ftu, Fpu, Su] = synth(vtu, vpu, lam .* Fu);', bad);

230+ if (source === model.source) throw new Error('grouped-call fixture no longer matches allencahn.m');

225231 let message = '';

226232 try {

227233 const session = await ModelSession.create({