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: 3
314 // synths + 2 analyses (the flux-form matvec's five Legendre transforms,
315 // docs/reduced-transforms.md Sec 4 with the dphig variation) + the
316 // grid-space phi-derivative + 3 coefficient-space shuffles plus 7
317 // generated kernels -- see test/modelChecks.ts's KERNELS_PER_ITERATION,
318 // which counts the kernels alone; this counts every op.
319 const perIteration = ops[1] - ops[0];
320 const want = 32;
321 check(
322 'loop: unrolling is exactly linear in the trip count',
323 perIteration === want && ops[2] - ops[0] === 4 * perIteration,
324 `${perIteration} ops per iteration (expected ${want}), ` +
325 `${ops[2] - ops[0]} for 4 iterations`,
326 );
328 // On the sphere lap_g = lap_s exactly, so the correction should compute
329 // (numerically) close to zero regardless of niter -- not bit-identical
330 // (it is a real computation now, through 8+ chained fp32 transforms per
331 // iteration, not the literal `0 * Un` placeholder that used to make this
332 // exact), but close. The tolerance is set by that chain's fp32 roundoff,
333 // not by the scheme: a real geometry-correction bug would miss by orders
334 // of magnitude more than this.
335 let worst = 0;
336 for (let k = 1; k < states.length; k++) {
337 for (let i = 0; i < states[0].length; i++) {
338 worst = Math.max(worst, Math.abs(states[k][i] - states[0][i]));
339 }
340 }
341 check(
342 'loop: on the sphere, the correction stays near zero across niter',
343 worst < 2e-3,
344 `states differ by up to ${worst.toExponential(2)} after ${STEPS} steps at ${counts.join('/')} iterations`,
345 );
346 }
348 // ---- on a curved surface, the correction actually changes the answer ----
349 {
350 const model = mModelByKey('schnakenberg')!;
351 const params = defaultParams(model);
352 const peanut = mGeometryByKey('peanut')!;
353 const peanutParams = defaultGeometryParams(peanut);
354 // niter 0 vs 1 only -- deliberately not the 4/8 the sweep below already
355 // documents as outside the Richardson iteration's convergence radius on
356 // this geometry. The point here is just that the correction is not a
357 // no-op, which a much smaller, still-converging niter already shows.
358 const states: Float32Array[] = [];
359 for (const niter of [0, 1]) {
360 const session = await ModelSession.create({
361 device, model, params, lmax: SWEEP_LMAX,
362 geometry: peanut, geometryParams: peanutParams, niter,
363 });
364 await session.seed(1);
365 session.step(STEPS);
366 states.push(await session.read('U'));
367 session.destroy();
368 }
369 let worst = 0;
370 for (let i = 0; i < states[0].length; i++) {
371 worst = Math.max(worst, Math.abs(states[1][i] - states[0][i]));
372 }
373 check(
374 'loop: on peanut, the correction measurably changes the answer',
375 worst > 1e-4 && states[1].every((v) => Number.isFinite(v)),
376 `states differ by ${worst.toExponential(2)} after ${STEPS} steps at niter 0 vs 1`,
377 );
378 }
380 // ---- niter x geometry sweep: catch a "doesn't run" regression early -----
381 // This is what actually turned up the two real issues found while building
382 // the correction: peanut diverging at niter >= 2 with schnak-spots'
383 // shipped default dt (the plain round-sphere preconditioner's convergence
384 // radius -- the mean-J preconditioner has since lifted it; see the control
385 // check after the sweep), and a since-fixed compiler bug where a loop-body
386 // statement could silently reuse a *different* statement's compiled kernel
387 // (test/modelChecks.ts's pipeline-cache check guards that one directly).
388 // Every shipped geometry x every niter the app's <select> actually offers,
389 // so a regression anywhere in that grid is caught.
390 //
391 // SWEEP_LMAX stays at the app's default: the divergence the control check
392 // pins down is lmax-dependent (at lmax 31 or 15 even the plain
393 // preconditioner stays finite on peanut), so a smaller grid would stop
394 // testing the thing the control exists to demonstrate.
395 if (!runSweep) {
396 log(
397 ' sweep: skipped — run `npm run test:node` (desktop Dawn, ~3 s) or ' +
398 '`npm run test:gpu -- --sweep` for the niter x geometry sweep.',
399 );
400 } else {
401 const model = mModelByKey('schnakenberg')!;
402 const params = defaultParams(model);
403 const SWEEP_NITER = [0, 1, 2, 4, 8];
404 // Empty since the symbol-based preconditioner: its high-degree
405 // contraction rate (muMax - muMin)/(muMax + muMin) < 1 on any surface
406 // (mu = the symbol eigenvalues, i.e. inverse squared principal
407 // stretches), where the plain preconditioner diverges wherever mu > 2
408 // -- which is exactly what used to make peanut/2, /4 and /8 diverge.
409 // The mechanism stays: a regression lands here with its evidence.
410 const KNOWN_DIVERGENT = new Set<string>([]);
412 for (const geomSpec of mGeometries) {
413 for (const niter of SWEEP_NITER) {
414 const session = await ModelSession.create({
415 device, model, params, lmax: SWEEP_LMAX,
416 geometry: geomSpec, geometryParams: defaultGeometryParams(geomSpec),
417 niter,
418 });
419 await session.seed(1);
420 session.step(STEPS);
421 const values = await session.read('u');
422 const finite = values.every((v) => Number.isFinite(v));
423 session.destroy();
425 const key = `${geomSpec.key}/${niter}`;
426 const expectDivergent = KNOWN_DIVERGENT.has(key);
427 check(
428 expectDivergent
429 ? `sweep: ${key} is known to diverge (outside the Richardson convergence radius)`
430 : `sweep: ${key} stays finite after ${STEPS} steps`,
431 expectDivergent ? !finite : finite,
432 expectDivergent
433 ? finite
434 ? 'now finite -- the convergence radius may have improved; update KNOWN_DIVERGENT'
435 : 'diverged as expected'
436 : finite
437 ? 'finite'
438 : 'NOT FINITE -- unexpected divergence, investigate before treating this as another known case',
439 );
440 }
441 }
443 // ---- the mean-J control: what the sweep's health is owed to ----------
444 // peanut at niter 4 was the canonical divergent case before the mean-J
445 // preconditioner. Pinning jhat to 1 reproduces the plain round-sphere
446 // preconditioner on today's code, so this asserts both directions at
447 // once: mean-J converges where plain diverges, on the same operator,
448 // same surface, same dt. If this check ever finds jhat = 1 finite, the
449 // sweep above has stopped exercising the regime the preconditioner
450 // exists for (e.g. someone lowered SWEEP_LMAX or dt).
451 {
452 const peanut = mGeometryByKey('peanut')!;
453 const outcomes: boolean[] = [];
454 let jstats = '';
455 for (const jhat of [undefined, 1]) {
456 const session = await ModelSession.create({
457 device, model,
458 params: jhat === undefined ? params : { ...params, jhat },
459 lmax: SWEEP_LMAX,
460 geometry: peanut, geometryParams: defaultGeometryParams(peanut),
461 niter: 4,
462 });
463 if (jhat === undefined) {
464 const g = session.geometry;
465 jstats =
466 `mu in [${g.muMin.toFixed(3)}, ${g.muMax.toFixed(3)}] ` +
467 `(J in [${g.Jmin.toFixed(3)}, ${g.Jmax.toFixed(3)}]), ` +
468 `Jhat ${g.Jhat.toFixed(3)}, ` +
469 `rate ${((g.muMax - g.muMin) / (g.muMax + g.muMin)).toFixed(3)} ` +
470 `vs plain ${(g.muMax - 1).toFixed(2)}`;
471 }
472 await session.seed(1);
473 session.step(STEPS);
474 const values = await session.read('u');
475 outcomes.push(values.every((v) => Number.isFinite(v)));
476 session.destroy();
477 }
478 log(` mean-J on peanut: ${jstats}`);
479 check(
480 'mean-J: converges on peanut/4 where the plain preconditioner diverges',
481 outcomes[0] && !outcomes[1],
482 `mean-J finite: ${outcomes[0]}, jhat=1 finite: ${outcomes[1]}`,
483 );
484 }
485 }
487 // ---- a loop whose length is not known at compile time is refused --------
488 {
489 const model = mModelByKey('allencahn')!;
490 // `dt` is a tunable parameter, so it reaches the compiler with no value:
491 // the plan cannot know how many iterations to emit.
492 const bad = model.source.replace('for k = 1:niter', 'for k = 1:dt');
493 let message = '';
494 try {
495 const session = await ModelSession.create({
496 device, model, params: defaultParams(model), lmax: LMAX, source: bad, niter: 1,
497 });
498 session.destroy();
499 } catch (e) {
500 message = e instanceof ModelCompileError ? e.message : `wrong error type: ${e}`;
501 }
502 check(
503 'loop: a runtime loop bound is refused at compile time',
504 message.includes('known when the model is compiled'),
505 message ? `refused: ${message.slice(0, 72)}…` : 'compiled anyway',
506 );
507 }
509 // ---- swapping the surface leaves the simulation alone ------------------
510 {
511 const model = mModelByKey('schnakenberg')!;
512 const session = await ModelSession.create({
513 device, model, params: defaultParams(model), lmax: LMAX,
514 });
515 await session.seed(1);
516 session.step(STEPS);
517 const before = await session.read('U');
519 const peanut = mGeometryByKey('peanut')!;
520 await session.setGeometry(peanut, defaultGeometryParams(peanut));
521 const after = await session.read('U');
523 let survived = before.length === after.length;
524 for (let i = 0; survived && i < before.length; i++) {
525 if (before[i] !== after[i]) survived = false;
526 }
527 const { lo, hi } = session.geometry.radiusRange();
528 check(
529 'geometry: swapping the surface mid-run does not disturb the state',
530 survived && session.geometryModel.key === 'peanut' && hi - lo > 0.1,
531 survived
532 ? `state identical, now on ${session.geometryModel.key} (radius ${lo.toFixed(3)}–${hi.toFixed(3)})`
533 : 'state changed',
534 );
535 session.destroy();
536 }
538 await randnfun3Checks(device, check, log);
539}
541/**
542 * The seeded initial condition: chebfun's randnfun3, drawn on the host and
543 * summed on the GPU (src/mgpu/randnfun3.ts).
544 *
545 * The split is the thing worth testing. The draw is MATLAB whose distribution
546 * is checked directly, and the sum is a WGSL kernel checked against the same
547 * modes evaluated in f64 on the CPU — if the kernel's indexing into the packed
548 * mode table were wrong it would still produce a smooth random-looking field,
549 * which is exactly the kind of wrong no "looks patterned" check would catch.
550 */
551async function randnfun3Checks(
552 device: GPUDevice,
553 check: Check,
554 log: Log,
555): Promise<void> {
556 const model = mModelByKey('schnakenberg')!;
557 const params = defaultParams(model);
558 const make = (lam3: number): Promise<ModelSession> =>
559 ModelSession.create({ device, model, params, lmax: LMAX, lam3 });
561 // ---- the GPU sum matches the same modes evaluated on the CPU -----------
562 {
563 const session = await make(DEFAULT_LAMBDA);
564 await session.seed(3);
565 // `u` after init is the steady state plus 0.01*f, so the field is
566 // recovered by removing the model's own uniform offset.
567 const u = await session.read('u');
568 const g = session.geometry;
569 const modes = drawModes(
570 DEFAULT_LAMBDA,
571 boundingBox(g.x, g.y, g.z),
572 3,
573 g.x.length,
574 );
575 const nmodes = modes[0];
577 // The same sum in f64, straight from the packed table the GPU read.
578 let maxErr = 0;
579 let amp = 0;
580 const us = params.a + params.b;
581 for (let i = 0; i < g.x.length; i++) {
582 let f = 0;
583 for (let j = 0; j < nmodes; j++) {
584 const b = 4 + 5 * j;
585 const t = modes[b] * g.x[i] + modes[b + 1] * g.y[i] + modes[b + 2] * g.z[i];
586 f += modes[b + 3] * Math.cos(t) - modes[b + 4] * Math.sin(t);
587 }
588 const want = us + 0.01 * f;
589 maxErr = Math.max(maxErr, Math.abs(u[i] - want));
590 amp = Math.max(amp, Math.abs(0.01 * f));
591 }
592 log(` randnfun3: ${nmodes} modes at lambda ${DEFAULT_LAMBDA}, |perturbation| up to ${amp.toExponential(2)}`);
593 check(
594 'randnfun3: the GPU sum matches the same modes summed on the CPU',
595 // fp32 over ~1400 terms against f64, on a field of amplitude ~1e-2.
596 maxErr < 2e-6 && amp > 1e-3,
597 `max |GPU - CPU| = ${maxErr.toExponential(2)}, perturbation amplitude ${amp.toExponential(2)}`,
598 );
599 session.destroy();
600 }
602 // ---- a seed reproduces, a different seed does not ----------------------
603 {
604 const a = await make(DEFAULT_LAMBDA);
605 await a.seed(11);
606 const first = await a.read('u');
607 await a.seed(11);
608 const again = await a.read('u');
609 await a.seed(12);
610 const other = await a.read('u');
611 let same = true;
612 let differs = false;
613 for (let i = 0; i < first.length; i++) {
614 if (first[i] !== again[i]) same = false;
615 if (first[i] !== other[i]) differs = true;
616 }
617 check(
618 'randnfun3: the same seed redraws the same field, a different one does not',
619 same && differs,
620 same ? (differs ? 'reproducible and seed-dependent' : 'seed 12 gave seed 11 back') : 'not reproducible',
621 );
622 a.destroy();
623 }
625 // ---- the field is smooth, and lambda sets how smooth -------------------
626 //
627 // This is what randnfun3 buys over the white noise it replaced: the seed is
628 // band-limited, so it is fully resolved by the grid instead of being
629 // whatever the grid happened to alias. Measured as the share of spectral
630 // energy above degree 20 — near zero for a smooth field, and larger for a
631 // shorter wavelength, which is the direction lambda is supposed to move it.
632 {
633 const tail = async (lam3: number): Promise<number> => {
634 const session = await make(lam3);
635 await session.seed(5);
636 const U = await session.read('U');
637 let lo = 0;
638 let hi = 0;
639 for (let m = 0; m <= LMAX; m++) {
640 for (let l = m; l <= LMAX; l++) {
641 const i = lmIndex(LMAX, l, m);
642 const e = U[2 * i] ** 2 + U[2 * i + 1] ** 2;
643 if (l > 20) hi += e;
644 else lo += e;
645 }
646 }
647 session.destroy();
648 return hi / (lo + hi);
649 };
650 const coarse = await tail(1);
651 const fine = await tail(0.4);
652 log(` randnfun3: energy above l=20 is ${coarse.toExponential(2)} at lambda 1, ${fine.toExponential(2)} at lambda 0.4`);
653 check(
654 'randnfun3: the seed is band-limited, and lambda sets its scale',
655 coarse < 1e-3 && fine > coarse,
656 `tail ${coarse.toExponential(2)} (lambda 1) < ${fine.toExponential(2)} (lambda 0.4)`,
657 );
658 }
660 // ---- a finer wavelength grows the table rather than being capped -------
661 //
662 // The mode table is sized to the wavelength asked for, so going finer
663 // reallocates it and rebinds the dispatch. Getting that wrong would leave
664 // the kernel reading a destroyed buffer or a stale one, so check that a
665 // fine field is actually there and actually different.
666 {
667 const session = await make(DEFAULT_LAMBDA);
668 await session.seed(21);
669 const coarse = await session.read('u');
670 session.setLam3(0.12);
671 await session.seed(21);
672 const fine = await session.read('u');
673 let differs = false;
674 let finite = true;
675 for (let i = 0; i < fine.length; i++) {
676 if (!Number.isFinite(fine[i])) finite = false;
677 if (fine[i] !== coarse[i]) differs = true;
678 }
679 check(
680 'randnfun3: a finer wavelength grows the mode table and rebinds',
681 finite && differs,
682 finite ? 'redrew finer, buffer rebound' : 'field went non-finite after resize',
683 );
684 session.destroy();
685 }
687 // ---- a wavelength past the cost budget is refused, not truncated -------
688 {
689 const session = await make(DEFAULT_LAMBDA);
690 let message = '';
691 try {
692 session.setLam3(1e-4);
693 await session.seed(1);
694 } catch (e) {
695 message = e instanceof Error ? e.message : String(e);
696 }
697 check(
698 'randnfun3: a wavelength whose table could not be built is refused',
699 message.includes('Fourier modes on this surface'),
700 message ? `refused: ${message.slice(0, 62)}…` : 'drew it anyway',
701 );
702 session.destroy();
703 }
704}
706/** Index of the entry minimizing `score`, over the first `n` entries. */
707function argMin(
708 xs: Float64Array | Float32Array,
709 n: number,
710 score: (v: number) => number,
711): number {
712 let best = 0;
713 let bestScore = Infinity;
714 for (let i = 0; i < n; i++) {
715 const s = score(xs[i]);
716 if (s < bestScore) {
717 bestScore = s;
718 best = i;
719 }
720 }
721 return best;
722}