Precondition with the operator's symbol; project the correction onto the band
Two fixes to the implicit solve's Richardson iteration, found from a
report that raising niter (8 -> 16) or lmax (63 -> 127) could make it
diverge, and that surviving runs drifted qualitatively higher-frequency.
1. The models precondition with lam/jhat, where jhat = 2/(muMin + muMax)
is the minimax constant over the eigenvalues mu(x) of the operator's
principal symbol -- the inverse squared principal stretches of the
embedding, read pointwise off the flux-metric arrays as eigenvalues of
S = (1/J)[[p1,p2],[p2,q2]]. This contracts every mode and every
direction at rate (muMax - muMin)/(muMax + muMin) < 1 on any surface;
the plain preconditioner (jhat = 1) diverges wherever mu > 2, which is
what put peanut (mu up to 6.2) outside the convergence radius at
niter >= 2, and a det-based mean of the area factor J (mu's geometric
mean, exact only conformally) still left directional high-degree bands
above 1 on the ellipsoid. jhat is a per-geometry uniform folded into
every params upload -- geometry swaps update it with no recompile, the
fixed point never depends on it, and jhat: 1 pins the old behavior for
A/B. The whole niter x geometry sweep now converges; KNOWN_DIVERGENT
is empty, with a jhat=1 control asserting the divergence it retired.
2. The correction is projected onto the band: dLu = (... ) .* filt,
algos.tex Algorithm 5's zeroing of the l >= L-2 coefficients, which
the port had dropped. Without it the top two degrees iterate toward
the *undiffused* right-hand side at species-dependent rates -- a
manufactured Turing band at the band edge, growing ~3%/step on the
round sphere at lmax 127, niter 8.
With both, the pattern's spectral centroid is resolution-independent
(l ~ 26 at lmax 63 and 127 alike, schnak-spots on ellipsoid) where it
previously climbed to the top of the band, and every formerly divergent
configuration runs: peanut at lmax 63/niter 16 and lmax 127/niter 8
included. On mildly deformed surfaces iteration counts do not drop --
one iteration already sits at the fp32 accumulation floor -- so the
payoff is reach and spectral correctness, recorded honestly in
docs/reduced-transforms.md Sec 10 and docs/richardson-iteration.md.
12 changed files+297−61
README.mdmodified+40−9View file
@@ -105,25 +105,56 @@ solve (see [docs/richardson-iteration.md](docs/richardson-iteration.md)). One
105105 species of [`models/schnakenberg.m`](models/schnakenberg.m)'s solve loop:
106106
107107 ```matlab
108+lamJ = lam ./ jhat; % mean-J preconditioner eigenvalues (below)
109+...
108110 for k = 1:niter
109111 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+ vtu = dthetac(Fu);
113+ vpu = dphic(Fu);
114+ [Ftu, Fpu] = synth(vtu, vpu); % sin(theta)*dtheta(u), dphi(u) -- smooth on
115+ % the sphere, one batched dispatch
112116 Pu = p1 .* Ftu + p2 .* Fpu; % the two fluxes, also smooth: the precomputed
113117 Qu = p2 .* Ftu + q2 .* Fpu; % weights carry every 1/sin(theta) there is
114- Pcu = analys(Pu) .* filt;
115- Qcu = analys(Qu) .* filt;
118+ [PAu, QAu] = analys(Pu, Qu);
119+ Pcu = PAu .* filt;
120+ Qcu = QAu .* filt;
116121 scu = dthetac(Pcu) + dphic(Qcu); % divergence, in coefficient space
117122 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);
123+ dLu = (analys(lapu) + lamJ .* Un) .* filt; % dlap, projected onto the band
124+ Un = (Bu + (dt * D1) * dLu) ./ (1 + (dt * D1) * lamJ);
120125 end
121126 ```
122127
123128 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.
129+`r = 1/sin²θ`, `jhat = 1`, and the composition collapses to `lap_s` — so the
130+sphere case reproduces turing-sphere to fp32 round-off, and the tests assert
131+the state stays put across 0, 1 and 4 iterations.
132+
133+**The preconditioner folds in the symbol of the operator.** `jhat` is the
134+host's minimax scale `2/(μmin + μmax)` over the eigenvalues `μ(x)` of the
135+operator's principal symbol — the inverse squared principal stretches of
136+the embedding, direction included, read straight off the flux-metric
137+arrays (`S = (1/J)·[[p1,p2],[p2,q2]]`). Preconditioning with `lam/jhat`
138+then contracts every mode *and every direction* at rate
139+`(μmax − μmin)/(μmax + μmin) < 1` on any surface, where the plain `lam`
140+diverges wherever `μ > 2` — peanut reaches `μ = 6.2`. A det-based mean of
141+the area factor (μ's geometric mean, exact only for conformal surfaces) is
142+not enough: it under-corrects anisotropic stretching and leaves directional
143+high-degree bands with amplification > 1, which surfaced as patterns going
144+high-frequency and diverging as `niter` or `lmax` grew. The answer never
145+depends on `jhat` — the `lamJ` term added inside `dLu` is the term divided
146+back out — only the convergence rate does.
147+
148+**The correction is projected onto the band** (`.* filt` on `dLu`,
149+matching algos.tex Algorithm 5's zeroing of the top coefficients). Without
150+it the top two degrees iterate toward the *undiffused* `Bu` — each solve
151+iteration strips a bit more of their implicit diffusion, at species-
152+dependent rates, which manufactures a spurious Turing band at the band
153+edge: visible on the round sphere as top-degree energy growing ~3%/step at
154+`lmax 127, niter 8`. With both fixes the whole niter × geometry sweep
155+converges, the spectral centroid of the pattern is resolution-independent
156+(l ≈ 26 at lmax 63 and 127 alike), and `jhat: 1` is kept as the divergent
157+control in the tests.
127158
128159 ### The geometry in the operator
129160
docs/reduced-transforms.mdmodified+29−0View file
@@ -386,3 +386,32 @@ and the hand-rolled $\alpha^\pm$ recurrences on this path are no longer needed.
386386 The other lever is iteration count rather than cost per iteration. Since $M^{-1}A$ approaches
387387 multiplication by $1/J$ at high $\ell$, folding a mean or smoothed $J$ into the preconditioner could
388388 reduce GMRES iterations by more than any of the above reduces transforms.
389+
390+### Measured (2026-08-05, Richardson iteration, fp32)
391+
392+Implemented, with two corrections the measurements forced.
393+
394+**The right constant is the minimax over the symbol, not over $J$.** The high-$\ell$ per-mode
395+factor is governed by the full principal symbol: in the orthonormal frame the symbol matrix is
396+$S = (1/J)\begin{pmatrix} p_1 & p_2 \\ p_2 & q_2\end{pmatrix}$, whose eigenvalues $\mu(x)$ are the
397+inverse squared principal stretches — direction matters. Preconditioning with $\lambda/\hat J$
398+contracts every mode and direction iff $\hat J\mu \in (0,2)$, so
399+$$\hat J = 2/(\mu_{\min} + \mu_{\max}), \qquad \text{rate} = (\mu_{\max}-\mu_{\min})/(\mu_{\max}+\mu_{\min}) < 1.$$
400+The det-based mean of $J$ (this section's original suggestion; $\mu$'s geometric mean, exact only
401+for conformal surfaces) is insufficient: on the shipped ellipsoid it leaves directional
402+high-degree bands with amplification $> 1$ — patterns went qualitatively high-frequency at
403+moderate settings and diverged as niter or $L$ grew. With the symbol-based constant every
404+niter/geometry combination in the test sweep converges (peanut: $\mu \in [0.44, 6.2]$, plain rate
405+5.2, preconditioned rate 0.87).
406+
407+**The correction must be band-projected.** Algorithm 5's "zero $\ell \ge L-2$" is load-bearing:
408+without applying the same mask to the correction $d\Delta u$, the top two degrees iterate toward
409+the *undiffused* right-hand side — each Richardson iteration strips more of their implicit
410+diffusion, at species-dependent rates, manufacturing a spurious Turing band at the band edge
411+(observed on the round sphere: top-degree energy growing $\sim 3\%$/step at $L=127$, 8 iterations).
412+
413+**Payoff shape:** on mildly deformed surfaces one iteration already reaches the $\sim10^{-4}$ fp32
414+accumulation floor, so iteration counts do not drop — the speculation above does not hold at fp32.
415+The gain is reach and correctness: stiff geometries and high niter/$L$ combinations that
416+previously diverged (or silently shifted the pattern's wavelength) now converge with a
417+resolution-independent spectrum.
docs/richardson-iteration.mdmodified+17−0View file
@@ -86,6 +86,23 @@ diverges over many steps — rather than being caught the way algos.tex's
8686 `solve_step` catches it (its `info != 0` return, logged when GMRES fails to
8787 reach `tol` within `maxiter`).
8888
89+**Update (symbol-based preconditioning).** The models now precondition
90+with `M = I + dt*D*lam/jhat` where `jhat = 2/(muMin + muMax)` is the host's
91+minimax scale over the eigenvalues of the operator's principal symbol — the
92+inverse squared principal stretches of the embedding, direction included
93+(see docs/reduced-transforms.md Sec 10). At high degree the iteration then
94+contracts at rate `(muMax - muMin)/(muMax + muMin) < 1` on any surface,
95+where the plain `M` diverges wherever `mu > 2` — which is what used to put
96+peanut outside the convergence radius at niter >= 2, and what made patterns
97+drift high-frequency on the ellipsoid as niter or lmax grew. The correction
98+is also projected onto the band (`.* filt` on `dLu`, algos.tex Algorithm
99+5's zeroing), without which the top two degrees iterate toward an
100+undiffused fixed point. The silent-failure caveat above still stands for
101+what a constant scale cannot capture (strong *spatial* variation of the
102+symbol at low degree, or `dt*D` beyond the correction's reach), but the
103+sweep's previously divergent cases all converge now, and `jhat: 1`
104+reproduces the old behavior for A/B.
105+
89106 That tradeoff is deliberate, not an oversight, and it comes from where the
90107 two projects run. algos.tex's GMRES needs, every iteration: a dot product
91108 across the whole spectral state (Arnoldi orthogonalization) and a residual
models/allencahn.mmodified+7−4View file
@@ -9,11 +9,14 @@ function [U, u] = init(noise)
99 u = synth(U);
1010 end
1111
12-function [Un, u] = step(U, lam, filt, gx, gy, gz, p1, p2, q2, r, eps2, dt, niter)
12+function [Un, u] = step(U, lam, filt, gx, gy, gz, p1, p2, q2, r, jhat, eps2, dt, niter)
1313 u = synth(U);
1414
1515 Bu = U + dt * analys(u - u.^3);
16- Un = Bu ./ (1 + (dt * eps2) * lam);
16+
17+ % Mean-J preconditioning -- see models/schnakenberg.m.
18+ lamJ = lam ./ jhat;
19+ Un = Bu ./ (1 + (dt * eps2) * lamJ);
1720
1821 for k = 1:niter
1922 % dlap = lap_g - lap_s, evaluated at the current iterate in flux form
@@ -31,8 +34,8 @@ function [Un, u] = step(U, lam, filt, gx, gy, gz, p1, p2, q2, r, eps2, dt, niter
3134 Qcu = QAu .* filt;
3235 scu = dthetac(Pcu) + dphic(Qcu);
3336 lapu = r .* synth(scu);
34- dLu = analys(lapu) + lam .* Un;
37+ dLu = (analys(lapu) + lamJ .* Un) .* filt;
3538
36- Un = (Bu + (dt * eps2) * dLu) ./ (1 + (dt * eps2) * lam);
39+ Un = (Bu + (dt * eps2) * dLu) ./ (1 + (dt * eps2) * lamJ);
3740 end
3841 end
models/brusselator.mmodified+9−7View file
@@ -11,7 +11,7 @@ function [U, V, u, v] = init(noise, A, B)
1111 [u, v] = synth(U, V);
1212 end
1313
14-function [Un, Vn, u, v] = step(U, V, lam, filt, gx, gy, gz, p1, p2, q2, r, A, B, D1, D2, dt, niter)
14+function [Un, Vn, u, v] = step(U, V, lam, filt, gx, gy, gz, p1, p2, q2, r, jhat, A, B, D1, D2, dt, niter)
1515 [u, v] = synth(U, V);
1616 uuv = u .* u .* v;
1717
@@ -21,8 +21,10 @@ function [Un, Vn, u, v] = step(U, V, lam, filt, gx, gy, gz, p1, p2, q2, r, A, B,
2121 Bu = U + dt * Ru;
2222 Bv = V + dt * Rv;
2323
24- Un = Bu ./ (1 + (dt * D1) * lam);
25- Vn = Bv ./ (1 + (dt * D2) * lam);
24+ % Mean-J preconditioning -- see models/schnakenberg.m.
25+ lamJ = lam ./ jhat;
26+ Un = Bu ./ (1 + (dt * D1) * lamJ);
27+ Vn = Bv ./ (1 + (dt * D2) * lamJ);
2628
2729 for k = 1:niter
2830 % dlap = lap_g - lap_s, evaluated at the current iterate in flux form
@@ -50,10 +52,10 @@ function [Un, Vn, u, v] = step(U, V, lam, filt, gx, gy, gz, p1, p2, q2, r, A, B,
5052 lapu = r .* Lu;
5153 lapv = r .* Lv;
5254 [LAu, LAv] = analys(lapu, lapv);
53- dLu = LAu + lam .* Un;
54- dLv = LAv + lam .* Vn;
55+ dLu = (LAu + lamJ .* Un) .* filt;
56+ dLv = (LAv + lamJ .* Vn) .* filt;
5557
56- Un = (Bu + (dt * D1) * dLu) ./ (1 + (dt * D1) * lam);
57- Vn = (Bv + (dt * D2) * dLv) ./ (1 + (dt * D2) * lam);
58+ Un = (Bu + (dt * D1) * dLu) ./ (1 + (dt * D1) * lamJ);
59+ Vn = (Bv + (dt * D2) * dLv) ./ (1 + (dt * D2) * lamJ);
5860 end
5961 end
models/schnakenberg.mmodified+26−11View file
@@ -22,7 +22,7 @@ function [U, V, u, v] = init(noise, a, b)
2222 [u, v] = synth(U, V);
2323 end
2424
25-function [Un, Vn, u, v] = step(U, V, lam, filt, gx, gy, gz, p1, p2, q2, r, a, b, D1, D2, dt, niter)
25+function [Un, Vn, u, v] = step(U, V, lam, filt, gx, gy, gz, p1, p2, q2, r, jhat, a, b, D1, D2, dt, niter)
2626 % Grouped transforms -- [a, b] = synth(x, y) -- are explicit batching:
2727 % output k is the transform of input k, and the whole group runs as one
2828 % batched Legendre dispatch, or as many as the device's lane width allows
@@ -38,9 +38,19 @@ function [Un, Vn, u, v] = step(U, V, lam, filt, gx, gy, gz, p1, p2, q2, r, a, b,
3838 Bu = U + dt * Ru;
3939 Bv = V + dt * Rv;
4040
41- % Round-sphere solve, then iterate the geometric correction.
42- Un = Bu ./ (1 + (dt * D1) * lam);
43- Vn = Bv ./ (1 + (dt * D2) * lam);
41+ % Preconditioned solve, then iterate the geometric correction. jhat is
42+ % the host's minimax scale over the operator's symbol eigenvalues mu(x)
43+ % -- the inverse squared principal stretches of the embedding, direction
44+ % included (src/geom/geometry.ts, Jhat): preconditioning with lam/jhat
45+ % contracts every mode and direction at rate
46+ % (muMax - muMin)/(muMax + muMin) < 1 on any surface, where the plain
47+ % lam (jhat = 1) diverges wherever mu > 2 -- docs/reduced-transforms.md
48+ % Sec 10. The answer never depends on jhat (the lamJ term added inside
49+ % dLu is the term divided back out); only the convergence rate does. On
50+ % the sphere mu = 1 and lamJ = lam.
51+ lamJ = lam ./ jhat;
52+ Un = Bu ./ (1 + (dt * D1) * lamJ);
53+ Vn = Bv ./ (1 + (dt * D2) * lamJ);
4454
4555 for k = 1:niter
4656 % dlap = lap_g - lap_s at the current iterate, in flux form
@@ -52,9 +62,14 @@ function [Un, Vn, u, v] = step(U, V, lam, filt, gx, gy, gz, p1, p2, q2, r, a, b,
5262 % through the *same* shuffles again and summed before the one synthesis
5363 % of the divergence, which r scales into lap_g(u). The only division by
5464 % sin(theta) anywhere is folded into p1,p2,q2,r at precompute time.
55- % lam.*Un adds back -lap_s(Un), since lam holds +l(l+1). filt zeroes the
56- % top two degrees, where the derivative recurrences cannot exactly
57- % represent a derivative.
65+ % lamJ.*Un adds back the preconditioner's -lap_s(Un)/jhat, since lam
66+ % holds +l(l+1). filt zeroes the top two degrees, where the derivative
67+ % recurrences cannot exactly represent a derivative -- and the correction
68+ % itself is projected onto the same band (algos.tex Algorithm 5 zeroes
69+ % the same coefficients): without that, each iteration replaces a bit
70+ % more of the top degrees' implicit diffusion with nothing (their fixed
71+ % point is the undiffused Bu), and the two species un-diffuse at
72+ % different rates -- a spurious Turing band at the band edge.
5873 %
5974 % The two species share each grouped call: the four gradient
6075 % syntheses, the four flux analyses, the two divergence syntheses and
@@ -81,10 +96,10 @@ function [Un, Vn, u, v] = step(U, V, lam, filt, gx, gy, gz, p1, p2, q2, r, a, b,
8196 lapu = r .* Lu;
8297 lapv = r .* Lv;
8398 [LAu, LAv] = analys(lapu, lapv);
84- dLu = LAu + lam .* Un;
85- dLv = LAv + lam .* Vn;
99+ dLu = (LAu + lamJ .* Un) .* filt;
100+ dLv = (LAv + lamJ .* Vn) .* filt;
86101
87- Un = (Bu + (dt * D1) * dLu) ./ (1 + (dt * D1) * lam);
88- Vn = (Bv + (dt * D2) * dLv) ./ (1 + (dt * D2) * lam);
102+ Un = (Bu + (dt * D1) * dLu) ./ (1 + (dt * D1) * lamJ);
103+ Vn = (Bv + (dt * D2) * dLv) ./ (1 + (dt * D2) * lamJ);
89104 end
90105 end
models/schnakenberg_alg4.mmodified+10−8View file
@@ -26,7 +26,7 @@ function [U, V, u, v] = init(noise, a, b)
2626 v = synth(V);
2727 end
2828
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)
29+function [Un, Vn, u, v] = step(U, V, lam, filt, gx, gy, gz, Vtx, Vty, Vtz, Vpx, Vpy, Vpz, jhat, a, b, D1, D2, dt, niter)
3030 u = synth(U);
3131 v = synth(V);
3232 uuv = u .* u .* v;
@@ -35,9 +35,11 @@ function [Un, Vn, u, v] = step(U, V, lam, filt, gx, gy, gz, Vtx, Vty, Vtz, Vpx,
3535 Bu = U + dt * analys(a - u + uuv);
3636 Bv = V + dt * analys(b - uuv);
3737
38- % Round-sphere solve, then iterate the geometric correction.
39- Un = Bu ./ (1 + (dt * D1) * lam);
40- Vn = Bv ./ (1 + (dt * D2) * lam);
38+ % Mean-J preconditioned solve (see models/schnakenberg.m), then iterate
39+ % the geometric correction.
40+ lamJ = lam ./ jhat;
41+ Un = Bu ./ (1 + (dt * D1) * lamJ);
42+ Vn = Bv ./ (1 + (dt * D2) * lamJ);
4143
4244 for k = 1:niter
4345 % dlap = lap_g - lap_s, evaluated at the current iterate (Algorithm 3 of
@@ -68,7 +70,7 @@ function [Un, Vn, u, v] = step(U, V, lam, filt, gx, gy, gz, Vtx, Vty, Vtz, Vpx,
6870 lapu = lapu + Fpcuy .* Vpy;
6971 lapu = lapu + Ftcuz .* Vtz;
7072 lapu = lapu + Fpcuz .* Vpz;
71- dLu = analys(lapu) + lam .* Un;
73+ dLu = (analys(lapu) + lamJ .* Un) .* filt;
7274
7375 Fv = Vn .* filt;
7476 Ftv = dtheta(Fv);
@@ -90,9 +92,9 @@ function [Un, Vn, u, v] = step(U, V, lam, filt, gx, gy, gz, Vtx, Vty, Vtz, Vpx,
9092 lapv = lapv + Fpcvy .* Vpy;
9193 lapv = lapv + Ftcvz .* Vtz;
9294 lapv = lapv + Fpcvz .* Vpz;
93- dLv = analys(lapv) + lam .* Vn;
95+ dLv = (analys(lapv) + lamJ .* Vn) .* filt;
9496
95- Un = (Bu + (dt * D1) * dLu) ./ (1 + (dt * D1) * lam);
96- Vn = (Bv + (dt * D2) * dLv) ./ (1 + (dt * D2) * lam);
97+ Un = (Bu + (dt * D1) * dLu) ./ (1 + (dt * D1) * lamJ);
98+ Vn = (Bv + (dt * D2) * dLv) ./ (1 + (dt * D2) * lamJ);
9799 end
98100 end
src/geom/geometry.tsmodified+65−0View file
@@ -85,6 +85,34 @@ export class Geometry {
8585 readonly p2: Float32Array;
8686 readonly q2: Float32Array;
8787 readonly r: Float32Array;
88+ /**
89+ * Preconditioner scale for the implicit solve (docs/reduced-transforms.md
90+ * Sec 10). At high degree the Richardson iteration's per-mode factor is
91+ * governed by the operator's principal symbol: in the orthonormal frame
92+ * the surface symbol matrix is S = (1/J)[[p1, p2], [p2, q2]], whose
93+ * eigenvalues mu(x) are the inverse squared principal stretches of the
94+ * embedding — the round sphere has mu = 1. Preconditioning with lam/Jhat
95+ * contracts every mode and every direction iff Jhat*mu stays in (0, 2),
96+ * so the minimax constant is the harmonic mean of the symbol extremes,
97+ *
98+ * Jhat = 2/(muMin + muMax), rate = (muMax - muMin)/(muMax + muMin) < 1.
99+ *
100+ * The direction dependence is the point: a det-based mean of the area
101+ * factor J (mu's geometric mean, exact only for conformal surfaces)
102+ * under-corrects anisotropic stretching — on the shipped ellipsoid it
103+ * leaves a band of directional high-degree modes with amplification > 1,
104+ * which inflates the pattern's spectrum at moderate niter/lmax and
105+ * diverges at larger ones. The plain scheme (Jhat = 1) diverges wherever
106+ * muMax > 2. The solve's fixed point never depends on Jhat; only the
107+ * convergence rate does.
108+ */
109+ readonly Jhat: number;
110+ /** Symbol-eigenvalue range over the grid (see Jhat), for diagnostics. */
111+ readonly muMin: number;
112+ readonly muMax: number;
113+ /** Area-factor range over the grid, for diagnostics. */
114+ readonly Jmin: number;
115+ readonly Jmax: number;
88116
89117 private constructor(init: {
90118 x: Float32Array; y: Float32Array; z: Float32Array;
@@ -92,6 +120,7 @@ export class Geometry {
92120 Vtx: Float32Array; Vty: Float32Array; Vtz: Float32Array;
93121 Vpx: Float32Array; Vpy: Float32Array; Vpz: Float32Array;
94122 p1: Float32Array; p2: Float32Array; q2: Float32Array; r: Float32Array;
123+ Jhat: number; muMin: number; muMax: number; Jmin: number; Jmax: number;
95124 }) {
96125 this.x = init.x;
97126 this.y = init.y;
@@ -109,6 +138,11 @@ export class Geometry {
109138 this.p2 = init.p2;
110139 this.q2 = init.q2;
111140 this.r = init.r;
141+ this.Jhat = init.Jhat;
142+ this.muMin = init.muMin;
143+ this.muMax = init.muMax;
144+ this.Jmin = init.Jmin;
145+ this.Jmax = init.Jmax;
112146 }
113147
114148 /**
@@ -189,12 +223,43 @@ export class Geometry {
189223 const sXtz = await deriv.sinDtheta(Z);
190224 const flux = computeFluxMetric(npts, sXtx, sXty, sXtz, Xp, Yp, Zp);
191225
226+ // The preconditioner scale — see the Jhat field comment. The symbol
227+ // matrix in the orthonormal frame is S = (1/J)[[p1,p2],[p2,q2]] with
228+ // 1/J = r sin^2(theta); its entries are the bounded quantities
229+ // g^tt, sin g^tp, sin^2 g^pp, so the eigenvalue extremes are clean to
230+ // take over the grid. det S = 1/J^2, so the area factor comes along
231+ // for free. f64 throughout.
232+ let muMin = Infinity;
233+ let muMax = 0;
234+ let Jmin = Infinity;
235+ let Jmax = 0;
236+ for (let i = 0; i < cfg.nlat; i++) {
237+ const ct = sht.cosTheta[i];
238+ const st2 = Math.max(0, 1 - ct * ct);
239+ for (let j = 0; j < cfg.nphi; j++) {
240+ const k = i * cfg.nphi + j;
241+ const invJ = flux.r[k] * st2;
242+ const s11 = flux.p1[k] * invJ;
243+ const s12 = flux.p2[k] * invJ;
244+ const s22 = flux.q2[k] * invJ;
245+ const mean = (s11 + s22) / 2;
246+ const disc = Math.sqrt(((s11 - s22) / 2) ** 2 + s12 * s12);
247+ if (mean - disc < muMin) muMin = mean - disc;
248+ if (mean + disc > muMax) muMax = mean + disc;
249+ const J = 1 / invJ;
250+ if (J < Jmin) Jmin = J;
251+ if (J > Jmax) Jmax = J;
252+ }
253+ }
254+ const Jhat = 2 / (muMin + muMax);
255+
192256 return new Geometry({
193257 x, y, z, X, Y, Z, Vtx, Vty, Vtz, Vpx, Vpy, Vpz,
194258 p1: new Float32Array(flux.p1),
195259 p2: new Float32Array(flux.p2),
196260 q2: new Float32Array(flux.q2),
197261 r: new Float32Array(flux.r),
262+ Jhat, muMin, muMax, Jmin, Jmax,
198263 });
199264 } finally {
200265 plan.destroy();
src/mgpu/model.tsmodified+21−3View file
@@ -83,6 +83,12 @@ export interface GeometryBuffers {
8383 p2: Float32Array;
8484 q2: Float32Array;
8585 r: Float32Array;
86+ /** Mean-J preconditioner scale (Geometry.Jhat): folded into every
87+ * setParams upload as the 'jhat' uniform, so a .m that takes jhat is
88+ * never left with the zero a missing parameter would default to. An
89+ * explicit jhat in the params wins (jhat: 1 pins the plain round-sphere
90+ * preconditioner, for A/B). */
91+ Jhat: number;
8692 }
8793
8894 /** Names the .m may take for the grid coordinates and for their coefficients. */
@@ -140,6 +146,8 @@ export class GpuModel {
140146 #host: HostBuffers;
141147 #initPlan: ModelPlan;
142148 #stepPlan: ModelPlan;
149+ /** Current geometry's mean-J scale; 1 with no geometry (the sphere). */
150+ #jhat = 1;
143151 #readback: GPUBuffer;
144152 /** Scratch holding a copy of the whole spectral state; see snapshotState. */
145153 #stash: GPUBuffer;
@@ -196,6 +204,10 @@ export class GpuModel {
196204 for (const g of GEOMETRY_SPECTRAL_NAMES) bindings[g] = { kind: 'tensor', shape: [2, nlm] };
197205 for (const g of METRIC_GRID_NAMES) bindings[g] = { kind: 'tensor', shape: [npts, 1] };
198206 for (const g of FLUX_METRIC_GRID_NAMES) bindings[g] = { kind: 'tensor', shape: [npts, 1] };
207+ // Mean-J preconditioner scale (Geometry.Jhat): a uniform, not a const,
208+ // so swapping the surface updates it with no recompile. The session
209+ // folds the current geometry's value into every setParams call.
210+ bindings['jhat'] = { kind: 'param' };
199211 }
200212 for (const s of state) bindings[s] = { kind: 'tensor', shape: [2, nlm] };
201213 for (const p of paramNames) bindings[p] = { kind: 'param' };
@@ -266,15 +278,18 @@ export class GpuModel {
266278 usage: GPUBufferUsage.COPY_SRC | GPUBufferUsage.COPY_DST,
267279 });
268280
269- return new GpuModel({
281+ const gpu = new GpuModel({
270282 device, host, initPlan, stepPlan, readback, stash,
271283 paramNames, state, view, npts, nlm,
272284 });
285+ if (geometry) gpu.#jhat = geometry.Jhat;
286+ return gpu;
273287 }
274288
275289 setParams(params: ModelParams): void {
276- this.#initPlan.setParams(params);
277- this.#stepPlan.setParams(params);
290+ const merged = { jhat: this.#jhat, ...params };
291+ this.#initPlan.setParams(merged);
292+ this.#stepPlan.setParams(merged);
278293 }
279294
280295 /**
@@ -304,6 +319,9 @@ export class GpuModel {
304319 for (const [name, data] of fields) {
305320 if (this.#host.get(name)) this.#host.upload(name, data);
306321 }
322+ // The new surface's preconditioner scale takes effect on the next
323+ // setParams (the session re-applies its params after a swap).
324+ this.#jhat = geometry.Jhat;
307325 }
308326
309327 /** Upload the seeded perturbation and run `init`. */
src/mgpu/session.tsmodified+3−0View file
@@ -205,6 +205,9 @@ export class ModelSession {
205205 this.#geometry = next;
206206 this.#geometryModel = geometryModel;
207207 this.gpu.uploadGeometry(next);
208+ // The new surface brings a new preconditioner scale (GpuModel folds its
209+ // current geometry's jhat into every params upload).
210+ this.gpu.setParams(this.#params);
208211 }
209212
210213 /** The plan whose grid `readSpecies` samples on — the display plan when
test/geometryChecks.tsmodified+62−14View file
@@ -379,20 +379,19 @@ export async function geometryChecks(
379379
380380 // ---- niter x geometry sweep: catch a "doesn't run" regression early -----
381381 // This is what actually turned up the two real issues found while building
382- // the correction: peanut diverging at niter >= 4 with schnak-spots'
383- // shipped default dt (a genuine Richardson-convergence-radius limit, not a
384- // bug -- see docs/richardson-iteration.md), and a since-fixed compiler bug
385- // where a loop-body statement could silently reuse a *different*
386- // statement's compiled kernel (test/modelChecks.ts's pipeline-cache check
387- // guards that one directly). Every shipped geometry x every niter the
388- // app's <select> actually offers, so a regression anywhere in that grid is
389- // caught -- without asserting away the one combination already known to be
390- // outside the convergence radius.
382+ // the correction: peanut diverging at niter >= 2 with schnak-spots'
383+ // shipped default dt (the plain round-sphere preconditioner's convergence
384+ // radius -- the mean-J preconditioner has since lifted it; see the control
385+ // check after the sweep), and a since-fixed compiler bug where a loop-body
386+ // statement could silently reuse a *different* statement's compiled kernel
387+ // (test/modelChecks.ts's pipeline-cache check guards that one directly).
388+ // Every shipped geometry x every niter the app's <select> actually offers,
389+ // so a regression anywhere in that grid is caught.
391390 //
392- // SWEEP_LMAX cannot be lowered to make this cheaper: at lmax 31 or 15 the
393- // peanut/2, /4 and /8 combinations below all stay finite, so a smaller grid
394- // would quietly turn KNOWN_DIVERGENT into three failures and stop testing
395- // the thing this sweep exists to pin down.
391+ // SWEEP_LMAX stays at the app's default: the divergence the control check
392+ // pins down is lmax-dependent (at lmax 31 or 15 even the plain
393+ // preconditioner stays finite on peanut), so a smaller grid would stop
394+ // testing the thing the control exists to demonstrate.
396395 if (!runSweep) {
397396 log(
398397 ' sweep: skipped — run `npm run test:node` (desktop Dawn, ~3 s) or ' +
@@ -402,7 +401,13 @@ export async function geometryChecks(
402401 const model = mModelByKey('schnakenberg')!;
403402 const params = defaultParams(model);
404403 const SWEEP_NITER = [0, 1, 2, 4, 8];
405- const KNOWN_DIVERGENT = new Set(['peanut/2', 'peanut/4', 'peanut/8']);
404+ // Empty since the symbol-based preconditioner: its high-degree
405+ // contraction rate (muMax - muMin)/(muMax + muMin) < 1 on any surface
406+ // (mu = the symbol eigenvalues, i.e. inverse squared principal
407+ // stretches), where the plain preconditioner diverges wherever mu > 2
408+ // -- which is exactly what used to make peanut/2, /4 and /8 diverge.
409+ // The mechanism stays: a regression lands here with its evidence.
410+ const KNOWN_DIVERGENT = new Set<string>([]);
406411
407412 for (const geomSpec of mGeometries) {
408413 for (const niter of SWEEP_NITER) {
@@ -434,6 +439,49 @@ export async function geometryChecks(
434439 );
435440 }
436441 }
442+
443+ // ---- the mean-J control: what the sweep's health is owed to ----------
444+ // peanut at niter 4 was the canonical divergent case before the mean-J
445+ // preconditioner. Pinning jhat to 1 reproduces the plain round-sphere
446+ // preconditioner on today's code, so this asserts both directions at
447+ // once: mean-J converges where plain diverges, on the same operator,
448+ // same surface, same dt. If this check ever finds jhat = 1 finite, the
449+ // sweep above has stopped exercising the regime the preconditioner
450+ // exists for (e.g. someone lowered SWEEP_LMAX or dt).
451+ {
452+ const peanut = mGeometryByKey('peanut')!;
453+ const outcomes: boolean[] = [];
454+ let jstats = '';
455+ for (const jhat of [undefined, 1]) {
456+ const session = await ModelSession.create({
457+ device, model,
458+ params: jhat === undefined ? params : { ...params, jhat },
459+ lmax: SWEEP_LMAX,
460+ geometry: peanut, geometryParams: defaultGeometryParams(peanut),
461+ niter: 4,
462+ });
463+ if (jhat === undefined) {
464+ const g = session.geometry;
465+ jstats =
466+ `mu in [${g.muMin.toFixed(3)}, ${g.muMax.toFixed(3)}] ` +
467+ `(J in [${g.Jmin.toFixed(3)}, ${g.Jmax.toFixed(3)}]), ` +
468+ `Jhat ${g.Jhat.toFixed(3)}, ` +
469+ `rate ${((g.muMax - g.muMin) / (g.muMax + g.muMin)).toFixed(3)} ` +
470+ `vs plain ${(g.muMax - 1).toFixed(2)}`;
471+ }
472+ session.seed(1);
473+ session.step(STEPS);
474+ const values = await session.read('u');
475+ outcomes.push(values.every((v) => Number.isFinite(v)));
476+ session.destroy();
477+ }
478+ log(` mean-J on peanut: ${jstats}`);
479+ check(
480+ 'mean-J: converges on peanut/4 where the plain preconditioner diverges',
481+ outcomes[0] && !outcomes[1],
482+ `mean-J finite: ${outcomes[0]}, jhat=1 finite: ${outcomes[1]}`,
483+ );
484+ }
437485 }
438486
439487 // ---- a loop whose length is not known at compile time is refused --------
test/modelChecks.tsmodified+8−5View file
@@ -25,10 +25,10 @@ import type { Check, Log } from './analyticChecks.ts';
2525 * (it cannot fuse into an external call).
2626 */
2727 const EXPECTED_KERNELS: Record<string, number> = {
28- schnakenberg: 7,
29- brusselator: 7,
30- allencahn: 3,
31- 'schnakenberg-alg4': 7,
28+ schnakenberg: 8,
29+ brusselator: 8,
30+ allencahn: 4,
31+ 'schnakenberg-alg4': 8,
3232 };
3333
3434 /**
@@ -50,7 +50,10 @@ const KERNELS_PER_ITERATION: Record<string, number> = {
5050 schnakenberg: 18,
5151 brusselator: 18,
5252 allencahn: 9,
53- 'schnakenberg-alg4': 30,
53+ // 30 before the correction gained its band projection (.* filt on dLu):
54+ // that line fused into the state update in this model's expression shape,
55+ // and no longer does — one extra 2 x nlm kernel per species per iteration.
56+ 'schnakenberg-alg4': 32,
5457 };
5558
5659 const LMAX = 31;