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