1/**
2 * Import a reference HDF5 file — geometry, initial and final spherical-
3 * harmonic coefficients for a run of this repo's solver, in the format
4 * documented alongside the sibling test-data repo's case files (see
5 * ../turing-surface-test-data/cases/) — run this repo's own solver from that
6 * file's exact initial condition, and report the numerical error against
7 * its final state.
8 *
9 * This is the regression check for the surface Laplace-Beltrami correction:
10 * replay a saved-off run and see how far this repo's own output has drifted
11 * (or use --niter to probe how much the correction term itself matters).
12 *
13 * npm run ref -- --in data/schnak-spots.h5
14 * npm run ref -- --in data/schnak-spots.h5 --niter 0
15 */
16import { requestShtDevice, describeAdapter } from '../src/sht/sht.ts';
17import { ModelSession } from '../src/mgpu/session.ts';
18import { mModelByKey, defaultParams, type Params } from '../src/mgpu/registry.ts';
19import { mGeometryByKey, defaultGeometryParams } from '../src/geom/registry.ts';
20import { relL2, relLinf } from '../src/mgpu/digest.ts';
21import { installWebGpu, errMsg, NO_ADAPTER_HINT } from './nodeWebGpu.ts';
22import * as h5wasm from 'h5wasm/node';
24const USAGE = `usage: npm run ref -- --in <file> [options]
26 --in <file> the reference HDF5 file to check against (required)
27 --niter <n> override the solve iteration count (default: the file's own)
28 --tolerance <n> if given, exit 1 when any reported relL2 meets or exceeds it
29 --tolerance-linf <n> if given, exit 1 when any reported relLinf meets or exceeds it
30 --json machine-readable output
31 --help
33Runs this repo's solver from the file's exact initial spectral state, to the
34same physical end time, and reports the relative-L2 and relative-L-infinity
35(max-norm) error of the resulting state against the file's final state (and,
36as a sanity check, of the regenerated geometry against the file's own
37geometry coefficients). --tolerance and --tolerance-linf gate independently:
38either can fail the run on its own.`;
40function fail(msg: string, code = 1): never {
41 console.error(`ref: ${msg}`);
42 process.exit(code);
43}
45const argv = process.argv.slice(2);
46if (argv.includes('--help') || argv.includes('-h')) {
47 console.log(USAGE);
48 process.exit(0);
49}
50let inFile: string | null = null;
51let niterOverride: number | null = null;
52let tolerance: number | null = null;
53let toleranceLinf: number | null = null;
54const wantJson = argv.includes('--json');
55for (let i = 0; i < argv.length; i++) {
56 const a = argv[i];
57 if (a === '--json') continue;
58 const valued = (name: string): string | null => {
59 if (a === `--${name}`) return argv[++i];
60 if (a.startsWith(`--${name}=`)) return a.slice(name.length + 3);
61 return null;
62 };
63 const inv = valued('in');
64 if (inv !== null) {
65 inFile = inv;
66 continue;
67 }
68 const niterv = valued('niter');
69 if (niterv !== null) {
70 niterOverride = Number(niterv);
71 if (!Number.isInteger(niterOverride) || niterOverride < 0) {
72 fail(`--niter must be an integer >= 0 (got '${niterv}')`, 2);
73 }
74 continue;
75 }
76 const tolLinfv = valued('tolerance-linf');
77 if (tolLinfv !== null) {
78 toleranceLinf = Number(tolLinfv);
79 if (!Number.isFinite(toleranceLinf)) fail(`--tolerance-linf must be a number (got '${tolLinfv}')`, 2);
80 continue;
81 }
82 const tolv = valued('tolerance');
83 if (tolv !== null) {
84 tolerance = Number(tolv);
85 if (!Number.isFinite(tolerance)) fail(`--tolerance must be a number (got '${tolv}')`, 2);
86 continue;
87 }
88 fail(`unrecognized argument '${a}'\n\n${USAGE}`, 2);
89}
90if (!inFile) fail(`--in <file> is required\n\n${USAGE}`, 2);
92const attrsOf = (entity: { attrs: Record<string, { value: unknown }> }): Record<string, unknown> =>
93 Object.fromEntries(Object.entries(entity.attrs).map(([k, v]) => [k, v.value]));
95const numberAttrs = (entity: { attrs: Record<string, { value: unknown }> }): Params =>
96 Object.fromEntries(
97 Object.entries(attrsOf(entity)).map(([k, v]) => [k, Number(v)]),
98 );
100let device: GPUDevice | null = null;
101let session: ModelSession | null = null;
102let h5file: InstanceType<typeof h5wasm.File> | null = null;
104try {
105 await h5wasm.ready;
106 h5file = new h5wasm.File(inFile, 'r');
108 const rootAttrs = attrsOf(h5file);
109 const modelKey = String(rootAttrs.model);
110 const model = mModelByKey(modelKey);
111 if (!model) fail(`unknown model '${modelKey}' in ${inFile}`);
113 const specGroup = h5file.get('spec') as InstanceType<typeof h5wasm.Group>;
114 const specAttrs = attrsOf(specGroup);
115 const geometryKey = String(specAttrs.geometry);
116 const geometryModel = mGeometryByKey(geometryKey);
117 if (!geometryModel) fail(`unknown geometry '${geometryKey}' in ${inFile}`);
119 const lmax = Number(specAttrs.lmax);
120 const steps = Number(specAttrs.steps);
121 const niter = niterOverride ?? Number(specAttrs.niter);
123 const params: Params = {
124 ...defaultParams(model),
125 ...numberAttrs(specGroup.get('params') as InstanceType<typeof h5wasm.Group>),
126 };
127 const geometryParams: Params = {
128 ...defaultGeometryParams(geometryModel),
129 ...numberAttrs(specGroup.get('geometry_params') as InstanceType<typeof h5wasm.Group>),
130 };
132 const geomGroup = h5file.get('geometry') as InstanceType<typeof h5wasm.Group>;
133 const fileGeom = {
134 X: (geomGroup.get('Gx') as InstanceType<typeof h5wasm.Dataset>).value as Float32Array,
135 Y: (geomGroup.get('Gy') as InstanceType<typeof h5wasm.Dataset>).value as Float32Array,
136 Z: (geomGroup.get('Gz') as InstanceType<typeof h5wasm.Dataset>).value as Float32Array,
137 };
139 const initialGroup = h5file.get('initial') as InstanceType<typeof h5wasm.Group>;
140 const finalGroup = h5file.get('final') as InstanceType<typeof h5wasm.Group>;
141 const fileInitial: Record<string, Float32Array> = {};
142 const fileFinal: Record<string, Float32Array> = {};
143 for (const name of model.state) {
144 fileInitial[name] = (initialGroup.get(name) as InstanceType<typeof h5wasm.Dataset>)
145 .value as Float32Array;
146 fileFinal[name] = (finalGroup.get(name) as InstanceType<typeof h5wasm.Dataset>)
147 .value as Float32Array;
148 }
150 h5file.close();
151 h5file = null;
153 const runtime = await installWebGpu();
154 device = await requestShtDevice().catch((e: unknown) => {
155 throw new Error(`${errMsg(e)}\n${NO_ADAPTER_HINT}`);
156 });
157 const adapter = await describeAdapter(device);
159 session = await ModelSession.create({
160 device,
161 model,
162 params,
163 lmax,
164 geometry: geometryModel,
165 geometryParams,
166 niter,
167 });
169 const errorOf = (a: Float32Array, b: Float32Array) => ({ relL2: relL2(a, b), relLinf: relLinf(a, b) });
171 const geometryError = {
172 Gx: errorOf(session.geometry.X, fileGeom.X),
173 Gy: errorOf(session.geometry.Y, fileGeom.Y),
174 Gz: errorOf(session.geometry.Z, fileGeom.Z),
175 };
177 session.loadState(fileInitial);
178 session.step(steps);
180 const stateError: Record<string, { relL2: number; relLinf: number }> = {};
181 for (const name of model.state) {
182 // Sequential: GpuModel.read() shares one staging buffer across calls.
183 const ours = await session.read(name);
184 stateError[name] = errorOf(ours, fileFinal[name]);
185 }
187 const allErrors = [...Object.values(geometryError), ...Object.values(stateError)];
188 const worstL2 = Math.max(...allErrors.map((e) => e.relL2));
189 const worstLinf = Math.max(...allErrors.map((e) => e.relLinf));
190 const passL2 = tolerance === null ? null : worstL2 < tolerance;
191 const passLinf = toleranceLinf === null ? null : worstLinf < toleranceLinf;
192 const checks = [passL2, passLinf].filter((p): p is boolean => p !== null);
193 const pass = checks.length === 0 ? null : checks.every(Boolean);
195 if (wantJson) {
196 console.log(
197 JSON.stringify(
198 {
199 in: inFile,
200 model: model.key,
201 geometry: geometryModel.key,
202 grid: { lmax, nlm: session.sht.nlm },
203 niter,
204 steps,
205 dt: params.dt,
206 T: steps * (params.dt ?? 0),
207 backend: { adapter, runtime, precision: 'fp32' },
208 geometryError,
209 stateError,
210 worstL2,
211 worstLinf,
212 tolerance,
213 toleranceLinf,
214 passL2,
215 passLinf,
216 pass,
217 },
218 null,
219 2,
220 ),
221 );
222 } else {
223 console.log(`ref: ${inFile}`);
224 console.log(
225 ` model ${model.label} (${model.state.join(', ')})\n` +
226 ` geometry ${geometryModel.label} ` +
227 geometryModel.params.map((p) => `${p.key}=${geometryParams[p.key]}`).join(' ') +
228 `\n grid lmax ${lmax} · nlm ${session.sht.nlm}\n` +
229 ` niter ${niter}${niterOverride !== null ? ` (file: ${specAttrs.niter})` : ''}\n` +
230 ` run ${steps} steps, dt=${params.dt} (T=${(steps * (params.dt ?? 0)).toFixed(2)})\n`,
231 );
232 const fmtErr = (v: { relL2: number; relLinf: number }) =>
233 `relL2 ${v.relL2.toExponential(3)} relLinf ${v.relLinf.toExponential(3)}`;
234 console.log(` geometry check (regenerated vs file):`);
235 for (const [k, v] of Object.entries(geometryError)) console.log(` ${k} ${fmtErr(v)}`);
236 console.log(`\n final state (this run vs file):`);
237 for (const [k, v] of Object.entries(stateError)) console.log(` ${k} ${fmtErr(v)}`);
238 if (tolerance !== null) {
239 console.log(
240 `\n worst relL2 ${worstL2.toExponential(3)} vs tolerance ${tolerance.toExponential(3)}: ` +
241 (passL2 ? 'PASS' : 'FAIL'),
242 );
243 }
244 if (toleranceLinf !== null) {
245 console.log(
246 ` worst relLinf ${worstLinf.toExponential(3)} vs tolerance-linf ${toleranceLinf.toExponential(3)}: ` +
247 (passLinf ? 'PASS' : 'FAIL'),
248 );
249 }
250 }
252 session.destroy();
253 device.destroy();
254 process.exit(pass === false ? 1 : 0);
255} catch (e) {
256 h5file?.close();
257 session?.destroy();
258 device?.destroy();
259 fail(errMsg(e));
260}