Comparing changes
main is 13 commits ahead of reference_soln.
Create pull request
export standalone matlab
Jeremy Magland committed
9b71d5bMerge pull request #10 from concept-collection/compare-interface
Jeremy Magland committed
71b6492Add short descriptions for each mode
Owen Melia committed
2bb4f78Adding a way to reset simulation from same random IC.
Owen Melia committed
ca37955Fixing some bugs in the UI
Owen Melia committed
c8e230aWIP: First draft at new interface with multiple comparison modes
Owen Melia committed
9389d73Update ids in index.html for easier refactor
Owen Melia committed
a158c73Merge pull request #9 from concept-collection/reference-compare-mode
Jeremy Magland committed
c88cad7Generate numbl's stdlib bundle in the preview workflow too
Jeremy Magland committed
1a01a2fGenerate numbl's stdlib bundle in the preview workflow too
Jeremy Magland committed
90d9678Open the reference comparison in one click
Jeremy Magland committed
3210e7dCheck reference files in the browser's compare mode
Jeremy Magland committed
c90d0e2Merge pull request #6 from concept-collection/reference_soln
Jeremy Magland committed
327b4ec18 changed files+2137−195
.github/workflows/preview.ymlmodified+6−1View file
@@ -93,7 +93,11 @@ jobs:
9393 # numbl is a `file:../../numbl` dependency: we use its compiler internals
9494 # (parser, lowerer, IR, inline pass), which its published package
9595 # `exports` do not expose. Clone it where that relative path expects it,
96- # pinned to the same ref as ci.yml and deploy.yml.
96+ # pinned to the same ref as ci.yml and deploy.yml. One file in it is
97+ # generated rather than committed — the interpreter's stdlib bundle,
98+ # which numbl gitignores — so a bare checkout is missing it and
99+ # `executeCode.ts` fails to resolve it; its generator only reads .m files
100+ # off disk, so plain `node` runs it without installing anything.
97101 - name: Check out numbl (sibling dependency)
98102 env:
99103 NUMBL_REF: main
@@ -101,6 +105,7 @@ jobs:
101105 git clone --filter=blob:none --no-checkout \
102106 https://github.com/flatironinstitute/numbl.git "$GITHUB_WORKSPACE/../../numbl"
103107 git -C "$GITHUB_WORKSPACE/../../numbl" checkout --quiet "$NUMBL_REF"
108+ node "$GITHUB_WORKSPACE/../../numbl/scripts/bundle-stdlib.ts"
104109 # --ignore-scripts: npm runs a linked package's `prepare` script, and
105110 # numbl's is husky, which is not installed here.
106111 - run: npm ci --ignore-scripts
.gitignoremodified+1−0View file
@@ -1,3 +1,4 @@
1+tmp/
12 node_modules/
23 dist/
34 *.log
README.mdmodified+20−0View file
@@ -497,6 +497,12 @@ package alone. Its binaries need glibc 2.29+. Other flags: `--steps`,
497497 `--warmup`, `--batch`, `--json`, `--help`; `DAWN_FLAGS='backend=vulkan'`
498498 (`;`-separated) passes Dawn options through.
499499
500+### The same run in MATLAB
501+
502+A run in the page needs a browser and a GPU; further analysis usually wants neither. The app therefore exports the run on screen as one self-contained MATLAB function file: **The same run as a standalone MATLAB script**, under the benchmark command, shows the script for copying and downloads it as `turing_surface_run.m`. The current model and geometry `.m` go in verbatim, edits in the page included, with the parameter values baked in; around them the file carries double-precision ports of everything the host provides: the transforms and their derivative shuffles, the metric weights, the seeded random field, and the run loop ([`src/export/`](src/export/)). It needs base MATLAB only, R2020b or newer, no toolboxes.
503+
504+Two deliberate differences from the page are stated in the script's own header: it runs in f64 where the GPU path is f32, and random draws use MATLAB's own `rng`, so a seed value picks a different member of the same random ensemble than the same value in the app. The script plots the pattern live and writes its initial and final spectral state to HDF5 in the reference-run layout of [docs/ellipsoid-reference-spec.md](docs/ellipsoid-reference-spec.md), so a MATLAB run can be loaded back into the page (**Compare against uploaded data**) or checked with `npm run ref -- --in turing_surface_run.h5`. Exported at the defaults, a 60-step Schnakenberg run on the ellipsoid replayed that way agrees with the app to relative L2 of about 1e-7, which is fp32 accumulation; the exported flux-form and Algorithm-4 models track each other to about 3e-10 in f64.
505+
500506 ## Tests
501507
502508 There is no second implementation of the solver to diff against, so the `.m`
@@ -609,6 +615,20 @@ the surface-correction iteration count independent of the file, and
609615 `--tolerance`/`--tolerance-linf` each independently turn their metric into a
610616 pass/fail for CI.
611617
618+The same check runs in the page: **Compare to reference…** picks a `.h5` and
619+opens the comparison in one step — the file's own settings (its recorded
620+niter, its band, its dt) as the single variant, paused at the file's exact
621+initial state, ready to Run. The file defines the whole problem — model,
622+parameters, geometry, initial state — and the run stops at the file's end
623+time, measured against one extra static row showing its final state on its
624+own surface. Watching *where* a variant leaves the reference (rather than
625+just reading one number per run) is the point. To widen the study, stop
626+comparing, pick more chips, and press Compare — the file stays loaded, with
627+the lmax choices floored at its band, since a narrower one could not hold
628+its initial state. Reading the file uses
629+[h5wasm](https://github.com/usnistgov/h5wasm)'s wasm build, loaded lazily on
630+the first file opened.
631+
612632 ## Development
613633
614634 ```
docs/ellipsoid-reference-spec.mdmodified+14−0View file
@@ -107,6 +107,20 @@ how much that correction term actually matters for a given run. `--tolerance
107107 <n>` and `--tolerance-linf <n>` each independently turn their metric into a
108108 pass/fail (nonzero exit code on failure), for use in CI.
109109
110+The browser demo runs the same check visually: **Compare to reference…**
111+picks a reference file and opens the comparison in one step, with the file's
112+own settings as the single variant, paused at its exact initial state. Run
113+takes it to the file's end time and stops; its final state shows as one
114+extra static row — on the file's own surface, with each variant's
115+relative-L2 distance to it updating live — and more variants can be added
116+from the compare bar's chips. Both readers share one parser
117+(`src/compare/referenceCase.ts`), so the layout above is interpreted
118+identically on the CLI and in the page.
119+
120+Note the files record only the two endpoint states (`initial/`, `final/`) —
121+no intermediate snapshots — so the comparison is meaningful at the end time;
122+the live Δ before that reads as "distance still to the final state".
123+
110124 ## Caveat
111125
112126 This repo runs fp32 on GPU; expect ~1e-4–1e-6 relative floating-point noise
index.htmlmodified+125−60View file
@@ -51,6 +51,19 @@
5151 }
5252 /* display:flex above would otherwise override the UA's [hidden] rule */
5353 .controls[hidden] { display: none; }
54+ .modes {
55+ display: flex; flex-wrap: wrap; gap: 8px;
56+ padding: 4px 0 10px; margin-bottom: 4px;
57+ border-bottom: 1px solid var(--line);
58+ }
59+ .modes .chip { font-size: 13px; padding: 4px 12px; }
60+ .mode-desc {
61+ color: var(--ink); margin: 0 0 12px; font-size: 13.5px;
62+ padding: 10px 14px; border-radius: 6px;
63+ border-left: 3px solid var(--accent);
64+ background: color-mix(in srgb, var(--accent) 10%, var(--bg));
65+ }
66+ .mode-desc:empty { display: none; }
5467 .controls label { color: var(--ink-2); font-size: 13px; white-space: nowrap; }
5568 select, input[type="number"], button {
5669 font: inherit; font-size: 13px;
@@ -104,6 +117,14 @@
104117 font: 12px/1.6 ui-monospace, SFMono-Regular, Menlo, Consolas, monospace;
105118 color: var(--ink); user-select: all;
106119 }
120+ details.cli > summary { cursor: pointer; }
121+ details.cli > summary::before { content: '▸'; font-size: 10px; color: var(--ink-2); }
122+ details.cli[open] > summary::before { content: '▾'; }
123+ #matlabscript {
124+ margin: 0; padding: 8px 10px; max-height: 30em; overflow: auto;
125+ font: 12px/1.5 ui-monospace, SFMono-Regular, Menlo, Consolas, monospace;
126+ color: var(--ink); white-space: pre; user-select: text;
127+ }
107128 #blurb { margin-top: 4px; font-size: 13px; color: var(--ink-2); }
108129 #err { color: #b35900; white-space: pre-wrap; font-size: 13px; }
109130 .editor {
@@ -116,10 +137,14 @@
116137 background: var(--sphere-bg); border-bottom: 1px solid var(--line);
117138 }
118139 .editor-head button { padding: 2px 10px; font-size: 12px; }
119- /* Source and compiled-op list side by side, so the editor gets the
120- height rather than sharing it with the list below. */
121- .editor-body { display: flex; align-items: stretch; }
122- .editor-code { position: relative; flex: 1 1 62%; min-width: 0; height: 34em; }
140+ /* Source and compiled-op list side by side. The fixed height lives on
141+ the row itself, not on either child: a child's own `height` would
142+ only be a *hint* stretch fills when unset, and `#compiled` sets its
143+ own smaller font, so the same em value would resolve to a shorter
144+ box than .editor-code's. Sizing the row instead makes both children
145+ stretch to one shared pixel height regardless of either's font. */
146+ .editor-body { display: flex; align-items: stretch; height: 34em; }
147+ .editor-code { position: relative; flex: 1 1 62%; min-width: 0; }
123148 /* The overlay and the textarea must agree on every metric that affects
124149 where a character lands. Keep these two rules together. */
125150 .editor-code > pre,
@@ -150,7 +175,9 @@
150175 color: var(--ink-2); white-space: pre;
151176 }
152177 @media (max-width: 860px) {
153- .editor-body { flex-direction: column; }
178+ /* Stacked, so each panel goes back to managing its own height rather
179+ than sharing the row's fixed one. */
180+ .editor-body { flex-direction: column; height: auto; }
154181 .editor-code { flex: none; height: 26em; }
155182 #compiled { border-left: 0; border-top: 1px solid var(--line); max-height: 12em; }
156183 }
@@ -229,16 +256,54 @@
229256 browser by <a href="https://numbl.org">numbl</a>. Edit either and watch
230257 it change. Drag to rotate.
231258 </p>
232- <div class="controls">
233- <label>preset
259+ <div class="modes" id="modebar">
260+ <button type="button" class="chip" id="mode-simulate" aria-pressed="true">Simulate</button>
261+ <button type="button" class="chip" id="mode-effort"
262+ title="Run several solver settings side by side on one clock">Compare computational effort</button>
263+ <button type="button" class="chip" id="mode-vs-sphere" disabled title="Coming soon">Compare against sphere (Coming soon)</button>
264+ <button type="button" class="chip" id="mode-vs-upload"
265+ title="Check this solver against a saved reference run">
266+ Compare against uploaded data</button>
267+ <input type="file" id="cmp-file" accept=".h5" hidden>
268+ </div>
269+ <p class="mode-desc" id="mode-desc"></p>
270+ <div class="controls ctrl-group" data-group="surface">
271+ <label>reaction-diffusion model
234272 <select id="model"></select>
235273 </label>
236274 <label title="The surface the pattern is solved on. Swapping it does not recompile the solver or restart the run.">geometry
237275 <select id="geometry"></select>
238276 </label>
239- <label title="Blend between the sphere (0) and the surface (1). Display only.">morph
240- <input type="range" id="morph" min="0" max="1" step="0.01" value="1" />
277+ </div>
278+ <div class="ctrl-group" data-group="surface-params">
279+ <div class="controls" id="params"></div>
280+ <div class="controls" id="geomparams"></div>
281+ </div>
282+ <div class="controls" id="comparebar" hidden>
283+ <div class="cmp-axes">
284+ <div class="cmp-axis">
285+ <span title="Iterations of the implicit diffusion solve">solve iters</span>
286+ <span id="cmp-niter" class="chips"></span>
287+ </div>
288+ <div class="cmp-axis">
289+ <span>lmax</span>
290+ <span id="cmp-lmax" class="chips"></span>
291+ </div>
292+ <div class="cmp-axis">
293+ <span title="Timestep, as a divisor of the model's dt. Divisors keep every variant on the same clock exactly.">dt</span>
294+ <span id="cmp-dt" class="chips"></span>
295+ </div>
296+ </div>
297+ <label title="The run everything else is measured against">reference
298+ <select id="cmp-ref"></select>
241299 </label>
300+ <span id="cmp-fileinfo" class="stats" hidden></span>
301+ <button id="cmp-fileclear" hidden
302+ title="Drop the reference file and compare the variants against each other again">×</button>
303+ <button id="cmp-start" class="primary">Compile comparison</button>
304+ <span id="cmp-count" class="stats"></span>
305+ </div>
306+ <div class="controls ctrl-group" data-group="solver">
242307 <label title="Iterations of the implicit diffusion solve. Changing it recompiles.">solve iters
243308 <select id="niter">
244309 <option value="0">0</option>
@@ -260,6 +325,8 @@
260325 <option value="255">255</option>
261326 </select>
262327 </label>
328+ </div>
329+ <div class="controls ctrl-group" data-group="display">
263330 <label title="Render on a finer grid. Display only.">display oversampling
264331 <select id="oversample">
265332 <option value="auto" selected>auto</option>
@@ -272,65 +339,53 @@
272339 <label>colormap
273340 <select id="colormap"></select>
274341 </label>
342+ <label title="Blend between the sphere (0) and the surface (1). Display only.">morph
343+ <input type="range" id="morph" min="0" max="1" step="0.01" value="1" />
344+ </label>
345+ <button id="resetview">Reset view</button>
346+ </div>
347+ <div class="controls ctrl-group" data-group="playback">
275348 <button id="runpause" class="primary">Run</button>
349+ <button id="restart" title="Rewind to the initial condition this run started from, without drawing a new one">Restart</button>
350+ </div>
351+ <div class="controls ctrl-group" data-group="benchmark">
276352 <button id="benchmark">Benchmark</button>
277- <label title="Wavelength of the smooth random field the initial condition is seeded from (chebfun's randnfun3, restricted to the surface). An absolute length in the surface's own units, as in chebfun — not a fraction of its size — so smaller means finer features to grow from. Nothing caps it but memory and patience — 0.02 takes about 3 s to seed, 0.01 about 25 s — but it is only *useful* down to about 2*pi/lmax (0.1 at lmax 63): below that the field carries more detail than the grid holds, init's analys discards it, and the seed gets weaker rather than finer. Raise lmax to go finer.">seed λ
353+ </div>
354+ <div class="controls ctrl-group" data-group="seed">
355+ <label title="Wavelength of the smooth random field generating the initial condition">Initial condition wavelength λ
278356 <input id="lam3" type="number" min="0" step="0.05" value="0.5">
279357 </label>
280358 <button id="reseed">Re-seed</button>
281- <button id="resetview">Reset view</button>
282- <button id="movietoggle" title="Export the run as an MP4 movie">Export movie</button>
283- <button id="comparetoggle" title="Run several solver settings side by side on one clock">Compare</button>
284359 </div>
285- <div class="controls" id="comparebar" hidden>
286- <div class="cmp-axes">
287- <div class="cmp-axis">
288- <span title="Iterations of the implicit diffusion solve">solve iters</span>
289- <span id="cmp-niter" class="chips"></span>
290- </div>
291- <div class="cmp-axis">
292- <span>lmax</span>
293- <span id="cmp-lmax" class="chips"></span>
294- </div>
295- <div class="cmp-axis">
296- <span title="Timestep, as a divisor of the model's dt. Divisors keep every variant on the same clock exactly.">dt</span>
297- <span id="cmp-dt" class="chips"></span>
298- </div>
360+ <div class="controls ctrl-group" data-group="movie">
361+ <button id="movietoggle" title="Export the run as an MP4 movie">Export movie</button>
362+ <div class="controls" id="moviebar" hidden>
363+ <label title="Simulation-time units per second of video">movie speed
364+ <select id="moviespeed">
365+ <option value="0.1">0.1×</option>
366+ <option value="0.5">0.5×</option>
367+ <option value="1">1×</option>
368+ <option value="3">3×</option>
369+ <option value="5">5×</option>
370+ <option value="10" selected>10×</option>
371+ <option value="20">20×</option>
372+ </select>
373+ </label>
374+ <label title="Size of each sphere panel in the video, in pixels">resolution
375+ <select id="movieres">
376+ <option value="480">480</option>
377+ <option value="640">640</option>
378+ <option value="768" selected>768</option>
379+ <option value="1080">1080</option>
380+ <option value="1440">1440</option>
381+ </select>
382+ </label>
383+ <label title="Slowly orbit the camera during the movie">
384+ <input type="checkbox" id="movierotate" checked /> auto-rotate
385+ </label>
386+ <button id="movie" title="Replay the run from t = 0 and download an MP4">Export</button>
299387 </div>
300- <label title="The run everything else is measured against">reference
301- <select id="cmp-ref"></select>
302- </label>
303- <button id="cmp-start" class="primary">Compare</button>
304- <span id="cmp-count" class="stats"></span>
305388 </div>
306- <div class="controls" id="moviebar" hidden>
307- <label title="Simulation-time units per second of video">movie speed
308- <select id="moviespeed">
309- <option value="0.1">0.1×</option>
310- <option value="0.5">0.5×</option>
311- <option value="1">1×</option>
312- <option value="3">3×</option>
313- <option value="5">5×</option>
314- <option value="10" selected>10×</option>
315- <option value="20">20×</option>
316- </select>
317- </label>
318- <label title="Size of each sphere panel in the video, in pixels">resolution
319- <select id="movieres">
320- <option value="480">480</option>
321- <option value="640">640</option>
322- <option value="768" selected>768</option>
323- <option value="1080">1080</option>
324- <option value="1440">1440</option>
325- </select>
326- </label>
327- <label title="Slowly orbit the camera during the movie">
328- <input type="checkbox" id="movierotate" checked /> auto-rotate
329- </label>
330- <button id="movie" title="Replay the run from t = 0 and download an MP4">Export</button>
331- </div>
332- <div class="controls" id="params"></div>
333- <div class="controls" id="geomparams"></div>
334389 <p id="geomnote" class="stats"></p>
335390 <div id="panels"></div>
336391 <p class="stats" id="stats"></p>
@@ -367,6 +422,16 @@
367422 </div>
368423 <code id="cmd"></code>
369424 </div>
425+ <details class="cli" id="matlab">
426+ <summary class="cli-head">
427+ <span>The same run as a standalone MATLAB script</span>
428+ <span>
429+ <button id="copymatlab" type="button">Copy</button>
430+ <button id="downloadmatlab" type="button">Download .m</button>
431+ </span>
432+ </summary>
433+ <pre id="matlabscript"></pre>
434+ </details>
370435 <p id="blurb"></p>
371436 <p id="err"></p>
372437 </main>
scripts/ref.tsmodified+9−54View file
@@ -15,8 +15,7 @@
1515 */
1616 import { requestShtDevice, describeAdapter } from '../src/sht/sht.ts';
1717 import { ModelSession } from '../src/mgpu/session.ts';
18-import { mModelByKey, defaultParams, type Params } from '../src/mgpu/registry.ts';
19-import { mGeometryByKey, defaultGeometryParams } from '../src/geom/registry.ts';
18+import { extractReferenceCase, type H5Node } from '../src/compare/referenceCase.ts';
2019 import { relL2, relLinf } from '../src/mgpu/digest.ts';
2120 import { installWebGpu, errMsg, NO_ADAPTER_HINT } from './nodeWebGpu.ts';
2221 import * as h5wasm from 'h5wasm/node';
@@ -89,14 +88,6 @@ for (let i = 0; i < argv.length; i++) {
8988 }
9089 if (!inFile) fail(`--in <file> is required\n\n${USAGE}`, 2);
9190
92-const attrsOf = (entity: { attrs: Record<string, { value: unknown }> }): Record<string, unknown> =>
93- Object.fromEntries(Object.entries(entity.attrs).map(([k, v]) => [k, v.value]));
94-
95-const numberAttrs = (entity: { attrs: Record<string, { value: unknown }> }): Params =>
96- Object.fromEntries(
97- Object.entries(attrsOf(entity)).map(([k, v]) => [k, Number(v)]),
98- );
99-
10091 let device: GPUDevice | null = null;
10192 let session: ModelSession | null = null;
10293 let h5file: InstanceType<typeof h5wasm.File> | null = null;
@@ -104,52 +95,16 @@ let h5file: InstanceType<typeof h5wasm.File> | null = null;
10495 try {
10596 await h5wasm.ready;
10697 h5file = new h5wasm.File(inFile, 'r');
107-
108- const rootAttrs = attrsOf(h5file);
109- const modelKey = String(rootAttrs.model);
110- const model = mModelByKey(modelKey);
111- if (!model) fail(`unknown model '${modelKey}' in ${inFile}`);
112-
113- const specGroup = h5file.get('spec') as InstanceType<typeof h5wasm.Group>;
114- const specAttrs = attrsOf(specGroup);
115- const geometryKey = String(specAttrs.geometry);
116- const geometryModel = mGeometryByKey(geometryKey);
117- if (!geometryModel) fail(`unknown geometry '${geometryKey}' in ${inFile}`);
118-
119- const lmax = Number(specAttrs.lmax);
120- const steps = Number(specAttrs.steps);
121- const niter = niterOverride ?? Number(specAttrs.niter);
122-
123- const params: Params = {
124- ...defaultParams(model),
125- ...numberAttrs(specGroup.get('params') as InstanceType<typeof h5wasm.Group>),
126- };
127- const geometryParams: Params = {
128- ...defaultGeometryParams(geometryModel),
129- ...numberAttrs(specGroup.get('geometry_params') as InstanceType<typeof h5wasm.Group>),
130- };
131-
132- const geomGroup = h5file.get('geometry') as InstanceType<typeof h5wasm.Group>;
133- const fileGeom = {
134- X: (geomGroup.get('Gx') as InstanceType<typeof h5wasm.Dataset>).value as Float32Array,
135- Y: (geomGroup.get('Gy') as InstanceType<typeof h5wasm.Dataset>).value as Float32Array,
136- Z: (geomGroup.get('Gz') as InstanceType<typeof h5wasm.Dataset>).value as Float32Array,
137- };
138-
139- const initialGroup = h5file.get('initial') as InstanceType<typeof h5wasm.Group>;
140- const finalGroup = h5file.get('final') as InstanceType<typeof h5wasm.Group>;
141- const fileInitial: Record<string, Float32Array> = {};
142- const fileFinal: Record<string, Float32Array> = {};
143- for (const name of model.state) {
144- fileInitial[name] = (initialGroup.get(name) as InstanceType<typeof h5wasm.Dataset>)
145- .value as Float32Array;
146- fileFinal[name] = (finalGroup.get(name) as InstanceType<typeof h5wasm.Dataset>)
147- .value as Float32Array;
148- }
149-
98+ const rc = extractReferenceCase(h5file as H5Node, inFile);
15099 h5file.close();
151100 h5file = null;
152101
102+ const { model, geometry: geometryModel, params, geometryParams, lmax, steps } = rc;
103+ const niter = niterOverride ?? rc.niter;
104+ const fileGeom = rc.geometryCoeffs;
105+ const fileInitial = rc.initial;
106+ const fileFinal = rc.final;
107+
153108 const runtime = await installWebGpu();
154109 device = await requestShtDevice().catch((e: unknown) => {
155110 throw new Error(`${errMsg(e)}\n${NO_ADAPTER_HINT}`);
@@ -226,7 +181,7 @@ try {
226181 ` geometry ${geometryModel.label} ` +
227182 geometryModel.params.map((p) => `${p.key}=${geometryParams[p.key]}`).join(' ') +
228183 `\n grid lmax ${lmax} · nlm ${session.sht.nlm}\n` +
229- ` niter ${niter}${niterOverride !== null ? ` (file: ${specAttrs.niter})` : ''}\n` +
184+ ` niter ${niter}${niterOverride !== null ? ` (file: ${rc.niter})` : ''}\n` +
230185 ` run ${steps} steps, dt=${params.dt} (T=${(steps * (params.dt ?? 0)).toFixed(2)})\n`,
231186 );
232187 const fmtErr = (v: { relL2: number; relLinf: number }) =>
scripts/test-node.tsmodified+12−0View file
@@ -8,6 +8,9 @@
88 *
99 * npm run test:node
1010 */
11+import { tmpdir } from 'node:os';
12+import { join } from 'node:path';
13+import * as h5wasm from 'h5wasm/node';
1114 import { requestShtDevice } from '../src/sht/sht.ts';
1215 import { installWebGpu, errMsg, NO_ADAPTER_HINT } from './nodeWebGpu.ts';
1316 import { transformChecks } from '../test/transformChecks.ts';
@@ -16,6 +19,8 @@ import { modelChecks } from '../test/modelChecks.ts';
1619 import { geometryChecks } from '../test/geometryChecks.ts';
1720 import { fluxChecks } from '../test/fluxChecks.ts';
1821 import { compareChecks } from '../test/compareChecks.ts';
22+import { referenceChecks, type H5Rt } from '../test/referenceChecks.ts';
23+import { matlabExportChecks } from '../test/matlabExportChecks.ts';
1924
2025 let failures = 0;
2126 const check = (name: string, ok: boolean, detail: string): void => {
@@ -55,6 +60,13 @@ await modelChecks(device, check, log);
5560 await geometryChecks(device, check, log);
5661 await fluxChecks(device, check, log);
5762 await compareChecks(device, check, log);
63+matlabExportChecks(check, log);
64+await referenceChecks(
65+ h5wasm as unknown as H5Rt,
66+ (name) => join(tmpdir(), `turing-surface-${process.pid}-${name}`),
67+ check,
68+ log,
69+);
5870
5971 console.log(failures === 0 ? '\nAll tests passed.' : `\n${failures} failed.`);
6072 process.exit(failures === 0 ? 0 : 1);
src/compare/compareRun.tsmodified+310−32View file
@@ -41,8 +41,9 @@ import {
4141 import { SphereScene } from '../render/SphereScene.ts';
4242 import { colormaps } from '../render/colormaps.ts';
4343 import { fmtValue, floorRange } from '../render/colorbar.ts';
44-import { sharedModes, sharedNoise } from './sharedStart.ts';
44+import { prolongCoeffs, sharedModes, sharedNoise } from './sharedStart.ts';
4545 import { variantLabel, VARIANT_COLORS, type Variant } from './variants.ts';
46+import type { ReferenceCase } from './referenceCase.ts';
4647
4748 /**
4849 * Latitudes of the shared display grid. 256 is the same target the single-run
@@ -76,8 +77,19 @@ export interface CompareOptions {
7677 geometryParams: Params;
7778 geometrySource: string;
7879 variants: Variant[];
79- /** Index into `variants` of the run everything else is measured against. */
80+ /** Index into `variants` of the run everything else is measured against.
81+ * Ignored when `refFile` is given — the file is the reference then. */
8082 reference: number;
83+ /**
84+ * Check against a reference file instead of against each other: its exact
85+ * initial state seeds every variant (so `seed` and `lam3` go unused), a
86+ * static extra row shows its final state, every Δ is measured against that
87+ * row, and the clock stops at the file's end time. Every variant's lmax must
88+ * be >= the file's — a narrower band could not hold the initial state.
89+ */
90+ refFile?: ReferenceCase;
91+ /** Called when a refFile run reaches the file's end time and stops. */
92+ onFinished?: () => void;
8193 seed: number;
8294 /** Wavelength of the seeded random field, shared by every variant — one
8395 * initial condition means one wavelength as much as one seed. */
@@ -111,9 +123,37 @@ interface Row {
111123 statEl: HTMLElement;
112124 }
113125
126+/**
127+ * The reference file's final state, as one more row of panels — with no
128+ * session behind it: its surface and fields are the file's coefficients
129+ * synthesized once on the shared display grid, fixed for the whole run. Only
130+ * its coloring changes, with the shared range.
131+ */
132+interface FileRow {
133+ coords: Float32Array;
134+ posBuf: Float32Array;
135+ scenes: SphereScene[];
136+ valueBufs: Float32Array[];
137+ colorBufs: Float32Array[];
138+ /** The file's final state on the shared grid, one per species. */
139+ fields: Float32Array[];
140+ /** Its extent, precomputed — a candidate for the shared color range. */
141+ bounds: (Bounds | null)[];
142+}
143+
114144 export class CompareRun {
115145 #opts: CompareOptions;
116146 #rows: Row[] = [];
147+ #fileRow: FileRow | null = null;
148+ /** What restart() reloads: the file's fixed state if opts.refFile is set,
149+ * otherwise a snapshot of the coarsest variant's state as of the last
150+ * (re-)seed — see the capture in create() and in reseed()'s plain branch. */
151+ #initial: Record<string, Float32Array>;
152+ #initialLmax: number;
153+ /** Base steps taken since the initial state — the refFile clock. */
154+ #stepsDone = 0;
155+ /** True once a refFile run has reached the file's end time. */
156+ #finished = false;
117157 #topo: SphereMeshTopology;
118158 /** Quadrature weight per grid point of the shared grid, for the L2 norm. */
119159 #weights: Float64Array;
@@ -137,14 +177,18 @@ export class CompareRun {
137177 private constructor(init: {
138178 opts: CompareOptions;
139179 rows: Row[];
180+ fileRow: FileRow | null;
140181 topo: SphereMeshTopology;
141182 weights: Float64Array;
142183 rangeBars: { fill: (lo: number, hi: number) => void }[];
143184 frameSteps: number;
144185 note: string;
186+ initial: Record<string, Float32Array>;
187+ initialLmax: number;
145188 }) {
146189 this.#opts = init.opts;
147190 this.#rows = init.rows;
191+ this.#fileRow = init.fileRow;
148192 this.#topo = init.topo;
149193 this.#weights = init.weights;
150194 this.#rangeBars = init.rangeBars;
@@ -152,6 +196,8 @@ export class CompareRun {
152196 this.#note = init.note;
153197 this.#morph = init.opts.morph;
154198 this.#ranges = init.opts.model.species.map(() => ({ lo: NaN, hi: NaN }));
199+ this.#initial = init.initial;
200+ this.#initialLmax = init.initialLmax;
155201 }
156202
157203 get variants(): Variant[] {
@@ -168,6 +214,11 @@ export class CompareRun {
168214 return this.#opts.reference;
169215 }
170216
217+ /** The reference file this study is checking against, if any. */
218+ get refFile(): ReferenceCase | null {
219+ return this.#opts.refFile ?? null;
220+ }
221+
171222 /** The base timestep a variant's dtDiv divides. */
172223 static baseDt(params: Params): number {
173224 return params.dt ?? 0;
@@ -182,6 +233,7 @@ export class CompareRun {
182233 // after the grid is up has to take them down explicitly — removing their
183234 // canvases from the DOM would leave both running.
184235 let built: Row[] = [];
236+ let builtFile: FileRow | null = null;
185237
186238 try {
187239 for (let i = 0; i < variants.length; i++) {
@@ -211,7 +263,7 @@ export class CompareRun {
211263
212264 // ---- the shared display grid ----------------------------------------
213265 const maxLmax = Math.max(...variants.map((v) => v.lmax));
214- const panels = variants.length * model.species.length;
266+ const panels = (variants.length + (opts.refFile ? 1 : 0)) * model.species.length;
215267 const target = panels > CROWDED_PANELS ? RENDER_NLAT_CROWDED : RENDER_NLAT;
216268 // Never below what the finest band needs to be representable at all
217269 // (ShtPlan requires nlat > lmax), whatever the panel count says.
@@ -221,12 +273,36 @@ export class CompareRun {
221273 for (const s of sessions) await s.setDisplayGrid(nlat, nphi);
222274
223275 // ---- 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);
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);
276+ // Also what restart() reloads later — the file's fixed state, or (for
277+ // the plain case) a snapshot of the coarsest variant's own state,
278+ // taken after seeding it: the same lowest-lmax session sharedNoise
279+ // itself draws from, so prolonging it up to any other variant later is
280+ // always widening a band, never narrowing one.
281+ let initial: Record<string, Float32Array>;
282+ let initialLmax: number;
283+ if (opts.refFile) {
284+ // The file's exact spectral state, prolonged into each variant's band.
285+ // Exact, not approximate: the state is band-limited at the file's lmax
286+ // and every variant's band contains it, so each session starts from
287+ // the very field the reference run started from.
288+ opts.onStatus('loading the initial state from the reference file…');
289+ for (const s of sessions) {
290+ s.loadState(prolongState(opts.refFile.initial, model.state, opts.refFile.lmax, s.cfg.lmax));
291+ }
292+ initial = opts.refFile.initial;
293+ initialLmax = opts.refFile.lmax;
294+ } else {
295+ opts.onStatus('seeding all variants from one band-limited perturbation…');
296+ const noise = await sharedNoise(sessions, model.seedAmp, opts.seed);
297+ const modes = await sharedModes(sessions[opts.reference] ?? sessions[0], opts.seed);
298+ // One at a time: a seed submits its whole mode sum in pieces, and there
299+ // is nothing to gain from interleaving several variants' worth of it.
300+ for (let i = 0; i < sessions.length; i++) await sessions[i].seedWith(noise[i], modes);
301+ let coarsest = sessions[0];
302+ for (const s of sessions) if (s.cfg.lmax < coarsest.cfg.lmax) coarsest = s;
303+ initial = await coarsest.readState();
304+ initialLmax = coarsest.cfg.lmax;
305+ }
230306
231307 // ---- the mesh, shared; the surface, per variant ---------------------
232308 const view = sessions[0].viewSht;
@@ -256,12 +332,14 @@ export class CompareRun {
256332 frameSteps = Math.max(1, frameSteps);
257333
258334 // ---- the grid of panels ---------------------------------------------
259- const { rows, rangeBars } = await buildGrid(opts, sessions, topo, showDt);
335+ const { rows, fileRow, rangeBars } = await buildGrid(opts, sessions, topo, showDt);
260336 built = rows;
337+ builtFile = fileRow;
261338
262339 const solverGrid = sessions.map((s) => `${s.cfg.nlat}×${s.cfg.nphi}`);
263340 const note =
264- `${variants.length} variants · display grid ${nlat}×${nphi}` +
341+ `${variants.length} variant${variants.length === 1 ? '' : 's'} · ` +
342+ `display grid ${nlat}×${nphi}` +
265343 (sessions.some((s) => s.cfg.nlat > nlat)
266344 ? ` (below the finest solver grid ${solverGrid[solverGrid.length - 1]} — display only)`
267345 : '') +
@@ -269,7 +347,7 @@ export class CompareRun {
269347 ` · ops/step ${ops.join(', ')}`;
270348
271349 const run = new CompareRun({
272- opts, rows, topo, weights, rangeBars, frameSteps, note,
350+ opts, rows, fileRow, topo, weights, rangeBars, frameSteps, note, initial, initialLmax,
273351 });
274352 await run.draw();
275353 run.#observeResize();
@@ -277,6 +355,7 @@ export class CompareRun {
277355 return run;
278356 } catch (e) {
279357 for (const r of built) for (const s of r.scenes) s.dispose();
358+ for (const s of builtFile?.scenes ?? []) s.dispose();
280359 for (const s of sessions) s.destroy();
281360 opts.container.replaceChildren();
282361 opts.container.classList.remove('compare');
@@ -294,27 +373,72 @@ export class CompareRun {
294373 return this.#running;
295374 }
296375
297- /** Re-seed every variant from one new shared perturbation. */
376+ /** Re-seed every variant from one new shared perturbation — or, against a
377+ * reference file, restart from its initial state (there is nothing to
378+ * draw; the seed is ignored). */
298379 async reseed(seed: number): Promise<void> {
299380 const wasRunning = this.#running;
300381 this.#running = false;
301382 while (this.#pumping) await nextFrame();
302383 if (this.#disposed) return;
303384 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++) {
385+ const refFile = this.#opts.refFile;
386+ if (refFile) {
387+ for (const s of sessions) {
388+ s.loadState(prolongState(refFile.initial, this.#opts.model.state, refFile.lmax, s.cfg.lmax));
389+ }
390+ } else {
391+ const noise = await sharedNoise(sessions, this.#opts.model.seedAmp, seed);
392+ const modes = await sharedModes(this.referenceSession ?? sessions[0], seed);
393+ // Checked per variant, not once: a seed awaits its own submission, so a
394+ // dispose can land between two of them and destroy the sessions left.
395+ for (let i = 0; i < sessions.length; i++) {
396+ if (this.#disposed) return;
397+ await sessions[i].seedWith(noise[i], modes);
398+ }
309399 if (this.#disposed) return;
310- await sessions[i].seedWith(noise[i], modes);
400+ // This draw becomes what restart() rewinds to from now on — see the
401+ // identical selection in create(). Recaptured here rather than left
402+ // pointing at the pre-reseed field.
403+ let coarsest = sessions[0];
404+ for (const s of sessions) if (s.cfg.lmax < coarsest.cfg.lmax) coarsest = s;
405+ this.#initial = await coarsest.readState();
406+ this.#initialLmax = coarsest.cfg.lmax;
311407 }
312408 this.#t = 0;
409+ this.#stepsDone = 0;
410+ this.#finished = false;
313411 for (const r of this.#ranges) {
314412 r.lo = NaN;
315413 r.hi = NaN;
316414 }
317415 await this.draw();
416+ this.#status();
417+ if (!this.#disposed && wasRunning) this.setRunning(true);
418+ }
419+
420+ /** Rewind every variant to the saved initial condition — the file's fixed
421+ * state, or (for the plain case) the last (re-)seed, not necessarily the
422+ * very first one — without drawing anything new. */
423+ async restart(): Promise<void> {
424+ const wasRunning = this.#running;
425+ this.#running = false;
426+ while (this.#pumping) await nextFrame();
427+ if (this.#disposed) return;
428+ for (const r of this.#rows) {
429+ r.session.loadState(
430+ prolongState(this.#initial, this.#opts.model.state, this.#initialLmax, r.session.cfg.lmax),
431+ );
432+ }
433+ this.#t = 0;
434+ this.#stepsDone = 0;
435+ this.#finished = false;
436+ for (const r of this.#ranges) {
437+ r.lo = NaN;
438+ r.hi = NaN;
439+ }
440+ await this.draw();
441+ this.#status();
318442 if (!this.#disposed && wasRunning) this.setRunning(true);
319443 }
320444
@@ -333,6 +457,10 @@ export class CompareRun {
333457
334458 /** Model parameters changed. Each variant keeps its own dt. */
335459 setParams(params: Params): void {
460+ // Against a reference file the parameters *are* the file's — they define
461+ // the problem being checked — and the page's parameter panel edits the
462+ // page's own model, which need not even be this one. Nothing to apply.
463+ if (this.#opts.refFile) return;
336464 this.#opts.params = params;
337465 const baseDt = CompareRun.baseDt(params);
338466 for (const r of this.#rows) {
@@ -346,10 +474,15 @@ export class CompareRun {
346474 fillPositions(r.posBuf, r.coords, this.#topo, morph);
347475 for (const s of r.scenes) s.updatePositions(r.posBuf);
348476 }
477+ const f = this.#fileRow;
478+ if (f) {
479+ fillPositions(f.posBuf, f.coords, this.#topo, morph);
480+ for (const s of f.scenes) s.updatePositions(f.posBuf);
481+ }
349482 }
350483
351484 resetView(): void {
352- for (const r of this.#rows) for (const s of r.scenes) s.resetCamera();
485+ for (const s of this.#allScenes()) s.resetCamera();
353486 }
354487
355488 dispose(): void {
@@ -361,11 +494,17 @@ export class CompareRun {
361494 for (const s of r.scenes) s.dispose();
362495 r.session.destroy();
363496 }
497+ for (const s of this.#fileRow?.scenes ?? []) s.dispose();
364498 this.#rows = [];
499+ this.#fileRow = null;
365500 this.#opts.container.replaceChildren();
366501 this.#opts.container.classList.remove('compare');
367502 }
368503
504+ #allScenes(): SphereScene[] {
505+ return [...this.#rows.flatMap((r) => r.scenes), ...(this.#fileRow?.scenes ?? [])];
506+ }
507+
369508 // ----------------------------------------------------------------- drawing
370509 /**
371510 * One frame's readback: every variant's every species, on the shared grid.
@@ -429,7 +568,14 @@ export class CompareRun {
429568 });
430569
431570 for (let k = 0; k < species.length; k++) {
432- const anchor = leastPeak(this.#rows.map((r, i) => (r.healthy ? bounds[i][k] : null)));
571+ // The file row, when there is one, is a candidate like any healthy
572+ // variant: early on the variants' small fields set the scale (it merely
573+ // clips), and if every variant diverges it is the row that keeps the
574+ // grid readable.
575+ const anchor = leastPeak([
576+ ...this.#rows.map((r, i) => (r.healthy ? bounds[i][k] : null)),
577+ this.#fileRow?.bounds[k] ?? null,
578+ ]);
433579 const range = this.#ranges[k];
434580 if (anchor) {
435581 if (!Number.isFinite(range.lo)) {
@@ -456,6 +602,12 @@ export class CompareRun {
456602 fillColors(r.colorBufs[k], r.valueBufs[k], shown.lo, shown.hi, cmap);
457603 r.scenes[k]?.updateColors(r.colorBufs[k]);
458604 }
605+ const f = this.#fileRow;
606+ if (f) {
607+ // Its values never change; only its coloring follows the shared range.
608+ fillColors(f.colorBufs[k], f.valueBufs[k], shown.lo, shown.hi, cmap);
609+ f.scenes[k]?.updateColors(f.colorBufs[k]);
610+ }
459611 }
460612
461613 this.#measureDifference();
@@ -470,8 +622,11 @@ export class CompareRun {
470622 * a physical quantity, which is all it is used for.
471623 */
472624 #measureDifference(): void {
473- const ref = this.#rows[this.#opts.reference];
474- if (!ref) return;
625+ // Against a reference file, every row is measured against its final state;
626+ // otherwise against the chosen reference variant, whose own Δ is zero.
627+ const ref = this.#fileRow ? null : this.#rows[this.#opts.reference];
628+ const refFields = this.#fileRow?.fields ?? ref?.fields;
629+ if (!refFields) return;
475630 const species = this.#opts.model.species;
476631 for (const r of this.#rows) {
477632 for (let k = 0; k < species.length; k++) {
@@ -480,7 +635,7 @@ export class CompareRun {
480635 continue;
481636 }
482637 const a = r.fields[k];
483- const b = ref.fields[k];
638+ const b = refFields[k];
484639 if (!a || !b || a.length !== b.length) {
485640 r.err[k] = NaN;
486641 continue;
@@ -507,7 +662,7 @@ export class CompareRun {
507662 */
508663 #updateRowStats(): void {
509664 const species = this.#opts.model.species;
510- const ref = this.#rows[this.#opts.reference];
665+ const ref = this.#fileRow ? null : this.#rows[this.#opts.reference];
511666 for (const r of this.#rows) {
512667 const per = species
513668 .map((s, k) => `${s} ${Number.isFinite(r.err[k]) ? r.err[k].toExponential(2) : '—'}`)
@@ -525,15 +680,22 @@ export class CompareRun {
525680 }
526681
527682 #status(): void {
683+ const refFile = this.#opts.refFile;
684+ const clock = refFile
685+ ? `<b>t = ${this.#t.toFixed(2)} / ${(refFile.steps * CompareRun.baseDt(this.#opts.params)).toFixed(2)}</b>` +
686+ (this.#finished
687+ ? ` — <b>at the file's end time</b>: Δ is the final comparison against its final state`
688+ : ` · Δ is the distance still to the file's <i>final</i> state — read it at the end time`)
689+ : `<b>t = ${this.#t.toFixed(2)}</b> (same for every variant)`;
528690 this.#opts.onStatus(
529- `<b>t = ${this.#t.toFixed(2)}</b> (same for every variant) · ` +
691+ `${clock} · ` +
530692 (this.#frameMs > 0 ? `${this.#frameMs.toFixed(1)} ms/frame · ` : '') +
531693 this.#note,
532694 );
533695 }
534696
535697 #observeResize(): void {
536- const scenes = this.#rows.flatMap((r) => r.scenes);
698+ const scenes = this.#allScenes();
537699 this.#resizeObs = new ResizeObserver(() => {
538700 for (const s of scenes) {
539701 const box = s.canvas.parentElement;
@@ -560,13 +722,34 @@ export class CompareRun {
560722 this.#pumping = true;
561723 try {
562724 while (this.#running && !this.#disposed) {
725+ // Against a reference file the run is finite: the last frame takes
726+ // however many base steps remain, so every variant lands exactly on
727+ // the file's end time — where Δ against its final state is the
728+ // comparison — and stops there rather than drifting past it.
729+ const refFile = this.#opts.refFile;
730+ const n = refFile
731+ ? Math.min(this.#frameSteps, refFile.steps - this.#stepsDone)
732+ : this.#frameSteps;
733+ if (n <= 0) {
734+ this.#running = false;
735+ this.#opts.onFinished?.();
736+ break;
737+ }
563738 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);
739+ for (const r of this.#rows) r.session.step(n * r.variant.dtDiv);
740+ this.#stepsDone += n;
741+ this.#t += n * CompareRun.baseDt(this.#opts.params);
566742 await this.draw();
567743 if (this.#disposed) break;
568744 const dt = performance.now() - t0;
569745 this.#frameMs = this.#frameMs === 0 ? dt : this.#frameMs + 0.05 * (dt - this.#frameMs);
746+ if (refFile && this.#stepsDone >= refFile.steps) {
747+ this.#finished = true;
748+ this.#running = false;
749+ this.#status();
750+ this.#opts.onFinished?.();
751+ break;
752+ }
570753 this.#status();
571754 await nextFrame();
572755 }
@@ -582,6 +765,19 @@ export class CompareRun {
582765
583766 const nextFrame = (): Promise<number> => new Promise(requestAnimationFrame);
584767
768+/** A whole spectral state re-indexed into a (wider) band's layout — the
769+ * reference file's initial condition, in the form loadState takes. */
770+function prolongState(
771+ coeffs: Record<string, Float32Array>,
772+ names: string[],
773+ lmaxFrom: number,
774+ lmaxTo: number,
775+): Record<string, Float32Array> {
776+ const out: Record<string, Float32Array> = {};
777+ for (const name of names) out[name] = prolongCoeffs(coeffs[name], lmaxFrom, lmaxTo);
778+ return out;
779+}
780+
585781 /** Whether every entry is an ordinary number — false once a variant has left
586782 * its convergence radius and saturated to infinity or NaN. */
587783 function allFinite(f: Float32Array | undefined): boolean {
@@ -626,12 +822,20 @@ function finiteRange(f: Float32Array | undefined): { lo: number; hi: number } |
626822 * the same bar repeated, and would suggest each panel had its own scaling,
627823 * which is exactly the thing that would make the comparison a lie.
628824 */
825+/** The file row's label color — none of the variant palette, since it is not
826+ * a variant: it is the thing they are all measured against. */
827+const FILE_ROW_COLOR = '#57606a';
828+
629829 async function buildGrid(
630830 opts: CompareOptions,
631831 sessions: ModelSession[],
632832 topo: SphereMeshTopology,
633833 showDt: boolean,
634-): Promise<{ rows: Row[]; rangeBars: { fill: (lo: number, hi: number) => void }[] }> {
834+): Promise<{
835+ rows: Row[];
836+ fileRow: FileRow | null;
837+ rangeBars: { fill: (lo: number, hi: number) => void }[];
838+}> {
635839 const { container, model } = opts;
636840 container.replaceChildren();
637841 container.classList.add('compare');
@@ -732,10 +936,84 @@ async function buildGrid(
732936 });
733937 }
734938
939+ // ---- the reference file's final state, as one more (static) row ---------
940+ let fileRow: FileRow | null = null;
941+ if (opts.refFile) {
942+ const rf = opts.refFile;
943+ // Synthesized through the coarsest session's display plan — exact, like
944+ // every other use of the shared grid: the file's coefficients are
945+ // band-limited at its lmax, which every variant's band contains.
946+ const view = sessions[0].viewSht;
947+ const lmaxTo = sessions[0].cfg.lmax;
948+ const on = (q: Float32Array): Promise<Float32Array> =>
949+ view.synth(prolongCoeffs(q, rf.lmax, lmaxTo));
950+ const [gx, gy, gz] = [
951+ await on(rf.geometryCoeffs.X),
952+ await on(rf.geometryCoeffs.Y),
953+ await on(rf.geometryCoeffs.Z),
954+ ];
955+ // The file's own surface, not a regeneration of it — interleaved xyz, the
956+ // same layout renderPositions() hands back.
957+ const coords = new Float32Array(3 * gx.length);
958+ for (let i = 0; i < gx.length; i++) {
959+ coords[3 * i] = gx[i];
960+ coords[3 * i + 1] = gy[i];
961+ coords[3 * i + 2] = gz[i];
962+ }
963+ const posBuf = new Float32Array(topo.numVertices * 3);
964+ fillPositions(posBuf, coords, topo, opts.morph);
965+
966+ const rowEl = document.createElement('div');
967+ rowEl.className = 'cmp-row';
968+ const labelEl = document.createElement('div');
969+ labelEl.className = 'cmp-rowlabel';
970+ labelEl.style.setProperty('--c', FILE_ROW_COLOR);
971+ const nameEl = document.createElement('div');
972+ nameEl.className = 'cmp-rowname';
973+ nameEl.textContent = 'reference file';
974+ nameEl.title = rf.label;
975+ const statEl = document.createElement('div');
976+ statEl.className = 'cmp-rowstat';
977+ statEl.innerHTML = `${rf.steps.toLocaleString()} steps<br><b>final state</b>`;
978+ labelEl.append(nameEl, statEl);
979+ const colsEl = document.createElement('div');
980+ colsEl.className = 'cmp-cols';
981+ rowEl.append(labelEl, colsEl);
982+ container.append(rowEl);
983+
984+ const scenes: SphereScene[] = [];
985+ const valueBufs: Float32Array[] = [];
986+ const colorBufs: Float32Array[] = [];
987+ const fields: Float32Array[] = [];
988+ const bounds: (Bounds | null)[] = [];
989+ for (let k = 0; k < model.species.length; k++) {
990+ const box = document.createElement('div');
991+ box.className = 'sphere-box cmp-box';
992+ colsEl.append(box);
993+ const scene = new SphereScene(
994+ box,
995+ topo.numVertices,
996+ topo.indices,
997+ Float32Array.from(posBuf),
998+ sphereBg || undefined,
999+ );
1000+ scene.fitCamera();
1001+ scenes.push(scene);
1002+ const field = await on(rf.final[model.state[k]]);
1003+ fields.push(field);
1004+ bounds.push(finiteRange(field));
1005+ const valueBuf = new Float32Array(topo.numVertices);
1006+ fillFieldValues(valueBuf, field, topo);
1007+ valueBufs.push(valueBuf);
1008+ colorBufs.push(new Float32Array(topo.numVertices * 3));
1009+ }
1010+ fileRow = { coords, posBuf, scenes, valueBufs, colorBufs, fields, bounds };
1011+ }
1012+
7351013 // Every panel shares one camera: the study is about the fields, and looking
7361014 // at two of them from different angles is not comparing them.
737- const all = rows.flatMap((r) => r.scenes);
1015+ const all = [...rows.flatMap((r) => r.scenes), ...(fileRow?.scenes ?? [])];
7381016 for (let i = 1; i < all.length; i++) all[0].syncCamerasWith(all[i]);
7391017
740- return { rows, rangeBars };
1018+ return { rows, fileRow, rangeBars };
7411019 }
src/compare/referenceCase.tsadded+124−0View file
@@ -0,0 +1,124 @@
1+/**
2+ * Reading a reference HDF5 file into the pieces a replay needs.
3+ *
4+ * A reference file is a saved run from an independently-implemented solver —
5+ * geometry, initial and final spherical-harmonic coefficients, and the run's
6+ * parameters — in the layout documented in docs/ellipsoid-reference-spec.md.
7+ * Two things read it: the `npm run ref` CLI (through `h5wasm/node`) and the
8+ * browser's compare mode (through `h5wasm`, lazily loaded — see
9+ * referenceFile.ts). Both hand this module the same object shape, so the
10+ * format knowledge lives once.
11+ */
12+import { mModelByKey, defaultParams, type MModel, type Params } from '../mgpu/registry.ts';
13+import { mGeometryByKey, defaultGeometryParams, type MGeometry } from '../geom/registry.ts';
14+import { nlmCalc } from '../sht/layout.ts';
15+
16+/** The slice of h5wasm's File/Group/Dataset API this reader touches — enough
17+ * that the node and browser builds both satisfy it structurally. */
18+export interface H5Node {
19+ attrs: Record<string, { value: unknown }>;
20+ get(name: string): unknown;
21+}
22+
23+export interface ReferenceCase {
24+ /** Where it came from — the file name, for labels and messages. */
25+ label: string;
26+ model: MModel;
27+ geometry: MGeometry;
28+ /** The model's defaults overlaid with the file's own — `dt` included, so
29+ * `steps * params.dt` is the file's end time. */
30+ params: Params;
31+ geometryParams: Params;
32+ lmax: number;
33+ /** The solve-iteration count recorded in the file — the replay's default. */
34+ niter: number;
35+ /** Steps at `params.dt` from the initial state to the final one. */
36+ steps: number;
37+ /** The band-limited surface's own coefficients, [re, im] per (l, m). The
38+ * reference solver ran on this exact surface, not the analytic shape. */
39+ geometryCoeffs: { X: Float32Array; Y: Float32Array; Z: Float32Array };
40+ /** Spectral state per species (keyed by `model.state` name) at t = 0. */
41+ initial: Record<string, Float32Array>;
42+ /** The same, at the end time. */
43+ final: Record<string, Float32Array>;
44+}
45+
46+const attrsOf = (node: H5Node): Record<string, unknown> =>
47+ Object.fromEntries(Object.entries(node.attrs).map(([k, v]) => [k, v.value]));
48+
49+/** Attributes as numbers — h5wasm hands back number or BigInt by dtype. */
50+const numberAttrs = (node: H5Node): Params =>
51+ Object.fromEntries(Object.entries(attrsOf(node)).map(([k, v]) => [k, Number(v)]));
52+
53+function groupOf(node: H5Node, name: string): H5Node {
54+ const g = node.get(name) as H5Node | null;
55+ if (!g || typeof g.get !== 'function') {
56+ throw new Error(`no '${name}/' group — is this a reference file?`);
57+ }
58+ return g;
59+}
60+
61+function coeffsOf(group: H5Node, groupName: string, name: string, nlm: number): Float32Array {
62+ const v = (group.get(name) as { value?: unknown } | null)?.value;
63+ if (!(v instanceof Float32Array)) {
64+ throw new Error(`'${groupName}/${name}' is not a float32 dataset`);
65+ }
66+ if (v.length !== 2 * nlm) {
67+ throw new Error(`'${groupName}/${name}' has ${v.length} values, expected 2*nlm = ${2 * nlm}`);
68+ }
69+ return v;
70+}
71+
72+/** Read an open reference file. Throws with a plain message on anything the
73+ * replay could not act on — unknown model or geometry, missing or misshapen
74+ * coefficients — so both the CLI and the page can just show it. */
75+export function extractReferenceCase(file: H5Node, label: string): ReferenceCase {
76+ const modelKey = String(attrsOf(file).model);
77+ const model = mModelByKey(modelKey);
78+ if (!model) throw new Error(`unknown model '${modelKey}'`);
79+
80+ const spec = groupOf(file, 'spec');
81+ const specAttrs = attrsOf(spec);
82+ const geometryKey = String(specAttrs.geometry);
83+ const geometry = mGeometryByKey(geometryKey);
84+ if (!geometry) throw new Error(`unknown geometry '${geometryKey}'`);
85+
86+ const lmax = Number(specAttrs.lmax);
87+ const steps = Number(specAttrs.steps);
88+ const niter = Number(specAttrs.niter);
89+ if (!Number.isInteger(lmax) || lmax < 1) throw new Error(`bad lmax '${String(specAttrs.lmax)}'`);
90+ if (!Number.isInteger(steps) || steps < 1) throw new Error(`bad steps '${String(specAttrs.steps)}'`);
91+ if (!Number.isInteger(niter) || niter < 0) throw new Error(`bad niter '${String(specAttrs.niter)}'`);
92+ const nlm = nlmCalc(lmax, lmax);
93+
94+ const params: Params = {
95+ ...defaultParams(model),
96+ ...numberAttrs(groupOf(spec, 'params')),
97+ };
98+ if (!(params.dt! > 0)) throw new Error(`bad dt '${params.dt}'`);
99+ const geometryParams: Params = {
100+ ...defaultGeometryParams(geometry),
101+ ...numberAttrs(groupOf(spec, 'geometry_params')),
102+ };
103+
104+ const geom = groupOf(file, 'geometry');
105+ const geometryCoeffs = {
106+ X: coeffsOf(geom, 'geometry', 'Gx', nlm),
107+ Y: coeffsOf(geom, 'geometry', 'Gy', nlm),
108+ Z: coeffsOf(geom, 'geometry', 'Gz', nlm),
109+ };
110+
111+ const initialGroup = groupOf(file, 'initial');
112+ const finalGroup = groupOf(file, 'final');
113+ const initial: Record<string, Float32Array> = {};
114+ const final: Record<string, Float32Array> = {};
115+ for (const name of model.state) {
116+ initial[name] = coeffsOf(initialGroup, 'initial', name, nlm);
117+ final[name] = coeffsOf(finalGroup, 'final', name, nlm);
118+ }
119+
120+ return {
121+ label, model, geometry, params, geometryParams,
122+ lmax, niter, steps, geometryCoeffs, initial, final,
123+ };
124+}
src/compare/referenceFile.tsadded+29−0View file
@@ -0,0 +1,29 @@
1+/**
2+ * Reading a reference .h5 in the page.
3+ *
4+ * h5wasm's browser build carries the whole HDF5 library as embedded wasm —
5+ * about 4 MB — so it is imported here, dynamically, and nowhere else: the page
6+ * pays for it on the first file actually loaded, never on startup. The bytes
7+ * are written into the wasm module's in-memory filesystem under a fixed
8+ * scratch name (loads are sequential — there is one file input), opened,
9+ * extracted, and unlinked.
10+ */
11+import { extractReferenceCase, type H5Node, type ReferenceCase } from './referenceCase.ts';
12+
13+const SCRATCH = '/loaded-reference.h5';
14+
15+export async function loadReferenceFile(file: File): Promise<ReferenceCase> {
16+ const bytes = new Uint8Array(await file.arrayBuffer());
17+ const h5 = await import('h5wasm');
18+ const { FS } = (await h5.ready) as unknown as {
19+ FS: { writeFile(path: string, data: Uint8Array): void; unlink(path: string): void };
20+ };
21+ FS.writeFile(SCRATCH, bytes);
22+ const opened = new h5.File(SCRATCH, 'r');
23+ try {
24+ return extractReferenceCase(opened as unknown as H5Node, file.name);
25+ } finally {
26+ opened.close();
27+ FS.unlink(SCRATCH);
28+ }
29+}
src/export/matlabScript.tsadded+313−0View file
@@ -0,0 +1,313 @@
1+/**
2+ * The run on screen, as one standalone MATLAB script.
3+ *
4+ * The models and geometries are already MATLAB; what the app supplies around
5+ * them — the transforms, the geometry weights, the seeded field, the driver —
6+ * exists only as TypeScript and WGSL. This module assembles a single function
7+ * file carrying all of it: the current model and geometry sources verbatim as
8+ * local functions, double-precision MATLAB ports of the host-provided
9+ * operations (support.m), and a generated driver with the run's settings
10+ * baked in.
11+ *
12+ * Fidelity is method-for-method, not bit-for-bit: the ports run in f64 where
13+ * the GPU path is f32, and random draws use MATLAB's own rng, so a seed value
14+ * selects a different member of the same random ensemble than the same value
15+ * in the app. The script's results file uses the app's reference-run layout
16+ * (docs/ellipsoid-reference-spec.md), so a MATLAB run can be loaded back into
17+ * the page or checked with `npm run ref`.
18+ */
19+import supportSource from './support.m?raw';
20+import randnfun3Source from '../../tools/randnfun3.m?raw';
21+import randnfunsphereSource from '../../tools/randnfunsphere.m?raw';
22+import type { MModel, Params } from '../mgpu/registry.ts';
23+import type { MGeometry } from '../geom/registry.ts';
24+
25+/** The generated function's name — and therefore the file name to save as. */
26+export const MATLAB_SCRIPT_NAME = 'turing_surface_run';
27+
28+export interface MatlabExportSpec {
29+ model: MModel;
30+ /** Model source as running — the editor's working copy when edited. */
31+ modelSource: string;
32+ params: Params;
33+ geometry: MGeometry;
34+ geometrySource: string;
35+ geometryParams: Params;
36+ lmax: number;
37+ niter: number;
38+ /** Wavelength of the seeded random field. */
39+ lam3: number;
40+ seed: number;
41+ /** Preset key, recorded in the results file's /spec. */
42+ preset: string;
43+ /** The equivalent `npm run bench` command, recorded for provenance. */
44+ command: string;
45+ /** The generated script's own run controls; app defaults when omitted. */
46+ controls?: { nsteps?: number; plotEvery?: number; outFile?: string };
47+}
48+
49+interface Signature {
50+ outputs: string[];
51+ params: string[];
52+}
53+
54+/** First `function [outs] = name(args)` line in a .m — the same contract the
55+ * compiler applies, minus everything it checks later. */
56+function parseSignature(source: string, name: string, file: string): Signature {
57+ const re = new RegExp(
58+ String.raw`^[ \t]*function\s+(?:\[([^\]]*)\]|([A-Za-z]\w*))\s*=\s*${name}\s*\(([^)]*)\)`,
59+ 'm',
60+ );
61+ const m = re.exec(source);
62+ if (!m) {
63+ throw new Error(`cannot export: ${file} defines no function named '${name}'`);
64+ }
65+ const split = (s: string): string[] =>
66+ s.split(',').map((t) => t.trim()).filter((t) => t.length > 0);
67+ return {
68+ outputs: m[1] !== undefined ? split(m[1]) : [m[2]],
69+ params: split(m[3]),
70+ };
71+}
72+
73+/** A number as MATLAB source. JS stringification round-trips doubles exactly
74+ * and every form it produces (0.0004, 1e-21, -3) is a MATLAB literal. */
75+const num = (v: number): string => (Number.isFinite(v) ? String(v) : '0');
76+
77+/** A string as a MATLAB char literal. */
78+const str = (s: string): string => `'${s.replace(/'/g, "''")}'`;
79+
80+const banner = (title: string): string => {
81+ const line = `% ${'='.repeat(72)}`;
82+ return `${line}\n% ${title}\n${line}`;
83+};
84+
85+export function generateMatlabScript(spec: MatlabExportSpec): string {
86+ const { model, geometry } = spec;
87+ const controls = {
88+ nsteps: spec.controls?.nsteps ?? 2000,
89+ plotEvery: spec.controls?.plotEvery ?? 10,
90+ outFile: spec.controls?.outFile ?? `${MATLAB_SCRIPT_NAME}.h5`,
91+ };
92+
93+ const init = parseSignature(spec.modelSource, 'init', `models/${model.key}.m`);
94+ const step = parseSignature(spec.modelSource, 'step', `models/${model.key}.m`);
95+ const shape = parseSignature(spec.geometrySource, 'shape', `geometries/${geometry.key}.m`);
96+
97+ // The driver defines every host-provided name the .m may ask for (lam,
98+ // filt, the geometry fields, jhat, niter, ...) under its canonical name,
99+ // so a model call is its own signature read back. Only the tunable
100+ // parameters live elsewhere — in the mp/gp structs, where the person
101+ // running the script edits them — so those names are mapped.
102+ const modelParamKeys = new Set(model.params.map((p) => p.key));
103+ const modelArg = (a: string): string => (modelParamKeys.has(a) ? `mp.${a}` : a);
104+ const shapeArg = (a: string): string =>
105+ a === 'theta' || a === 'phi' ? a : `gp.${a}`;
106+
107+ const stateOuts = [...model.state, ...model.species];
108+ const outs = `[${stateOuts.join(', ')}]`;
109+ const initCall = `${outs} = init(${init.params.map(modelArg).join(', ')});`;
110+ const stepCall = `${outs} = step(${step.params.map(modelArg).join(', ')});`;
111+ const shapeCall = `[gxr, gyr, gzr] = shape(${shape.params.map(shapeArg).join(', ')});`;
112+
113+ const speciesCell = `{${model.species.join(', ')}}`;
114+ const namesCell = `{${model.species.map((s) => str(s)).join(', ')}}`;
115+
116+ // `noise` is the plain seeded grid perturbation, for a .m that takes it
117+ // instead of calling randnfun3 (none of the shipped models do).
118+ const takesNoise = init.params.includes('noise') || step.params.includes('noise');
119+
120+ const mpBlock = model.params
121+ .map((p) => `mp.${p.key} = ${num(spec.params[p.key] ?? p.value)};`)
122+ .join('\n');
123+ const gpBlock = geometry.params
124+ .map((p) => `gp.${p.key} = ${num(spec.geometryParams[p.key] ?? p.value)};`)
125+ .join('\n');
126+
127+ const driver = `function ${MATLAB_SCRIPT_NAME}()
128+% ${model.label} on ${geometry.label} -- a run captured from the
129+% turing-surface app as one standalone MATLAB script.
130+%
131+% The model and geometry .m below are the app's own, verbatim; around them
132+% this file carries double-precision MATLAB ports of everything the app
133+% provides from the host side: the spherical-harmonic transforms and their
134+% derivative shuffles, the metric weights of the surface Laplace-Beltrami
135+% operator, the seeded random field, and the run loop (src/sht and src/geom
136+% in the repository). The scheme is the app's: IMEX Euler, implicit
137+% diffusion preconditioned on the round sphere, the geometric correction
138+% iterated niter times per step.
139+%
140+% Two deliberate differences from the page. Everything here runs in double
141+% precision, where the app's GPU path is single. And random draws use
142+% MATLAB's own rng, so a seed value selects a different member of the same
143+% random ensemble than the same value in the app.
144+%
145+% Save as ${MATLAB_SCRIPT_NAME}.m and run it. The run plots live, and the
146+% final state is written to an HDF5 file in the app's reference-run layout
147+% (docs/ellipsoid-reference-spec.md in the repository), so it can be loaded
148+% back into the page ("Compare against uploaded data") or checked on the
149+% desktop with \`npm run ref -- --in ${controls.outFile}\`.
150+% Needs base MATLAB, R2020b or newer; no toolboxes.
151+
152+% ---- run controls --------------------------------------------------------
153+nsteps = ${controls.nsteps}; % timesteps to run
154+plot_every = ${controls.plotEvery}; % live-plot interval, in steps; 0 disables plotting
155+out_file = ${str(controls.outFile)}; % results file; '' disables
156+seed = ${num(spec.seed)}; % rng seed for the initial condition
157+
158+% ---- captured from the app -----------------------------------------------
159+lmax = ${spec.lmax}; % spherical-harmonic truncation degree
160+niter = ${spec.niter}; % iterations of the implicit solve's geometric correction
161+lam3 = ${num(spec.lam3)}; % wavelength of the seeded random field
162+${model.params.length ? `% ${model.label} parameters\n${mpBlock}` : `% ${model.label} has no parameters`}
163+${geometry.params.length ? `% ${geometry.label} parameters\n${gpBlock}` : `% ${geometry.label} has no parameters`}
164+
165+% ---- grid and transforms -------------------------------------------------
166+% nlat/nphi follow lmax by the app's dealiasing rule (src/sht/layout.ts) for
167+% a reaction of polynomial degree pdeg.
168+pdeg = ${model.pdeg};
169+mmax = lmax;
170+nlat = 2 * ceil(max(lmax + 1, ((pdeg + 1) * lmax + 1) / 2) / 2);
171+nphi = 2 ^ nextpow2((pdeg + 1) * lmax + 1);
172+npts = nlat * nphi;
173+sht_tables(sht_setup(lmax, mmax, nlat, nphi));
174+S = sht_tables();
175+nlm = S.nlm;
176+lam = S.lam;
177+filt = S.filt;
178+theta = S.theta;
179+phi = S.phi;
180+
181+% ---- the surface ---------------------------------------------------------
182+${shapeCall}
183+% A constant coordinate comes back scalar; spread it over the grid.
184+gxr = gxr + zeros(npts, 1);
185+gyr = gyr + zeros(npts, 1);
186+gzr = gzr + zeros(npts, 1);
187+G = surface_tables(gxr, gyr, gzr);
188+gx = G.gx; gy = G.gy; gz = G.gz;
189+Gx = G.Gx; Gy = G.Gy; Gz = G.Gz;
190+p1 = G.p1; p2 = G.p2; q2 = G.q2; r = G.r;
191+dp1 = G.dp1; dq2 = G.dq2; jinv = G.jinv;
192+Vtx = G.Vtx; Vty = G.Vty; Vtz = G.Vtz;
193+Vpx = G.Vpx; Vpy = G.Vpy; Vpz = G.Vpz;
194+jhat = G.Jhat;
195+radius = sqrt(gx.^2 + gy.^2 + gz.^2);
196+fprintf('grid %d x %d, nlm %d, radius %.3f-%.3f, Jhat %.3f\\n', ...
197+ nlat, nphi, nlm, min(radius), max(radius), jhat);
198+
199+% ---- initial condition ---------------------------------------------------
200+rng(seed);
201+${takesNoise ? `noise = ${num(model.seedAmp)} * randn(npts, 1);\n` : ''}${initCall}
202+${model.state.map((s) => `${s}0 = ${s};`).join('\n')}
203+
204+% ---- time loop -----------------------------------------------------------
205+if plot_every > 0
206+ ph = plot_setup(gx, gy, gz, ${speciesCell}, ${namesCell});
207+ plot_update(ph, ${speciesCell}, 0, 0, nsteps);
208+end
209+report_every = max(1, round(nsteps / 10));
210+t = 0;
211+tstart = tic;
212+for k = 1:nsteps
213+ ${stepCall}
214+ t = t + mp.dt;
215+ if plot_every > 0 && (mod(k, plot_every) == 0 || k == nsteps)
216+ plot_update(ph, ${speciesCell}, t, k, nsteps);
217+ end
218+ if mod(k, report_every) == 0 || k == nsteps
219+ fprintf('step %d/%d t = %.3f (%.1f s)\\n', k, nsteps, t, toc(tstart));
220+ end
221+end
222+
223+% ---- results file --------------------------------------------------------
224+% The app's reference-run layout, plus a /fields group with the final grid
225+% fields, the surface and the grid angles (each field stored nphi x nlat,
226+% ring by ring from the north pole).
227+if ~isempty(out_file)
228+ if exist(out_file, 'file') == 2
229+ delete(out_file);
230+ end
231+${['Gx', 'Gy', 'Gz']
232+ .map((c) => ` write_coeffs(out_file, '/geometry/${c}', ${c});`)
233+ .join('\n')}
234+${model.state
235+ .map((s) => ` write_coeffs(out_file, '/initial/${s}', ${s}0);`)
236+ .join('\n')}
237+${model.state
238+ .map((s) => ` write_coeffs(out_file, '/final/${s}', ${s});`)
239+ .join('\n')}
240+${[...model.species.map((s) => [s, s] as const), (['x', 'gx'] as const), (['y', 'gy'] as const), (['z', 'gz'] as const)]
241+ .map(
242+ ([name, v]) =>
243+ ` h5create(out_file, '/fields/${name}', [nphi nlat]);\n` +
244+ ` h5write(out_file, '/fields/${name}', reshape(${v}, nphi, nlat));`,
245+ )
246+ .join('\n')}
247+ h5create(out_file, '/fields/theta', nlat);
248+ h5write(out_file, '/fields/theta', acos(min(1, max(-1, S.ct))));
249+ h5create(out_file, '/fields/phi', nphi);
250+ h5write(out_file, '/fields/phi', 2*pi*(0:nphi-1)'/nphi);
251+ make_group(out_file, '/backend');
252+ make_group(out_file, '/spec');
253+ make_group(out_file, '/spec/params');
254+ make_group(out_file, '/spec/geometry_params');
255+ make_group(out_file, '/grid');
256+ h5writeatt(out_file, '/', 'model', ${str(model.key)});
257+ h5writeatt(out_file, '/', 'species', [${model.state.map((s) => `"${s}"`).join(' ')}]);
258+ h5writeatt(out_file, '/', 'command', ${str(spec.command)});
259+ h5writeatt(out_file, '/backend', 'runtime', 'matlab');
260+ h5writeatt(out_file, '/backend', 'adapter', ['MATLAB ' version]);
261+ h5writeatt(out_file, '/backend', 'precision', 'double');
262+ h5writeatt(out_file, '/spec', 'preset', ${str(spec.preset)});
263+ h5writeatt(out_file, '/spec', 'geometry', ${str(geometry.key)});
264+ h5writeatt(out_file, '/spec', 'lmax', lmax);
265+ h5writeatt(out_file, '/spec', 'seed', seed);
266+ h5writeatt(out_file, '/spec', 'steps', nsteps);
267+ h5writeatt(out_file, '/spec', 'warmup', 0);
268+ h5writeatt(out_file, '/spec', 'niter', niter);
269+ h5writeatt(out_file, '/spec', 'lam3', lam3);
270+${model.params
271+ .map((p) => ` h5writeatt(out_file, '/spec/params', ${str(p.key)}, mp.${p.key});`)
272+ .join('\n')}
273+${geometry.params
274+ .map((p) => ` h5writeatt(out_file, '/spec/geometry_params', ${str(p.key)}, gp.${p.key});`)
275+ .join('\n')}
276+ h5writeatt(out_file, '/grid', 'lmax', lmax);
277+ h5writeatt(out_file, '/grid', 'mmax', mmax);
278+ h5writeatt(out_file, '/grid', 'nlat', nlat);
279+ h5writeatt(out_file, '/grid', 'nphi', nphi);
280+ h5writeatt(out_file, '/grid', 'nlm', nlm);
281+ fprintf('wrote %s\\n', out_file);
282+end
283+end`;
284+
285+ // tools/randnfun3.m verbatim, renamed: the models call the app's builtin
286+ // `randnfun3(lam3, gx, gy, gz)`, which support.m provides as a dispatcher
287+ // over this mode draw.
288+ const modesSource = randnfun3Source.replace(
289+ /function\s*\[\s*k\s*,\s*c\s*\]\s*=\s*randnfun3\s*\(/,
290+ 'function [k, c] = randnfun3_modes(',
291+ );
292+ if (modesSource === randnfun3Source) {
293+ throw new Error('cannot export: tools/randnfun3.m no longer matches the expected signature');
294+ }
295+
296+ const usesSphere = /\brandnfunsphere\b/.test(spec.geometrySource + spec.modelSource);
297+
298+ const parts = [
299+ driver,
300+ banner(`models/${model.key}.m -- the model, verbatim`),
301+ spec.modelSource.trim(),
302+ banner(`geometries/${geometry.key}.m -- the surface, verbatim`),
303+ spec.geometrySource.trim(),
304+ banner('tools/randnfun3.m -- the random-field mode draw, verbatim'),
305+ modesSource.trim(),
306+ ...(usesSphere
307+ ? [banner('tools/randnfunsphere.m -- verbatim'), randnfunsphereSource.trim()]
308+ : []),
309+ banner('host-provided operations, ported from src/sht and src/geom'),
310+ supportSource.trim(),
311+ ];
312+ return parts.join('\n\n') + '\n';
313+}
src/export/support.madded+388−0View file
@@ -0,0 +1,388 @@
1+% ---------------------------------------------------------------- transforms
2+%
3+% Double-precision MATLAB ports of the operations the app provides to a .m
4+% around its compiled GPU pipeline. Conventions follow src/sht/layout.ts:
5+% orthonormal spherical harmonics with the Condon-Shortley phase, coefficients
6+% stored for m >= 0 only in m-major order (m = 0..mmax, l = m..lmax within
7+% each m) -- here as complex nlm x 1 column vectors where the GPU carries
8+% interleaved [re, im] pairs. Grid fields are npts x 1 columns, phi-fastest:
9+% point (itheta, iphi) sits at row (itheta-1)*nphi + iphi, north row first.
10+
11+% Holds the precomputed tables between calls: set once from the top of the
12+% run, read back by every transform below.
13+function S = sht_tables(S)
14+ persistent stored
15+ if nargin > 0
16+ stored = S;
17+ end
18+ S = stored;
19+end
20+
21+% Everything the transforms need for one grid: Gauss nodes and weights,
22+% per-m Legendre tables, the coefficient layout, the derivative shuffles,
23+% and the eigenvalue/filter vectors the models take as `lam` and `filt`.
24+function S = sht_setup(lmax, mmax, nlat, nphi)
25+ S.lmax = lmax;
26+ S.mmax = mmax;
27+ S.nlat = nlat;
28+ S.nphi = nphi;
29+ S.npts = nlat * nphi;
30+ [ct, wg] = gauss_legendre(nlat);
31+ S.ct = ct;
32+ S.st = sqrt(1 - ct.^2);
33+ S.wg = wg;
34+ S.nlm = (mmax + 1) * (lmax + 1) - mmax * (mmax + 1) / 2;
35+
36+ % The grid angles as npts x 1 fields, phi-fastest like everything else.
37+ S.theta = repelem(acos(min(1, max(-1, ct))), nphi);
38+ S.phi = repmat(2*pi*(0:nphi-1)'/nphi, nlat, 1);
39+ S.stpt = repelem(S.st, nphi);
40+
41+ % Degree and order of each coefficient, and each m block's start.
42+ off = zeros(mmax + 1, 1);
43+ lv = zeros(S.nlm, 1);
44+ mv = zeros(S.nlm, 1);
45+ pos = 1;
46+ for m = 0:mmax
47+ n = lmax - m + 1;
48+ off(m + 1) = pos;
49+ lv(pos:pos + n - 1) = (m:lmax)';
50+ mv(pos:pos + n - 1) = m;
51+ pos = pos + n;
52+ end
53+ S.off = off;
54+ S.lv = lv;
55+ S.mv = mv;
56+ % Laplace-Beltrami eigenvalues l(l+1) and the top-mode filter: 1 below
57+ % lmax-2, 0 at the top two degrees, where the derivative recurrences cannot
58+ % exactly represent a derivative (src/mgpu/model.ts).
59+ S.lam = lv .* (lv + 1);
60+ S.filt = double(lv < lmax - 2);
61+
62+ % Orthonormal Legendre tables ytilde_l^m(theta_i), one nlat x (lmax-m+1)
63+ % block per m, by the standard three-term recurrence (src/sht/coeffs.ts;
64+ % SHTNS normalization, Condon-Shortley phase carried in the seed's sign).
65+ S.Y = cell(mmax + 1, 1);
66+ t = 1 / (4*pi);
67+ amm = sqrt(t);
68+ for m = 0:mmax
69+ if m > 0
70+ t = t * (2*m + 1) / (2*m);
71+ amm = (-1)^m * sqrt(t);
72+ end
73+ n = lmax - m + 1;
74+ Y = zeros(nlat, n);
75+ y0 = amm * S.st.^m;
76+ Y(:, 1) = y0;
77+ if n > 1
78+ y1 = sqrt(2*m + 3) * ct .* y0;
79+ Y(:, 2) = y1;
80+ for l = m + 2:lmax
81+ t1 = (l + m) * (l - m);
82+ a = sqrt((2*l + 1) * (2*l - 1) / t1);
83+ b = -sqrt(((2*l + 1) / (2*l - 3)) * ((l - 1 + m) * (l - 1 - m) / t1));
84+ y2 = a * ct .* y1 + b * y0;
85+ Y(:, l - m + 1) = y2;
86+ y0 = y1;
87+ y1 = y2;
88+ end
89+ end
90+ S.Y{m + 1} = Y;
91+ end
92+
93+ % sin(theta)*dtheta in coefficient space: v_l^m = ap(lm) u_{l-1}^m +
94+ % am(lm) u_{l+1}^m (src/sht/derivCoeffs.ts). Neighbors sit at +-1 within
95+ % each m block; ap/am are zero at the block edges, so the clamped index
96+ % vectors never read across a boundary.
97+ l = lv;
98+ m = mv;
99+ ap = (l - 1) .* sqrt(max(0, (l - m) .* (l + m)) ./ ((2*l - 1) .* (2*l + 1)));
100+ ap(l <= m) = 0;
101+ am = -(l + 2) .* sqrt((l + 1 - m) .* (l + 1 + m) ./ ((2*l + 1) .* (2*l + 3)));
102+ am(l >= lmax) = 0;
103+ S.ap = ap;
104+ S.am = am;
105+ S.iprev = max((1:S.nlm)' - 1, 1);
106+ S.inext = min((1:S.nlm)' + 1, S.nlm);
107+
108+ % dphig's Fourier multiplier: i*m on fft's frequency layout, masked past
109+ % the filter's reach (mcut = lmax-3), mirroring src/sht/wgsl/deriv.ts.
110+ freq = [(0:nphi/2)'; (-nphi/2 + 1:-1)'];
111+ S.dmul = 1i * freq .* (abs(freq) <= max(0, lmax - 3));
112+end
113+
114+% Gauss-Legendre nodes cos(theta), in decreasing order (north pole first),
115+% and weights for integration over cos(theta) -- Newton iteration on P_n,
116+% as src/sht/gauss.ts.
117+function [x, w] = gauss_legendre(n)
118+ x = zeros(n, 1);
119+ w = zeros(n, 1);
120+ half = floor((n + 1) / 2);
121+ for i = 1:half
122+ z = cos(pi * (i - 0.25) / (n + 0.5));
123+ pp = 0;
124+ for it = 1:100
125+ p1 = 1;
126+ p2 = 0;
127+ for j = 1:n
128+ p3 = p2;
129+ p2 = p1;
130+ p1 = ((2*j - 1) * z * p2 - (j - 1) * p3) / j;
131+ end
132+ pp = n * (z * p1 - p2) / (z^2 - 1);
133+ dz = p1 / pp;
134+ z = z - dz;
135+ if abs(dz) < 1e-15 * abs(z) + 1e-300
136+ p1 = 1;
137+ p2 = 0;
138+ for j = 1:n
139+ p3 = p2;
140+ p2 = p1;
141+ p1 = ((2*j - 1) * z * p2 - (j - 1) * p3) / j;
142+ end
143+ pp = n * (z * p1 - p2) / (z^2 - 1);
144+ z = z - p1 / pp;
145+ break;
146+ end
147+ end
148+ x(i) = z;
149+ x(n + 1 - i) = -z;
150+ wi = 2 / ((1 - z^2) * pp^2);
151+ w(i) = wi;
152+ w(n + 1 - i) = wi;
153+ end
154+ if mod(n, 2) == 1
155+ x(half) = 0;
156+ end
157+end
158+
159+% Synthesis, spectral -> grid. Grouped calls -- [a, b] = synth(x, y) -- are
160+% the app's batching hint; here each member simply runs in turn.
161+function varargout = synth(varargin)
162+ S = sht_tables();
163+ varargout = cell(1, nargin);
164+ for k = 1:nargin
165+ varargout{k} = synth_one(S, varargin{k});
166+ end
167+end
168+
169+function f = synth_one(S, Q)
170+ % Legendre stage per m, then one inverse FFT per latitude ring with the
171+ % m < 0 modes filled in by conjugate symmetry (the field is real).
172+ G = zeros(S.nphi, S.nlat);
173+ for m = 0:S.mmax
174+ Fm = (S.Y{m + 1} * Q(S.off(m + 1):S.off(m + 1) + S.lmax - m)).';
175+ G(m + 1, :) = Fm;
176+ if m > 0
177+ G(S.nphi + 1 - m, :) = conj(Fm);
178+ end
179+ end
180+ f = S.nphi * real(ifft(G, [], 1));
181+ f = f(:);
182+end
183+
184+% Analysis, grid -> spectral: forward FFT per ring, then Gauss quadrature
185+% against the same Legendre tables.
186+function varargout = analys(varargin)
187+ S = sht_tables();
188+ varargout = cell(1, nargin);
189+ for k = 1:nargin
190+ varargout{k} = analys_one(S, varargin{k});
191+ end
192+end
193+
194+function Q = analys_one(S, f)
195+ F = fft(reshape(f, S.nphi, S.nlat), [], 1) * (2*pi/S.nphi);
196+ Q = complex(zeros(S.nlm, 1));
197+ for m = 0:S.mmax
198+ Q(S.off(m + 1):S.off(m + 1) + S.lmax - m) = S.Y{m + 1}.' * (S.wg .* F(m + 1, :).');
199+ end
200+end
201+
202+% The coefficients of sin(theta)*dtheta(u): the alpha^+/alpha^- shift by one
203+% degree within each m block.
204+function V = dthetac(Q)
205+ S = sht_tables();
206+ V = S.ap .* Q(S.iprev) + S.am .* Q(S.inext);
207+end
208+
209+% The coefficients of dphi(u): i*m, diagonal.
210+function V = dphic(Q)
211+ S = sht_tables();
212+ V = 1i * (S.mv .* Q);
213+end
214+
215+% Grid-space derivatives, coefficients in: compositions of the shuffles and
216+% the synthesis (src/sht/deriv.ts). dtheta divides by sin(theta) afterwards.
217+function f = dtheta(Q)
218+ S = sht_tables();
219+ f = synth(dthetac(Q)) ./ S.stpt;
220+end
221+
222+function f = dphi(Q)
223+ f = synth(dphic(Q));
224+end
225+
226+% Grid-space phi derivative, grid in: two FFT stages and a pointwise i*m,
227+% no Legendre work -- d/dphi is diagonal in the Fourier index.
228+function g = dphig(f)
229+ S = sht_tables();
230+ F = fft(reshape(f, S.nphi, S.nlat), [], 1);
231+ g = real(ifft(S.dmul .* F, [], 1));
232+ g = g(:);
233+end
234+
235+% ---------------------------------------------------------------- the surface
236+%
237+% What the app precomputes from a shape's raw grid values: the band-limited
238+% embedding and both metric formulations built on it (src/geom/geometry.ts,
239+% src/geom/metric.ts). The solver runs on the synthesis of the coefficients,
240+% not on the raw values -- for a shape with sharp features the two differ.
241+function G = surface_tables(gxr, gyr, gzr)
242+ S = sht_tables();
243+ [G.Gx, G.Gy, G.Gz] = analys(gxr, gyr, gzr);
244+ [G.gx, G.gy, G.gz] = synth(G.Gx, G.Gy, G.Gz);
245+
246+ % Flux-form metric weights, from the sin-weighted theta tangent
247+ % sin(theta)*X_theta and X_phi, both smooth on the sphere:
248+ % gtt~ = sin^2 g_tt, gtp~ = sin g_tp, D = J sin^2(theta).
249+ [sXtx, sXty, sXtz] = synth(dthetac(G.Gx), dthetac(G.Gy), dthetac(G.Gz));
250+ [Xpx, Xpy, Xpz] = synth(dphic(G.Gx), dphic(G.Gy), dphic(G.Gz));
251+ gtt = sXtx.^2 + sXty.^2 + sXtz.^2;
252+ gtp = sXtx.*Xpx + sXty.*Xpy + sXtz.*Xpz;
253+ gpp = Xpx.^2 + Xpy.^2 + Xpz.^2;
254+ D = sqrt(gtt .* gpp - gtp.^2);
255+ G.p1 = gpp ./ D;
256+ G.p2 = -gtp ./ D;
257+ G.q2 = gtt ./ D;
258+ G.r = 1 ./ D;
259+
260+ % The sphere-subtracted weights and the bounded 1/J = r sin^2(theta) --
261+ % what keeps the concentrated division off the round sphere's share of the
262+ % flux divergence. Formed here in f64, as the app forms them.
263+ G.jinv = G.r .* S.stpt.^2;
264+ G.dp1 = G.p1 - 1;
265+ G.dq2 = G.q2 - 1;
266+
267+ % Preconditioner scale Jhat = 2/(muMin + muMax) over the eigenvalues of
268+ % the operator's symbol S = (1/J)[[p1, p2], [p2, q2]].
269+ s11 = G.p1 .* G.jinv;
270+ s12 = G.p2 .* G.jinv;
271+ s22 = G.q2 .* G.jinv;
272+ mn = (s11 + s22) / 2;
273+ disc = sqrt(((s11 - s22) / 2).^2 + s12.^2);
274+ G.Jhat = 2 / (min(mn - disc) + max(mn + disc));
275+
276+ % Inverse metric quantities V_theta/V_phi, for the Algorithm-4 models.
277+ Xtx = sXtx ./ S.stpt;
278+ Xty = sXty ./ S.stpt;
279+ Xtz = sXtz ./ S.stpt;
280+ g11 = Xtx.^2 + Xty.^2 + Xtz.^2;
281+ g12 = Xtx.*Xpx + Xty.*Xpy + Xtz.*Xpz;
282+ g22 = gpp;
283+ det = g11 .* g22 - g12.^2;
284+ G.Vtx = (g22 .* Xtx - g12 .* Xpx) ./ det;
285+ G.Vty = (g22 .* Xty - g12 .* Xpy) ./ det;
286+ G.Vtz = (g22 .* Xtz - g12 .* Xpz) ./ det;
287+ G.Vpx = (g11 .* Xpx - g12 .* Xtx) ./ det;
288+ G.Vpy = (g11 .* Xpy - g12 .* Xty) ./ det;
289+ G.Vpz = (g11 .* Xpz - g12 .* Xtz) ./ det;
290+end
291+
292+% ---------------------------------------------------------------- random field
293+%
294+% chebfun-style smooth random field in 3D, restricted to the surface by
295+% evaluating it at the grid points -- the way surfacefun seeds a run. Two
296+% signatures, as in the app:
297+% [k, c] = randnfun3(lambda, dom) the Fourier-mode draw (tools/randnfun3.m)
298+% f = randnfun3(lambda, gx, gy, gz) that draw, summed at the surface points
299+% Seed with rng(...) before calling.
300+function varargout = randnfun3(lambda, varargin)
301+ if nargin == 2
302+ [k, c] = randnfun3_modes(lambda, varargin{1});
303+ varargout = {k, c};
304+ return;
305+ end
306+ [gx, gy, gz] = deal(varargin{1:3});
307+ dom = [min(gx) max(gx) min(gy) max(gy) min(gz) max(gz)];
308+ [k, c] = randnfun3_modes(lambda, dom);
309+ % Summed in blocks of modes: the full npts x nmodes phase matrix can reach
310+ % hundreds of MB at a fine wavelength.
311+ f = zeros(numel(gx), 1);
312+ blk = 2048;
313+ for j0 = 1:blk:size(k, 1)
314+ j1 = min(j0 + blk - 1, size(k, 1));
315+ t = gx * k(j0:j1, 1)' + gy * k(j0:j1, 2)' + gz * k(j0:j1, 3)';
316+ f = f + cos(t) * c(j0:j1, 1) - sin(t) * c(j0:j1, 2);
317+ end
318+ varargout = {f};
319+end
320+
321+% ---------------------------------------------------------------- display
322+%
323+% The pattern on the surface, one panel per species. The solver grid has no
324+% pole rows and an open phi seam; wrap_grid closes both for display, capping
325+% each pole with the mean of its nearest ring.
326+function h = plot_setup(gx, gy, gz, fields, names)
327+ fig = figure('Name', 'turing-surface', 'Color', 'w');
328+ Xs = wrap_grid(gx);
329+ Ys = wrap_grid(gy);
330+ Zs = wrap_grid(gz);
331+ n = numel(fields);
332+ h.surf = gobjects(1, n);
333+ h.ax = gobjects(1, n);
334+ for k = 1:n
335+ ax = subplot(1, n, k, 'Parent', fig);
336+ h.surf(k) = surf(ax, Xs, Ys, Zs, wrap_grid(fields{k}), 'EdgeColor', 'none');
337+ shading(ax, 'interp');
338+ axis(ax, 'equal');
339+ axis(ax, 'off');
340+ colormap(ax, 'jet');
341+ colorbar(ax);
342+ h.ax(k) = ax;
343+ end
344+ h.names = names;
345+end
346+
347+function plot_update(h, fields, t, k, nsteps)
348+ for i = 1:numel(fields)
349+ C = wrap_grid(fields{i});
350+ set(h.surf(i), 'CData', C);
351+ lo = min(C(:));
352+ hi = max(C(:));
353+ if ~(hi > lo)
354+ hi = lo + 1;
355+ end
356+ caxis(h.ax(i), [lo hi]);
357+ title(h.ax(i), sprintf('%s t = %.3f (step %d/%d)', h.names{i}, t, k, nsteps));
358+ end
359+ drawnow;
360+end
361+
362+function M = wrap_grid(f)
363+ S = sht_tables();
364+ M = reshape(f, S.nphi, S.nlat).';
365+ M = [M, M(:, 1)];
366+ M = [mean(M(1, :)) * ones(1, S.nphi + 1); M; mean(M(end, :)) * ones(1, S.nphi + 1)];
367+end
368+
369+% ---------------------------------------------------------------- results file
370+%
371+% Complex coefficients -> flat float32 [re, im] per (l, m), the layout the
372+% app's reference-file reader expects (docs/ellipsoid-reference-spec.md).
373+function write_coeffs(fname, path, Q)
374+ flat = zeros(2 * numel(Q), 1);
375+ flat(1:2:end) = real(Q);
376+ flat(2:2:end) = imag(Q);
377+ h5create(fname, path, numel(flat), 'Datatype', 'single');
378+ h5write(fname, path, single(flat));
379+end
380+
381+% h5writeatt cannot create a bare group, so the attribute-only groups of the
382+% reference layout are made through the low-level API.
383+function make_group(fname, path)
384+ fid = H5F.open(fname, 'H5F_ACC_RDWR', 'H5P_DEFAULT');
385+ gid = H5G.create(fid, path, 'H5P_DEFAULT', 'H5P_DEFAULT', 'H5P_DEFAULT');
386+ H5G.close(gid);
387+ H5F.close(fid);
388+end
src/main.tsmodified+426−41View file
@@ -39,6 +39,9 @@ import {
3939 variantLabel,
4040 type Variant,
4141 } from './compare/variants.ts';
42+import { loadReferenceFile } from './compare/referenceFile.ts';
43+import type { ReferenceCase } from './compare/referenceCase.ts';
44+import { generateMatlabScript, MATLAB_SCRIPT_NAME } from './export/matlabScript.ts';
4245
4346 const $ = <T extends HTMLElement>(id: string): T =>
4447 document.getElementById(id) as T;
@@ -51,6 +54,7 @@ const elLmax = $<HTMLSelectElement>('lmax');
5154 const elOversample = $<HTMLSelectElement>('oversample');
5255 const elColormap = $<HTMLSelectElement>('colormap');
5356 const elRunPause = $<HTMLButtonElement>('runpause');
57+const elRestart = $<HTMLButtonElement>('restart');
5458 const elBenchmark = $<HTMLButtonElement>('benchmark');
5559 const elReseed = $<HTMLButtonElement>('reseed');
5660 const elLam3 = $<HTMLInputElement>('lam3');
@@ -61,12 +65,18 @@ const elMovieSpeed = $<HTMLSelectElement>('moviespeed');
6165 const elMovieRes = $<HTMLSelectElement>('movieres');
6266 const elMovieRotate = $<HTMLInputElement>('movierotate');
6367 const elMovie = $<HTMLButtonElement>('movie');
64-const elCompareToggle = $<HTMLButtonElement>('comparetoggle');
68+const elModeSimulate = $<HTMLButtonElement>('mode-simulate');
69+const elModeEffort = $<HTMLButtonElement>('mode-effort');
70+const elModeVsUpload = $<HTMLButtonElement>('mode-vs-upload');
71+const elModeDesc = $('mode-desc');
6572 const elCompareBar = $('comparebar');
6673 const elCmpNiter = $('cmp-niter');
6774 const elCmpLmax = $('cmp-lmax');
6875 const elCmpDt = $('cmp-dt');
6976 const elCmpRef = $<HTMLSelectElement>('cmp-ref');
77+const elCmpFile = $<HTMLInputElement>('cmp-file');
78+const elCmpFileInfo = $('cmp-fileinfo');
79+const elCmpFileClear = $<HTMLButtonElement>('cmp-fileclear');
7080 const elCmpStart = $<HTMLButtonElement>('cmp-start');
7181 const elCmpCount = $('cmp-count');
7282 const elParams = $('params');
@@ -77,6 +87,10 @@ const elStats = $('stats');
7787 const elBenchResult = $('benchresult');
7888 const elCmd = $('cmd');
7989 const elCopyCmd = $<HTMLButtonElement>('copycmd');
90+const elMatlab = $<HTMLDetailsElement>('matlab');
91+const elMatlabScript = $('matlabscript');
92+const elCopyMatlab = $<HTMLButtonElement>('copymatlab');
93+const elDownloadMatlab = $<HTMLButtonElement>('downloadmatlab');
8094 const elBlurb = $('blurb');
8195 const elErr = $('err');
8296 const elSource = $<HTMLTextAreaElement>('source');
@@ -87,6 +101,21 @@ const elEditorFile = $<HTMLSelectElement>('editor-file');
87101 const elRecompile = $<HTMLButtonElement>('recompile');
88102 const elRevert = $<HTMLButtonElement>('revert');
89103
104+/** The named groups the control area is organized into (index.html's
105+ * `.ctrl-group[data-group]` wrappers). Each mode shows a declared subset of
106+ * these — see MODE_GROUPS and applyModeVisibility below. */
107+const GROUP_NAMES = [
108+ 'surface', 'surface-params', 'solver', 'display',
109+ 'playback', 'benchmark', 'seed', 'movie',
110+] as const;
111+type GroupName = (typeof GROUP_NAMES)[number];
112+const groupEls: Record<GroupName, HTMLElement> = Object.fromEntries(
113+ GROUP_NAMES.map((name) => [
114+ name,
115+ document.querySelector(`.ctrl-group[data-group="${name}"]`) as HTMLElement,
116+ ]),
117+) as Record<GroupName, HTMLElement>;
118+
90119 for (const p of presets) {
91120 const o = document.createElement('option');
92121 o.value = p.key;
@@ -254,6 +283,11 @@ let posBuf: Float32Array | null = null;
254283 * mode. While it is non-null there is no `session`: the study owns one per
255284 * variant, and the panels area is its grid. */
256285 let compareRun: CompareRun | null = null;
286+/** `session`'s spectral state as of the last (re-)seed — what "Restart"
287+ * rewinds to. Captured fresh each time a new field is actually established
288+ * (rebuild/reseed), not just once, so Restart reflects the run's current
289+ * starting point rather than permanently the very first draw. */
290+let initialState: Record<string, Float32Array> | null = null;
257291
258292 const source = (): string => editedSource ?? model.source;
259293 const geomSource = (): string => editedGeomSource ?? geometry.source;
@@ -261,6 +295,10 @@ const geomSource = (): string => editedGeomSource ?? geometry.source;
261295 // ---------------------------------------------------------------- UI wiring
262296 function buildParamInputs(): void {
263297 elParams.replaceChildren();
298+ if (model.params.length === 0) return;
299+ const tag = document.createElement('label');
300+ tag.textContent = 'model parameters';
301+ elParams.append(tag);
264302 for (const spec of model.params) {
265303 const label = document.createElement('label');
266304 label.textContent = `${spec.label} `;
@@ -296,7 +334,7 @@ function buildGeomParamInputs(): void {
296334 elGeomParams.replaceChildren();
297335 if (geometry.params.length === 0) return;
298336 const tag = document.createElement('label');
299- tag.textContent = `${geometry.key}.m`;
337+ tag.textContent = 'geometry parameters';
300338 elGeomParams.append(tag);
301339 for (const spec of geometry.params) {
302340 // A random seed picks a draw and means nothing on its own, so it gets a
@@ -316,6 +354,7 @@ function buildGeomParamInputs(): void {
316354 next = spec.min + Math.floor(Math.random() * (span + 1));
317355 }
318356 geomParams[spec.key] = next;
357+ updateCommand();
319358 viewChange = viewChange.then(() => applyGeometry());
320359 });
321360 elGeomParams.append(button);
@@ -332,6 +371,7 @@ function buildGeomParamInputs(): void {
332371 input.addEventListener('change', () => {
333372 const v = Number(input.value);
334373 if (Number.isFinite(v)) geomParams[spec.key] = v;
374+ updateCommand();
335375 viewChange = viewChange.then(() => applyGeometry());
336376 });
337377 label.append(input);
@@ -366,6 +406,8 @@ function applyGeometryChoice(key: string): void {
366406 editedGeomSource = null;
367407 buildGeomParamInputs();
368408 showEditorFile();
409+ // The command line and the MATLAB export both bake the surface in.
410+ updateCommand();
369411 }
370412
371413 /** Load the chosen file into the editor, keeping any unsaved edit to it. */
@@ -403,7 +445,47 @@ function currentSpec(): RunSpec {
403445 }
404446
405447 function updateCommand(): void {
406- elCmd.textContent = formatCommand(currentSpec());
448+ // A study against a reference file replays the file, so its desktop
449+ // equivalent is the ref checker, not the benchmark.
450+ if (compareRun?.refFile) {
451+ elCmd.textContent = `npm run ref -- --in ${compareRun.refFile.label}`;
452+ } else {
453+ elCmd.textContent = formatCommand(currentSpec());
454+ }
455+ if (elMatlab.open) refreshMatlabScript();
456+}
457+
458+/** The run on screen as one standalone .m — in a study, its reference
459+ * variant, the same choice `currentSpec` makes. Throws on a working copy
460+ * the export cannot parse (no init/step/shape function). */
461+function matlabScriptText(): string {
462+ const spec = currentSpec();
463+ return generateMatlabScript({
464+ model,
465+ modelSource: source(),
466+ params: spec.params,
467+ geometry,
468+ geometrySource: geomSource(),
469+ geometryParams: geomParams,
470+ lmax: spec.lmax,
471+ niter: spec.niter,
472+ lam3: Number(elLam3.value),
473+ seed,
474+ preset: spec.preset,
475+ command: formatCommand(spec),
476+ });
477+}
478+
479+/** Regenerate the visible script text; on failure, show why in its place. */
480+function refreshMatlabScript(): string | null {
481+ try {
482+ const text = matlabScriptText();
483+ elMatlabScript.textContent = text;
484+ return text;
485+ } catch (e) {
486+ elMatlabScript.textContent = e instanceof Error ? e.message : String(e);
487+ return null;
488+ }
407489 }
408490
409491 elModel.addEventListener('change', () => {
@@ -436,6 +518,8 @@ elGeometry.addEventListener('change', () => {
436518 elLam3.addEventListener('change', () => {
437519 const v = Number(elLam3.value);
438520 if (!Number.isFinite(v) || v <= 0) return;
521+ // Not in the bench command, but the MATLAB export bakes lam3 in.
522+ updateCommand();
439523 // Changing the wavelength redraws the field, which restarts the run — so
440524 // pause first, exactly as the Re-seed button does. Without it the reseed's
441525 // readback races the pump's own, and the two collide on the staging buffer.
@@ -491,6 +575,10 @@ elReseed.addEventListener('click', () => {
491575 updateCommand();
492576 void reseed();
493577 });
578+elRestart.addEventListener('click', () => {
579+ setRunning(false);
580+ void restart();
581+});
494582 elResetView.addEventListener('click', () => {
495583 compareRun?.resetView();
496584 for (const s of scenes) s.resetCamera();
@@ -535,6 +623,46 @@ elCopyCmd.addEventListener('click', () => {
535623 navigator.clipboard.writeText(text).then(() => flash('Copied'), selectCommand);
536624 });
537625
626+// The MATLAB export: the same run as one self-contained .m. Regenerated from
627+// the current UI state whenever it is shown, copied or downloaded, so the
628+// text always matches the run on screen.
629+elMatlab.addEventListener('toggle', () => {
630+ if (elMatlab.open) refreshMatlabScript();
631+});
632+elCopyMatlab.addEventListener('click', (e) => {
633+ // The buttons live inside the <summary>; without this a click also toggles.
634+ e.preventDefault();
635+ e.stopPropagation();
636+ const text = refreshMatlabScript();
637+ if (text === null || !navigator.clipboard) {
638+ // Open to show the error, or the text to select by hand.
639+ elMatlab.open = true;
640+ return;
641+ }
642+ navigator.clipboard.writeText(text).then(
643+ () => {
644+ elCopyMatlab.textContent = 'Copied';
645+ setTimeout(() => (elCopyMatlab.textContent = 'Copy'), 1200);
646+ },
647+ () => (elMatlab.open = true),
648+ );
649+});
650+elDownloadMatlab.addEventListener('click', (e) => {
651+ e.preventDefault();
652+ e.stopPropagation();
653+ const text = refreshMatlabScript();
654+ if (text === null) {
655+ elMatlab.open = true;
656+ return;
657+ }
658+ const url = URL.createObjectURL(new Blob([text], { type: 'text/x-matlab' }));
659+ const a = document.createElement('a');
660+ a.href = url;
661+ a.download = `${MATLAB_SCRIPT_NAME}.m`;
662+ a.click();
663+ URL.revokeObjectURL(url);
664+});
665+
538666 // ---------------------------------------------------------------- setup
539667 function disposeView(): void {
540668 for (const s of scenes) s.dispose();
@@ -751,6 +879,9 @@ async function rebuild(): Promise<void> {
751879 if (gen !== generation) return;
752880
753881 await session.seed(seed);
882+ if (gen !== generation) return;
883+ initialState = await session.readState();
884+ if (gen !== generation) return;
754885
755886 const plan = session.describe();
756887 elCompiled.textContent =
@@ -785,6 +916,24 @@ async function reseed(): Promise<void> {
785916 const gen = generation;
786917 await session.seed(seed);
787918 if (gen !== generation) return;
919+ initialState = await session.readState();
920+ if (gen !== generation) return;
921+ for (const r of ranges) {
922+ r.lo = NaN;
923+ r.hi = NaN;
924+ }
925+ await draw();
926+ updateStats();
927+}
928+
929+/** Rewind to the field this run is currently starting from — the last
930+ * (re-)seed, not necessarily the very first one — without drawing a new
931+ * one. Unlike reseed(), the seed value and lam3 are untouched, so nothing
932+ * the CLI command line encodes changes. */
933+async function restart(): Promise<void> {
934+ if (compareRun) return compareRun.restart();
935+ if (!session || !initialState) return;
936+ session.loadState(initialState);
788937 for (const r of ranges) {
789938 r.lo = NaN;
790939 r.hi = NaN;
@@ -986,7 +1135,7 @@ function submitSteps(n: number): void {
9861135 function setMovieUi(on: boolean): void {
9871136 const locked = [
9881137 elModel, elGeometry, elMorph, elNiter, elLmax, elOversample, elColormap,
989- elRunPause, elBenchmark, elReseed, elRecompile, elRevert, elEditorFile,
1138+ elRunPause, elRestart, elBenchmark, elReseed, elRecompile, elRevert, elEditorFile,
9901139 elMovieSpeed, elMovieRes, elMovieRotate, elMovieToggle,
9911140 ];
9921141 for (const el of locked) el.disabled = on;
@@ -1203,6 +1352,15 @@ function buildChips(host: HTMLElement, values: number[], selected: Set<number>,
12031352 const cmpVariants = (): Variant[] =>
12041353 crossProduct([...cmpSelected.niter], [...cmpSelected.lmax], [...cmpSelected.dt]);
12051354
1355+/**
1356+ * A loaded reference file, or null. While one is loaded the study checks the
1357+ * variants against it instead of against each other: the file defines the
1358+ * whole problem (model, parameters, geometry, initial state, end time), so
1359+ * the page's own model and geometry choices do not enter the study at all —
1360+ * only the solver knobs above do.
1361+ */
1362+let refCase: ReferenceCase | null = null;
1363+
12061364 /** The reference the user picked, clamped to the current variant list. */
12071365 let cmpRefKey = '';
12081366
@@ -1215,19 +1373,32 @@ function compareRefIndex(): number {
12151373 function refreshVariants(): void {
12161374 const variants = cmpVariants();
12171375 const showDt = cmpSelected.dt.size > 1;
1218- const panels = variants.length * model.species.length;
1376+ // With a file loaded the study's model is the file's, and its final state
1377+ // is one more row of panels.
1378+ const cmpModel = refCase?.model ?? model;
1379+ const rowCount = variants.length + (refCase ? 1 : 0);
1380+ const panels = rowCount * cmpModel.species.length;
12191381
12201382 const prev = cmpRefKey;
12211383 elCmpRef.replaceChildren();
1222- for (const v of variants) {
1384+ if (refCase) {
1385+ // The file is the reference; the pick among variants means nothing here.
12231386 const o = document.createElement('option');
1224- o.value = variantKey(v);
1225- o.textContent = variantLabel(v, showDt);
1387+ o.textContent = `the file's final state`;
12261388 elCmpRef.append(o);
1389+ elCmpRef.disabled = true;
1390+ } else {
1391+ elCmpRef.disabled = false;
1392+ for (const v of variants) {
1393+ const o = document.createElement('option');
1394+ o.value = variantKey(v);
1395+ o.textContent = variantLabel(v, showDt);
1396+ elCmpRef.append(o);
1397+ }
1398+ const keys = variants.map(variantKey);
1399+ cmpRefKey = keys.includes(prev) ? prev : keys[mostResolved(variants)];
1400+ elCmpRef.value = cmpRefKey;
12271401 }
1228- const keys = variants.map(variantKey);
1229- cmpRefKey = keys.includes(prev) ? prev : keys[mostResolved(variants)];
1230- elCmpRef.value = cmpRefKey;
12311402
12321403 const tooMany =
12331404 variants.length > MAX_VARIANTS
@@ -1237,23 +1408,53 @@ function refreshVariants(): void {
12371408 : '';
12381409 elCmpCount.textContent = tooMany
12391410 ? `too many: ${tooMany}`
1240- : `${variants.length} variants × ${model.species.length} species = ${panels} panels`;
1411+ : `${variants.length} variant${variants.length === 1 ? '' : 's'}` +
1412+ `${refCase ? ' + the file' : ''} × ` +
1413+ `${cmpModel.species.length} species = ${panels} panels`;
12411414 elCmpCount.style.color = tooMany ? '#b35900' : '';
12421415 elCmpStart.disabled = tooMany !== '' && compareRun === null;
12431416 }
12441417
1245-buildChips(
1246- elCmpNiter,
1247- [...elNiter.options].map((o) => Number(o.value)),
1248- cmpSelected.niter,
1249- String,
1250-);
1251-buildChips(
1252- elCmpLmax,
1253- [...elLmax.options].map((o) => Number(o.value)),
1254- cmpSelected.lmax,
1255- String,
1256-);
1418+/**
1419+ * The niter chips on offer. A loaded reference file adds its own recorded
1420+ * iteration count if the standard list lacks it, so the file's settings are
1421+ * always selectable; clearing the file drops any selection outside the
1422+ * standard list again.
1423+ */
1424+function rebuildNiterChips(): void {
1425+ const all = [...elNiter.options].map((o) => Number(o.value));
1426+ let values = all;
1427+ if (refCase && !all.includes(refCase.niter)) {
1428+ values = [...all, refCase.niter].sort((a, b) => a - b);
1429+ }
1430+ if (!refCase) {
1431+ for (const v of [...cmpSelected.niter]) if (!values.includes(v)) cmpSelected.niter.delete(v);
1432+ if (cmpSelected.niter.size === 0) cmpSelected.niter.add(DEFAULT_NITER);
1433+ }
1434+ buildChips(elCmpNiter, values, cmpSelected.niter, String);
1435+}
1436+
1437+/**
1438+ * The lmax chips on offer. A loaded reference file floors them at its own
1439+ * band: a variant below it could not even hold the file's initial state
1440+ * (prolongation only widens), so those values are not offered rather than
1441+ * offered and refused.
1442+ */
1443+function rebuildLmaxChips(): void {
1444+ const all = [...elLmax.options].map((o) => Number(o.value));
1445+ let values = all;
1446+ if (refCase) {
1447+ const floor = refCase.lmax;
1448+ values = all.filter((v) => v >= floor);
1449+ if (!values.includes(floor)) values = [floor, ...values];
1450+ for (const v of [...cmpSelected.lmax]) if (!values.includes(v)) cmpSelected.lmax.delete(v);
1451+ if (cmpSelected.lmax.size === 0) cmpSelected.lmax.add(floor);
1452+ }
1453+ buildChips(elCmpLmax, values, cmpSelected.lmax, String);
1454+}
1455+
1456+rebuildNiterChips();
1457+rebuildLmaxChips();
12571458 buildChips(elCmpDt, DT_DIVISORS, cmpSelected.dt, (v) => (v === 1 ? 'dt' : `dt/${v}`));
12581459 refreshVariants();
12591460
@@ -1262,8 +1463,171 @@ elCmpRef.addEventListener('change', () => {
12621463 if (compareRun) void rebuildCompare();
12631464 });
12641465
1265-elCompareToggle.addEventListener('click', () => {
1266- elCompareBar.hidden = !elCompareBar.hidden;
1466+/** Reflect the loaded (or cleared) reference file in the compare bar. */
1467+function applyRefUi(): void {
1468+ elCmpFileInfo.hidden = elCmpFileClear.hidden = refCase === null;
1469+ if (refCase) {
1470+ const rc = refCase;
1471+ const geomParamText = rc.geometry.params
1472+ .map((p) => `${p.key}=${rc.geometryParams[p.key]}`)
1473+ .join(' ');
1474+ const name = document.createElement('b');
1475+ name.textContent = rc.label;
1476+ const info = document.createElement('span');
1477+ info.textContent =
1478+ ` — ${rc.model.label} on ${rc.geometry.label.toLowerCase()}` +
1479+ (geomParamText ? ` (${geomParamText})` : '') +
1480+ `, lmax ${rc.lmax}, T = ${(rc.steps * (rc.params.dt ?? 0)).toFixed(2)}` +
1481+ ` (${rc.steps} × dt ${rc.params.dt})`;
1482+ elCmpFileInfo.replaceChildren(name, info);
1483+ }
1484+ rebuildNiterChips();
1485+ rebuildLmaxChips();
1486+ refreshVariants();
1487+}
1488+
1489+/**
1490+ * The four top-level modes and which control groups each shows (see
1491+ * GROUP_NAMES/groupEls above; `.ctrl-group` wrappers in index.html).
1492+ * `currentMode` tracks which configuration is on screen — the compare bar
1493+ * being open, and in which flavor — not whether a study has actually been
1494+ * started inside it. That match matters: without it, opening the bar
1495+ * (which already shows the right groups) leaves its top-row button
1496+ * unhighlighted until a study happens to start, which is inconsistent with
1497+ * `vs-upload`'s one-click flow and reads as broken.
1498+ */
1499+type Mode = 'simulate' | 'compute-effort' | 'vs-sphere' | 'vs-upload';
1500+let currentMode: Mode = 'simulate';
1501+
1502+const MODE_GROUPS: Record<Mode, readonly GroupName[]> = {
1503+ simulate: ['surface', 'surface-params', 'solver', 'display', 'playback', 'benchmark', 'seed', 'movie'],
1504+ 'compute-effort': ['surface', 'surface-params', 'display', 'playback', 'seed'],
1505+ 'vs-sphere': [], // unreachable — the button is disabled, no listener ever calls setMode with this
1506+ // No `seed` here: nothing in that group does anything useful against a
1507+ // loaded file (lam3 is silently absorbed, and Restart already covers what
1508+ // Re-seed would otherwise be doing — reloading the file's fixed initial
1509+ // state) — see CompareRun.restart().
1510+ 'vs-upload': ['display', 'playback'],
1511+};
1512+
1513+const MODE_DESCRIPTIONS: Record<Mode, string> = {
1514+ simulate:
1515+ 'This mode runs one standalone reaction-diffusion solver.',
1516+ 'compute-effort':
1517+ 'When we change the computational effort of the solver by varying solve iterations, lmax, or timestep, ' +
1518+ 'how does the solution change? Find out by running several ' +
1519+ 'so you can see how each setting trades accuracy for speed.',
1520+ 'vs-sphere': '',
1521+ 'vs-upload':
1522+ 'Load a saved reference run (an .h5 file) and run this solver to the ' +
1523+ 'same physical end time from the same initial condition, to check how ' +
1524+ 'closely it reproduces the reference. You can adjust the solver settings ' +
1525+ 'to see how they affect the outcome.',
1526+};
1527+
1528+function setModeButtons(mode: Mode): void {
1529+ elModeSimulate.setAttribute('aria-pressed', String(mode === 'simulate'));
1530+ elModeEffort.setAttribute('aria-pressed', String(mode === 'compute-effort'));
1531+ elModeVsUpload.setAttribute('aria-pressed', String(mode === 'vs-upload'));
1532+ elModeDesc.textContent = MODE_DESCRIPTIONS[mode];
1533+}
1534+
1535+/** Show exactly the groups `mode` declares; hide the rest. */
1536+function applyModeVisibility(mode: Mode): void {
1537+ currentMode = mode;
1538+ const shown = new Set<GroupName>(MODE_GROUPS[mode]);
1539+ for (const name of GROUP_NAMES) groupEls[name].hidden = !shown.has(name);
1540+ setModeButtons(mode);
1541+}
1542+
1543+/**
1544+ * Enter `mode`: groups, top-row buttons, and the compare bar's own
1545+ * visibility (open for the two compare flavors, closed for Simulate).
1546+ * Doesn't touch `compareRun`/`refCase` or start/stop a study — callers
1547+ * decide that; this only decides what's on screen, and it decides it
1548+ * immediately, so the button you clicked lights up right away rather than
1549+ * waiting on a study that may not exist yet (or may never start, if the
1550+ * bar's own Compare is never pressed).
1551+ */
1552+function enterMode(mode: Mode): void {
1553+ applyModeVisibility(mode);
1554+ elCompareBar.hidden = mode === 'simulate';
1555+}
1556+
1557+/** Entering a mode from the top row. */
1558+function setMode(mode: Mode): void {
1559+ if (mode === 'vs-sphere') return; // unreachable — button is disabled
1560+ if (mode === 'simulate') {
1561+ if (compareRun) void stopCompare();
1562+ enterMode('simulate');
1563+ return;
1564+ }
1565+ if (mode === 'compute-effort') {
1566+ // Tear down whatever study is running first (mirrors Simulate above) —
1567+ // stopCompare's synchronous prefix disposes it and nulls `compareRun`
1568+ // before its first `await`, so `refCase` is safe to drop right after.
1569+ if (compareRun) void stopCompare();
1570+ if (refCase) {
1571+ refCase = null;
1572+ applyRefUi();
1573+ }
1574+ enterMode('compute-effort');
1575+ return;
1576+ }
1577+ // vs-upload: opens the file picker; entering the mode itself happens once
1578+ // a file is actually chosen (elCmpFile's change handler below) — not here,
1579+ // since cancelling the dialog must leave the current mode untouched.
1580+ elCmpFile.click();
1581+}
1582+
1583+elModeSimulate.addEventListener('click', () => setMode('simulate'));
1584+elModeEffort.addEventListener('click', () => setMode('compute-effort'));
1585+elModeVsUpload.addEventListener('click', () => setMode('vs-upload'));
1586+
1587+elCmpFile.addEventListener('change', () => {
1588+ const file = elCmpFile.files?.[0];
1589+ // Cleared so picking the same file again still fires a change event.
1590+ elCmpFile.value = '';
1591+ if (!file) return;
1592+ void (async () => {
1593+ try {
1594+ refCase = await loadReferenceFile(file);
1595+ elErr.textContent = '';
1596+ } catch (e) {
1597+ refCase = null;
1598+ elErr.textContent = `reference file ${file.name}: ${e instanceof Error ? e.message : e}`;
1599+ applyRefUi();
1600+ return;
1601+ }
1602+ // One click, one study: the file's own settings become the single
1603+ // variant — its recorded niter, its band, its dt undivided — and the
1604+ // comparison opens on them, paused at the initial state so what runs is
1605+ // the user's choice. (Widening it is: teardown the comparison, pick more
1606+ // chips, compile it again — the file stays loaded.)
1607+ cmpSelected.niter.clear();
1608+ cmpSelected.niter.add(refCase.niter);
1609+ cmpSelected.lmax.clear();
1610+ cmpSelected.lmax.add(refCase.lmax);
1611+ cmpSelected.dt.clear();
1612+ cmpSelected.dt.add(1);
1613+ applyRefUi();
1614+ enterMode('vs-upload');
1615+ if (compareRun) {
1616+ // A study is already up (this one loaded over it): same teardown as
1617+ // rebuildCompare, then the new file's study takes its place.
1618+ compareRun.dispose();
1619+ compareRun = null;
1620+ setCompareUi(false);
1621+ }
1622+ await startCompare();
1623+ })();
1624+});
1625+elCmpFileClear.addEventListener('click', () => {
1626+ refCase = null;
1627+ applyRefUi();
1628+ // The bar stays open — this only drops back to the plain chip comparison.
1629+ // Only reachable while idle (elCmpFileClear is disabled during a study).
1630+ enterMode('compute-effort');
12671631 });
12681632
12691633 elCmpStart.addEventListener('click', () => {
@@ -1271,23 +1635,38 @@ elCmpStart.addEventListener('click', () => {
12711635 else void startCompare();
12721636 });
12731637
1274-/** Controls the study supersedes or cannot honour while it is running. */
1638+/**
1639+ * Controls the study supersedes or cannot honour while it is running.
1640+ * Mode/group/button state is not this function's job — that's set the
1641+ * moment a mode is entered (enterMode, above), independent of whether a
1642+ * study inside it has actually started or stopped.
1643+ */
12751644 function setCompareUi(on: boolean): void {
1276- for (const el of [elNiter, elLmax, elOversample, elBenchmark, elMovieToggle]) {
1277- el.disabled = on;
1278- }
1645+ // A study picks its own display grid, so oversample stays individually
1646+ // disabled inside the still-visible display group; and clearing a loaded
1647+ // file out from under a running study would leave it checking against one
1648+ // that no longer exists.
1649+ elOversample.disabled = on;
1650+ elCmpFileClear.disabled = on;
12791651 elCmpNiter.querySelectorAll('button').forEach((b) => (b.disabled = on));
12801652 elCmpLmax.querySelectorAll('button').forEach((b) => (b.disabled = on));
12811653 elCmpDt.querySelectorAll('button').forEach((b) => (b.disabled = on));
1282- elCmpStart.textContent = on ? 'Stop comparing' : 'Compare';
1283- elCompareToggle.textContent = on ? 'Comparing' : 'Compare';
1654+ elCmpStart.textContent = on ? 'Teardown comparison' : 'Compile comparison';
1655+ // The movie bar's own hidden flag is independent of the movie *group's* —
1656+ // force it closed so it doesn't reappear open once the group is shown
1657+ // again on returning to Simulate.
12841658 if (on) elMovieBar.hidden = true;
12851659 }
12861660
12871661 async function startCompare(): Promise<void> {
12881662 if (compareRun || !device) return;
1663+ // Snapshotted for the whole study: `refCase` only changes with no study up
1664+ // (clearing is disabled during one, and loading tears it down first).
1665+ const rc = refCase;
1666+ const cmpModel = rc?.model ?? model;
12891667 const variants = cmpVariants();
1290- if (variants.length > MAX_VARIANTS || variants.length * model.species.length > MAX_PANELS) {
1668+ const rowCount = variants.length + (rc ? 1 : 0);
1669+ if (variants.length > MAX_VARIANTS || rowCount * cmpModel.species.length > MAX_PANELS) {
12911670 return;
12921671 }
12931672 // Take down the single run first: its pump, its scenes, its session. The
@@ -1303,18 +1682,23 @@ async function startCompare(): Promise<void> {
13031682 setCompareUi(true);
13041683
13051684 try {
1685+ // Against a reference file, the problem is the file's — its model,
1686+ // parameters and geometry, from the registry sources (the editor's
1687+ // working copies describe the page's run, not the file's).
13061688 compareRun = await CompareRun.create({
13071689 device,
1308- model,
1309- params,
1310- source: source(),
1311- geometry,
1312- geometryParams: geomParams,
1313- geometrySource: geomSource(),
1690+ model: cmpModel,
1691+ params: rc ? rc.params : params,
1692+ source: rc ? rc.model.source : source(),
1693+ geometry: rc ? rc.geometry : geometry,
1694+ geometryParams: rc ? rc.geometryParams : geomParams,
1695+ geometrySource: rc ? rc.geometry.source : geomSource(),
13141696 variants,
1315- reference: compareRefIndex(),
1697+ reference: rc ? 0 : compareRefIndex(),
1698+ refFile: rc ?? undefined,
1699+ onFinished: () => setRunning(false),
13161700 seed,
1317- lam3: Number(elLam3.value),
1701+ lam3: rc ? undefined : Number(elLam3.value),
13181702 morph,
13191703 colormapName: () => elColormap.value,
13201704 container: elPanels,
@@ -1356,6 +1740,7 @@ async function rebuildCompare(): Promise<void> {
13561740
13571741 // ---------------------------------------------------------------- boot
13581742 async function boot(): Promise<void> {
1743+ enterMode('simulate');
13591744 elModel.value = presets[0].key;
13601745 // The iteration count is one default shared with the benchmark, like the
13611746 // rest of the RunSpec's — take it from there rather than from the markup, so
src/mgpu/session.tsmodified+21−7View file
@@ -398,18 +398,32 @@ export class ModelSession {
398398 }
399399
400400 /**
401- * Read species `k` at render resolution (`viewSht`'s grid). Without
402- * oversampling this is the grid field the .m returned. With oversampling the
403- * spectral state is synthesized on the finer grid instead — the same field,
404- * since the models define each species as synth of its state, evaluated
405- * exactly on more points.
401+ * Read every state species back to the CPU, in the shape `loadState`
402+ * consumes — the capture side of that method's reload, so a caller can
403+ * hold onto the current spectral state and restore it exactly later
404+ * (e.g. a "restart to this run's initial condition" control). One at a
405+ * time, not `Promise.all`: every `read()` copies into the same shared
406+ * readback buffer (`GpuModel#readback`, model.ts:448-452), so two in
407+ * flight at once race `mapAsync` against each other's `unmap`.
408+ */
409+ async readState(): Promise<Record<string, Float32Array>> {
410+ const out: Record<string, Float32Array> = {};
411+ for (const name of this.model.state) out[name] = await this.read(name);
412+ return out;
413+ }
414+
415+ /**
416+ * Read species `k` at render resolution (`viewSht`'s grid): the spectral
417+ * state synthesized there. The models define each species as synth of its
418+ * state, so this is the field the .m returned — evaluated exactly, whatever
419+ * the grid — and it is current however the state last changed, including a
420+ * `loadState`, which runs no kernel that would write the grid-space fields.
406421 */
407422 readSpecies(k: number): Promise<Float32Array> {
408- if (!this.#displaySht) return this.read(this.model.species[k]);
409423 const state = this.model.state[k];
410424 const buf = this.gpu.valueBuffer(state);
411425 if (!buf) throw new Error(`readSpecies: no buffer for state '${state}'`);
412- return this.#displaySht.synthFrom(buf);
426+ return this.viewSht.synthFrom(buf);
413427 }
414428
415429 describe(): { init: string[]; step: string[] } {
test/compareChecks.tsmodified+67−0View file
@@ -75,6 +75,73 @@ export async function compareChecks(
7575 );
7676 }
7777
78+ // ---- a file's exact state loads onto every grid --------------------------
79+ // What a reference-file study does instead of seeding: the file's spectral
80+ // state pushed into each variant by loadState, prolonged into its band. The
81+ // load is a plain upload, so the state must come back bit-exact; and read on
82+ // one shared grid the variants must then show one field, because synthesis
83+ // of the same band-limited coefficients is evaluation, not resampling.
84+ {
85+ const model = mModelByKey('allencahn')!;
86+ const params = defaultParams(model);
87+ const sessions: ModelSession[] = [];
88+ try {
89+ for (const lmax of [COARSE, FINE]) {
90+ sessions.push(await ModelSession.create({ device, model, params, lmax, niter: 0 }));
91+ }
92+ const [coarse, fine] = sessions;
93+ // A deterministic band-limited state, decaying like a real spectrum;
94+ // m = 0 imaginary parts stay zero (the state is a real field).
95+ const q = new Float32Array(2 * nlmCalc(COARSE, COARSE));
96+ for (let m = 0; m <= COARSE; m++) {
97+ for (let l = m; l <= COARSE; l++) {
98+ const i = 2 * lmIndex(COARSE, l, m);
99+ const amp = Math.exp(-l / 6);
100+ q[i] = amp * Math.sin(1 + 3 * l + 7 * m);
101+ q[i + 1] = m === 0 ? 0 : amp * Math.cos(2 + 5 * l + 11 * m);
102+ }
103+ }
104+ coarse.loadState({ U: q });
105+ fine.loadState({ U: prolongCoeffs(q, COARSE, FINE) });
106+
107+ const back = await coarse.read('U');
108+ let exact = back.length === q.length;
109+ if (exact) {
110+ for (let i = 0; i < q.length; i++) {
111+ if (back[i] !== q[i]) {
112+ exact = false;
113+ break;
114+ }
115+ }
116+ }
117+ check(
118+ 'compare: loadState puts the exact coefficients in the state',
119+ exact,
120+ `${q.length} float32 values round-tripped bit-exact at lmax ${COARSE}`,
121+ );
122+
123+ // The coarse session's own solver grid, so its display plan is the
124+ // solver's — the branch a crowded study lands on.
125+ for (const s of sessions) await s.setDisplayGrid(64, 128);
126+ const cu = await coarse.readSpecies(0);
127+ const fu = await fine.readSpecies(0);
128+ let maxd = 0;
129+ let scale = 0;
130+ for (let i = 0; i < cu.length; i++) {
131+ maxd = Math.max(maxd, Math.abs(cu[i] - fu[i]));
132+ scale = Math.max(scale, Math.abs(cu[i]));
133+ }
134+ check(
135+ 'compare: one loaded state reads back as one field on a shared grid',
136+ maxd < 1e-4 * scale,
137+ `max |du| = ${maxd.toExponential(2)} vs max |u| = ${scale.toExponential(2)} ` +
138+ `across lmax ${COARSE} vs ${FINE}`,
139+ );
140+ } finally {
141+ for (const s of sessions) s.destroy();
142+ }
143+ }
144+
78145 // ---- one random field across lmax: the shipped models' seeding -----------
79146 {
80147 const model = mModelByKey('schnakenberg')!;
test/matlabExportChecks.tsadded+134−0View file
@@ -0,0 +1,134 @@
1+/**
2+ * The MATLAB export (src/export/matlabScript.ts) is string assembly, so these
3+ * checks are cheap and need no GPU: every preset x geometry combination must
4+ * generate, the assembled file must keep its local-function namespace free of
5+ * collisions, and the driver's calls must match the signatures the .m files
6+ * declare. Whether the generated MATLAB actually reproduces a run is checked
7+ * against MATLAB itself, not here: a run exported at defaults and executed in
8+ * MATLAB R2026b lands within fp32 accumulation error of the app's own replay
9+ * (relL2 ~1e-7 over 60 steps via `npm run ref`), and the flux and Algorithm-4
10+ * exports track each other to ~3e-10 in f64.
11+ */
12+import { generateMatlabScript, MATLAB_SCRIPT_NAME } from '../src/export/matlabScript.ts';
13+import { presets, mModelByKey } from '../src/mgpu/registry.ts';
14+import { mGeometries } from '../src/geom/registry.ts';
15+import { formatCommand, resolvePreset, DEFAULT_WARMUP } from '../src/bench/runSpec.ts';
16+
17+type Check = (name: string, ok: boolean, detail: string) => void;
18+type Log = (s: string) => void;
19+
20+export function matlabExportChecks(check: Check, log: Log): void {
21+ log('--- MATLAB export ---');
22+ for (const preset of presets) {
23+ const { model, params } = resolvePreset(preset.key);
24+ for (const geometry of mGeometries) {
25+ const geometryParams = Object.fromEntries(
26+ geometry.params.map((p) => [p.key, p.value]),
27+ );
28+ const spec = {
29+ preset: preset.key,
30+ lmax: 63,
31+ seed: 1,
32+ steps: 2000,
33+ warmup: DEFAULT_WARMUP,
34+ params,
35+ geometry: geometry.key,
36+ geometryParams,
37+ niter: 8,
38+ };
39+ const name = `matlab-export ${preset.key} on ${geometry.key}`;
40+ let text: string;
41+ try {
42+ text = generateMatlabScript({
43+ model,
44+ modelSource: model.source,
45+ params,
46+ geometry,
47+ geometrySource: geometry.source,
48+ geometryParams,
49+ lmax: 63,
50+ niter: 8,
51+ lam3: 0.5,
52+ seed: 1,
53+ preset: preset.key,
54+ command: formatCommand(spec),
55+ });
56+ } catch (e) {
57+ check(name, false, e instanceof Error ? e.message : String(e));
58+ continue;
59+ }
60+
61+ // One file, one namespace: every local function name must be unique,
62+ // or MATLAB silently shadows one definition with another.
63+ const fnNames = [...text.matchAll(/^[ \t]*function\s+(?:\[[^\]]*\]|\w+)\s*=\s*(\w+)\s*\(/gm)]
64+ .map((m) => m[1]);
65+ const dupes = fnNames.filter((n, i) => fnNames.indexOf(n) !== i);
66+
67+ // The driver must define what it calls: the state it steps, the model
68+ // call mapped through the mp struct, and the transform setup.
69+ const wants = [
70+ `function ${MATLAB_SCRIPT_NAME}()`,
71+ 'sht_tables(sht_setup(lmax, mmax, nlat, nphi));',
72+ `= init(`,
73+ `= step(${model.state.join(', ')}, `,
74+ 'surface_tables(gxr, gyr, gzr)',
75+ `'/final/${model.state[0]}'`,
76+ ];
77+ const missing = wants.filter((w) => !text.includes(w));
78+
79+ // randnfunsphere rides along exactly when the geometry draws on it.
80+ const wantsSphereTool = /\brandnfunsphere\b/.test(geometry.source);
81+ const carriesSphereTool = /function f = randnfunsphere\(/.test(text);
82+
83+ const problems = [
84+ ...(dupes.length ? [`duplicate local functions: ${[...new Set(dupes)].join(', ')}`] : []),
85+ ...(missing.length ? [`missing: ${missing.join(' | ')}`] : []),
86+ ...(wantsSphereTool !== carriesSphereTool
87+ ? [`randnfunsphere ${wantsSphereTool ? 'missing' : 'included needlessly'}`]
88+ : []),
89+ ];
90+ check(name, problems.length === 0, problems.join('; ') || `${fnNames.length} local functions`);
91+ }
92+ }
93+
94+ // An edited working copy that dropped a required function is refused with a
95+ // message naming the file, not exported broken.
96+ const { model, params } = resolvePreset(presets[0].key);
97+ const geometry = mGeometries[0];
98+ try {
99+ generateMatlabScript({
100+ model,
101+ modelSource: '% nothing here',
102+ params,
103+ geometry,
104+ geometrySource: geometry.source,
105+ geometryParams: {},
106+ lmax: 63,
107+ niter: 8,
108+ lam3: 0.5,
109+ seed: 1,
110+ preset: presets[0].key,
111+ command: '',
112+ });
113+ check('matlab-export refuses a source without init', false, 'no error thrown');
114+ } catch (e) {
115+ const msg = e instanceof Error ? e.message : String(e);
116+ check(
117+ 'matlab-export refuses a source without init',
118+ msg.includes("'init'") && msg.includes(model.key),
119+ msg,
120+ );
121+ }
122+ // Guard the assumption the model registry makes for stateFor: state names
123+ // are used as `<name>0` initial-capture variables, which must not collide
124+ // with the species names.
125+ for (const p of presets) {
126+ const m = mModelByKey(p.modelKey)!;
127+ const all = new Set([...m.state, ...m.species]);
128+ check(
129+ `matlab-export names disjoint for ${m.key}`,
130+ all.size === m.state.length + m.species.length,
131+ [...all].join(', '),
132+ );
133+ }
134+}
test/referenceChecks.tsadded+132−0View file
@@ -0,0 +1,132 @@
1+/**
2+ * The reference-file reader, against a file this test writes itself.
3+ *
4+ * No GPU: this is about the format — that what h5wasm writes in the
5+ * documented layout (docs/ellipsoid-reference-spec.md) comes back through
6+ * `extractReferenceCase` with nothing renamed, rescaled or truncated, and
7+ * that a file the replay could not act on is refused with a message rather
8+ * than half-read. The h5wasm module is injected: the node harness passes
9+ * `h5wasm/node` (real files), the browser harness `h5wasm` (in-memory wasm
10+ * filesystem) — so the browser run also proves the wasm build actually ships.
11+ */
12+import { extractReferenceCase, type H5Node } from '../src/compare/referenceCase.ts';
13+import { nlmCalc } from '../src/sht/layout.ts';
14+
15+type Check = (name: string, ok: boolean, detail: string) => void;
16+type Log = (line: string) => void;
17+
18+/** The slice of h5wasm's writing API these checks touch — the node and
19+ * browser builds both satisfy it structurally. */
20+interface H5Out {
21+ create_group(name: string): H5Out;
22+ create_attribute(name: string, data: unknown): void;
23+ create_dataset(args: { name: string; data: unknown; dtype?: string }): unknown;
24+}
25+export interface H5Rt {
26+ ready: Promise<unknown>;
27+ File: new (path: string, mode?: string) => H5Out & H5Node & { close(): unknown };
28+}
29+
30+const LMAX = 3;
31+const STEPS = 8;
32+
33+export async function referenceChecks(
34+ h5: H5Rt,
35+ /** Where a named scratch file may live: a temp dir on node, '/' in the
36+ * browser's in-memory filesystem. */
37+ pathFor: (name: string) => string,
38+ check: Check,
39+ log: Log,
40+): Promise<void> {
41+ log('\nreference files (HDF5 layout):');
42+ const mod = (await h5.ready) as { FS?: { unlink(path: string): void } };
43+ const nlm = nlmCalc(LMAX, LMAX);
44+ const series = (offset: number): Float32Array =>
45+ Float32Array.from({ length: 2 * nlm }, (_, i) => offset + i / 16);
46+ const arrays = {
47+ Gx: series(100), Gy: series(200), Gz: series(300),
48+ initialU: series(1), finalU: series(2),
49+ };
50+
51+ // ---- write the documented layout, read it back ---------------------------
52+ const goodPath = pathFor('ref-roundtrip.h5');
53+ {
54+ const f = new h5.File(goodPath, 'w');
55+ f.create_attribute('model', 'allencahn');
56+ f.create_attribute('species', ['U']);
57+ const spec = f.create_group('spec');
58+ spec.create_attribute('geometry', 'ellipsoid');
59+ spec.create_attribute('lmax', LMAX);
60+ spec.create_attribute('steps', STEPS);
61+ spec.create_attribute('niter', 2);
62+ spec.create_attribute('seed', 1);
63+ spec.create_attribute('warmup', 0);
64+ const params = spec.create_group('params');
65+ params.create_attribute('dt', 0.0625);
66+ params.create_attribute('eps2', 0.5);
67+ const geomParams = spec.create_group('geometry_params');
68+ geomParams.create_attribute('ax', 2.5);
69+ geomParams.create_attribute('ay', 1.25);
70+ geomParams.create_attribute('az', 0.75);
71+ const geom = f.create_group('geometry');
72+ geom.create_dataset({ name: 'Gx', data: arrays.Gx, dtype: '<f4' });
73+ geom.create_dataset({ name: 'Gy', data: arrays.Gy, dtype: '<f4' });
74+ geom.create_dataset({ name: 'Gz', data: arrays.Gz, dtype: '<f4' });
75+ f.create_group('initial').create_dataset({ name: 'U', data: arrays.initialU, dtype: '<f4' });
76+ f.create_group('final').create_dataset({ name: 'U', data: arrays.finalU, dtype: '<f4' });
77+ f.close();
78+ }
79+ {
80+ const f = new h5.File(goodPath, 'r');
81+ const rc = extractReferenceCase(f, 'ref-roundtrip.h5');
82+ f.close();
83+ mod.FS?.unlink(goodPath);
84+
85+ check(
86+ 'reference: the run identity survives the round trip',
87+ rc.model.key === 'allencahn' && rc.geometry.key === 'ellipsoid' &&
88+ rc.lmax === LMAX && rc.steps === STEPS && rc.niter === 2,
89+ `${rc.model.key} on ${rc.geometry.key}, lmax ${rc.lmax}, ` +
90+ `${rc.steps} steps, niter ${rc.niter}`,
91+ );
92+ check(
93+ 'reference: the file’s parameters override the defaults',
94+ rc.params.dt === 0.0625 && rc.params.eps2 === 0.5 &&
95+ rc.geometryParams.ax === 2.5 && rc.geometryParams.ay === 1.25 &&
96+ rc.geometryParams.az === 0.75,
97+ `dt ${rc.params.dt}, eps2 ${rc.params.eps2}, ` +
98+ `ax/ay/az ${rc.geometryParams.ax}/${rc.geometryParams.ay}/${rc.geometryParams.az}`,
99+ );
100+ const same = (a: Float32Array, b: Float32Array): boolean =>
101+ a.length === b.length && a.every((v, i) => v === b[i]);
102+ check(
103+ 'reference: every coefficient array comes back bit-exact',
104+ same(rc.geometryCoeffs.X, arrays.Gx) && same(rc.geometryCoeffs.Y, arrays.Gy) &&
105+ same(rc.geometryCoeffs.Z, arrays.Gz) && same(rc.initial.U, arrays.initialU) &&
106+ same(rc.final.U, arrays.finalU),
107+ `5 arrays x ${2 * nlm} float32 values`,
108+ );
109+ }
110+
111+ // ---- a file the replay cannot act on is refused, not half-read -----------
112+ {
113+ const badPath = pathFor('ref-unknown-model.h5');
114+ const f = new h5.File(badPath, 'w');
115+ f.create_attribute('model', 'nosuchmodel');
116+ f.close();
117+ const r = new h5.File(badPath, 'r');
118+ let message = '';
119+ try {
120+ extractReferenceCase(r, 'ref-unknown-model.h5');
121+ } catch (e) {
122+ message = e instanceof Error ? e.message : String(e);
123+ }
124+ r.close();
125+ mod.FS?.unlink(badPath);
126+ check(
127+ 'reference: an unknown model is refused with its name',
128+ message.includes('nosuchmodel'),
129+ message || 'no error thrown',
130+ );
131+ }
132+}
test/test-page.tsmodified+6−0View file
@@ -23,12 +23,15 @@ import {
2323 defaultGeometryParams,
2424 DEFAULT_GEOMETRY_KEY,
2525 } from '../src/geom/registry.ts';
26+import * as h5wasm from 'h5wasm';
2627 import { transformChecks } from './transformChecks.ts';
2728 import { analyticChecks } from './analyticChecks.ts';
2829 import { modelChecks } from './modelChecks.ts';
2930 import { geometryChecks } from './geometryChecks.ts';
3031 import { fluxChecks } from './fluxChecks.ts';
3132 import { compareChecks } from './compareChecks.ts';
33+import { referenceChecks, type H5Rt } from './referenceChecks.ts';
34+import { matlabExportChecks } from './matlabExportChecks.ts';
3235
3336 declare global {
3437 interface Window {
@@ -222,6 +225,9 @@ async function main(): Promise<void> {
222225 await geometryChecks(device, check, log, { sweep: q.has('sweep') });
223226 await fluxChecks(device, check, log, { ab: q.has('sweep') });
224227 await compareChecks(device, check, log);
228+ // '/' is the wasm module's in-memory filesystem — nothing touches disk.
229+ await referenceChecks(h5wasm as unknown as H5Rt, (name) => `/${name}`, check, log);
230+ matlabExportChecks(check, log);
225231
226232 window.__RESULTS__ = { ok: failures === 0, lines };
227233 log(failures === 0 ? 'ALL PASS' : `${failures} FAILURE(S)`);