0af3386Reaction-diffusion on spherical-harmonic surfacesJeremy Magland 1import { requestShtDevice, describeAdapter } from './sht/sht.ts';
2import { gridForLmax } from './sht/layout.ts';
3import { ModelSession } from './mgpu/session.ts';
4import { mModelByKey, presets, type MModel, type Params } from './mgpu/registry.ts';
5import { ModelCompileError, formatFailure } from './mgpu/errors.ts';
6import { EXTERNAL_OPS } from './mgpu/externals.ts';
7import { CodeEditor } from './editor/codeEditor.ts';
8import {
9 formatCommand,
10 resolvePreset,
11 DEFAULT_STEPS,
12 DEFAULT_WARMUP,
13 type RunSpec,
14} from './bench/runSpec.ts';
15import {
16 mGeometries,
17 mGeometryByKey,
18 defaultGeometryParams,
19 SPHERE_KEY,
20 DEFAULT_GEOMETRY_KEY,
21 type MGeometry,
22} from './geom/registry.ts';
23import {
24 buildTopology,
25 fillFieldValues,
26 fillPositions,
27 fillColors,
28 type SphereMeshTopology,
29} from './render/sphereMesh.ts';
30import { SphereScene } from './render/SphereScene.ts';
31import { Colorbar, fmtValue } from './render/colorbar.ts';
32import { colormaps, colormapNames } from './render/colormaps.ts';
33import { MovieRecorder } from './render/movie.ts';
35const $ = <T extends HTMLElement>(id: string): T =>
36 document.getElementById(id) as T;
38const elModel = $<HTMLSelectElement>('model');
39const elGeometry = $<HTMLSelectElement>('geometry');
40const elMorph = $<HTMLInputElement>('morph');
41const elNiter = $<HTMLSelectElement>('niter');
42const elLmax = $<HTMLSelectElement>('lmax');
43const elOversample = $<HTMLSelectElement>('oversample');
44const elColormap = $<HTMLSelectElement>('colormap');
45const elRunPause = $<HTMLButtonElement>('runpause');
46const elBenchmark = $<HTMLButtonElement>('benchmark');
47const elReseed = $<HTMLButtonElement>('reseed');
48const elResetView = $<HTMLButtonElement>('resetview');
49const elMovieToggle = $<HTMLButtonElement>('movietoggle');
50const elMovieBar = $('moviebar');
51const elMovieSpeed = $<HTMLSelectElement>('moviespeed');
52const elMovieRes = $<HTMLSelectElement>('movieres');
53const elMovieRotate = $<HTMLInputElement>('movierotate');
54const elMovie = $<HTMLButtonElement>('movie');
55const elParams = $('params');
56const elGeomParams = $('geomparams');
57const elGeomNote = $('geomnote');
58const elPanels = $('panels');
59const elStats = $('stats');
60const elBenchResult = $('benchresult');
61const elCmd = $('cmd');
62const elCopyCmd = $<HTMLButtonElement>('copycmd');
63const elBlurb = $('blurb');
64const elErr = $('err');
65const elSource = $<HTMLTextAreaElement>('source');
66const elHighlight = $('highlight');
67const elCompiled = $('compiled');
68const elEditorTitle = $('editor-title');
69const elEditorFile = $<HTMLSelectElement>('editor-file');
70const elRecompile = $<HTMLButtonElement>('recompile');
71const elRevert = $<HTMLButtonElement>('revert');
73for (const p of presets) {
74 const o = document.createElement('option');
75 o.value = p.key;
76 o.textContent = p.label;
77 elModel.append(o);
78}
79for (const g of mGeometries) {
80 const o = document.createElement('option');
81 o.value = g.key;
82 o.textContent = g.label;
83 elGeometry.append(o);
84}
85for (const [value, label] of [['model', 'the solver'], ['geometry', 'the surface']]) {
86 const o = document.createElement('option');
87 o.value = value;
88 o.textContent = label;
89 elEditorFile.append(o);
90}
91for (const name of colormapNames) {
92 const o = document.createElement('option');
93 o.value = name;
94 o.textContent = name;
95 elColormap.append(o);
96}
97elColormap.value = 'jet';
99/** Whichever .m is open: the solver or the surface. Both are MATLAB, compiled
100 * by the same backend, so one editor serves both. The host-provided operations
101 * are marked so the boundary between the file and what it is given is
102 * visible. */
103const editor = new CodeEditor({
104 textarea: elSource,
105 overlay: elHighlight,
106 external: EXTERNAL_OPS,
107 onInput: (value) => {
108 if (editing === 'geometry') editedGeomSource = value;
109 else editedSource = value;
110 elRecompile.textContent = 'Recompile *';
111 },
112});
114/** Timesteps submitted per rendered frame. Nothing is read back between them,
115 * so the batch costs one submit and one readback regardless of size. */
116const STEPS_PER_FRAME = 4;
118/**
119 * Steps in a solver-timing burst, and how often to run one.
120 *
121 * Timing the solver needs a `queue.onSubmittedWorkDone()` to know the work
122 * finished, and in a browser that is an IPC round trip into the GPU process — a
123 * fixed cost of a few milliseconds. Spread over one frame's four steps it would
124 * swamp them on a fast GPU and make the solver look far slower than it is. So the
125 * rate is measured in an occasional larger batch, where the single sync is
126 * amortized the way the desktop benchmark amortizes its own. The state is
127 * snapshotted and restored around the batch, so measuring never advances the
128 * simulation — otherwise the pattern would visibly lurch forward at every
129 * measurement.
130 */
131const MEASURE_BURST = 32;
132const MEASURE_EVERY_MS = 2000;
134/**
135 * 'auto' display oversampling targets this many render latitudes: the factor is
136 * the smallest power of two (up to 4) that reaches it. A solver grid already
137 * this fine gains nothing visually and is not oversampled.
138 */
139const AUTO_RENDER_NLAT = 256;
141/** The display oversampling factor the UI currently asks for. */
142function resolveOversample(): number {
143 if (elOversample.value !== 'auto') return Number(elOversample.value);
144 const { nlat } = gridForLmax(Number(elLmax.value), model.pdeg);
145 let os = 1;
146 while (os < 4 && os * nlat < AUTO_RENDER_NLAT) os *= 2;
147 return os;
148}
150/**
151 * Movie frame rate, and a cap on frames per movie. Playback speed comes from
152 * the UI, in simulation-time units per second of video; the movie's length is
153 * the run's t at that speed, and the frame count follows from it — capped by
154 * the run's own step count (a step is at most one frame) and by
155 * MOVIE_MAX_FRAMES to bound encode time and file size. Frame timestamps are
156 * derived from simulation time, so a capped movie keeps its duration and
157 * speed exactly, at a lower effective frame rate.
158 */
159const MOVIE_FPS = 30;
160const MOVIE_MAX_FRAMES = 3600;
162/** Movie auto-rotation: camera revolutions per second of video. Measured in
163 * video time, so the orbit pace on screen is the same at every export speed. */
164const MOVIE_ROTATE_RPS = 1 / 120;
166// ---------------------------------------------------------------- state
167let device: GPUDevice | null = null;
168let session: ModelSession | null = null;
169let topo: SphereMeshTopology | null = null;
170let scenes: SphereScene[] = [];
171let colorbars: Colorbar[] = [];
172let valueBufs: Float32Array[] = [];
173let colorBufs: Float32Array[] = [];
174let ranges: { lo: number; hi: number }[] = [];
175let resizeObs: ResizeObserver | null = null;
177const initial = resolvePreset(presets[0].key);
178let model: MModel = mModelByKey(initial.model.key)!;
179let params: Params = initial.params;
180let geometry: MGeometry = mGeometryByKey(DEFAULT_GEOMETRY_KEY)!;
181let geomParams: Params = defaultGeometryParams(geometry);
182/** Which file the editor is showing. */
183let editing: 'model' | 'geometry' = 'model';
184/** Each .m as edited in the page; `null` while it matches the file. */
185let editedSource: string | null = null;
186let editedGeomSource: string | null = null;
187/** Sphere (0) to surface (1). Display only; does not touch the solver. */
188let morph = 1;
189let seed = 1;
190let running = false;
191let adapterName = '';
192let pumping = false;
193let movieBusy = false;
194let movieCancel = false;
195let solverMs = 0;
196let frameMs = 0;
197let lastMeasure = 0;
198let generation = 0; // bumped on every rebuild to cancel stale pumps
199/** Surface coordinates on the render grid, interleaved xyz; null before the
200 * first build. Kept so the morph slider can re-fill positions without
201 * re-synthesizing. */
202let coords: Float32Array | null = null;
203let posBuf: Float32Array | null = null;
205const source = (): string => editedSource ?? model.source;
206const geomSource = (): string => editedGeomSource ?? geometry.source;
208// ---------------------------------------------------------------- UI wiring
209function buildParamInputs(): void {
210 elParams.replaceChildren();
211 for (const spec of model.params) {
212 const label = document.createElement('label');
213 label.textContent = `${spec.label} `;
214 const input = document.createElement('input');
215 input.type = 'number';
216 input.min = String(spec.min);
217 input.max = String(spec.max);
218 input.step = String(spec.step);
219 input.value = String(params[spec.key]);
220 input.addEventListener('change', () => {
221 const v = Number(input.value);
222 if (Number.isFinite(v)) params[spec.key] = v;
223 // Parameters are uniforms, not constants baked into the kernels, so a
224 // change costs an upload rather than a recompile.
225 session?.setParams(params);
226 updateCommand();
227 });
228 label.append(input);
229 elParams.append(label);
230 }
231}
233/**
234 * The shape's own parameters. Unlike the model's, these are NOT uniforms: the
235 * surface is evaluated once at build time and reduced to coefficients, so
236 * moving one rebuilds the geometry (and with it the mesh), though not the
237 * simulation's compiled step.
238 */
239function buildGeomParamInputs(): void {
240 elGeomParams.replaceChildren();
241 if (geometry.params.length === 0) return;
242 const tag = document.createElement('label');
243 tag.textContent = `${geometry.key}.m`;
244 elGeomParams.append(tag);
245 for (const spec of geometry.params) {
246 const label = document.createElement('label');
247 label.textContent = `${spec.label} `;
248 const input = document.createElement('input');
249 input.type = 'number';
250 input.min = String(spec.min);
251 input.max = String(spec.max);
252 input.step = String(spec.step);
253 input.value = String(geomParams[spec.key]);
254 input.addEventListener('change', () => {
255 const v = Number(input.value);
256 if (Number.isFinite(v)) geomParams[spec.key] = v;
257 viewChange = viewChange.then(() => applyGeometry());
258 });
259 label.append(input);
260 elGeomParams.append(label);
261 }
262}
264function applyPreset(presetKey: string): void {
265 const resolved = resolvePreset(presetKey);
266 const next = mModelByKey(resolved.model.key);
267 if (!next) {
268 elErr.textContent = `No .m model for '${resolved.model.key}'`;
269 return;
270 }
271 model = next;
272 params = resolved.params;
273 editedSource = null;
274 buildParamInputs();
275 elBlurb.textContent = model.blurb;
276 showEditorFile();
277 updateCommand();
278}
280function applyGeometryChoice(key: string): void {
281 const next = mGeometryByKey(key);
282 if (!next) {
283 elErr.textContent = `No .m geometry for '${key}'`;
284 return;
285 }
286 geometry = next;
287 geomParams = defaultGeometryParams(geometry);
288 editedGeomSource = null;
289 buildGeomParamInputs();
290 showEditorFile();
291}
293/** Load the chosen file into the editor, keeping any unsaved edit to it. */
294function showEditorFile(): void {
295 editing = elEditorFile.value === 'geometry' ? 'geometry' : 'model';
296 if (editing === 'geometry') {
297 editor.value = geomSource();
298 elEditorTitle.textContent = `geometries/${geometry.key}.m — shape(), compiled to WebGPU`;
299 } else {
300 editor.value = source();
301 elEditorTitle.textContent = `models/${model.key}.m — init() and step(), compiled to WebGPU`;
302 }
303}
305/** The run currently on screen, as the benchmark's RunSpec. */
306function currentSpec(): RunSpec {
307 return {
308 preset: elModel.value,
309 lmax: Number(elLmax.value),
310 seed,
311 steps: DEFAULT_STEPS,
312 warmup: DEFAULT_WARMUP,
313 params,
314 geometry: geometry.key,
315 geometryParams: geomParams,
316 niter: Number(elNiter.value),
317 };
318}
320function updateCommand(): void {
321 elCmd.textContent = formatCommand(currentSpec());
322}
324elModel.addEventListener('change', () => {
325 applyPreset(elModel.value);
326 void rebuild();
327});
328elLmax.addEventListener('change', () => void rebuild());
329// The solve iteration count is unrolled into the compiled step, so unlike a
330// parameter it cannot be changed without recompiling.
331elNiter.addEventListener('change', () => void rebuild());
332// Oversampling and geometry are display-or-data changes, not code ones, so
333// they swap things in place rather than rebuilding the run. Serialized through
334// one chain: a rapid second change waits its turn.
335let viewChange = Promise.resolve();
336elOversample.addEventListener('change', () => {
337 viewChange = viewChange.then(() => applyOversample());
338});
339elGeometry.addEventListener('change', () => {
340 applyGeometryChoice(elGeometry.value);
341 viewChange = viewChange.then(() => applyGeometry());
342});
343// Morph is pure rendering: no readback, no GPU work, just the vertex buffer.
344elMorph.addEventListener('input', () => {
345 morph = Number(elMorph.value);
346 applyMorph();
347});
348elColormap.addEventListener('change', () => void draw());
349elEditorFile.addEventListener('change', () => showEditorFile());
351function setRunning(next: boolean): void {
352 running = next;
353 elRunPause.textContent = running ? 'Pause' : 'Run';
354 if (running) void pump();
355}
357elRunPause.addEventListener('click', () => setRunning(!running));
358elBenchmark.addEventListener('click', () => void benchmark());
359elReseed.addEventListener('click', () => {
360 seed = (Math.random() * 2 ** 31) >>> 0;
361 setRunning(false);
362 updateCommand();
363 void reseed();
364});
365elResetView.addEventListener('click', () => {
366 for (const s of scenes) s.resetCamera();
367});
368elMovieToggle.addEventListener('click', () => {
369 elMovieBar.hidden = !elMovieBar.hidden;
370});
371elMovie.addEventListener('click', () => {
372 if (movieBusy) movieCancel = true;
373 else void recordMovie();
374});
376elRecompile.addEventListener('click', () => {
377 if (editing === 'geometry') editedGeomSource = editor.value;
378 else editedSource = editor.value;
379 void rebuild();
380});
381elRevert.addEventListener('click', () => {
382 if (editing === 'geometry') editedGeomSource = null;
383 else editedSource = null;
384 showEditorFile();
385 void rebuild();
386});
388// The command reproduces this run's parameters on the desktop; keep it
389// selectable even where the clipboard API is unavailable.
390elCopyCmd.addEventListener('click', () => {
391 const text = elCmd.textContent ?? '';
392 const flash = (msg: string): void => {
393 elCopyCmd.textContent = msg;
394 setTimeout(() => (elCopyCmd.textContent = 'Copy'), 1200);
395 };
396 const selectCommand = (): void => {
397 const range = document.createRange();
398 range.selectNodeContents(elCmd);
399 const sel = getSelection();
400 sel?.removeAllRanges();
401 sel?.addRange(range);
402 flash('Selected');
403 };
404 if (!navigator.clipboard) return selectCommand();
405 navigator.clipboard.writeText(text).then(() => flash('Copied'), selectCommand);
406});
408// ---------------------------------------------------------------- setup
409function disposeView(): void {
410 for (const s of scenes) s.dispose();
411 scenes = [];
412 colorbars = [];
413 topo = null;
414 coords = null;
415 posBuf = null;
416 resizeObs?.disconnect();
417 resizeObs = null;
418 elPanels.replaceChildren();
419}
421/**
422 * Build the mesh, scenes, colorbars and per-species buffers on the current
423 * render grid, from surface coordinates already synthesized there. Call
424 * disposeView() first. The color ranges are kept if present, so a display-only
425 * rebuild (an oversampling change) does not pop the shading; a full rebuild
426 * clears `ranges` beforehand.
427 */
428function buildView(surface: Float32Array): void {
429 if (!session) return;
430 const view = session.viewSht;
431 const { nphi } = view.cfg;
432 const phi = new Float64Array(nphi);
433 for (let j = 0; j < nphi; j++) phi[j] = (2 * Math.PI * j) / nphi;
434 topo = buildTopology(view.cosTheta, phi);
435 coords = surface;
436 posBuf = new Float32Array(topo.numVertices * 3);
437 fillPositions(posBuf, coords, topo, morph);
439 const sphereBg = getComputedStyle(document.documentElement)
440 .getPropertyValue('--sphere-bg')
441 .trim();
442 for (let k = 0; k < model.species.length; k++) {
443 const panel = document.createElement('div');
444 panel.className = 'panel';
445 const box = document.createElement('div');
446 box.className = 'sphere-box';
447 const tag = document.createElement('div');
448 tag.className = 'species-tag';
449 tag.textContent = model.species[k];
450 box.append(tag);
451 const side = document.createElement('div');
452 panel.append(box, side);
453 elPanels.append(panel);
455 const scene = new SphereScene(
456 box,
457 topo.numVertices,
458 topo.indices,
459 // Each scene owns its position buffer: three.js uploads from it, and the
460 // morph rewrites all of them from the one shared `coords`.
461 Float32Array.from(posBuf),
462 sphereBg || undefined,
463 );
464 scene.fitCamera();
465 scenes.push(scene);
466 colorbars.push(new Colorbar(side));
467 valueBufs[k] = new Float32Array(topo.numVertices);
468 colorBufs[k] = new Float32Array(topo.numVertices * 3);
469 if (!ranges[k]) ranges[k] = { lo: NaN, hi: NaN };
470 }
471 for (let k = 1; k < scenes.length; k++) scenes[0].syncCamerasWith(scenes[k]);
473 resizeObs = new ResizeObserver(() => {
474 const boxes = elPanels.querySelectorAll<HTMLElement>('.sphere-box');
475 boxes.forEach((box, i) => {
476 scenes[i]?.resize(box.clientWidth, box.clientHeight);
477 });
478 });
479 elPanels
480 .querySelectorAll<HTMLElement>('.sphere-box')
481 .forEach((box) => resizeObs!.observe(box));
482}
484/**
485 * Apply the UI's oversampling choice to the running session. Display-only: the
486 * session and its state survive; only the display plan, mesh and scenes are
487 * rebuilt, keeping the camera pose and color ranges. The pump is drained first
488 * so no readback is in flight on the plan being replaced.
489 */
490async function applyOversample(): Promise<void> {
491 if (!session) return;
492 const gen = generation;
493 const os = resolveOversample();
494 if (os === session.oversample) return;
495 const wasRunning = running;
496 setRunning(false);
497 while (pumping) await nextFrame();
498 if (gen !== generation || !session) return;
499 await session.setOversample(os);
500 if (gen !== generation || !session) return;
501 const surface = await session.renderPositions();
502 if (gen !== generation || !session) return;
503 const cam = scenes[0]?.cameraState();
504 disposeView();
505 buildView(surface);
506 if (cam) for (const s of scenes) s.setCameraState(cam);
507 await draw();
508 updateStats();
509 if (wasRunning) setRunning(true);
510}
512/**
513 * Re-evaluate the surface and swap it in. Data, not code: the compiled step is
514 * untouched and the simulation keeps its state and its model time, so a shape
515 * can be changed mid-run. Only the mesh is rebuilt.
516 */
517async function applyGeometry(): Promise<void> {
518 if (!session) return;
519 const gen = generation;
520 const wasRunning = running;
521 setRunning(false);
522 while (pumping) await nextFrame();
523 if (gen !== generation || !session) return;
524 try {
525 await session.setGeometry(geometry, geomParams, geomSource());
526 } catch (e) {
527 reportCompileError(e);
528 return;
529 }
530 if (gen !== generation || !session) return;
531 const surface = await session.renderPositions();
532 if (gen !== generation || !session) return;
533 const cam = scenes[0]?.cameraState();
534 disposeView();
535 buildView(surface);
536 if (cam) for (const s of scenes) s.setCameraState(cam);
537 elErr.textContent = '';
538 await draw();
539 updateGeomNote();
540 updateStats();
541 if (wasRunning) setRunning(true);
542}
544/** Re-place the vertices for the current morph. No GPU work and no readback —
545 * the surface is already on the CPU, so this is a buffer fill per panel. */
546function applyMorph(): void {
547 if (!topo || !coords || !posBuf) return;
548 fillPositions(posBuf, coords, topo, morph);
549 for (const s of scenes) s.updatePositions(posBuf);
550}
552/** What the surface is, and the standing caveat about where it is not. */
553function updateGeomNote(): void {
554 if (!session) {
555 elGeomNote.textContent = '';
556 return;
557 }
558 const { lo, hi } = session.geometry.radiusRange();
559 const isSphere = session.geometryModel.key === SPHERE_KEY;
560 elGeomNote.innerHTML =
561 `<b>${session.geometryModel.label}</b> — ${session.geometryModel.blurb} ` +
562 `Radius ${lo.toFixed(3)}–${hi.toFixed(3)}. ` +
563 (isSphere
564 ? 'This is the round-sphere case, so the solver is exact here.'
565 : '<b>Rendered only:</b> the Laplace–Beltrami operator in the .m is still ' +
566 'the round sphere\'s, so the pattern is the sphere\'s pattern painted ' +
567 'onto this shape.');
568}
570/** Report a compile failure, and select the offending text in the editor. */
571function reportCompileError(e: unknown): void {
572 elErr.textContent = formatFailure(e, source());
573 elCompiled.textContent = '';
574 if (e instanceof ModelCompileError && e.start !== undefined) {
575 editor.select(e.start, e.end ?? e.start);
576 }
577}
579async function rebuild(): Promise<void> {
580 generation++;
581 const gen = generation;
582 setRunning(false);
583 disposeView();
584 session?.destroy();
585 session = null;
586 solverMs = 0;
587 frameMs = 0;
588 lastMeasure = 0;
589 elErr.textContent = '';
590 updateCommand();
591 if (!device) return;
593 try {
594 session = await ModelSession.create({
595 device,
596 model,
597 params,
598 lmax: Number(elLmax.value),
599 source: source(),
600 oversample: resolveOversample(),
601 geometry,
602 geometryParams: geomParams,
603 geometrySource: geomSource(),
604 niter: Number(elNiter.value),
605 });
606 } catch (e) {
607 reportCompileError(e);
608 return;
609 }
610 if (gen !== generation) return;
612 session.seed(seed);
614 const plan = session.describe();
615 elCompiled.textContent =
616 `one step compiled to ${plan.step.length} GPU operations:\n` +
617 plan.step.map((l) => ` ${l}`).join('\n');
618 elRecompile.textContent = 'Recompile';
620 const surface = await session.renderPositions();
621 if (gen !== generation) return;
623 ranges = [];
624 buildView(surface);
626 await draw();
627 updateGeomNote();
628 updateStats();
629 void pump();
630}
632async function reseed(): Promise<void> {
633 if (!session) return;
634 const gen = generation;
635 session.seed(seed);
636 if (gen !== generation) return;
637 for (const r of ranges) {
638 r.lo = NaN;
639 r.hi = NaN;
640 }
641 await draw();
642 updateStats();
643}
645// ---------------------------------------------------------------- drawing
646async function draw(): Promise<void> {
647 if (!session || !topo) return;
648 const gen = generation;
649 const cmap = colormaps[elColormap.value] ?? colormaps.viridis;
650 for (let k = 0; k < model.species.length; k++) {
651 // The one readback per frame — the loop is otherwise entirely on the GPU.
652 // A rebuild can land while this is in flight and destroy the buffer being
653 // mapped, which rejects the map; that result is stale anyway, so drop it.
654 let field: Float32Array;
655 try {
656 field = await session.readSpecies(k);
657 } catch (e) {
658 if (gen !== generation) return;
659 throw e;
660 }
661 if (gen !== generation || !topo) return;
662 fillFieldValues(valueBufs[k], field, topo);
663 let lo = Infinity;
664 let hi = -Infinity;
665 for (const v of valueBufs[k]) {
666 if (v < lo) lo = v;
667 if (v > hi) hi = v;
668 }
669 // smooth the color range in both directions so the shading evolves
670 // gently as the pattern grows (out-of-range values clamp meanwhile)
671 const r = ranges[k];
672 if (!Number.isFinite(r.lo)) {
673 r.lo = lo;
674 r.hi = hi;
675 } else {
676 const a = 0.15;
677 r.lo += a * (lo - r.lo);
678 r.hi += a * (hi - r.hi);
679 }
680 if (r.hi - r.lo < 1e-9) {
681 const mid = (r.hi + r.lo) / 2;
682 r.lo = mid - 5e-10;
683 r.hi = mid + 5e-10;
684 }
685 fillColors(colorBufs[k], valueBufs[k], r.lo, r.hi, cmap);
686 scenes[k]?.updateColors(colorBufs[k]);
687 colorbars[k]?.update(cmap, r.lo, r.hi);
688 }
689}
691function updateStats(): void {
692 if (!session) return;
693 const { nlat, nphi } = session.cfg;
694 const kind = `WebGPU fp32${adapterName ? ` — ${adapterName}` : ''}`;
695 const solver =
696 solverMs > 0
697 ? `<b>${solverMs.toFixed(2)} ms/step</b> (${(1000 / solverMs).toFixed(0)} steps/s, ` +
698 `batch of ${MEASURE_BURST}, no readback)`
699 : '—';
700 const frame =
701 frameMs > 0
702 ? `${frameMs.toFixed(1)} ms/frame (${STEPS_PER_FRAME} steps + readback + render)`
703 : '—';
704 const view = session.viewSht.cfg;
705 const render =
706 session.oversample > 1
707 ? ` (display ${view.nlat}×${view.nphi}, ${session.oversample}×)`
708 : '';
709 elStats.innerHTML =
710 `<b>${kind}</b> · grid ${nlat}×${nphi}${render} · nlm ${session.sht.nlm.toLocaleString()} · ` +
711 `${session.sht.fourierMode.toUpperCase()} · ${session.geometryModel.key} · ` +
712 `${session.niter} solve iter${session.niter === 1 ? '' : 's'} · ` +
713 `solver ${solver} · ${frame} · ` +
714 `t = <b>${session.t.toFixed(2)}</b> (${session.steps} steps)`;
715}
717// ---------------------------------------------------------------- sim loop
718const nextFrame = () => new Promise<number>(requestAnimationFrame);
720async function pump(): Promise<void> {
721 if (pumping) return;
722 pumping = true;
723 const gen = generation;
724 try {
725 while (running && session && gen === generation) {
726 // Occasionally, a burst purely to measure the solver rate: many steps,
727 // one sync, nothing read back — directly comparable to the desktop
728 // benchmark's throughput number. State-preserving: the display and
729 // model time are unaffected.
730 if (performance.now() - lastMeasure > MEASURE_EVERY_MS) {
731 const ms = await session.measure(MEASURE_BURST);
732 if (gen !== generation) break;
733 solverMs = ms;
734 lastMeasure = performance.now();
735 }
737 // The frame itself. No explicit sync here — draw()'s readback already
738 // waits for the steps, so asking twice would only add a round trip.
739 const t0 = performance.now();
740 session.step(STEPS_PER_FRAME);
741 await draw();
742 if (gen !== generation) break;
743 frameMs = frameMs === 0
744 ? performance.now() - t0
745 : frameMs + 0.05 * (performance.now() - t0 - frameMs);
746 updateStats();
747 await nextFrame();
748 }
749 if (gen === generation) {
750 await draw();
751 updateStats();
752 }
753 } finally {
754 pumping = false;
755 }
756}
758/**
759 * Sustained solver benchmark, in the page.
760 *
761 * The same measurement `npm run bench` makes: batches of steps submitted
762 * together, waited for, never read back, with no rendering and no animation
763 * pacing in between. That makes it directly comparable to the terminal number,
764 * which is the only way to tell a genuinely slower browser GPU stack apart from
765 * the costs the app adds on top.
766 *
767 * It also reports the ramp — the first third of the run against the last. GPUs
768 * downclock when idle, and an animation-paced loop leaves them idle most of every
769 * frame, so a large ramp means the app's steady-state number is limited by clocks
770 * rather than by the work.
771 *
772 * These are ordinary steps: the simulation advances by them.
773 */
774async function benchmark(): Promise<void> {
775 if (!session || movieBusy) return;
776 setRunning(false);
777 const BATCH = 32;
778 const DURATION_MS = 2000;
779 elBenchResult.textContent = 'benchmarking…';
780 // A movie started mid-benchmark would replay while this loop still steps.
781 elMovie.disabled = true;
782 try {
783 await nextFrame();
785 const gen = generation;
786 const perStep: number[] = [];
787 const t0 = performance.now();
788 while (performance.now() - t0 < DURATION_MS) {
789 const b0 = performance.now();
790 session.step(BATCH);
791 await session.sync();
792 if (gen !== generation) return;
793 perStep.push((performance.now() - b0) / BATCH);
794 }
796 const mean = (xs: number[]): number => xs.reduce((a, b) => a + b, 0) / xs.length;
797 const all = mean(perStep);
798 const best = Math.min(...perStep);
799 const third = Math.max(1, Math.floor(perStep.length / 3));
800 const first = mean(perStep.slice(0, third));
801 const last = mean(perStep.slice(-third));
802 const steps = perStep.length * BATCH;
804 elBenchResult.innerHTML =
805 `sustained solver: <b>${all.toFixed(2)} ms/step</b> ` +
806 `(${(1000 / all).toFixed(0)} steps/s) · best ${best.toFixed(2)} · ` +
807 `ramp ${(first / last).toFixed(2)}× (${first.toFixed(2)} → ${last.toFixed(2)}) · ` +
808 `${steps} steps in batches of ${BATCH} · ` +
809 `compare with <code>npm run bench -- --lmax ${session.cfg.lmax}</code>`;
810 await draw();
811 updateStats();
812 } finally {
813 elMovie.disabled = false;
814 }
815}
817// ---------------------------------------------------------------- movie
818function saveBlob(blob: Blob, filename: string): void {
819 const url = URL.createObjectURL(blob);
820 const a = document.createElement('a');
821 a.href = url;
822 a.download = filename;
823 a.click();
824 setTimeout(() => URL.revokeObjectURL(url), 10_000);
825}
827/** Submit `n` steps in bounded command buffers — a single buffer encoding
828 * many thousands of steps can exhaust the encoder. */
829function submitSteps(n: number): void {
830 while (n > 0 && session) {
831 const chunk = Math.min(512, n);
832 session.step(chunk);
833 n -= chunk;
834 }
835}
837/** While recording, lock everything that could change the run mid-replay;
838 * the Movie button itself becomes the cancel button. */
839function setMovieUi(on: boolean): void {
840 const locked = [
841 elModel, elGeometry, elMorph, elNiter, elLmax, elOversample, elColormap,
842 elRunPause, elBenchmark, elReseed, elRecompile, elRevert, elEditorFile,
843 elMovieSpeed, elMovieRes, elMovieRotate, elMovieToggle,
844 ];
845 for (const el of locked) el.disabled = on;
846 elParams.querySelectorAll('input').forEach((input) => (input.disabled = on));
847 elGeomParams.querySelectorAll('input').forEach((input) => (input.disabled = on));
848 elMovie.textContent = on ? 'Cancel · 0%' : 'Export';
849}
851/**
852 * Recompute the run from t = 0 and download it as an MP4.
853 *
854 * The movie is not a recording of what already happened — it is the same
855 * trajectory recomputed: same seed, same source, and the *current* parameters
856 * and colormap throughout. Determinism makes this exact: after the replay the
857 * state is where it was, so the one session is reused and the app resumes as
858 * if nothing happened. Frames are composited from the live panels, so the
859 * movie shows the spheres at the current camera orientation — and the replay
860 * doubles as the progress display, since it is visible on screen.
861 */
862async function recordMovie(): Promise<void> {
863 if (!session || movieBusy) return;
864 if (session.steps === 0) {
865 elMovie.textContent = 'run first';
866 setTimeout(() => (elMovie.textContent = 'Export'), 1200);
867 return;
868 }
869 movieBusy = true;
870 movieCancel = false;
871 const gen = generation;
872 setMovieUi(true);
873 let wasRunning = false;
874 let total = 0;
875 let done = 0;
876 let seeded = false;
877 let camBefore: ReturnType<SphereScene['cameraState']> | undefined;
878 try {
879 // An in-flight display-grid swap replaces the scenes whose canvases the
880 // recorder captures, and resumes the run when it lands — let it finish.
881 await viewChange;
882 if (gen !== generation || !session) return;
883 wasRunning = running;
884 setRunning(false);
885 while (pumping) await nextFrame(); // let an in-flight live frame drain
886 if (gen !== generation || !session) return;
887 total = session.steps;
888 const speed = Number(elMovieSpeed.value) || 10;
889 const sphere = Number(elMovieRes.value) || 768;
890 const rotate = elMovieRotate.checked;
891 if (rotate) camBefore = scenes[0]?.cameraState();
892 // Render the scenes at exactly the chosen resolution for the recording —
893 // independent of the window size — and restore afterwards.
894 for (const s of scenes) s.captureSize(sphere);
895 const durationS = Math.max(session.t / speed, 2 / MOVIE_FPS);
896 const frames = Math.max(
897 2,
898 Math.min(Math.round(durationS * MOVIE_FPS) + 1, total + 1, MOVIE_MAX_FRAMES),
899 );
900 /** The step index captured as frame `i`; both endpoints land exactly. */
901 const stepAt = (i: number): number => Math.round((i * total) / (frames - 1));
903 const title =
904 (presets.find((p) => p.key === elModel.value)?.label ?? model.label) +
905 ` on ${geometry.label.toLowerCase()}` +
906 (editedSource !== null || editedGeomSource !== null ? ' (edited)' : '');
907 const subtitle = model.params
908 .map((spec) => `${spec.label} ${fmtValue(params[spec.key])}`)
909 .join(' · ');
910 const rec = await MovieRecorder.create({
911 panels: model.species.map((label, k) => ({ canvas: scenes[k].canvas, label })),
912 title,
913 subtitle,
914 speed,
915 fps: (frames - 1) / durationS,
916 sphere,
917 });
919 let finished = false;
920 try {
921 // Reset the color-range smoothing as a re-seed does, so the shading
922 // evolves in the movie the way it did live.
923 session.seed(seed);
924 seeded = true;
925 for (const r of ranges) {
926 r.lo = NaN;
927 r.hi = NaN;
928 }
929 const cmap = colormaps[elColormap.value] ?? colormaps.viridis;
930 let lastVideoS = 0;
931 for (let frame = 0; ; ) {
932 await draw();
933 if (gen !== generation) return;
934 if (movieCancel) break;
935 if (rotate) {
936 // Advance the orbit by this frame's share of video time; siblings
937 // follow scenes[0] through the usual camera sync.
938 const videoS = session.t / speed;
939 scenes[0]?.orbitBy(2 * Math.PI * MOVIE_ROTATE_RPS * (videoS - lastVideoS));
940 lastVideoS = videoS;
941 }
942 for (const s of scenes) s.renderNow();
943 await rec.addFrame(
944 session.t,
945 model.species.map((_, k) => ({ cmap, lo: ranges[k].lo, hi: ranges[k].hi })),
946 );
947 if (++frame >= frames) {
948 finished = true;
949 break;
950 }
951 const target = stepAt(frame);
952 submitSteps(target - done);
953 done = target;
954 elMovie.textContent = `Cancel · ${Math.round((100 * done) / total)}%`;
955 }
956 if (finished) {
957 const blob = await rec.finish();
958 saveBlob(
959 blob,
960 `turing-surface-${model.key}-${geometry.key}-` +
961 `t${session.t.toFixed(2)}-${speed}x.mp4`,
962 );
963 }
964 } finally {
965 if (!finished) rec.cancel();
966 }
967 } catch (e) {
968 elErr.textContent = `movie: ${e instanceof Error ? e.message : e}`;
969 } finally {
970 // A cancelled replay stopped short of where the run was; step the
971 // remainder — determinism makes this land exactly there.
972 if (seeded && gen === generation && session) {
973 while (done < total && gen === generation && session) {
974 const n = Math.min(4096, total - done);
975 submitSteps(n);
976 done += n;
977 elMovie.textContent = `restoring · ${Math.round((100 * done) / total)}%`;
978 await session.sync();
979 }
980 await draw();
981 updateStats();
982 }
983 if (gen === generation) {
984 for (const s of scenes) s.restoreSize();
985 }
986 if (camBefore && gen === generation) {
987 for (const s of scenes) s.setCameraState(camBefore);
988 }
989 movieBusy = false;
990 setMovieUi(false);
991 if (gen === generation) setRunning(wasRunning);
992 }
993}
995// ---------------------------------------------------------------- boot
996async function boot(): Promise<void> {
997 elModel.value = presets[0].key;
998 elGeometry.value = DEFAULT_GEOMETRY_KEY;
999 elMorph.value = String(morph);
1000 applyGeometryChoice(DEFAULT_GEOMETRY_KEY);
1001 applyPreset(presets[0].key);
1002 try {
1003 device = await requestShtDevice();
1004 adapterName = await describeAdapter(device);
1005 } catch (e) {
1006 device = null;
1007 elErr.textContent =
1008 `WebGPU is not available (${e instanceof Error ? e.message : e}). ` +
1009 `This demo compiles the MATLAB solver to WebGPU compute shaders, so it ` +
1010 `needs a WebGPU-capable browser (Chrome/Edge 113+).`;
1011 return;
1012 }
1013 device.lost.then((info) => {
1014 if (info.reason !== 'destroyed') {
1015 elErr.textContent = `WebGPU device lost: ${info.message}`;
1016 }
1017 });
1018 await rebuild();
1019}
1021void boot();