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 { relL2 } from '../src/mgpu/digest.ts';
41import type { Check, Log } from './analyticChecks.ts';
43const LMAX = 31;
44const STEPS = 20;
46/**
47 * The shipped Schnakenberg model with its solver calls switched from
48 * richardson to another shipped solver — the same edit a user makes in the
49 * page, which is the point: same operator, same model, different solver.
50 * (gmres additionally takes `nlm`, for sizing its Krylov basis bank.)
51 * Throws if the model text drifted from what this rewrites, so the test
52 * fails loudly rather than silently comparing richardson with itself.
53 */
54function solverSchnak(source: string, solver: 'bicgstab' | 'gmres'): string {
55 const out = source
56 .replace('step(U, V, lam, filt, ', 'step(U, V, lam, filt, wlm, ')
57 .replace('D2, dt, niter)', solver === 'gmres' ? 'D2, dt, nlm, niter)' : 'D2, dt, niter)')
58 .replaceAll('richardson(Bu, dt * D1, lam, filt, ', `${solver}(Bu, dt * D1, lam, filt, wlm, `)
59 .replaceAll('richardson(Bv, dt * D2, lam, filt, ', `${solver}(Bv, dt * D2, lam, filt, wlm, `)
60 .replaceAll('Vpz, niter);', solver === 'gmres' ? 'Vpz, nlm, niter);' : 'Vpz, niter);');
61 if (!out.includes('wlm,') || !out.includes(`${solver}(Bu`) || !out.includes(`${solver}(Bv`)) {
62 throw new Error('solverSchnak: the model source no longer matches the rewrite');
63 }
64 return out;
65}
66/** The app's actual default lmax (README: "at the default lmax 63 that is a
67 * 128x256 grid"), used for the niter/geometry sweep below and the peanut
68 * check next to it -- the divergence they're both about is a real, lmax-
69 * dependent numerical property of the Richardson iteration, not one this
70 * file's other, smaller LMAX happens to reproduce. */
71const SWEEP_LMAX = 63;
73/** Build one geometry on its own transform plan, for inspection. */
74async function buildGeometry(device: GPUDevice, key: string) {
75 const g = mGeometryByKey(key)!;
76 const { nlat, nphi } = gridForLmax(LMAX, 3);
77 const cfg = { lmax: LMAX, mmax: LMAX, nlat, nphi };
78 const sht = await ShtPlan.create(device, cfg);
79 const deriv = await DerivPlan.create(device, sht);
80 const geometry = await Geometry.create({
81 device,
82 sht,
83 cfg,
84 source: g.source,
85 paramNames: g.params.map((p) => p.key),
86 params: defaultGeometryParams(g),
87 deriv,
88 });
89 return { g, sht, deriv, cfg, geometry };
90}
92export interface GeometryCheckOptions {
93 /**
94 * Run the niter x geometry sweep at the end. On by default, and nearly free
95 * on desktop Dawn (~3 s for all 20 combinations), but in a browser it
96 * dominates the whole suite: every session recompiles its unrolled step from
97 * scratch — there is no pipeline cache across sessions — so the sweep costs
98 * ~6 minutes on software WebGPU against ~30 s for every other check here put
99 * together. The browser page therefore leaves it out unless asked (?sweep=1),
100 * which is what keeps CI short.
101 */
102 sweep?: boolean;
103}
105export async function geometryChecks(
106 device: GPUDevice,
107 check: Check,
108 log: Log,
109 opts: GeometryCheckOptions = {},
110): Promise<void> {
111 const runSweep = opts.sweep ?? true;
112 // ---- every geometry compiles and closes ---------------------------------
113 for (const spec of mGeometries) {
114 const { sht, deriv, geometry } = await buildGeometry(device, spec.key);
115 let finite = true;
116 for (const a of [geometry.x, geometry.y, geometry.z]) {
117 for (const v of a) if (!Number.isFinite(v)) finite = false;
118 }
119 for (const a of [geometry.Vtx, geometry.Vty, geometry.Vtz, geometry.Vpx, geometry.Vpy, geometry.Vpz]) {
120 for (const v of a) if (!Number.isFinite(v)) finite = false;
121 }
122 const { lo, hi } = geometry.radiusRange();
123 check(
124 `geometry: ${spec.key}.m evaluates to a finite surface`,
125 finite && lo > 1e-3,
126 `radius ${lo.toFixed(4)}–${hi.toFixed(4)}`,
127 );
128 deriv.destroy();
129 sht.destroy();
130 }
132 // ---- the sphere is the unit sphere, exactly, and is degree 1 ------------
133 {
134 const { sht, deriv, geometry } = await buildGeometry(device, SPHERE_KEY);
136 let maxRadiusErr = 0;
137 for (let i = 0; i < geometry.x.length; i++) {
138 const r = Math.hypot(geometry.x[i], geometry.y[i], geometry.z[i]);
139 maxRadiusErr = Math.max(maxRadiusErr, Math.abs(r - 1));
140 }
141 // Tolerance is fp32 through a full analysis/synthesis round trip, not the
142 // geometry: the exact answer is representable, and what is measured here
143 // is the transforms' own round-off. It is set by the loosest stack this
144 // runs on — SwiftShader in CI is an order of magnitude worse than Dawn on
145 // real hardware (4e-4 against 2e-5). A geometry that was actually wrong
146 // would miss by O(1), so the slack costs nothing.
147 check(
148 'geometry: sphere.m has radius 1 everywhere',
149 maxRadiusErr < 2e-3,
150 `max |r - 1| = ${maxRadiusErr.toExponential(2)}`,
151 );
153 // x, y, z of the unit sphere are the three degree-1 harmonics and nothing
154 // else, so analysing them must leave every other coefficient at zero.
155 // This is what makes the sphere case exact rather than merely accurate:
156 // there is no content for the band limit to throw away.
157 const degreeOne = new Set([
158 lmIndex(LMAX, 1, 0),
159 lmIndex(LMAX, 1, 1),
160 ]);
161 let leak = 0;
162 for (const coeffs of [geometry.X, geometry.Y, geometry.Z]) {
163 for (let i = 0; i < coeffs.length / 2; i++) {
164 if (degreeOne.has(i)) continue;
165 leak = Math.max(leak, Math.abs(coeffs[2 * i]), Math.abs(coeffs[2 * i + 1]));
166 }
167 }
168 check(
169 'geometry: sphere.m is exactly degree 1 in the harmonics',
170 leak < 1e-3,
171 `max |coefficient| outside l = 1 is ${leak.toExponential(2)}`,
172 );
174 // The inverse metric quantities have a closed form on the unit sphere:
175 // V_theta = (cos(theta)cos(phi), cos(theta)sin(phi), -sin(theta)),
176 // V_phi = (-sin(phi)/sin(theta), cos(phi)/sin(theta), 0). Checking these
177 // pins the sign convention of computeMetric (src/geom/metric.ts) before
178 // it is buried under the Laplace-Beltrami operator built on top of it.
179 //
180 // Split by latitude, because the accuracy available here is not uniform.
181 // computeMetric divides by det = g_tt*g_pp - g_tp^2, and on the sphere
182 // g_pp and det are both O(sin^2 theta) -- 1.4e-3 at the outermost Gauss
183 // latitude of this grid. The transforms deliver X_theta/X_phi with an
184 // *absolute* fp32 error of ~1e-6, which is a large *relative* error once
185 // it is squared into quantities that small, so the error in both V's grows
186 // like 1/sin^2 theta toward the poles. That is conditioning, not a wrong
187 // formula: measured worst cases are
188 //
189 // sin(theta) >= 0.2 sin(theta) < 0.2 (8 of 64 latitudes)
190 // Dawn 3.2e-4 1.7e-3
191 // SwiftShader 1.5e-3 8.8e-3
192 //
193 // and a sign or formula error would be O(1) in either band, so tolerances
194 // a few times the looser stack still catch one.
195 const POLE_SIN = 0.2;
196 let maxMetricErr = 0;
197 let maxPoleErr = 0;
198 for (let i = 0; i < sht.cfg.nlat; i++) {
199 const ct = sht.cosTheta[i];
200 const st = Math.sqrt(Math.max(0, 1 - ct * ct));
201 for (let j = 0; j < sht.cfg.nphi; j++) {
202 const phi = (2 * Math.PI * j) / sht.cfg.nphi;
203 const k = i * sht.cfg.nphi + j;
204 const cphi = Math.cos(phi);
205 const sphi = Math.sin(phi);
206 const wantVtx = ct * cphi;
207 const wantVty = ct * sphi;
208 const wantVtz = -st;
209 const wantVpx = -sphi / st;
210 const wantVpy = cphi / st;
211 const wantVpz = 0;
212 const worst = Math.max(
213 Math.abs(geometry.Vtx[k] - wantVtx),
214 Math.abs(geometry.Vty[k] - wantVty),
215 Math.abs(geometry.Vtz[k] - wantVtz),
216 Math.abs(geometry.Vpx[k] - wantVpx),
217 Math.abs(geometry.Vpy[k] - wantVpy),
218 Math.abs(geometry.Vpz[k] - wantVpz),
219 );
220 if (st < POLE_SIN) maxPoleErr = Math.max(maxPoleErr, worst);
221 else maxMetricErr = Math.max(maxMetricErr, worst);
222 }
223 }
224 check(
225 'geometry: sphere.m has the closed-form inverse metric quantities',
226 maxMetricErr < 4e-3,
227 `max |V - closed form| = ${maxMetricErr.toExponential(2)} ` +
228 `away from the poles (sin theta >= ${POLE_SIN})`,
229 );
230 check(
231 'geometry: the polar caps stay within their conditioning',
232 maxPoleErr < 2e-2,
233 `max |V - closed form| = ${maxPoleErr.toExponential(2)} at sin theta < ${POLE_SIN}`,
234 );
236 deriv.destroy();
237 sht.destroy();
238 }
240 // ---- a deformed surface matches its own formula, on any grid ------------
241 {
242 const { g, sht, deriv, cfg, geometry } = await buildGeometry(device, 'peanut');
243 const p = defaultGeometryParams(g);
245 // peanut.m written out: r = 1 - waist*sin(theta)^2 scales the unit sphere,
246 // and z is then stretched, so the distance from the origin depends on
247 // theta alone. Checking every point against this closed form checks the
248 // whole path at once — the compiled shape kernel, the analysis into
249 // coefficients, the synthesis back — and, because the formula has no phi
250 // in it, that the surface really is a surface of revolution.
251 const peanutRadius = (ct: number): number => {
252 const st2 = Math.max(0, 1 - ct * ct);
253 const r = 1 - p.waist * st2;
254 return r * Math.hypot(Math.sqrt(st2), (1 + p.stretch) * ct);
255 };
257 const onGrid = (
258 cosTheta: Float64Array,
259 nlat: number,
260 nphi: number,
261 at: (i: number) => number,
262 ): number => {
263 let worst = 0;
264 for (let i = 0; i < nlat; i++) {
265 const want = peanutRadius(cosTheta[i]);
266 for (let j = 0; j < nphi; j++) {
267 worst = Math.max(worst, Math.abs(at(i * nphi + j) - want));
268 }
269 }
270 return worst;
271 };
273 const coarse = onGrid(sht.cosTheta, cfg.nlat, cfg.nphi, (k) =>
274 Math.hypot(geometry.x[k], geometry.y[k], geometry.z[k]),
275 );
276 check(
277 'geometry: peanut.m matches its own radial formula on the solver grid',
278 coarse < 1e-3,
279 `max |dr| = ${coarse.toExponential(2)}`,
280 );
282 // And the same on a finer grid, from the same coefficients. This is what
283 // "the rendered surface is the surface being solved on" means: display
284 // oversampling evaluates the embedding at more points, it does not
285 // subdivide or smooth it. The 2x Gauss latitudes share no point with the
286 // 1x ones, so agreeing here is agreeing everywhere, not at samples.
287 const fine = await ShtPlan.create(device, {
288 lmax: cfg.lmax,
289 mmax: cfg.mmax,
290 nlat: 2 * cfg.nlat,
291 nphi: 2 * cfg.nphi,
292 });
293 const finePos = await geometry.positionsOn(fine);
294 const refined = onGrid(fine.cosTheta, 2 * cfg.nlat, 2 * cfg.nphi, (k) =>
295 Math.hypot(finePos[3 * k], finePos[3 * k + 1], finePos[3 * k + 2]),
296 );
297 check(
298 'geometry: the same coefficients give the same surface on a 2x grid',
299 refined < 1e-3,
300 `max |dr| = ${refined.toExponential(2)} at ${2 * cfg.nlat}×${2 * cfg.nphi} points`,
301 );
302 fine.destroy();
303 deriv.destroy();
304 sht.destroy();
305 }
307 // ---- the unrolled loop: more ops, identical answer ----------------------
308 {
309 const model = mModelByKey('schnakenberg')!;
310 const params = defaultParams(model);
311 const counts = [0, 1, 4];
312 const ops: number[] = [];
313 const states: Float32Array[] = [];
315 for (const niter of counts) {
316 const session = await ModelSession.create({
317 device, model, params, lmax: LMAX, niter,
318 });
319 ops.push(session.describe().step.length);
320 session.seed(1);
321 session.step(STEPS);
322 states.push(await session.read('U'));
323 session.destroy();
324 }
326 log(` schnakenberg.m ops/step by solve iterations: ${
327 counts.map((n, i) => `${n} -> ${ops[i]}`).join(', ')
328 }`);
329 check(
330 'loop: each solve iteration adds GPU operations',
331 ops[0] < ops[1] && ops[1] < ops[2],
332 `${ops.join(' < ')} ops for ${counts.join(', ')} iterations`,
333 );
334 // Unrolling has to be exactly linear in the trip count: the body planned
335 // once per iteration, no more and no less. Per species per iteration: 8
336 // dtheta/dphi + 4 analys transforms (Algorithm 3's cost, applied to the
337 // field and to each of its three Cartesian gradient components) plus 14
338 // generated kernels -- see test/modelChecks.ts's KERNELS_PER_ITERATION,
339 // which counts the kernels alone; this counts every op, transforms
340 // included.
341 const perIteration = ops[1] - ops[0];
342 const want = 52;
343 check(
344 'loop: unrolling is exactly linear in the trip count',
345 perIteration === want && ops[2] - ops[0] === 4 * perIteration,
346 `${perIteration} ops per iteration (expected ${want}), ` +
347 `${ops[2] - ops[0]} for 4 iterations`,
348 );
350 // On the sphere lap_g = lap_s exactly, so the correction should compute
351 // (numerically) close to zero regardless of niter -- not bit-identical
352 // (it is a real computation now, through 8+ chained fp32 transforms per
353 // iteration, not the literal `0 * Un` placeholder that used to make this
354 // exact), but close. The tolerance is set by that chain's fp32 roundoff,
355 // not by the scheme: a real geometry-correction bug would miss by orders
356 // of magnitude more than this.
357 let worst = 0;
358 for (let k = 1; k < states.length; k++) {
359 for (let i = 0; i < states[0].length; i++) {
360 worst = Math.max(worst, Math.abs(states[k][i] - states[0][i]));
361 }
362 }
363 check(
364 'loop: on the sphere, the correction stays near zero across niter',
365 worst < 2e-3,
366 `states differ by up to ${worst.toExponential(2)} after ${STEPS} steps at ${counts.join('/')} iterations`,
367 );
368 }
370 // ---- on a curved surface, the correction actually changes the answer ----
371 {
372 const model = mModelByKey('schnakenberg')!;
373 const params = defaultParams(model);
374 const peanut = mGeometryByKey('peanut')!;
375 const peanutParams = defaultGeometryParams(peanut);
376 // niter 0 vs 1 only -- deliberately not the 4/8 the sweep below already
377 // documents as outside the Richardson iteration's convergence radius on
378 // this geometry. The point here is just that the correction is not a
379 // no-op, which a much smaller, still-converging niter already shows.
380 const states: Float32Array[] = [];
381 for (const niter of [0, 1]) {
382 const session = await ModelSession.create({
383 device, model, params, lmax: SWEEP_LMAX,
384 geometry: peanut, geometryParams: peanutParams, niter,
385 });
386 session.seed(1);
387 session.step(STEPS);
388 states.push(await session.read('U'));
389 session.destroy();
390 }
391 let worst = 0;
392 for (let i = 0; i < states[0].length; i++) {
393 worst = Math.max(worst, Math.abs(states[1][i] - states[0][i]));
394 }
395 check(
396 'loop: on peanut, the correction measurably changes the answer',
397 worst > 1e-4 && states[1].every((v) => Number.isFinite(v)),
398 `states differ by ${worst.toExponential(2)} after ${STEPS} steps at niter 0 vs 1`,
399 );
400 }
402 // ---- two solvers, one operator ------------------------------------------
403 // solvers/bicgstab.m against solvers/richardson.m on the same implicit
404 // system: a Krylov iteration converges superlinearly where the stationary
405 // one converges linearly, so at equal niter it must land much closer to the
406 // converged answer. The comparison is a ratio against the same reference,
407 // which keeps it meaningful on SwiftShader's looser fp32 too.
408 {
409 const model = mModelByKey('schnakenberg')!;
410 const params = defaultParams(model);
411 const ellipsoid = mGeometryByKey('ellipsoid')!;
412 const run = async (source: string | undefined, niter: number): Promise<Float32Array> => {
413 const session = await ModelSession.create({
414 device, model, params, lmax: LMAX,
415 geometry: ellipsoid, geometryParams: defaultGeometryParams(ellipsoid),
416 niter, ...(source ? { source } : {}),
417 });
418 session.seed(1);
419 session.step(STEPS);
420 const U = await session.read('U');
421 session.destroy();
422 return U;
423 };
424 const ref = await run(undefined, 8); // richardson, effectively converged
425 const rich = await run(undefined, 2);
426 const bicg = await run(solverSchnak(model.source, 'bicgstab'), 2);
427 const gmres = await run(solverSchnak(model.source, 'gmres'), 2);
428 const relRich = relL2(rich, ref);
429 const relBicg = relL2(bicg, ref);
430 const relGmres = relL2(gmres, ref);
431 check(
432 'solvers: bicgstab(2) converges far past richardson(2) on the same operator',
433 bicg.every((v) => Number.isFinite(v)) && relBicg < relRich / 5 && relBicg < 1e-4,
434 `relL2 vs richardson(8): bicgstab ${relBicg.toExponential(2)}, ` +
435 `richardson ${relRich.toExponential(2)}`,
436 );
437 check(
438 'solvers: gmres(2) converges far past richardson(2) on the same operator',
439 gmres.every((v) => Number.isFinite(v)) && relGmres < relRich / 5 && relGmres < 1e-3,
440 `relL2 vs richardson(8): gmres ${relGmres.toExponential(2)}, ` +
441 `richardson ${relRich.toExponential(2)}`,
442 );
443 }
445 // ---- niter x geometry sweep: catch a "doesn't run" regression early -----
446 // This is what actually turned up the two real issues found while building
447 // the correction: peanut diverging at niter >= 4 with schnak-spots'
448 // shipped default dt (a genuine Richardson-convergence-radius limit, not a
449 // bug -- see docs/richardson-iteration.md), and a since-fixed compiler bug
450 // where a loop-body statement could silently reuse a *different*
451 // statement's compiled kernel (test/modelChecks.ts's pipeline-cache check
452 // guards that one directly). Every shipped geometry x every niter the
453 // app's <select> actually offers, so a regression anywhere in that grid is
454 // caught -- without asserting away the one combination already known to be
455 // outside the convergence radius.
456 //
457 // SWEEP_LMAX cannot be lowered to make this cheaper: at lmax 31 or 15 the
458 // peanut/2, /4 and /8 combinations below all stay finite, so a smaller grid
459 // would quietly turn KNOWN_DIVERGENT into three failures and stop testing
460 // the thing this sweep exists to pin down.
461 if (!runSweep) {
462 log(
463 ' sweep: skipped — run `npm run test:node` (desktop Dawn, ~3 s) or ' +
464 '`npm run test:gpu -- --sweep` for the niter x geometry sweep.',
465 );
466 } else {
467 const model = mModelByKey('schnakenberg')!;
468 const params = defaultParams(model);
469 const SWEEP_NITER = [0, 1, 2, 4, 8];
470 const KNOWN_DIVERGENT = new Set(['peanut/2', 'peanut/4', 'peanut/8']);
472 for (const geomSpec of mGeometries) {
473 for (const niter of SWEEP_NITER) {
474 const session = await ModelSession.create({
475 device, model, params, lmax: SWEEP_LMAX,
476 geometry: geomSpec, geometryParams: defaultGeometryParams(geomSpec),
477 niter,
478 });
479 session.seed(1);
480 session.step(STEPS);
481 const values = await session.read('u');
482 const finite = values.every((v) => Number.isFinite(v));
483 session.destroy();
485 const key = `${geomSpec.key}/${niter}`;
486 const expectDivergent = KNOWN_DIVERGENT.has(key);
487 check(
488 expectDivergent
489 ? `sweep: ${key} is known to diverge (outside the Richardson convergence radius)`
490 : `sweep: ${key} stays finite after ${STEPS} steps`,
491 expectDivergent ? !finite : finite,
492 expectDivergent
493 ? finite
494 ? 'now finite -- the convergence radius may have improved; update KNOWN_DIVERGENT'
495 : 'diverged as expected'
496 : finite
497 ? 'finite'
498 : 'NOT FINITE -- unexpected divergence, investigate before treating this as another known case',
499 );
500 }
501 }
503 // The other side of KNOWN_DIVERGENT: on the very combinations where the
504 // Richardson iteration leaves its convergence radius, the Krylov solvers
505 // — same operator, same preconditioner — keep converging in niter. This
506 // is what having the solver as its own .m is for.
507 {
508 const model = mModelByKey('schnakenberg')!;
509 const params = defaultParams(model);
510 const peanut = mGeometryByKey('peanut')!;
511 const run = async (solver: 'bicgstab' | 'gmres', niter: number): Promise<Float32Array> => {
512 const session = await ModelSession.create({
513 device, model, params, lmax: SWEEP_LMAX,
514 geometry: peanut, geometryParams: defaultGeometryParams(peanut),
515 niter, source: solverSchnak(model.source, solver),
516 });
517 session.seed(1);
518 session.step(STEPS);
519 const U = await session.read('U');
520 session.destroy();
521 return U;
522 };
523 const finiteAll = (U: Float32Array): boolean => U.every((v) => Number.isFinite(v));
525 const bicg = new Map<number, Float32Array>();
526 for (const niter of [1, 2, 4, 8]) bicg.set(niter, await run('bicgstab', niter));
527 const bref = bicg.get(8)!;
528 const brel = (n: number): number => relL2(bicg.get(n)!, bref);
529 check(
530 'sweep: bicgstab converges on the peanut combinations richardson cannot',
531 [...bicg.values()].every(finiteAll) && brel(4) < brel(2) && brel(2) < brel(1),
532 `relL2 vs bicgstab(8): niter 1 -> ${brel(1).toExponential(2)}, ` +
533 `2 -> ${brel(2).toExponential(2)}, 4 -> ${brel(4).toExponential(2)}`,
534 );
536 // gmres exercises the whole indexed-access machinery (the basis bank,
537 // the Hessenberg updates, the triangular inner loops) at the sweep's
538 // full lmax, on the operator's hardest shipped case.
539 const g1 = await run('gmres', 1);
540 const g4 = await run('gmres', 4);
541 const grel1 = relL2(g1, bref);
542 const grel4 = relL2(g4, bref);
543 check(
544 'sweep: gmres converges there too',
545 finiteAll(g1) && finiteAll(g4) && grel4 < grel1,
546 `relL2 vs bicgstab(8): niter 1 -> ${grel1.toExponential(2)}, ` +
547 `4 -> ${grel4.toExponential(2)}`,
548 );
549 }
550 }
552 // ---- a loop whose length is not known at compile time is refused --------
553 {
554 const model = mModelByKey('allencahn')!;
555 // `dt` is a tunable parameter, so it reaches the compiler with no value:
556 // the plan cannot know how many iterations to emit. The model's own loop
557 // lives in solvers/richardson.m now, so the bad loop is written out here.
558 const bad = `
559function [U, u] = init(noise)
560 U = analys(noise);
561 u = synth(U);
562end
564function [Un, u] = step(U, lam, eps2, dt, niter)
565 u = synth(U);
566 Bu = U + dt * analys(u - u.^3);
567 Un = Bu ./ (1 + (dt * eps2) * lam);
568 for k = 1:dt
569 Un = Un + 0 * Un;
570 end
571end
572`;
573 let message = '';
574 try {
575 const session = await ModelSession.create({
576 device, model, params: defaultParams(model), lmax: LMAX, source: bad, niter: 1,
577 });
578 session.destroy();
579 } catch (e) {
580 message = e instanceof ModelCompileError ? e.message : `wrong error type: ${e}`;
581 }
582 check(
583 'loop: a runtime loop bound is refused at compile time',
584 message.includes('known when the model is compiled'),
585 message ? `refused: ${message.slice(0, 72)}…` : 'compiled anyway',
586 );
587 }
589 // ---- swapping the surface leaves the simulation alone ------------------
590 {
591 const model = mModelByKey('schnakenberg')!;
592 const session = await ModelSession.create({
593 device, model, params: defaultParams(model), lmax: LMAX,
594 });
595 session.seed(1);
596 session.step(STEPS);
597 const before = await session.read('U');
599 const peanut = mGeometryByKey('peanut')!;
600 await session.setGeometry(peanut, defaultGeometryParams(peanut));
601 const after = await session.read('U');
603 let survived = before.length === after.length;
604 for (let i = 0; survived && i < before.length; i++) {
605 if (before[i] !== after[i]) survived = false;
606 }
607 const { lo, hi } = session.geometry.radiusRange();
608 check(
609 'geometry: swapping the surface mid-run does not disturb the state',
610 survived && session.geometryModel.key === 'peanut' && hi - lo > 0.1,
611 survived
612 ? `state identical, now on ${session.geometryModel.key} (radius ${lo.toFixed(3)}–${hi.toFixed(3)})`
613 : 'state changed',
614 );
615 session.destroy();
616 }
617}
619/** Index of the entry minimizing `score`, over the first `n` entries. */
620function argMin(
621 xs: Float64Array | Float32Array,
622 n: number,
623 score: (v: number) => number,
624): number {
625 let best = 0;
626 let bestScore = Infinity;
627 for (let i = 0; i < n; i++) {
628 const s = score(xs[i]);
629 if (s < bestScore) {
630 bestScore = s;
631 best = i;
632 }
633 }
634 return best;
635}