WIP: added code for Laplace-Beltrami operator evaluation on smooth genus-0 surface
19 changed files+1569−179
models/allencahn.mmodified+26−3View file
@@ -9,15 +9,38 @@ function [U, u] = init(noise)
99 u = synth(U);
1010 end
1111
12-function [Un, u] = step(U, lam, gx, gy, gz, eps2, dt, niter)
12+function [Un, u] = step(U, lam, filt, gx, gy, gz, Vtx, Vty, Vtz, Vpx, Vpy, Vpz, 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- % Placeholder: dlap = lap_g - lap_s (see models/schnakenberg.m).
20- dLu = 0 * Un;
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).
22+ 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;
42+ dLu = analys(lapu) + lam .* Un;
43+
2144 Un = (Bu + (dt * eps2) * dLu) ./ (1 + (dt * eps2) * lam);
2245 end
2346 end
models/brusselator.mmodified+48−4View 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, gx, gy, gz, A, B, D1, D2, dt, niter)
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)
1616 u = synth(U);
1717 v = synth(V);
1818 uuv = u .* u .* v;
@@ -24,9 +24,53 @@ function [Un, Vn, u, v] = step(U, V, lam, gx, gy, gz, A, B, D1, D2, dt, niter)
2424 Vn = Bv ./ (1 + (dt * D2) * lam);
2525
2626 for k = 1:niter
27- % Placeholder: dlap = lap_g - lap_s (see models/schnakenberg.m).
28- dLu = 0 * Un;
29- dLv = 0 * Vn;
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).
30+ 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;
50+ dLu = analys(lapu) + lam .* Un;
51+
52+ 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;
72+ dLv = analys(lapv) + lam .* Vn;
73+
3074 Un = (Bu + (dt * D1) * dLu) ./ (1 + (dt * D1) * lam);
3175 Vn = (Bv + (dt * D2) * dLv) ./ (1 + (dt * D2) * lam);
3276 end
models/schnakenberg.mmodified+54−7View file
@@ -7,7 +7,7 @@
77 % splits lap_g = lap_s + dlap: the round-sphere part lap_s is diagonal in
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;
10-% spectral fields are real 2 x nlm.
10+% spectral fields are real 2 x nlm. See docs/richardson-iteration.md.
1111
1212 function [U, V, u, v] = init(noise, a, b)
1313 us = a + b;
@@ -18,7 +18,7 @@ function [U, V, u, v] = init(noise, a, b)
1818 v = synth(V);
1919 end
2020
21-function [Un, Vn, u, v] = step(U, V, lam, gx, gy, gz, a, b, D1, D2, dt, niter)
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)
2222 u = synth(U);
2323 v = synth(V);
2424 uuv = u .* u .* v;
@@ -32,11 +32,58 @@ function [Un, Vn, u, v] = step(U, V, lam, gx, gy, gz, a, b, D1, D2, dt, niter)
3232 Vn = Bv ./ (1 + (dt * D2) * lam);
3333
3434 for k = 1:niter
35- % Placeholder: dlap = lap_g - lap_s is still identically zero, so this
36- % is exactly the round-sphere solver and the geometry is only drawn.
37- % See the README.
38- dLu = 0 * Un;
39- dLv = 0 * Vn;
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.
43+ 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;
63+ dLu = analys(lapu) + lam .* Un;
64+
65+ 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;
85+ dLv = analys(lapv) + lam .* Vn;
86+
4087 Un = (Bu + (dt * D1) * dLu) ./ (1 + (dt * D1) * lam);
4188 Vn = (Bv + (dt * D2) * dLv) ./ (1 + (dt * D2) * lam);
4289 end
package-lock.jsonmodified+446−96View file
@@ -112,31 +112,6 @@
112112 "dev": true,
113113 "license": "Apache-2.0"
114114 },
115- "node_modules/@emnapi/core": {
116- "version": "2.0.0-alpha.3",
117- "resolved": "https://registry.npmjs.org/@emnapi/core/-/core-2.0.0-alpha.3.tgz",
118- "integrity": "sha512-AZypUeJ/yByuxyS7BlSNRDOMLMlROYtjYdIAuBmJssVz1UJDSeYxLrdizhXCFYhedC5bqd/ASy8EuNXbVVXp9g==",
119- "dev": true,
120- "license": "MIT",
121- "optional": true,
122- "peer": true,
123- "dependencies": {
124- "@emnapi/wasi-threads": "2.0.1",
125- "tslib": "^2.4.0"
126- }
127- },
128- "node_modules/@emnapi/runtime": {
129- "version": "2.0.0-alpha.3",
130- "resolved": "https://registry.npmjs.org/@emnapi/runtime/-/runtime-2.0.0-alpha.3.tgz",
131- "integrity": "sha512-hFPAhMUjJD9BSyCANEISPOogeXC9Zo9ZQl7L6vKnaVsMkCtzznaW/naYypeyl0Gv5rYfWYsZbpixTMpjDJzQeA==",
132- "dev": true,
133- "license": "MIT",
134- "optional": true,
135- "peer": true,
136- "dependencies": {
137- "tslib": "^2.4.0"
138- }
139- },
140115 "node_modules/@emnapi/wasi-threads": {
141116 "version": "2.0.1",
142117 "resolved": "https://registry.npmjs.org/@emnapi/wasi-threads/-/wasi-threads-2.0.1.tgz",
@@ -144,7 +119,6 @@
144119 "dev": true,
145120 "license": "MIT",
146121 "optional": true,
147- "peer": true,
148122 "dependencies": {
149123 "tslib": "^2.4.0"
150124 }
@@ -438,6 +412,23 @@
438412 "node": ">=12"
439413 }
440414 },
415+ "node_modules/@esbuild/netbsd-arm64": {
416+ "version": "0.28.1",
417+ "resolved": "https://registry.npmjs.org/@esbuild/netbsd-arm64/-/netbsd-arm64-0.28.1.tgz",
418+ "integrity": "sha512-oks0DYbLwWMmaakTsCb+zL4E+aHRVLom9IJZOAthMQEPiQmydXHkziYEsGYRx0uNV/IjEKGAV941JzH02pflqw==",
419+ "cpu": [
420+ "arm64"
421+ ],
422+ "dev": true,
423+ "license": "MIT",
424+ "optional": true,
425+ "os": [
426+ "netbsd"
427+ ],
428+ "engines": {
429+ "node": ">=18"
430+ }
431+ },
441432 "node_modules/@esbuild/netbsd-x64": {
442433 "version": "0.21.5",
443434 "resolved": "https://registry.npmjs.org/@esbuild/netbsd-x64/-/netbsd-x64-0.21.5.tgz",
@@ -455,6 +446,23 @@
455446 "node": ">=12"
456447 }
457448 },
449+ "node_modules/@esbuild/openbsd-arm64": {
450+ "version": "0.28.1",
451+ "resolved": "https://registry.npmjs.org/@esbuild/openbsd-arm64/-/openbsd-arm64-0.28.1.tgz",
452+ "integrity": "sha512-MEFJe5C3R8pwXdZ5Y21oo6m7ePiS0d9pWucn99O/wvyJZChoIQKrQDxKrGeW8F5+T0okTHesAmDeiHDTIq0V/Q==",
453+ "cpu": [
454+ "arm64"
455+ ],
456+ "dev": true,
457+ "license": "MIT",
458+ "optional": true,
459+ "os": [
460+ "openbsd"
461+ ],
462+ "engines": {
463+ "node": ">=18"
464+ }
465+ },
458466 "node_modules/@esbuild/openbsd-x64": {
459467 "version": "0.21.5",
460468 "resolved": "https://registry.npmjs.org/@esbuild/openbsd-x64/-/openbsd-x64-0.21.5.tgz",
@@ -472,6 +480,23 @@
472480 "node": ">=12"
473481 }
474482 },
483+ "node_modules/@esbuild/openharmony-arm64": {
484+ "version": "0.28.1",
485+ "resolved": "https://registry.npmjs.org/@esbuild/openharmony-arm64/-/openharmony-arm64-0.28.1.tgz",
486+ "integrity": "sha512-ge+Z7EXFNt2BO1oAMsVpiQ8EwndV9i1xXerAeTIK7AtPs3bKFXQM7nlRxDSIUIMeueR1CNXxqztLzdNeReKBJg==",
487+ "cpu": [
488+ "arm64"
489+ ],
490+ "dev": true,
491+ "license": "MIT",
492+ "optional": true,
493+ "os": [
494+ "openharmony"
495+ ],
496+ "engines": {
497+ "node": ">=18"
498+ }
499+ },
475500 "node_modules/@esbuild/sunos-x64": {
476501 "version": "0.21.5",
477502 "resolved": "https://registry.npmjs.org/@esbuild/sunos-x64/-/sunos-x64-0.21.5.tgz",
@@ -688,9 +713,6 @@
688713 "arm64"
689714 ],
690715 "dev": true,
691- "libc": [
692- "glibc"
693- ],
694716 "license": "MIT",
695717 "optional": true,
696718 "os": [
@@ -708,9 +730,6 @@
708730 "arm64"
709731 ],
710732 "dev": true,
711- "libc": [
712- "musl"
713- ],
714733 "license": "MIT",
715734 "optional": true,
716735 "os": [
@@ -728,9 +747,6 @@
728747 "ppc64"
729748 ],
730749 "dev": true,
731- "libc": [
732- "glibc"
733- ],
734750 "license": "MIT",
735751 "optional": true,
736752 "os": [
@@ -748,9 +764,6 @@
748764 "s390x"
749765 ],
750766 "dev": true,
751- "libc": [
752- "glibc"
753- ],
754767 "license": "MIT",
755768 "optional": true,
756769 "os": [
@@ -768,9 +781,6 @@
768781 "x64"
769782 ],
770783 "dev": true,
771- "libc": [
772- "glibc"
773- ],
774784 "license": "MIT",
775785 "optional": true,
776786 "os": [
@@ -788,9 +798,6 @@
788798 "x64"
789799 ],
790800 "dev": true,
791- "libc": [
792- "musl"
793- ],
794801 "license": "MIT",
795802 "optional": true,
796803 "os": [
@@ -1003,9 +1010,6 @@
10031010 "arm"
10041011 ],
10051012 "dev": true,
1006- "libc": [
1007- "glibc"
1008- ],
10091013 "license": "MIT",
10101014 "optional": true,
10111015 "os": [
@@ -1020,9 +1024,6 @@
10201024 "arm"
10211025 ],
10221026 "dev": true,
1023- "libc": [
1024- "musl"
1025- ],
10261027 "license": "MIT",
10271028 "optional": true,
10281029 "os": [
@@ -1037,9 +1038,6 @@
10371038 "arm64"
10381039 ],
10391040 "dev": true,
1040- "libc": [
1041- "glibc"
1042- ],
10431041 "license": "MIT",
10441042 "optional": true,
10451043 "os": [
@@ -1054,9 +1052,6 @@
10541052 "arm64"
10551053 ],
10561054 "dev": true,
1057- "libc": [
1058- "musl"
1059- ],
10601055 "license": "MIT",
10611056 "optional": true,
10621057 "os": [
@@ -1071,9 +1066,6 @@
10711066 "loong64"
10721067 ],
10731068 "dev": true,
1074- "libc": [
1075- "glibc"
1076- ],
10771069 "license": "MIT",
10781070 "optional": true,
10791071 "os": [
@@ -1088,9 +1080,6 @@
10881080 "loong64"
10891081 ],
10901082 "dev": true,
1091- "libc": [
1092- "musl"
1093- ],
10941083 "license": "MIT",
10951084 "optional": true,
10961085 "os": [
@@ -1105,9 +1094,6 @@
11051094 "ppc64"
11061095 ],
11071096 "dev": true,
1108- "libc": [
1109- "glibc"
1110- ],
11111097 "license": "MIT",
11121098 "optional": true,
11131099 "os": [
@@ -1122,9 +1108,6 @@
11221108 "ppc64"
11231109 ],
11241110 "dev": true,
1125- "libc": [
1126- "musl"
1127- ],
11281111 "license": "MIT",
11291112 "optional": true,
11301113 "os": [
@@ -1139,9 +1122,6 @@
11391122 "riscv64"
11401123 ],
11411124 "dev": true,
1142- "libc": [
1143- "glibc"
1144- ],
11451125 "license": "MIT",
11461126 "optional": true,
11471127 "os": [
@@ -1156,9 +1136,6 @@
11561136 "riscv64"
11571137 ],
11581138 "dev": true,
1159- "libc": [
1160- "musl"
1161- ],
11621139 "license": "MIT",
11631140 "optional": true,
11641141 "os": [
@@ -1173,9 +1150,6 @@
11731150 "s390x"
11741151 ],
11751152 "dev": true,
1176- "libc": [
1177- "glibc"
1178- ],
11791153 "license": "MIT",
11801154 "optional": true,
11811155 "os": [
@@ -1190,9 +1164,6 @@
11901164 "x64"
11911165 ],
11921166 "dev": true,
1193- "libc": [
1194- "glibc"
1195- ],
11961167 "license": "MIT",
11971168 "optional": true,
11981169 "os": [
@@ -1207,9 +1178,6 @@
12071178 "x64"
12081179 ],
12091180 "dev": true,
1210- "libc": [
1211- "musl"
1212- ],
12131181 "license": "MIT",
12141182 "optional": true,
12151183 "os": [
@@ -1733,7 +1701,8 @@
17331701 "resolved": "https://registry.npmjs.org/devtools-protocol/-/devtools-protocol-0.0.1367902.tgz",
17341702 "integrity": "sha512-XxtPuC3PGakY6PD7dG66/o8KwJ/LkH2/EKe19Dcw58w53dv4/vSQEkn/SzuyhHE2q4zPgCkxQBxus3VV4ql+Pg==",
17351703 "dev": true,
1736- "license": "BSD-3-Clause"
1704+ "license": "BSD-3-Clause",
1705+ "peer": true
17371706 },
17381707 "node_modules/emoji-regex": {
17391708 "version": "8.0.0",
@@ -2068,6 +2037,7 @@
20682037 "integrity": "sha512-WkUDrojuJs0xkgGf2udWxa3yGBRxPtxUkB79i6aCZLRgc7PM8fZe9TosfPDcvEpQZbuFASnHYmRLBLUbmLOIIA==",
20692038 "dev": true,
20702039 "license": "MPL-2.0",
2040+ "peer": true,
20712041 "dependencies": {
20722042 "detect-libc": "^2.0.3"
20732043 },
@@ -2205,9 +2175,6 @@
22052175 "arm64"
22062176 ],
22072177 "dev": true,
2208- "libc": [
2209- "glibc"
2210- ],
22112178 "license": "MPL-2.0",
22122179 "optional": true,
22132180 "os": [
@@ -2229,9 +2196,6 @@
22292196 "arm64"
22302197 ],
22312198 "dev": true,
2232- "libc": [
2233- "musl"
2234- ],
22352199 "license": "MPL-2.0",
22362200 "optional": true,
22372201 "os": [
@@ -2253,9 +2217,6 @@
22532217 "x64"
22542218 ],
22552219 "dev": true,
2256- "libc": [
2257- "glibc"
2258- ],
22592220 "license": "MPL-2.0",
22602221 "optional": true,
22612222 "os": [
@@ -2277,9 +2238,6 @@
22772238 "x64"
22782239 ],
22792240 "dev": true,
2280- "libc": [
2281- "musl"
2282- ],
22832241 "license": "MPL-2.0",
22842242 "optional": true,
22852243 "os": [
@@ -2495,6 +2453,7 @@
24952453 "integrity": "sha512-RvwwcruNjI1ncT5xRakeyS9Lf8lcItv34KD+aif+VH9kduAyfYBipGh12274xtenIPZ119/R9BdTBa8gAwSh0A==",
24962454 "dev": true,
24972455 "license": "MIT",
2456+ "peer": true,
24982457 "engines": {
24992458 "node": ">=12"
25002459 },
@@ -3008,6 +2967,397 @@
30082967 "url": "https://opencollective.com/antfu"
30092968 }
30102969 },
2970+ "node_modules/vite-node/node_modules/@esbuild/aix-ppc64": {
2971+ "version": "0.28.1",
2972+ "resolved": "https://registry.npmjs.org/@esbuild/aix-ppc64/-/aix-ppc64-0.28.1.tgz",
2973+ "integrity": "sha512-Svl7tq8k/08+p6CXPpRjQ1fKX+1odH/BQbb48fV6fj3CWHhsoIOoY87w1oHXm0qEpkIK3ZfVgp0hed3XBXzXMQ==",
2974+ "cpu": [
2975+ "ppc64"
2976+ ],
2977+ "dev": true,
2978+ "license": "MIT",
2979+ "optional": true,
2980+ "os": [
2981+ "aix"
2982+ ],
2983+ "engines": {
2984+ "node": ">=18"
2985+ }
2986+ },
2987+ "node_modules/vite-node/node_modules/@esbuild/android-arm": {
2988+ "version": "0.28.1",
2989+ "resolved": "https://registry.npmjs.org/@esbuild/android-arm/-/android-arm-0.28.1.tgz",
2990+ "integrity": "sha512-0k2F129Xdio1TdJfzJ8sy1Q47vUD2NnwdhiAf7drUN1EBTfPf4hsFCtmMgu/6m8JSzsBrlmVjudMBQqOfG8usQ==",
2991+ "cpu": [
2992+ "arm"
2993+ ],
2994+ "dev": true,
2995+ "license": "MIT",
2996+ "optional": true,
2997+ "os": [
2998+ "android"
2999+ ],
3000+ "engines": {
3001+ "node": ">=18"
3002+ }
3003+ },
3004+ "node_modules/vite-node/node_modules/@esbuild/android-arm64": {
3005+ "version": "0.28.1",
3006+ "resolved": "https://registry.npmjs.org/@esbuild/android-arm64/-/android-arm64-0.28.1.tgz",
3007+ "integrity": "sha512-34EGEbCIAgosYz6goLcopX6Mo7NyGv9tfwEM2/7Ce2VcVRk568iSvniGWcUXIy7wEDR1wzolcxcriFVrWYcwBg==",
3008+ "cpu": [
3009+ "arm64"
3010+ ],
3011+ "dev": true,
3012+ "license": "MIT",
3013+ "optional": true,
3014+ "os": [
3015+ "android"
3016+ ],
3017+ "engines": {
3018+ "node": ">=18"
3019+ }
3020+ },
3021+ "node_modules/vite-node/node_modules/@esbuild/android-x64": {
3022+ "version": "0.28.1",
3023+ "resolved": "https://registry.npmjs.org/@esbuild/android-x64/-/android-x64-0.28.1.tgz",
3024+ "integrity": "sha512-dbwY7ltSMDWsRatcRpCnES4F+im88OCUgGZjy52shC7GqHRE/cYlxNbB4Z4UpJswpcc4Qxd2oE/ufM0p61IKng==",
3025+ "cpu": [
3026+ "x64"
3027+ ],
3028+ "dev": true,
3029+ "license": "MIT",
3030+ "optional": true,
3031+ "os": [
3032+ "android"
3033+ ],
3034+ "engines": {
3035+ "node": ">=18"
3036+ }
3037+ },
3038+ "node_modules/vite-node/node_modules/@esbuild/darwin-arm64": {
3039+ "version": "0.28.1",
3040+ "resolved": "https://registry.npmjs.org/@esbuild/darwin-arm64/-/darwin-arm64-0.28.1.tgz",
3041+ "integrity": "sha512-TZbWkQY7kvTAXbXUT7uVACR5cMHsDiSz9z7ZKAX/RTq/WJEk3QyRr0wZpNhBDX+/0CtdqUIJlOiodQcta6tY3Q==",
3042+ "cpu": [
3043+ "arm64"
3044+ ],
3045+ "dev": true,
3046+ "license": "MIT",
3047+ "optional": true,
3048+ "os": [
3049+ "darwin"
3050+ ],
3051+ "engines": {
3052+ "node": ">=18"
3053+ }
3054+ },
3055+ "node_modules/vite-node/node_modules/@esbuild/darwin-x64": {
3056+ "version": "0.28.1",
3057+ "resolved": "https://registry.npmjs.org/@esbuild/darwin-x64/-/darwin-x64-0.28.1.tgz",
3058+ "integrity": "sha512-zfdzgK9ACBNZLI/CyHTOx81SyNbM6YXn7rxSgX97VjyiPl9W1i4Ka4fgKECEoFCKGpvBj5qArWIGgQjOwkgskQ==",
3059+ "cpu": [
3060+ "x64"
3061+ ],
3062+ "dev": true,
3063+ "license": "MIT",
3064+ "optional": true,
3065+ "os": [
3066+ "darwin"
3067+ ],
3068+ "engines": {
3069+ "node": ">=18"
3070+ }
3071+ },
3072+ "node_modules/vite-node/node_modules/@esbuild/freebsd-arm64": {
3073+ "version": "0.28.1",
3074+ "resolved": "https://registry.npmjs.org/@esbuild/freebsd-arm64/-/freebsd-arm64-0.28.1.tgz",
3075+ "integrity": "sha512-wG2EA8ENdEI0qhkSZMjfqrdY+ziCYCPMmtZjjIwOmXFjmyzEHn+UUxk5of+SYsjtfs3VpnlC7QLzSI5hY/rOAw==",
3076+ "cpu": [
3077+ "arm64"
3078+ ],
3079+ "dev": true,
3080+ "license": "MIT",
3081+ "optional": true,
3082+ "os": [
3083+ "freebsd"
3084+ ],
3085+ "engines": {
3086+ "node": ">=18"
3087+ }
3088+ },
3089+ "node_modules/vite-node/node_modules/@esbuild/freebsd-x64": {
3090+ "version": "0.28.1",
3091+ "resolved": "https://registry.npmjs.org/@esbuild/freebsd-x64/-/freebsd-x64-0.28.1.tgz",
3092+ "integrity": "sha512-i7dZ9vQgnvSCzi/rYCXNgtF/U+eKZNJBzu3eTQbRgHnM7tNSizLOkRFAl3qzVc/Op/u5YkHHa4pf/3DOYHthLQ==",
3093+ "cpu": [
3094+ "x64"
3095+ ],
3096+ "dev": true,
3097+ "license": "MIT",
3098+ "optional": true,
3099+ "os": [
3100+ "freebsd"
3101+ ],
3102+ "engines": {
3103+ "node": ">=18"
3104+ }
3105+ },
3106+ "node_modules/vite-node/node_modules/@esbuild/linux-arm": {
3107+ "version": "0.28.1",
3108+ "resolved": "https://registry.npmjs.org/@esbuild/linux-arm/-/linux-arm-0.28.1.tgz",
3109+ "integrity": "sha512-qVXBOHQS+d5Y722GwJzJUtOLlX7km3CraOaGormF1pDtPd2C/l1SHRPgjLunLGe51Sh5YYWKMFDyV4SxgMQYTQ==",
3110+ "cpu": [
3111+ "arm"
3112+ ],
3113+ "dev": true,
3114+ "license": "MIT",
3115+ "optional": true,
3116+ "os": [
3117+ "linux"
3118+ ],
3119+ "engines": {
3120+ "node": ">=18"
3121+ }
3122+ },
3123+ "node_modules/vite-node/node_modules/@esbuild/linux-arm64": {
3124+ "version": "0.28.1",
3125+ "resolved": "https://registry.npmjs.org/@esbuild/linux-arm64/-/linux-arm64-0.28.1.tgz",
3126+ "integrity": "sha512-yHs+0uc8+nvEAfAfxrWQKK5peSNzBc4PegcMO0EJ2hT71uA7vB8Ihg2e77R2P7SG5uYjPbHlLLmve4LLLRCf0g==",
3127+ "cpu": [
3128+ "arm64"
3129+ ],
3130+ "dev": true,
3131+ "license": "MIT",
3132+ "optional": true,
3133+ "os": [
3134+ "linux"
3135+ ],
3136+ "engines": {
3137+ "node": ">=18"
3138+ }
3139+ },
3140+ "node_modules/vite-node/node_modules/@esbuild/linux-ia32": {
3141+ "version": "0.28.1",
3142+ "resolved": "https://registry.npmjs.org/@esbuild/linux-ia32/-/linux-ia32-0.28.1.tgz",
3143+ "integrity": "sha512-d1z4ZuP0ajrfz/FhGT4vv278rX8KnPPJx8i5+AtK7TYbx9Le9F1hyzurZpkEyjkGa9dUGhQow4C1NmeGvqxN2w==",
3144+ "cpu": [
3145+ "ia32"
3146+ ],
3147+ "dev": true,
3148+ "license": "MIT",
3149+ "optional": true,
3150+ "os": [
3151+ "linux"
3152+ ],
3153+ "engines": {
3154+ "node": ">=18"
3155+ }
3156+ },
3157+ "node_modules/vite-node/node_modules/@esbuild/linux-loong64": {
3158+ "version": "0.28.1",
3159+ "resolved": "https://registry.npmjs.org/@esbuild/linux-loong64/-/linux-loong64-0.28.1.tgz",
3160+ "integrity": "sha512-M5sRjUVZrkm1OAPR3dlOYzNmN+loZKGVi1VUQGrwuqLcbR6qeAz+famMhjASeH3YVKvZz+zT1jlh/keC3Rj/lg==",
3161+ "cpu": [
3162+ "loong64"
3163+ ],
3164+ "dev": true,
3165+ "license": "MIT",
3166+ "optional": true,
3167+ "os": [
3168+ "linux"
3169+ ],
3170+ "engines": {
3171+ "node": ">=18"
3172+ }
3173+ },
3174+ "node_modules/vite-node/node_modules/@esbuild/linux-mips64el": {
3175+ "version": "0.28.1",
3176+ "resolved": "https://registry.npmjs.org/@esbuild/linux-mips64el/-/linux-mips64el-0.28.1.tgz",
3177+ "integrity": "sha512-mRObBZeHh2OxcBFPWE/FjylkRgZdYuiTR3vaTozquCGOH14iP9oN4x4Ge81CoIDYQrXmIxpFumJBu5MtZpnQJQ==",
3178+ "cpu": [
3179+ "mips64el"
3180+ ],
3181+ "dev": true,
3182+ "license": "MIT",
3183+ "optional": true,
3184+ "os": [
3185+ "linux"
3186+ ],
3187+ "engines": {
3188+ "node": ">=18"
3189+ }
3190+ },
3191+ "node_modules/vite-node/node_modules/@esbuild/linux-ppc64": {
3192+ "version": "0.28.1",
3193+ "resolved": "https://registry.npmjs.org/@esbuild/linux-ppc64/-/linux-ppc64-0.28.1.tgz",
3194+ "integrity": "sha512-slScBsMAb3GFDcdrCgLwZtPYRoH2H/youv10QiZyRjmsP48fznoveWytSgCI/R0ZcUgpc0ZhIUEx6LHts8yrfQ==",
3195+ "cpu": [
3196+ "ppc64"
3197+ ],
3198+ "dev": true,
3199+ "license": "MIT",
3200+ "optional": true,
3201+ "os": [
3202+ "linux"
3203+ ],
3204+ "engines": {
3205+ "node": ">=18"
3206+ }
3207+ },
3208+ "node_modules/vite-node/node_modules/@esbuild/linux-riscv64": {
3209+ "version": "0.28.1",
3210+ "resolved": "https://registry.npmjs.org/@esbuild/linux-riscv64/-/linux-riscv64-0.28.1.tgz",
3211+ "integrity": "sha512-kw0owk1o0GFETUJyW0jc0G4Yzs0BHZn0JDZ8JRT088vjJYX777BAs1fDGxAC+q831qOs2DTC96mNsG2opdfyyQ==",
3212+ "cpu": [
3213+ "riscv64"
3214+ ],
3215+ "dev": true,
3216+ "license": "MIT",
3217+ "optional": true,
3218+ "os": [
3219+ "linux"
3220+ ],
3221+ "engines": {
3222+ "node": ">=18"
3223+ }
3224+ },
3225+ "node_modules/vite-node/node_modules/@esbuild/linux-s390x": {
3226+ "version": "0.28.1",
3227+ "resolved": "https://registry.npmjs.org/@esbuild/linux-s390x/-/linux-s390x-0.28.1.tgz",
3228+ "integrity": "sha512-/lAIjX8aYFRByhh6L5rYtPEDRqa9de/4V/juOXcta5frjvzXO4/sqEtyytse0g3zZFuWu5cDN0MkLz2qRDD2Ag==",
3229+ "cpu": [
3230+ "s390x"
3231+ ],
3232+ "dev": true,
3233+ "license": "MIT",
3234+ "optional": true,
3235+ "os": [
3236+ "linux"
3237+ ],
3238+ "engines": {
3239+ "node": ">=18"
3240+ }
3241+ },
3242+ "node_modules/vite-node/node_modules/@esbuild/linux-x64": {
3243+ "version": "0.28.1",
3244+ "resolved": "https://registry.npmjs.org/@esbuild/linux-x64/-/linux-x64-0.28.1.tgz",
3245+ "integrity": "sha512-u/anNYF2mmVOEDwLtnQ1wOr3EZ9sTNGLWrsYGYwHWzGA3Si84IOkHXlbWTD1NB+9/1lcnweYKO54uhxZydNzfA==",
3246+ "cpu": [
3247+ "x64"
3248+ ],
3249+ "dev": true,
3250+ "license": "MIT",
3251+ "optional": true,
3252+ "os": [
3253+ "linux"
3254+ ],
3255+ "engines": {
3256+ "node": ">=18"
3257+ }
3258+ },
3259+ "node_modules/vite-node/node_modules/@esbuild/netbsd-x64": {
3260+ "version": "0.28.1",
3261+ "resolved": "https://registry.npmjs.org/@esbuild/netbsd-x64/-/netbsd-x64-0.28.1.tgz",
3262+ "integrity": "sha512-aeL6lAnN89Hz43Mlh1G8ARasbuoYvSITDEx0tHh5b7jJnHcssqgjy9Yx430GDpmCa6OyrKoS0aNRjKundRizGg==",
3263+ "cpu": [
3264+ "x64"
3265+ ],
3266+ "dev": true,
3267+ "license": "MIT",
3268+ "optional": true,
3269+ "os": [
3270+ "netbsd"
3271+ ],
3272+ "engines": {
3273+ "node": ">=18"
3274+ }
3275+ },
3276+ "node_modules/vite-node/node_modules/@esbuild/openbsd-x64": {
3277+ "version": "0.28.1",
3278+ "resolved": "https://registry.npmjs.org/@esbuild/openbsd-x64/-/openbsd-x64-0.28.1.tgz",
3279+ "integrity": "sha512-i/ZLIOafE0Z8cI/XANJAixoJL/uRAoS2xOA3rb0xN+KK0K177cMAsQYkzHtBrtMXAKuAc7HGgcWiZ/sRC1Nxgw==",
3280+ "cpu": [
3281+ "x64"
3282+ ],
3283+ "dev": true,
3284+ "license": "MIT",
3285+ "optional": true,
3286+ "os": [
3287+ "openbsd"
3288+ ],
3289+ "engines": {
3290+ "node": ">=18"
3291+ }
3292+ },
3293+ "node_modules/vite-node/node_modules/@esbuild/sunos-x64": {
3294+ "version": "0.28.1",
3295+ "resolved": "https://registry.npmjs.org/@esbuild/sunos-x64/-/sunos-x64-0.28.1.tgz",
3296+ "integrity": "sha512-BEjgtECkL3vY+SaSQ6nzVfiALUeFxpawyp8Jmf5PtYhf1Ug40N1h/hxlhts+f1FvSvarEigdxS3BlSMI2PJLcQ==",
3297+ "cpu": [
3298+ "x64"
3299+ ],
3300+ "dev": true,
3301+ "license": "MIT",
3302+ "optional": true,
3303+ "os": [
3304+ "sunos"
3305+ ],
3306+ "engines": {
3307+ "node": ">=18"
3308+ }
3309+ },
3310+ "node_modules/vite-node/node_modules/@esbuild/win32-arm64": {
3311+ "version": "0.28.1",
3312+ "resolved": "https://registry.npmjs.org/@esbuild/win32-arm64/-/win32-arm64-0.28.1.tgz",
3313+ "integrity": "sha512-lCv9eK/H6ZJWbE7bh2nw54CZ9M2nupBxJcTsdk/QQnWkdSjKGuxmmH8/GWrlT1eMmZfn4dGcCjRte397WqfQXA==",
3314+ "cpu": [
3315+ "arm64"
3316+ ],
3317+ "dev": true,
3318+ "license": "MIT",
3319+ "optional": true,
3320+ "os": [
3321+ "win32"
3322+ ],
3323+ "engines": {
3324+ "node": ">=18"
3325+ }
3326+ },
3327+ "node_modules/vite-node/node_modules/@esbuild/win32-ia32": {
3328+ "version": "0.28.1",
3329+ "resolved": "https://registry.npmjs.org/@esbuild/win32-ia32/-/win32-ia32-0.28.1.tgz",
3330+ "integrity": "sha512-zvb/mB2bSCoJOpoCBgYKKpX6YM6mJBlBUVUtVj41DlZJVEB6/0CKlRYxP5wWl1C1ILiCoAU5wZZ4q1P3qeS6Eg==",
3331+ "cpu": [
3332+ "ia32"
3333+ ],
3334+ "dev": true,
3335+ "license": "MIT",
3336+ "optional": true,
3337+ "os": [
3338+ "win32"
3339+ ],
3340+ "engines": {
3341+ "node": ">=18"
3342+ }
3343+ },
3344+ "node_modules/vite-node/node_modules/@esbuild/win32-x64": {
3345+ "version": "0.28.1",
3346+ "resolved": "https://registry.npmjs.org/@esbuild/win32-x64/-/win32-x64-0.28.1.tgz",
3347+ "integrity": "sha512-bm4Mowrv+GXMlpWX++EcXw/iLyd1o3+bJkC2DkWXYVvgZCqD/bSj9ctZeAMC3cIxgjRVR2Dufaiu4YPxr5gW1A==",
3348+ "cpu": [
3349+ "x64"
3350+ ],
3351+ "dev": true,
3352+ "license": "MIT",
3353+ "optional": true,
3354+ "os": [
3355+ "win32"
3356+ ],
3357+ "engines": {
3358+ "node": ">=18"
3359+ }
3360+ },
30113361 "node_modules/vite-node/node_modules/vite": {
30123362 "version": "8.1.5",
30133363 "resolved": "https://registry.npmjs.org/vite/-/vite-8.1.5.tgz",
src/geom/geometry.tsmodified+38−2View file
@@ -30,6 +30,8 @@
3030 */
3131 import { ShtPlan } from '../sht/sht.ts';
3232 import type { ShtConfig } from '../sht/layout.ts';
33+import type { DerivPlan } from '../sht/deriv.ts';
34+import { computeMetric } from './metric.ts';
3335 import { HostBuffers, ModelPlan } from '../mgpu/plan.ts';
3436 import { CompiledModel, type Binding } from '../mgpu/compile.ts';
3537 import { inFunction, inFunctionAsync, inModel } from '../mgpu/errors.ts';
@@ -48,6 +50,8 @@ export interface GeometryOptions {
4850 /** Parameter names the .m may take beyond `theta` and `phi`. */
4951 paramNames: string[];
5052 params: ModelParams;
53+ /** Computes the theta/phi derivatives the inverse metric quantities need. */
54+ deriv: DerivPlan;
5155 }
5256
5357 export class Geometry {
@@ -59,10 +63,23 @@ export class Geometry {
5963 readonly X: Float32Array;
6064 readonly Y: Float32Array;
6165 readonly Z: Float32Array;
66+ /**
67+ * Inverse metric quantities (src/geom/metric.ts), grid space, npts each.
68+ * 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.
70+ */
71+ readonly Vtx: Float32Array;
72+ readonly Vty: Float32Array;
73+ readonly Vtz: Float32Array;
74+ readonly Vpx: Float32Array;
75+ readonly Vpy: Float32Array;
76+ readonly Vpz: Float32Array;
6277
6378 private constructor(init: {
6479 x: Float32Array; y: Float32Array; z: Float32Array;
6580 X: Float32Array; Y: Float32Array; Z: Float32Array;
81+ Vtx: Float32Array; Vty: Float32Array; Vtz: Float32Array;
82+ Vpx: Float32Array; Vpy: Float32Array; Vpz: Float32Array;
6683 }) {
6784 this.x = init.x;
6885 this.y = init.y;
@@ -70,6 +87,12 @@ export class Geometry {
7087 this.X = init.X;
7188 this.Y = init.Y;
7289 this.Z = init.Z;
90+ this.Vtx = init.Vtx;
91+ this.Vty = init.Vty;
92+ this.Vtz = init.Vtz;
93+ this.Vpx = init.Vpx;
94+ this.Vpy = init.Vpy;
95+ this.Vpz = init.Vpz;
7396 }
7497
7598 /**
@@ -78,7 +101,7 @@ export class Geometry {
78101 * takes part in the timestep — so it reads back through the CPU freely.
79102 */
80103 static async create(opts: GeometryOptions): Promise<Geometry> {
81- const { device, sht, cfg, source, paramNames, params } = opts;
104+ const { device, sht, cfg, source, paramNames, params, deriv } = opts;
82105 const npts = cfg.nlat * cfg.nphi;
83106 const nlm = sht.nlm;
84107
@@ -128,7 +151,20 @@ export class Geometry {
128151 await sht.synth(Y),
129152 await sht.synth(Z),
130153 ];
131- return new Geometry({ x, y, z, X, Y, Z });
154+
155+ // Inverse metric quantities (algos.tex Algorithm 2): theta/phi
156+ // derivatives of the embedding's coefficients, contracted through the
157+ // inverse first fundamental form. Depends only on the geometry, so
158+ // this is a one-off alongside x,y,z above, not per-step work.
159+ const Xt = await deriv.dtheta(X);
160+ const Xp = await deriv.dphi(X);
161+ const Yt = await deriv.dtheta(Y);
162+ const Yp = await deriv.dphi(Y);
163+ const Zt = await deriv.dtheta(Z);
164+ const Zp = await deriv.dphi(Z);
165+ const { Vtx, Vty, Vtz, Vpx, Vpy, Vpz } = computeMetric(npts, Xt, Xp, Yt, Yp, Zt, Zp);
166+
167+ return new Geometry({ x, y, z, X, Y, Z, Vtx, Vty, Vtz, Vpx, Vpy, Vpz });
132168 } finally {
133169 plan.destroy();
134170 host.destroy();
src/geom/metric.tsadded+69−0View file
@@ -0,0 +1,69 @@
1+/**
2+ * Inverse metric quantities V_theta, V_phi of a surface embedding X=(x,y,z)
3+ * (evolving_surface/notes/algos.tex Algorithm 2 / SurfaceDiffOperator.
4+ * _precompute_metric_quantities, clear_denominators=False branch): six grid
5+ * scalar fields depending only on the geometry, used by the surface
6+ * Laplace-Beltrami operator (Algorithm 3) to contract a field's theta/phi
7+ * derivatives into a tangential gradient/divergence.
8+ *
9+ * g_tt = Xt.Xt, g_tp = Xt.Xp, g_pp = Xp.Xp (first fundamental form)
10+ * det = g_tt*g_pp - g_tp^2
11+ * V_theta = ( g_pp*Xt - g_tp*Xp ) / det
12+ * V_phi = ( g_tt*Xp - g_tp*Xt ) / det
13+ */
14+
15+export interface MetricFields {
16+ /** V_theta, Cartesian components, npts each. */
17+ Vtx: Float32Array;
18+ Vty: Float32Array;
19+ Vtz: Float32Array;
20+ /** V_phi, Cartesian components, npts each. */
21+ Vpx: Float32Array;
22+ Vpy: Float32Array;
23+ Vpz: Float32Array;
24+}
25+
26+/**
27+ * Xt/Xp (etc) are the theta/phi derivatives of each Cartesian embedding
28+ * component, grid space, npts each -- the tangent vectors X_theta, X_phi of
29+ * algos.tex Sec 4.1, one component per array.
30+ */
31+export function computeMetric(
32+ npts: number,
33+ Xt: Float32Array,
34+ Xp: Float32Array,
35+ Yt: Float32Array,
36+ Yp: Float32Array,
37+ Zt: Float32Array,
38+ Zp: Float32Array,
39+): MetricFields {
40+ const Vtx = new Float32Array(npts);
41+ const Vty = new Float32Array(npts);
42+ const Vtz = new Float32Array(npts);
43+ const Vpx = new Float32Array(npts);
44+ const Vpy = new Float32Array(npts);
45+ const Vpz = new Float32Array(npts);
46+
47+ for (let i = 0; i < npts; i++) {
48+ const xt = Xt[i];
49+ const xp = Xp[i];
50+ const yt = Yt[i];
51+ const yp = Yp[i];
52+ const zt = Zt[i];
53+ const zp = Zp[i];
54+
55+ const gtt = xt * xt + yt * yt + zt * zt;
56+ const gtp = xt * xp + yt * yp + zt * zp;
57+ const gpp = xp * xp + yp * yp + zp * zp;
58+ const det = gtt * gpp - gtp * gtp;
59+
60+ Vtx[i] = (gpp * xt - gtp * xp) / det;
61+ Vty[i] = (gpp * yt - gtp * yp) / det;
62+ Vtz[i] = (gpp * zt - gtp * zp) / det;
63+ Vpx[i] = (gtt * xp - gtp * xt) / det;
64+ Vpy[i] = (gtt * yp - gtp * yt) / det;
65+ Vpz[i] = (gtt * zp - gtp * zt) / det;
66+ }
67+
68+ return { Vtx, Vty, Vtz, Vpx, Vpy, Vpz };
69+}
src/main.tsmodified+56−8View file
@@ -111,9 +111,25 @@ const editor = new CodeEditor({
111111 },
112112 });
113113
114-/** Timesteps submitted per rendered frame. Nothing is read back between them,
115- * so the batch costs one submit and one readback regardless of size. */
116-const STEPS_PER_FRAME = 4;
114+/**
115+ * Timesteps submitted per rendered frame, at most. Nothing is read back
116+ * between them, so the batch costs one submit and one readback regardless of
117+ * size — but a compute pass is still real GPU work, and a browser's GPU
118+ * process enforces a watchdog timeout a headless desktop run does not: a
119+ * submission with enough dispatches in it can trip "device lost" outright,
120+ * on weak-enough hardware, well before it would ever show up as merely slow.
121+ * The `for k = 1:niter` correction loop makes a step's dispatch count scale
122+ * with niter (each iteration is ~15 dispatches per species — see
123+ * models/schnakenberg.m), so a fixed per-frame step count that was safe when
124+ * every model's step was a handful of dispatches is not safe once niter is
125+ * large. `stepsPerFrame`/`measureBurst` below scale it down — never up, so
126+ * the common case does not change — to keep one submission's total dispatch
127+ * count under DISPATCH_BUDGET regardless of how expensive the compiled step
128+ * is.
129+ */
130+const STEPS_PER_FRAME_BASE = 4;
131+/** See STEPS_PER_FRAME_BASE. Recomputed per rebuild in `rebuild()`. */
132+let stepsPerFrame = STEPS_PER_FRAME_BASE;
117133
118134 /**
119135 * Steps in a solver-timing burst, and how often to run one.
@@ -128,9 +144,25 @@ const STEPS_PER_FRAME = 4;
128144 * simulation — otherwise the pattern would visibly lurch forward at every
129145 * measurement.
130146 */
131-const MEASURE_BURST = 32;
147+const MEASURE_BURST_BASE = 32;
148+/** See STEPS_PER_FRAME_BASE — the measurement burst is one submission too,
149+ * and a bigger one: 32 steps is the single largest batch this app ever
150+ * submits, so it is the first thing to cross DISPATCH_BUDGET as niter grows. */
151+let measureBurst = MEASURE_BURST_BASE;
132152 const MEASURE_EVERY_MS = 2000;
133153
154+/**
155+ * Upper bound on dispatches in one submission — the frame batch and the
156+ * measurement burst are both scaled down to stay under this, never up, so
157+ * a cheap model's pacing is unchanged. Chosen well under what this project's
158+ * own desktop benchmark measures as trivially fast (single-digit ms even at
159+ * niter=8's ~450 dispatches/step), because the risk here is not GPU time on
160+ * capable hardware — it is a browser's GPU-process watchdog on weak
161+ * (integrated-graphics) hardware, which a headless desktop run never
162+ * exercises and this project has no way to benchmark directly.
163+ */
164+const DISPATCH_BUDGET = 1000;
165+
134166 /**
135167 * 'auto' display oversampling targets this many render latitudes: the factor is
136168 * the smallest power of two (up to 4) that reaches it. A solver grid already
@@ -581,7 +613,12 @@ async function rebuild(): Promise<void> {
581613 session = null;
582614 solverMs = 0;
583615 frameMs = 0;
584- lastMeasure = 0;
616+ // Not 0: with a large niter's dispatch count not yet known (that needs the
617+ // compiled plan below), the first measurement burst should wait for the
618+ // ordinary per-frame batch — already sized to this model — to prove itself
619+ // first, rather than firing a possibly-oversized burst before a single
620+ // frame has run.
621+ lastMeasure = performance.now();
585622 elErr.textContent = '';
586623 updateCommand();
587624 if (!device) return;
@@ -613,6 +650,13 @@ async function rebuild(): Promise<void> {
613650 plan.step.map((l) => ` ${l}`).join('\n');
614651 elRecompile.textContent = 'Recompile';
615652
653+ // Scale the frame batch and the measurement burst down — never up — so
654+ // neither submission's total dispatch count exceeds DISPATCH_BUDGET, no
655+ // matter how expensive niter has made one step. See STEPS_PER_FRAME_BASE.
656+ const opsPerStep = Math.max(1, plan.step.length);
657+ stepsPerFrame = Math.max(1, Math.min(STEPS_PER_FRAME_BASE, Math.floor(DISPATCH_BUDGET / opsPerStep)));
658+ measureBurst = Math.max(1, Math.min(MEASURE_BURST_BASE, Math.floor(DISPATCH_BUDGET / opsPerStep)));
659+
616660 const surface = await session.renderPositions();
617661 if (gen !== generation) return;
618662
@@ -718,7 +762,7 @@ async function pump(): Promise<void> {
718762 // benchmark's throughput number. State-preserving: the display and
719763 // model time are unaffected.
720764 if (performance.now() - lastMeasure > MEASURE_EVERY_MS) {
721- const ms = await session.measure(MEASURE_BURST);
765+ const ms = await session.measure(measureBurst);
722766 if (gen !== generation) break;
723767 solverMs = ms;
724768 lastMeasure = performance.now();
@@ -727,7 +771,7 @@ async function pump(): Promise<void> {
727771 // The frame itself. No explicit sync here — draw()'s readback already
728772 // waits for the steps, so asking twice would only add a round trip.
729773 const t0 = performance.now();
730- session.step(STEPS_PER_FRAME);
774+ session.step(stepsPerFrame);
731775 await draw();
732776 if (gen !== generation) break;
733777 frameMs = frameMs === 0
@@ -764,7 +808,11 @@ async function pump(): Promise<void> {
764808 async function benchmark(): Promise<void> {
765809 if (!session || movieBusy) return;
766810 setRunning(false);
767- const BATCH = 32;
811+ // Same base size and the same DISPATCH_BUDGET scaling as the automatic
812+ // measurement burst (see STEPS_PER_FRAME_BASE) — this is a user-triggered
813+ // 32-step submission, exactly the shape of thing that risks a browser's
814+ // GPU-process watchdog on weak hardware once niter makes a step expensive.
815+ const BATCH = measureBurst;
768816 const DURATION_MS = 2000;
769817 elBenchResult.textContent = 'benchmarking…';
770818 // A movie started mid-benchmark would replay while this loop still steps.
src/mgpu/externals.tsmodified+16−2View file
@@ -77,7 +77,13 @@ exports.cBody = function () {
7777 `;
7878 }
7979
80-/** Workspace files that make `synth` / `analys` resolvable during lowering. */
80+/**
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.
86+ */
8187 export function externalOpFiles(g: GridSizes): { name: string; source: string }[] {
8288 return [
8389 {
@@ -88,8 +94,16 @@ export function externalOpFiles(g: GridSizes): { name: string; source: string }[
8894 name: 'analys.mtoc2.js',
8995 source: transformSource('analys', g.npts, 1, 2, g.nlm),
9096 },
97+ {
98+ name: 'dtheta.mtoc2.js',
99+ source: transformSource('dtheta', 2, g.nlm, g.npts, 1),
100+ },
101+ {
102+ name: 'dphi.mtoc2.js',
103+ source: transformSource('dphi', 2, g.nlm, g.npts, 1),
104+ },
91105 ];
92106 }
93107
94108 /** Names the WGSL backend must implement as GPU encodes rather than kernels. */
95-export const EXTERNAL_OPS = new Set(['synth', 'analys']);
109+export const EXTERNAL_OPS = new Set(['synth', 'analys', 'dtheta', 'dphi']);
src/mgpu/model.tsmodified+53−5View file
@@ -18,6 +18,7 @@
1818 * name, so the file documents its own interface.
1919 */
2020 import { ShtPlan } from '../sht/sht.ts';
21+import type { DerivPlan } from '../sht/deriv.ts';
2122 import { lmIndex, type ShtConfig } from '../sht/layout.ts';
2223 import { HostBuffers, ModelPlan } from './plan.ts';
2324 import { inFunction, inFunctionAsync, inModel } from './errors.ts';
@@ -46,6 +47,9 @@ export interface GpuModelOptions {
4647 * sphere, where the .m has no geometry to take.
4748 */
4849 geometry?: GeometryBuffers;
50+ /** Computes `dtheta`/`dphi` for the .m's surface Laplace-Beltrami
51+ * correction. Omitted for a bare unit sphere, same as `geometry`. */
52+ deriv?: DerivPlan;
4953 /**
5054 * Iterations of the implicit solve the .m's `for` loop runs. A fixed scalar
5155 * rather than a tunable one: the loop is unrolled into the op sequence, so
@@ -64,11 +68,20 @@ export interface GeometryBuffers {
6468 X: Float32Array;
6569 Y: Float32Array;
6670 Z: Float32Array;
71+ /** Inverse metric quantities (src/geom/metric.ts), grid space, npts each. */
72+ Vtx: Float32Array;
73+ Vty: Float32Array;
74+ Vtz: Float32Array;
75+ Vpx: Float32Array;
76+ Vpy: Float32Array;
77+ Vpz: Float32Array;
6778 }
6879
6980 /** Names the .m may take for the grid coordinates and for their coefficients. */
7081 export const GEOMETRY_GRID_NAMES = ['gx', 'gy', 'gz'] as const;
7182 export const GEOMETRY_SPECTRAL_NAMES = ['Gx', 'Gy', 'Gz'] as const;
83+/** Names the .m may take for the inverse metric quantities. */
84+export const METRIC_GRID_NAMES = ['Vtx', 'Vty', 'Vtz', 'Vpx', 'Vpy', 'Vpz'] as const;
7285
7386 /** Laplace-Beltrami eigenvalues l(l+1), duplicated across re/im so the array
7487 * matches the 2 x nlm spectral layout element for element. */
@@ -84,6 +97,27 @@ export function eigenvalues(cfg: ShtConfig, nlm: number): Float32Array {
8497 return lam;
8598 }
8699
100+/**
101+ * 1 where l < lmax-2, else 0, duplicated across re/im like `lam`. The
102+ * theta/phi derivative recurrences (src/sht/derivCoeffs.ts) cannot exactly
103+ * represent a derivative at the top two degrees of the band limit, so the
104+ * surface Laplace-Beltrami correction filters them out wherever it
105+ * re-differentiates a field (evolving_surface/notes/algos.tex Sec 6,
106+ * "Miscellaneous implementation details").
107+ */
108+export function filterMask(cfg: ShtConfig, nlm: number): Float32Array {
109+ const filt = new Float32Array(2 * nlm);
110+ for (let m = 0; m <= cfg.mmax; m++) {
111+ for (let l = m; l <= cfg.lmax; l++) {
112+ const i = lmIndex(cfg.lmax, l, m);
113+ const keep = l < cfg.lmax - 2 ? 1 : 0;
114+ filt[2 * i] = keep;
115+ filt[2 * i + 1] = keep;
116+ }
117+ }
118+ return filt;
119+}
120+
87121 export class GpuModel {
88122 readonly paramNames: string[];
89123 readonly state: string[];
@@ -130,15 +164,17 @@ export class GpuModel {
130164 }
131165
132166 static async create(opts: GpuModelOptions): Promise<GpuModel> {
133- const { device, sht, cfg, source, paramNames, state, view, geometry } = opts;
167+ const { device, sht, cfg, source, paramNames, state, view, geometry, deriv } = opts;
134168 const npts = cfg.nlat * cfg.nphi;
135169 const nlm = sht.nlm;
136170 const niter = opts.niter ?? 0;
137171
138- // What the .m may ask for by parameter name. Spectral state and the
139- // eigenvalues are 2 x nlm; the seeded perturbation is a grid field.
172+ // What the .m may ask for by parameter name. Spectral state, the
173+ // eigenvalues and the top-mode filter are 2 x nlm; the seeded
174+ // perturbation is a grid field.
140175 const bindings: Record<string, Binding> = {
141176 lam: { kind: 'tensor', shape: [2, nlm] },
177+ filt: { kind: 'tensor', shape: [2, nlm] },
142178 noise: { kind: 'tensor', shape: [npts, 1] },
143179 npts: { kind: 'const', value: npts },
144180 nlm: { kind: 'const', value: nlm },
@@ -147,6 +183,7 @@ export class GpuModel {
147183 if (geometry) {
148184 for (const g of GEOMETRY_GRID_NAMES) bindings[g] = { kind: 'tensor', shape: [npts, 1] };
149185 for (const g of GEOMETRY_SPECTRAL_NAMES) bindings[g] = { kind: 'tensor', shape: [2, nlm] };
186+ for (const g of METRIC_GRID_NAMES) bindings[g] = { kind: 'tensor', shape: [npts, 1] };
150187 }
151188 for (const s of state) bindings[s] = { kind: 'tensor', shape: [2, nlm] };
152189 for (const p of paramNames) bindings[p] = { kind: 'param' };
@@ -169,20 +206,23 @@ export class GpuModel {
169206 // it writes it, and `step` reads it back.
170207 for (const s of state) host.ensure(s, 2 * nlm);
171208 host.ensure('lam', 2 * nlm);
209+ host.ensure('filt', 2 * nlm);
172210 host.ensure('noise', npts);
173211 if (geometry) {
174212 for (const g of GEOMETRY_GRID_NAMES) host.ensure(g, npts);
175213 for (const g of GEOMETRY_SPECTRAL_NAMES) host.ensure(g, 2 * nlm);
214+ for (const g of METRIC_GRID_NAMES) host.ensure(g, npts);
176215 }
177216
178217 const initPlan = await inFunctionAsync('init', () =>
179- ModelPlan.create(device, sht, { fn: initFn, feedback }, host),
218+ ModelPlan.create(device, sht, { fn: initFn, feedback }, host, deriv),
180219 );
181220 const stepPlan = await inFunctionAsync('step', () =>
182- ModelPlan.create(device, sht, { fn: stepFn, feedback }, host),
221+ ModelPlan.create(device, sht, { fn: stepFn, feedback }, host, deriv),
183222 );
184223
185224 host.upload('lam', eigenvalues(cfg, nlm));
225+ host.upload('filt', filterMask(cfg, nlm));
186226 if (geometry) {
187227 host.upload('gx', geometry.x);
188228 host.upload('gy', geometry.y);
@@ -190,6 +230,12 @@ export class GpuModel {
190230 host.upload('Gx', geometry.X);
191231 host.upload('Gy', geometry.Y);
192232 host.upload('Gz', geometry.Z);
233+ host.upload('Vtx', geometry.Vtx);
234+ host.upload('Vty', geometry.Vty);
235+ host.upload('Vtz', geometry.Vtz);
236+ host.upload('Vpx', geometry.Vpx);
237+ host.upload('Vpy', geometry.Vpy);
238+ host.upload('Vpz', geometry.Vpz);
193239 }
194240
195241 const readback = device.createBuffer({
@@ -233,6 +279,8 @@ export class GpuModel {
233279 const fields: [string, Float32Array][] = [
234280 ['gx', geometry.x], ['gy', geometry.y], ['gz', geometry.z],
235281 ['Gx', geometry.X], ['Gy', geometry.Y], ['Gz', geometry.Z],
282+ ['Vtx', geometry.Vtx], ['Vty', geometry.Vty], ['Vtz', geometry.Vtz],
283+ ['Vpx', geometry.Vpx], ['Vpy', geometry.Vpy], ['Vpz', geometry.Vpz],
236284 ];
237285 for (const [name, data] of fields) {
238286 if (this.#host.get(name)) this.#host.upload(name, data);
src/mgpu/plan.tsmodified+51−14View file
@@ -12,6 +12,7 @@ import { isMultiElement, scalarDouble } from 'numbl-src/numbl-core/jit/lowering/
1212 import type { Assign, For, IRExpr, IRStmt } from 'numbl-src/numbl-core/jit/lowering/ir.ts';
1313 import type { NumericType, Type } from 'numbl-src/numbl-core/jit/lowering/types.ts';
1414 import { ShtPlan, type ShtBinding } from '../sht/sht.ts';
15+import { DerivPlan, type DerivBinding } from '../sht/deriv.ts';
1516 import type { CompiledFunction } from './compile.ts';
1617 import { EXTERNAL_OPS } from './externals.ts';
1718 import {
@@ -119,6 +120,7 @@ type Op =
119120 copyBack?: { from: GPUBuffer; to: GPUBuffer; bytes: number };
120121 }
121122 | { kind: 'synth' | 'analys'; binding: ShtBinding; label: string }
123+ | { kind: 'dtheta' | 'dphi'; binding: DerivBinding; label: string }
122124 | { kind: 'copy'; from: GPUBuffer; to: GPUBuffer; bytes: number; label: string };
123125
124126 export interface PlanSpec {
@@ -192,6 +194,7 @@ export class ModelPlan {
192194
193195 #device: GPUDevice;
194196 #sht: ShtPlan;
197+ #deriv?: DerivPlan;
195198 #ops: Op[];
196199 #owned: GPUBuffer[];
197200 #paramBuf: GPUBuffer;
@@ -202,6 +205,7 @@ export class ModelPlan {
202205 private constructor(init: {
203206 device: GPUDevice;
204207 sht: ShtPlan;
208+ deriv?: DerivPlan;
205209 ops: Op[];
206210 byName: Map<string, Slot>;
207211 owned: GPUBuffer[];
@@ -211,6 +215,7 @@ export class ModelPlan {
211215 }) {
212216 this.#device = init.device;
213217 this.#sht = init.sht;
218+ this.#deriv = init.deriv;
214219 this.#ops = init.ops;
215220 this.#byName = init.byName;
216221 this.#owned = init.owned;
@@ -224,6 +229,8 @@ export class ModelPlan {
224229 sht: ShtPlan,
225230 spec: PlanSpec,
226231 host: HostBuffers,
232+ /** Computes dtheta/dphi — only needed if the .m calls them. */
233+ deriv?: DerivPlan,
227234 ): Promise<ModelPlan> {
228235 const { fn } = spec;
229236
@@ -297,7 +304,7 @@ export class ModelPlan {
297304 });
298305
299306 return new ModelPlan({
300- device, sht, ops, byName, owned, paramBuf, paramData, paramNames,
307+ device, sht, deriv, ops, byName, owned, paramBuf, paramData, paramNames,
301308 });
302309
303310 async function planStatement(stmt: IRStmt): Promise<void> {
@@ -348,19 +355,35 @@ export class ModelPlan {
348355 stmt.span,
349356 );
350357 }
351- ops.push(
352- ext.name === 'synth'
353- ? {
354- kind: 'synth',
355- binding: sht.createSynthBinding(argSlot.buffer, dest.buffer),
356- label: `${stmt.name} = synth(${ext.argName})`,
357- }
358- : {
359- kind: 'analys',
360- binding: sht.createAnalysBinding(argSlot.buffer, dest.buffer),
361- label: `${stmt.name} = analys(${ext.argName})`,
362- },
363- );
358+ const label = `${stmt.name} = ${ext.name}(${ext.argName})`;
359+ if (ext.name === 'synth') {
360+ ops.push({
361+ kind: 'synth',
362+ binding: sht.createSynthBinding(argSlot.buffer, dest.buffer),
363+ label,
364+ });
365+ } else if (ext.name === 'analys') {
366+ ops.push({
367+ kind: 'analys',
368+ binding: sht.createAnalysBinding(argSlot.buffer, dest.buffer),
369+ label,
370+ });
371+ } else if (ext.name === 'dtheta' || ext.name === 'dphi') {
372+ if (!deriv) {
373+ throw new UnsupportedOnGpu(
374+ `'${ext.name}' needs the surface's derivative transforms, ` +
375+ `which this plan was not given`,
376+ stmt.span,
377+ );
378+ }
379+ ops.push(
380+ ext.name === 'dtheta'
381+ ? { kind: 'dtheta', binding: deriv.createDthetaBinding(argSlot.buffer, dest.buffer), label }
382+ : { kind: 'dphi', binding: deriv.createDphiBinding(argSlot.buffer, dest.buffer), label },
383+ );
384+ } else {
385+ throw new UnsupportedOnGpu(`unknown external op '${ext.name}'`, stmt.span);
386+ }
364387 return;
365388 }
366389
@@ -382,6 +405,7 @@ export class ModelPlan {
382405 count,
383406 label,
384407 );
408+
385409 const bindGroupLayout = kernelLayout(device, tensors.size);
386410 const pipeline = await makePipeline(device, kernel.code, label, bindGroupLayout);
387411
@@ -542,6 +566,12 @@ export class ModelPlan {
542566 case 'analys':
543567 this.#shtInto(inPass(), op);
544568 break;
569+ case 'dtheta':
570+ this.#derivInto(inPass(), op);
571+ break;
572+ case 'dphi':
573+ this.#derivInto(inPass(), op);
574+ break;
545575 case 'copy':
546576 endPass();
547577 encoder.copyBufferToBuffer(op.from, 0, op.to, 0, op.bytes);
@@ -557,6 +587,13 @@ export class ModelPlan {
557587 else this.#sht.encodeAnalysInto(pass, op.binding);
558588 }
559589
590+ #derivInto(pass: GPUComputePassEncoder, op: Op & { kind: 'dtheta' | 'dphi' }): void {
591+ // planStatement already refused to plan a dtheta/dphi op without a
592+ // DerivPlan, so #deriv is guaranteed set whenever an op of this kind exists.
593+ if (op.kind === 'dtheta') this.#deriv!.encodeDthetaInto(pass, op.binding);
594+ else this.#deriv!.encodeDphiInto(pass, op.binding);
595+ }
596+
560597 /** Human-readable op sequence — what the .m actually compiled to. */
561598 describe(): string[] {
562599 return this.#ops.map((op) => `${op.kind.padEnd(7)} ${op.label}`);
src/mgpu/session.tsmodified+17−3View file
@@ -7,6 +7,7 @@
77 * browser-specific beyond needing a GPUDevice.
88 */
99 import { ShtPlan } from '../sht/sht.ts';
10+import { DerivPlan } from '../sht/deriv.ts';
1011 import { gridForLmax, type ShtConfig } from '../sht/layout.ts';
1112 import { GpuModel, type ModelParams } from './model.ts';
1213 import { seededNoise } from './noise.ts';
@@ -50,6 +51,8 @@ export class ModelSession {
5051 /** The surface being solved on, as spherical-harmonic coefficients. */
5152 #geometry: Geometry;
5253 #geometryModel: MGeometry;
54+ /** Computes the theta/phi derivatives a geometry's metric quantities need. */
55+ #deriv: DerivPlan;
5356
5457 /** Model time and step count since the last seeding. */
5558 t = 0;
@@ -71,6 +74,7 @@ export class ModelSession {
7174 oversample: number;
7275 geometry: Geometry;
7376 geometryModel: MGeometry;
77+ deriv: DerivPlan;
7478 niter: number;
7579 }) {
7680 this.device = init.device;
@@ -84,6 +88,7 @@ export class ModelSession {
8488 this.#displaySht = init.displaySht;
8589 this.#geometry = init.geometry;
8690 this.#geometryModel = init.geometryModel;
91+ this.#deriv = init.deriv;
8792 this.niter = init.niter;
8893 }
8994
@@ -110,6 +115,7 @@ export class ModelSession {
110115 const cfg = { lmax, mmax: lmax, nlat, nphi };
111116 const sht = await ShtPlan.create(device, cfg);
112117 let displaySht: ShtPlan | null = null;
118+ let deriv: DerivPlan | null = null;
113119 try {
114120 // The display plan shares nothing with the solver's beyond the
115121 // coefficients copied into it per readback; its grid is the solver's
@@ -123,10 +129,13 @@ export class ModelSession {
123129 nphi: oversample * nphi,
124130 });
125131 }
132+ // Computes the theta/phi derivatives the geometry's metric quantities
133+ // (and, per step, the surface Laplace-Beltrami correction) need.
134+ deriv = await DerivPlan.create(device, sht);
126135 // The surface is built before the model, because the model takes it as
127136 // an argument. It is a one-off: compiled, evaluated, read back, and its
128- // plan discarded — nothing of it survives into the timestep but six
129- // buffers of numbers.
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).
130139 const geometry = await Geometry.create({
131140 device,
132141 sht,
@@ -134,6 +143,7 @@ export class ModelSession {
134143 source: opts.geometrySource ?? geometryModel.source,
135144 paramNames: geometryModel.params.map((p) => p.key),
136145 params: geometryParams,
146+ deriv,
137147 });
138148 const gpu = await GpuModel.create({
139149 device,
@@ -144,15 +154,17 @@ export class ModelSession {
144154 state: model.state,
145155 view: model.species,
146156 geometry,
157+ deriv,
147158 niter,
148159 });
149160 gpu.setParams(params);
150161 return new ModelSession({
151162 device, model, cfg, sht, displaySht, gpu, params, oversample,
152- geometry, geometryModel, niter,
163+ geometry, geometryModel, deriv, niter,
153164 });
154165 } catch (e) {
155166 // The transform plans own GPU buffers; do not leak them on a compile error.
167+ deriv?.destroy();
156168 displaySht?.destroy();
157169 sht.destroy();
158170 throw e;
@@ -187,6 +199,7 @@ export class ModelSession {
187199 source: source ?? geometryModel.source,
188200 paramNames: geometryModel.params.map((p) => p.key),
189201 params,
202+ deriv: this.#deriv,
190203 });
191204 this.#geometry = next;
192205 this.#geometryModel = geometryModel;
@@ -296,6 +309,7 @@ export class ModelSession {
296309
297310 destroy(): void {
298311 this.gpu.destroy();
312+ this.#deriv.destroy();
299313 this.#displaySht?.destroy();
300314 this.sht.destroy();
301315 }
src/sht/deriv.tsadded+234−0View file
@@ -0,0 +1,234 @@
1+/**
2+ * First derivatives of a scalar field, coefficients -> grid: dtheta and dphi
3+ * (evolving_surface/notes/algos.tex Algorithm 1, theta/phi branches only --
4+ * the Laplace-Beltrami operator built on these never needs the second-
5+ * derivative/curvature branches, so they are not ported).
6+ *
7+ * Both derivatives start with a shuffle in coefficient space (the theta
8+ * branch's +-1 index gather via the alpha recurrence, the phi branch's i*m
9+ * row-swap) and then reuse the *existing* Legendre+Fourier synthesis
10+ * pipeline (ShtPlan.createSynthBinding/encodeSynthInto) unchanged -- neither
11+ * derivative touches the Legendre recurrence stage itself. dtheta
12+ * additionally divides by sin(theta) on the grid afterwards.
13+ */
14+import type { ShtPlan, ShtBinding } from './sht.ts';
15+import { derivCoeffs } from './derivCoeffs.ts';
16+import { dthetaShuffleWGSL, dphiShuffleWGSL, divideSinThetaWGSL } from './wgsl/deriv.ts';
17+
18+const WG = 64;
19+
20+async function makePipeline(
21+ device: GPUDevice,
22+ code: string,
23+ entryPoint: string,
24+): Promise<GPUComputePipeline> {
25+ device.pushErrorScope('validation');
26+ const module = device.createShaderModule({ code, label: entryPoint });
27+ const info = await module.getCompilationInfo();
28+ const errors = info.messages.filter((m) => m.type === 'error');
29+ if (errors.length) {
30+ throw new Error(
31+ `WGSL compile error in ${entryPoint}:\n` +
32+ errors.map((e) => ` ${e.lineNum}:${e.linePos} ${e.message}`).join('\n'),
33+ );
34+ }
35+ const pipeline = await device.createComputePipelineAsync({
36+ layout: 'auto',
37+ compute: { module, entryPoint },
38+ label: entryPoint,
39+ });
40+ const err = await device.popErrorScope();
41+ if (err) throw new Error(`pipeline ${entryPoint}: ${err.message}`);
42+ return pipeline;
43+}
44+
45+/** Bindings for one dtheta/dphi call against caller-supplied buffers. */
46+export interface DerivBinding {
47+ readonly shuffle: GPUBindGroup;
48+ readonly sht: ShtBinding;
49+ /** Only present for dtheta: the post-synthesis divide by sin(theta). */
50+ readonly divide?: GPUBindGroup;
51+}
52+
53+export class DerivPlan {
54+ private device: GPUDevice;
55+ private sht: ShtPlan;
56+ private nlm: number;
57+ private npts: number;
58+
59+ private bufAPlus!: GPUBuffer;
60+ private bufAMinus!: GPUBuffer;
61+ private bufMOf!: GPUBuffer;
62+ private bufSinTheta!: GPUBuffer;
63+ /** Scratch coefficient buffer for the shuffled input to synth -- shared
64+ * sequentially like ShtPlan's fmBuf, since ops within one pass execute
65+ * in submission order. */
66+ private scratch!: GPUBuffer;
67+
68+ private pipeDtheta!: GPUComputePipeline;
69+ private pipeDphi!: GPUComputePipeline;
70+ private pipeDivide!: GPUComputePipeline;
71+
72+ private constructor(device: GPUDevice, sht: ShtPlan) {
73+ this.device = device;
74+ this.sht = sht;
75+ this.nlm = sht.nlm;
76+ this.npts = sht.cfg.nlat * sht.cfg.nphi;
77+ }
78+
79+ static async create(device: GPUDevice, sht: ShtPlan): Promise<DerivPlan> {
80+ const plan = new DerivPlan(device, sht);
81+ await plan.init();
82+ return plan;
83+ }
84+
85+ private async init(): Promise<void> {
86+ const { nlat, nphi } = this.sht.cfg;
87+ const dev = this.device;
88+
89+ const { aPlus, aMinus, mOf } = derivCoeffs(this.sht.cfg.lmax, this.sht.cfg.mmax);
90+ const sinTheta = new Float32Array(nlat);
91+ for (let i = 0; i < nlat; i++) {
92+ const ct = this.sht.cosTheta[i];
93+ sinTheta[i] = Math.sqrt(Math.max(0, 1 - ct * ct));
94+ }
95+
96+ const mk = (label: string, size: number, usage = GPUBufferUsage.STORAGE | GPUBufferUsage.COPY_DST) =>
97+ dev.createBuffer({ label, size, usage });
98+ this.bufAPlus = mk('deriv-aplus', 4 * this.nlm);
99+ this.bufAMinus = mk('deriv-aminus', 4 * this.nlm);
100+ this.bufMOf = mk('deriv-mof', 4 * this.nlm);
101+ this.bufSinTheta = mk('deriv-sintheta', 4 * nlat);
102+ this.scratch = mk(
103+ 'deriv-scratch',
104+ 8 * this.nlm,
105+ GPUBufferUsage.STORAGE | GPUBufferUsage.COPY_DST | GPUBufferUsage.COPY_SRC,
106+ );
107+
108+ dev.queue.writeBuffer(this.bufAPlus, 0, new Float32Array(aPlus));
109+ dev.queue.writeBuffer(this.bufAMinus, 0, new Float32Array(aMinus));
110+ dev.queue.writeBuffer(this.bufMOf, 0, mOf as Uint32Array<ArrayBuffer>);
111+ dev.queue.writeBuffer(this.bufSinTheta, 0, sinTheta);
112+
113+ const [pDtheta, pDphi, pDivide] = await Promise.all([
114+ makePipeline(dev, dthetaShuffleWGSL({ nlm: this.nlm }), 'dtheta_shuffle'),
115+ makePipeline(dev, dphiShuffleWGSL({ nlm: this.nlm }), 'dphi_shuffle'),
116+ makePipeline(dev, divideSinThetaWGSL({ nlat, nphi }), 'divide_sin_theta'),
117+ ]);
118+ this.pipeDtheta = pDtheta;
119+ this.pipeDphi = pDphi;
120+ this.pipeDivide = pDivide;
121+ }
122+
123+ /** Bindings for dtheta(qlmIn) -> spatOut, against caller-owned buffers. */
124+ createDthetaBinding(qlmIn: GPUBuffer, spatOut: GPUBuffer): DerivBinding {
125+ const shuffle = this.device.createBindGroup({
126+ layout: this.pipeDtheta.getBindGroupLayout(0),
127+ entries: [
128+ { binding: 0, resource: { buffer: this.bufAPlus } },
129+ { binding: 1, resource: { buffer: this.bufAMinus } },
130+ { binding: 2, resource: { buffer: qlmIn } },
131+ { binding: 3, resource: { buffer: this.scratch } },
132+ ],
133+ });
134+ const sht = this.sht.createSynthBinding(this.scratch, spatOut);
135+ const divide = this.device.createBindGroup({
136+ layout: this.pipeDivide.getBindGroupLayout(0),
137+ entries: [
138+ { binding: 0, resource: { buffer: this.bufSinTheta } },
139+ { binding: 1, resource: { buffer: spatOut } },
140+ ],
141+ });
142+ return { shuffle, sht, divide };
143+ }
144+
145+ /** Bindings for dphi(qlmIn) -> spatOut, against caller-owned buffers. */
146+ 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+ });
155+ const sht = this.sht.createSynthBinding(this.scratch, spatOut);
156+ return { shuffle, sht };
157+ }
158+
159+ /** Record dtheta into an existing compute pass. */
160+ 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);
165+ pass.setPipeline(this.pipeDivide);
166+ pass.setBindGroup(0, b.divide!);
167+ pass.dispatchWorkgroups(Math.ceil(this.npts / WG));
168+ }
169+
170+ /** Record dphi into an existing compute pass. */
171+ 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));
175+ this.sht.encodeSynthInto(pass, b.sht);
176+ }
177+
178+ /** CPU convenience: qlm (interleaved [re,im], length 2*nlm) -> grid field. */
179+ async dtheta(qlm: Float32Array): Promise<Float32Array> {
180+ return this.#runToGrid(qlm, true);
181+ }
182+
183+ /** CPU convenience: qlm (interleaved [re,im], length 2*nlm) -> grid field. */
184+ async dphi(qlm: Float32Array): Promise<Float32Array> {
185+ return this.#runToGrid(qlm, false);
186+ }
187+
188+ async #runToGrid(qlm: Float32Array, withDivide: boolean): Promise<Float32Array> {
189+ if (qlm.length !== 2 * this.nlm) throw new Error(`qlm must have length ${2 * this.nlm}`);
190+ const dev = this.device;
191+ const qlmIn = dev.createBuffer({
192+ label: 'deriv-qlm-in',
193+ size: 8 * this.nlm,
194+ usage: GPUBufferUsage.STORAGE | GPUBufferUsage.COPY_DST,
195+ });
196+ const spatOut = dev.createBuffer({
197+ label: 'deriv-spat-out',
198+ size: 4 * this.npts,
199+ usage: GPUBufferUsage.STORAGE | GPUBufferUsage.COPY_SRC,
200+ });
201+ const stage = dev.createBuffer({
202+ label: 'deriv-stage',
203+ size: 4 * this.npts,
204+ usage: GPUBufferUsage.MAP_READ | GPUBufferUsage.COPY_DST,
205+ });
206+ try {
207+ dev.queue.writeBuffer(qlmIn, 0, qlm as Float32Array<ArrayBuffer>);
208+ const binding = withDivide
209+ ? this.createDthetaBinding(qlmIn, spatOut)
210+ : this.createDphiBinding(qlmIn, spatOut);
211+ const enc = dev.createCommandEncoder({ label: 'deriv-run' });
212+ const pass = enc.beginComputePass({ label: 'deriv-run' });
213+ if (withDivide) this.encodeDthetaInto(pass, binding);
214+ else this.encodeDphiInto(pass, binding);
215+ pass.end();
216+ enc.copyBufferToBuffer(spatOut, 0, stage, 0, 4 * this.npts);
217+ dev.queue.submit([enc.finish()]);
218+ await stage.mapAsync(GPUMapMode.READ);
219+ const out = new Float32Array(stage.getMappedRange().slice(0));
220+ stage.unmap();
221+ return out;
222+ } finally {
223+ qlmIn.destroy();
224+ spatOut.destroy();
225+ stage.destroy();
226+ }
227+ }
228+
229+ destroy(): void {
230+ for (const b of [
231+ this.bufAPlus, this.bufAMinus, this.bufMOf, this.bufSinTheta, this.scratch,
232+ ]) b?.destroy();
233+ }
234+}
src/sht/derivCoeffs.tsadded+49−0View file
@@ -0,0 +1,49 @@
1+/**
2+ * Recurrence coefficients for the first theta-derivative of orthonormal
3+ * associated Legendre functions (Condon-Shortley phase included), matching
4+ * the alpha^+/alpha^- recurrence in evolving_surface/notes/algos.tex Sec 2.1:
5+ *
6+ * sin(theta) d/dtheta Y_l^m = alpha^+(l,m) Y_{l+1}^m + alpha^-(l,m) Y_{l-1}^m
7+ *
8+ * so the coefficients of sin(theta)*dtheta(u), by degree, are
9+ *
10+ * v_l^m = alpha^+(l-1,m) u_{l-1}^m + alpha^-(l+1,m) u_{l+1}^m
11+ *
12+ * dropping any term referring to a degree outside 0 <= l <= lmax. Baked to
13+ * zero at each m-block's first/last element (rather than left undefined), so
14+ * a consuming WGSL kernel needs only an in-bounds check, not a validity check.
15+ */
16+import { lmIndex, nlmCalc } from './layout.ts';
17+
18+export interface DerivCoeffs {
19+ /** aPlus[lm] = alpha^+(l-1,m) when l>m, else 0 -- multiplies u_{l-1}^m. */
20+ aPlus: Float64Array;
21+ /** aMinus[lm] = alpha^-(l+1,m) when l<lmax, else 0 -- multiplies u_{l+1}^m. */
22+ aMinus: Float64Array;
23+ /** m of the coefficient at flat index lm (the phi-derivative needs only this). */
24+ mOf: Uint32Array;
25+}
26+
27+export function alphaPlus(l: number, m: number): number {
28+ return l * Math.sqrt(((l - m + 1) * (l + m + 1)) / ((2 * l + 1) * (2 * l + 3)));
29+}
30+
31+export function alphaMinus(l: number, m: number): number {
32+ return -(l + 1) * Math.sqrt(((l - m) * (l + m)) / ((2 * l - 1) * (2 * l + 1)));
33+}
34+
35+export function derivCoeffs(lmax: number, mmax: number): DerivCoeffs {
36+ const nlm = nlmCalc(lmax, mmax);
37+ const aPlus = new Float64Array(nlm);
38+ const aMinus = new Float64Array(nlm);
39+ const mOf = new Uint32Array(nlm);
40+ for (let m = 0; m <= mmax; m++) {
41+ for (let l = m; l <= lmax; l++) {
42+ const lm = lmIndex(lmax, l, m);
43+ mOf[lm] = m;
44+ if (l - 1 >= m) aPlus[lm] = alphaPlus(l - 1, m);
45+ if (l + 1 <= lmax) aMinus[lm] = alphaMinus(l + 1, m);
46+ }
47+ }
48+ return { aPlus, aMinus, mOf };
49+}
src/sht/reference.tsmodified+51−0View file
@@ -8,6 +8,7 @@
88 */
99 import { gaussNodesWeights } from './gauss.ts';
1010 import { legendreCoeffs, legendreRow, type LegendreCoeffs } from './coeffs.ts';
11+import { alphaPlus, alphaMinus } from './derivCoeffs.ts';
1112 import { lmIndex, nlmCalc, validateConfig, type ShtConfig } from './layout.ts';
1213
1314 export class ShtReference {
@@ -116,6 +117,56 @@ export class ShtReference {
116117 }
117118 return qlm;
118119 }
120+
121+ /**
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).
124+ */
125+ dtheta(qlm: ArrayLike<number>): Float64Array {
126+ const { lmax, mmax, nlat, nphi } = this.cfg;
127+ const v = new Float64Array(2 * this.nlm);
128+ for (let m = 0; m <= mmax; m++) {
129+ for (let l = m; l <= lmax; l++) {
130+ const lm = lmIndex(lmax, l, m);
131+ let re = 0;
132+ let im = 0;
133+ if (l - 1 >= m) {
134+ const lm1 = lmIndex(lmax, l - 1, m);
135+ const a = alphaPlus(l - 1, m);
136+ re += a * qlm[2 * lm1];
137+ im += a * qlm[2 * lm1 + 1];
138+ }
139+ if (l + 1 <= lmax) {
140+ const lm2 = lmIndex(lmax, l + 1, m);
141+ const a = alphaMinus(l + 1, m);
142+ re += a * qlm[2 * lm2];
143+ im += a * qlm[2 * lm2 + 1];
144+ }
145+ v[2 * lm] = re;
146+ v[2 * lm + 1] = im;
147+ }
148+ }
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;
155+ }
156+
157+ /** Phi-derivative, f64: (dphi u)_l^m = i*m*u_l^m, then synthesize. */
158+ dphi(qlm: ArrayLike<number>): Float64Array {
159+ const { lmax, mmax } = this.cfg;
160+ const v = new Float64Array(2 * this.nlm);
161+ for (let m = 0; m <= mmax; m++) {
162+ for (let l = m; l <= lmax; l++) {
163+ const lm = lmIndex(lmax, l, m);
164+ v[2 * lm] = -m * qlm[2 * lm + 1];
165+ v[2 * lm + 1] = m * qlm[2 * lm];
166+ }
167+ }
168+ return this.synth(v);
169+ }
119170 }
120171
121172 /** Random band-limited spectrum for testing (m=0 imaginary parts zeroed). */
src/sht/wgsl/deriv.tsadded+101−0View file
@@ -0,0 +1,101 @@
1+/**
2+ * WGSL kernels for the coefficient-space step of the theta/phi first-
3+ * derivative algorithm (evolving_surface/notes/algos.tex, Algorithm 1, theta
4+ * and phi branches only -- the Laplace-Beltrami operator built on these never
5+ * needs the second-derivative/curvature branches).
6+ *
7+ * Both derivatives are a shuffle across nearby spectral coefficients --
8+ * independent of latitude/longitude, so neither touches the Legendre
9+ * recurrence or Fourier stages in leg.ts/fourier.ts -- followed by the
10+ * *existing*, unchanged scalar synthesis pipeline. dtheta additionally
11+ * divides the synthesized grid field by sin(theta) afterwards.
12+ */
13+
14+const WG = 64;
15+
16+export interface DerivCoeffParams {
17+ nlm: number;
18+}
19+
20+/**
21+ * v_l^m = alpha^+(l-1,m) * u_{l-1}^m + alpha^-(l+1,m) * u_{l+1}^m, the
22+ * coefficients of sin(theta) * dtheta(u) (algos.tex eq. v_coeffs). aPlus/
23+ * aMinus are precomputed zero at each m-block's boundary
24+ * (src/sht/derivCoeffs.ts), so the multiply is always mathematically
25+ * correct; the bounds checks below exist only to avoid reading past the ends
26+ * of the qlm array (an m-block-internal +-1 step never leaves the array, so
27+ * this is the only place it could).
28+ */
29+export function dthetaShuffleWGSL(p: DerivCoeffParams): string {
30+ return /* wgsl */ `
31+const NLM: u32 = ${p.nlm}u;
32+
33+@group(0) @binding(0) var<storage, read> aPlus: array<f32>;
34+@group(0) @binding(1) var<storage, read> aMinus: array<f32>;
35+@group(0) @binding(2) var<storage, read> qlmIn: array<vec2f>;
36+@group(0) @binding(3) var<storage, read_write> vOut: array<vec2f>;
37+
38+@compute @workgroup_size(${WG})
39+fn dtheta_shuffle(@builtin(global_invocation_id) gid: vec3u) {
40+ let lm = gid.x;
41+ if (lm >= NLM) { return; }
42+ var v = vec2f(0.0);
43+ if (lm > 0u) { v += aPlus[lm] * qlmIn[lm - 1u]; }
44+ if (lm + 1u < NLM) { v += aMinus[lm] * qlmIn[lm + 1u]; }
45+ vOut[lm] = v;
46+}
47+`;
48+}
49+
50+/**
51+ * (dphi u)_l^m = i*m*u_l^m: in the [re, im] row layout this swaps and
52+ * negates, re' = -m*im, im' = m*re (algos.tex eq. dYdphi).
53+ */
54+export function dphiShuffleWGSL(p: DerivCoeffParams): string {
55+ return /* wgsl */ `
56+const NLM: u32 = ${p.nlm}u;
57+
58+@group(0) @binding(0) var<storage, read> mOf: array<u32>;
59+@group(0) @binding(1) var<storage, read> qlmIn: array<vec2f>;
60+@group(0) @binding(2) var<storage, read_write> vOut: array<vec2f>;
61+
62+@compute @workgroup_size(${WG})
63+fn dphi_shuffle(@builtin(global_invocation_id) gid: vec3u) {
64+ let lm = gid.x;
65+ if (lm >= NLM) { return; }
66+ let m = f32(mOf[lm]);
67+ let c = qlmIn[lm];
68+ vOut[lm] = vec2f(-m * c.y, m * c.x);
69+}
70+`;
71+}
72+
73+export interface DivideParams {
74+ nlat: number;
75+ nphi: number;
76+}
77+
78+/**
79+ * Elementwise divide by sin(theta): the grid-space finish of Algorithm 1's
80+ * dtheta branch (dtheta(u) = synth(v_l^m) / sin(theta)). Gauss nodes never
81+ * sit at the poles, so this never divides by zero.
82+ */
83+export function divideSinThetaWGSL(p: DivideParams): string {
84+ const npts = p.nlat * p.nphi;
85+ return /* wgsl */ `
86+const NLAT: u32 = ${p.nlat}u;
87+const NPHI: u32 = ${p.nphi}u;
88+const NPTS: u32 = ${npts}u;
89+
90+@group(0) @binding(0) var<storage, read> sinTheta: array<f32>;
91+@group(0) @binding(1) var<storage, read_write> spat: array<f32>;
92+
93+@compute @workgroup_size(${WG})
94+fn divide_sin_theta(@builtin(global_invocation_id) gid: vec3u) {
95+ let i = gid.x;
96+ if (i >= NPTS) { return; }
97+ let ilat = i / NPHI;
98+ spat[i] = spat[i] / sinTheta[ilat];
99+}
100+`;
101+}
test/analyticChecks.tsmodified+12−5View file
@@ -21,6 +21,7 @@
2121 * (~1e-7 relative) rather than by the scheme.
2222 */
2323 import { ShtPlan } from '../src/sht/sht.ts';
24+import { DerivPlan } from '../src/sht/deriv.ts';
2425 import { gridForLmax, lmIndex, nlmCalc, type ShtConfig } from '../src/sht/layout.ts';
2526 import { GpuModel } from '../src/mgpu/model.ts';
2627 import { mModelByKey, defaultParams, type MModel, type ParamSpec } from '../src/mgpu/registry.ts';
@@ -60,8 +61,9 @@ async function makeModel(
6061 model: MModel,
6162 cfg: ShtConfig,
6263 niter = 1,
63-): Promise<{ sht: ShtPlan; gpu: GpuModel }> {
64+): Promise<{ sht: ShtPlan; deriv: DerivPlan; gpu: GpuModel }> {
6465 const sht = await ShtPlan.create(device, cfg);
66+ const deriv = await DerivPlan.create(device, sht);
6567 const geometry = await Geometry.create({
6668 device,
6769 sht,
@@ -69,6 +71,7 @@ async function makeModel(
6971 source: mGeometryByKey(SPHERE_KEY)!.source,
7072 paramNames: [],
7173 params: {},
74+ deriv,
7275 });
7376 const gpu = await GpuModel.create({
7477 device,
@@ -79,9 +82,10 @@ async function makeModel(
7982 state: model.state,
8083 view: model.species,
8184 geometry,
85+ deriv,
8286 niter,
8387 });
84- return { sht, gpu };
88+ return { sht, deriv, gpu };
8589 }
8690
8791 export async function analyticChecks(
@@ -101,7 +105,7 @@ export async function analyticChecks(
101105 const nsteps = 20;
102106
103107 const model = testModel('linear', linearSource, ['c', 'D', 'dt']);
104- const { sht, gpu } = await makeModel(device, model, cfg);
108+ const { sht, deriv, gpu } = await makeModel(device, model, cfg);
105109 gpu.setParams({ c, D, dt });
106110
107111 // A single (l, m) mode, written straight into the spectral state.
@@ -137,6 +141,7 @@ export async function analyticChecks(
137141 check('A: no leakage into other modes', leak < 2e-6, `max |other| ${leak.toExponential(2)}`);
138142
139143 gpu.destroy();
144+ deriv.destroy();
140145 sht.destroy();
141146 }
142147
@@ -153,7 +158,7 @@ export async function analyticChecks(
153158 const u0 = 0.3;
154159
155160 const model = testModel('logistic', logisticSource, ['r', 'D', 'dt']);
156- const { sht, gpu } = await makeModel(device, model, cfg);
161+ const { sht, deriv, gpu } = await makeModel(device, model, cfg);
157162 gpu.setParams({ r, D, dt });
158163
159164 // Uniform initial field: stays uniform, and diffusion cannot touch it.
@@ -196,6 +201,7 @@ export async function analyticChecks(
196201 );
197202
198203 gpu.destroy();
204+ deriv.destroy();
199205 sht.destroy();
200206 }
201207
@@ -209,7 +215,7 @@ export async function analyticChecks(
209215 const nlm = nlmCalc(lmax, lmax);
210216 const npts = nlat * nphi;
211217
212- const { sht, gpu } = await makeModel(device, model, cfg);
218+ const { sht, deriv, gpu } = await makeModel(device, model, cfg);
213219 gpu.setParams(p);
214220
215221 // Seed the exact homogeneous fixed point by handing init a zero
@@ -271,6 +277,7 @@ export async function analyticChecks(
271277 log(` C: growth over ${nsteps} steps = ${(Math.abs(cu) / eps).toFixed(3)}x (predicted)`);
272278
273279 gpu.destroy();
280+ deriv.destroy();
274281 sht.destroy();
275282 }
276283 }
test/geometryChecks.tsmodified+170−21View file
@@ -8,14 +8,24 @@
88 * same surface on a finer grid.
99 *
1010 * The loop is checked for the property the whole design rests on: it is
11- * unrolled into the fixed op sequence, so more iterations means more GPU ops —
12- * and, while the geometry correction inside it is identically zero, the answer
13- * must be *bit for bit* independent of how many times it runs. That is a
14- * stronger statement than "close enough": if the placeholder were ever
15- * something that merely rounds to zero, or if the loop were miscompiled to
16- * read a stale buffer, these would differ in the last bits and this fails.
11+ * unrolled into the fixed op sequence, so more iterations means more GPU ops.
12+ * On the sphere, where the surface Laplace-Beltrami correction is
13+ * mathematically zero (lap_g = lap_s exactly), the answer must stay close
14+ * across niter to fp32 tolerance — not bit-identical, since the correction is
15+ * now a real (if numerically near-zero) computation rather than the literal
16+ * `0 * Un` placeholder, so the op sequence differs even though the answer
17+ * shouldn't move much. On a genuinely curved surface the correction must
18+ * actually change the answer, and — since the Richardson iteration only
19+ * converges while the correction stays small relative to what the
20+ * round-sphere solve inverts (docs/richardson-iteration.md) — a niter/dt/
21+ * geometry combination outside that radius is expected to diverge. The
22+ * niter x geometry sweep below documents which shipped combinations that
23+ * currently affects, so a regression that makes a *currently-healthy*
24+ * combination diverge is caught without this file silently asserting away a
25+ * real, known numerical limit.
1726 */
1827 import { ShtPlan } from '../src/sht/sht.ts';
28+import { DerivPlan } from '../src/sht/deriv.ts';
1929 import { gridForLmax, lmIndex } from '../src/sht/layout.ts';
2030 import { ModelSession } from '../src/mgpu/session.ts';
2131 import { mModelByKey, defaultParams } from '../src/mgpu/registry.ts';
@@ -31,6 +41,12 @@ import type { Check, Log } from './analyticChecks.ts';
3141
3242 const LMAX = 31;
3343 const STEPS = 20;
44+/** The app's actual default lmax (README: "at the default lmax 63 that is a
45+ * 128x256 grid"), used for the niter/geometry sweep below and the peanut
46+ * check next to it -- the divergence they're both about is a real, lmax-
47+ * dependent numerical property of the Richardson iteration, not one this
48+ * file's other, smaller LMAX happens to reproduce. */
49+const SWEEP_LMAX = 63;
3450
3551 /** Build one geometry on its own transform plan, for inspection. */
3652 async function buildGeometry(device: GPUDevice, key: string) {
@@ -38,6 +54,7 @@ async function buildGeometry(device: GPUDevice, key: string) {
3854 const { nlat, nphi } = gridForLmax(LMAX, 3);
3955 const cfg = { lmax: LMAX, mmax: LMAX, nlat, nphi };
4056 const sht = await ShtPlan.create(device, cfg);
57+ const deriv = await DerivPlan.create(device, sht);
4158 const geometry = await Geometry.create({
4259 device,
4360 sht,
@@ -45,8 +62,9 @@ async function buildGeometry(device: GPUDevice, key: string) {
4562 source: g.source,
4663 paramNames: g.params.map((p) => p.key),
4764 params: defaultGeometryParams(g),
65+ deriv,
4866 });
49- return { g, sht, cfg, geometry };
67+ return { g, sht, deriv, cfg, geometry };
5068 }
5169
5270 export async function geometryChecks(
@@ -56,23 +74,27 @@ export async function geometryChecks(
5674 ): Promise<void> {
5775 // ---- every geometry compiles and closes ---------------------------------
5876 for (const spec of mGeometries) {
59- const { sht, geometry } = await buildGeometry(device, spec.key);
77+ const { sht, deriv, geometry } = await buildGeometry(device, spec.key);
6078 let finite = true;
6179 for (const a of [geometry.x, geometry.y, geometry.z]) {
6280 for (const v of a) if (!Number.isFinite(v)) finite = false;
6381 }
82+ for (const a of [geometry.Vtx, geometry.Vty, geometry.Vtz, geometry.Vpx, geometry.Vpy, geometry.Vpz]) {
83+ for (const v of a) if (!Number.isFinite(v)) finite = false;
84+ }
6485 const { lo, hi } = geometry.radiusRange();
6586 check(
6687 `geometry: ${spec.key}.m evaluates to a finite surface`,
6788 finite && lo > 1e-3,
6889 `radius ${lo.toFixed(4)}–${hi.toFixed(4)}`,
6990 );
91+ deriv.destroy();
7092 sht.destroy();
7193 }
7294
7395 // ---- the sphere is the unit sphere, exactly, and is degree 1 ------------
7496 {
75- const { sht, geometry } = await buildGeometry(device, SPHERE_KEY);
97+ const { sht, deriv, geometry } = await buildGeometry(device, SPHERE_KEY);
7698
7799 let maxRadiusErr = 0;
78100 for (let i = 0; i < geometry.x.length; i++) {
@@ -111,12 +133,51 @@ export async function geometryChecks(
111133 leak < 1e-3,
112134 `max |coefficient| outside l = 1 is ${leak.toExponential(2)}`,
113135 );
136+
137+ // The inverse metric quantities have a closed form on the unit sphere:
138+ // V_theta = (cos(theta)cos(phi), cos(theta)sin(phi), -sin(theta)),
139+ // V_phi = (-sin(phi)/sin(theta), cos(phi)/sin(theta), 0). Checking these
140+ // pins the sign convention of computeMetric (src/geom/metric.ts) before
141+ // it is buried under the Laplace-Beltrami operator built on top of it.
142+ let maxMetricErr = 0;
143+ for (let i = 0; i < sht.cfg.nlat; i++) {
144+ const ct = sht.cosTheta[i];
145+ const st = Math.sqrt(Math.max(0, 1 - ct * ct));
146+ for (let j = 0; j < sht.cfg.nphi; j++) {
147+ const phi = (2 * Math.PI * j) / sht.cfg.nphi;
148+ const k = i * sht.cfg.nphi + j;
149+ const cphi = Math.cos(phi);
150+ const sphi = Math.sin(phi);
151+ const wantVtx = ct * cphi;
152+ const wantVty = ct * sphi;
153+ const wantVtz = -st;
154+ const wantVpx = -sphi / st;
155+ const wantVpy = cphi / st;
156+ const wantVpz = 0;
157+ maxMetricErr = Math.max(
158+ maxMetricErr,
159+ Math.abs(geometry.Vtx[k] - wantVtx),
160+ Math.abs(geometry.Vty[k] - wantVty),
161+ Math.abs(geometry.Vtz[k] - wantVtz),
162+ Math.abs(geometry.Vpx[k] - wantVpx),
163+ Math.abs(geometry.Vpy[k] - wantVpy),
164+ Math.abs(geometry.Vpz[k] - wantVpz),
165+ );
166+ }
167+ }
168+ check(
169+ 'geometry: sphere.m has the closed-form inverse metric quantities',
170+ maxMetricErr < 2e-3,
171+ `max |V - closed form| = ${maxMetricErr.toExponential(2)}`,
172+ );
173+
174+ deriv.destroy();
114175 sht.destroy();
115176 }
116177
117178 // ---- a deformed surface matches its own formula, on any grid ------------
118179 {
119- const { g, sht, cfg, geometry } = await buildGeometry(device, 'peanut');
180+ const { g, sht, deriv, cfg, geometry } = await buildGeometry(device, 'peanut');
120181 const p = defaultGeometryParams(g);
121182
122183 // peanut.m written out: r = 1 - waist*sin(theta)^2 scales the unit sphere,
@@ -177,6 +238,7 @@ export async function geometryChecks(
177238 `max |dr| = ${refined.toExponential(2)} at ${2 * cfg.nlat}×${2 * cfg.nphi} points`,
178239 );
179240 fine.destroy();
241+ deriv.destroy();
180242 sht.destroy();
181243 }
182244
@@ -208,10 +270,14 @@ export async function geometryChecks(
208270 `${ops.join(' < ')} ops for ${counts.join(', ')} iterations`,
209271 );
210272 // Unrolling has to be exactly linear in the trip count: the body planned
211- // once per iteration, no more and no less. Two dispatches per species per
212- // iteration — the placeholder line and the update that reads it.
273+ // once per iteration, no more and no less. Per species per iteration: 8
274+ // dtheta/dphi + 4 analys transforms (Algorithm 3's cost, applied to the
275+ // field and to each of its three Cartesian gradient components) plus 15
276+ // generated kernels -- see test/modelChecks.ts's KERNELS_PER_ITERATION,
277+ // which counts the kernels alone; this counts every op, transforms
278+ // included.
213279 const perIteration = ops[1] - ops[0];
214- const want = 2 * model.species.length;
280+ const want = 54;
215281 check(
216282 'loop: unrolling is exactly linear in the trip count',
217283 perIteration === want && ops[2] - ops[0] === 4 * perIteration,
@@ -219,24 +285,107 @@ export async function geometryChecks(
219285 `${ops[2] - ops[0]} for 4 iterations`,
220286 );
221287
222- let identical = true;
288+ // On the sphere lap_g = lap_s exactly, so the correction should compute
289+ // (numerically) close to zero regardless of niter -- not bit-identical
290+ // (it is a real computation now, through 8+ chained fp32 transforms per
291+ // iteration, not the literal `0 * Un` placeholder that used to make this
292+ // exact), but close. The tolerance is set by that chain's fp32 roundoff,
293+ // not by the scheme: a real geometry-correction bug would miss by orders
294+ // of magnitude more than this.
223295 let worst = 0;
224296 for (let k = 1; k < states.length; k++) {
225- if (states[k].length !== states[0].length) identical = false;
226297 for (let i = 0; i < states[0].length; i++) {
227- if (states[k][i] !== states[0][i]) identical = false;
228298 worst = Math.max(worst, Math.abs(states[k][i] - states[0][i]));
229299 }
230300 }
231301 check(
232- 'loop: the geometry correction is exactly zero, so the answer does not move',
233- identical,
234- identical
235- ? `bit-identical after ${STEPS} steps at ${counts.join('/')} iterations`
236- : `states differ by up to ${worst.toExponential(2)}`,
302+ 'loop: on the sphere, the correction stays near zero across niter',
303+ worst < 2e-3,
304+ `states differ by up to ${worst.toExponential(2)} after ${STEPS} steps at ${counts.join('/')} iterations`,
237305 );
238306 }
239307
308+ // ---- on a curved surface, the correction actually changes the answer ----
309+ {
310+ const model = mModelByKey('schnakenberg')!;
311+ const params = defaultParams(model);
312+ const peanut = mGeometryByKey('peanut')!;
313+ const peanutParams = defaultGeometryParams(peanut);
314+ // niter 0 vs 1 only -- deliberately not the 4/8 the sweep below already
315+ // documents as outside the Richardson iteration's convergence radius on
316+ // this geometry. The point here is just that the correction is not a
317+ // no-op, which a much smaller, still-converging niter already shows.
318+ const states: Float32Array[] = [];
319+ for (const niter of [0, 1]) {
320+ const session = await ModelSession.create({
321+ device, model, params, lmax: SWEEP_LMAX,
322+ geometry: peanut, geometryParams: peanutParams, niter,
323+ });
324+ session.seed(1);
325+ session.step(STEPS);
326+ states.push(await session.read('U'));
327+ session.destroy();
328+ }
329+ let worst = 0;
330+ for (let i = 0; i < states[0].length; i++) {
331+ worst = Math.max(worst, Math.abs(states[1][i] - states[0][i]));
332+ }
333+ check(
334+ 'loop: on peanut, the correction measurably changes the answer',
335+ worst > 1e-4 && states[1].every((v) => Number.isFinite(v)),
336+ `states differ by ${worst.toExponential(2)} after ${STEPS} steps at niter 0 vs 1`,
337+ );
338+ }
339+
340+ // ---- niter x geometry sweep: catch a "doesn't run" regression early -----
341+ // This is what actually turned up the two real issues found while building
342+ // the correction: peanut diverging at niter >= 4 with schnak-spots'
343+ // shipped default dt (a genuine Richardson-convergence-radius limit, not a
344+ // bug -- see docs/richardson-iteration.md), and a since-fixed compiler bug
345+ // where a loop-body statement could silently reuse a *different*
346+ // statement's compiled kernel (test/modelChecks.ts's pipeline-cache check
347+ // guards that one directly). Every shipped geometry x every niter the
348+ // app's <select> actually offers, so a regression anywhere in that grid is
349+ // caught -- without asserting away the one combination already known to be
350+ // outside the convergence radius.
351+ {
352+ const model = mModelByKey('schnakenberg')!;
353+ const params = defaultParams(model);
354+ const SWEEP_NITER = [0, 1, 2, 4, 8];
355+ const KNOWN_DIVERGENT = new Set(['peanut/2', 'peanut/4', 'peanut/8']);
356+
357+ for (const geomSpec of mGeometries) {
358+ for (const niter of SWEEP_NITER) {
359+ const session = await ModelSession.create({
360+ device, model, params, lmax: SWEEP_LMAX,
361+ geometry: geomSpec, geometryParams: defaultGeometryParams(geomSpec),
362+ niter,
363+ });
364+ session.seed(1);
365+ session.step(STEPS);
366+ const values = await session.read('u');
367+ const finite = values.every((v) => Number.isFinite(v));
368+ session.destroy();
369+
370+ const key = `${geomSpec.key}/${niter}`;
371+ const expectDivergent = KNOWN_DIVERGENT.has(key);
372+ check(
373+ expectDivergent
374+ ? `sweep: ${key} is known to diverge (outside the Richardson convergence radius)`
375+ : `sweep: ${key} stays finite after ${STEPS} steps`,
376+ expectDivergent ? !finite : finite,
377+ expectDivergent
378+ ? finite
379+ ? 'now finite -- the convergence radius may have improved; update KNOWN_DIVERGENT'
380+ : 'diverged as expected'
381+ : finite
382+ ? 'finite'
383+ : 'NOT FINITE -- unexpected divergence, investigate before treating this as another known case',
384+ );
385+ }
386+ }
387+ }
388+
240389 // ---- a loop whose length is not known at compile time is refused --------
241390 {
242391 const model = mModelByKey('allencahn')!;
test/modelChecks.tsmodified+16−9View file
@@ -31,13 +31,22 @@ const EXPECTED_KERNELS: Record<string, number> = {
3131 };
3232
3333 /**
34- * And what one unrolled iteration of the solve loop adds, per species: the
35- * placeholder line that will become the geometry correction, and the update
36- * that consumes it. Two rather than one because the correction does not fuse
37- * into its consumer — which is right, since the operator that replaces it will
38- * be transforms and kernels of its own, not an expression.
34+ * What one unrolled iteration of the solve loop adds, total (not per
35+ * species — the surface Laplace-Beltrami correction's per-species kernel
36+ * count is a byproduct of exactly how its expression tree happens to fuse,
37+ * not a clean per-species multiple, so this is measured per model rather
38+ * 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.
3944 */
40-const KERNELS_PER_ITERATION = 2;
45+const KERNELS_PER_ITERATION: Record<string, number> = {
46+ schnakenberg: 30,
47+ brusselator: 30,
48+ allencahn: 14,
49+};
4150
4251 const LMAX = 31;
4352 const STEPS = 40;
@@ -91,9 +100,7 @@ export async function modelChecks(
91100 const xforms = plan.step.filter(
92101 (l) => l.startsWith('synth') || l.startsWith('analys'),
93102 ).length;
94- const expected =
95- EXPECTED_KERNELS[model.key] +
96- NITER * KERNELS_PER_ITERATION * model.species.length;
103+ const expected = EXPECTED_KERNELS[model.key] + NITER * KERNELS_PER_ITERATION[model.key];
97104 log(
98105 ` ${model.key}.m -> ${plan.step.length} ops/step ` +
99106 `(${kernels} generated kernels, ${xforms} transforms, ${NITER} solve iter)`,
test/transformChecks.tsmodified+62−0View file
@@ -8,6 +8,7 @@
88 */
99 import { ShtPlan } from '../src/sht/sht.ts';
1010 import { ShtReference, randomSpectrum } from '../src/sht/reference.ts';
11+import { DerivPlan } from '../src/sht/deriv.ts';
1112 import { gridForLmax } from '../src/sht/layout.ts';
1213 import type { Check, Log } from './analyticChecks.ts';
1314
@@ -51,5 +52,66 @@ export async function transformChecks(
5152 `synth ${errSynth.toExponential(2)}, analys ${errAnalys.toExponential(2)}`,
5253 );
5354
55+ // ---- f64 reference dtheta/dphi vs an independent closed form -----------
56+ // x(theta,phi) = sin(theta)*cos(phi) is exactly degree 1, so quadrature
57+ // recovers it to f64 round-off; comparing its dtheta/dphi against the
58+ // grid-space analytic derivatives (not derived from the same recurrence
59+ // being tested) catches a sign or indexing error the random-spectrum check
60+ // below, which compares two implementations of the same formula, would not.
61+ {
62+ const x = new Float64Array(nlat * nphi);
63+ for (let i = 0; i < nlat; i++) {
64+ const st = ref.st[i];
65+ for (let j = 0; j < nphi; j++) {
66+ const phi = (2 * Math.PI * j) / nphi;
67+ x[i * nphi + j] = st * Math.cos(phi);
68+ }
69+ }
70+ const X = ref.analys(x);
71+ const dThetaX = ref.dtheta(X);
72+ const dPhiX = ref.dphi(X);
73+
74+ let errNum = 0;
75+ let norm = 0;
76+ for (let i = 0; i < nlat; i++) {
77+ const ct = ref.ct[i];
78+ const st = ref.st[i];
79+ for (let j = 0; j < nphi; j++) {
80+ const phi = (2 * Math.PI * j) / nphi;
81+ const k = i * nphi + j;
82+ const wantTheta = ct * Math.cos(phi);
83+ const wantPhi = -st * Math.sin(phi);
84+ errNum += (dThetaX[k] - wantTheta) ** 2 + (dPhiX[k] - wantPhi) ** 2;
85+ norm += wantTheta * wantTheta + wantPhi * wantPhi;
86+ }
87+ }
88+ const relErr = Math.sqrt(errNum / Math.max(norm, 1e-300));
89+ check(
90+ 'deriv: f64 reference dtheta/dphi match the closed form on x = sin(theta)cos(phi)',
91+ relErr < 1e-6,
92+ `rel L2 error ${relErr.toExponential(2)}`,
93+ );
94+ }
95+
96+ // ---- WGSL fp32 dtheta/dphi vs the (now closed-form-verified) f64 reference
97+ {
98+ const deriv = await DerivPlan.create(device, plan);
99+
100+ const dThetaGpu = await deriv.dtheta(new Float32Array(q64));
101+ const dThetaCpu = ref.dtheta(q64);
102+ const errDtheta = relL2(dThetaGpu, dThetaCpu);
103+
104+ const dPhiGpu = await deriv.dphi(new Float32Array(q64));
105+ const dPhiCpu = ref.dphi(q64);
106+ const errDphi = relL2(dPhiGpu, dPhiCpu);
107+
108+ check(
109+ 'deriv: WGSL fp32 dtheta/dphi vs f64 CPU reference',
110+ errDtheta < 1e-4 && errDphi < 1e-4,
111+ `dtheta ${errDtheta.toExponential(2)}, dphi ${errDphi.toExponential(2)}`,
112+ );
113+ deriv.destroy();
114+ }
115+
54116 plan.destroy();
55117 }