forked from github.com/concept-collection/turing-surface
Comparing changes
main is 7 commits ahead of compare-interface.
Create pull request
Merge pull request #11 from concept-collection/compare-against-sphere
Jeremy Magland committed
9321562Merge main into compare-against-sphere
Jeremy Magland committed
9916f17Accept a surface handed over as coefficients (reharm's export button)
Dan Fortunato committed
77f627eexport standalone matlab
Jeremy Magland committed
9b71d5bWIP: comparison with unit sphere
Owen Melia committed
cf4af12Update comparison UI
Owen Melia committed
02e7a36Merge pull request #10 from concept-collection/compare-interface
Jeremy Magland committed
71b649216 changed files+2119−93
.gitignoremodified+1−0View file
@@ -1,3 +1,4 @@
1+tmp/
12 node_modules/
23 dist/
34 *.log
README.mdmodified+16−0View file
@@ -68,6 +68,16 @@ solver and does not disturb the run: the geometry is data whose shape in the
6868 bindings depends only on the grid, so a swap is sixteen buffer writes and the
6969 pattern carries straight on.
7070
71+A sixth kind of geometry arrives from outside: reharm (the sibling
72+spherical-harmonic fitting app) can hand a fitted surface — a scanned mesh, a
73+brain, a point cloud — to this page with its **Export to turing-surface**
74+button. The surface arrives as coefficients by `postMessage` (there is no .m;
75+`Geometry` is built from them directly, re-indexed onto whatever band is
76+selected), lands as an **imported** entry in the geometry dropdown at its own
77+lmax, and behaves like any other shape from there — except that the editor,
78+the benchmark command line and the MATLAB export, all of which need a .m, say
79+so instead.
80+
7181 A **morph** slider blends the drawn surface back to the unit sphere. The
7282 parametrization is the sphere's either way, so sweeping it shows which point
7383 went where.
@@ -497,6 +507,12 @@ package alone. Its binaries need glibc 2.29+. Other flags: `--steps`,
497507 `--warmup`, `--batch`, `--json`, `--help`; `DAWN_FLAGS='backend=vulkan'`
498508 (`;`-separated) passes Dawn options through.
499509
510+### The same run in MATLAB
511+
512+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.
513+
514+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.
515+
500516 ## Tests
501517
502518 There is no second implementation of the solver to diff against, so the `.m`
index.htmlmodified+71−3View file
@@ -117,6 +117,14 @@
117117 font: 12px/1.6 ui-monospace, SFMono-Regular, Menlo, Consolas, monospace;
118118 color: var(--ink); user-select: all;
119119 }
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+ }
120128 #blurb { margin-top: 4px; font-size: 13px; color: var(--ink-2); }
121129 #err { color: #b35900; white-space: pre-wrap; font-size: 13px; }
122130 .editor {
@@ -177,6 +185,12 @@
177185 /* The bar's own layout: three chip rows stacked, then the reference
178186 picker and the button beside them. */
179187 .cmp-axes { display: flex; flex-direction: column; gap: 4px; }
188+ /* Without this, JS setting the `hidden` attribute (vs-sphere mode,
189+ which has no chip grid) does nothing: [hidden] and .cmp-axes tie on
190+ specificity, and the author rule above wins the tie over the
191+ browser's default — the exact opposite of what `hidden` is supposed
192+ to mean. */
193+ .cmp-axes[hidden] { display: none; }
180194 .cmp-axis { display: flex; align-items: center; gap: 8px; }
181195 .cmp-axis > span:first-child {
182196 color: var(--ink-2); font-size: 13px; width: 5.5em; text-align: right;
@@ -225,6 +239,47 @@
225239 border: 1px solid var(--line); border-radius: 2px;
226240 }
227241 .cmp-rangelab { font-variant-numeric: tabular-nums; white-space: nowrap; }
242+ /* The diff row sits right under its variant's value row: same label
243+ width so the columns still line up, but a lighter tone since it names
244+ no new variant, only the difference from one already named above. */
245+ .cmp-diffrow .cmp-rowlabel { color: var(--ink-2); font-style: italic; }
246+ .cmp-diffbox { position: relative; }
247+ .cmp-diffcap {
248+ position: absolute; bottom: 4px; right: 6px;
249+ font-size: 11px; font-variant-numeric: tabular-nums;
250+ color: #fff; background: rgba(0, 0, 0, 0.45);
251+ padding: 1px 5px; border-radius: 4px; pointer-events: none;
252+ }
253+ .cmp-chart {
254+ margin-top: 8px; padding: 8px 10px;
255+ border: 1px solid var(--line); border-radius: 8px;
256+ }
257+ .cmp-chart-legend {
258+ display: flex; flex-wrap: wrap; gap: 6px 18px;
259+ font-size: 12px; color: var(--ink-2); margin-bottom: 6px;
260+ }
261+ .cmp-chart-legend-group { display: flex; flex-wrap: wrap; gap: 12px; }
262+ /* The species group (dash styles) reads as a separate cluster from the
263+ variant group (colors) — a rule between them, not just a gap, since
264+ both groups otherwise look like more of the same kind of chip. */
265+ .cmp-chart-legend-species {
266+ padding-left: 14px; border-left: 1px solid var(--line);
267+ }
268+ .cmp-chart-legend-item { display: inline-flex; align-items: center; gap: 5px; }
269+ .cmp-chart-legend-item i {
270+ display: inline-block; width: 10px; height: 10px; border-radius: 2px;
271+ }
272+ /* Species swatches reuse the color chip's box model but draw a dashed
273+ line instead of a filled square, so the two legend groups still read
274+ as one visual language (a small icon + a name). */
275+ .cmp-chart-legend-species i.cmp-chart-dash {
276+ width: 16px; height: 0; border-radius: 0;
277+ border-top-width: 2px; border-top-color: var(--ink-2);
278+ }
279+ .cmp-chart-dash-0 { border-top-style: solid; }
280+ .cmp-chart-dash-1 { border-top-style: dashed; }
281+ .cmp-chart-dash-2 { border-top-style: dotted; }
282+ .cmp-chart-canvas { width: 100%; height: 240px; display: block; }
228283 @media (max-width: 860px) {
229284 .cmp-row { flex-direction: column; }
230285 .cmp-rowlabel { width: auto; flex-direction: row; gap: 10px; }
@@ -252,7 +307,9 @@
252307 <button type="button" class="chip" id="mode-simulate" aria-pressed="true">Simulate</button>
253308 <button type="button" class="chip" id="mode-effort"
254309 title="Run several solver settings side by side on one clock">Compare computational effort</button>
255- <button type="button" class="chip" id="mode-vs-sphere" disabled title="Coming soon">Compare against sphere (Coming soon)</button>
310+ <button type="button" class="chip" id="mode-vs-sphere"
311+ title="Run this model once on the true geometry and once on the unit sphere, from the same starting state">
312+ Compare against sphere</button>
256313 <button type="button" class="chip" id="mode-vs-upload"
257314 title="Check this solver against a saved reference run">
258315 Compare against uploaded data</button>
@@ -272,7 +329,7 @@
272329 <div class="controls" id="geomparams"></div>
273330 </div>
274331 <div class="controls" id="comparebar" hidden>
275- <div class="cmp-axes">
332+ <div class="cmp-axes" id="cmp-axes">
276333 <div class="cmp-axis">
277334 <span title="Iterations of the implicit diffusion solve">solve iters</span>
278335 <span id="cmp-niter" class="chips"></span>
@@ -286,7 +343,7 @@
286343 <span id="cmp-dt" class="chips"></span>
287344 </div>
288345 </div>
289- <label title="The run everything else is measured against">reference
346+ <label title="The run everything else is measured against" id="cmp-reflabel">reference
290347 <select id="cmp-ref"></select>
291348 </label>
292349 <span id="cmp-fileinfo" class="stats" hidden></span>
@@ -380,6 +437,7 @@
380437 </div>
381438 <p id="geomnote" class="stats"></p>
382439 <div id="panels"></div>
440+ <div id="cmp-chart" class="cmp-chart" hidden></div>
383441 <p class="stats" id="stats"></p>
384442 <p class="stats" id="benchresult"></p>
385443 <div class="editor">
@@ -414,6 +472,16 @@
414472 </div>
415473 <code id="cmd"></code>
416474 </div>
475+ <details class="cli" id="matlab">
476+ <summary class="cli-head">
477+ <span>The same run as a standalone MATLAB script</span>
478+ <span>
479+ <button id="copymatlab" type="button">Copy</button>
480+ <button id="downloadmatlab" type="button">Download .m</button>
481+ </span>
482+ </summary>
483+ <pre id="matlabscript"></pre>
484+ </details>
417485 <p id="blurb"></p>
418486 <p id="err"></p>
419487 </main>
scripts/test-node.tsmodified+2−0View file
@@ -20,6 +20,7 @@ import { geometryChecks } from '../test/geometryChecks.ts';
2020 import { fluxChecks } from '../test/fluxChecks.ts';
2121 import { compareChecks } from '../test/compareChecks.ts';
2222 import { referenceChecks, type H5Rt } from '../test/referenceChecks.ts';
23+import { matlabExportChecks } from '../test/matlabExportChecks.ts';
2324
2425 let failures = 0;
2526 const check = (name: string, ok: boolean, detail: string): void => {
@@ -59,6 +60,7 @@ await modelChecks(device, check, log);
5960 await geometryChecks(device, check, log);
6061 await fluxChecks(device, check, log);
6162 await compareChecks(device, check, log);
63+matlabExportChecks(check, log);
6264 await referenceChecks(
6365 h5wasm as unknown as H5Rt,
6466 (name) => join(tmpdir(), `turing-surface-${process.pid}-${name}`),
src/compare/compareRun.tsmodified+297−53View file
@@ -44,6 +44,7 @@ import { fmtValue, floorRange } from '../render/colorbar.ts';
4444 import { prolongCoeffs, sharedModes, sharedNoise } from './sharedStart.ts';
4545 import { variantLabel, VARIANT_COLORS, type Variant } from './variants.ts';
4646 import type { ReferenceCase } from './referenceCase.ts';
47+import { ErrorChart, type ErrorChartRow } from '../render/errorChart.ts';
4748
4849 /**
4950 * Latitudes of the shared display grid. 256 is the same target the single-run
@@ -76,6 +77,21 @@ export interface CompareOptions {
7677 geometry: MGeometry;
7778 geometryParams: Params;
7879 geometrySource: string;
80+ /** Per-row geometry override, parallel-indexed to `variants` — the
81+ * "several rows on several geometries" case (vs-sphere mode), as opposed
82+ * to every other mode's "several rows, one geometry." Rows without an
83+ * entry (or when this is omitted entirely) fall back to the single
84+ * `geometry`/`geometryParams`/`geometrySource` above. */
85+ geometries?: { geometry: MGeometry; geometryParams: Params; geometrySource: string }[];
86+ /** Render every row — value panels and diff panels alike — on the
87+ * *reference* row's own surface (`coords`/`posBuf`) instead of each row's
88+ * own. For comparing two different geometries where only the *field*
89+ * should read as different, not the displayed shape. */
90+ renderOnReferenceGeometry?: boolean;
91+ /** Row names, overriding `variantLabel(variant, showDt)` — for a mode
92+ * where every row shares the same niter/lmax/dt, so that label alone
93+ * wouldn't tell the rows apart. */
94+ rowLabels?: string[];
7995 variants: Variant[];
8096 /** Index into `variants` of the run everything else is measured against.
8197 * Ignored when `refFile` is given — the file is the reference then. */
@@ -98,6 +114,8 @@ export interface CompareOptions {
98114 colormapName: () => string;
99115 /** Where the variant grid goes (the app's #panels). */
100116 container: HTMLElement;
117+ /** Where the error-vs-time chart goes (the app's #cmp-chart). */
118+ chartContainer: HTMLElement;
101119 /** Progress and, afterwards, the standing description of the study. */
102120 onStatus: (html: string) => void;
103121 }
@@ -114,6 +132,21 @@ interface Row {
114132 colorBufs: Float32Array[];
115133 /** Fields read this frame, one per species, on the shared grid. */
116134 fields: Float32Array[];
135+ /** Pointwise difference from the reference, one per species, on the shared
136+ * grid — filled by #measureDifference as it accumulates the norm below.
137+ * Unused (and never touched) for the reference row itself. */
138+ diffFields: Float32Array[];
139+ /** The diff row's own panels, one per species — empty for the reference
140+ * row, whose diff against itself is trivially zero. */
141+ diffScenes: SphereScene[];
142+ diffValueBufs: Float32Array[];
143+ diffColorBufs: Float32Array[];
144+ /** This row's own smoothed symmetric color range per species — each row
145+ * is scaled to its own diff extent, not a range shared across the
146+ * column, so one row's diff panels never affect another's coloring. */
147+ diffRanges: { lo: number; hi: number }[];
148+ /** Each diff panel's "± magnitude" caption, parallel to diffScenes. */
149+ diffCaps: HTMLElement[];
117150 /** Relative difference from the reference, one per species. */
118151 err: number[];
119152 /** False once any species has left the floating-point numbers — the shape a
@@ -161,6 +194,7 @@ export class CompareRun {
161194 /** Smoothed color range per species, shared by every variant so the panels
162195 * in a column are directly comparable by eye and not just by number. */
163196 #ranges: { lo: number; hi: number }[] = [];
197+ #errorChart: ErrorChart | null = null;
164198 #resizeObs: ResizeObserver | null = null;
165199
166200 #running = false;
@@ -234,10 +268,12 @@ export class CompareRun {
234268 // canvases from the DOM would leave both running.
235269 let built: Row[] = [];
236270 let builtFile: FileRow | null = null;
271+ let errorChart: ErrorChart | null = null;
237272
238273 try {
239274 for (let i = 0; i < variants.length; i++) {
240275 const v = variants[i];
276+ const g = opts.geometries?.[i];
241277 opts.onStatus(
242278 `compiling ${i + 1}/${variants.length} — ${variantLabel(v, showDt)} ` +
243279 `(a solve iteration is ~15 kernels per species, and there is no ` +
@@ -252,9 +288,9 @@ export class CompareRun {
252288 params: { ...opts.params, dt: baseDt / v.dtDiv },
253289 lmax: v.lmax,
254290 source: opts.source,
255- geometry: opts.geometry,
256- geometryParams: opts.geometryParams,
257- geometrySource: opts.geometrySource,
291+ geometry: g?.geometry ?? opts.geometry,
292+ geometryParams: g?.geometryParams ?? opts.geometryParams,
293+ geometrySource: g?.geometrySource ?? opts.geometrySource,
258294 niter: v.niter,
259295 lam3: opts.lam3,
260296 }),
@@ -298,10 +334,30 @@ export class CompareRun {
298334 // One at a time: a seed submits its whole mode sum in pieces, and there
299335 // is nothing to gain from interleaving several variants' worth of it.
300336 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;
337+ if (opts.geometries) {
338+ // Several geometries, not several lmax bands: the seeding above
339+ // drew a spatially shared field, but evaluated it on each row's
340+ // own (geometry-dependent) points — a different field per row in
341+ // coefficient space. The reference's exact resulting coefficients
342+ // replace that for every other row, so every row starts from the
343+ // identical spectral state and only the operator applied to it
344+ // differs from then on.
345+ const refSession = sessions[opts.reference];
346+ const refState = await refSession.readState();
347+ for (let i = 0; i < sessions.length; i++) {
348+ if (i === opts.reference) continue;
349+ sessions[i].loadState(
350+ prolongState(refState, model.state, refSession.cfg.lmax, sessions[i].cfg.lmax),
351+ );
352+ }
353+ initial = refState;
354+ initialLmax = refSession.cfg.lmax;
355+ } else {
356+ let coarsest = sessions[0];
357+ for (const s of sessions) if (s.cfg.lmax < coarsest.cfg.lmax) coarsest = s;
358+ initial = await coarsest.readState();
359+ initialLmax = coarsest.cfg.lmax;
360+ }
305361 }
306362
307363 // ---- the mesh, shared; the surface, per variant ---------------------
@@ -336,6 +392,17 @@ export class CompareRun {
336392 built = rows;
337393 builtFile = fileRow;
338394
395+ // ---- the error-vs-time chart, one line per row with a diff panel ----
396+ const chartRows = rows.filter((r) => r.diffScenes.length > 0);
397+ errorChart = new ErrorChart(
398+ opts.chartContainer,
399+ model.species,
400+ chartRows.map((r): ErrorChartRow => ({ label: variantLabel(r.variant, showDt), color: r.color })),
401+ );
402+ // Nothing to chart with a single variant and no reference file — the
403+ // one row present is the reference itself.
404+ opts.chartContainer.hidden = chartRows.length === 0;
405+
339406 const solverGrid = sessions.map((s) => `${s.cfg.nlat}×${s.cfg.nphi}`);
340407 const note =
341408 `${variants.length} variant${variants.length === 1 ? '' : 's'} · ` +
@@ -349,14 +416,20 @@ export class CompareRun {
349416 const run = new CompareRun({
350417 opts, rows, fileRow, topo, weights, rangeBars, frameSteps, note, initial, initialLmax,
351418 });
419+ run.#errorChart = errorChart;
352420 await run.draw();
353421 run.#observeResize();
354422 run.#status();
355423 return run;
356424 } catch (e) {
357- for (const r of built) for (const s of r.scenes) s.dispose();
425+ for (const r of built) {
426+ for (const s of r.scenes) s.dispose();
427+ for (const s of r.diffScenes) s.dispose();
428+ }
358429 for (const s of builtFile?.scenes ?? []) s.dispose();
359430 for (const s of sessions) s.destroy();
431+ errorChart?.dispose();
432+ opts.chartContainer.hidden = true;
360433 opts.container.replaceChildren();
361434 opts.container.classList.remove('compare');
362435 throw e;
@@ -400,10 +473,26 @@ export class CompareRun {
400473 // This draw becomes what restart() rewinds to from now on — see the
401474 // identical selection in create(). Recaptured here rather than left
402475 // 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;
476+ if (this.#opts.geometries) {
477+ // Mirrors create()'s IC block: copy the reference's exact resulting
478+ // coefficients into every other row rather than trusting their own
479+ // (geometry-dependent) seeding to have landed on the same state.
480+ const refSession = sessions[this.#opts.reference];
481+ const refState = await refSession.readState();
482+ for (let i = 0; i < sessions.length; i++) {
483+ if (i === this.#opts.reference) continue;
484+ sessions[i].loadState(
485+ prolongState(refState, this.#opts.model.state, refSession.cfg.lmax, sessions[i].cfg.lmax),
486+ );
487+ }
488+ this.#initial = refState;
489+ this.#initialLmax = refSession.cfg.lmax;
490+ } else {
491+ let coarsest = sessions[0];
492+ for (const s of sessions) if (s.cfg.lmax < coarsest.cfg.lmax) coarsest = s;
493+ this.#initial = await coarsest.readState();
494+ this.#initialLmax = coarsest.cfg.lmax;
495+ }
407496 }
408497 this.#t = 0;
409498 this.#stepsDone = 0;
@@ -412,6 +501,13 @@ export class CompareRun {
412501 r.lo = NaN;
413502 r.hi = NaN;
414503 }
504+ for (const r of this.#rows) {
505+ for (const rr of r.diffRanges) {
506+ rr.lo = NaN;
507+ rr.hi = NaN;
508+ }
509+ }
510+ this.#errorChart?.reset();
415511 await this.draw();
416512 this.#status();
417513 if (!this.#disposed && wasRunning) this.setRunning(true);
@@ -437,6 +533,13 @@ export class CompareRun {
437533 r.lo = NaN;
438534 r.hi = NaN;
439535 }
536+ for (const r of this.#rows) {
537+ for (const rr of r.diffRanges) {
538+ rr.lo = NaN;
539+ rr.hi = NaN;
540+ }
541+ }
542+ this.#errorChart?.reset();
440543 await this.draw();
441544 this.#status();
442545 if (!this.#disposed && wasRunning) this.setRunning(true);
@@ -473,6 +576,11 @@ export class CompareRun {
473576 for (const r of this.#rows) {
474577 fillPositions(r.posBuf, r.coords, this.#topo, morph);
475578 for (const s of r.scenes) s.updatePositions(r.posBuf);
579+ // The diff row sits on the exact same mesh as the value row above it
580+ // (same coords, same posBuf) — it just never got told to re-render
581+ // when this method was first written, so it stayed fixed at whatever
582+ // shape the study was compiled with.
583+ for (const s of r.diffScenes) s.updatePositions(r.posBuf);
476584 }
477585 const f = this.#fileRow;
478586 if (f) {
@@ -492,17 +600,24 @@ export class CompareRun {
492600 this.#resizeObs = null;
493601 for (const r of this.#rows) {
494602 for (const s of r.scenes) s.dispose();
603+ for (const s of r.diffScenes) s.dispose();
495604 r.session.destroy();
496605 }
497606 for (const s of this.#fileRow?.scenes ?? []) s.dispose();
498607 this.#rows = [];
499608 this.#fileRow = null;
609+ this.#errorChart?.dispose();
610+ this.#errorChart = null;
611+ this.#opts.chartContainer.hidden = true;
500612 this.#opts.container.replaceChildren();
501613 this.#opts.container.classList.remove('compare');
502614 }
503615
504616 #allScenes(): SphereScene[] {
505- return [...this.#rows.flatMap((r) => r.scenes), ...(this.#fileRow?.scenes ?? [])];
617+ return [
618+ ...this.#rows.flatMap((r) => [...r.scenes, ...r.diffScenes]),
619+ ...(this.#fileRow?.scenes ?? []),
620+ ];
506621 }
507622
508623 // ----------------------------------------------------------------- drawing
@@ -611,7 +726,54 @@ export class CompareRun {
611726 }
612727
613728 this.#measureDifference();
729+ this.#colorDiffPanels();
614730 this.#updateRowStats();
731+ this.#errorChart?.push(
732+ this.#t,
733+ this.#rows.filter((r) => r.diffScenes.length > 0).map((r) => r.err),
734+ );
735+ }
736+
737+ /**
738+ * Color every diff panel from #measureDifference's pointwise diffFields, on
739+ * a *symmetric* range that is each row's own — not shared across the
740+ * column. A diverging row's diff explodes, but with a per-row range that
741+ * only saturates its own panels; it can no longer affect how any other
742+ * row's diff panels are scaled, which is a simpler fix for exactly the
743+ * flooding problem the value panels' shared #ranges/leastPeak logic exists
744+ * to manage there. The cost: diff-panel color alone no longer tells you
745+ * which row has more error than another — that comparison now lives in the
746+ * error chart, which has actual numbers. Always drawn with a diverging
747+ * colormap regardless of the user's chosen (sequential) one — a signed
748+ * quantity centered at zero needs a diverging map to read correctly, which
749+ * coolwarm is and viridis etc. are not.
750+ */
751+ #colorDiffPanels(): void {
752+ const species = this.#opts.model.species;
753+ for (const r of this.#rows) {
754+ if (r.diffScenes.length === 0) continue;
755+ for (let k = 0; k < species.length; k++) {
756+ const m = r.healthy ? maxAbsFinite(r.diffFields[k]) : null;
757+ const range = r.diffRanges[k];
758+ if (m !== null) {
759+ if (!Number.isFinite(range.lo)) {
760+ range.lo = -m;
761+ range.hi = m;
762+ } else {
763+ const a = 0.15;
764+ range.lo += a * (-m - range.lo);
765+ range.hi += a * (m - range.hi);
766+ }
767+ }
768+ if (!Number.isFinite(range.lo) || !Number.isFinite(range.hi)) continue;
769+ const shown = floorRange(range.lo, range.hi);
770+ fillFieldValues(r.diffValueBufs[k], r.diffFields[k], this.#topo);
771+ fillColors(r.diffColorBufs[k], r.diffValueBufs[k], shown.lo, shown.hi, colormaps.coolwarm);
772+ r.diffScenes[k]?.updateColors(r.diffColorBufs[k]);
773+ const cap = r.diffCaps[k];
774+ if (cap) cap.textContent = `± ${fmtValue(shown.hi)}`;
775+ }
776+ }
615777 }
616778
617779 /**
@@ -640,11 +802,13 @@ export class CompareRun {
640802 r.err[k] = NaN;
641803 continue;
642804 }
805+ const diff = r.diffFields[k];
643806 let num = 0;
644807 let den = 0;
645808 for (let i = 0; i < a.length; i++) {
646809 const w = this.#weights[i];
647810 const d = a[i] - b[i];
811+ diff[i] = d;
648812 num += w * d * d;
649813 den += w * b[i] * b[i];
650814 }
@@ -661,21 +825,19 @@ export class CompareRun {
661825 * usually the one that has converged and the fast one the one that has not.
662826 */
663827 #updateRowStats(): void {
664- const species = this.#opts.model.species;
665828 const ref = this.#fileRow ? null : this.#rows[this.#opts.reference];
666829 for (const r of this.#rows) {
667- const per = species
668- .map((s, k) => `${s} ${Number.isFinite(r.err[k]) ? r.err[k].toExponential(2) : '—'}`)
669- .join('<br>');
670830 // Divergence is said, not implied. Scaled to a healthy row, a blown-up
671831 // variant is a flat saturated panel, which on its own is easy to misread
672- // as a converged uniform state.
832+ // as a converged uniform state. Δ itself now lives in the error chart,
833+ // not here.
673834 const body = !r.healthy
674835 ? '<b class="cmp-diverged">diverged</b>'
675836 : r === ref
676837 ? '<b>reference</b>'
677- : `Δ ${per}`;
678- r.statEl.innerHTML = `${r.session.steps.toLocaleString()} steps<br>${body}`;
838+ : '';
839+ r.statEl.innerHTML =
840+ `${r.session.steps.toLocaleString()} steps` + (body ? `<br>${body}` : '');
679841 }
680842 }
681843
@@ -683,9 +845,7 @@ export class CompareRun {
683845 const refFile = this.#opts.refFile;
684846 const clock = refFile
685847 ? `<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`)
848+ ` · NOTE: vertical axis measures difference from uploaded simulation's end state.`
689849 : `<b>t = ${this.#t.toFixed(2)}</b> (same for every variant)`;
690850 this.#opts.onStatus(
691851 `${clock} · ` +
@@ -701,11 +861,13 @@ export class CompareRun {
701861 const box = s.canvas.parentElement;
702862 if (box) s.resize(box.clientWidth, box.clientHeight);
703863 }
864+ this.#errorChart?.redraw();
704865 });
705866 for (const s of scenes) {
706867 const box = s.canvas.parentElement;
707868 if (box) this.#resizeObs.observe(box);
708869 }
870+ this.#resizeObs.observe(this.#opts.chartContainer);
709871 }
710872
711873 // -------------------------------------------------------------- the clock
@@ -801,6 +963,22 @@ function leastPeak(all: (Bounds | null)[]): Bounds | null {
801963 return best;
802964 }
803965
966+/** Max |value| over the finite entries of a field; null if none are finite —
967+ * the diff panels' analogue of finiteRange below, since a symmetric range
968+ * only needs the one number. */
969+function maxAbsFinite(f: Float32Array): number | null {
970+ let m = -Infinity;
971+ let any = false;
972+ for (let i = 0; i < f.length; i++) {
973+ const v = f[i];
974+ if (!Number.isFinite(v)) continue;
975+ any = true;
976+ const a = Math.abs(v);
977+ if (a > m) m = a;
978+ }
979+ return any ? m : null;
980+}
981+
804982 /** Min and max over the finite entries only; null when there are none. */
805983 function finiteRange(f: Float32Array | undefined): { lo: number; hi: number } | null {
806984 if (!f) return null;
@@ -826,6 +1004,38 @@ function finiteRange(f: Float32Array | undefined): { lo: number; hi: number } |
8261004 * a variant: it is the thing they are all measured against. */
8271005 const FILE_ROW_COLOR = '#57606a';
8281006
1007+/**
1008+ * One sphere panel: a boxed SphereScene plus the value/color buffers that
1009+ * feed it. Shared by the value row, the diff row, and the file row — all
1010+ * three build a panel the same way, only differing in the class on the box
1011+ * (for styling) and in what fills the buffers afterward.
1012+ */
1013+function makeSpherePanel(
1014+ colsEl: HTMLElement,
1015+ topo: SphereMeshTopology,
1016+ posBuf: Float32Array,
1017+ background: string | undefined,
1018+ extraClass = '',
1019+): { box: HTMLElement; scene: SphereScene; valueBuf: Float32Array; colorBuf: Float32Array } {
1020+ const box = document.createElement('div');
1021+ box.className = extraClass ? `sphere-box cmp-box ${extraClass}` : 'sphere-box cmp-box';
1022+ colsEl.append(box);
1023+ const scene = new SphereScene(
1024+ box,
1025+ topo.numVertices,
1026+ topo.indices,
1027+ Float32Array.from(posBuf),
1028+ background,
1029+ );
1030+ scene.fitCamera();
1031+ return {
1032+ box,
1033+ scene,
1034+ valueBuf: new Float32Array(topo.numVertices),
1035+ colorBuf: new Float32Array(topo.numVertices * 3),
1036+ };
1037+}
1038+
8291039 async function buildGrid(
8301040 opts: CompareOptions,
8311041 sessions: ModelSession[],
@@ -884,13 +1094,28 @@ async function buildGrid(
8841094 .getPropertyValue('--sphere-bg')
8851095 .trim();
8861096
1097+ // When every row should be drawn on the same shape (vs-sphere: the point
1098+ // is to isolate the field, not the surface), fetch that shape once from
1099+ // the reference row's session and hand the identical array to every row.
1100+ // fillPositions/fillFieldValues already treat "whose mesh this is" and
1101+ // "whose field this is" as fully independent buffers, so this is the only
1102+ // place that needs to know about it — every method downstream that reads
1103+ // r.coords/r.posBuf just sees one row's shape reused by every other.
1104+ const sharedCoords = opts.renderOnReferenceGeometry
1105+ ? await sessions[opts.reference].renderPositions()
1106+ : null;
1107+
8871108 const rows: Row[] = [];
8881109 for (let i = 0; i < sessions.length; i++) {
8891110 const session = sessions[i];
8901111 const variant = opts.variants[i];
8911112 const color = VARIANT_COLORS[i % VARIANT_COLORS.length];
1113+ // The reference itself never gets a diff row — its diff against itself
1114+ // is trivially zero, and a flat zero panel would just spend a WebGL
1115+ // context for nothing.
1116+ const isRef = !opts.refFile && i === opts.reference;
8921117
893- const coords = await session.renderPositions();
1118+ const coords = sharedCoords ?? (await session.renderPositions());
8941119 const posBuf = new Float32Array(topo.numVertices * 3);
8951120 fillPositions(posBuf, coords, topo, opts.morph);
8961121
@@ -901,7 +1126,7 @@ async function buildGrid(
9011126 labelEl.style.setProperty('--c', color);
9021127 const nameEl = document.createElement('div');
9031128 nameEl.className = 'cmp-rowname';
904- nameEl.textContent = variantLabel(variant, showDt);
1129+ nameEl.textContent = opts.rowLabels?.[i] ?? variantLabel(variant, showDt);
9051130 const statEl = document.createElement('div');
9061131 statEl.className = 'cmp-rowstat';
9071132 labelEl.append(nameEl, statEl);
@@ -914,25 +1139,52 @@ async function buildGrid(
9141139 const valueBufs: Float32Array[] = [];
9151140 const colorBufs: Float32Array[] = [];
9161141 for (let k = 0; k < model.species.length; k++) {
917- const box = document.createElement('div');
918- box.className = 'sphere-box cmp-box';
919- colsEl.append(box);
920- const scene = new SphereScene(
921- box,
922- topo.numVertices,
923- topo.indices,
924- Float32Array.from(posBuf),
925- sphereBg || undefined,
926- );
927- scene.fitCamera();
1142+ const { scene, valueBuf, colorBuf } = makeSpherePanel(colsEl, topo, posBuf, sphereBg || undefined);
9281143 scenes.push(scene);
929- valueBufs.push(new Float32Array(topo.numVertices));
930- colorBufs.push(new Float32Array(topo.numVertices * 3));
1144+ valueBufs.push(valueBuf);
1145+ colorBufs.push(colorBuf);
1146+ }
1147+
1148+ // ---- its diff row, right underneath ----------------------------------
1149+ const diffFields = model.species.map(() => new Float32Array(topo.nlat * topo.nphi));
1150+ const diffScenes: SphereScene[] = [];
1151+ const diffValueBufs: Float32Array[] = [];
1152+ const diffColorBufs: Float32Array[] = [];
1153+ const diffCaps: HTMLElement[] = [];
1154+ if (!isRef) {
1155+ const diffRowEl = document.createElement('div');
1156+ diffRowEl.className = 'cmp-row cmp-diffrow';
1157+ const diffLabelEl = document.createElement('div');
1158+ diffLabelEl.className = 'cmp-rowlabel';
1159+ diffLabelEl.style.setProperty('--c', color);
1160+ const diffNameEl = document.createElement('div');
1161+ diffNameEl.className = 'cmp-rowname';
1162+ diffNameEl.textContent = 'Δ vs reference';
1163+ diffLabelEl.append(diffNameEl);
1164+ const diffColsEl = document.createElement('div');
1165+ diffColsEl.className = 'cmp-cols';
1166+ diffRowEl.append(diffLabelEl, diffColsEl);
1167+ container.append(diffRowEl);
1168+
1169+ for (let k = 0; k < model.species.length; k++) {
1170+ const { box, scene, valueBuf, colorBuf } = makeSpherePanel(
1171+ diffColsEl, topo, posBuf, sphereBg || undefined, 'cmp-diffbox',
1172+ );
1173+ diffScenes.push(scene);
1174+ diffValueBufs.push(valueBuf);
1175+ diffColorBufs.push(colorBuf);
1176+ const cap = document.createElement('span');
1177+ cap.className = 'cmp-diffcap';
1178+ box.append(cap);
1179+ diffCaps.push(cap);
1180+ }
9311181 }
9321182
1183+ const diffRanges = model.species.map(() => ({ lo: NaN, hi: NaN }));
9331184 rows.push({
9341185 variant, session, color, coords, posBuf, scenes, valueBufs, colorBufs,
935- fields: [], err: model.species.map(() => 0), healthy: true, statEl,
1186+ fields: [], diffFields, diffScenes, diffValueBufs, diffColorBufs, diffCaps, diffRanges,
1187+ err: model.species.map(() => 0), healthy: true, statEl,
9361188 });
9371189 }
9381190
@@ -987,32 +1239,24 @@ async function buildGrid(
9871239 const fields: Float32Array[] = [];
9881240 const bounds: (Bounds | null)[] = [];
9891241 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();
1242+ const { scene, valueBuf, colorBuf } = makeSpherePanel(colsEl, topo, posBuf, sphereBg || undefined);
10011243 scenes.push(scene);
10021244 const field = await on(rf.final[model.state[k]]);
10031245 fields.push(field);
10041246 bounds.push(finiteRange(field));
1005- const valueBuf = new Float32Array(topo.numVertices);
10061247 fillFieldValues(valueBuf, field, topo);
10071248 valueBufs.push(valueBuf);
1008- colorBufs.push(new Float32Array(topo.numVertices * 3));
1249+ colorBufs.push(colorBuf);
10091250 }
10101251 fileRow = { coords, posBuf, scenes, valueBufs, colorBufs, fields, bounds };
10111252 }
10121253
10131254 // Every panel shares one camera: the study is about the fields, and looking
10141255 // at two of them from different angles is not comparing them.
1015- const all = [...rows.flatMap((r) => r.scenes), ...(fileRow?.scenes ?? [])];
1256+ const all = [
1257+ ...rows.flatMap((r) => [...r.scenes, ...r.diffScenes]),
1258+ ...(fileRow?.scenes ?? []),
1259+ ];
10161260 for (let i = 1; i < all.length; i++) all[0].syncCamerasWith(all[i]);
10171261
10181262 return { rows, fileRow, rangeBars };
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/geom/geometry.tsmodified+28−9View file
@@ -40,8 +40,9 @@ import {
4040 type RuntimeValue,
4141 } from 'numbl-src/numbl-core/runtime/types.ts';
4242 import { ShtPlan } from '../sht/sht.ts';
43-import type { ShtConfig } from '../sht/layout.ts';
43+import { relayoutCoeffs, type ShtConfig } from '../sht/layout.ts';
4444 import type { DerivPlan } from '../sht/deriv.ts';
45+import type { ImportedCoeffs } from './registry.ts';
4546 import { computeMetric, computeFluxMetric } from './metric.ts';
4647 import { toolFiles } from '../tools.ts';
4748 import { inFunction, inModel, ModelCompileError } from '../mgpu/errors.ts';
@@ -61,6 +62,13 @@ export interface GeometryOptions {
6162 params: ModelParams;
6263 /** Computes the theta/phi derivatives the inverse metric quantities need. */
6364 deriv: DerivPlan;
65+ /**
66+ * The surface as coefficients instead of a .m. When set, `source` is never
67+ * evaluated: the coefficients are re-indexed onto this plan's band
68+ * (truncation is exact projection, padding is exact) and everything from
69+ * synthesis onward runs unchanged.
70+ */
71+ coeffs?: ImportedCoeffs;
6472 }
6573
6674 export class Geometry {
@@ -189,16 +197,27 @@ export class Geometry {
189197 const { sht, cfg, source, paramNames, params, deriv } = opts;
190198 const npts = cfg.nlat * cfg.nphi;
191199
192- const { theta, phi } = gridAngles(sht, cfg);
193- const raw = evaluateShape(source, paramNames, params, theta, phi, npts);
194-
195200 // Coefficients first, then back to the grid: what the solver and the
196201 // renderer both see is the band-limited surface, not the raw .m output.
197- const [X, Y, Z] = [
198- await sht.analys(raw[0]),
199- await sht.analys(raw[1]),
200- await sht.analys(raw[2]),
201- ];
202+ // An imported surface arrives already as coefficients, so its .m stage is
203+ // a re-indexing onto this plan's band rather than an evaluation.
204+ let X: Float32Array;
205+ let Y: Float32Array;
206+ let Z: Float32Array;
207+ if (opts.coeffs) {
208+ const from = { lmax: opts.coeffs.lmax, mmax: opts.coeffs.mmax };
209+ X = relayoutCoeffs(opts.coeffs.X, from, cfg);
210+ Y = relayoutCoeffs(opts.coeffs.Y, from, cfg);
211+ Z = relayoutCoeffs(opts.coeffs.Z, from, cfg);
212+ } else {
213+ const { theta, phi } = gridAngles(sht, cfg);
214+ const raw = evaluateShape(source, paramNames, params, theta, phi, npts);
215+ [X, Y, Z] = [
216+ await sht.analys(raw[0]),
217+ await sht.analys(raw[1]),
218+ await sht.analys(raw[2]),
219+ ];
220+ }
202221 const [x, y, z] = [
203222 await sht.synth(X),
204223 await sht.synth(Y),
src/geom/registry.tsmodified+46−0View file
@@ -23,6 +23,23 @@ export interface MGeometry {
2323 params: ParamSpec[];
2424 /** MATLAB source — the shape itself. */
2525 source: string;
26+ /**
27+ * When set, the shape IS these spherical-harmonic coefficients and `source`
28+ * is a display-only stub: Geometry.create re-indexes them onto the session's
29+ * band instead of evaluating any .m. How a surface handed over by another
30+ * app (reharm's "Export to turing-surface") enters the pipeline.
31+ */
32+ coeffs?: ImportedCoeffs;
33+}
34+
35+/** A surface as coefficients of x, y, z — the same three arrays a reference
36+ * file's /geometry group holds, in the shared SHTNS layout (src/sht/layout.ts). */
37+export interface ImportedCoeffs {
38+ lmax: number;
39+ mmax: number;
40+ X: Float32Array;
41+ Y: Float32Array;
42+ Z: Float32Array;
2643 }
2744
2845 const sphere: MGeometry = {
@@ -92,6 +109,35 @@ export const defaultGeometryParams = (g: MGeometry): Params =>
92109 * none gets: the case where the solver is exact. */
93110 export const SPHERE_KEY = 'sphere';
94111
112+/** The dropdown key an imported surface lives under. Not in `mGeometries`:
113+ * the entry exists only after an import arrives, and the page owns it. */
114+export const IMPORTED_GEOMETRY_KEY = 'imported';
115+
116+/**
117+ * Wrap a handed-over surface as a geometry the rest of the app can hold like
118+ * any registry entry. No parameters — the shape is whatever the coefficients
119+ * say — and the .m slot carries a note rather than code, since the import
120+ * path in Geometry.create never evaluates it.
121+ */
122+export function makeImportedGeometry(
123+ name: string,
124+ blurb: string,
125+ coeffs: ImportedCoeffs,
126+): MGeometry {
127+ return {
128+ key: IMPORTED_GEOMETRY_KEY,
129+ label: name,
130+ blurb,
131+ params: [],
132+ source:
133+ `% ${name}\n` +
134+ `% This surface was handed over as spherical-harmonic coefficients\n` +
135+ `% (lmax ${coeffs.lmax}), so there is no .m to edit: the solver uses the\n` +
136+ `% coefficients directly, and edits here have no effect.\n`,
137+ coeffs,
138+ };
139+}
140+
95141 /** What the app and the benchmark start on. Not the sphere: this project
96142 * exists for the other shapes, and opening on the reference case would hide
97143 * the one thing it adds. */
src/main.tsmodified+343−28View file
@@ -1,5 +1,5 @@
11 import { requestShtDevice, describeAdapter } from './sht/sht.ts';
2-import { gridForLmax } from './sht/layout.ts';
2+import { gridForLmax, nlmCalc } from './sht/layout.ts';
33 import { ModelSession } from './mgpu/session.ts';
44 import { mModelByKey, presets, type MModel, type Params } from './mgpu/registry.ts';
55 import { ModelCompileError, formatFailure } from './mgpu/errors.ts';
@@ -18,6 +18,10 @@ import {
1818 mGeometryByKey,
1919 defaultGeometryParams,
2020 DEFAULT_GEOMETRY_KEY,
21+ SPHERE_KEY,
22+ IMPORTED_GEOMETRY_KEY,
23+ makeImportedGeometry,
24+ type ImportedCoeffs,
2125 type MGeometry,
2226 } from './geom/registry.ts';
2327 import {
@@ -41,6 +45,7 @@ import {
4145 } from './compare/variants.ts';
4246 import { loadReferenceFile } from './compare/referenceFile.ts';
4347 import type { ReferenceCase } from './compare/referenceCase.ts';
48+import { generateMatlabScript, MATLAB_SCRIPT_NAME } from './export/matlabScript.ts';
4449
4550 const $ = <T extends HTMLElement>(id: string): T =>
4651 document.getElementById(id) as T;
@@ -66,12 +71,15 @@ const elMovieRotate = $<HTMLInputElement>('movierotate');
6671 const elMovie = $<HTMLButtonElement>('movie');
6772 const elModeSimulate = $<HTMLButtonElement>('mode-simulate');
6873 const elModeEffort = $<HTMLButtonElement>('mode-effort');
74+const elModeVsSphere = $<HTMLButtonElement>('mode-vs-sphere');
6975 const elModeVsUpload = $<HTMLButtonElement>('mode-vs-upload');
7076 const elModeDesc = $('mode-desc');
7177 const elCompareBar = $('comparebar');
78+const elCmpAxes = $('cmp-axes');
7279 const elCmpNiter = $('cmp-niter');
7380 const elCmpLmax = $('cmp-lmax');
7481 const elCmpDt = $('cmp-dt');
82+const elCmpRefLabel = $('cmp-reflabel');
7583 const elCmpRef = $<HTMLSelectElement>('cmp-ref');
7684 const elCmpFile = $<HTMLInputElement>('cmp-file');
7785 const elCmpFileInfo = $('cmp-fileinfo');
@@ -82,10 +90,15 @@ const elParams = $('params');
8290 const elGeomParams = $('geomparams');
8391 const elGeomNote = $('geomnote');
8492 const elPanels = $('panels');
93+const elCmpChart = $('cmp-chart');
8594 const elStats = $('stats');
8695 const elBenchResult = $('benchresult');
8796 const elCmd = $('cmd');
8897 const elCopyCmd = $<HTMLButtonElement>('copycmd');
98+const elMatlab = $<HTMLDetailsElement>('matlab');
99+const elMatlabScript = $('matlabscript');
100+const elCopyMatlab = $<HTMLButtonElement>('copymatlab');
101+const elDownloadMatlab = $<HTMLButtonElement>('downloadmatlab');
89102 const elBlurb = $('blurb');
90103 const elErr = $('err');
91104 const elSource = $<HTMLTextAreaElement>('source');
@@ -278,6 +291,9 @@ let posBuf: Float32Array | null = null;
278291 * mode. While it is non-null there is no `session`: the study owns one per
279292 * variant, and the panels area is its grid. */
280293 let compareRun: CompareRun | null = null;
294+/** The surface handed over by another app (reharm), when one has arrived —
295+ * what the geometry dropdown's 'imported' entry resolves to. */
296+let importedGeometry: MGeometry | null = null;
281297 /** `session`'s spectral state as of the last (re-)seed — what "Restart"
282298 * rewinds to. Captured fresh each time a new field is actually established
283299 * (rebuild/reseed), not just once, so Restart reflects the run's current
@@ -349,6 +365,7 @@ function buildGeomParamInputs(): void {
349365 next = spec.min + Math.floor(Math.random() * (span + 1));
350366 }
351367 geomParams[spec.key] = next;
368+ updateCommand();
352369 viewChange = viewChange.then(() => applyGeometry());
353370 });
354371 elGeomParams.append(button);
@@ -365,6 +382,7 @@ function buildGeomParamInputs(): void {
365382 input.addEventListener('change', () => {
366383 const v = Number(input.value);
367384 if (Number.isFinite(v)) geomParams[spec.key] = v;
385+ updateCommand();
368386 viewChange = viewChange.then(() => applyGeometry());
369387 });
370388 label.append(input);
@@ -389,7 +407,10 @@ function applyPreset(presetKey: string): void {
389407 }
390408
391409 function applyGeometryChoice(key: string): void {
392- const next = mGeometryByKey(key);
410+ // The imported entry is the page's own, not the registry's — it exists only
411+ // after a surface has been handed over (see the reharm import block below).
412+ const next =
413+ key === IMPORTED_GEOMETRY_KEY ? importedGeometry : mGeometryByKey(key);
393414 if (!next) {
394415 elErr.textContent = `No .m geometry for '${key}'`;
395416 return;
@@ -399,6 +420,8 @@ function applyGeometryChoice(key: string): void {
399420 editedGeomSource = null;
400421 buildGeomParamInputs();
401422 showEditorFile();
423+ // The command line and the MATLAB export both bake the surface in.
424+ updateCommand();
402425 }
403426
404427 /** Load the chosen file into the editor, keeping any unsaved edit to it. */
@@ -440,9 +463,54 @@ function updateCommand(): void {
440463 // equivalent is the ref checker, not the benchmark.
441464 if (compareRun?.refFile) {
442465 elCmd.textContent = `npm run ref -- --in ${compareRun.refFile.label}`;
443- return;
466+ } else if (geometry.coeffs) {
467+ // An imported surface exists only as coefficients in this page, so no
468+ // command line can regenerate the run.
469+ elCmd.textContent = '(imported surface — no desktop equivalent)';
470+ } else {
471+ elCmd.textContent = formatCommand(currentSpec());
472+ }
473+ if (elMatlab.open) refreshMatlabScript();
474+}
475+
476+/** The run on screen as one standalone .m — in a study, its reference
477+ * variant, the same choice `currentSpec` makes. Throws on a working copy
478+ * the export cannot parse (no init/step/shape function). */
479+function matlabScriptText(): string {
480+ if (geometry.coeffs) {
481+ throw new Error(
482+ 'The imported surface has no .m, so the MATLAB export cannot bake it in. ' +
483+ 'Export the surface from reharm as .h5 and load its /geometry ' +
484+ 'coefficients with h5read instead.',
485+ );
486+ }
487+ const spec = currentSpec();
488+ return generateMatlabScript({
489+ model,
490+ modelSource: source(),
491+ params: spec.params,
492+ geometry,
493+ geometrySource: geomSource(),
494+ geometryParams: geomParams,
495+ lmax: spec.lmax,
496+ niter: spec.niter,
497+ lam3: Number(elLam3.value),
498+ seed,
499+ preset: spec.preset,
500+ command: formatCommand(spec),
501+ });
502+}
503+
504+/** Regenerate the visible script text; on failure, show why in its place. */
505+function refreshMatlabScript(): string | null {
506+ try {
507+ const text = matlabScriptText();
508+ elMatlabScript.textContent = text;
509+ return text;
510+ } catch (e) {
511+ elMatlabScript.textContent = e instanceof Error ? e.message : String(e);
512+ return null;
444513 }
445- elCmd.textContent = formatCommand(currentSpec());
446514 }
447515
448516 elModel.addEventListener('change', () => {
@@ -475,6 +543,8 @@ elGeometry.addEventListener('change', () => {
475543 elLam3.addEventListener('change', () => {
476544 const v = Number(elLam3.value);
477545 if (!Number.isFinite(v) || v <= 0) return;
546+ // Not in the bench command, but the MATLAB export bakes lam3 in.
547+ updateCommand();
478548 // Changing the wavelength redraws the field, which restarts the run — so
479549 // pause first, exactly as the Re-seed button does. Without it the reseed's
480550 // readback races the pump's own, and the two collide on the staging buffer.
@@ -578,6 +648,46 @@ elCopyCmd.addEventListener('click', () => {
578648 navigator.clipboard.writeText(text).then(() => flash('Copied'), selectCommand);
579649 });
580650
651+// The MATLAB export: the same run as one self-contained .m. Regenerated from
652+// the current UI state whenever it is shown, copied or downloaded, so the
653+// text always matches the run on screen.
654+elMatlab.addEventListener('toggle', () => {
655+ if (elMatlab.open) refreshMatlabScript();
656+});
657+elCopyMatlab.addEventListener('click', (e) => {
658+ // The buttons live inside the <summary>; without this a click also toggles.
659+ e.preventDefault();
660+ e.stopPropagation();
661+ const text = refreshMatlabScript();
662+ if (text === null || !navigator.clipboard) {
663+ // Open to show the error, or the text to select by hand.
664+ elMatlab.open = true;
665+ return;
666+ }
667+ navigator.clipboard.writeText(text).then(
668+ () => {
669+ elCopyMatlab.textContent = 'Copied';
670+ setTimeout(() => (elCopyMatlab.textContent = 'Copy'), 1200);
671+ },
672+ () => (elMatlab.open = true),
673+ );
674+});
675+elDownloadMatlab.addEventListener('click', (e) => {
676+ e.preventDefault();
677+ e.stopPropagation();
678+ const text = refreshMatlabScript();
679+ if (text === null) {
680+ elMatlab.open = true;
681+ return;
682+ }
683+ const url = URL.createObjectURL(new Blob([text], { type: 'text/x-matlab' }));
684+ const a = document.createElement('a');
685+ a.href = url;
686+ a.download = `${MATLAB_SCRIPT_NAME}.m`;
687+ a.click();
688+ URL.revokeObjectURL(url);
689+});
690+
581691 // ---------------------------------------------------------------- setup
582692 function disposeView(): void {
583693 for (const s of scenes) s.dispose();
@@ -1286,13 +1396,26 @@ function compareRefIndex(): number {
12861396 }
12871397
12881398 function refreshVariants(): void {
1399+ // vs-sphere has no chip grid feeding cmpVariants() — always exactly two
1400+ // fixed rows plus one diff row, regardless of MAX_VARIANTS/MAX_PANELS —
1401+ // so the grid-count logic below doesn't apply here at all.
1402+ if (currentMode === 'vs-sphere') {
1403+ elCmpCount.textContent = `2 rows × ${model.species.length} species + 1 diff row`;
1404+ elCmpCount.style.color = '';
1405+ elCmpStart.disabled = false;
1406+ return;
1407+ }
12891408 const variants = cmpVariants();
12901409 const showDt = cmpSelected.dt.size > 1;
12911410 // With a file loaded the study's model is the file's, and its final state
12921411 // is one more row of panels.
12931412 const cmpModel = refCase?.model ?? model;
12941413 const rowCount = variants.length + (refCase ? 1 : 0);
1295- const panels = rowCount * cmpModel.species.length;
1414+ // Every row but the reference variant (or, against a file, every variant)
1415+ // gets a second diff row underneath it — each one more WebGL context per
1416+ // species, so MAX_PANELS has to bound the real total, not just the values.
1417+ const diffRowCount = refCase ? variants.length : Math.max(0, variants.length - 1);
1418+ const panels = (rowCount + diffRowCount) * cmpModel.species.length;
12961419
12971420 const prev = cmpRefKey;
12981421 elCmpRef.replaceChildren();
@@ -1325,7 +1448,9 @@ function refreshVariants(): void {
13251448 ? `too many: ${tooMany}`
13261449 : `${variants.length} variant${variants.length === 1 ? '' : 's'}` +
13271450 `${refCase ? ' + the file' : ''} × ` +
1328- `${cmpModel.species.length} species = ${panels} panels`;
1451+ `${cmpModel.species.length} species` +
1452+ (diffRowCount > 0 ? ` + ${diffRowCount} diff row${diffRowCount === 1 ? '' : 's'}` : '') +
1453+ ` = ${panels} panels`;
13291454 elCmpCount.style.color = tooMany ? '#b35900' : '';
13301455 elCmpStart.disabled = tooMany !== '' && compareRun === null;
13311456 }
@@ -1368,6 +1493,25 @@ function rebuildLmaxChips(): void {
13681493 buildChips(elCmpLmax, values, cmpSelected.lmax, String);
13691494 }
13701495
1496+/**
1497+ * The four top-level modes and which control groups each shows (see
1498+ * GROUP_NAMES/groupEls above; `.ctrl-group` wrappers in index.html).
1499+ * `currentMode` tracks which configuration is on screen — the compare bar
1500+ * being open, and in which flavor — not whether a study has actually been
1501+ * started inside it. That match matters: without it, opening the bar
1502+ * (which already shows the right groups) leaves its top-row button
1503+ * unhighlighted until a study happens to start, which is inconsistent with
1504+ * `vs-upload`'s one-click flow and reads as broken.
1505+ *
1506+ * Declared here, ahead of the top-level `refreshVariants()` call just below
1507+ * — that call reads `currentMode` (to know whether to skip its chip-grid
1508+ * counting for vs-sphere), so the declaration has to be in scope by the time
1509+ * this file's top-level code actually runs, not merely by the time
1510+ * `refreshVariants` is later invoked from an event handler.
1511+ */
1512+type Mode = 'simulate' | 'compute-effort' | 'vs-sphere' | 'vs-upload';
1513+let currentMode: Mode = 'simulate';
1514+
13711515 rebuildNiterChips();
13721516 rebuildLmaxChips();
13731517 buildChips(elCmpDt, DT_DIVISORS, cmpSelected.dt, (v) => (v === 1 ? 'dt' : `dt/${v}`));
@@ -1401,23 +1545,13 @@ function applyRefUi(): void {
14011545 refreshVariants();
14021546 }
14031547
1404-/**
1405- * The four top-level modes and which control groups each shows (see
1406- * GROUP_NAMES/groupEls above; `.ctrl-group` wrappers in index.html).
1407- * `currentMode` tracks which configuration is on screen — the compare bar
1408- * being open, and in which flavor — not whether a study has actually been
1409- * started inside it. That match matters: without it, opening the bar
1410- * (which already shows the right groups) leaves its top-row button
1411- * unhighlighted until a study happens to start, which is inconsistent with
1412- * `vs-upload`'s one-click flow and reads as broken.
1413- */
1414-type Mode = 'simulate' | 'compute-effort' | 'vs-sphere' | 'vs-upload';
1415-let currentMode: Mode = 'simulate';
1416-
14171548 const MODE_GROUPS: Record<Mode, readonly GroupName[]> = {
14181549 simulate: ['surface', 'surface-params', 'solver', 'display', 'playback', 'benchmark', 'seed', 'movie'],
14191550 'compute-effort': ['surface', 'surface-params', 'display', 'playback', 'seed'],
1420- 'vs-sphere': [], // unreachable — the button is disabled, no listener ever calls setMode with this
1551+ // Unlike compute-effort, there is no separate chip grid for this mode —
1552+ // both rows always mirror whatever niter/lmax the Simulate controls have,
1553+ // so `solver` has to stay visible; it's the only place those get set.
1554+ 'vs-sphere': ['surface', 'surface-params', 'solver', 'display', 'playback', 'seed'],
14211555 // No `seed` here: nothing in that group does anything useful against a
14221556 // loaded file (lam3 is silently absorbed, and Restart already covers what
14231557 // Re-seed would otherwise be doing — reloading the file's fixed initial
@@ -1432,7 +1566,11 @@ const MODE_DESCRIPTIONS: Record<Mode, string> = {
14321566 'When we change the computational effort of the solver by varying solve iterations, lmax, or timestep, ' +
14331567 'how does the solution change? Find out by running several ' +
14341568 'so you can see how each setting trades accuracy for speed.',
1435- 'vs-sphere': '',
1569+ 'vs-sphere':
1570+ 'Run this model once on the selected geometry and once on the plain unit sphere, from the exact same ' +
1571+ 'starting state and the same solve iterations, lmax and timestep, so geometry is the only thing that ' +
1572+ 'differs. See how much resolving the true shape actually changes the pattern, versus approximating it ' +
1573+ 'as a sphere.',
14361574 'vs-upload':
14371575 'Load a saved reference run (an .h5 file) and run this solver to the ' +
14381576 'same physical end time from the same initial condition, to check how ' +
@@ -1443,6 +1581,7 @@ const MODE_DESCRIPTIONS: Record<Mode, string> = {
14431581 function setModeButtons(mode: Mode): void {
14441582 elModeSimulate.setAttribute('aria-pressed', String(mode === 'simulate'));
14451583 elModeEffort.setAttribute('aria-pressed', String(mode === 'compute-effort'));
1584+ elModeVsSphere.setAttribute('aria-pressed', String(mode === 'vs-sphere'));
14461585 elModeVsUpload.setAttribute('aria-pressed', String(mode === 'vs-upload'));
14471586 elModeDesc.textContent = MODE_DESCRIPTIONS[mode];
14481587 }
@@ -1467,26 +1606,33 @@ function applyModeVisibility(mode: Mode): void {
14671606 function enterMode(mode: Mode): void {
14681607 applyModeVisibility(mode);
14691608 elCompareBar.hidden = mode === 'simulate';
1609+ // vs-sphere has no chip grid to pick from — both rows mirror the Simulate
1610+ // panel's own niter/lmax — and no reference to pick among variants either,
1611+ // since the reference is always "the selected geometry." Just the
1612+ // description and the Compile button apply.
1613+ elCmpAxes.hidden = mode === 'vs-sphere';
1614+ elCmpRefLabel.hidden = mode === 'vs-sphere';
14701615 }
14711616
14721617 /** Entering a mode from the top row. */
14731618 function setMode(mode: Mode): void {
1474- if (mode === 'vs-sphere') return; // unreachable — button is disabled
14751619 if (mode === 'simulate') {
14761620 if (compareRun) void stopCompare();
14771621 enterMode('simulate');
14781622 return;
14791623 }
1480- if (mode === 'compute-effort') {
1624+ if (mode === 'compute-effort' || mode === 'vs-sphere') {
14811625 // Tear down whatever study is running first (mirrors Simulate above) —
14821626 // stopCompare's synchronous prefix disposes it and nulls `compareRun`
14831627 // before its first `await`, so `refCase` is safe to drop right after.
1628+ // vs-sphere never uses a reference file either — its "reference" is
1629+ // always the selected geometry — so the same drop applies there too.
14841630 if (compareRun) void stopCompare();
14851631 if (refCase) {
14861632 refCase = null;
14871633 applyRefUi();
14881634 }
1489- enterMode('compute-effort');
1635+ enterMode(mode);
14901636 return;
14911637 }
14921638 // vs-upload: opens the file picker; entering the mode itself happens once
@@ -1497,6 +1643,7 @@ function setMode(mode: Mode): void {
14971643
14981644 elModeSimulate.addEventListener('click', () => setMode('simulate'));
14991645 elModeEffort.addEventListener('click', () => setMode('compute-effort'));
1646+elModeVsSphere.addEventListener('click', () => setMode('vs-sphere'));
15001647 elModeVsUpload.addEventListener('click', () => setMode('vs-upload'));
15011648
15021649 elCmpFile.addEventListener('change', () => {
@@ -1577,11 +1724,42 @@ async function startCompare(): Promise<void> {
15771724 if (compareRun || !device) return;
15781725 // Snapshotted for the whole study: `refCase` only changes with no study up
15791726 // (clearing is disabled during one, and loading tears it down first).
1580- const rc = refCase;
1727+ // vs-sphere never has one — its "reference" is always the selected
1728+ // geometry, not a file.
1729+ const rc = currentMode === 'vs-sphere' ? null : refCase;
15811730 const cmpModel = rc?.model ?? model;
1582- const variants = cmpVariants();
1731+
1732+ // vs-sphere: exactly two rows — the selected geometry (the reference) and
1733+ // the plain unit sphere — mirroring whatever niter/lmax the Simulate
1734+ // controls have, so geometry is the only thing that differs between them.
1735+ // Every other mode still drives its rows from the chip grid.
1736+ const isVsSphere = currentMode === 'vs-sphere';
1737+ let variants: Variant[];
1738+ let reference: number;
1739+ let geometries: { geometry: MGeometry; geometryParams: Params; geometrySource: string }[] | undefined;
1740+ let renderOnReferenceGeometry: boolean | undefined;
1741+ let rowLabels: string[] | undefined;
1742+ if (isVsSphere) {
1743+ const niter = Number(elNiter.value);
1744+ const lmax = Number(elLmax.value);
1745+ variants = [{ niter, lmax, dtDiv: 1 }, { niter, lmax, dtDiv: 1 }];
1746+ reference = 0;
1747+ const sphereGeom = mGeometryByKey(SPHERE_KEY)!;
1748+ geometries = [
1749+ { geometry, geometryParams: geomParams, geometrySource: geomSource() },
1750+ { geometry: sphereGeom, geometryParams: defaultGeometryParams(sphereGeom), geometrySource: sphereGeom.source },
1751+ ];
1752+ renderOnReferenceGeometry = true;
1753+ rowLabels = [geometry.label, 'unit sphere'];
1754+ } else {
1755+ variants = cmpVariants();
1756+ reference = rc ? 0 : compareRefIndex();
1757+ }
1758+
15831759 const rowCount = variants.length + (rc ? 1 : 0);
1584- if (variants.length > MAX_VARIANTS || rowCount * cmpModel.species.length > MAX_PANELS) {
1760+ const diffRowCount = rc ? variants.length : Math.max(0, variants.length - 1);
1761+ const panels = (rowCount + diffRowCount) * cmpModel.species.length;
1762+ if (!isVsSphere && (variants.length > MAX_VARIANTS || panels > MAX_PANELS)) {
15851763 return;
15861764 }
15871765 // Take down the single run first: its pump, its scenes, its session. The
@@ -1609,7 +1787,10 @@ async function startCompare(): Promise<void> {
16091787 geometryParams: rc ? rc.geometryParams : geomParams,
16101788 geometrySource: rc ? rc.geometry.source : geomSource(),
16111789 variants,
1612- reference: rc ? 0 : compareRefIndex(),
1790+ reference,
1791+ geometries,
1792+ renderOnReferenceGeometry,
1793+ rowLabels,
16131794 refFile: rc ?? undefined,
16141795 onFinished: () => setRunning(false),
16151796 seed,
@@ -1617,6 +1798,7 @@ async function startCompare(): Promise<void> {
16171798 morph,
16181799 colormapName: () => elColormap.value,
16191800 container: elPanels,
1801+ chartContainer: elCmpChart,
16201802 onStatus: (html) => (elStats.innerHTML = html),
16211803 });
16221804 } catch (e) {
@@ -1653,6 +1835,118 @@ async function rebuildCompare(): Promise<void> {
16531835 await startCompare();
16541836 }
16551837
1838+// ------------------------------------------------------ import from reharm
1839+// reharm's "Export to turing-surface" button opens this page with ?import and
1840+// hands the fitted surface over by postMessage — the two apps live on
1841+// different origins, so a message between opener and opened tab is the only
1842+// channel, and structured clone carries the Float32Arrays without encoding.
1843+// Protocol: the opener pings {type:'reharm-import-ping'} until this page
1844+// answers {type:'reharm-import-ready'} (only once boot has a live session to
1845+// swap the surface into), then posts {type:'reharm-geometry', version, lmax,
1846+// mmax, Gx, Gy, Gz, name, provenance}; this page validates, installs it as
1847+// the dropdown's 'imported' entry, and acks. The payload is data, not code:
1848+// three coefficient arrays and label strings, the strings stripped of markup
1849+// before they reach any innerHTML (updateGeomNote).
1850+
1851+/** null until boot resolves; then '' when imports can proceed, or the reason
1852+ * they cannot (no WebGPU). */
1853+let importGate: string | null = null;
1854+
1855+const cleanLabel = (v: unknown, max: number): string =>
1856+ String(v ?? '')
1857+ .replace(/[<>&"'`]/g, ' ')
1858+ .replace(/\s+/g, ' ')
1859+ .trim()
1860+ .slice(0, max);
1861+
1862+/** Validate a handed-over payload into coefficients, or say what is wrong. */
1863+function importedCoeffsOf(d: Record<string, unknown>): ImportedCoeffs | string {
1864+ const version = Number(d.version);
1865+ if (version !== 1) return `unknown import payload version '${String(d.version)}'`;
1866+ const lmax = Number(d.lmax);
1867+ const mmax = Number(d.mmax);
1868+ if (!Number.isInteger(lmax) || lmax < 1 || lmax > 4096) return `bad lmax '${String(d.lmax)}'`;
1869+ if (!Number.isInteger(mmax) || mmax < 0 || mmax > lmax) return `bad mmax '${String(d.mmax)}'`;
1870+ const nlm = nlmCalc(lmax, mmax);
1871+ const arrays: Float32Array[] = [];
1872+ for (const name of ['Gx', 'Gy', 'Gz']) {
1873+ const q = d[name];
1874+ if (!(q instanceof Float32Array)) return `'${name}' is not a Float32Array`;
1875+ if (q.length !== 2 * nlm) return `'${name}' has ${q.length} values, expected 2*nlm = ${2 * nlm}`;
1876+ if (!q.every(Number.isFinite)) return `'${name}' has non-finite values`;
1877+ arrays.push(q);
1878+ }
1879+ return { lmax, mmax, X: arrays[0], Y: arrays[1], Z: arrays[2] };
1880+}
1881+
1882+/** Adopt the surface: it becomes (or replaces) the dropdown's imported entry
1883+ * and is selected. The run also adopts the surface's own band — truncating a
1884+ * fit to a coarser default would silently smooth the shape that was just
1885+ * exported — so the lmax dropdown gains the exact value when it is not a
1886+ * stock one, and a band change goes through the full rebuild (new grid)
1887+ * rather than the in-place swap. */
1888+function installImportedGeometry(g: MGeometry): void {
1889+ importedGeometry = g;
1890+ let opt = elGeometry.querySelector<HTMLOptionElement>(
1891+ `option[value="${IMPORTED_GEOMETRY_KEY}"]`,
1892+ );
1893+ if (!opt) {
1894+ opt = document.createElement('option');
1895+ opt.value = IMPORTED_GEOMETRY_KEY;
1896+ elGeometry.append(opt);
1897+ }
1898+ opt.textContent = g.label;
1899+ elGeometry.value = IMPORTED_GEOMETRY_KEY;
1900+ applyGeometryChoice(IMPORTED_GEOMETRY_KEY);
1901+ const lmax = g.coeffs!.lmax;
1902+ if (Number(elLmax.value) !== lmax) {
1903+ const options = [...elLmax.options];
1904+ if (!options.some((o) => Number(o.value) === lmax)) {
1905+ const o = document.createElement('option');
1906+ o.value = String(lmax);
1907+ o.textContent = String(lmax);
1908+ elLmax.add(o, options.find((x) => Number(x.value) > lmax) ?? null);
1909+ }
1910+ elLmax.value = String(lmax);
1911+ void rebuild();
1912+ } else {
1913+ viewChange = viewChange.then(() => applyGeometry());
1914+ }
1915+}
1916+
1917+window.addEventListener('message', (ev: MessageEvent) => {
1918+ const d = ev.data as Record<string, unknown> | null;
1919+ if (!d || typeof d !== 'object' || typeof d.type !== 'string') return;
1920+ const reply = (msg: Record<string, unknown>): void => {
1921+ (ev.source as Window | null)?.postMessage(msg, { targetOrigin: ev.origin });
1922+ };
1923+ if (d.type === 'reharm-import-ping') {
1924+ // Silence before boot resolves is fine — the opener keeps pinging.
1925+ if (importGate === '') reply({ type: 'reharm-import-ready' });
1926+ else if (importGate) reply({ type: 'reharm-import-error', message: importGate });
1927+ return;
1928+ }
1929+ if (d.type !== 'reharm-geometry') return;
1930+ if (importGate !== '') {
1931+ reply({ type: 'reharm-import-error', message: importGate ?? 'the page is still starting up' });
1932+ return;
1933+ }
1934+ const coeffs = importedCoeffsOf(d);
1935+ if (typeof coeffs === 'string') {
1936+ reply({ type: 'reharm-import-error', message: coeffs });
1937+ return;
1938+ }
1939+ const p = (d.provenance ?? {}) as Record<string, unknown>;
1940+ const name = `imported: ${cleanLabel(d.name, 40) || 'surface'}`;
1941+ const blurb =
1942+ `Handed over by reharm (lmax ${coeffs.lmax}` +
1943+ (p.map ? `, map ${cleanLabel(p.map, 60)}` : '') +
1944+ (p.filter ? `, filter ${cleanLabel(p.filter, 60)}` : '') +
1945+ `).`;
1946+ installImportedGeometry(makeImportedGeometry(name, blurb, coeffs));
1947+ reply({ type: 'reharm-import-ack' });
1948+});
1949+
16561950 // ---------------------------------------------------------------- boot
16571951 async function boot(): Promise<void> {
16581952 enterMode('simulate');
@@ -1673,6 +1967,8 @@ async function boot(): Promise<void> {
16731967 elErr.textContent =
16741968 `WebGPU is not available (${e instanceof Error ? e.message : e}). ` +
16751969 `Use a WebGPU-capable browser such as Chrome or Edge.`;
1970+ importGate = 'WebGPU is not available in the opened tab';
1971+ announceImportGate();
16761972 return;
16771973 }
16781974 device.lost.then((info) => {
@@ -1681,6 +1977,25 @@ async function boot(): Promise<void> {
16811977 }
16821978 });
16831979 await rebuild();
1980+ // Imports wait for a live run to swap the surface into; a tab opened by
1981+ // reharm's export button (?import) is told without waiting to be pinged.
1982+ importGate = '';
1983+ announceImportGate();
1984+}
1985+
1986+/** Tell the opener that navigated here with ?import how imports stand. Only
1987+ * a signal — no data — so the opener's origin need not be known. */
1988+function announceImportGate(): void {
1989+ if (!window.opener || !new URLSearchParams(location.search).has('import')) return;
1990+ const msg =
1991+ importGate === ''
1992+ ? { type: 'reharm-import-ready' }
1993+ : { type: 'reharm-import-error', message: importGate };
1994+ try {
1995+ (window.opener as Window).postMessage(msg, { targetOrigin: '*' });
1996+ } catch {
1997+ // The opener may already be gone; the ping path covers the rest.
1998+ }
16841999 }
16852000
16862001 void boot();
src/mgpu/session.tsmodified+2−0View file
@@ -155,6 +155,7 @@ export class ModelSession {
155155 paramNames: geometryModel.params.map((p) => p.key),
156156 params: geometryParams,
157157 deriv,
158+ coeffs: geometryModel.coeffs,
158159 });
159160 const gpu = await GpuModel.create({
160161 device,
@@ -211,6 +212,7 @@ export class ModelSession {
211212 paramNames: geometryModel.params.map((p) => p.key),
212213 params,
213214 deriv: this.#deriv,
215+ coeffs: geometryModel.coeffs,
214216 });
215217 this.#geometry = next;
216218 this.#geometryModel = geometryModel;
src/render/errorChart.tsadded+328−0View file
@@ -0,0 +1,328 @@
1+/**
2+ * Live Δ(t) chart for a comparison study: one canvas, one line per
3+ * (variant, species) pair, sharing a log-scaled y-axis and an x-axis in
4+ * model time `t`. A variant's lines are colored with that variant's own
5+ * accent color (the same color shown on the left edge of its row of
6+ * panels); species are told apart within one color by line style (solid,
7+ * dashed, dotted, …), matched by the legend's dash swatches.
8+ *
9+ * Deliberately dumb — no library, just a canvas redrawn from scratch on every
10+ * push(). History is unbounded (a study is a "does it converge" question, and
11+ * capping the window would hide exactly the slow drift that question is
12+ * about); only the running log-min/max is maintained incrementally so a
13+ * redraw stays O(points on screen), not O(history) beyond what it already
14+ * draws.
15+ */
16+
17+import { fmtValue } from './colorbar.ts';
18+
19+export interface ErrorChartRow {
20+ label: string;
21+ color: string;
22+}
23+
24+/** Floor for the log scale — well below any error this app ever measures,
25+ * just far enough that log10 never sees zero or a negative. */
26+const MIN_LOG_VALUE = 1e-12;
27+
28+/** Canvas line-dash patterns by species index — solid, dashed, dotted,
29+ * repeating for a fourth-plus species (no shipped model has one). Mirrored
30+ * in index.html's `.cmp-chart-dash-{k}` border-style classes for the
31+ * legend swatches. */
32+const DASH_PATTERNS: number[][] = [[], [7, 4], [2, 3]];
33+
34+const MARGIN_LEFT = 44;
35+const MARGIN_BOTTOM = 38;
36+const MARGIN_TOP = 8;
37+const MARGIN_RIGHT = 10;
38+
39+export class ErrorChart {
40+ #container: HTMLElement;
41+ #species: string[];
42+ #rows: ErrorChartRow[];
43+ #canvas: HTMLCanvasElement | null = null;
44+
45+ #ts: number[] = [];
46+ /** #errs[k][i] is one species' one row's history, same length as #ts. */
47+ #errs: number[][][] = [];
48+ /** Running log10 extent across every (species, row) pushed — one shared
49+ * y-axis now, not one per species. */
50+ #minLog = Infinity;
51+ #maxLog = -Infinity;
52+
53+ constructor(container: HTMLElement, species: string[], rows: ErrorChartRow[]) {
54+ this.#container = container;
55+ this.#species = species;
56+ this.#rows = rows;
57+ this.#resetSeries();
58+ this.#build();
59+ }
60+
61+ #resetSeries(): void {
62+ this.#ts = [];
63+ this.#errs = this.#species.map(() => this.#rows.map(() => []));
64+ this.#minLog = Infinity;
65+ this.#maxLog = -Infinity;
66+ }
67+
68+ #build(): void {
69+ this.#container.replaceChildren();
70+ this.#canvas = null;
71+ if (this.#rows.length === 0 || this.#species.length === 0) return;
72+
73+ const legend = document.createElement('div');
74+ legend.className = 'cmp-chart-legend';
75+
76+ const variantGroup = document.createElement('div');
77+ variantGroup.className = 'cmp-chart-legend-group';
78+ for (const r of this.#rows) {
79+ const item = document.createElement('span');
80+ item.className = 'cmp-chart-legend-item';
81+ const swatch = document.createElement('i');
82+ swatch.style.background = r.color;
83+ item.append(swatch, document.createTextNode(r.label));
84+ variantGroup.append(item);
85+ }
86+
87+ const speciesGroup = document.createElement('div');
88+ speciesGroup.className = 'cmp-chart-legend-group cmp-chart-legend-species';
89+ this.#species.forEach((name, k) => {
90+ const item = document.createElement('span');
91+ item.className = 'cmp-chart-legend-item';
92+ const swatch = document.createElement('i');
93+ swatch.className = `cmp-chart-dash cmp-chart-dash-${k % DASH_PATTERNS.length}`;
94+ item.append(swatch, document.createTextNode(name));
95+ speciesGroup.append(item);
96+ });
97+
98+ legend.append(variantGroup, speciesGroup);
99+ this.#container.append(legend);
100+
101+ const canvas = document.createElement('canvas');
102+ canvas.className = 'cmp-chart-canvas';
103+ this.#container.append(canvas);
104+ this.#canvas = canvas;
105+ }
106+
107+ /** One frame's sample: `t` shared by every row, `perRowErr[i][k]` the
108+ * relative-L2 error of row i's species k against the reference. Rows here
109+ * are exactly the ones passed to the constructor, in the same order — the
110+ * reference row is never included (its Δ against itself is always 0). */
111+ push(t: number, perRowErr: number[][]): void {
112+ this.#ts.push(t);
113+ for (let k = 0; k < this.#species.length; k++) {
114+ for (let i = 0; i < this.#rows.length; i++) {
115+ const v = perRowErr[i]?.[k];
116+ this.#errs[k][i].push(v === undefined ? NaN : v);
117+ if (Number.isFinite(v) && v! > 0) {
118+ const lv = Math.log10(v!);
119+ if (lv < this.#minLog) this.#minLog = lv;
120+ if (lv > this.#maxLog) this.#maxLog = lv;
121+ }
122+ }
123+ }
124+ this.#draw();
125+ }
126+
127+ /** Back to no history — called wherever the study's clock itself resets
128+ * (restart, reseed), so the chart never shows a curve spanning a rewind. */
129+ reset(): void {
130+ this.#resetSeries();
131+ this.#draw();
132+ }
133+
134+ /** Redraw with the data as it stands — for a container resize, where
135+ * nothing new has been measured but the canvas backing buffer has to
136+ * change to match. */
137+ redraw(): void {
138+ this.#draw();
139+ }
140+
141+ dispose(): void {
142+ this.#container.replaceChildren();
143+ this.#canvas = null;
144+ }
145+
146+ #draw(): void {
147+ const canvas = this.#canvas;
148+ if (!canvas) return;
149+ const rect = canvas.getBoundingClientRect();
150+ const dpr = window.devicePixelRatio || 1;
151+ const cssW = Math.max(1, rect.width);
152+ const cssH = Math.max(1, rect.height || 240);
153+ const pxW = Math.max(1, Math.round(cssW * dpr));
154+ const pxH = Math.max(1, Math.round(cssH * dpr));
155+ if (canvas.width !== pxW) canvas.width = pxW;
156+ if (canvas.height !== pxH) canvas.height = pxH;
157+ const ctx = canvas.getContext('2d');
158+ if (!ctx) return;
159+ // Draw in CSS-pixel coordinates throughout; the transform alone accounts
160+ // for devicePixelRatio, rather than scaling every margin/font by hand.
161+ ctx.setTransform(dpr, 0, 0, dpr, 0, 0);
162+ ctx.clearRect(0, 0, cssW, cssH);
163+
164+ const n = this.#ts.length;
165+ if (n < 2) return;
166+
167+ const rawLo = this.#minLog;
168+ const rawHi = this.#maxLog;
169+ if (!Number.isFinite(rawLo) || !Number.isFinite(rawHi)) return;
170+ let lo = rawLo;
171+ let hi = rawHi;
172+ if (hi - lo < 1e-6) {
173+ lo -= 0.5;
174+ hi += 0.5;
175+ }
176+ const pad = (hi - lo) * 0.08;
177+ lo -= pad;
178+ hi += pad;
179+
180+ const style = getComputedStyle(document.documentElement);
181+ const inkColor = style.getPropertyValue('--ink-2').trim() || '#666';
182+ const gridColor = style.getPropertyValue('--line').trim() || '#ccc';
183+
184+ const plotX = MARGIN_LEFT;
185+ const plotY = MARGIN_TOP;
186+ const plotW = Math.max(1, cssW - MARGIN_LEFT - MARGIN_RIGHT);
187+ const plotH = Math.max(1, cssH - MARGIN_TOP - MARGIN_BOTTOM);
188+
189+ const t0 = this.#ts[0];
190+ const t1 = this.#ts[n - 1];
191+ const tSpan = t1 - t0 || 1;
192+ const xAt = (t: number): number => plotX + ((t - t0) / tSpan) * plotW;
193+ const yAt = (v: number): number => {
194+ const lv = Math.max(lo, Math.min(hi, Math.log10(Math.max(v, MIN_LOG_VALUE))));
195+ return plotY + (1 - (lv - lo) / (hi - lo)) * plotH;
196+ };
197+
198+ ctx.font = '11px sans-serif';
199+
200+ // ---- y-axis: decade gridlines + ticks + labels -----------------------
201+ const yTicks = decadeTicks(lo, hi, rawLo, rawHi);
202+ ctx.textAlign = 'right';
203+ ctx.textBaseline = 'middle';
204+ for (const tick of yTicks) {
205+ const y = yAt(tick.value);
206+ ctx.strokeStyle = gridColor;
207+ ctx.globalAlpha = 0.3;
208+ ctx.lineWidth = 1;
209+ ctx.beginPath();
210+ ctx.moveTo(plotX, y);
211+ ctx.lineTo(plotX + plotW, y);
212+ ctx.stroke();
213+ ctx.globalAlpha = 1;
214+ ctx.beginPath();
215+ ctx.moveTo(plotX - 4, y);
216+ ctx.lineTo(plotX, y);
217+ ctx.stroke();
218+ ctx.fillStyle = inkColor;
219+ ctx.fillText(tick.label, plotX - 6, y);
220+ }
221+
222+ // ---- x-axis: ticks + labels --------------------------------------------
223+ const xTicks = evenTicks(t0, t1, 5);
224+ ctx.textAlign = 'center';
225+ ctx.textBaseline = 'top';
226+ ctx.strokeStyle = gridColor;
227+ for (const t of xTicks) {
228+ const x = xAt(t);
229+ ctx.beginPath();
230+ ctx.moveTo(x, plotY + plotH);
231+ ctx.lineTo(x, plotY + plotH + 4);
232+ ctx.stroke();
233+ ctx.fillStyle = inkColor;
234+ ctx.fillText(fmtValue(t), x, plotY + plotH + 6);
235+ }
236+
237+ // ---- axis frame ---------------------------------------------------------
238+ ctx.strokeStyle = gridColor;
239+ ctx.beginPath();
240+ ctx.moveTo(plotX, plotY);
241+ ctx.lineTo(plotX, plotY + plotH);
242+ ctx.lineTo(plotX + plotW, plotY + plotH);
243+ ctx.stroke();
244+
245+ // ---- axis titles ---------------------------------------------------------
246+ ctx.fillStyle = inkColor;
247+ ctx.textAlign = 'center';
248+ ctx.textBaseline = 'bottom';
249+ ctx.fillText('t', plotX + plotW / 2, cssH - 2);
250+
251+ ctx.save();
252+ ctx.translate(12, plotY + plotH / 2);
253+ ctx.rotate(-Math.PI / 2);
254+ ctx.textAlign = 'center';
255+ ctx.textBaseline = 'alphabetic';
256+ ctx.fillText('relative L2 error', 0, 0);
257+ ctx.restore();
258+
259+ // ---- the data: one line per (species, row) -----------------------------
260+ for (let k = 0; k < this.#species.length; k++) {
261+ const dash = DASH_PATTERNS[k % DASH_PATTERNS.length];
262+ for (let i = 0; i < this.#rows.length; i++) {
263+ const series = this.#errs[k][i];
264+ ctx.setLineDash(dash);
265+ ctx.strokeStyle = this.#rows[i].color;
266+ ctx.lineWidth = 1.5;
267+ ctx.beginPath();
268+ let started = false;
269+ for (let j = 0; j < n; j++) {
270+ const v = series[j];
271+ // A non-finite or non-positive sample (a diverged row, or a norm
272+ // with zero denominator) breaks the line rather than drawing a
273+ // spurious segment through a point that has no place on a log axis.
274+ if (!Number.isFinite(v) || v <= 0) {
275+ started = false;
276+ continue;
277+ }
278+ const x = xAt(this.#ts[j]);
279+ const y = yAt(v);
280+ if (started) ctx.lineTo(x, y);
281+ else ctx.moveTo(x, y);
282+ started = true;
283+ }
284+ ctx.stroke();
285+ }
286+ }
287+ ctx.setLineDash([]);
288+ }
289+}
290+
291+/**
292+ * Y-axis ticks: one per decade spanned by [lo, hi] (the padded plotting
293+ * range), labeled `1e{n}`. If the visible range covers less than a full
294+ * decade — a study that hasn't had time to spread out yet — decade ticks
295+ * would give zero or one of them, so fall back to two ticks at the actual
296+ * (unpadded) data extent instead, labeled with their real value.
297+ */
298+function decadeTicks(
299+ lo: number,
300+ hi: number,
301+ rawLo: number,
302+ rawHi: number,
303+): { value: number; label: string }[] {
304+ const start = Math.ceil(lo);
305+ const end = Math.floor(hi);
306+ const decades: number[] = [];
307+ for (let d = start; d <= end; d++) decades.push(d);
308+ if (decades.length >= 2) {
309+ return decades.map((d) => ({
310+ value: 10 ** d,
311+ label: (10 ** d).toExponential(0).replace('e+', 'e'),
312+ }));
313+ }
314+ const loVal = 10 ** rawLo;
315+ const hiVal = 10 ** rawHi;
316+ if (hiVal <= loVal) return [{ value: loVal, label: fmtValue(loVal) }];
317+ return [
318+ { value: loVal, label: fmtValue(loVal) },
319+ { value: hiVal, label: fmtValue(hiVal) },
320+ ];
321+}
322+
323+/** `count` evenly spaced ticks between t0 and t1 inclusive; just t0 if the
324+ * span is degenerate (a single point pushed so far). */
325+function evenTicks(t0: number, t1: number, count: number): number[] {
326+ if (t1 <= t0) return [t0];
327+ return Array.from({ length: count }, (_, i) => t0 + (i * (t1 - t0)) / (count - 1));
328+}
src/sht/layout.tsmodified+27−0View file
@@ -46,6 +46,33 @@ export function isPowerOfTwo(n: number): boolean {
4646 return n > 0 && (n & (n - 1)) === 0;
4747 }
4848
49+/**
50+ * Re-index coefficients from one (lmax, mmax) layout into another, dropping
51+ * the modes the target cannot hold and zero-filling the ones the source does
52+ * not have. Both truncation and padding are exact spectral operations —
53+ * truncation is orthogonal projection onto the smaller band, padding changes
54+ * nothing — so a surface carried as coefficients can be moved onto any plan.
55+ */
56+export function relayoutCoeffs(
57+ q: Float32Array,
58+ from: { lmax: number; mmax: number },
59+ to: { lmax: number; mmax: number },
60+): Float32Array {
61+ if (from.lmax === to.lmax && from.mmax === to.mmax) return q;
62+ const out = new Float32Array(2 * nlmCalc(to.lmax, to.mmax));
63+ const mTop = Math.min(from.mmax, to.mmax);
64+ const lTop = Math.min(from.lmax, to.lmax);
65+ for (let m = 0; m <= mTop; m++) {
66+ for (let l = m; l <= lTop; l++) {
67+ const src = 2 * lmIndex(from.lmax, l, m);
68+ const dst = 2 * lmIndex(to.lmax, l, m);
69+ out[dst] = q[src];
70+ out[dst + 1] = q[src + 1];
71+ }
72+ }
73+ return out;
74+}
75+
4976 /** Grid sizes for a given lmax, dealiased for a reaction of polynomial degree
5077 * `pdeg` (the rule from websph's reference implementation):
5178 * nlat >= ((pdeg+1)*lmax+1)/2, nphi >= (pdeg+1)*lmax+1. nphi is rounded up to
test/geometryChecks.tsmodified+121−0View file
@@ -282,6 +282,127 @@ export async function geometryChecks(
282282 sht.destroy();
283283 }
284284
285+ // ---- a surface handed over as coefficients is the same surface -----------
286+ // The import path (reharm's "Export to turing-surface"): Geometry.create
287+ // with `coeffs` set evaluates no .m — the coefficients are re-indexed onto
288+ // the session's band and everything from synthesis onward runs unchanged.
289+ // Checked in all three band relations: same (must reproduce the original
290+ // geometry, arrays and all), wider (zero-padding; the surface must still
291+ // match its closed form on the finer plan's own grid), and narrower
292+ // (truncation; exact when the content fits the smaller band, as the
293+ // degree-1 sphere does).
294+ {
295+ const importGeometry = async (
296+ coeffs: { lmax: number; mmax: number; X: Float32Array; Y: Float32Array; Z: Float32Array },
297+ lmax: number,
298+ ) => {
299+ const { nlat, nphi } = gridForLmax(lmax, 3);
300+ const cfg = { lmax, mmax: lmax, nlat, nphi };
301+ const sht = await ShtPlan.create(device, cfg);
302+ const deriv = await DerivPlan.create(device, sht);
303+ const geometry = await Geometry.create({
304+ sht, cfg, source: '', paramNames: [], params: {}, deriv, coeffs,
305+ });
306+ return { sht, deriv, cfg, geometry };
307+ };
308+
309+ const peanut = await buildGeometry(device, 'peanut');
310+ const coeffs = {
311+ lmax: peanut.cfg.lmax,
312+ mmax: peanut.cfg.mmax,
313+ X: peanut.geometry.X,
314+ Y: peanut.geometry.Y,
315+ Z: peanut.geometry.Z,
316+ };
317+
318+ // Same band: the identical pipeline on identical inputs, so the surface
319+ // and both metric formulations must come back numerically unchanged.
320+ const same = await importGeometry(coeffs, peanut.cfg.lmax);
321+ let worst = 0;
322+ const pairs: [Float32Array, Float32Array][] = [
323+ [same.geometry.x, peanut.geometry.x],
324+ [same.geometry.y, peanut.geometry.y],
325+ [same.geometry.z, peanut.geometry.z],
326+ [same.geometry.p1, peanut.geometry.p1],
327+ [same.geometry.p2, peanut.geometry.p2],
328+ [same.geometry.q2, peanut.geometry.q2],
329+ [same.geometry.r, peanut.geometry.r],
330+ [same.geometry.Vtx, peanut.geometry.Vtx],
331+ [same.geometry.Vpx, peanut.geometry.Vpx],
332+ ];
333+ for (const [a, b] of pairs) {
334+ for (let i = 0; i < a.length; i++) {
335+ worst = Math.max(worst, Math.abs(a[i] - b[i]));
336+ }
337+ }
338+ check(
339+ 'import: coefficients rebuild the identical geometry on the same band',
340+ worst <= 1e-7,
341+ `max |difference| = ${worst.toExponential(2)} over surface and metric`,
342+ );
343+ same.deriv.destroy();
344+ same.sht.destroy();
345+
346+ // Wider band: zero-padded onto a finer plan, then checked against the
347+ // .m's own closed form on that plan's grid — points the original build
348+ // never evaluated (Gauss latitudes do not nest).
349+ const p = defaultGeometryParams(peanut.g);
350+ const peanutRadius = (ct: number): number => {
351+ const st2 = Math.max(0, 1 - ct * ct);
352+ const r = 1 - p.waist * st2;
353+ return r * Math.hypot(Math.sqrt(st2), (1 + p.stretch) * ct);
354+ };
355+ const wide = await importGeometry(coeffs, 2 * peanut.cfg.lmax + 1);
356+ let dr = 0;
357+ for (let i = 0; i < wide.cfg.nlat; i++) {
358+ const want = peanutRadius(wide.sht.cosTheta[i]);
359+ for (let j = 0; j < wide.cfg.nphi; j++) {
360+ const k = i * wide.cfg.nphi + j;
361+ const got = Math.hypot(wide.geometry.x[k], wide.geometry.y[k], wide.geometry.z[k]);
362+ dr = Math.max(dr, Math.abs(got - want));
363+ }
364+ }
365+ check(
366+ 'import: zero-padding onto a wider band keeps the surface',
367+ dr < 1e-3,
368+ `max |dr| = ${dr.toExponential(2)} at lmax ${wide.cfg.lmax}`,
369+ );
370+ wide.deriv.destroy();
371+ wide.sht.destroy();
372+ peanut.deriv.destroy();
373+ peanut.sht.destroy();
374+
375+ // Narrower band: truncation is orthogonal projection, so a surface whose
376+ // content fits the smaller band comes through exactly — the sphere is
377+ // degree 1 and must still have radius 1 after import at any lmax.
378+ const sphere = await buildGeometry(device, SPHERE_KEY);
379+ const narrow = await importGeometry(
380+ {
381+ lmax: sphere.cfg.lmax,
382+ mmax: sphere.cfg.mmax,
383+ X: sphere.geometry.X,
384+ Y: sphere.geometry.Y,
385+ Z: sphere.geometry.Z,
386+ },
387+ Math.floor(sphere.cfg.lmax / 2),
388+ );
389+ let radiusErr = 0;
390+ for (let k = 0; k < narrow.geometry.x.length; k++) {
391+ const r = Math.hypot(narrow.geometry.x[k], narrow.geometry.y[k], narrow.geometry.z[k]);
392+ radiusErr = Math.max(radiusErr, Math.abs(r - 1));
393+ }
394+ check(
395+ 'import: truncation onto a narrower band is exact when the content fits',
396+ radiusErr < 2e-3,
397+ `sphere from lmax ${sphere.cfg.lmax} coefficients at lmax ${narrow.cfg.lmax}: ` +
398+ `max |r - 1| = ${radiusErr.toExponential(2)}`,
399+ );
400+ narrow.deriv.destroy();
401+ narrow.sht.destroy();
402+ sphere.deriv.destroy();
403+ sphere.sht.destroy();
404+ }
405+
285406 // ---- the unrolled loop: more ops, identical answer ----------------------
286407 {
287408 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/test-page.tsmodified+2−0View file
@@ -31,6 +31,7 @@ import { geometryChecks } from './geometryChecks.ts';
3131 import { fluxChecks } from './fluxChecks.ts';
3232 import { compareChecks } from './compareChecks.ts';
3333 import { referenceChecks, type H5Rt } from './referenceChecks.ts';
34+import { matlabExportChecks } from './matlabExportChecks.ts';
3435
3536 declare global {
3637 interface Window {
@@ -226,6 +227,7 @@ async function main(): Promise<void> {
226227 await compareChecks(device, check, log);
227228 // '/' is the wasm module's in-memory filesystem — nothing touches disk.
228229 await referenceChecks(h5wasm as unknown as H5Rt, (name) => `/${name}`, check, log);
230+ matlabExportChecks(check, log);
229231
230232 window.__RESULTS__ = { ok: failures === 0, lines };
231233 log(failures === 0 ? 'ALL PASS' : `${failures} FAILURE(S)`);