1/**
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';
28type Check = (name: string, ok: boolean, detail: string) => void;
29type Log = (line: string) => void;
31const COARSE = 31;
32const FINE = 63;
34export async function compareChecks(
35 device: GPUDevice,
36 check: Check,
37 log: Log,
38): Promise<void> {
39 log('\ncompare mode (convergence study):');
41 // ---- prolongCoeffs: every (l, m) lands on itself -------------------------
42 {
43 const src = new Float32Array(2 * nlmCalc(COARSE, COARSE));
44 for (let m = 0; m <= COARSE; m++) {
45 for (let l = m; l <= COARSE; l++) {
46 const i = 2 * lmIndex(COARSE, l, m);
47 src[i] = l + m / 100;
48 src[i + 1] = -l - m / 100;
49 }
50 }
51 const out = prolongCoeffs(src, COARSE, FINE);
52 let moved = 0;
53 let leaked = 0;
54 for (let m = 0; m <= FINE; m++) {
55 for (let l = m; l <= FINE; l++) {
56 const j = 2 * lmIndex(FINE, l, m);
57 if (l <= COARSE && m <= COARSE) {
58 const i = 2 * lmIndex(COARSE, l, m);
59 if (out[j] !== src[i] || out[j + 1] !== src[i + 1]) moved++;
60 } else if (out[j] !== 0 || out[j + 1] !== 0) {
61 leaked++;
62 }
63 }
64 }
65 check(
66 'compare: prolongation puts every coefficient at its own (l, m)',
67 moved === 0 && leaked === 0,
68 `${moved} misplaced, ${leaked} non-zero above the source band ` +
69 `(${nlmCalc(COARSE, COARSE)} -> ${nlmCalc(FINE, FINE)} coefficients)`,
70 );
71 }
73 // ---- one initial condition across lmax, and the control ------------------
74 {
75 const model = mModelByKey('schnakenberg')!;
76 const params = defaultParams(model);
77 const sessions: ModelSession[] = [];
78 try {
79 for (const lmax of [COARSE, FINE]) {
80 sessions.push(await ModelSession.create({ device, model, params, lmax, niter: 1 }));
81 }
82 const [coarse, fine] = sessions;
84 // What the study does: one band-limited field, evaluated on each grid.
85 const noise = await sharedNoise(sessions, model.seedAmp, 1);
86 sessions.forEach((s, i) => s.seedWith(noise[i]));
87 const shared = compareStates(
88 prolongCoeffs(await coarse.read('U'), COARSE, FINE),
89 await fine.read('U'),
90 );
91 check(
92 'compare: one shared perturbation gives both grids the same state',
93 shared.rel < 5e-5,
94 `max |dU| = ${shared.abs.toExponential(2)} ` +
95 `(${(100 * shared.rel).toFixed(4)}% of max |U| = ${shared.scale.toFixed(3)}) ` +
96 `across lmax ${COARSE} vs ${FINE}`,
97 );
99 // The control: the ordinary per-grid seeding these two would otherwise
100 // get. One deviate per grid point, and the grids differ, so the same
101 // integer seed is two different initial conditions -- comparing runs
102 // started this way would report a difference that is entirely the seed.
103 sessions.forEach((s) => s.seed(1));
104 const plain = compareStates(
105 prolongCoeffs(await coarse.read('U'), COARSE, FINE),
106 await fine.read('U'),
107 );
108 check(
109 'compare: control — the same integer seed alone does not do it',
110 plain.abs > 20 * shared.abs,
111 `per-grid seeding differs by ${plain.abs.toExponential(2)}, ` +
112 `${(plain.abs / Math.max(shared.abs, 1e-30)).toExponential(1)}x the shared start's ` +
113 `(seed amplitude ${model.seedAmp})`,
114 );
115 log(
116 ` shared start: ${shared.abs.toExponential(2)}, ` +
117 `per-grid seeds: ${plain.abs.toExponential(2)}`,
118 );
119 } finally {
120 for (const s of sessions) s.destroy();
121 }
122 }
124 // ---- one clock: steps * dt is bit-identical across the divisors ----------
125 {
126 const divisors = [1, 2, 4, 8];
127 const steps = 4;
128 let worst = 0;
129 const cases: string[] = [];
130 for (const model of ['schnakenberg', 'brusselator', 'allencahn']) {
131 const dt = defaultParams(mModelByKey(model)!).dt;
132 for (const div of divisors) {
133 // A variant at dt/div takes div times as many steps to cover the same
134 // span. Powers of two only touch the exponent, so both the divide and
135 // the multiply back are exact and the two spans are the same float.
136 const span = (steps * div) * (dt / div);
137 const ulps = Math.abs(span - steps * dt);
138 if (ulps > worst) worst = ulps;
139 if (div === divisors[divisors.length - 1]) {
140 cases.push(`${model} dt ${dt} -> ${dt / div}`);
141 }
142 }
143 }
144 check(
145 'compare: a power-of-two dt divisor keeps every variant on one clock',
146 worst === 0,
147 `exact for every shipped dt x ${divisors.join('/')} (${cases.join(', ')})`,
148 );
149 }
151 // ---- the variant grid and its reference ---------------------------------
152 {
153 const variants = crossProduct([1, 4], [31, 63], [1, 2]);
154 const ref = variants[mostResolved(variants)];
155 check(
156 'compare: the reference is the most-resolved corner of the grid',
157 variants.length === 8 &&
158 ref.niter === 4 && ref.lmax === 63 && ref.dtDiv === 2 &&
159 new Set(variants.map((v) => `${v.niter}/${v.lmax}/${v.dtDiv}`)).size === 8,
160 `${variants.length} distinct variants, reference niter ${ref.niter} · ` +
161 `lmax ${ref.lmax} · dt/${ref.dtDiv}`,
162 );
163 }
164}
166/** Max absolute difference of two equal-length spectral states, and that
167 * difference relative to the scale of the reference. */
168function compareStates(
169 a: Float32Array,
170 b: Float32Array,
171): { abs: number; rel: number; scale: number } {
172 let abs = 0;
173 let scale = 0;
174 const n = Math.min(a.length, b.length);
175 for (let i = 0; i < n; i++) {
176 abs = Math.max(abs, Math.abs(a[i] - b[i]));
177 scale = Math.max(scale, Math.abs(b[i]));
178 }
179 return { abs, rel: scale > 0 ? abs / scale : Infinity, scale };
180}