Compare several solver settings side by side, on one clock - concept-collection/turing-surface

concept-collection / turing-surface

Compare several solver settings side by side, on one clock

A convergence study in the page: pick a set of niter / lmax / dt values and run them all at once, one row per variant, so "does this setting matter?" is visible rather than argued. It is a mode, not a widening of the ordinary controls. With the Compare bar closed nothing about using the page changes; opening it and pressing Compare tears down the single session and hands the panels area to a CompareRun that owns one session per variant, and pressing it again puts the single run back. Three things the variants are forced to share, each of which the comparison would be meaningless without: - One initial condition. The host's seeded perturbation is one deviate per grid point, so two sessions at different lmax seeded from the same integer start from unrelated fields. src/compare/sharedStart.ts builds the field once, band-limited at the coarsest variant's lmax, and evaluates it on each variant's own grid -- exact, since every band contains that one. - One clock. dt varies only by a power-of-two divisor, and a frame advances each variant by frameSteps * dtDiv of its own steps. Every variant lands on the same model time at the end of every frame, so nothing is ever compared across a time offset. - One grid to look at. Each session's display plan is pointed at a common grid (the new ModelSession.setDisplayGrid), which is evaluation rather than resampling because the state is band-limited. The fields come back directly comparable point by point, one mesh serves every panel, and the relative L2 difference from the reference is an ordinary weighted sum. Columns share a colour range, or the panels would each be scaled to themselves and hide the difference. The range is set by whichever variant reaches least far from zero: a comparison between rows rather than a threshold on any of them, so a variant that leaves its convergence radius stops setting the scale from its first bad step instead of dragging the whole grid flat until some limit trips. A row that has gone non-finite says so in its label. Also: the default solve-iteration count is now 8, taken from DEFAULT_NITER so the page and `npm run bench` cannot start out disagreeing about it.

Jeremy Magland <jmagland@flatironinstitute.org> committed 14 days ago commit f467ffaaf82b parent ae74084 Browse files

10 changed files+1475−26

index.htmlmodified+81−2View file

@@ -171,6 +171,63 @@

171171 .editor-code { flex: none; height: 26em; }

172172 #compiled { border-left: 0; border-top: 1px solid var(--line); max-height: 12em; }

173173 }

174+ /* ---- compare mode ------------------------------------------------ */

175+ /* The bar's own layout: three chip rows stacked, then the reference

176+ picker and the button beside them. */

177+ .cmp-axes { display: flex; flex-direction: column; gap: 4px; }

178+ .cmp-axis { display: flex; align-items: center; gap: 8px; }

179+ .cmp-axis > span:first-child {

180+ color: var(--ink-2); font-size: 13px; width: 5.5em; text-align: right;

181+ }

182+ .chips { display: flex; flex-wrap: wrap; gap: 4px; }

183+ .chip {

184+ font: inherit; font-size: 12px; padding: 2px 8px;

185+ border: 1px solid var(--line); border-radius: 999px;

186+ background: var(--bg); color: var(--ink-2); cursor: pointer;

187+ }

188+ .chip:hover { border-color: var(--accent); }

189+ .chip[aria-pressed="true"] {

190+ border-color: var(--accent); color: var(--accent);

191+ background: color-mix(in srgb, var(--accent) 12%, transparent);

192+ font-weight: 600;

193+ }

194+ /* The panel area becomes a stack of labelled rows. Overrides the flex

195+ wrap the single-run view uses. */

196+ #panels.compare { flex-direction: column; gap: 8px; }

197+ .cmp-row { display: flex; align-items: stretch; gap: 8px; }

198+ .cmp-rowlabel {

199+ flex: none; width: 13em; padding: 6px 8px;

200+ border-left: 4px solid var(--c, var(--line));

201+ font-size: 12px; color: var(--ink-2);

202+ display: flex; flex-direction: column; justify-content: center; gap: 3px;

203+ }

204+ .cmp-rowname { color: var(--ink); font-weight: 600; }

205+ .cmp-rowstat { font-variant-numeric: tabular-nums; line-height: 1.35; }

206+ .cmp-diverged { color: var(--warn-edge); }

207+ /* The header row's label cell is a spacer, not a variant — no swatch. */

208+ .cmp-head .cmp-rowlabel { border-left-color: transparent; padding: 0 8px; }

209+ .cmp-cols { flex: 1; display: flex; gap: 8px; min-width: 0; }

210+ .cmp-box {

211+ flex: 1 1 0; min-width: 0;

212+ border: 1px solid var(--line); border-radius: 8px; overflow: hidden;

213+ max-height: 42vh;

214+ }

215+ .cmp-head { align-items: flex-end; }

216+ .cmp-colhead {

217+ flex: 1 1 0; min-width: 0;

218+ display: flex; align-items: center; gap: 8px;

219+ font-size: 12px; color: var(--ink-2);

220+ }

221+ .cmp-rangebar {

222+ flex: 1 1 auto; min-width: 0; height: 8px;

223+ border: 1px solid var(--line); border-radius: 2px;

224+ }

225+ .cmp-rangelab { font-variant-numeric: tabular-nums; white-space: nowrap; }

226+ @media (max-width: 860px) {

227+ .cmp-row { flex-direction: column; }

228+ .cmp-rowlabel { width: auto; flex-direction: row; gap: 10px; }

229+ .cmp-head { display: none; }

230+ }

174231 .tok-com { color: var(--tok-com); }

175232 .tok-str { color: var(--tok-str); }

176233 .tok-num { color: var(--tok-num); }

@@ -207,10 +264,10 @@

207264 <label title="Iterations of the implicit diffusion solve. Changing it recompiles.">solve iters

208265 <select id="niter">

209266 <option value="0">0</option>

210- <option value="1" selected>1</option>

267+ <option value="1">1</option>

211268 <option value="2">2</option>

212269 <option value="4">4</option>

213- <option value="8">8</option>

270+ <option value="8" selected>8</option>

214271 <option value="16">16</option>

215272 <option value="32">32</option>

216273 <option value="64">64</option>

@@ -242,6 +299,28 @@

242299 <button id="reseed">Re-seed</button>

243300 <button id="resetview">Reset view</button>

244301 <button id="movietoggle" title="Export the run as an MP4 movie">Export movie</button>

302+ <button id="comparetoggle" title="Run several solver settings side by side on one clock">Compare</button>

303+ </div>

304+ <div class="controls" id="comparebar" hidden>

305+ <div class="cmp-axes">

306+ <div class="cmp-axis">

307+ <span title="Iterations of the implicit diffusion solve">solve iters</span>

308+ <span id="cmp-niter" class="chips"></span>

309+ </div>

310+ <div class="cmp-axis">

311+ <span>lmax</span>

312+ <span id="cmp-lmax" class="chips"></span>

313+ </div>

314+ <div class="cmp-axis">

315+ <span title="Timestep, as a divisor of the model's dt. Divisors keep every variant on the same clock exactly.">dt</span>

316+ <span id="cmp-dt" class="chips"></span>

317+ </div>

318+ </div>

319+ <label title="The run everything else is measured against">reference

320+ <select id="cmp-ref"></select>

321+ </label>

322+ <button id="cmp-start" class="primary">Compare</button>

323+ <span id="cmp-count" class="stats"></span>

245324 </div>

246325 <div class="controls" id="moviebar" hidden>

247326 <label title="Simulation-time units per second of video">movie speed

scripts/test-node.tsmodified+2−0View file

@@ -15,6 +15,7 @@ import { analyticChecks } from '../test/analyticChecks.ts';

1515 import { modelChecks } from '../test/modelChecks.ts';

1616 import { geometryChecks } from '../test/geometryChecks.ts';

1717 import { fluxChecks } from '../test/fluxChecks.ts';

18+import { compareChecks } from '../test/compareChecks.ts';

1819

1920 let failures = 0;

