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