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 */
31import { ModelSession } from '../mgpu/session.ts';
32import type { MModel, Params } from '../mgpu/registry.ts';
33import type { MGeometry } from '../geom/registry.ts';
34import {
35 buildTopology,
36 fillFieldValues,
37 fillPositions,
38 fillColors,
39 type SphereMeshTopology,
40} from '../render/sphereMesh.ts';
41import { SphereScene } from '../render/SphereScene.ts';
42import { colormaps } from '../render/colormaps.ts';
ef3ae33Do not stretch the colormap across a constant field's roundoffJeremy Magland 43import { fmtValue, floorRange } from '../render/colorbar.ts';
c90d0e2Check reference files in the browser's compare modeJeremy Magland 44import { prolongCoeffs, sharedModes, sharedNoise } from './sharedStart.ts';
f467ffaCompare several solver settings side by side, on one clockJeremy Magland 45import { variantLabel, VARIANT_COLORS, type Variant } from './variants.ts';
c90d0e2Check reference files in the browser's compare modeJeremy Magland 46import type { ReferenceCase } from './referenceCase.ts';
48/**
49 * Latitudes of the shared display grid. 256 is the same target the single-run
50 * view uses for 'auto' oversampling, and for the same reason — beyond it a
51 * finer mesh costs vertices without showing anything.
52 *
53 * Here it is a ceiling as well as a target, in two directions. At lmax 255 the
54 * solver grid is finer than this, so the panels sample the (exact) state more
55 * coarsely than the solver carries it; and past a handful of panels the mesh is
56 * paid for once per panel, in vertices, normals and a WebGL context each, so it
57 * halves. Both are display choices, both are reported in the status line, and
58 * neither touches the difference norm's meaning: that is computed on this same
59 * grid for every variant, so it stays a consistent comparison whatever the grid.
60 */
61const RENDER_NLAT = 256;
62const RENDER_NLAT_CROWDED = 128;
63const CROWDED_PANELS = 6;
65/** See main.ts's DISPATCH_BUDGET — the same watchdog argument, per variant. */
66const DISPATCH_BUDGET = 1000;
67const STEPS_PER_FRAME_BASE = 4;
69export interface CompareOptions {
70 device: GPUDevice;
71 model: MModel;
72 /** The model's parameters, with `dt` read as the *base* timestep that each
73 * variant's dtDiv divides. */
74 params: Params;
75 source: string;
76 geometry: MGeometry;
77 geometryParams: Params;
78 geometrySource: string;
79 variants: Variant[];
c90d0e2Check reference files in the browser's compare modeJeremy Magland 80 /** Index into `variants` of the run everything else is measured against.
81 * Ignored when `refFile` is given — the file is the reference then. */
f467ffaCompare several solver settings side by side, on one clockJeremy Magland 82 reference: number;
84 * Check against a reference file instead of against each other: its exact
85 * initial state seeds every variant (so `seed` and `lam3` go unused), a
86 * static extra row shows its final state, every Δ is measured against that
87 * row, and the clock stops at the file's end time. Every variant's lmax must
88 * be >= the file's — a narrower band could not hold the initial state.
89 */
90 refFile?: ReferenceCase;
91 /** Called when a refFile run reaches the file's end time and stops. */
92 onFinished?: () => void;
beac00aMerge main into random-fieldsJeremy Magland 94 /** Wavelength of the seeded random field, shared by every variant — one
95 * initial condition means one wavelength as much as one seed. */
96 lam3?: number;
98 colormapName: () => string;
99 /** Where the variant grid goes (the app's #panels). */
100 container: HTMLElement;
101 /** Progress and, afterwards, the standing description of the study. */
102 onStatus: (html: string) => void;
103}
105interface Row {
106 variant: Variant;
107 session: ModelSession;
108 color: string;
109 /** Surface coordinates on the shared render grid — this variant's own. */
110 coords: Float32Array;
111 posBuf: Float32Array;
112 scenes: SphereScene[];
113 valueBufs: Float32Array[];
114 colorBufs: Float32Array[];
115 /** Fields read this frame, one per species, on the shared grid. */
116 fields: Float32Array[];
117 /** Relative difference from the reference, one per species. */
118 err: number[];
119 /** False once any species has left the floating-point numbers — the shape a
120 * variant outside the convergence radius eventually fails in. Such a row is
121 * never used to scale a column, and its label says so. */
122 healthy: boolean;
123 statEl: HTMLElement;
124}
127 * The reference file's final state, as one more row of panels — with no
128 * session behind it: its surface and fields are the file's coefficients
129 * synthesized once on the shared display grid, fixed for the whole run. Only
130 * its coloring changes, with the shared range.
131 */
132interface FileRow {
133 coords: Float32Array;
134 posBuf: Float32Array;
135 scenes: SphereScene[];
136 valueBufs: Float32Array[];
137 colorBufs: Float32Array[];
138 /** The file's final state on the shared grid, one per species. */
139 fields: Float32Array[];
140 /** Its extent, precomputed — a candidate for the shared color range. */
141 bounds: (Bounds | null)[];
142}
f467ffaCompare several solver settings side by side, on one clockJeremy Magland 144export class CompareRun {
145 #opts: CompareOptions;
146 #rows: Row[] = [];
c90d0e2Check reference files in the browser's compare modeJeremy Magland 147 #fileRow: FileRow | null = null;
ca37955Adding a way to reset simulation from same random IC.Owen Melia 148 /** What restart() reloads: the file's fixed state if opts.refFile is set,
149 * otherwise a snapshot of the coarsest variant's state as of the last
150 * (re-)seed — see the capture in create() and in reseed()'s plain branch. */
151 #initial: Record<string, Float32Array>;
152 #initialLmax: number;
c90d0e2Check reference files in the browser's compare modeJeremy Magland 153 /** Base steps taken since the initial state — the refFile clock. */
154 #stepsDone = 0;
155 /** True once a refFile run has reached the file's end time. */
156 #finished = false;
f467ffaCompare several solver settings side by side, on one clockJeremy Magland 157 #topo: SphereMeshTopology;
158 /** Quadrature weight per grid point of the shared grid, for the L2 norm. */
159 #weights: Float64Array;
160 #rangeBars: { fill: (lo: number, hi: number) => void }[] = [];
161 /** Smoothed color range per species, shared by every variant so the panels
162 * in a column are directly comparable by eye and not just by number. */
163 #ranges: { lo: number; hi: number }[] = [];
164 #resizeObs: ResizeObserver | null = null;
166 #running = false;
167 #pumping = false;
168 #disposed = false;
169 #morph: number;
170 /** Base steps per frame; variant i takes this times its dtDiv. */
171 #frameSteps = STEPS_PER_FRAME_BASE;
172 /** Model time all variants are at — one number, by construction. */
173 #t = 0;
174 #frameMs = 0;
175 #note: string;
177 private constructor(init: {
178 opts: CompareOptions;
179 rows: Row[];
c90d0e2Check reference files in the browser's compare modeJeremy Magland 180 fileRow: FileRow | null;
f467ffaCompare several solver settings side by side, on one clockJeremy Magland 181 topo: SphereMeshTopology;
182 weights: Float64Array;
183 rangeBars: { fill: (lo: number, hi: number) => void }[];
184 frameSteps: number;
185 note: string;
ca37955Adding a way to reset simulation from same random IC.Owen Melia 186 initial: Record<string, Float32Array>;
187 initialLmax: number;
189 this.#opts = init.opts;
190 this.#rows = init.rows;
c90d0e2Check reference files in the browser's compare modeJeremy Magland 191 this.#fileRow = init.fileRow;
f467ffaCompare several solver settings side by side, on one clockJeremy Magland 192 this.#topo = init.topo;
193 this.#weights = init.weights;
194 this.#rangeBars = init.rangeBars;
195 this.#frameSteps = init.frameSteps;
196 this.#note = init.note;
197 this.#morph = init.opts.morph;
198 this.#ranges = init.opts.model.species.map(() => ({ lo: NaN, hi: NaN }));
ca37955Adding a way to reset simulation from same random IC.Owen Melia 199 this.#initial = init.initial;
200 this.#initialLmax = init.initialLmax;
203 get variants(): Variant[] {
204 return this.#rows.map((r) => r.variant);
205 }
207 /** The variant everything else is measured against — the one whose numbers
208 * stand on their own, so the one the app quotes when it has to quote one. */
209 get referenceSession(): ModelSession | null {
210 return this.#rows[this.#opts.reference]?.session ?? null;
211 }
213 get referenceIndex(): number {
214 return this.#opts.reference;
215 }
c90d0e2Check reference files in the browser's compare modeJeremy Magland 217 /** The reference file this study is checking against, if any. */
218 get refFile(): ReferenceCase | null {
219 return this.#opts.refFile ?? null;
220 }
f467ffaCompare several solver settings side by side, on one clockJeremy Magland 222 /** The base timestep a variant's dtDiv divides. */
223 static baseDt(params: Params): number {
224 return params.dt ?? 0;
225 }
227 static async create(opts: CompareOptions): Promise<CompareRun> {
228 const { device, model, variants } = opts;
229 const baseDt = CompareRun.baseDt(opts.params);
230 const showDt = variants.some((v) => v.dtDiv !== variants[0].dtDiv);
231 const sessions: ModelSession[] = [];
232 // Scenes own a WebGL context and an animation frame each, so a failure
233 // after the grid is up has to take them down explicitly — removing their
234 // canvases from the DOM would leave both running.
235 let built: Row[] = [];
c90d0e2Check reference files in the browser's compare modeJeremy Magland 236 let builtFile: FileRow | null = null;
238 try {
239 for (let i = 0; i < variants.length; i++) {
240 const v = variants[i];
241 opts.onStatus(
242 `compiling ${i + 1}/${variants.length} — ${variantLabel(v, showDt)} ` +
243 `(a solve iteration is ~15 kernels per species, and there is no ` +
244 `pipeline cache across sessions)`,
245 );
246 // Yield, so the status actually paints before the compile blocks.
247 await new Promise<number>(requestAnimationFrame);
248 sessions.push(
249 await ModelSession.create({
250 device,
251 model,
252 params: { ...opts.params, dt: baseDt / v.dtDiv },
253 lmax: v.lmax,
254 source: opts.source,
255 geometry: opts.geometry,
256 geometryParams: opts.geometryParams,
257 geometrySource: opts.geometrySource,
258 niter: v.niter,
261 );
262 }
264 // ---- the shared display grid ----------------------------------------
265 const maxLmax = Math.max(...variants.map((v) => v.lmax));
c90d0e2Check reference files in the browser's compare modeJeremy Magland 266 const panels = (variants.length + (opts.refFile ? 1 : 0)) * model.species.length;
f467ffaCompare several solver settings side by side, on one clockJeremy Magland 267 const target = panels > CROWDED_PANELS ? RENDER_NLAT_CROWDED : RENDER_NLAT;
268 // Never below what the finest band needs to be representable at all
269 // (ShtPlan requires nlat > lmax), whatever the panel count says.
270 const nlat = Math.max(target, 2 * Math.ceil((maxLmax + 2) / 2));
271 let nphi = 1;
272 while (nphi < Math.max(2 * nlat, 2 * maxLmax + 1)) nphi *= 2;
273 for (const s of sessions) await s.setDisplayGrid(nlat, nphi);
275 // ---- one initial condition, on every grid ---------------------------
ca37955Adding a way to reset simulation from same random IC.Owen Melia 276 // Also what restart() reloads later — the file's fixed state, or (for
277 // the plain case) a snapshot of the coarsest variant's own state,
278 // taken after seeding it: the same lowest-lmax session sharedNoise
279 // itself draws from, so prolonging it up to any other variant later is
280 // always widening a band, never narrowing one.
281 let initial: Record<string, Float32Array>;
282 let initialLmax: number;
284 // The file's exact spectral state, prolonged into each variant's band.
285 // Exact, not approximate: the state is band-limited at the file's lmax
286 // and every variant's band contains it, so each session starts from
287 // the very field the reference run started from.
288 opts.onStatus('loading the initial state from the reference file…');
289 for (const s of sessions) {
290 s.loadState(prolongState(opts.refFile.initial, model.state, opts.refFile.lmax, s.cfg.lmax));
291 }
ca37955Adding a way to reset simulation from same random IC.Owen Melia 292 initial = opts.refFile.initial;
293 initialLmax = opts.refFile.lmax;
295 opts.onStatus('seeding all variants from one band-limited perturbation…');
296 const noise = await sharedNoise(sessions, model.seedAmp, opts.seed);
297 const modes = await sharedModes(sessions[opts.reference] ?? sessions[0], opts.seed);
298 // One at a time: a seed submits its whole mode sum in pieces, and there
299 // is nothing to gain from interleaving several variants' worth of it.
300 for (let i = 0; i < sessions.length; i++) await sessions[i].seedWith(noise[i], modes);
ca37955Adding a way to reset simulation from same random IC.Owen Melia 301 let coarsest = sessions[0];
302 for (const s of sessions) if (s.cfg.lmax < coarsest.cfg.lmax) coarsest = s;
303 initial = await coarsest.readState();
304 initialLmax = coarsest.cfg.lmax;
307 // ---- the mesh, shared; the surface, per variant ---------------------
308 const view = sessions[0].viewSht;
309 const phi = new Float64Array(nphi);
310 for (let j = 0; j < nphi; j++) phi[j] = (2 * Math.PI * j) / nphi;
311 const topo = buildTopology(view.cosTheta, phi);
312 // Gauss weights carry the sin(theta) of the area element; the constant
313 // 2*pi/nphi is common to every point and cancels in the relative norm.
314 const weights = new Float64Array(nlat * nphi);
315 for (let i = 0; i < nlat; i++) {
316 for (let j = 0; j < nphi; j++) weights[i * nphi + j] = view.gaussWeights[i];
317 }
319 // ---- how many steps a frame may submit ------------------------------
320 // Per variant: its own unrolled step size times its dtDiv, since a ÷K
321 // variant takes K times as many steps to reach the same time.
322 let frameSteps = STEPS_PER_FRAME_BASE;
323 const ops: number[] = [];
324 for (let i = 0; i < sessions.length; i++) {
325 const n = Math.max(1, sessions[i].describe().step.length);
326 ops.push(n);
327 frameSteps = Math.min(
328 frameSteps,
329 Math.max(1, Math.floor(DISPATCH_BUDGET / (n * variants[i].dtDiv))),
330 );
331 }
332 frameSteps = Math.max(1, frameSteps);
334 // ---- the grid of panels ---------------------------------------------
c90d0e2Check reference files in the browser's compare modeJeremy Magland 335 const { rows, fileRow, rangeBars } = await buildGrid(opts, sessions, topo, showDt);
339 const solverGrid = sessions.map((s) => `${s.cfg.nlat}×${s.cfg.nphi}`);
340 const note =
3210e7dOpen the reference comparison in one clickJeremy Magland 341 `${variants.length} variant${variants.length === 1 ? '' : 's'} · ` +
342 `display grid ${nlat}×${nphi}` +
f467ffaCompare several solver settings side by side, on one clockJeremy Magland 343 (sessions.some((s) => s.cfg.nlat > nlat)
344 ? ` (below the finest solver grid ${solverGrid[solverGrid.length - 1]} — display only)`
345 : '') +
346 ` · ${frameSteps} base step${frameSteps === 1 ? '' : 's'}/frame` +
347 ` · ops/step ${ops.join(', ')}`;
349 const run = new CompareRun({
ca37955Adding a way to reset simulation from same random IC.Owen Melia 350 opts, rows, fileRow, topo, weights, rangeBars, frameSteps, note, initial, initialLmax,
352 await run.draw();
353 run.#observeResize();
354 run.#status();
355 return run;
356 } catch (e) {
357 for (const r of built) for (const s of r.scenes) s.dispose();
c90d0e2Check reference files in the browser's compare modeJeremy Magland 358 for (const s of builtFile?.scenes ?? []) s.dispose();
f467ffaCompare several solver settings side by side, on one clockJeremy Magland 359 for (const s of sessions) s.destroy();
360 opts.container.replaceChildren();
361 opts.container.classList.remove('compare');
362 throw e;
363 }
364 }
366 // ------------------------------------------------------------------ state
367 setRunning(next: boolean): void {
368 this.#running = next;
369 if (next) void this.#pump();
370 }
372 get running(): boolean {
373 return this.#running;
374 }
c90d0e2Check reference files in the browser's compare modeJeremy Magland 376 /** Re-seed every variant from one new shared perturbation — or, against a
377 * reference file, restart from its initial state (there is nothing to
378 * draw; the seed is ignored). */
f467ffaCompare several solver settings side by side, on one clockJeremy Magland 379 async reseed(seed: number): Promise<void> {
380 const wasRunning = this.#running;
381 this.#running = false;
382 while (this.#pumping) await nextFrame();
383 if (this.#disposed) return;
beac00aMerge main into random-fieldsJeremy Magland 384 const sessions = this.#rows.map((r) => r.session);
c90d0e2Check reference files in the browser's compare modeJeremy Magland 385 const refFile = this.#opts.refFile;
386 if (refFile) {
387 for (const s of sessions) {
388 s.loadState(prolongState(refFile.initial, this.#opts.model.state, refFile.lmax, s.cfg.lmax));
389 }
390 } else {
391 const noise = await sharedNoise(sessions, this.#opts.model.seedAmp, seed);
392 const modes = await sharedModes(this.referenceSession ?? sessions[0], seed);
393 // Checked per variant, not once: a seed awaits its own submission, so a
394 // dispose can land between two of them and destroy the sessions left.
395 for (let i = 0; i < sessions.length; i++) {
396 if (this.#disposed) return;
397 await sessions[i].seedWith(noise[i], modes);
398 }
ca37955Adding a way to reset simulation from same random IC.Owen Melia 399 if (this.#disposed) return;
400 // This draw becomes what restart() rewinds to from now on — see the
401 // identical selection in create(). Recaptured here rather than left
402 // pointing at the pre-reseed field.
403 let coarsest = sessions[0];
404 for (const s of sessions) if (s.cfg.lmax < coarsest.cfg.lmax) coarsest = s;
405 this.#initial = await coarsest.readState();
406 this.#initialLmax = coarsest.cfg.lmax;
407 }
408 this.#t = 0;
409 this.#stepsDone = 0;
410 this.#finished = false;
411 for (const r of this.#ranges) {
412 r.lo = NaN;
413 r.hi = NaN;
414 }
415 await this.draw();
416 this.#status();
417 if (!this.#disposed && wasRunning) this.setRunning(true);
418 }
420 /** Rewind every variant to the saved initial condition — the file's fixed
421 * state, or (for the plain case) the last (re-)seed, not necessarily the
422 * very first one — without drawing anything new. */
423 async restart(): Promise<void> {
424 const wasRunning = this.#running;
425 this.#running = false;
426 while (this.#pumping) await nextFrame();
427 if (this.#disposed) return;
428 for (const r of this.#rows) {
429 r.session.loadState(
430 prolongState(this.#initial, this.#opts.model.state, this.#initialLmax, r.session.cfg.lmax),
431 );
435 this.#finished = false;
f467ffaCompare several solver settings side by side, on one clockJeremy Magland 436 for (const r of this.#ranges) {
437 r.lo = NaN;
438 r.hi = NaN;
439 }
440 await this.draw();
f467ffaCompare several solver settings side by side, on one clockJeremy Magland 442 if (!this.#disposed && wasRunning) this.setRunning(true);
443 }
beac00aMerge main into random-fieldsJeremy Magland 445 /** Wavelength of the seeded random field. One number for the study: every
446 * variant seeds from the same field, so they seed at the same wavelength. */
447 get lam3(): number {
448 return this.#rows[0]?.session.lam3 ?? 0;
449 }
451 /** Change it on every variant. Like the single run's, this only takes effect
452 * on the next reseed, which is where the field is drawn. */
453 setLam3(lambda: number): void {
454 this.#opts.lam3 = lambda;
455 for (const r of this.#rows) r.session.setLam3(lambda);
456 }
f467ffaCompare several solver settings side by side, on one clockJeremy Magland 458 /** Model parameters changed. Each variant keeps its own dt. */
459 setParams(params: Params): void {
c90d0e2Check reference files in the browser's compare modeJeremy Magland 460 // Against a reference file the parameters *are* the file's — they define
461 // the problem being checked — and the page's parameter panel edits the
462 // page's own model, which need not even be this one. Nothing to apply.
463 if (this.#opts.refFile) return;
f467ffaCompare several solver settings side by side, on one clockJeremy Magland 464 this.#opts.params = params;
465 const baseDt = CompareRun.baseDt(params);
466 for (const r of this.#rows) {
467 r.session.setParams({ ...params, dt: baseDt / r.variant.dtDiv });
468 }
469 }
471 setMorph(morph: number): void {
472 this.#morph = morph;
473 for (const r of this.#rows) {
474 fillPositions(r.posBuf, r.coords, this.#topo, morph);
475 for (const s of r.scenes) s.updatePositions(r.posBuf);
476 }
c90d0e2Check reference files in the browser's compare modeJeremy Magland 477 const f = this.#fileRow;
478 if (f) {
479 fillPositions(f.posBuf, f.coords, this.#topo, morph);
480 for (const s of f.scenes) s.updatePositions(f.posBuf);
481 }
484 resetView(): void {
c90d0e2Check reference files in the browser's compare modeJeremy Magland 485 for (const s of this.#allScenes()) s.resetCamera();
488 dispose(): void {
489 this.#disposed = true;
490 this.#running = false;
491 this.#resizeObs?.disconnect();
492 this.#resizeObs = null;
493 for (const r of this.#rows) {
494 for (const s of r.scenes) s.dispose();
495 r.session.destroy();
496 }
c90d0e2Check reference files in the browser's compare modeJeremy Magland 497 for (const s of this.#fileRow?.scenes ?? []) s.dispose();
f467ffaCompare several solver settings side by side, on one clockJeremy Magland 498 this.#rows = [];
f467ffaCompare several solver settings side by side, on one clockJeremy Magland 500 this.#opts.container.replaceChildren();
501 this.#opts.container.classList.remove('compare');
502 }
c90d0e2Check reference files in the browser's compare modeJeremy Magland 504 #allScenes(): SphereScene[] {
505 return [...this.#rows.flatMap((r) => r.scenes), ...(this.#fileRow?.scenes ?? [])];
506 }
f467ffaCompare several solver settings side by side, on one clockJeremy Magland 508 // ----------------------------------------------------------------- drawing
509 /**
510 * One frame's readback: every variant's every species, on the shared grid.
511 * Read first, then color — the range is shared down a column, so no panel can
512 * be filled until the column's range is known.
513 */
514 async draw(): Promise<void> {
515 if (this.#disposed) return;
516 const species = this.#opts.model.species;
517 // Sessions are independent, so their readbacks can be in flight together;
518 // within one session they must not be (they share its staging buffers).
519 await Promise.all(
520 this.#rows.map(async (r) => {
521 for (let k = 0; k < species.length; k++) {
522 r.fields[k] = await r.session.readSpecies(k);
523 }
524 }),
525 );
526 if (this.#disposed) return;
528 const cmap = colormaps[this.#opts.colormapName()] ?? colormaps.viridis;
530 /**
531 * What scales a column is the whole question, and it has three wrong
532 * answers.
533 *
534 * Per panel is wrong: a range each rescales every variant to itself and
535 * hides exactly the difference the grid exists to show. The union over
536 * variants is wrong for the opposite reason: a variant outside the
537 * iteration's convergence radius runs away to 1e20 and then to NaN, and a
538 * union range rescales the *whole column* to it, flattening every panel to
539 * one colour — which reads as "they all blew up" when only one did.
540 *
541 * The reference alone is wrong too, less obviously, and it is the case that
542 * actually bites: outside the convergence radius *more* Richardson
543 * iterations diverge *faster*, so the row that goes first is usually the
544 * highest-niter one — which is the reference.
545 *
546 * So the column is scaled by whichever variant **reaches least far from
547 * zero** — the least-blown-up one. That is a comparison between the rows,
548 * not a threshold on any of them, and the distinction is the whole point:
549 * any "is this value too big?" test has a window in which a diverging field
550 * is still under the limit, and for as long as that window lasts it drags
551 * the scale and flattens the grid, until it finally trips and everything
552 * springs back. A comparison has no such window — a run-away only has to be
553 * *larger* than a healthy row to stop setting the scale, which it is from
554 * its first bad step, and it stays larger no matter how many other rows go
555 * with it. One healthy variant is enough to keep the grid readable.
556 *
557 * The cost is a slight bias: among healthy variants the scale comes from
558 * the one with the smallest peak, so the others clip by however much they
559 * exceed it. They are approximations of the same solution, so that is a
560 * fraction of a percent, and the alternative is a display that a single
561 * divergence can take away.
562 */
563 const bounds = this.#rows.map((r) => species.map((_, k) => finiteRange(r.fields[k])));
564 this.#rows.forEach((r, i) => {
565 // A row with any non-finite value is out of the running entirely: its
566 // finite entries are whatever survived, and no rank over them means much.
567 r.healthy = species.every((_, k) => allFinite(r.fields[k]) && bounds[i][k] !== null);
568 });
570 for (let k = 0; k < species.length; k++) {
c90d0e2Check reference files in the browser's compare modeJeremy Magland 571 // The file row, when there is one, is a candidate like any healthy
572 // variant: early on the variants' small fields set the scale (it merely
573 // clips), and if every variant diverges it is the row that keeps the
574 // grid readable.
575 const anchor = leastPeak([
576 ...this.#rows.map((r, i) => (r.healthy ? bounds[i][k] : null)),
577 this.#fileRow?.bounds[k] ?? null,
578 ]);
f467ffaCompare several solver settings side by side, on one clockJeremy Magland 579 const range = this.#ranges[k];
580 if (anchor) {
581 if (!Number.isFinite(range.lo)) {
582 range.lo = anchor.lo;
583 range.hi = anchor.hi;
584 } else {
585 // Smooth in both directions so the shading evolves gently as the
586 // pattern grows, as the single-run view does.
587 const a = 0.15;
588 range.lo += a * (anchor.lo - range.lo);
589 range.hi += a * (anchor.hi - range.hi);
590 }
591 }
592 // With every row gone, the last good range is kept rather than replaced
593 // by nothing: the panels freeze at a readable scale and the row labels
594 // say what happened, instead of the grid going blank.
595 if (!Number.isFinite(range.lo) || !Number.isFinite(range.hi)) continue;
ef3ae33Do not stretch the colormap across a constant field's roundoffJeremy Magland 596 // The floor is applied to what is drawn, not to what is tracked, so it
597 // never feeds back into the smoothing above.
598 const shown = floorRange(range.lo, range.hi);
599 this.#rangeBars[k]?.fill(shown.lo, shown.hi);
f467ffaCompare several solver settings side by side, on one clockJeremy Magland 600 for (const r of this.#rows) {
601 fillFieldValues(r.valueBufs[k], r.fields[k], this.#topo);
ef3ae33Do not stretch the colormap across a constant field's roundoffJeremy Magland 602 fillColors(r.colorBufs[k], r.valueBufs[k], shown.lo, shown.hi, cmap);
f467ffaCompare several solver settings side by side, on one clockJeremy Magland 603 r.scenes[k]?.updateColors(r.colorBufs[k]);
604 }
c90d0e2Check reference files in the browser's compare modeJeremy Magland 605 const f = this.#fileRow;
606 if (f) {
607 // Its values never change; only its coloring follows the shared range.
608 fillColors(f.colorBufs[k], f.valueBufs[k], shown.lo, shown.hi, cmap);
609 f.scenes[k]?.updateColors(f.colorBufs[k]);
610 }
613 this.#measureDifference();
614 this.#updateRowStats();
615 }
617 /**
618 * Relative L2 difference from the reference, per species, on the shared
619 * grid. Weighted by the Gauss weights, so it is the norm on the parameter
620 * sphere — not on the embedded surface, which would weight by the area
621 * element. That makes it a consistent diagnostic across variants rather than
622 * a physical quantity, which is all it is used for.
623 */
624 #measureDifference(): void {
c90d0e2Check reference files in the browser's compare modeJeremy Magland 625 // Against a reference file, every row is measured against its final state;
626 // otherwise against the chosen reference variant, whose own Δ is zero.
627 const ref = this.#fileRow ? null : this.#rows[this.#opts.reference];
628 const refFields = this.#fileRow?.fields ?? ref?.fields;
629 if (!refFields) return;
f467ffaCompare several solver settings side by side, on one clockJeremy Magland 630 const species = this.#opts.model.species;
631 for (const r of this.#rows) {
632 for (let k = 0; k < species.length; k++) {
633 if (r === ref) {
634 r.err[k] = 0;
635 continue;
636 }
637 const a = r.fields[k];
c90d0e2Check reference files in the browser's compare modeJeremy Magland 638 const b = refFields[k];
f467ffaCompare several solver settings side by side, on one clockJeremy Magland 639 if (!a || !b || a.length !== b.length) {
640 r.err[k] = NaN;
641 continue;
642 }
643 let num = 0;
644 let den = 0;
645 for (let i = 0; i < a.length; i++) {
646 const w = this.#weights[i];
647 const d = a[i] - b[i];
648 num += w * d * d;
649 den += w * b[i] * b[i];
650 }
651 r.err[k] = den > 0 ? Math.sqrt(num / den) : NaN;
652 }
653 }
654 }
656 /**
657 * Each row's standing line: how many of its own steps it took to reach the
658 * common time, and how far it is from the reference right now, per species.
659 * Per species rather than a single worst-case number because the two are
660 * genuinely different questions on a two-species model — the slow species is
661 * usually the one that has converged and the fast one the one that has not.
662 */
663 #updateRowStats(): void {
664 const species = this.#opts.model.species;
c90d0e2Check reference files in the browser's compare modeJeremy Magland 665 const ref = this.#fileRow ? null : this.#rows[this.#opts.reference];
f467ffaCompare several solver settings side by side, on one clockJeremy Magland 666 for (const r of this.#rows) {
667 const per = species
668 .map((s, k) => `${s} ${Number.isFinite(r.err[k]) ? r.err[k].toExponential(2) : '—'}`)
669 .join('<br>');
670 // Divergence is said, not implied. Scaled to a healthy row, a blown-up
671 // variant is a flat saturated panel, which on its own is easy to misread
672 // as a converged uniform state.
673 const body = !r.healthy
674 ? '<b class="cmp-diverged">diverged</b>'
675 : r === ref
676 ? '<b>reference</b>'
677 : `Δ ${per}`;
678 r.statEl.innerHTML = `${r.session.steps.toLocaleString()} steps<br>${body}`;
679 }
680 }
682 #status(): void {
c90d0e2Check reference files in the browser's compare modeJeremy Magland 683 const refFile = this.#opts.refFile;
684 const clock = refFile
685 ? `<b>t = ${this.#t.toFixed(2)} / ${(refFile.steps * CompareRun.baseDt(this.#opts.params)).toFixed(2)}</b>` +
686 (this.#finished
687 ? ` — <b>at the file's end time</b>: Δ is the final comparison against its final state`
688 : ` · Δ is the distance still to the file's <i>final</i> state — read it at the end time`)
689 : `<b>t = ${this.#t.toFixed(2)}</b> (same for every variant)`;
f467ffaCompare several solver settings side by side, on one clockJeremy Magland 690 this.#opts.onStatus(
f467ffaCompare several solver settings side by side, on one clockJeremy Magland 692 (this.#frameMs > 0 ? `${this.#frameMs.toFixed(1)} ms/frame · ` : '') +
693 this.#note,
694 );
695 }
697 #observeResize(): void {
c90d0e2Check reference files in the browser's compare modeJeremy Magland 698 const scenes = this.#allScenes();
f467ffaCompare several solver settings side by side, on one clockJeremy Magland 699 this.#resizeObs = new ResizeObserver(() => {
700 for (const s of scenes) {
701 const box = s.canvas.parentElement;
702 if (box) s.resize(box.clientWidth, box.clientHeight);
703 }
704 });
705 for (const s of scenes) {
706 const box = s.canvas.parentElement;
707 if (box) this.#resizeObs.observe(box);
708 }
709 }
711 // -------------------------------------------------------------- the clock
712 /**
713 * One frame advances every variant by the *same model time*: `frameSteps`
714 * base steps, which a ÷K variant covers in K times as many of its own. That
715 * is the whole reason dt varies by an integer divisor — the alternative is
716 * rounding each variant to the nearest step and comparing fields that are a
717 * fraction of a timestep apart, which would show up as a difference and be
718 * indistinguishable from a real one.
719 */
720 async #pump(): Promise<void> {
721 if (this.#pumping) return;
722 this.#pumping = true;
723 try {
724 while (this.#running && !this.#disposed) {
c90d0e2Check reference files in the browser's compare modeJeremy Magland 725 // Against a reference file the run is finite: the last frame takes
726 // however many base steps remain, so every variant lands exactly on
727 // the file's end time — where Δ against its final state is the
728 // comparison — and stops there rather than drifting past it.
729 const refFile = this.#opts.refFile;
730 const n = refFile
731 ? Math.min(this.#frameSteps, refFile.steps - this.#stepsDone)
732 : this.#frameSteps;
733 if (n <= 0) {
734 this.#running = false;
735 this.#opts.onFinished?.();
736 break;
737 }
f467ffaCompare several solver settings side by side, on one clockJeremy Magland 738 const t0 = performance.now();
c90d0e2Check reference files in the browser's compare modeJeremy Magland 739 for (const r of this.#rows) r.session.step(n * r.variant.dtDiv);
740 this.#stepsDone += n;
741 this.#t += n * CompareRun.baseDt(this.#opts.params);
f467ffaCompare several solver settings side by side, on one clockJeremy Magland 742 await this.draw();
743 if (this.#disposed) break;
744 const dt = performance.now() - t0;
745 this.#frameMs = this.#frameMs === 0 ? dt : this.#frameMs + 0.05 * (dt - this.#frameMs);
c90d0e2Check reference files in the browser's compare modeJeremy Magland 746 if (refFile && this.#stepsDone >= refFile.steps) {
747 this.#finished = true;
748 this.#running = false;
749 this.#status();
750 this.#opts.onFinished?.();
751 break;
752 }
754 await nextFrame();
755 }
756 if (!this.#disposed) {
757 await this.draw();
758 this.#status();
759 }
760 } finally {
761 this.#pumping = false;
762 }
763 }
764}
766const nextFrame = (): Promise<number> => new Promise(requestAnimationFrame);
c90d0e2Check reference files in the browser's compare modeJeremy Magland 768/** A whole spectral state re-indexed into a (wider) band's layout — the
769 * reference file's initial condition, in the form loadState takes. */
770function prolongState(
771 coeffs: Record<string, Float32Array>,
772 names: string[],
773 lmaxFrom: number,
774 lmaxTo: number,
775): Record<string, Float32Array> {
776 const out: Record<string, Float32Array> = {};
777 for (const name of names) out[name] = prolongCoeffs(coeffs[name], lmaxFrom, lmaxTo);
778 return out;
779}
f467ffaCompare several solver settings side by side, on one clockJeremy Magland 781/** Whether every entry is an ordinary number — false once a variant has left
782 * its convergence radius and saturated to infinity or NaN. */
783function allFinite(f: Float32Array | undefined): boolean {
784 if (!f) return false;
785 for (let i = 0; i < f.length; i++) if (!Number.isFinite(f[i])) return false;
786 return true;
787}
789type Bounds = { lo: number; hi: number };
791/** How far a field reaches from zero — the one number the rows are ranked by
792 * when deciding which of them sets a column's scale. */
793const peak = (b: Bounds): number => Math.max(Math.abs(b.lo), Math.abs(b.hi));
795/** Whichever of the given bounds reaches least far from zero; null if none. */
796function leastPeak(all: (Bounds | null)[]): Bounds | null {
797 let best: Bounds | null = null;
798 for (const b of all) {
799 if (b !== null && (best === null || peak(b) < peak(best))) best = b;
800 }
801 return best;
802}
804/** Min and max over the finite entries only; null when there are none. */
805function finiteRange(f: Float32Array | undefined): { lo: number; hi: number } | null {
806 if (!f) return null;
807 let lo = Infinity;
808 let hi = -Infinity;
809 for (let i = 0; i < f.length; i++) {
810 const v = f[i];
811 if (!Number.isFinite(v)) continue;
812 if (v < lo) lo = v;
813 if (v > hi) hi = v;
814 }
815 return lo <= hi ? { lo, hi } : null;
816}
818/**
819 * The DOM: a header row naming each species and carrying that column's shared
820 * color range, then one row per variant. The colorbar is per *column* rather
821 * than per panel because the range is shared — a bar on every panel would be
822 * the same bar repeated, and would suggest each panel had its own scaling,
823 * which is exactly the thing that would make the comparison a lie.
824 */
c90d0e2Check reference files in the browser's compare modeJeremy Magland 825/** The file row's label color — none of the variant palette, since it is not
826 * a variant: it is the thing they are all measured against. */
827const FILE_ROW_COLOR = '#57606a';
f467ffaCompare several solver settings side by side, on one clockJeremy Magland 829async function buildGrid(
830 opts: CompareOptions,
831 sessions: ModelSession[],
832 topo: SphereMeshTopology,
833 showDt: boolean,
835 rows: Row[];
836 fileRow: FileRow | null;
837 rangeBars: { fill: (lo: number, hi: number) => void }[];
838}> {
f467ffaCompare several solver settings side by side, on one clockJeremy Magland 839 const { container, model } = opts;
840 container.replaceChildren();
841 container.classList.add('compare');
843 const head = document.createElement('div');
844 head.className = 'cmp-row cmp-head';
845 const headSpacer = document.createElement('div');
846 headSpacer.className = 'cmp-rowlabel';
847 const headCols = document.createElement('div');
848 headCols.className = 'cmp-cols';
849 head.append(headSpacer, headCols);
850 container.append(head);
852 const rangeBars = model.species.map((name) => {
853 const col = document.createElement('div');
854 col.className = 'cmp-colhead';
855 const tag = document.createElement('b');
856 tag.textContent = name;
857 const canvas = document.createElement('canvas');
858 canvas.width = 160;
859 canvas.height = 8;
860 canvas.className = 'cmp-rangebar';
861 const lab = document.createElement('span');
862 lab.className = 'cmp-rangelab';
863 col.append(tag, canvas, lab);
864 headCols.append(col);
865 let painted = false;
866 return {
867 fill: (lo: number, hi: number): void => {
868 const ctx = canvas.getContext('2d');
869 if (ctx && !painted) {
870 painted = true;
871 const cmap = colormaps[opts.colormapName()] ?? colormaps.viridis;
872 for (let x = 0; x < canvas.width; x++) {
873 const [r, g, b] = cmap(x / (canvas.width - 1));
874 ctx.fillStyle = `rgb(${r},${g},${b})`;
875 ctx.fillRect(x, 0, 1, canvas.height);
876 }
877 }
878 lab.textContent = `${fmtValue(lo)} … ${fmtValue(hi)}`;
879 },
880 };
881 });
883 const sphereBg = getComputedStyle(document.documentElement)
884 .getPropertyValue('--sphere-bg')
885 .trim();
887 const rows: Row[] = [];
888 for (let i = 0; i < sessions.length; i++) {
889 const session = sessions[i];
890 const variant = opts.variants[i];
891 const color = VARIANT_COLORS[i % VARIANT_COLORS.length];
893 const coords = await session.renderPositions();
894 const posBuf = new Float32Array(topo.numVertices * 3);
895 fillPositions(posBuf, coords, topo, opts.morph);
897 const rowEl = document.createElement('div');
898 rowEl.className = 'cmp-row';
899 const labelEl = document.createElement('div');
900 labelEl.className = 'cmp-rowlabel';
901 labelEl.style.setProperty('--c', color);
902 const nameEl = document.createElement('div');
903 nameEl.className = 'cmp-rowname';
904 nameEl.textContent = variantLabel(variant, showDt);
905 const statEl = document.createElement('div');
906 statEl.className = 'cmp-rowstat';
907 labelEl.append(nameEl, statEl);
908 const colsEl = document.createElement('div');
909 colsEl.className = 'cmp-cols';
910 rowEl.append(labelEl, colsEl);
911 container.append(rowEl);
913 const scenes: SphereScene[] = [];
914 const valueBufs: Float32Array[] = [];
915 const colorBufs: Float32Array[] = [];
916 for (let k = 0; k < model.species.length; k++) {
917 const box = document.createElement('div');
918 box.className = 'sphere-box cmp-box';
919 colsEl.append(box);
920 const scene = new SphereScene(
921 box,
922 topo.numVertices,
923 topo.indices,
924 Float32Array.from(posBuf),
925 sphereBg || undefined,
926 );
927 scene.fitCamera();
928 scenes.push(scene);
929 valueBufs.push(new Float32Array(topo.numVertices));
930 colorBufs.push(new Float32Array(topo.numVertices * 3));
931 }
933 rows.push({
934 variant, session, color, coords, posBuf, scenes, valueBufs, colorBufs,
935 fields: [], err: model.species.map(() => 0), healthy: true, statEl,
936 });
937 }
c90d0e2Check reference files in the browser's compare modeJeremy Magland 939 // ---- the reference file's final state, as one more (static) row ---------
940 let fileRow: FileRow | null = null;
941 if (opts.refFile) {
942 const rf = opts.refFile;
943 // Synthesized through the coarsest session's display plan — exact, like
944 // every other use of the shared grid: the file's coefficients are
945 // band-limited at its lmax, which every variant's band contains.
946 const view = sessions[0].viewSht;
947 const lmaxTo = sessions[0].cfg.lmax;
948 const on = (q: Float32Array): Promise<Float32Array> =>
949 view.synth(prolongCoeffs(q, rf.lmax, lmaxTo));
950 const [gx, gy, gz] = [
951 await on(rf.geometryCoeffs.X),
952 await on(rf.geometryCoeffs.Y),
953 await on(rf.geometryCoeffs.Z),
954 ];
955 // The file's own surface, not a regeneration of it — interleaved xyz, the
956 // same layout renderPositions() hands back.
957 const coords = new Float32Array(3 * gx.length);
958 for (let i = 0; i < gx.length; i++) {
959 coords[3 * i] = gx[i];
960 coords[3 * i + 1] = gy[i];
961 coords[3 * i + 2] = gz[i];
962 }
963 const posBuf = new Float32Array(topo.numVertices * 3);
964 fillPositions(posBuf, coords, topo, opts.morph);
966 const rowEl = document.createElement('div');
967 rowEl.className = 'cmp-row';
968 const labelEl = document.createElement('div');
969 labelEl.className = 'cmp-rowlabel';
970 labelEl.style.setProperty('--c', FILE_ROW_COLOR);
971 const nameEl = document.createElement('div');
972 nameEl.className = 'cmp-rowname';
973 nameEl.textContent = 'reference file';
974 nameEl.title = rf.label;
975 const statEl = document.createElement('div');
976 statEl.className = 'cmp-rowstat';
977 statEl.innerHTML = `${rf.steps.toLocaleString()} steps<br><b>final state</b>`;
978 labelEl.append(nameEl, statEl);
979 const colsEl = document.createElement('div');
980 colsEl.className = 'cmp-cols';
981 rowEl.append(labelEl, colsEl);
982 container.append(rowEl);
984 const scenes: SphereScene[] = [];
985 const valueBufs: Float32Array[] = [];
986 const colorBufs: Float32Array[] = [];
987 const fields: Float32Array[] = [];
988 const bounds: (Bounds | null)[] = [];
989 for (let k = 0; k < model.species.length; k++) {
990 const box = document.createElement('div');
991 box.className = 'sphere-box cmp-box';
992 colsEl.append(box);
993 const scene = new SphereScene(
994 box,
995 topo.numVertices,
996 topo.indices,
997 Float32Array.from(posBuf),
998 sphereBg || undefined,
999 );
1000 scene.fitCamera();
1001 scenes.push(scene);
1002 const field = await on(rf.final[model.state[k]]);
1003 fields.push(field);
1004 bounds.push(finiteRange(field));
1005 const valueBuf = new Float32Array(topo.numVertices);
1006 fillFieldValues(valueBuf, field, topo);
1007 valueBufs.push(valueBuf);
1008 colorBufs.push(new Float32Array(topo.numVertices * 3));
1009 }
1010 fileRow = { coords, posBuf, scenes, valueBufs, colorBufs, fields, bounds };
1011 }
f467ffaCompare several solver settings side by side, on one clockJeremy Magland 1013 // Every panel shares one camera: the study is about the fields, and looking
1014 // at two of them from different angles is not comparing them.
c90d0e2Check reference files in the browser's compare modeJeremy Magland 1015 const all = [...rows.flatMap((r) => r.scenes), ...(fileRow?.scenes ?? [])];
f467ffaCompare several solver settings side by side, on one clockJeremy Magland 1016 for (let i = 1; i < all.length; i++) all[0].syncCamerasWith(all[i]);
c90d0e2Check reference files in the browser's compare modeJeremy Magland 1018 return { rows, fileRow, rangeBars };