2021 const check = (name: string, ok: boolean, detail: string): void => {

@@ -53,6 +54,7 @@ await analyticChecks(device, check, log);

5354 await modelChecks(device, check, log);

5455 await geometryChecks(device, check, log);

5556 await fluxChecks(device, check, log);

57+await compareChecks(device, check, log);

5658

5759 console.log(failures === 0 ? '\nAll tests passed.' : `\n${failures} failed.`);

5860 process.exit(failures === 0 ? 0 : 1);

src/bench/runSpec.tsmodified+3−1View file

@@ -46,7 +46,9 @@ export interface RunSpec {

4646 niter: number;

4747 }

4848

49-export const DEFAULT_NITER = 1;

49+/** Iterations of the implicit solve, everywhere that does not say otherwise:

50+ * the app's `solve iters` control, `npm run bench`, and the soak. */

51+export const DEFAULT_NITER = 8;

5052

5153 /** Geometry + starting parameters of a geometry key. */

5254 export function resolveGeometry(key: string): { geometry: MGeometry; params: Params } {

src/compare/compareRun.tsadded+721−0View file

@@ -0,0 +1,721 @@

1+/**

2+ * Several solver settings, one problem, one clock.

3+ *

4+ * A convergence study of the knobs that decide how well the implicit solve is

5+ * resolved — `niter`, `lmax`, `dt` — run side by side so the answer to "does it

6+ * matter?" is visible rather than argued. Every variant is its own

7+ * `ModelSession` (both `niter` and `lmax` are structural: they change the

8+ * compiled step and the grid), and what makes the set a comparison rather than

9+ * a collection is three things they are forced to share:

10+ *

11+ * - **One initial condition.** Band-limited at the coarsest variant's lmax and

12+ * evaluated on each variant's own grid, so every session starts from the same

13+ * *function* rather than from the same random seed — see sharedStart.ts for

14+ * why the seed alone is not enough.

15+ *

16+ * - **One clock.** Variants differ in dt only by an integer power-of-two

17+ * divisor, and a frame advances each of them by `frameSteps * dtDiv` steps.

18+ * Every variant therefore lands on exactly the same model time at the end of

19+ * every frame, having taken a different number of steps to get there. Nothing

20+ * is ever compared across a time offset.

21+ *

22+ * - **One grid to look at.** Each session's *display* plan is pointed at a

23+ * common grid (ModelSession.setDisplayGrid), which is exact evaluation rather

24+ * than resampling because the state is band-limited. So the fields come back

25+ * directly comparable point by point, one mesh topology serves every panel,

26+ * and the difference norm is an ordinary weighted sum.

27+ *

28+ * What is *not* shared is the surface: each variant carries the geometry

29+ * band-limited at its own lmax, and renders the surface it actually solves on.

30+ */

31+import { ModelSession } from '../mgpu/session.ts';

32+import type { MModel, Params } from '../mgpu/registry.ts';

33+import type { MGeometry } from '../geom/registry.ts';

34+import {

35+ buildTopology,

36+ fillFieldValues,

37+ fillPositions,

38+ fillColors,

39+ type SphereMeshTopology,

40+} from '../render/sphereMesh.ts';

41+import { SphereScene } from '../render/SphereScene.ts';

42+import { colormaps } from '../render/colormaps.ts';

43+import { fmtValue } from '../render/colorbar.ts';

44+import { sharedNoise } from './sharedStart.ts';

45+import { variantLabel, VARIANT_COLORS, type Variant } from './variants.ts';

46+

47+/**

48+ * Latitudes of the shared display grid. 256 is the same target the single-run

49+ * view uses for 'auto' oversampling, and for the same reason — beyond it a

50+ * finer mesh costs vertices without showing anything.

51+ *

52+ * Here it is a ceiling as well as a target, in two directions. At lmax 255 the

53+ * solver grid is finer than this, so the panels sample the (exact) state more

54+ * coarsely than the solver carries it; and past a handful of panels the mesh is

55+ * paid for once per panel, in vertices, normals and a WebGL context each, so it

56+ * halves. Both are display choices, both are reported in the status line, and

57+ * neither touches the difference norm's meaning: that is computed on this same

58+ * grid for every variant, so it stays a consistent comparison whatever the grid.

59+ */

60+const RENDER_NLAT = 256;

61+const RENDER_NLAT_CROWDED = 128;

62+const CROWDED_PANELS = 6;

63+

64+/** See main.ts's DISPATCH_BUDGET — the same watchdog argument, per variant. */

65+const DISPATCH_BUDGET = 1000;

66+const STEPS_PER_FRAME_BASE = 4;

67+

68+export interface CompareOptions {

69+ device: GPUDevice;

70+ model: MModel;

71+ /** The model's parameters, with `dt` read as the *base* timestep that each

72+ * variant's dtDiv divides. */

73+ params: Params;

74+ source: string;

75+ geometry: MGeometry;

76+ geometryParams: Params;

77+ geometrySource: string;

78+ variants: Variant[];

79+ /** Index into `variants` of the run everything else is measured against. */

80+ reference: number;

81+ seed: number;

82+ morph: number;

83+ colormapName: () => string;

84+ /** Where the variant grid goes (the app's #panels). */

85+ container: HTMLElement;

86+ /** Progress and, afterwards, the standing description of the study. */

87+ onStatus: (html: string) => void;

88+}

89+

90+interface Row {

91+ variant: Variant;

92+ session: ModelSession;

93+ color: string;

94+ /** Surface coordinates on the shared render grid — this variant's own. */

95+ coords: Float32Array;

96+ posBuf: Float32Array;

97+ scenes: SphereScene[];

98+ valueBufs: Float32Array[];

99+ colorBufs: Float32Array[];

100+ /** Fields read this frame, one per species, on the shared grid. */

101+ fields: Float32Array[];

102+ /** Relative difference from the reference, one per species. */

103+ err: number[];

104+ /** False once any species has left the floating-point numbers — the shape a

105+ * variant outside the convergence radius eventually fails in. Such a row is

106+ * never used to scale a column, and its label says so. */

107+ healthy: boolean;

108+ statEl: HTMLElement;

109+}

110+

111+export class CompareRun {

112+ #opts: CompareOptions;

113+ #rows: Row[] = [];

114+ #topo: SphereMeshTopology;

115+ /** Quadrature weight per grid point of the shared grid, for the L2 norm. */

116+ #weights: Float64Array;

117+ #rangeBars: { fill: (lo: number, hi: number) => void }[] = [];

118+ /** Smoothed color range per species, shared by every variant so the panels

119+ * in a column are directly comparable by eye and not just by number. */

120+ #ranges: { lo: number; hi: number }[] = [];

121+ #resizeObs: ResizeObserver | null = null;

122+

123+ #running = false;

124+ #pumping = false;

125+ #disposed = false;

126+ #morph: number;

127+ /** Base steps per frame; variant i takes this times its dtDiv. */

128+ #frameSteps = STEPS_PER_FRAME_BASE;

129+ /** Model time all variants are at — one number, by construction. */

130+ #t = 0;

131+ #frameMs = 0;

132+ #note: string;

133+

134+ private constructor(init: {

135+ opts: CompareOptions;

136+ rows: Row[];

137+ topo: SphereMeshTopology;

138+ weights: Float64Array;

139+ rangeBars: { fill: (lo: number, hi: number) => void }[];

140+ frameSteps: number;

141+ note: string;

142+ }) {

143+ this.#opts = init.opts;

144+ this.#rows = init.rows;

145+ this.#topo = init.topo;

146+ this.#weights = init.weights;

147+ this.#rangeBars = init.rangeBars;

148+ this.#frameSteps = init.frameSteps;

149+ this.#note = init.note;

150+ this.#morph = init.opts.morph;

151+ this.#ranges = init.opts.model.species.map(() => ({ lo: NaN, hi: NaN }));

152+ }

153+

154+ get variants(): Variant[] {

155+ return this.#rows.map((r) => r.variant);

156+ }

157+

158+ /** The variant everything else is measured against — the one whose numbers

159+ * stand on their own, so the one the app quotes when it has to quote one. */

160+ get referenceSession(): ModelSession | null {

161+ return this.#rows[this.#opts.reference]?.session ?? null;

162+ }

163+

164+ get referenceIndex(): number {

165+ return this.#opts.reference;

166+ }

167+

168+ /** The base timestep a variant's dtDiv divides. */

169+ static baseDt(params: Params): number {

170+ return params.dt ?? 0;

171+ }

172+

173+ static async create(opts: CompareOptions): Promise<CompareRun> {

174+ const { device, model, variants } = opts;

175+ const baseDt = CompareRun.baseDt(opts.params);

176+ const showDt = variants.some((v) => v.dtDiv !== variants[0].dtDiv);

177+ const sessions: ModelSession[] = [];

178+ // Scenes own a WebGL context and an animation frame each, so a failure

179+ // after the grid is up has to take them down explicitly — removing their

180+ // canvases from the DOM would leave both running.

181+ let built: Row[] = [];

182+

183+ try {

184+ for (let i = 0; i < variants.length; i++) {

185+ const v = variants[i];

186+ opts.onStatus(

187+ `compiling ${i + 1}/${variants.length} — ${variantLabel(v, showDt)} ` +

188+ `(a solve iteration is ~15 kernels per species, and there is no ` +

189+ `pipeline cache across sessions)`,

190+ );

191+ // Yield, so the status actually paints before the compile blocks.

192+ await new Promise<number>(requestAnimationFrame);

193+ sessions.push(

194+ await ModelSession.create({

195+ device,

196+ model,

197+ params: { ...opts.params, dt: baseDt / v.dtDiv },

198+ lmax: v.lmax,

199+ source: opts.source,

200+ geometry: opts.geometry,

201+ geometryParams: opts.geometryParams,

202+ geometrySource: opts.geometrySource,

203+ niter: v.niter,

204+ }),

205+ );

206+ }

207+

208+ // ---- the shared display grid ----------------------------------------

209+ const maxLmax = Math.max(...variants.map((v) => v.lmax));

210+ const panels = variants.length * model.species.length;

211+ const target = panels > CROWDED_PANELS ? RENDER_NLAT_CROWDED : RENDER_NLAT;

212+ // Never below what the finest band needs to be representable at all

213+ // (ShtPlan requires nlat > lmax), whatever the panel count says.

214+ const nlat = Math.max(target, 2 * Math.ceil((maxLmax + 2) / 2));

215+ let nphi = 1;

216+ while (nphi < Math.max(2 * nlat, 2 * maxLmax + 1)) nphi *= 2;

217+ for (const s of sessions) await s.setDisplayGrid(nlat, nphi);

218+

219+ // ---- one initial condition, on every grid ---------------------------

220+ opts.onStatus('seeding all variants from one band-limited perturbation…');

221+ const noise = await sharedNoise(sessions, model.seedAmp, opts.seed);

222+ sessions.forEach((s, i) => s.seedWith(noise[i]));

223+

224+ // ---- the mesh, shared; the surface, per variant ---------------------

225+ const view = sessions[0].viewSht;

226+ const phi = new Float64Array(nphi);

227+ for (let j = 0; j < nphi; j++) phi[j] = (2 * Math.PI * j) / nphi;

228+ const topo = buildTopology(view.cosTheta, phi);

229+ // Gauss weights carry the sin(theta) of the area element; the constant

230+ // 2*pi/nphi is common to every point and cancels in the relative norm.

231+ const weights = new Float64Array(nlat * nphi);

232+ for (let i = 0; i < nlat; i++) {

233+ for (let j = 0; j < nphi; j++) weights[i * nphi + j] = view.gaussWeights[i];

234+ }

235+

236+ // ---- how many steps a frame may submit ------------------------------

237+ // Per variant: its own unrolled step size times its dtDiv, since a ÷K

238+ // variant takes K times as many steps to reach the same time.

239+ let frameSteps = STEPS_PER_FRAME_BASE;

240+ const ops: number[] = [];

241+ for (let i = 0; i < sessions.length; i++) {

242+ const n = Math.max(1, sessions[i].describe().step.length);

243+ ops.push(n);

244+ frameSteps = Math.min(

245+ frameSteps,

246+ Math.max(1, Math.floor(DISPATCH_BUDGET / (n * variants[i].dtDiv))),

247+ );

248+ }

249+ frameSteps = Math.max(1, frameSteps);

250+

251+ // ---- the grid of panels ---------------------------------------------

252+ const { rows, rangeBars } = await buildGrid(opts, sessions, topo, showDt);

253+ built = rows;

254+

255+ const solverGrid = sessions.map((s) => `${s.cfg.nlat}×${s.cfg.nphi}`);

256+ const note =

257+ `${variants.length} variants · display grid ${nlat}×${nphi}` +

258+ (sessions.some((s) => s.cfg.nlat > nlat)

259+ ? ` (below the finest solver grid ${solverGrid[solverGrid.length - 1]} — display only)`

260+ : '') +

261+ ` · ${frameSteps} base step${frameSteps === 1 ? '' : 's'}/frame` +

262+ ` · ops/step ${ops.join(', ')}`;

263+

264+ const run = new CompareRun({

265+ opts, rows, topo, weights, rangeBars, frameSteps, note,

266+ });

267+ await run.draw();

268+ run.#observeResize();

269+ run.#status();

270+ return run;

271+ } catch (e) {

272+ for (const r of built) for (const s of r.scenes) s.dispose();

273+ for (const s of sessions) s.destroy();

274+ opts.container.replaceChildren();

275+ opts.container.classList.remove('compare');

276+ throw e;

277+ }

278+ }

279+

280+ // ------------------------------------------------------------------ state

281+ setRunning(next: boolean): void {

282+ this.#running = next;

283+ if (next) void this.#pump();

284+ }

285+

286+ get running(): boolean {

287+ return this.#running;

288+ }

289+

290+ /** Re-seed every variant from one new shared perturbation. */

291+ async reseed(seed: number): Promise<void> {

292+ const wasRunning = this.#running;

293+ this.#running = false;

294+ while (this.#pumping) await nextFrame();

295+ if (this.#disposed) return;

296+ const noise = await sharedNoise(

297+ this.#rows.map((r) => r.session),

298+ this.#opts.model.seedAmp,

299+ seed,

300+ );

301+ if (this.#disposed) return;

302+ this.#rows.forEach((r, i) => r.session.seedWith(noise[i]));

303+ this.#t = 0;

304+ for (const r of this.#ranges) {

305+ r.lo = NaN;

306+ r.hi = NaN;

307+ }

308+ await this.draw();

309+ if (!this.#disposed && wasRunning) this.setRunning(true);

310+ }

311+

312+ /** Model parameters changed. Each variant keeps its own dt. */

313+ setParams(params: Params): void {

314+ this.#opts.params = params;

315+ const baseDt = CompareRun.baseDt(params);

316+ for (const r of this.#rows) {

317+ r.session.setParams({ ...params, dt: baseDt / r.variant.dtDiv });

318+ }

319+ }

320+

321+ setMorph(morph: number): void {

322+ this.#morph = morph;

323+ for (const r of this.#rows) {

324+ fillPositions(r.posBuf, r.coords, this.#topo, morph);

325+ for (const s of r.scenes) s.updatePositions(r.posBuf);

326+ }

327+ }

328+

329+ resetView(): void {

330+ for (const r of this.#rows) for (const s of r.scenes) s.resetCamera();

331+ }

332+

333+ dispose(): void {

334+ this.#disposed = true;

335+ this.#running = false;

336+ this.#resizeObs?.disconnect();

337+ this.#resizeObs = null;

338+ for (const r of this.#rows) {

339+ for (const s of r.scenes) s.dispose();

340+ r.session.destroy();

341+ }

342+ this.#rows = [];

343+ this.#opts.container.replaceChildren();

344+ this.#opts.container.classList.remove('compare');

345+ }

346+

347+ // ----------------------------------------------------------------- drawing

348+ /**

349+ * One frame's readback: every variant's every species, on the shared grid.

350+ * Read first, then color — the range is shared down a column, so no panel can

351+ * be filled until the column's range is known.

352+ */

353+ async draw(): Promise<void> {

354+ if (this.#disposed) return;

355+ const species = this.#opts.model.species;

356+ // Sessions are independent, so their readbacks can be in flight together;

357+ // within one session they must not be (they share its staging buffers).

358+ await Promise.all(

359+ this.#rows.map(async (r) => {

360+ for (let k = 0; k < species.length; k++) {

361+ r.fields[k] = await r.session.readSpecies(k);

362+ }

363+ }),

364+ );

365+ if (this.#disposed) return;

366+

367+ const cmap = colormaps[this.#opts.colormapName()] ?? colormaps.viridis;

368+

369+ /**

370+ * What scales a column is the whole question, and it has three wrong

371+ * answers.

372+ *

373+ * Per panel is wrong: a range each rescales every variant to itself and

374+ * hides exactly the difference the grid exists to show. The union over

375+ * variants is wrong for the opposite reason: a variant outside the

376+ * iteration's convergence radius runs away to 1e20 and then to NaN, and a

377+ * union range rescales the *whole column* to it, flattening every panel to

378+ * one colour — which reads as "they all blew up" when only one did.

379+ *

380+ * The reference alone is wrong too, less obviously, and it is the case that

381+ * actually bites: outside the convergence radius *more* Richardson

382+ * iterations diverge *faster*, so the row that goes first is usually the

383+ * highest-niter one — which is the reference.

384+ *

385+ * So the column is scaled by whichever variant **reaches least far from

386+ * zero** — the least-blown-up one. That is a comparison between the rows,

387+ * not a threshold on any of them, and the distinction is the whole point:

388+ * any "is this value too big?" test has a window in which a diverging field

389+ * is still under the limit, and for as long as that window lasts it drags

390+ * the scale and flattens the grid, until it finally trips and everything

391+ * springs back. A comparison has no such window — a run-away only has to be

392+ * *larger* than a healthy row to stop setting the scale, which it is from

393+ * its first bad step, and it stays larger no matter how many other rows go

394+ * with it. One healthy variant is enough to keep the grid readable.

395+ *

396+ * The cost is a slight bias: among healthy variants the scale comes from

397+ * the one with the smallest peak, so the others clip by however much they

398+ * exceed it. They are approximations of the same solution, so that is a

399+ * fraction of a percent, and the alternative is a display that a single

400+ * divergence can take away.

401+ */

402+ const bounds = this.#rows.map((r) => species.map((_, k) => finiteRange(r.fields[k])));

403+ this.#rows.forEach((r, i) => {

404+ // A row with any non-finite value is out of the running entirely: its

405+ // finite entries are whatever survived, and no rank over them means much.

406+ r.healthy = species.every((_, k) => allFinite(r.fields[k]) && bounds[i][k] !== null);

407+ });

408+

409+ for (let k = 0; k < species.length; k++) {

410+ const anchor = leastPeak(this.#rows.map((r, i) => (r.healthy ? bounds[i][k] : null)));

411+ const range = this.#ranges[k];

412+ if (anchor) {

413+ if (!Number.isFinite(range.lo)) {

414+ range.lo = anchor.lo;

415+ range.hi = anchor.hi;

416+ } else {

417+ // Smooth in both directions so the shading evolves gently as the

418+ // pattern grows, as the single-run view does.

419+ const a = 0.15;

420+ range.lo += a * (anchor.lo - range.lo);

421+ range.hi += a * (anchor.hi - range.hi);

422+ }

423+ }

424+ // With every row gone, the last good range is kept rather than replaced

425+ // by nothing: the panels freeze at a readable scale and the row labels

426+ // say what happened, instead of the grid going blank.

427+ if (!Number.isFinite(range.lo) || !Number.isFinite(range.hi)) continue;

428+ if (range.hi - range.lo < 1e-9) {

429+ const mid = (range.hi + range.lo) / 2;

430+ range.lo = mid - 5e-10;

431+ range.hi = mid + 5e-10;

432+ }

433+ this.#rangeBars[k]?.fill(range.lo, range.hi);

434+ for (const r of this.#rows) {

435+ fillFieldValues(r.valueBufs[k], r.fields[k], this.#topo);

436+ fillColors(r.colorBufs[k], r.valueBufs[k], range.lo, range.hi, cmap);

437+ r.scenes[k]?.updateColors(r.colorBufs[k]);

438+ }

439+ }

440+

441+ this.#measureDifference();

442+ this.#updateRowStats();

443+ }

444+

445+ /**

446+ * Relative L2 difference from the reference, per species, on the shared

447+ * grid. Weighted by the Gauss weights, so it is the norm on the parameter

448+ * sphere — not on the embedded surface, which would weight by the area

449+ * element. That makes it a consistent diagnostic across variants rather than

450+ * a physical quantity, which is all it is used for.

451+ */

452+ #measureDifference(): void {

453+ const ref = this.#rows[this.#opts.reference];

454+ if (!ref) return;

455+ const species = this.#opts.model.species;

456+ for (const r of this.#rows) {

457+ for (let k = 0; k < species.length; k++) {

458+ if (r === ref) {

459+ r.err[k] = 0;

460+ continue;

461+ }

462+ const a = r.fields[k];

463+ const b = ref.fields[k];

464+ if (!a || !b || a.length !== b.length) {

465+ r.err[k] = NaN;

466+ continue;

467+ }

468+ let num = 0;

469+ let den = 0;

470+ for (let i = 0; i < a.length; i++) {

471+ const w = this.#weights[i];

472+ const d = a[i] - b[i];

473+ num += w * d * d;

474+ den += w * b[i] * b[i];

475+ }

476+ r.err[k] = den > 0 ? Math.sqrt(num / den) : NaN;

477+ }

478+ }

479+ }

480+

481+ /**

482+ * Each row's standing line: how many of its own steps it took to reach the

483+ * common time, and how far it is from the reference right now, per species.

484+ * Per species rather than a single worst-case number because the two are

485+ * genuinely different questions on a two-species model — the slow species is

486+ * usually the one that has converged and the fast one the one that has not.

487+ */

488+ #updateRowStats(): void {

489+ const species = this.#opts.model.species;

490+ const ref = this.#rows[this.#opts.reference];

491+ for (const r of this.#rows) {

492+ const per = species

493+ .map((s, k) => `${s} ${Number.isFinite(r.err[k]) ? r.err[k].toExponential(2) : '—'}`)

494+ .join('<br>');

495+ // Divergence is said, not implied. Scaled to a healthy row, a blown-up

496+ // variant is a flat saturated panel, which on its own is easy to misread

497+ // as a converged uniform state.

498+ const body = !r.healthy

499+ ? '<b class="cmp-diverged">diverged</b>'

500+ : r === ref

501+ ? '<b>reference</b>'

502+ : `Δ ${per}`;

503+ r.statEl.innerHTML = `${r.session.steps.toLocaleString()} steps<br>${body}`;

504+ }

505+ }

506+

507+ #status(): void {

508+ this.#opts.onStatus(

509+ `<b>t = ${this.#t.toFixed(2)}</b> (same for every variant) · ` +

510+ (this.#frameMs > 0 ? `${this.#frameMs.toFixed(1)} ms/frame · ` : '') +

511+ this.#note,

512+ );

513+ }

514+

515+ #observeResize(): void {

516+ const scenes = this.#rows.flatMap((r) => r.scenes);

517+ this.#resizeObs = new ResizeObserver(() => {

518+ for (const s of scenes) {

519+ const box = s.canvas.parentElement;

520+ if (box) s.resize(box.clientWidth, box.clientHeight);

521+ }

522+ });

523+ for (const s of scenes) {

524+ const box = s.canvas.parentElement;

525+ if (box) this.#resizeObs.observe(box);

526+ }

527+ }

528+

529+ // -------------------------------------------------------------- the clock

530+ /**

531+ * One frame advances every variant by the *same model time*: `frameSteps`

532+ * base steps, which a ÷K variant covers in K times as many of its own. That

533+ * is the whole reason dt varies by an integer divisor — the alternative is

534+ * rounding each variant to the nearest step and comparing fields that are a

535+ * fraction of a timestep apart, which would show up as a difference and be

536+ * indistinguishable from a real one.

537+ */

538+ async #pump(): Promise<void> {

539+ if (this.#pumping) return;

540+ this.#pumping = true;

541+ try {

542+ while (this.#running && !this.#disposed) {

543+ const t0 = performance.now();

544+ for (const r of this.#rows) r.session.step(this.#frameSteps * r.variant.dtDiv);

545+ this.#t += this.#frameSteps * CompareRun.baseDt(this.#opts.params);

546+ await this.draw();

547+ if (this.#disposed) break;

548+ const dt = performance.now() - t0;

549+ this.#frameMs = this.#frameMs === 0 ? dt : this.#frameMs + 0.05 * (dt - this.#frameMs);

550+ this.#status();

551+ await nextFrame();

552+ }

553+ if (!this.#disposed) {

554+ await this.draw();

555+ this.#status();

556+ }

557+ } finally {

558+ this.#pumping = false;

559+ }

560+ }

561+}

562+

563+const nextFrame = (): Promise<number> => new Promise(requestAnimationFrame);

564+

565+/** Whether every entry is an ordinary number — false once a variant has left

566+ * its convergence radius and saturated to infinity or NaN. */

567+function allFinite(f: Float32Array | undefined): boolean {

568+ if (!f) return false;

569+ for (let i = 0; i < f.length; i++) if (!Number.isFinite(f[i])) return false;

570+ return true;

571+}

572+

573+type Bounds = { lo: number; hi: number };

574+

575+/** How far a field reaches from zero — the one number the rows are ranked by

576+ * when deciding which of them sets a column's scale. */

577+const peak = (b: Bounds): number => Math.max(Math.abs(b.lo), Math.abs(b.hi));

578+

579+/** Whichever of the given bounds reaches least far from zero; null if none. */

580+function leastPeak(all: (Bounds | null)[]): Bounds | null {

581+ let best: Bounds | null = null;

582+ for (const b of all) {

583+ if (b !== null && (best === null || peak(b) < peak(best))) best = b;

584+ }

585+ return best;

586+}

587+

588+/** Min and max over the finite entries only; null when there are none. */

589+function finiteRange(f: Float32Array | undefined): { lo: number; hi: number } | null {

590+ if (!f) return null;

591+ let lo = Infinity;

592+ let hi = -Infinity;

593+ for (let i = 0; i < f.length; i++) {

594+ const v = f[i];

595+ if (!Number.isFinite(v)) continue;

596+ if (v < lo) lo = v;

597+ if (v > hi) hi = v;

598+ }

599+ return lo <= hi ? { lo, hi } : null;

600+}

601+

602+/**

603+ * The DOM: a header row naming each species and carrying that column's shared

604+ * color range, then one row per variant. The colorbar is per *column* rather

605+ * than per panel because the range is shared — a bar on every panel would be

606+ * the same bar repeated, and would suggest each panel had its own scaling,

607+ * which is exactly the thing that would make the comparison a lie.

608+ */

609+async function buildGrid(

610+ opts: CompareOptions,

611+ sessions: ModelSession[],

612+ topo: SphereMeshTopology,

613+ showDt: boolean,

614+): Promise<{ rows: Row[]; rangeBars: { fill: (lo: number, hi: number) => void }[] }> {

615+ const { container, model } = opts;

616+ container.replaceChildren();

617+ container.classList.add('compare');

618+

619+ const head = document.createElement('div');

620+ head.className = 'cmp-row cmp-head';

621+ const headSpacer = document.createElement('div');

622+ headSpacer.className = 'cmp-rowlabel';

623+ const headCols = document.createElement('div');

624+ headCols.className = 'cmp-cols';

625+ head.append(headSpacer, headCols);

626+ container.append(head);

627+

628+ const rangeBars = model.species.map((name) => {

629+ const col = document.createElement('div');

630+ col.className = 'cmp-colhead';

631+ const tag = document.createElement('b');

632+ tag.textContent = name;

633+ const canvas = document.createElement('canvas');

634+ canvas.width = 160;

635+ canvas.height = 8;

636+ canvas.className = 'cmp-rangebar';

637+ const lab = document.createElement('span');

638+ lab.className = 'cmp-rangelab';

639+ col.append(tag, canvas, lab);

640+ headCols.append(col);

641+ let painted = false;

642+ return {

643+ fill: (lo: number, hi: number): void => {

644+ const ctx = canvas.getContext('2d');

645+ if (ctx && !painted) {

646+ painted = true;

647+ const cmap = colormaps[opts.colormapName()] ?? colormaps.viridis;

648+ for (let x = 0; x < canvas.width; x++) {

649+ const [r, g, b] = cmap(x / (canvas.width - 1));

650+ ctx.fillStyle = `rgb(${r},${g},${b})`;

651+ ctx.fillRect(x, 0, 1, canvas.height);

652+ }

653+ }

654+ lab.textContent = `${fmtValue(lo)} … ${fmtValue(hi)}`;

655+ },

656+ };

657+ });

658+

659+ const sphereBg = getComputedStyle(document.documentElement)

660+ .getPropertyValue('--sphere-bg')

661+ .trim();

662+

663+ const rows: Row[] = [];

664+ for (let i = 0; i < sessions.length; i++) {

665+ const session = sessions[i];

666+ const variant = opts.variants[i];

667+ const color = VARIANT_COLORS[i % VARIANT_COLORS.length];

668+

669+ const coords = await session.renderPositions();

670+ const posBuf = new Float32Array(topo.numVertices * 3);

671+ fillPositions(posBuf, coords, topo, opts.morph);

672+

673+ const rowEl = document.createElement('div');

674+ rowEl.className = 'cmp-row';

675+ const labelEl = document.createElement('div');

676+ labelEl.className = 'cmp-rowlabel';

677+ labelEl.style.setProperty('--c', color);

678+ const nameEl = document.createElement('div');

679+ nameEl.className = 'cmp-rowname';

680+ nameEl.textContent = variantLabel(variant, showDt);

681+ const statEl = document.createElement('div');

682+ statEl.className = 'cmp-rowstat';

683+ labelEl.append(nameEl, statEl);

684+ const colsEl = document.createElement('div');

685+ colsEl.className = 'cmp-cols';

686+ rowEl.append(labelEl, colsEl);

687+ container.append(rowEl);

688+

689+ const scenes: SphereScene[] = [];

690+ const valueBufs: Float32Array[] = [];

691+ const colorBufs: Float32Array[] = [];

692+ for (let k = 0; k < model.species.length; k++) {

693+ const box = document.createElement('div');

694+ box.className = 'sphere-box cmp-box';

695+ colsEl.append(box);

696+ const scene = new SphereScene(

697+ box,

698+ topo.numVertices,

699+ topo.indices,

700+ Float32Array.from(posBuf),

701+ sphereBg || undefined,

702+ );

703+ scene.fitCamera();

704+ scenes.push(scene);

705+ valueBufs.push(new Float32Array(topo.numVertices));

706+ colorBufs.push(new Float32Array(topo.numVertices * 3));

707+ }

708+

709+ rows.push({

710+ variant, session, color, coords, posBuf, scenes, valueBufs, colorBufs,

711+ fields: [], err: model.species.map(() => 0), healthy: true, statEl,

712+ });

713+ }

714+

715+ // Every panel shares one camera: the study is about the fields, and looking

716+ // at two of them from different angles is not comparing them.

717+ const all = rows.flatMap((r) => r.scenes);

718+ for (let i = 1; i < all.length; i++) all[0].syncCamerasWith(all[i]);

719+

720+ return { rows, rangeBars };

721+}

src/compare/sharedStart.tsadded+76−0View file

@@ -0,0 +1,76 @@

1+/**

2+ * One initial condition, on every variant's grid.

3+ *

4+ * The host's seeded perturbation is one normal deviate per *grid point*

5+ * (src/mgpu/noise.ts), so two sessions at different lmax seeded from the same

6+ * integer do not start from the same field — they start from unrelated fields

7+ * that merely share a random seed. Comparing them would compare two different

8+ * problems, and every number the comparison produced would be meaningless.

9+ *

10+ * So the field is built once, band-limited at the *coarsest* variant's lmax,

11+ * and evaluated on each variant's own grid:

12+ *

13+ * 1. white noise on the coarsest grid

14+ * 2. analysed there -> coefficients up to lmax_min

15+ * 3. zero-padded into each variant's coefficient layout

16+ * 4. synthesized on that variant's grid

17+ *

18+ * Steps 3 and 4 are exact: the field is band-limited at lmax_min, and every

19+ * variant's band contains that, so each one receives the *same function*

20+ * sampled where it needs it. Running each model's own `init` on it then leaves

21+ * every session holding the identical spectral state (zero-padded), which is

22+ * what makes a pointwise comparison at later times mean something.

23+ *

24+ * The coarsest variant gets the projected field too, not the raw white noise

25+ * it was analysed from — otherwise it alone would start somewhere slightly

26+ * different from the others.

27+ */

28+import { lmIndex, nlmCalc } from '../sht/layout.ts';

29+import { seededNoise } from '../mgpu/noise.ts';

30+import type { ModelSession } from '../mgpu/session.ts';

31+

32+/**

33+ * Re-index coefficients from a band limit into a wider one's layout, zero-

34+ * filling the degrees the source does not have. Both layouts are SHTNS

35+ * m-major with mmax = lmax, so nothing but the index mapping changes.

36+ */

37+export function prolongCoeffs(

38+ q: Float32Array,

39+ lmaxFrom: number,

40+ lmaxTo: number,

41+): Float32Array {

42+ if (lmaxTo === lmaxFrom) return q;

43+ if (lmaxTo < lmaxFrom) {

44+ throw new Error(`prolongCoeffs: cannot widen ${lmaxFrom} into a smaller ${lmaxTo}`);

45+ }

46+ const out = new Float32Array(2 * nlmCalc(lmaxTo, lmaxTo));

47+ for (let m = 0; m <= lmaxFrom; m++) {

48+ for (let l = m; l <= lmaxFrom; l++) {

49+ const from = 2 * lmIndex(lmaxFrom, l, m);

50+ const to = 2 * lmIndex(lmaxTo, l, m);

51+ out[to] = q[from];

52+ out[to + 1] = q[from + 1];

53+ }

54+ }

55+ return out;

56+}

57+

58+/**

59+ * The same band-limited perturbation, sampled on each session's grid. Order

60+ * follows `sessions`. Nothing may be in flight on any session's transform

61+ * plan — the one-off analys/synth here use the plan's own scratch buffers.

62+ */

63+export async function sharedNoise(

64+ sessions: ModelSession[],

65+ amp: number,

66+ seed: number,

67+): Promise<Float32Array[]> {

68+ let base = sessions[0];

69+ for (const s of sessions) if (s.cfg.lmax < base.cfg.lmax) base = s;

70+ const coeffs = await base.sht.analys(seededNoise(base.npts, amp, seed));

71+ const out: Float32Array[] = [];

72+ for (const s of sessions) {

73+ out.push(await s.sht.synth(prolongCoeffs(coeffs, base.cfg.lmax, s.cfg.lmax)));

74+ }

75+ return out;

76+}

src/compare/variants.tsadded+81−0View file

@@ -0,0 +1,81 @@

1+/**

2+ * One point of a convergence study: a choice of the three knobs that decide

3+ * *how well* the same problem is being solved, rather than what the problem is.

4+ *

5+ * niter iterations of the implicit solve (structural — it unrolls into the

6+ * compiled step, so each value is its own compiled session)

7+ * lmax the spectral band, and with it the grid (also structural)

8+ * dtDiv the timestep, as an integer divisor of the model's own dt

9+ *

10+ * dt is a *divisor* rather than a free value on purpose, and it is the whole

11+ * reason the comparison can be trusted: variants have to be compared at the

12+ * same model time, and with dt = dtBase/K every variant lands exactly on the

13+ * same t after K times as many steps — no rounding, no drift, no interpolation

14+ * in time. A free dt would put each variant on its own timeline and every

15+ * difference reported would be part real and part "these are 0.003 apart".

16+ */

17+

18+export interface Variant {

19+ /** Iterations of the implicit solve. */

20+ niter: number;

21+ /** Spectral band limit. */

22+ lmax: number;

23+ /** Timestep divisor: this variant runs at dtBase / dtDiv. */

24+ dtDiv: number;

25+}

26+

27+/** Stable identity of a variant, for keying maps and the reference <select>. */

28+export const variantKey = (v: Variant): string => `${v.niter}/${v.lmax}/${v.dtDiv}`;

29+

30+/** Human label. The dt term is dropped when nothing varies it, so the common

31+ * case (niter x lmax) reads as just those two. */

32+export const variantLabel = (v: Variant, showDt: boolean): string =>

33+ `niter ${v.niter} · lmax ${v.lmax}` + (showDt ? ` · dt/${v.dtDiv}` : '');

34+

35+/**

36+ * Every combination of the selected values, in a stable order: coarsest first,

37+ * so the grid reads top-to-bottom from least to most resolved and the

38+ * reference (the last row) is the one everything is measured against.

39+ */

40+export function crossProduct(

41+ niters: number[],

42+ lmaxes: number[],

43+ dtDivs: number[],

44+): Variant[] {

45+ const out: Variant[] = [];

46+ for (const lmax of [...lmaxes].sort((a, b) => a - b)) {

47+ for (const dtDiv of [...dtDivs].sort((a, b) => a - b)) {

48+ for (const niter of [...niters].sort((a, b) => a - b)) {

49+ out.push({ niter, lmax, dtDiv });

50+ }

51+ }

52+ }

53+ return out;

54+}

55+

56+/**

57+ * Index of the most-resolved variant: the natural reference, since it is the

58+ * one every other choice is an approximation of. Finer band first (it bounds

59+ * what can be represented at all), then more solve iterations, then smaller

60+ * timestep.

61+ */

62+export function mostResolved(variants: Variant[]): number {

63+ let best = 0;

64+ for (let i = 1; i < variants.length; i++) {

65+ const a = variants[i];

66+ const b = variants[best];

67+ if (

68+ a.lmax > b.lmax ||

69+ (a.lmax === b.lmax && a.niter > b.niter) ||

70+ (a.lmax === b.lmax && a.niter === b.niter && a.dtDiv > b.dtDiv)

71+ ) {

72+ best = i;

73+ }

74+ }

75+ return best;

76+}

77+

78+/** Distinguishable line/label colors, one per variant row. */

79+export const VARIANT_COLORS = [

80+ '#0969da', '#bf8700', '#1a7f37', '#cf222e', '#8250df', '#0f7c8a',

81+];

src/main.tsmodified+287−11View file

@@ -8,6 +8,7 @@ import { CodeEditor } from './editor/codeEditor.ts';

88 import {

99 formatCommand,

1010 resolvePreset,

11+ DEFAULT_NITER,

1112 DEFAULT_STEPS,

1213 DEFAULT_WARMUP,

1314 type RunSpec,

@@ -31,6 +32,14 @@ import { SphereScene } from './render/SphereScene.ts';

3132 import { Colorbar, fmtValue } from './render/colorbar.ts';

3233 import { colormaps, colormapNames } from './render/colormaps.ts';

3334 import { MovieRecorder } from './render/movie.ts';

35+import { CompareRun } from './compare/compareRun.ts';

36+import {

37+ crossProduct,

38+ mostResolved,

39+ variantKey,

40+ variantLabel,

41+ type Variant,

42+} from './compare/variants.ts';

3443

3544 const $ = <T extends HTMLElement>(id: string): T =>

3645 document.getElementById(id) as T;

@@ -52,6 +61,14 @@ const elMovieSpeed = $<HTMLSelectElement>('moviespeed');

5261 const elMovieRes = $<HTMLSelectElement>('movieres');

5362 const elMovieRotate = $<HTMLInputElement>('movierotate');

5463 const elMovie = $<HTMLButtonElement>('movie');

64+const elCompareToggle = $<HTMLButtonElement>('comparetoggle');

65+const elCompareBar = $('comparebar');

66+const elCmpNiter = $('cmp-niter');

67+const elCmpLmax = $('cmp-lmax');

68+const elCmpDt = $('cmp-dt');

69+const elCmpRef = $<HTMLSelectElement>('cmp-ref');

70+const elCmpStart = $<HTMLButtonElement>('cmp-start');

71+const elCmpCount = $('cmp-count');

5572 const elParams = $('params');

5673 const elGeomParams = $('geomparams');

5774 const elGeomNote = $('geomnote');

@@ -233,6 +250,10 @@ let generation = 0; // bumped on every rebuild to cancel stale pumps

233250 * re-synthesizing. */

234251 let coords: Float32Array | null = null;

235252 let posBuf: Float32Array | null = null;

253+/** The convergence study, when one is running; null in ordinary single-run

254+ * mode. While it is non-null there is no `session`: the study owns one per

255+ * variant, and the panels area is its grid. */

256+let compareRun: CompareRun | null = null;

236257

237258 const source = (): string => editedSource ?? model.source;

238259 const geomSource = (): string => editedGeomSource ?? geometry.source;

@@ -253,8 +274,11 @@ function buildParamInputs(): void {

253274 const v = Number(input.value);

254275 if (Number.isFinite(v)) params[spec.key] = v;

255276 // Parameters are uniforms, not constants baked into the kernels, so a

256- // change costs an upload rather than a recompile.

277+ // change costs an upload rather than a recompile. In compare mode `dt`

278+ // is the *base* timestep each variant's divisor divides, so the study

279+ // re-derives every variant's dt from it.

257280 session?.setParams(params);

281+ compareRun?.setParams(params);

258282 updateCommand();

259283 });

260284 label.append(input);

@@ -334,18 +358,25 @@ function showEditorFile(): void {

334358 }

335359 }

336360

337-/** The run currently on screen, as the benchmark's RunSpec. */

361+/**

362+ * The run currently on screen, as the benchmark's RunSpec. While a study is

363+ * running there is no single run, so this describes its *reference* variant —

364+ * the one the other rows are measured against, and the only one of them whose

365+ * numbers mean anything on their own.

366+ */

338367 function currentSpec(): RunSpec {

368+ const ref = compareRun?.variants[compareRefIndex()];

369+ const dt = ref ? { dt: (params.dt ?? 0) / ref.dtDiv } : null;

339370 return {

340371 preset: elModel.value,

341- lmax: Number(elLmax.value),

372+ lmax: ref ? ref.lmax : Number(elLmax.value),

342373 seed,

343374 steps: DEFAULT_STEPS,

344375 warmup: DEFAULT_WARMUP,

345- params,

376+ params: dt ? { ...params, ...dt } : params,

346377 geometry: geometry.key,

347378 geometryParams: geomParams,

348- niter: Number(elNiter.value),

379+ niter: ref ? ref.niter : Number(elNiter.value),

349380 };

350381 }

351382

@@ -366,6 +397,10 @@ elNiter.addEventListener('change', () => void rebuild());

366397 // one chain: a rapid second change waits its turn.

367398 let viewChange = Promise.resolve();

368399 elOversample.addEventListener('change', () => {

400+ // The study picks its own display grid — one grid common to every variant is

401+ // what makes their fields comparable — so this control is inert (and

402+ // disabled) while one is running.

403+ if (compareRun) return;

369404 viewChange = viewChange.then(() => applyOversample());

370405 });

371406 elGeometry.addEventListener('change', () => {

@@ -375,14 +410,22 @@ elGeometry.addEventListener('change', () => {

375410 // Morph is pure rendering: no readback, no GPU work, just the vertex buffer.

376411 elMorph.addEventListener('input', () => {

377412 morph = Number(elMorph.value);

378- applyMorph();

413+ if (compareRun) compareRun.setMorph(morph);

414+ else applyMorph();

415+});

416+elColormap.addEventListener('change', () => {

417+ if (compareRun) void compareRun.draw();

418+ else void draw();

379419 });

380-elColormap.addEventListener('change', () => void draw());

381420 elEditorFile.addEventListener('change', () => showEditorFile());

382421

383422 function setRunning(next: boolean): void {

384423 running = next;

385424 elRunPause.textContent = running ? 'Pause' : 'Run';

425+ if (compareRun) {

426+ compareRun.setRunning(next);

427+ return;

428+ }

386429 if (running) void pump();

387430 }

388431

@@ -395,6 +438,7 @@ elReseed.addEventListener('click', () => {

395438 void reseed();

396439 });

397440 elResetView.addEventListener('click', () => {

441+ compareRun?.resetView();

398442 for (const s of scenes) s.resetCamera();

399443 });

400444 elMovieToggle.addEventListener('click', () => {

@@ -547,6 +591,10 @@ async function applyOversample(): Promise<void> {

547591 * can be changed mid-run. Only the mesh is rebuilt.

548592 */

549593 async function applyGeometry(): Promise<void> {

594+ // The in-place swap below is a single session's trick. Each variant carries

595+ // the surface band-limited at its own lmax, and the study's meshes are built

596+ // from those, so a shape change goes through the full rebuild instead.

597+ if (compareRun) return rebuildCompare();

550598 if (!session) return;

551599 const gen = generation;

552600 const wasRunning = running;

@@ -583,14 +631,17 @@ function applyMorph(): void {

583631

584632 /** What the surface is, and the standing caveat about where it is not. */

585633 function updateGeomNote(): void {

586- if (!session) {

634+ // In compare mode each variant carries the surface band-limited at its own

635+ // lmax; the reference's is the one quoted, as everywhere else.

636+ const s = session ?? compareRun?.referenceSession ?? null;

637+ if (!s) {

587638 elGeomNote.textContent = '';

588639 return;

589640 }

590- const { lo, hi } = session.geometry.radiusRange();

591- const isSphere = session.geometryModel.key === SPHERE_KEY;

641+ const { lo, hi } = s.geometry.radiusRange();

642+ const isSphere = s.geometryModel.key === SPHERE_KEY;

592643 elGeomNote.innerHTML =

593- `<b>${session.geometryModel.label}</b> — ${session.geometryModel.blurb} ` +

644+ `<b>${s.geometryModel.label}</b> — ${s.geometryModel.blurb} ` +

594645 `Radius ${lo.toFixed(3)}–${hi.toFixed(3)}.` +

595646 (isSphere ? '' : ' <b>Rendered only</b> — not yet in the operator.');

596647 }

@@ -605,6 +656,10 @@ function reportCompileError(e: unknown): void {

605656 }

606657

607658 async function rebuild(): Promise<void> {

659+ // A study is several runs, so "rebuild the run" means rebuild all of them.

660+ // Everything that recompiles — a model or preset change, an edit to either

661+ // .m, a revert — arrives here, and none of it needs to know which mode is up.

662+ if (compareRun) return rebuildCompare();

608663 generation++;

609664 const gen = generation;

610665 setRunning(false);

@@ -670,6 +725,9 @@ async function rebuild(): Promise<void> {

670725 }

671726

672727 async function reseed(): Promise<void> {

728+ // One new perturbation for the whole study, band-limited at its coarsest

729+ // variant and evaluated on each grid — see src/compare/sharedStart.ts.

730+ if (compareRun) return compareRun.reseed(seed);

673731 if (!session) return;

674732 const gen = generation;

675733 session.seed(seed);

@@ -1029,9 +1087,227 @@ async function recordMovie(): Promise<void> {

10291087 }

10301088 }

10311089

1090+// ---------------------------------------------------------------- compare

1091+/**

1092+ * Comparing several solver settings at once.

1093+ *

1094+ * Deliberately a mode rather than a widening of the ordinary controls: the

1095+ * single-run path above is untouched, and with the bar closed nothing about

1096+ * using this page has changed. Opening it and pressing Compare tears down the

1097+ * one session and hands the panels area to a CompareRun, which owns a session

1098+ * per variant; pressing it again puts the single run back.

1099+ *

1100+ * The ceilings below are not arbitrary. Each variant compiles its whole

1101+ * unrolled step with no pipeline cache between sessions (a solve iteration is

1102+ * ~15 kernels per species), so the variant count is what you wait for; and

1103+ * each panel is a WebGL context and a full mesh, so the panel count is what

1104+ * the browser has to keep alive at once.

1105+ */

1106+const MAX_VARIANTS = 6;

1107+const MAX_PANELS = 12;

1108+/** dt divisors. Powers of two so that dtBase/K is exact in binary and every

1109+ * variant lands on the same model time with no accumulated drift. */

1110+const DT_DIVISORS = [1, 2, 4, 8];

1111+

1112+/**

1113+ * What the bar opens on: the default iteration count against the next step up,

1114+ * at the default band. Two variants, so the first study is quick to compile,

1115+ * and it asks the question the control exists for — is the default already

1116+ * converged? A flat, low curve says yes; one that climbs says the answer is

1117+ * still moving at niter 8 and the default is not enough for this shape.

1118+ */

1119+const cmpSelected = {

1120+ niter: new Set<number>([DEFAULT_NITER, 2 * DEFAULT_NITER]),

1121+ lmax: new Set<number>([63]),

1122+ dt: new Set<number>([1]),

1123+};

1124+

1125+/** A row of toggle chips backed by a Set. At least one stays selected — an

1126+ * empty axis has no meaning here, and silently falling back to a default

1127+ * would hide which values are actually being run. */

1128+function buildChips(host: HTMLElement, values: number[], selected: Set<number>, label: (v: number) => string): void {

1129+ host.replaceChildren();

1130+ for (const value of values) {

1131+ const chip = document.createElement('button');

1132+ chip.type = 'button';

1133+ chip.className = 'chip';

1134+ chip.textContent = label(value);

1135+ const paint = (): void => chip.setAttribute('aria-pressed', String(selected.has(value)));

1136+ paint();

1137+ chip.addEventListener('click', () => {

1138+ if (selected.has(value)) {

1139+ if (selected.size === 1) return;

1140+ selected.delete(value);

1141+ } else {

1142+ selected.add(value);

1143+ }

1144+ paint();

1145+ refreshVariants();

1146+ });

1147+ host.append(chip);

1148+ }

1149+}

1150+

1151+const cmpVariants = (): Variant[] =>

1152+ crossProduct([...cmpSelected.niter], [...cmpSelected.lmax], [...cmpSelected.dt]);

1153+

1154+/** The reference the user picked, clamped to the current variant list. */

1155+let cmpRefKey = '';

1156+

1157+/** Index of the reference in the current variant list, never negative. */

1158+function compareRefIndex(): number {

1159+ const i = cmpVariants().map(variantKey).indexOf(cmpRefKey);

1160+ return i < 0 ? 0 : i;

1161+}

1162+

1163+function refreshVariants(): void {

1164+ const variants = cmpVariants();

1165+ const showDt = cmpSelected.dt.size > 1;

1166+ const panels = variants.length * model.species.length;

1167+

1168+ const prev = cmpRefKey;

1169+ elCmpRef.replaceChildren();

1170+ for (const v of variants) {

1171+ const o = document.createElement('option');

1172+ o.value = variantKey(v);

1173+ o.textContent = variantLabel(v, showDt);

1174+ elCmpRef.append(o);

1175+ }

1176+ const keys = variants.map(variantKey);

1177+ cmpRefKey = keys.includes(prev) ? prev : keys[mostResolved(variants)];

1178+ elCmpRef.value = cmpRefKey;

1179+

1180+ const tooMany =

1181+ variants.length > MAX_VARIANTS

1182+ ? `${variants.length} variants — at most ${MAX_VARIANTS}`

1183+ : panels > MAX_PANELS

1184+ ? `${panels} panels — at most ${MAX_PANELS}`

1185+ : '';

1186+ elCmpCount.textContent = tooMany

1187+ ? `too many: ${tooMany}`

1188+ : `${variants.length} variants × ${model.species.length} species = ${panels} panels`;

1189+ elCmpCount.style.color = tooMany ? '#b35900' : '';

1190+ elCmpStart.disabled = tooMany !== '' && compareRun === null;

1191+}

1192+

1193+buildChips(

1194+ elCmpNiter,

1195+ [...elNiter.options].map((o) => Number(o.value)),

1196+ cmpSelected.niter,

1197+ String,

1198+);

1199+buildChips(

1200+ elCmpLmax,

1201+ [...elLmax.options].map((o) => Number(o.value)),

1202+ cmpSelected.lmax,

1203+ String,

1204+);

1205+buildChips(elCmpDt, DT_DIVISORS, cmpSelected.dt, (v) => (v === 1 ? 'dt' : `dt/${v}`));

1206+refreshVariants();

1207+

1208+elCmpRef.addEventListener('change', () => {

1209+ cmpRefKey = elCmpRef.value;

1210+ if (compareRun) void rebuildCompare();

1211+});

1212+

1213+elCompareToggle.addEventListener('click', () => {

1214+ elCompareBar.hidden = !elCompareBar.hidden;

1215+});

1216+

1217+elCmpStart.addEventListener('click', () => {

1218+ if (compareRun) void stopCompare();

1219+ else void startCompare();

1220+});

1221+

1222+/** Controls the study supersedes or cannot honour while it is running. */

1223+function setCompareUi(on: boolean): void {

1224+ for (const el of [elNiter, elLmax, elOversample, elBenchmark, elMovieToggle]) {

1225+ el.disabled = on;

1226+ }

1227+ elCmpNiter.querySelectorAll('button').forEach((b) => (b.disabled = on));

1228+ elCmpLmax.querySelectorAll('button').forEach((b) => (b.disabled = on));

1229+ elCmpDt.querySelectorAll('button').forEach((b) => (b.disabled = on));

1230+ elCmpStart.textContent = on ? 'Stop comparing' : 'Compare';

1231+ elCompareToggle.textContent = on ? 'Comparing' : 'Compare';

1232+ if (on) elMovieBar.hidden = true;

1233+}

1234+

1235+async function startCompare(): Promise<void> {

1236+ if (compareRun || !device) return;

1237+ const variants = cmpVariants();

1238+ if (variants.length > MAX_VARIANTS || variants.length * model.species.length > MAX_PANELS) {

1239+ return;

1240+ }

1241+ // Take down the single run first: its pump, its scenes, its session. The

1242+ // generation bump makes any readback already in flight drop its result.

1243+ generation++;

1244+ setRunning(false);

1245+ while (pumping) await nextFrame();

1246+ disposeView();

1247+ session?.destroy();

1248+ session = null;

1249+ elBenchResult.textContent = '';

1250+ elErr.textContent = '';

1251+ setCompareUi(true);

1252+

1253+ try {

1254+ compareRun = await CompareRun.create({

1255+ device,

1256+ model,

1257+ params,

1258+ source: source(),

1259+ geometry,

1260+ geometryParams: geomParams,

1261+ geometrySource: geomSource(),

1262+ variants,

1263+ reference: compareRefIndex(),

1264+ seed,

1265+ morph,

1266+ colormapName: () => elColormap.value,

1267+ container: elPanels,

1268+ onStatus: (html) => (elStats.innerHTML = html),

1269+ });

1270+ } catch (e) {

1271+ compareRun = null;

1272+ setCompareUi(false);

1273+ refreshVariants();

1274+ reportCompileError(e);

1275+ await rebuild();

1276+ return;

1277+ }

1278+ updateGeomNote();

1279+ // The command describes the reference variant, which only exists now.

1280+ updateCommand();

1281+ elRunPause.textContent = 'Run';

1282+}

1283+

1284+async function stopCompare(): Promise<void> {

1285+ if (!compareRun) return;

1286+ compareRun.dispose();

1287+ compareRun = null;

1288+ setCompareUi(false);

1289+ refreshVariants();

1290+ elStats.textContent = '';

1291+ await rebuild();

1292+}

1293+

1294+/** Rebuild the study in place — after a model, geometry, source or reference

1295+ * change. Same teardown as stopping, without leaving the mode. */

1296+async function rebuildCompare(): Promise<void> {

1297+ if (!compareRun) return;

1298+ compareRun.dispose();

1299+ compareRun = null;

1300+ setCompareUi(false);

1301+ await startCompare();

1302+}

1303+

10321304 // ---------------------------------------------------------------- boot

10331305 async function boot(): Promise<void> {

10341306 elModel.value = presets[0].key;

1307+ // The iteration count is one default shared with the benchmark, like the

1308+ // rest of the RunSpec's — take it from there rather than from the markup, so

1309+ // the page and `npm run bench` cannot start out disagreeing about it.

1310+ elNiter.value = String(DEFAULT_NITER);

10351311 elGeometry.value = DEFAULT_GEOMETRY_KEY;

10361312 elMorph.value = String(morph);

10371313 applyGeometryChoice(DEFAULT_GEOMETRY_KEY);

src/mgpu/session.tsmodified+42−12View file

@@ -224,25 +224,55 @@ export class ModelSession {

224224 */

225225 async setOversample(oversample: number): Promise<void> {

226226 const os = Math.max(1, Math.round(oversample));

227- if (os === this.#oversample) return;

228- const next =

229- os > 1

230- ? await ShtPlan.create(this.device, {

231- lmax: this.cfg.lmax,

232- mmax: this.cfg.mmax,

233- nlat: os * this.cfg.nlat,

234- nphi: os * this.cfg.nphi,

235- })

236- : null;

227+ await this.setDisplayGrid(os * this.cfg.nlat, os * this.cfg.nphi);

228+ }

229+

230+ /**

231+ * Point the display plan at an arbitrary grid, rather than an integer

232+ * multiple of the solver's. Same contract as setOversample — display-only,

233+ * no readback may be in flight — and the same exactness argument, which does

234+ * not care about the ratio: the state is band-limited at lmax, so

235+ * synthesizing it anywhere is evaluation, not resampling. What this adds is a

236+ * grid that need not be *finer*: several sessions at different lmax can be

237+ * put on one common grid, which is what makes their fields directly

238+ * comparable point by point and lets one mesh serve all of them.

239+ */

240+ async setDisplayGrid(nlat: number, nphi: number): Promise<void> {

241+ const view = this.viewSht.cfg;

242+ if (nlat === view.nlat && nphi === view.nphi) return;

243+ const onSolverGrid = nlat === this.cfg.nlat && nphi === this.cfg.nphi;

244+ const next = onSolverGrid

245+ ? null

246+ : await ShtPlan.create(this.device, {

247+ lmax: this.cfg.lmax,

248+ mmax: this.cfg.mmax,

249+ nlat,

250+ nphi,

251+ });

237252 const old = this.#displaySht;

238253 this.#displaySht = next;

239- this.#oversample = os;

254+ this.#oversample = nlat / this.cfg.nlat;

240255 old?.destroy();

241256 }

242257

243258 /** Run `init` from a seeded perturbation, resetting model time. */

244259 seed(seed: number): void {

245- this.gpu.init(seededNoise(this.npts, this.model.seedAmp, seed));

260+ this.seedWith(seededNoise(this.npts, this.model.seedAmp, seed));

261+ }

262+

263+ /**

264+ * Run `init` from a caller-supplied perturbation field, resetting model time.

265+ * `seed()` is this with the field the host's RNG produces on this session's

266+ * grid; supplying the field instead is how several sessions on *different*

267+ * grids can be started from the same band-limited initial condition, which is

268+ * the only way a comparison across lmax compares one problem rather than two

269+ * (see src/compare/sharedStart.ts).

270+ */

271+ seedWith(noise: Float32Array): void {

272+ if (noise.length !== this.npts) {

273+ throw new Error(`seedWith: noise must have length ${this.npts} (got ${noise.length})`);

274+ }

275+ this.gpu.init(noise);

246276 this.t = 0;

247277 this.steps = 0;

248278 }

test/compareChecks.tsadded+180−0View file

@@ -0,0 +1,180 @@

1+/**

2+ * The two things a side-by-side comparison of solver settings rests on.

3+ *

4+ * Both are silent when broken: the panels still animate, the difference norm

5+ * still produces a number, and the number is simply wrong — it reports a

6+ * disagreement between two runs that were never solving the same problem, or

7+ * that were never at the same time. Neither failure looks like a failure, which

8+ * is exactly why they are pinned here.

9+ *

10+ * 1. One initial condition. Sessions at different lmax seeded through

11+ * sharedNoise hold the *same* spectral state, zero-padded — and the

12+ * control shows what that is owed to: seeded the ordinary per-grid way,

13+ * the same integer seed gives two unrelated fields.

14+ *

15+ * 2. One clock. dt varies by a power-of-two divisor, so `steps * dt` is

16+ * bit-identical across variants and no comparison is ever made across a

17+ * fraction of a timestep.

18+ *

19+ * Deliberately small — two sessions at niter 1, lmax 31 and 63 — because a

20+ * session compiles its whole unrolled step and this suite has to stay short.

21+ */

22+import { ModelSession } from '../src/mgpu/session.ts';

23+import { mModelByKey, defaultParams } from '../src/mgpu/registry.ts';

24+import { prolongCoeffs, sharedNoise } from '../src/compare/sharedStart.ts';

25+import { lmIndex, nlmCalc } from '../src/sht/layout.ts';

26+import { crossProduct, mostResolved } from '../src/compare/variants.ts';

27+

28+type Check = (name: string, ok: boolean, detail: string) => void;

29+type Log = (line: string) => void;

30+

31+const COARSE = 31;

32+const FINE = 63;

33+

34+export async function compareChecks(

35+ device: GPUDevice,

36+ check: Check,

37+ log: Log,

38+): Promise<void> {

39+ log('\ncompare mode (convergence study):');

40+

41+ // ---- prolongCoeffs: every (l, m) lands on itself -------------------------

42+ {

43+ const src = new Float32Array(2 * nlmCalc(COARSE, COARSE));

44+ for (let m = 0; m <= COARSE; m++) {

45+ for (let l = m; l <= COARSE; l++) {

46+ const i = 2 * lmIndex(COARSE, l, m);

47+ src[i] = l + m / 100;

48+ src[i + 1] = -l - m / 100;

49+ }

50+ }

51+ const out = prolongCoeffs(src, COARSE, FINE);

52+ let moved = 0;

53+ let leaked = 0;

54+ for (let m = 0; m <= FINE; m++) {

55+ for (let l = m; l <= FINE; l++) {

56+ const j = 2 * lmIndex(FINE, l, m);

57+ if (l <= COARSE && m <= COARSE) {

58+ const i = 2 * lmIndex(COARSE, l, m);

59+ if (out[j] !== src[i] || out[j + 1] !== src[i + 1]) moved++;

60+ } else if (out[j] !== 0 || out[j + 1] !== 0) {

61+ leaked++;

62+ }

63+ }

64+ }

65+ check(

66+ 'compare: prolongation puts every coefficient at its own (l, m)',

67+ moved === 0 && leaked === 0,

68+ `${moved} misplaced, ${leaked} non-zero above the source band ` +

69+ `(${nlmCalc(COARSE, COARSE)} -> ${nlmCalc(FINE, FINE)} coefficients)`,

70+ );

71+ }

72+

73+ // ---- one initial condition across lmax, and the control ------------------

74+ {

75+ const model = mModelByKey('schnakenberg')!;

76+ const params = defaultParams(model);

77+ const sessions: ModelSession[] = [];

78+ try {

79+ for (const lmax of [COARSE, FINE]) {

80+ sessions.push(await ModelSession.create({ device, model, params, lmax, niter: 1 }));

81+ }

82+ const [coarse, fine] = sessions;

83+

84+ // What the study does: one band-limited field, evaluated on each grid.

85+ const noise = await sharedNoise(sessions, model.seedAmp, 1);

86+ sessions.forEach((s, i) => s.seedWith(noise[i]));

87+ const shared = compareStates(

88+ prolongCoeffs(await coarse.read('U'), COARSE, FINE),

89+ await fine.read('U'),

90+ );

91+ check(

92+ 'compare: one shared perturbation gives both grids the same state',

93+ shared.rel < 5e-5,

94+ `max |dU| = ${shared.abs.toExponential(2)} ` +

95+ `(${(100 * shared.rel).toFixed(4)}% of max |U| = ${shared.scale.toFixed(3)}) ` +

96+ `across lmax ${COARSE} vs ${FINE}`,

97+ );

98+

99+ // The control: the ordinary per-grid seeding these two would otherwise

100+ // get. One deviate per grid point, and the grids differ, so the same

101+ // integer seed is two different initial conditions -- comparing runs

102+ // started this way would report a difference that is entirely the seed.

103+ sessions.forEach((s) => s.seed(1));

104+ const plain = compareStates(

105+ prolongCoeffs(await coarse.read('U'), COARSE, FINE),

106+ await fine.read('U'),

107+ );

108+ check(

109+ 'compare: control — the same integer seed alone does not do it',

110+ plain.abs > 20 * shared.abs,

111+ `per-grid seeding differs by ${plain.abs.toExponential(2)}, ` +

112+ `${(plain.abs / Math.max(shared.abs, 1e-30)).toExponential(1)}x the shared start's ` +

113+ `(seed amplitude ${model.seedAmp})`,

114+ );

115+ log(

116+ ` shared start: ${shared.abs.toExponential(2)}, ` +

117+ `per-grid seeds: ${plain.abs.toExponential(2)}`,

118+ );

119+ } finally {

120+ for (const s of sessions) s.destroy();

121+ }

122+ }

123+

124+ // ---- one clock: steps * dt is bit-identical across the divisors ----------

125+ {

126+ const divisors = [1, 2, 4, 8];

127+ const steps = 4;

128+ let worst = 0;

129+ const cases: string[] = [];

130+ for (const model of ['schnakenberg', 'brusselator', 'allencahn']) {

131+ const dt = defaultParams(mModelByKey(model)!).dt;

132+ for (const div of divisors) {

133+ // A variant at dt/div takes div times as many steps to cover the same

134+ // span. Powers of two only touch the exponent, so both the divide and

135+ // the multiply back are exact and the two spans are the same float.

136+ const span = (steps * div) * (dt / div);

137+ const ulps = Math.abs(span - steps * dt);

138+ if (ulps > worst) worst = ulps;

139+ if (div === divisors[divisors.length - 1]) {

140+ cases.push(`${model} dt ${dt} -> ${dt / div}`);

141+ }

142+ }

143+ }

144+ check(

145+ 'compare: a power-of-two dt divisor keeps every variant on one clock',

146+ worst === 0,

147+ `exact for every shipped dt x ${divisors.join('/')} (${cases.join(', ')})`,

148+ );

149+ }

150+

151+ // ---- the variant grid and its reference ---------------------------------

152+ {

153+ const variants = crossProduct([1, 4], [31, 63], [1, 2]);

154+ const ref = variants[mostResolved(variants)];

155+ check(

156+ 'compare: the reference is the most-resolved corner of the grid',

157+ variants.length === 8 &&

158+ ref.niter === 4 && ref.lmax === 63 && ref.dtDiv === 2 &&

159+ new Set(variants.map((v) => `${v.niter}/${v.lmax}/${v.dtDiv}`)).size === 8,

160+ `${variants.length} distinct variants, reference niter ${ref.niter} · ` +

161+ `lmax ${ref.lmax} · dt/${ref.dtDiv}`,

162+ );

163+ }

164+}

165+

166+/** Max absolute difference of two equal-length spectral states, and that

167+ * difference relative to the scale of the reference. */

168+function compareStates(

169+ a: Float32Array,

170+ b: Float32Array,

171+): { abs: number; rel: number; scale: number } {

172+ let abs = 0;

173+ let scale = 0;

174+ const n = Math.min(a.length, b.length);

175+ for (let i = 0; i < n; i++) {

176+ abs = Math.max(abs, Math.abs(a[i] - b[i]));

177+ scale = Math.max(scale, Math.abs(b[i]));

178+ }

179+ return { abs, rel: scale > 0 ? abs / scale : Infinity, scale };

180+}

test/test-page.tsmodified+2−0View file

@@ -28,6 +28,7 @@ import { analyticChecks } from './analyticChecks.ts';

2828 import { modelChecks } from './modelChecks.ts';

2929 import { geometryChecks } from './geometryChecks.ts';

3030 import { fluxChecks } from './fluxChecks.ts';

31+import { compareChecks } from './compareChecks.ts';

3132

3233 declare global {

3334 interface Window {

@@ -220,6 +221,7 @@ async function main(): Promise<void> {

220221 // but minutes in a browser, where each session recompiles its unrolled step.

221222 await geometryChecks(device, check, log, { sweep: q.has('sweep') });

222223 await fluxChecks(device, check, log, { ab: q.has('sweep') });

224+ await compareChecks(device, check, log);

223225

224226 window.__RESULTS__ = { ok: failures === 0, lines };

225227 log(failures === 0 ? 'ALL PASS' : `${failures} FAILURE(S)`);