1/**
2 * The two things this project adds to turing-sphere: a surface, and a `for`
3 * loop in the compiled step.
4 *
5 * The surface is checked against what it is supposed to be — the sphere really
6 * is the unit sphere and really is degree 1, a deformed shape really has the
7 * radius profile its .m says, and the coefficients really do evaluate to the
8 * same surface on a finer grid.
9 *
10 * 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 * 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.
26 */
27import { ShtPlan } from '../src/sht/sht.ts';
28import { DerivPlan } from '../src/sht/deriv.ts';
29import { gridForLmax, lmIndex } from '../src/sht/layout.ts';
30import { ModelSession } from '../src/mgpu/session.ts';
31import { mModelByKey, defaultParams } from '../src/mgpu/registry.ts';
32import { Geometry } from '../src/geom/geometry.ts';
33import {
34 mGeometries,
35 mGeometryByKey,
36 defaultGeometryParams,
37 SPHERE_KEY,
38} from '../src/geom/registry.ts';
39import { ModelCompileError } from '../src/mgpu/errors.ts';
40import { boundingBox, drawModes, DEFAULT_LAMBDA } from '../src/mgpu/randnfun3.ts';
41import type { Check, Log } from './analyticChecks.ts';
43const LMAX = 31;
44const STEPS = 20;
45/** The app's actual default lmax (README: "at the default lmax 63 that is a
46 * 128x256 grid"), used for the niter/geometry sweep below and the peanut
47 * check next to it -- the divergence they're both about is a real, lmax-
48 * dependent numerical property of the Richardson iteration, not one this
49 * file's other, smaller LMAX happens to reproduce. */
50const SWEEP_LMAX = 63;
52/** Build one geometry on its own transform plan, for inspection. */
53async function buildGeometry(device: GPUDevice, key: string) {
54 const g = mGeometryByKey(key)!;
55 const { nlat, nphi } = gridForLmax(LMAX, 3);
56 const cfg = { lmax: LMAX, mmax: LMAX, nlat, nphi };
57 const sht = await ShtPlan.create(device, cfg);
58 const deriv = await DerivPlan.create(device, sht);
59 const geometry = await Geometry.create({
60 sht,
61 cfg,
62 source: g.source,
63 paramNames: g.params.map((p) => p.key),
64 params: defaultGeometryParams(g),
65 deriv,
66 });
67 return { g, sht, deriv, cfg, geometry };
68}
70export interface GeometryCheckOptions {
71 /**
72 * Run the niter x geometry sweep at the end. On by default, and nearly free
73 * on desktop Dawn (~3 s for all 20 combinations), but in a browser it
74 * dominates the whole suite: every session recompiles its unrolled step from
75 * scratch — there is no pipeline cache across sessions — so the sweep costs
76 * ~6 minutes on software WebGPU against ~30 s for every other check here put
77 * together. The browser page therefore leaves it out unless asked (?sweep=1),
78 * which is what keeps CI short.
79 */
80 sweep?: boolean;
81}
83export async function geometryChecks(
84 device: GPUDevice,
85 check: Check,
86 log: Log,
87 opts: GeometryCheckOptions = {},
88): Promise<void> {
89 const runSweep = opts.sweep ?? true;
90 // ---- every geometry compiles and closes ---------------------------------
91 for (const spec of mGeometries) {
92 const { sht, deriv, geometry } = await buildGeometry(device, spec.key);
93 let finite = true;
94 for (const a of [geometry.x, geometry.y, geometry.z]) {
95 for (const v of a) if (!Number.isFinite(v)) finite = false;
96 }
97 for (const a of [geometry.Vtx, geometry.Vty, geometry.Vtz, geometry.Vpx, geometry.Vpy, geometry.Vpz]) {
98 for (const v of a) if (!Number.isFinite(v)) finite = false;
99 }
100 const { lo, hi } = geometry.radiusRange();
101 check(
102 `geometry: ${spec.key}.m evaluates to a finite surface`,
103 finite && lo > 1e-3,
104 `radius ${lo.toFixed(4)}–${hi.toFixed(4)}`,
105 );
106 deriv.destroy();
107 sht.destroy();
108 }
110 // ---- the sphere is the unit sphere, exactly, and is degree 1 ------------
111 {
112 const { sht, deriv, geometry } = await buildGeometry(device, SPHERE_KEY);
114 let maxRadiusErr = 0;
115 for (let i = 0; i < geometry.x.length; i++) {
116 const r = Math.hypot(geometry.x[i], geometry.y[i], geometry.z[i]);
117 maxRadiusErr = Math.max(maxRadiusErr, Math.abs(r - 1));
118 }
119 // Tolerance is fp32 through a full analysis/synthesis round trip, not the
120 // geometry: the exact answer is representable, and what is measured here
121 // is the transforms' own round-off. It is set by the loosest stack this
122 // runs on — SwiftShader in CI is an order of magnitude worse than Dawn on
123 // real hardware (4e-4 against 2e-5). A geometry that was actually wrong
124 // would miss by O(1), so the slack costs nothing.
125 check(
126 'geometry: sphere.m has radius 1 everywhere',
127 maxRadiusErr < 2e-3,
128 `max |r - 1| = ${maxRadiusErr.toExponential(2)}`,
129 );
131 // x, y, z of the unit sphere are the three degree-1 harmonics and nothing
132 // else, so analysing them must leave every other coefficient at zero.
133 // This is what makes the sphere case exact rather than merely accurate:
134 // there is no content for the band limit to throw away.
135 const degreeOne = new Set([
136 lmIndex(LMAX, 1, 0),
137 lmIndex(LMAX, 1, 1),
138 ]);
139 let leak = 0;
140 for (const coeffs of [geometry.X, geometry.Y, geometry.Z]) {
141 for (let i = 0; i < coeffs.length / 2; i++) {
142 if (degreeOne.has(i)) continue;
143 leak = Math.max(leak, Math.abs(coeffs[2 * i]), Math.abs(coeffs[2 * i + 1]));
144 }
145 }
146 check(
147 'geometry: sphere.m is exactly degree 1 in the harmonics',
148 leak < 1e-3,
149 `max |coefficient| outside l = 1 is ${leak.toExponential(2)}`,
150 );
152 // The inverse metric quantities have a closed form on the unit sphere:
153 // V_theta = (cos(theta)cos(phi), cos(theta)sin(phi), -sin(theta)),
154 // V_phi = (-sin(phi)/sin(theta), cos(phi)/sin(theta), 0). Checking these
155 // pins the sign convention of computeMetric (src/geom/metric.ts) before
156 // it is buried under the Laplace-Beltrami operator built on top of it.
157 //
158 // Split by latitude, because the accuracy available here is not uniform.
159 // computeMetric divides by det = g_tt*g_pp - g_tp^2, and on the sphere
160 // g_pp and det are both O(sin^2 theta) -- 1.4e-3 at the outermost Gauss
161 // latitude of this grid. The transforms deliver X_theta/X_phi with an
162 // *absolute* fp32 error of ~1e-6, which is a large *relative* error once
163 // it is squared into quantities that small, so the error in both V's grows
164 // like 1/sin^2 theta toward the poles. That is conditioning, not a wrong
165 // formula: measured worst cases are
166 //
167 // sin(theta) >= 0.2 sin(theta) < 0.2 (8 of 64 latitudes)
168 // Dawn 3.2e-4 1.7e-3
169 // SwiftShader 1.5e-3 8.8e-3
170 //
171 // and a sign or formula error would be O(1) in either band, so tolerances
172 // a few times the looser stack still catch one.
173 const POLE_SIN = 0.2;
174 let maxMetricErr = 0;
175 let maxPoleErr = 0;
176 for (let i = 0; i < sht.cfg.nlat; i++) {
177 const ct = sht.cosTheta[i];
178 const st = Math.sqrt(Math.max(0, 1 - ct * ct));
179 for (let j = 0; j < sht.cfg.nphi; j++) {
180 const phi = (2 * Math.PI * j) / sht.cfg.nphi;
181 const k = i * sht.cfg.nphi + j;
182 const cphi = Math.cos(phi);
183 const sphi = Math.sin(phi);
184 const wantVtx = ct * cphi;
185 const wantVty = ct * sphi;
186 const wantVtz = -st;
187 const wantVpx = -sphi / st;
188 const wantVpy = cphi / st;
189 const wantVpz = 0;
190 const worst = Math.max(
191 Math.abs(geometry.Vtx[k] - wantVtx),
192 Math.abs(geometry.Vty[k] - wantVty),
193 Math.abs(geometry.Vtz[k] - wantVtz),
194 Math.abs(geometry.Vpx[k] - wantVpx),
195 Math.abs(geometry.Vpy[k] - wantVpy),
196 Math.abs(geometry.Vpz[k] - wantVpz),
197 );
198 if (st < POLE_SIN) maxPoleErr = Math.max(maxPoleErr, worst);
199 else maxMetricErr = Math.max(maxMetricErr, worst);
200 }
201 }
202 check(
203 'geometry: sphere.m has the closed-form inverse metric quantities',
204 maxMetricErr < 4e-3,
205 `max |V - closed form| = ${maxMetricErr.toExponential(2)} ` +
206 `away from the poles (sin theta >= ${POLE_SIN})`,
207 );
208 check(
209 'geometry: the polar caps stay within their conditioning',
210 maxPoleErr < 2e-2,
211 `max |V - closed form| = ${maxPoleErr.toExponential(2)} at sin theta < ${POLE_SIN}`,
212 );
214 deriv.destroy();
215 sht.destroy();
216 }
218 // ---- a deformed surface matches its own formula, on any grid ------------
219 {
220 const { g, sht, deriv, cfg, geometry } = await buildGeometry(device, 'peanut');
221 const p = defaultGeometryParams(g);
223 // peanut.m written out: r = 1 - waist*sin(theta)^2 scales the unit sphere,
224 // and z is then stretched, so the distance from the origin depends on
225 // theta alone. Checking every point against this closed form checks the
226 // whole path at once — the compiled shape kernel, the analysis into
227 // coefficients, the synthesis back — and, because the formula has no phi
228 // in it, that the surface really is a surface of revolution.
229 const peanutRadius = (ct: number): number => {
230 const st2 = Math.max(0, 1 - ct * ct);
231 const r = 1 - p.waist * st2;
232 return r * Math.hypot(Math.sqrt(st2), (1 + p.stretch) * ct);
233 };
235 const onGrid = (
236 cosTheta: Float64Array,
237 nlat: number,
238 nphi: number,
239 at: (i: number) => number,
240 ): number => {
241 let worst = 0;
242 for (let i = 0; i < nlat; i++) {
243 const want = peanutRadius(cosTheta[i]);
244 for (let j = 0; j < nphi; j++) {
245 worst = Math.max(worst, Math.abs(at(i * nphi + j) - want));
246 }
247 }
248 return worst;
249 };
251 const coarse = onGrid(sht.cosTheta, cfg.nlat, cfg.nphi, (k) =>
252 Math.hypot(geometry.x[k], geometry.y[k], geometry.z[k]),
253 );
254 check(
255 'geometry: peanut.m matches its own radial formula on the solver grid',
256 coarse < 1e-3,
257 `max |dr| = ${coarse.toExponential(2)}`,
258 );
260 // And the same on a finer grid, from the same coefficients. This is what
261 // "the rendered surface is the surface being solved on" means: display
262 // oversampling evaluates the embedding at more points, it does not
263 // subdivide or smooth it. The 2x Gauss latitudes share no point with the
264 // 1x ones, so agreeing here is agreeing everywhere, not at samples.
265 const fine = await ShtPlan.create(device, {
266 lmax: cfg.lmax,
267 mmax: cfg.mmax,
268 nlat: 2 * cfg.nlat,
269 nphi: 2 * cfg.nphi,
270 });
271 const finePos = await geometry.positionsOn(fine);
272 const refined = onGrid(fine.cosTheta, 2 * cfg.nlat, 2 * cfg.nphi, (k) =>
273 Math.hypot(finePos[3 * k], finePos[3 * k + 1], finePos[3 * k + 2]),
274 );
275 check(
276 'geometry: the same coefficients give the same surface on a 2x grid',
277 refined < 1e-3,
278 `max |dr| = ${refined.toExponential(2)} at ${2 * cfg.nlat}×${2 * cfg.nphi} points`,
279 );
280 fine.destroy();
281 deriv.destroy();
282 sht.destroy();
283 }
285 // ---- the unrolled loop: more ops, identical answer ----------------------
286 {
287 const model = mModelByKey('schnakenberg')!;
288 const params = defaultParams(model);
289 const counts = [0, 1, 4];
290 const ops: number[] = [];
291 const states: Float32Array[] = [];
293 for (const niter of counts) {
294 const session = await ModelSession.create({
295 device, model, params, lmax: LMAX, niter,
296 });
297 ops.push(session.describe().step.length);
298 await session.seed(1);
299 session.step(STEPS);
300 states.push(await session.read('U'));
301 session.destroy();
302 }
304 log(` schnakenberg.m ops/step by solve iterations: ${
305 counts.map((n, i) => `${n} -> ${ops[i]}`).join(', ')
306 }`);
307 check(
308 'loop: each solve iteration adds GPU operations',
309 ops[0] < ops[1] && ops[1] < ops[2],
310 `${ops.join(' < ')} ops for ${counts.join(', ')} iterations`,
311 );
312 // Unrolling has to be exactly linear in the trip count: the body planned
313 // once per iteration, no more and no less. Per species per iteration: 4
314 // synths + 2 analyses (the flux-form matvec's five Legendre transforms,
315 // docs/reduced-transforms.md Sec 4 with the dphig variation, plus the
316 // round-sphere synthesis of the divergence split) + the grid-space
317 // phi-derivative + 3 coefficient-space shuffles plus 8 generated kernels
318 // -- see test/modelChecks.ts's KERNELS_PER_ITERATION, which counts the
319 // kernels alone; this counts every op.
320 const perIteration = ops[1] - ops[0];
321 const want = 36;
322 check(
323 'loop: unrolling is exactly linear in the trip count',
324 perIteration === want && ops[2] - ops[0] === 4 * perIteration,
325 `${perIteration} ops per iteration (expected ${want}), ` +
326 `${ops[2] - ops[0]} for 4 iterations`,
327 );
329 // On the sphere lap_g = lap_s exactly, so the correction should compute
330 // (numerically) close to zero regardless of niter -- not bit-identical
331 // (it is a real computation now, through 8+ chained fp32 transforms per
332 // iteration, not the literal `0 * Un` placeholder that used to make this
333 // exact), but close. The tolerance is set by that chain's fp32 roundoff,
334 // not by the scheme: a real geometry-correction bug would miss by orders
335 // of magnitude more than this.
336 let worst = 0;
337 for (let k = 1; k < states.length; k++) {
338 for (let i = 0; i < states[0].length; i++) {
339 worst = Math.max(worst, Math.abs(states[k][i] - states[0][i]));
340 }
341 }
342 check(
343 'loop: on the sphere, the correction stays near zero across niter',
344 worst < 2e-3,
345 `states differ by up to ${worst.toExponential(2)} after ${STEPS} steps at ${counts.join('/')} iterations`,
346 );
347 }
349 // ---- on a curved surface, the correction actually changes the answer ----
350 {
351 const model = mModelByKey('schnakenberg')!;
352 const params = defaultParams(model);
353 const peanut = mGeometryByKey('peanut')!;
354 const peanutParams = defaultGeometryParams(peanut);
355 // niter 0 vs 1 only -- deliberately not the 4/8 the sweep below already
356 // documents as outside the Richardson iteration's convergence radius on
357 // this geometry. The point here is just that the correction is not a
358 // no-op, which a much smaller, still-converging niter already shows.
359 const states: Float32Array[] = [];
360 for (const niter of [0, 1]) {
361 const session = await ModelSession.create({
362 device, model, params, lmax: SWEEP_LMAX,
363 geometry: peanut, geometryParams: peanutParams, niter,
364 });
365 await session.seed(1);
366 session.step(STEPS);
367 states.push(await session.read('U'));
368 session.destroy();
369 }
370 let worst = 0;
371 for (let i = 0; i < states[0].length; i++) {
372 worst = Math.max(worst, Math.abs(states[1][i] - states[0][i]));
373 }
374 check(
375 'loop: on peanut, the correction measurably changes the answer',
376 worst > 1e-4 && states[1].every((v) => Number.isFinite(v)),
377 `states differ by ${worst.toExponential(2)} after ${STEPS} steps at niter 0 vs 1`,
378 );
379 }
381 // ---- niter x geometry sweep: catch a "doesn't run" regression early -----
382 // This is what actually turned up the two real issues found while building
383 // the correction: peanut diverging at niter >= 2 with schnak-spots'
384 // shipped default dt (the plain round-sphere preconditioner's convergence
385 // radius -- the mean-J preconditioner has since lifted it; see the control
386 // check after the sweep), and a since-fixed compiler bug where a loop-body
387 // statement could silently reuse a *different* statement's compiled kernel
388 // (test/modelChecks.ts's pipeline-cache check guards that one directly).
389 // Every shipped geometry x every niter the app's <select> actually offers,
390 // so a regression anywhere in that grid is caught.
391 //
392 // SWEEP_LMAX stays at the app's default: the divergence the control check
393 // pins down is lmax-dependent (at lmax 31 or 15 even the plain
394 // preconditioner stays finite on peanut), so a smaller grid would stop
395 // testing the thing the control exists to demonstrate.
396 if (!runSweep) {
397 log(
398 ' sweep: skipped — run `npm run test:node` (desktop Dawn, ~3 s) or ' +
399 '`npm run test:gpu -- --sweep` for the niter x geometry sweep.',
400 );
401 } else {
402 const model = mModelByKey('schnakenberg')!;
403 const params = defaultParams(model);
404 const SWEEP_NITER = [0, 1, 2, 4, 8];
405 // Empty since the symbol-based preconditioner: its high-degree
406 // contraction rate (muMax - muMin)/(muMax + muMin) < 1 on any surface
407 // (mu = the symbol eigenvalues, i.e. inverse squared principal
408 // stretches), where the plain preconditioner diverges wherever mu > 2
409 // -- which is exactly what used to make peanut/2, /4 and /8 diverge.
410 // The mechanism stays: a regression lands here with its evidence.
411 const KNOWN_DIVERGENT = new Set<string>([]);
413 for (const geomSpec of mGeometries) {
414 for (const niter of SWEEP_NITER) {
415 const session = await ModelSession.create({
416 device, model, params, lmax: SWEEP_LMAX,
417 geometry: geomSpec, geometryParams: defaultGeometryParams(geomSpec),
418 niter,
419 });
420 await session.seed(1);
421 session.step(STEPS);
422 const values = await session.read('u');
423 const finite = values.every((v) => Number.isFinite(v));
424 session.destroy();
426 const key = `${geomSpec.key}/${niter}`;
427 const expectDivergent = KNOWN_DIVERGENT.has(key);
428 check(
429 expectDivergent
430 ? `sweep: ${key} is known to diverge (outside the Richardson convergence radius)`
431 : `sweep: ${key} stays finite after ${STEPS} steps`,
432 expectDivergent ? !finite : finite,
433 expectDivergent
434 ? finite
435 ? 'now finite -- the convergence radius may have improved; update KNOWN_DIVERGENT'
436 : 'diverged as expected'
437 : finite
438 ? 'finite'
439 : 'NOT FINITE -- unexpected divergence, investigate before treating this as another known case',
440 );
441 }
442 }
444 // ---- the mean-J control: what the sweep's health is owed to ----------
445 // peanut at niter 4 was the canonical divergent case before the mean-J
446 // preconditioner. Pinning jhat to 1 reproduces the plain round-sphere
447 // preconditioner on today's code, so this asserts both directions at
448 // once: mean-J converges where plain diverges, on the same operator,
449 // same surface, same dt. If this check ever finds jhat = 1 finite, the
450 // sweep above has stopped exercising the regime the preconditioner
451 // exists for (e.g. someone lowered SWEEP_LMAX or dt).
452 {
453 const peanut = mGeometryByKey('peanut')!;
454 const outcomes: boolean[] = [];
455 let jstats = '';
456 for (const jhat of [undefined, 1]) {
457 const session = await ModelSession.create({
458 device, model,
459 params: jhat === undefined ? params : { ...params, jhat },
460 lmax: SWEEP_LMAX,
461 geometry: peanut, geometryParams: defaultGeometryParams(peanut),
462 niter: 4,
463 });
464 if (jhat === undefined) {
465 const g = session.geometry;
466 jstats =
467 `mu in [${g.muMin.toFixed(3)}, ${g.muMax.toFixed(3)}] ` +
468 `(J in [${g.Jmin.toFixed(3)}, ${g.Jmax.toFixed(3)}]), ` +
469 `Jhat ${g.Jhat.toFixed(3)}, ` +
470 `rate ${((g.muMax - g.muMin) / (g.muMax + g.muMin)).toFixed(3)} ` +
471 `vs plain ${(g.muMax - 1).toFixed(2)}`;
472 }
473 await session.seed(1);
474 session.step(STEPS);
475 const values = await session.read('u');
476 outcomes.push(values.every((v) => Number.isFinite(v)));
477 session.destroy();
478 }
479 log(` mean-J on peanut: ${jstats}`);
480 check(
481 'mean-J: converges on peanut/4 where the plain preconditioner diverges',
482 outcomes[0] && !outcomes[1],
483 `mean-J finite: ${outcomes[0]}, jhat=1 finite: ${outcomes[1]}`,
484 );
485 }
486 }
488 // ---- a loop whose length is not known at compile time is refused --------
489 {
490 const model = mModelByKey('allencahn')!;
491 // `dt` is a tunable parameter, so it reaches the compiler with no value:
492 // the plan cannot know how many iterations to emit.
493 const bad = model.source.replace('for k = 1:niter', 'for k = 1:dt');
494 let message = '';
495 try {
496 const session = await ModelSession.create({
497 device, model, params: defaultParams(model), lmax: LMAX, source: bad, niter: 1,
498 });
499 session.destroy();
500 } catch (e) {
501 message = e instanceof ModelCompileError ? e.message : `wrong error type: ${e}`;
502 }
503 check(
504 'loop: a runtime loop bound is refused at compile time',
505 message.includes('known when the model is compiled'),
506 message ? `refused: ${message.slice(0, 72)}…` : 'compiled anyway',
507 );
508 }
510 // ---- swapping the surface leaves the simulation alone ------------------
511 {
512 const model = mModelByKey('schnakenberg')!;
513 const session = await ModelSession.create({
514 device, model, params: defaultParams(model), lmax: LMAX,
515 });
516 await session.seed(1);
517 session.step(STEPS);
518 const before = await session.read('U');
520 const peanut = mGeometryByKey('peanut')!;
521 await session.setGeometry(peanut, defaultGeometryParams(peanut));
522 const after = await session.read('U');
524 let survived = before.length === after.length;
525 for (let i = 0; survived && i < before.length; i++) {
526 if (before[i] !== after[i]) survived = false;
527 }
528 const { lo, hi } = session.geometry.radiusRange();
529 check(
530 'geometry: swapping the surface mid-run does not disturb the state',
531 survived && session.geometryModel.key === 'peanut' && hi - lo > 0.1,
532 survived
533 ? `state identical, now on ${session.geometryModel.key} (radius ${lo.toFixed(3)}–${hi.toFixed(3)})`
534 : 'state changed',
535 );
536 session.destroy();
537 }
539 await randnfun3Checks(device, check, log);
540}
542/**
543 * The seeded initial condition: chebfun's randnfun3, drawn on the host and
544 * summed on the GPU (src/mgpu/randnfun3.ts).
545 *
546 * The split is the thing worth testing. The draw is MATLAB whose distribution
547 * is checked directly, and the sum is a WGSL kernel checked against the same
548 * modes evaluated in f64 on the CPU — if the kernel's indexing into the packed
549 * mode table were wrong it would still produce a smooth random-looking field,
550 * which is exactly the kind of wrong no "looks patterned" check would catch.
551 */
552async function randnfun3Checks(
553 device: GPUDevice,
554 check: Check,
555 log: Log,
556): Promise<void> {
557 const model = mModelByKey('schnakenberg')!;
558 const params = defaultParams(model);
559 const make = (lam3: number): Promise<ModelSession> =>
560 ModelSession.create({ device, model, params, lmax: LMAX, lam3 });
562 // ---- the GPU sum matches the same modes evaluated on the CPU -----------
563 {
564 const session = await make(DEFAULT_LAMBDA);
565 await session.seed(3);
566 // `u` after init is the steady state plus 0.01*f, so the field is
567 // recovered by removing the model's own uniform offset.
568 const u = await session.read('u');
569 const g = session.geometry;
570 const modes = drawModes(
571 DEFAULT_LAMBDA,
572 boundingBox(g.x, g.y, g.z),
573 3,
574 g.x.length,
575 );
576 const nmodes = modes[0];
578 // The same sum in f64, straight from the packed table the GPU read.
579 let maxErr = 0;
580 let amp = 0;
581 const us = params.a + params.b;
582 for (let i = 0; i < g.x.length; i++) {
583 let f = 0;
584 for (let j = 0; j < nmodes; j++) {
585 const b = 4 + 5 * j;
586 const t = modes[b] * g.x[i] + modes[b + 1] * g.y[i] + modes[b + 2] * g.z[i];
587 f += modes[b + 3] * Math.cos(t) - modes[b + 4] * Math.sin(t);
588 }
589 const want = us + 0.01 * f;
590 maxErr = Math.max(maxErr, Math.abs(u[i] - want));
591 amp = Math.max(amp, Math.abs(0.01 * f));
592 }
593 log(` randnfun3: ${nmodes} modes at lambda ${DEFAULT_LAMBDA}, |perturbation| up to ${amp.toExponential(2)}`);
594 check(
595 'randnfun3: the GPU sum matches the same modes summed on the CPU',
596 // fp32 over ~1400 terms against f64. What is being bounded is the
597 // summation floor, and its size is the backend's accumulation order:
598 // Metal lands at 2.0e-6, SwiftShader at 3.8e-6, so an absolute constant
599 // tuned on one is a coin flip on the other. Scale it to the field
600 // instead. The bug this exists to catch -- a mis-indexed read into the
601 // packed table, which would still look like a smooth random field -- is
602 // wrong by O(amp), a thousand times over the bound.
603 maxErr < 1e-3 * amp && amp > 1e-3,
604 `max |GPU - CPU| = ${maxErr.toExponential(2)}, perturbation amplitude ${amp.toExponential(2)}`,
605 );
606 session.destroy();
607 }
609 // ---- a seed reproduces, a different seed does not ----------------------
610 {
611 const a = await make(DEFAULT_LAMBDA);
612 await a.seed(11);
613 const first = await a.read('u');
614 await a.seed(11);
615 const again = await a.read('u');
616 await a.seed(12);
617 const other = await a.read('u');
618 let same = true;
619 let differs = false;
620 for (let i = 0; i < first.length; i++) {
621 if (first[i] !== again[i]) same = false;
622 if (first[i] !== other[i]) differs = true;
623 }
624 check(
625 'randnfun3: the same seed redraws the same field, a different one does not',
626 same && differs,
627 same ? (differs ? 'reproducible and seed-dependent' : 'seed 12 gave seed 11 back') : 'not reproducible',
628 );
629 a.destroy();
630 }
632 // ---- the field is smooth, and lambda sets how smooth -------------------
633 //
634 // This is what randnfun3 buys over the white noise it replaced: the seed is
635 // band-limited, so it is fully resolved by the grid instead of being
636 // whatever the grid happened to alias. Measured as the share of spectral
637 // energy above degree 20 — near zero for a smooth field, and larger for a
638 // shorter wavelength, which is the direction lambda is supposed to move it.
639 {
640 const tail = async (lam3: number): Promise<number> => {
641 const session = await make(lam3);
642 await session.seed(5);
643 const U = await session.read('U');
644 let lo = 0;
645 let hi = 0;
646 for (let m = 0; m <= LMAX; m++) {
647 for (let l = m; l <= LMAX; l++) {
648 const i = lmIndex(LMAX, l, m);
649 const e = U[2 * i] ** 2 + U[2 * i + 1] ** 2;
650 if (l > 20) hi += e;
651 else lo += e;
652 }
653 }
654 session.destroy();
655 return hi / (lo + hi);
656 };
657 const coarse = await tail(1);
658 const fine = await tail(0.4);
659 log(` randnfun3: energy above l=20 is ${coarse.toExponential(2)} at lambda 1, ${fine.toExponential(2)} at lambda 0.4`);
660 check(
661 'randnfun3: the seed is band-limited, and lambda sets its scale',
662 coarse < 1e-3 && fine > coarse,
663 `tail ${coarse.toExponential(2)} (lambda 1) < ${fine.toExponential(2)} (lambda 0.4)`,
664 );
665 }
667 // ---- a finer wavelength grows the table rather than being capped -------
668 //
669 // The mode table is sized to the wavelength asked for, so going finer
670 // reallocates it and rebinds the dispatch. Getting that wrong would leave
671 // the kernel reading a destroyed buffer or a stale one, so check that a
672 // fine field is actually there and actually different.
673 {
674 const session = await make(DEFAULT_LAMBDA);
675 await session.seed(21);
676 const coarse = await session.read('u');
677 session.setLam3(0.12);
678 await session.seed(21);
679 const fine = await session.read('u');
680 let differs = false;
681 let finite = true;
682 for (let i = 0; i < fine.length; i++) {
683 if (!Number.isFinite(fine[i])) finite = false;
684 if (fine[i] !== coarse[i]) differs = true;
685 }
686 check(
687 'randnfun3: a finer wavelength grows the mode table and rebinds',
688 finite && differs,
689 finite ? 'redrew finer, buffer rebound' : 'field went non-finite after resize',
690 );
691 session.destroy();
692 }
694 // ---- a wavelength past the cost budget is refused, not truncated -------
695 {
696 const session = await make(DEFAULT_LAMBDA);
697 let message = '';
698 try {
699 session.setLam3(1e-4);
700 await session.seed(1);
701 } catch (e) {
702 message = e instanceof Error ? e.message : String(e);
703 }
704 check(
705 'randnfun3: a wavelength whose table could not be built is refused',
706 message.includes('Fourier modes on this surface'),
707 message ? `refused: ${message.slice(0, 62)}…` : 'drew it anyway',
708 );
709 session.destroy();
710 }
711}
713/** Index of the entry minimizing `score`, over the first `n` entries. */
714function argMin(
715 xs: Float64Array | Float32Array,
716 n: number,
717 score: (v: number) => number,
718): number {
719 let best = 0;
720 let bestScore = Infinity;
721 for (let i = 0; i < n; i++) {
722 const s = score(xs[i]);
723 if (s < bestScore) {
724 bestScore = s;
725 best = i;
726 }
727 }
728 return best;
729}