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';
beac00aMerge main into random-fieldsJeremy Magland 44import { 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';
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 */
60const RENDER_NLAT = 256;
61const RENDER_NLAT_CROWDED = 128;
62const CROWDED_PANELS = 6;
64/** See main.ts's DISPATCH_BUDGET — the same watchdog argument, per variant. */
65const DISPATCH_BUDGET = 1000;
66const STEPS_PER_FRAME_BASE = 4;
68export 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;
beac00aMerge main into random-fieldsJeremy Magland 82 /** Wavelength of the seeded random field, shared by every variant — one
83 * initial condition means one wavelength as much as one seed. */
84 lam3?: number;
86 colormapName: () => string;
87 /** Where the variant grid goes (the app's #panels). */
88 container: HTMLElement;
89 /** Progress and, afterwards, the standing description of the study. */
90 onStatus: (html: string) => void;
91}
93interface Row {
94 variant: Variant;
95 session: ModelSession;
96 color: string;
97 /** Surface coordinates on the shared render grid — this variant's own. */
98 coords: Float32Array;
99 posBuf: Float32Array;
100 scenes: SphereScene[];
101 valueBufs: Float32Array[];
102 colorBufs: Float32Array[];
103 /** Fields read this frame, one per species, on the shared grid. */
104 fields: Float32Array[];
105 /** Relative difference from the reference, one per species. */
106 err: number[];
107 /** False once any species has left the floating-point numbers — the shape a
108 * variant outside the convergence radius eventually fails in. Such a row is
109 * never used to scale a column, and its label says so. */
110 healthy: boolean;
111 statEl: HTMLElement;
112}
114export class CompareRun {
115 #opts: CompareOptions;
116 #rows: Row[] = [];
117 #topo: SphereMeshTopology;
118 /** Quadrature weight per grid point of the shared grid, for the L2 norm. */
119 #weights: Float64Array;
120 #rangeBars: { fill: (lo: number, hi: number) => void }[] = [];
121 /** Smoothed color range per species, shared by every variant so the panels
122 * in a column are directly comparable by eye and not just by number. */
123 #ranges: { lo: number; hi: number }[] = [];
124 #resizeObs: ResizeObserver | null = null;
126 #running = false;
127 #pumping = false;
128 #disposed = false;
129 #morph: number;
130 /** Base steps per frame; variant i takes this times its dtDiv. */
131 #frameSteps = STEPS_PER_FRAME_BASE;
132 /** Model time all variants are at — one number, by construction. */
133 #t = 0;
134 #frameMs = 0;
135 #note: string;
137 private constructor(init: {
138 opts: CompareOptions;
139 rows: Row[];
140 topo: SphereMeshTopology;
141 weights: Float64Array;
142 rangeBars: { fill: (lo: number, hi: number) => void }[];
143 frameSteps: number;
144 note: string;
145 }) {
146 this.#opts = init.opts;
147 this.#rows = init.rows;
148 this.#topo = init.topo;
149 this.#weights = init.weights;
150 this.#rangeBars = init.rangeBars;
151 this.#frameSteps = init.frameSteps;
152 this.#note = init.note;
153 this.#morph = init.opts.morph;
154 this.#ranges = init.opts.model.species.map(() => ({ lo: NaN, hi: NaN }));
155 }
157 get variants(): Variant[] {
158 return this.#rows.map((r) => r.variant);
159 }
161 /** The variant everything else is measured against — the one whose numbers
162 * stand on their own, so the one the app quotes when it has to quote one. */
163 get referenceSession(): ModelSession | null {
164 return this.#rows[this.#opts.reference]?.session ?? null;
165 }
167 get referenceIndex(): number {
168 return this.#opts.reference;
169 }
171 /** The base timestep a variant's dtDiv divides. */
172 static baseDt(params: Params): number {
173 return params.dt ?? 0;
174 }
176 static async create(opts: CompareOptions): Promise<CompareRun> {
177 const { device, model, variants } = opts;
178 const baseDt = CompareRun.baseDt(opts.params);
179 const showDt = variants.some((v) => v.dtDiv !== variants[0].dtDiv);
180 const sessions: ModelSession[] = [];
181 // Scenes own a WebGL context and an animation frame each, so a failure
182 // after the grid is up has to take them down explicitly — removing their
183 // canvases from the DOM would leave both running.
184 let built: Row[] = [];
186 try {
187 for (let i = 0; i < variants.length; i++) {
188 const v = variants[i];
189 opts.onStatus(
190 `compiling ${i + 1}/${variants.length} — ${variantLabel(v, showDt)} ` +
191 `(a solve iteration is ~15 kernels per species, and there is no ` +
192 `pipeline cache across sessions)`,
193 );
194 // Yield, so the status actually paints before the compile blocks.
195 await new Promise<number>(requestAnimationFrame);
196 sessions.push(
197 await ModelSession.create({
198 device,
199 model,
200 params: { ...opts.params, dt: baseDt / v.dtDiv },
201 lmax: v.lmax,
202 source: opts.source,
203 geometry: opts.geometry,
204 geometryParams: opts.geometryParams,
205 geometrySource: opts.geometrySource,
206 niter: v.niter,
209 );
210 }
212 // ---- the shared display grid ----------------------------------------
213 const maxLmax = Math.max(...variants.map((v) => v.lmax));
214 const panels = variants.length * model.species.length;
215 const target = panels > CROWDED_PANELS ? RENDER_NLAT_CROWDED : RENDER_NLAT;
216 // Never below what the finest band needs to be representable at all
217 // (ShtPlan requires nlat > lmax), whatever the panel count says.
218 const nlat = Math.max(target, 2 * Math.ceil((maxLmax + 2) / 2));
219 let nphi = 1;
220 while (nphi < Math.max(2 * nlat, 2 * maxLmax + 1)) nphi *= 2;
221 for (const s of sessions) await s.setDisplayGrid(nlat, nphi);
223 // ---- one initial condition, on every grid ---------------------------
224 opts.onStatus('seeding all variants from one band-limited perturbation…');
225 const noise = await sharedNoise(sessions, model.seedAmp, opts.seed);
beac00aMerge main into random-fieldsJeremy Magland 226 const modes = await sharedModes(sessions[opts.reference] ?? sessions[0], opts.seed);
227 // One at a time: a seed submits its whole mode sum in pieces, and there
228 // is nothing to gain from interleaving several variants' worth of it.
229 for (let i = 0; i < sessions.length; i++) await sessions[i].seedWith(noise[i], modes);
231 // ---- the mesh, shared; the surface, per variant ---------------------
232 const view = sessions[0].viewSht;
233 const phi = new Float64Array(nphi);
234 for (let j = 0; j < nphi; j++) phi[j] = (2 * Math.PI * j) / nphi;
235 const topo = buildTopology(view.cosTheta, phi);
236 // Gauss weights carry the sin(theta) of the area element; the constant
237 // 2*pi/nphi is common to every point and cancels in the relative norm.
238 const weights = new Float64Array(nlat * nphi);
239 for (let i = 0; i < nlat; i++) {
240 for (let j = 0; j < nphi; j++) weights[i * nphi + j] = view.gaussWeights[i];
241 }
243 // ---- how many steps a frame may submit ------------------------------
244 // Per variant: its own unrolled step size times its dtDiv, since a ÷K
245 // variant takes K times as many steps to reach the same time.
246 let frameSteps = STEPS_PER_FRAME_BASE;
247 const ops: number[] = [];
248 for (let i = 0; i < sessions.length; i++) {
249 const n = Math.max(1, sessions[i].describe().step.length);
250 ops.push(n);
251 frameSteps = Math.min(
252 frameSteps,
253 Math.max(1, Math.floor(DISPATCH_BUDGET / (n * variants[i].dtDiv))),
254 );
255 }
256 frameSteps = Math.max(1, frameSteps);
258 // ---- the grid of panels ---------------------------------------------
259 const { rows, rangeBars } = await buildGrid(opts, sessions, topo, showDt);
260 built = rows;
262 const solverGrid = sessions.map((s) => `${s.cfg.nlat}×${s.cfg.nphi}`);
263 const note =
264 `${variants.length} variants · display grid ${nlat}×${nphi}` +
265 (sessions.some((s) => s.cfg.nlat > nlat)
266 ? ` (below the finest solver grid ${solverGrid[solverGrid.length - 1]} — display only)`
267 : '') +
268 ` · ${frameSteps} base step${frameSteps === 1 ? '' : 's'}/frame` +
269 ` · ops/step ${ops.join(', ')}`;
271 const run = new CompareRun({
272 opts, rows, topo, weights, rangeBars, frameSteps, note,
273 });
274 await run.draw();
275 run.#observeResize();
276 run.#status();
277 return run;
278 } catch (e) {
279 for (const r of built) for (const s of r.scenes) s.dispose();
280 for (const s of sessions) s.destroy();
281 opts.container.replaceChildren();
282 opts.container.classList.remove('compare');
283 throw e;
284 }
285 }
287 // ------------------------------------------------------------------ state
288 setRunning(next: boolean): void {
289 this.#running = next;
290 if (next) void this.#pump();
291 }
293 get running(): boolean {
294 return this.#running;
295 }
297 /** Re-seed every variant from one new shared perturbation. */
298 async reseed(seed: number): Promise<void> {
299 const wasRunning = this.#running;
300 this.#running = false;
301 while (this.#pumping) await nextFrame();
302 if (this.#disposed) return;
beac00aMerge main into random-fieldsJeremy Magland 303 const sessions = this.#rows.map((r) => r.session);
304 const noise = await sharedNoise(sessions, this.#opts.model.seedAmp, seed);
305 const modes = await sharedModes(this.referenceSession ?? sessions[0], seed);
306 // Checked per variant, not once: a seed awaits its own submission, so a
307 // dispose can land between two of them and destroy the sessions left.
308 for (let i = 0; i < sessions.length; i++) {
309 if (this.#disposed) return;
310 await sessions[i].seedWith(noise[i], modes);
311 }
313 for (const r of this.#ranges) {
314 r.lo = NaN;
315 r.hi = NaN;
316 }
317 await this.draw();
318 if (!this.#disposed && wasRunning) this.setRunning(true);
319 }
beac00aMerge main into random-fieldsJeremy Magland 321 /** Wavelength of the seeded random field. One number for the study: every
322 * variant seeds from the same field, so they seed at the same wavelength. */
323 get lam3(): number {
324 return this.#rows[0]?.session.lam3 ?? 0;
325 }
327 /** Change it on every variant. Like the single run's, this only takes effect
328 * on the next reseed, which is where the field is drawn. */
329 setLam3(lambda: number): void {
330 this.#opts.lam3 = lambda;
331 for (const r of this.#rows) r.session.setLam3(lambda);
332 }
f467ffaCompare several solver settings side by side, on one clockJeremy Magland 334 /** Model parameters changed. Each variant keeps its own dt. */
335 setParams(params: Params): void {
336 this.#opts.params = params;
337 const baseDt = CompareRun.baseDt(params);
338 for (const r of this.#rows) {
339 r.session.setParams({ ...params, dt: baseDt / r.variant.dtDiv });
340 }
341 }
343 setMorph(morph: number): void {
344 this.#morph = morph;
345 for (const r of this.#rows) {
346 fillPositions(r.posBuf, r.coords, this.#topo, morph);
347 for (const s of r.scenes) s.updatePositions(r.posBuf);
348 }
349 }
351 resetView(): void {
352 for (const r of this.#rows) for (const s of r.scenes) s.resetCamera();
353 }
355 dispose(): void {
356 this.#disposed = true;
357 this.#running = false;
358 this.#resizeObs?.disconnect();
359 this.#resizeObs = null;
360 for (const r of this.#rows) {
361 for (const s of r.scenes) s.dispose();
362 r.session.destroy();
363 }
364 this.#rows = [];
365 this.#opts.container.replaceChildren();
366 this.#opts.container.classList.remove('compare');
367 }
369 // ----------------------------------------------------------------- drawing
370 /**
371 * One frame's readback: every variant's every species, on the shared grid.
372 * Read first, then color — the range is shared down a column, so no panel can
373 * be filled until the column's range is known.
374 */
375 async draw(): Promise<void> {
376 if (this.#disposed) return;
377 const species = this.#opts.model.species;
378 // Sessions are independent, so their readbacks can be in flight together;
379 // within one session they must not be (they share its staging buffers).
380 await Promise.all(
381 this.#rows.map(async (r) => {
382 for (let k = 0; k < species.length; k++) {
383 r.fields[k] = await r.session.readSpecies(k);
384 }
385 }),
386 );
387 if (this.#disposed) return;
389 const cmap = colormaps[this.#opts.colormapName()] ?? colormaps.viridis;
391 /**
392 * What scales a column is the whole question, and it has three wrong
393 * answers.
394 *
395 * Per panel is wrong: a range each rescales every variant to itself and
396 * hides exactly the difference the grid exists to show. The union over
397 * variants is wrong for the opposite reason: a variant outside the
398 * iteration's convergence radius runs away to 1e20 and then to NaN, and a
399 * union range rescales the *whole column* to it, flattening every panel to
400 * one colour — which reads as "they all blew up" when only one did.
401 *
402 * The reference alone is wrong too, less obviously, and it is the case that
403 * actually bites: outside the convergence radius *more* Richardson
404 * iterations diverge *faster*, so the row that goes first is usually the
405 * highest-niter one — which is the reference.
406 *
407 * So the column is scaled by whichever variant **reaches least far from
408 * zero** — the least-blown-up one. That is a comparison between the rows,
409 * not a threshold on any of them, and the distinction is the whole point:
410 * any "is this value too big?" test has a window in which a diverging field
411 * is still under the limit, and for as long as that window lasts it drags
412 * the scale and flattens the grid, until it finally trips and everything
413 * springs back. A comparison has no such window — a run-away only has to be
414 * *larger* than a healthy row to stop setting the scale, which it is from
415 * its first bad step, and it stays larger no matter how many other rows go
416 * with it. One healthy variant is enough to keep the grid readable.
417 *
418 * The cost is a slight bias: among healthy variants the scale comes from
419 * the one with the smallest peak, so the others clip by however much they
420 * exceed it. They are approximations of the same solution, so that is a
421 * fraction of a percent, and the alternative is a display that a single
422 * divergence can take away.
423 */
424 const bounds = this.#rows.map((r) => species.map((_, k) => finiteRange(r.fields[k])));
425 this.#rows.forEach((r, i) => {
426 // A row with any non-finite value is out of the running entirely: its
427 // finite entries are whatever survived, and no rank over them means much.
428 r.healthy = species.every((_, k) => allFinite(r.fields[k]) && bounds[i][k] !== null);
429 });
431 for (let k = 0; k < species.length; k++) {
432 const anchor = leastPeak(this.#rows.map((r, i) => (r.healthy ? bounds[i][k] : null)));
433 const range = this.#ranges[k];
434 if (anchor) {
435 if (!Number.isFinite(range.lo)) {
436 range.lo = anchor.lo;
437 range.hi = anchor.hi;
438 } else {
439 // Smooth in both directions so the shading evolves gently as the
440 // pattern grows, as the single-run view does.
441 const a = 0.15;
442 range.lo += a * (anchor.lo - range.lo);
443 range.hi += a * (anchor.hi - range.hi);
444 }
445 }
446 // With every row gone, the last good range is kept rather than replaced
447 // by nothing: the panels freeze at a readable scale and the row labels
448 // say what happened, instead of the grid going blank.
449 if (!Number.isFinite(range.lo) || !Number.isFinite(range.hi)) continue;
ef3ae33Do not stretch the colormap across a constant field's roundoffJeremy Magland 450 // The floor is applied to what is drawn, not to what is tracked, so it
451 // never feeds back into the smoothing above.
452 const shown = floorRange(range.lo, range.hi);
453 this.#rangeBars[k]?.fill(shown.lo, shown.hi);
f467ffaCompare several solver settings side by side, on one clockJeremy Magland 454 for (const r of this.#rows) {
455 fillFieldValues(r.valueBufs[k], r.fields[k], this.#topo);
ef3ae33Do not stretch the colormap across a constant field's roundoffJeremy Magland 456 fillColors(r.colorBufs[k], r.valueBufs[k], shown.lo, shown.hi, cmap);
f467ffaCompare several solver settings side by side, on one clockJeremy Magland 457 r.scenes[k]?.updateColors(r.colorBufs[k]);
458 }
459 }
461 this.#measureDifference();
462 this.#updateRowStats();
463 }
465 /**
466 * Relative L2 difference from the reference, per species, on the shared
467 * grid. Weighted by the Gauss weights, so it is the norm on the parameter
468 * sphere — not on the embedded surface, which would weight by the area
469 * element. That makes it a consistent diagnostic across variants rather than
470 * a physical quantity, which is all it is used for.
471 */
472 #measureDifference(): void {
473 const ref = this.#rows[this.#opts.reference];
474 if (!ref) return;
475 const species = this.#opts.model.species;
476 for (const r of this.#rows) {
477 for (let k = 0; k < species.length; k++) {
478 if (r === ref) {
479 r.err[k] = 0;
480 continue;
481 }
482 const a = r.fields[k];
483 const b = ref.fields[k];
484 if (!a || !b || a.length !== b.length) {
485 r.err[k] = NaN;
486 continue;
487 }
488 let num = 0;
489 let den = 0;
490 for (let i = 0; i < a.length; i++) {
491 const w = this.#weights[i];
492 const d = a[i] - b[i];
493 num += w * d * d;
494 den += w * b[i] * b[i];
495 }
496 r.err[k] = den > 0 ? Math.sqrt(num / den) : NaN;
497 }
498 }
499 }
501 /**
502 * Each row's standing line: how many of its own steps it took to reach the
503 * common time, and how far it is from the reference right now, per species.
504 * Per species rather than a single worst-case number because the two are
505 * genuinely different questions on a two-species model — the slow species is
506 * usually the one that has converged and the fast one the one that has not.
507 */
508 #updateRowStats(): void {
509 const species = this.#opts.model.species;
510 const ref = this.#rows[this.#opts.reference];
511 for (const r of this.#rows) {
512 const per = species
513 .map((s, k) => `${s} ${Number.isFinite(r.err[k]) ? r.err[k].toExponential(2) : '—'}`)
514 .join('<br>');
515 // Divergence is said, not implied. Scaled to a healthy row, a blown-up
516 // variant is a flat saturated panel, which on its own is easy to misread
517 // as a converged uniform state.
518 const body = !r.healthy
519 ? '<b class="cmp-diverged">diverged</b>'
520 : r === ref
521 ? '<b>reference</b>'
522 : `Δ ${per}`;
523 r.statEl.innerHTML = `${r.session.steps.toLocaleString()} steps<br>${body}`;
524 }
525 }
527 #status(): void {
528 this.#opts.onStatus(
529 `<b>t = ${this.#t.toFixed(2)}</b> (same for every variant) · ` +
530 (this.#frameMs > 0 ? `${this.#frameMs.toFixed(1)} ms/frame · ` : '') +
531 this.#note,
532 );
533 }
535 #observeResize(): void {
536 const scenes = this.#rows.flatMap((r) => r.scenes);
537 this.#resizeObs = new ResizeObserver(() => {
538 for (const s of scenes) {
539 const box = s.canvas.parentElement;
540 if (box) s.resize(box.clientWidth, box.clientHeight);
541 }
542 });
543 for (const s of scenes) {
544 const box = s.canvas.parentElement;
545 if (box) this.#resizeObs.observe(box);
546 }
547 }
549 // -------------------------------------------------------------- the clock
550 /**
551 * One frame advances every variant by the *same model time*: `frameSteps`
552 * base steps, which a ÷K variant covers in K times as many of its own. That
553 * is the whole reason dt varies by an integer divisor — the alternative is
554 * rounding each variant to the nearest step and comparing fields that are a
555 * fraction of a timestep apart, which would show up as a difference and be
556 * indistinguishable from a real one.
557 */
558 async #pump(): Promise<void> {
559 if (this.#pumping) return;
560 this.#pumping = true;
561 try {
562 while (this.#running && !this.#disposed) {
563 const t0 = performance.now();
564 for (const r of this.#rows) r.session.step(this.#frameSteps * r.variant.dtDiv);
565 this.#t += this.#frameSteps * CompareRun.baseDt(this.#opts.params);
566 await this.draw();
567 if (this.#disposed) break;
568 const dt = performance.now() - t0;
569 this.#frameMs = this.#frameMs === 0 ? dt : this.#frameMs + 0.05 * (dt - this.#frameMs);
570 this.#status();
571 await nextFrame();
572 }
573 if (!this.#disposed) {
574 await this.draw();
575 this.#status();
576 }
577 } finally {
578 this.#pumping = false;
579 }
580 }
581}
583const nextFrame = (): Promise<number> => new Promise(requestAnimationFrame);
585/** Whether every entry is an ordinary number — false once a variant has left
586 * its convergence radius and saturated to infinity or NaN. */
587function allFinite(f: Float32Array | undefined): boolean {
588 if (!f) return false;
589 for (let i = 0; i < f.length; i++) if (!Number.isFinite(f[i])) return false;
590 return true;
591}
593type Bounds = { lo: number; hi: number };
595/** How far a field reaches from zero — the one number the rows are ranked by
596 * when deciding which of them sets a column's scale. */
597const peak = (b: Bounds): number => Math.max(Math.abs(b.lo), Math.abs(b.hi));
599/** Whichever of the given bounds reaches least far from zero; null if none. */
600function leastPeak(all: (Bounds | null)[]): Bounds | null {
601 let best: Bounds | null = null;
602 for (const b of all) {
603 if (b !== null && (best === null || peak(b) < peak(best))) best = b;
604 }
605 return best;
606}
608/** Min and max over the finite entries only; null when there are none. */
609function finiteRange(f: Float32Array | undefined): { lo: number; hi: number } | null {
610 if (!f) return null;
611 let lo = Infinity;
612 let hi = -Infinity;
613 for (let i = 0; i < f.length; i++) {
614 const v = f[i];
615 if (!Number.isFinite(v)) continue;
616 if (v < lo) lo = v;
617 if (v > hi) hi = v;
618 }
619 return lo <= hi ? { lo, hi } : null;
620}
622/**
623 * The DOM: a header row naming each species and carrying that column's shared
624 * color range, then one row per variant. The colorbar is per *column* rather
625 * than per panel because the range is shared — a bar on every panel would be
626 * the same bar repeated, and would suggest each panel had its own scaling,
627 * which is exactly the thing that would make the comparison a lie.
628 */
629async function buildGrid(
630 opts: CompareOptions,
631 sessions: ModelSession[],
632 topo: SphereMeshTopology,
633 showDt: boolean,
634): Promise<{ rows: Row[]; rangeBars: { fill: (lo: number, hi: number) => void }[] }> {
635 const { container, model } = opts;
636 container.replaceChildren();
637 container.classList.add('compare');
639 const head = document.createElement('div');
640 head.className = 'cmp-row cmp-head';
641 const headSpacer = document.createElement('div');
642 headSpacer.className = 'cmp-rowlabel';
643 const headCols = document.createElement('div');
644 headCols.className = 'cmp-cols';
645 head.append(headSpacer, headCols);
646 container.append(head);
648 const rangeBars = model.species.map((name) => {
649 const col = document.createElement('div');
650 col.className = 'cmp-colhead';
651 const tag = document.createElement('b');
652 tag.textContent = name;
653 const canvas = document.createElement('canvas');
654 canvas.width = 160;
655 canvas.height = 8;
656 canvas.className = 'cmp-rangebar';
657 const lab = document.createElement('span');
658 lab.className = 'cmp-rangelab';
659 col.append(tag, canvas, lab);
660 headCols.append(col);
661 let painted = false;
662 return {
663 fill: (lo: number, hi: number): void => {
664 const ctx = canvas.getContext('2d');
665 if (ctx && !painted) {
666 painted = true;
667 const cmap = colormaps[opts.colormapName()] ?? colormaps.viridis;
668 for (let x = 0; x < canvas.width; x++) {
669 const [r, g, b] = cmap(x / (canvas.width - 1));
670 ctx.fillStyle = `rgb(${r},${g},${b})`;
671 ctx.fillRect(x, 0, 1, canvas.height);
672 }
673 }
674 lab.textContent = `${fmtValue(lo)} … ${fmtValue(hi)}`;
675 },
676 };
677 });
679 const sphereBg = getComputedStyle(document.documentElement)
680 .getPropertyValue('--sphere-bg')
681 .trim();
683 const rows: Row[] = [];
684 for (let i = 0; i < sessions.length; i++) {
685 const session = sessions[i];
686 const variant = opts.variants[i];
687 const color = VARIANT_COLORS[i % VARIANT_COLORS.length];
689 const coords = await session.renderPositions();
690 const posBuf = new Float32Array(topo.numVertices * 3);
691 fillPositions(posBuf, coords, topo, opts.morph);
693 const rowEl = document.createElement('div');
694 rowEl.className = 'cmp-row';
695 const labelEl = document.createElement('div');
696 labelEl.className = 'cmp-rowlabel';
697 labelEl.style.setProperty('--c', color);
698 const nameEl = document.createElement('div');
699 nameEl.className = 'cmp-rowname';
700 nameEl.textContent = variantLabel(variant, showDt);
701 const statEl = document.createElement('div');
702 statEl.className = 'cmp-rowstat';
703 labelEl.append(nameEl, statEl);
704 const colsEl = document.createElement('div');
705 colsEl.className = 'cmp-cols';
706 rowEl.append(labelEl, colsEl);
707 container.append(rowEl);
709 const scenes: SphereScene[] = [];
710 const valueBufs: Float32Array[] = [];
711 const colorBufs: Float32Array[] = [];
712 for (let k = 0; k < model.species.length; k++) {
713 const box = document.createElement('div');
714 box.className = 'sphere-box cmp-box';
715 colsEl.append(box);
716 const scene = new SphereScene(
717 box,
718 topo.numVertices,
719 topo.indices,
720 Float32Array.from(posBuf),
721 sphereBg || undefined,
722 );
723 scene.fitCamera();
724 scenes.push(scene);
725 valueBufs.push(new Float32Array(topo.numVertices));
726 colorBufs.push(new Float32Array(topo.numVertices * 3));
727 }
729 rows.push({
730 variant, session, color, coords, posBuf, scenes, valueBufs, colorBufs,
731 fields: [], err: model.species.map(() => 0), healthy: true, statEl,
732 });
733 }
735 // Every panel shares one camera: the study is about the fields, and looking
736 // at two of them from different angles is not comparing them.
737 const all = rows.flatMap((r) => r.scenes);
738 for (let i = 1; i < all.length; i++) all[0].syncCamerasWith(all[i]);
740 return { rows, rangeBars };
741}