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 }
beac00aMerge main into random-fieldsJeremy Magland 78 // ---- one random field across lmax: the shipped models' seeding -----------
80 const model = mModelByKey('schnakenberg')!;
81 const params = defaultParams(model);
82 const sessions: ModelSession[] = [];
83 try {
84 for (const lmax of [COARSE, FINE]) {
85 sessions.push(await ModelSession.create({ device, model, params, lmax, niter: 1 }));
86 }
87 const [coarse, fine] = sessions;
beac00aMerge main into random-fieldsJeremy Magland 89 // What the study does: one coefficient table, drawn once, summed on each
90 // variant's own grid points. The residual is the coarse grid's analysis of
91 // a field with a little content above its band, not a difference in the
92 // field -- so it is bounded by the perturbation, not by |U|, which is why
93 // compareStates measures against the non-constant part.
f467ffaCompare several solver settings side by side, on one clockJeremy Magland 94 const noise = await sharedNoise(sessions, model.seedAmp, 1);
96 check(
97 'compare: a randnfun3 model seeds every variant from one drawn table',
98 modes !== null,
99 modes ? `${modes[0]} Fourier modes, one table for both grids` : 'no table drawn',
100 );
101 for (let i = 0; i < sessions.length; i++) await sessions[i].seedWith(noise[i], modes);
102 const shared = compareStates(
103 prolongCoeffs(await coarse.read('U'), COARSE, FINE),
104 await fine.read('U'),
105 );
106 check(
107 'compare: one shared random field gives both grids the same state',
108 shared.rel < 1e-3,
109 `max |dU| = ${shared.abs.toExponential(2)} ` +
110 `(${(100 * shared.rel).toFixed(3)}% of the perturbation, max ` +
111 `${shared.scale.toExponential(2)}) across lmax ${COARSE} vs ${FINE}`,
112 );
114 // The control, and the reason `sharedModes` exists: left to seed itself
115 // each session draws over *its own* bounding box, and a box is grid
116 // samples of the surface, so the two draws are near neighbours rather than
117 // one field. A tolerance is not what separates them — the shared table is
118 // simply closer, and would be however either number moved.
119 await coarse.seed(1);
120 await fine.seed(1);
121 const own = compareStates(
122 prolongCoeffs(await coarse.read('U'), COARSE, FINE),
123 await fine.read('U'),
124 );
125 check(
126 'compare: control — a per-session draw is not the same field',
127 own.abs > shared.abs,
128 `per-session draws differ by ${own.abs.toExponential(2)}, ` +
129 `${(own.abs / Math.max(shared.abs, 1e-30)).toFixed(1)}x the shared table's ` +
130 `${shared.abs.toExponential(2)}`,
131 );
132 } finally {
133 for (const s of sessions) s.destroy();
134 }
135 }
137 // ---- one grid-point perturbation across lmax, and the control ------------
138 // The other way a model can seed: `init(noise)` takes the host's field
139 // directly, one deviate per grid point (the test models here, and any .m
140 // edited to do it). Nothing about it is a function of space, so this is the
141 // case sharedNoise's projection is for — and the case where the same integer
142 // seed on two grids gives two entirely unrelated initial conditions.
143 {
144 const params = { c: 0, D: 1e-3, dt: 0.05 };
145 const model = noiseModel();
146 const sessions: ModelSession[] = [];
147 try {
148 for (const lmax of [COARSE, FINE]) {
149 sessions.push(await ModelSession.create({ device, model, params, lmax, niter: 1 }));
150 }
151 const [coarse, fine] = sessions;
153 const noise = await sharedNoise(sessions, model.seedAmp, 1);
154 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 155 const shared = compareStates(
156 prolongCoeffs(await coarse.read('U'), COARSE, FINE),
157 await fine.read('U'),
158 );
159 check(
160 'compare: one shared perturbation gives both grids the same state',
161 shared.rel < 5e-5,
162 `max |dU| = ${shared.abs.toExponential(2)} ` +
beac00aMerge main into random-fieldsJeremy Magland 163 `(${(100 * shared.rel).toFixed(4)}% of max |U| = ${shared.scale.toExponential(2)}) ` +
f467ffaCompare several solver settings side by side, on one clockJeremy Magland 164 `across lmax ${COARSE} vs ${FINE}`,
165 );
168 await fine.seed(1);
f467ffaCompare several solver settings side by side, on one clockJeremy Magland 169 const plain = compareStates(
170 prolongCoeffs(await coarse.read('U'), COARSE, FINE),
171 await fine.read('U'),
172 );
173 check(
174 'compare: control — the same integer seed alone does not do it',
175 plain.abs > 20 * shared.abs,
176 `per-grid seeding differs by ${plain.abs.toExponential(2)}, ` +
177 `${(plain.abs / Math.max(shared.abs, 1e-30)).toExponential(1)}x the shared start's ` +
178 `(seed amplitude ${model.seedAmp})`,
179 );
180 log(
181 ` shared start: ${shared.abs.toExponential(2)}, ` +
182 `per-grid seeds: ${plain.abs.toExponential(2)}`,
183 );
184 } finally {
185 for (const s of sessions) s.destroy();
186 }
187 }
189 // ---- one clock: steps * dt is bit-identical across the divisors ----------
190 {
191 const divisors = [1, 2, 4, 8];
192 const steps = 4;
193 let worst = 0;
194 const cases: string[] = [];
195 for (const model of ['schnakenberg', 'brusselator', 'allencahn']) {
196 const dt = defaultParams(mModelByKey(model)!).dt;
197 for (const div of divisors) {
198 // A variant at dt/div takes div times as many steps to cover the same
199 // span. Powers of two only touch the exponent, so both the divide and
200 // the multiply back are exact and the two spans are the same float.
201 const span = (steps * div) * (dt / div);
202 const ulps = Math.abs(span - steps * dt);
203 if (ulps > worst) worst = ulps;
204 if (div === divisors[divisors.length - 1]) {
205 cases.push(`${model} dt ${dt} -> ${dt / div}`);
206 }
207 }
208 }
209 check(
210 'compare: a power-of-two dt divisor keeps every variant on one clock',
211 worst === 0,
212 `exact for every shipped dt x ${divisors.join('/')} (${cases.join(', ')})`,
213 );
214 }
ef3ae33Do not stretch the colormap across a constant field's roundoffJeremy Magland 216 // ---- a uniform field is drawn uniform, on every grid --------------------
217 // Schnakenberg seeds v as a literal constant (`vs * ones(...)`), so its whole
218 // spread is the fp32 residue of the analys/synth round trip -- pole-localized
219 // and grid-dependent, so scaled to its own extremes it paints two unrelated
220 // pictures of the same constant, which is what a broken seeding would look
221 // like. The spans below are measured (worst |deviation| x 2, on vs = 0.9):
222 // lmax 63, 127, 255. See floorRange for where they come from.
223 {
224 const vs = 0.9;
225 const spans = [5.2e-5, 1.8e-4, 4.4e-4];
226 // Each must end up a small slice of the drawn range rather than all of it.
227 const shares = spans.map((sp) => sp / (floorRange(vs - sp / 2, vs + sp / 2).hi -
228 floorRange(vs - sp / 2, vs + sp / 2).lo));
229 // ...while real structure keeps its own range exactly. v once the spots
230 // have formed spans ~0.03 on the same 0.9, two orders above the residue.
231 const real = floorRange(0.895, 0.924);
232 check(
233 'compare: fp32 residue on a constant field does not become a picture',
234 shares.every((s) => s < 0.1) && real.lo === 0.895 && real.hi === 0.924,
235 `residue uses ${shares.map((s) => `${(100 * s).toFixed(1)}%`).join(', ')} ` +
236 `of the colormap at lmax 63/127/255; real pattern ` +
237 `[${real.lo}, ${real.hi}] left untouched`,
238 );
239 }
f467ffaCompare several solver settings side by side, on one clockJeremy Magland 241 // ---- the variant grid and its reference ---------------------------------
242 {
243 const variants = crossProduct([1, 4], [31, 63], [1, 2]);
244 const ref = variants[mostResolved(variants)];
245 check(
246 'compare: the reference is the most-resolved corner of the grid',
247 variants.length === 8 &&
248 ref.niter === 4 && ref.lmax === 63 && ref.dtDiv === 2 &&
249 new Set(variants.map((v) => `${v.niter}/${v.lmax}/${v.dtDiv}`)).size === 8,
250 `${variants.length} distinct variants, reference niter ${ref.niter} · ` +
251 `lmax ${ref.lmax} · dt/${ref.dtDiv}`,
252 );
253 }
254}
257 * The one-species linear test model, seeded from `noise` rather than from a
258 * random field — `init(noise)`, so the host's grid-point field is what reaches
259 * the state (test/models/linear.m). Never stepped here; the parameters exist
260 * because the .m names them.
261 */
262function noiseModel(): MModel {
263 const param = (key: string): ParamSpec => ({
264 key, label: key, value: 0, min: -1e9, max: 1e9, step: 1,
265 });
266 return {
267 key: 'linear',
268 label: 'linear',
269 blurb: '',
270 species: ['u'],
271 state: ['U'],
272 params: ['c', 'D', 'dt'].map(param),
273 pdeg: 1,
274 seedAmp: 1e-2,
275 source: linearSource,
276 };
277}
279/**
280 * Max absolute difference of two equal-length spectral states, and that
281 * difference relative to the scale of the reference's *non-constant* part —
282 * every coefficient but (l, m) = (0, 0), which is index 0 in either layout.
283 *
284 * Normalizing against the whole state would hide the question. A model seeded as
285 * a perturbation of a uniform steady state puts that state in (0, 0) alone, two
286 * orders above everything else, so |dU| / max |U| would report a comfortable
287 * fraction of the *background* however unrelated the two perturbations were —
288 * including when no perturbation arrived at all, which is what a table that
289 * never reaches a session looks like. Against the perturbation, that failure
290 * reads as a ratio of 1.
291 */
f467ffaCompare several solver settings side by side, on one clockJeremy Magland 292function compareStates(
293 a: Float32Array,
294 b: Float32Array,
295): { abs: number; rel: number; scale: number } {
296 let abs = 0;
297 let scale = 0;
298 const n = Math.min(a.length, b.length);
299 for (let i = 0; i < n; i++) {
300 abs = Math.max(abs, Math.abs(a[i] - b[i]));
beac00aMerge main into random-fieldsJeremy Magland 301 if (i >= 2) scale = Math.max(scale, Math.abs(b[i]));
303 return { abs, rel: scale > 0 ? abs / scale : Infinity, scale };
304}