2 * The two things a side-by-side comparison of solver settings rests on.
3 *
4 * Both are silent when broken: the panels still animate, the difference norm
5 * still produces a number, and the number is simply wrong — it reports a
6 * disagreement between two runs that were never solving the same problem, or
7 * that were never at the same time. Neither failure looks like a failure, which
8 * is exactly why they are pinned here.
9 *
beac00aMerge main into random-fieldsJeremy Magland 10 * 1. One initial condition, in both of the ways a model can seed. A model
11 * that calls `randnfun3` — every shipped one does — gets a field in space,
12 * so one table drawn once is one field on every grid; the control is the
13 * per-session draw the study must not do. A model that takes `noise` gets
14 * one deviate per grid point, which sharedNoise has to project; the
15 * control there is starker, since the same integer seed on two grids is
16 * simply two unrelated fields.
18 * 2. One clock. dt varies by a power-of-two divisor, so `steps * dt` is
19 * bit-identical across variants and no comparison is ever made across a
20 * fraction of a timestep.
21 *
beac00aMerge main into random-fieldsJeremy Magland 22 * Deliberately small — pairs of sessions at niter 1, lmax 31 and 63 — because a
f467ffaCompare several solver settings side by side, on one clockJeremy Magland 23 * session compiles its whole unrolled step and this suite has to stay short.
24 */
25import { ModelSession } from '../src/mgpu/session.ts';
beac00aMerge main into random-fieldsJeremy Magland 26import { mModelByKey, defaultParams, type MModel, type ParamSpec } from '../src/mgpu/registry.ts';
27import { prolongCoeffs, sharedNoise, sharedModes } from '../src/compare/sharedStart.ts';
28import linearSource from './models/linear.m?raw';
f467ffaCompare several solver settings side by side, on one clockJeremy Magland 29import { lmIndex, nlmCalc } from '../src/sht/layout.ts';
30import { crossProduct, mostResolved } from '../src/compare/variants.ts';
ef3ae33Do not stretch the colormap across a constant field's roundoffJeremy Magland 31import { floorRange } from '../src/render/colorbar.ts';
33type Check = (name: string, ok: boolean, detail: string) => void;
34type Log = (line: string) => void;
36const COARSE = 31;
37const FINE = 63;
39export async function compareChecks(
40 device: GPUDevice,
41 check: Check,
42 log: Log,
43): Promise<void> {
44 log('\ncompare mode (convergence study):');
46 // ---- prolongCoeffs: every (l, m) lands on itself -------------------------
47 {
48 const src = new Float32Array(2 * nlmCalc(COARSE, COARSE));
49 for (let m = 0; m <= COARSE; m++) {
50 for (let l = m; l <= COARSE; l++) {
51 const i = 2 * lmIndex(COARSE, l, m);
52 src[i] = l + m / 100;
53 src[i + 1] = -l - m / 100;
54 }
55 }
56 const out = prolongCoeffs(src, COARSE, FINE);
57 let moved = 0;
58 let leaked = 0;
59 for (let m = 0; m <= FINE; m++) {
60 for (let l = m; l <= FINE; l++) {
61 const j = 2 * lmIndex(FINE, l, m);
62 if (l <= COARSE && m <= COARSE) {
63 const i = 2 * lmIndex(COARSE, l, m);
64 if (out[j] !== src[i] || out[j + 1] !== src[i + 1]) moved++;
65 } else if (out[j] !== 0 || out[j + 1] !== 0) {
66 leaked++;
67 }
68 }
69 }
70 check(
71 'compare: prolongation puts every coefficient at its own (l, m)',
72 moved === 0 && leaked === 0,
73 `${moved} misplaced, ${leaked} non-zero above the source band ` +
74 `(${nlmCalc(COARSE, COARSE)} -> ${nlmCalc(FINE, FINE)} coefficients)`,
75 );
76 }
c90d0e2Check reference files in the browser's compare modeJeremy Magland 78 // ---- a file's exact state loads onto every grid --------------------------
79 // What a reference-file study does instead of seeding: the file's spectral
80 // state pushed into each variant by loadState, prolonged into its band. The
81 // load is a plain upload, so the state must come back bit-exact; and read on
82 // one shared grid the variants must then show one field, because synthesis
83 // of the same band-limited coefficients is evaluation, not resampling.
84 {
85 const model = mModelByKey('allencahn')!;
86 const params = defaultParams(model);
87 const sessions: ModelSession[] = [];
88 try {
89 for (const lmax of [COARSE, FINE]) {
90 sessions.push(await ModelSession.create({ device, model, params, lmax, niter: 0 }));
91 }
92 const [coarse, fine] = sessions;
93 // A deterministic band-limited state, decaying like a real spectrum;
94 // m = 0 imaginary parts stay zero (the state is a real field).
95 const q = new Float32Array(2 * nlmCalc(COARSE, COARSE));
96 for (let m = 0; m <= COARSE; m++) {
97 for (let l = m; l <= COARSE; l++) {
98 const i = 2 * lmIndex(COARSE, l, m);
99 const amp = Math.exp(-l / 6);
100 q[i] = amp * Math.sin(1 + 3 * l + 7 * m);
101 q[i + 1] = m === 0 ? 0 : amp * Math.cos(2 + 5 * l + 11 * m);
102 }
103 }
104 coarse.loadState({ U: q });
105 fine.loadState({ U: prolongCoeffs(q, COARSE, FINE) });
107 const back = await coarse.read('U');
108 let exact = back.length === q.length;
109 if (exact) {
110 for (let i = 0; i < q.length; i++) {
111 if (back[i] !== q[i]) {
112 exact = false;
113 break;
114 }
115 }
116 }
117 check(
118 'compare: loadState puts the exact coefficients in the state',
119 exact,
120 `${q.length} float32 values round-tripped bit-exact at lmax ${COARSE}`,
121 );
123 // The coarse session's own solver grid, so its display plan is the
124 // solver's — the branch a crowded study lands on.
125 for (const s of sessions) await s.setDisplayGrid(64, 128);
126 const cu = await coarse.readSpecies(0);
127 const fu = await fine.readSpecies(0);
128 let maxd = 0;
129 let scale = 0;
130 for (let i = 0; i < cu.length; i++) {
131 maxd = Math.max(maxd, Math.abs(cu[i] - fu[i]));
132 scale = Math.max(scale, Math.abs(cu[i]));
133 }
134 check(
135 'compare: one loaded state reads back as one field on a shared grid',
136 maxd < 1e-4 * scale,
137 `max |du| = ${maxd.toExponential(2)} vs max |u| = ${scale.toExponential(2)} ` +
138 `across lmax ${COARSE} vs ${FINE}`,
139 );
140 } finally {
141 for (const s of sessions) s.destroy();
142 }
143 }
beac00aMerge main into random-fieldsJeremy Magland 145 // ---- one random field across lmax: the shipped models' seeding -----------
147 const model = mModelByKey('schnakenberg')!;
148 const params = defaultParams(model);
149 const sessions: ModelSession[] = [];
150 try {
151 for (const lmax of [COARSE, FINE]) {
152 sessions.push(await ModelSession.create({ device, model, params, lmax, niter: 1 }));
153 }
154 const [coarse, fine] = sessions;
beac00aMerge main into random-fieldsJeremy Magland 156 // What the study does: one coefficient table, drawn once, summed on each
157 // variant's own grid points. The residual is the coarse grid's analysis of
158 // a field with a little content above its band, not a difference in the
159 // field -- so it is bounded by the perturbation, not by |U|, which is why
160 // compareStates measures against the non-constant part.
f467ffaCompare several solver settings side by side, on one clockJeremy Magland 161 const noise = await sharedNoise(sessions, model.seedAmp, 1);
163 check(
164 'compare: a randnfun3 model seeds every variant from one drawn table',
165 modes !== null,
166 modes ? `${modes[0]} Fourier modes, one table for both grids` : 'no table drawn',
167 );
168 for (let i = 0; i < sessions.length; i++) await sessions[i].seedWith(noise[i], modes);
169 const shared = compareStates(
170 prolongCoeffs(await coarse.read('U'), COARSE, FINE),
171 await fine.read('U'),
172 );
173 check(
174 'compare: one shared random field gives both grids the same state',
175 shared.rel < 1e-3,
176 `max |dU| = ${shared.abs.toExponential(2)} ` +
177 `(${(100 * shared.rel).toFixed(3)}% of the perturbation, max ` +
178 `${shared.scale.toExponential(2)}) across lmax ${COARSE} vs ${FINE}`,
179 );
181 // The control, and the reason `sharedModes` exists: left to seed itself
182 // each session draws over *its own* bounding box, and a box is grid
183 // samples of the surface, so the two draws are near neighbours rather than
184 // one field. A tolerance is not what separates them — the shared table is
185 // simply closer, and would be however either number moved.
186 await coarse.seed(1);
187 await fine.seed(1);
188 const own = compareStates(
189 prolongCoeffs(await coarse.read('U'), COARSE, FINE),
190 await fine.read('U'),
191 );
192 check(
193 'compare: control — a per-session draw is not the same field',
194 own.abs > shared.abs,
195 `per-session draws differ by ${own.abs.toExponential(2)}, ` +
196 `${(own.abs / Math.max(shared.abs, 1e-30)).toFixed(1)}x the shared table's ` +
197 `${shared.abs.toExponential(2)}`,
198 );
199 } finally {
200 for (const s of sessions) s.destroy();
201 }
202 }
204 // ---- one grid-point perturbation across lmax, and the control ------------
205 // The other way a model can seed: `init(noise)` takes the host's field
206 // directly, one deviate per grid point (the test models here, and any .m
207 // edited to do it). Nothing about it is a function of space, so this is the
208 // case sharedNoise's projection is for — and the case where the same integer
209 // seed on two grids gives two entirely unrelated initial conditions.
210 {
211 const params = { c: 0, D: 1e-3, dt: 0.05 };
212 const model = noiseModel();
213 const sessions: ModelSession[] = [];
214 try {
215 for (const lmax of [COARSE, FINE]) {
216 sessions.push(await ModelSession.create({ device, model, params, lmax, niter: 1 }));
217 }
218 const [coarse, fine] = sessions;
220 const noise = await sharedNoise(sessions, model.seedAmp, 1);
221 for (let i = 0; i < sessions.length; i++) await sessions[i].seedWith(noise[i]);
f467ffaCompare several solver settings side by side, on one clockJeremy Magland 222 const shared = compareStates(
223 prolongCoeffs(await coarse.read('U'), COARSE, FINE),
224 await fine.read('U'),
225 );
226 check(
227 'compare: one shared perturbation gives both grids the same state',
228 shared.rel < 5e-5,
229 `max |dU| = ${shared.abs.toExponential(2)} ` +
beac00aMerge main into random-fieldsJeremy Magland 230 `(${(100 * shared.rel).toFixed(4)}% of max |U| = ${shared.scale.toExponential(2)}) ` +
f467ffaCompare several solver settings side by side, on one clockJeremy Magland 231 `across lmax ${COARSE} vs ${FINE}`,
232 );
235 await fine.seed(1);
f467ffaCompare several solver settings side by side, on one clockJeremy Magland 236 const plain = compareStates(
237 prolongCoeffs(await coarse.read('U'), COARSE, FINE),
238 await fine.read('U'),
239 );
240 check(
241 'compare: control — the same integer seed alone does not do it',
242 plain.abs > 20 * shared.abs,
243 `per-grid seeding differs by ${plain.abs.toExponential(2)}, ` +
244 `${(plain.abs / Math.max(shared.abs, 1e-30)).toExponential(1)}x the shared start's ` +
245 `(seed amplitude ${model.seedAmp})`,
246 );
247 log(
248 ` shared start: ${shared.abs.toExponential(2)}, ` +
249 `per-grid seeds: ${plain.abs.toExponential(2)}`,
250 );
251 } finally {
252 for (const s of sessions) s.destroy();
253 }
254 }
256 // ---- one clock: steps * dt is bit-identical across the divisors ----------
257 {
258 const divisors = [1, 2, 4, 8];
259 const steps = 4;
260 let worst = 0;
261 const cases: string[] = [];
262 for (const model of ['schnakenberg', 'brusselator', 'allencahn']) {
263 const dt = defaultParams(mModelByKey(model)!).dt;
264 for (const div of divisors) {
265 // A variant at dt/div takes div times as many steps to cover the same
266 // span. Powers of two only touch the exponent, so both the divide and
267 // the multiply back are exact and the two spans are the same float.
268 const span = (steps * div) * (dt / div);
269 const ulps = Math.abs(span - steps * dt);
270 if (ulps > worst) worst = ulps;
271 if (div === divisors[divisors.length - 1]) {
272 cases.push(`${model} dt ${dt} -> ${dt / div}`);
273 }
274 }
275 }
276 check(
277 'compare: a power-of-two dt divisor keeps every variant on one clock',
278 worst === 0,
279 `exact for every shipped dt x ${divisors.join('/')} (${cases.join(', ')})`,
280 );
281 }
ef3ae33Do not stretch the colormap across a constant field's roundoffJeremy Magland 283 // ---- a uniform field is drawn uniform, on every grid --------------------
284 // Schnakenberg seeds v as a literal constant (`vs * ones(...)`), so its whole
285 // spread is the fp32 residue of the analys/synth round trip -- pole-localized
286 // and grid-dependent, so scaled to its own extremes it paints two unrelated
287 // pictures of the same constant, which is what a broken seeding would look
288 // like. The spans below are measured (worst |deviation| x 2, on vs = 0.9):
289 // lmax 63, 127, 255. See floorRange for where they come from.
290 {
291 const vs = 0.9;
292 const spans = [5.2e-5, 1.8e-4, 4.4e-4];
293 // Each must end up a small slice of the drawn range rather than all of it.
294 const shares = spans.map((sp) => sp / (floorRange(vs - sp / 2, vs + sp / 2).hi -
295 floorRange(vs - sp / 2, vs + sp / 2).lo));
296 // ...while real structure keeps its own range exactly. v once the spots
297 // have formed spans ~0.03 on the same 0.9, two orders above the residue.
298 const real = floorRange(0.895, 0.924);
299 check(
300 'compare: fp32 residue on a constant field does not become a picture',
301 shares.every((s) => s < 0.1) && real.lo === 0.895 && real.hi === 0.924,
302 `residue uses ${shares.map((s) => `${(100 * s).toFixed(1)}%`).join(', ')} ` +
303 `of the colormap at lmax 63/127/255; real pattern ` +
304 `[${real.lo}, ${real.hi}] left untouched`,
305 );
306 }
f467ffaCompare several solver settings side by side, on one clockJeremy Magland 308 // ---- the variant grid and its reference ---------------------------------
309 {
310 const variants = crossProduct([1, 4], [31, 63], [1, 2]);
311 const ref = variants[mostResolved(variants)];
312 check(
313 'compare: the reference is the most-resolved corner of the grid',
314 variants.length === 8 &&
315 ref.niter === 4 && ref.lmax === 63 && ref.dtDiv === 2 &&
316 new Set(variants.map((v) => `${v.niter}/${v.lmax}/${v.dtDiv}`)).size === 8,
317 `${variants.length} distinct variants, reference niter ${ref.niter} · ` +
318 `lmax ${ref.lmax} · dt/${ref.dtDiv}`,
319 );
320 }
321}
324 * The one-species linear test model, seeded from `noise` rather than from a
325 * random field — `init(noise)`, so the host's grid-point field is what reaches
326 * the state (test/models/linear.m). Never stepped here; the parameters exist
327 * because the .m names them.
328 */
329function noiseModel(): MModel {
330 const param = (key: string): ParamSpec => ({
331 key, label: key, value: 0, min: -1e9, max: 1e9, step: 1,
332 });
333 return {
334 key: 'linear',
335 label: 'linear',
336 blurb: '',
337 species: ['u'],
338 state: ['U'],
339 params: ['c', 'D', 'dt'].map(param),
340 pdeg: 1,
341 seedAmp: 1e-2,
342 source: linearSource,
343 };
344}
346/**
347 * Max absolute difference of two equal-length spectral states, and that
348 * difference relative to the scale of the reference's *non-constant* part —
349 * every coefficient but (l, m) = (0, 0), which is index 0 in either layout.
350 *
351 * Normalizing against the whole state would hide the question. A model seeded as
352 * a perturbation of a uniform steady state puts that state in (0, 0) alone, two
353 * orders above everything else, so |dU| / max |U| would report a comfortable
354 * fraction of the *background* however unrelated the two perturbations were —
355 * including when no perturbation arrived at all, which is what a table that
356 * never reaches a session looks like. Against the perturbation, that failure
357 * reads as a ratio of 1.
358 */
f467ffaCompare several solver settings side by side, on one clockJeremy Magland 359function compareStates(
360 a: Float32Array,
361 b: Float32Array,
362): { abs: number; rel: number; scale: number } {
363 let abs = 0;
364 let scale = 0;
365 const n = Math.min(a.length, b.length);
366 for (let i = 0; i < n; i++) {
367 abs = Math.max(abs, Math.abs(a[i] - b[i]));
beac00aMerge main into random-fieldsJeremy Magland 368 if (i >= 2) scale = Math.max(scale, Math.abs(b[i]));
370 return { abs, rel: scale > 0 ? abs / scale : Infinity, scale };
371}