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 *
10 * 1. One initial condition. Sessions at different lmax seeded through
11 * sharedNoise hold the *same* spectral state, zero-padded — and the
12 * control shows what that is owed to: seeded the ordinary per-grid way,
13 * the same integer seed gives two unrelated fields.
14 *
15 * 2. One clock. dt varies by a power-of-two divisor, so `steps * dt` is
16 * bit-identical across variants and no comparison is ever made across a
17 * fraction of a timestep.
18 *
19 * Deliberately small — two sessions at niter 1, lmax 31 and 63 — because a
20 * session compiles its whole unrolled step and this suite has to stay short.
21 */
22import { ModelSession } from '../src/mgpu/session.ts';
23import { mModelByKey, defaultParams } from '../src/mgpu/registry.ts';
24import { prolongCoeffs, sharedNoise } from '../src/compare/sharedStart.ts';
25import { lmIndex, nlmCalc } from '../src/sht/layout.ts';
26import { crossProduct, mostResolved } from '../src/compare/variants.ts';
ef3ae33Do not stretch the colormap across a constant field's roundoffJeremy Magland 27import { floorRange } from '../src/render/colorbar.ts';
29type Check = (name: string, ok: boolean, detail: string) => void;
30type Log = (line: string) => void;
32const COARSE = 31;
33const FINE = 63;
35export async function compareChecks(
36 device: GPUDevice,
37 check: Check,
38 log: Log,
39): Promise<void> {
40 log('\ncompare mode (convergence study):');
42 // ---- prolongCoeffs: every (l, m) lands on itself -------------------------
43 {
44 const src = new Float32Array(2 * nlmCalc(COARSE, COARSE));
45 for (let m = 0; m <= COARSE; m++) {
46 for (let l = m; l <= COARSE; l++) {
47 const i = 2 * lmIndex(COARSE, l, m);
48 src[i] = l + m / 100;
49 src[i + 1] = -l - m / 100;
50 }
51 }
52 const out = prolongCoeffs(src, COARSE, FINE);
53 let moved = 0;
54 let leaked = 0;
55 for (let m = 0; m <= FINE; m++) {
56 for (let l = m; l <= FINE; l++) {
57 const j = 2 * lmIndex(FINE, l, m);
58 if (l <= COARSE && m <= COARSE) {
59 const i = 2 * lmIndex(COARSE, l, m);
60 if (out[j] !== src[i] || out[j + 1] !== src[i + 1]) moved++;
61 } else if (out[j] !== 0 || out[j + 1] !== 0) {
62 leaked++;
63 }
64 }
65 }
66 check(
67 'compare: prolongation puts every coefficient at its own (l, m)',
68 moved === 0 && leaked === 0,
69 `${moved} misplaced, ${leaked} non-zero above the source band ` +
70 `(${nlmCalc(COARSE, COARSE)} -> ${nlmCalc(FINE, FINE)} coefficients)`,
71 );
72 }
74 // ---- one initial condition across lmax, and the control ------------------
75 {
76 const model = mModelByKey('schnakenberg')!;
77 const params = defaultParams(model);
78 const sessions: ModelSession[] = [];
79 try {
80 for (const lmax of [COARSE, FINE]) {
81 sessions.push(await ModelSession.create({ device, model, params, lmax, niter: 1 }));
82 }
83 const [coarse, fine] = sessions;
85 // What the study does: one band-limited field, evaluated on each grid.
86 const noise = await sharedNoise(sessions, model.seedAmp, 1);
87 sessions.forEach((s, i) => s.seedWith(noise[i]));
88 const shared = compareStates(
89 prolongCoeffs(await coarse.read('U'), COARSE, FINE),
90 await fine.read('U'),
91 );
92 check(
93 'compare: one shared perturbation gives both grids the same state',
94 shared.rel < 5e-5,
95 `max |dU| = ${shared.abs.toExponential(2)} ` +
96 `(${(100 * shared.rel).toFixed(4)}% of max |U| = ${shared.scale.toFixed(3)}) ` +
97 `across lmax ${COARSE} vs ${FINE}`,
98 );
100 // The control: the ordinary per-grid seeding these two would otherwise
101 // get. One deviate per grid point, and the grids differ, so the same
102 // integer seed is two different initial conditions -- comparing runs
103 // started this way would report a difference that is entirely the seed.
104 sessions.forEach((s) => s.seed(1));
105 const plain = compareStates(
106 prolongCoeffs(await coarse.read('U'), COARSE, FINE),
107 await fine.read('U'),
108 );
109 check(
110 'compare: control — the same integer seed alone does not do it',
111 plain.abs > 20 * shared.abs,
112 `per-grid seeding differs by ${plain.abs.toExponential(2)}, ` +
113 `${(plain.abs / Math.max(shared.abs, 1e-30)).toExponential(1)}x the shared start's ` +
114 `(seed amplitude ${model.seedAmp})`,
115 );
116 log(
117 ` shared start: ${shared.abs.toExponential(2)}, ` +
118 `per-grid seeds: ${plain.abs.toExponential(2)}`,
119 );
120 } finally {
121 for (const s of sessions) s.destroy();
122 }
123 }
125 // ---- one clock: steps * dt is bit-identical across the divisors ----------
126 {
127 const divisors = [1, 2, 4, 8];
128 const steps = 4;
129 let worst = 0;
130 const cases: string[] = [];
131 for (const model of ['schnakenberg', 'brusselator', 'allencahn']) {
132 const dt = defaultParams(mModelByKey(model)!).dt;
133 for (const div of divisors) {
134 // A variant at dt/div takes div times as many steps to cover the same
135 // span. Powers of two only touch the exponent, so both the divide and
136 // the multiply back are exact and the two spans are the same float.
137 const span = (steps * div) * (dt / div);
138 const ulps = Math.abs(span - steps * dt);
139 if (ulps > worst) worst = ulps;
140 if (div === divisors[divisors.length - 1]) {
141 cases.push(`${model} dt ${dt} -> ${dt / div}`);
142 }
143 }
144 }
145 check(
146 'compare: a power-of-two dt divisor keeps every variant on one clock',
147 worst === 0,
148 `exact for every shipped dt x ${divisors.join('/')} (${cases.join(', ')})`,
149 );
150 }
ef3ae33Do not stretch the colormap across a constant field's roundoffJeremy Magland 152 // ---- a uniform field is drawn uniform, on every grid --------------------
153 // Schnakenberg seeds v as a literal constant (`vs * ones(...)`), so its whole
154 // spread is the fp32 residue of the analys/synth round trip -- pole-localized
155 // and grid-dependent, so scaled to its own extremes it paints two unrelated
156 // pictures of the same constant, which is what a broken seeding would look
157 // like. The spans below are measured (worst |deviation| x 2, on vs = 0.9):
158 // lmax 63, 127, 255. See floorRange for where they come from.
159 {
160 const vs = 0.9;
161 const spans = [5.2e-5, 1.8e-4, 4.4e-4];
162 // Each must end up a small slice of the drawn range rather than all of it.
163 const shares = spans.map((sp) => sp / (floorRange(vs - sp / 2, vs + sp / 2).hi -
164 floorRange(vs - sp / 2, vs + sp / 2).lo));
165 // ...while real structure keeps its own range exactly. v once the spots
166 // have formed spans ~0.03 on the same 0.9, two orders above the residue.
167 const real = floorRange(0.895, 0.924);
168 check(
169 'compare: fp32 residue on a constant field does not become a picture',
170 shares.every((s) => s < 0.1) && real.lo === 0.895 && real.hi === 0.924,
171 `residue uses ${shares.map((s) => `${(100 * s).toFixed(1)}%`).join(', ')} ` +
172 `of the colormap at lmax 63/127/255; real pattern ` +
173 `[${real.lo}, ${real.hi}] left untouched`,
174 );
175 }
f467ffaCompare several solver settings side by side, on one clockJeremy Magland 177 // ---- the variant grid and its reference ---------------------------------
178 {
179 const variants = crossProduct([1, 4], [31, 63], [1, 2]);
180 const ref = variants[mostResolved(variants)];
181 check(
182 'compare: the reference is the most-resolved corner of the grid',
183 variants.length === 8 &&
184 ref.niter === 4 && ref.lmax === 63 && ref.dtDiv === 2 &&
185 new Set(variants.map((v) => `${v.niter}/${v.lmax}/${v.dtDiv}`)).size === 8,
186 `${variants.length} distinct variants, reference niter ${ref.niter} · ` +
187 `lmax ${ref.lmax} · dt/${ref.dtDiv}`,
188 );
189 }
190}
192/** Max absolute difference of two equal-length spectral states, and that
193 * difference relative to the scale of the reference. */
194function compareStates(
195 a: Float32Array,
196 b: Float32Array,
197): { abs: number; rel: number; scale: number } {
198 let abs = 0;
199 let scale = 0;
200 const n = Math.min(a.length, b.length);
201 for (let i = 0; i < n; i++) {
202 abs = Math.max(abs, Math.abs(a[i] - b[i]));
203 scale = Math.max(scale, Math.abs(b[i]));
204 }
205 return { abs, rel: scale > 0 ? abs / scale : Infinity, scale };
206}