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 } 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 --json machine-readable output
30 --help
32Runs this repo's solver from the file's exact initial spectral state, to the
33same physical end time, and reports the relative-L2 error of the resulting
34state against the file's final state (and, as a sanity check, of the
35regenerated geometry against the file's own geometry coefficients).`;
37function fail(msg: string, code = 1): never {
38 console.error(`ref: ${msg}`);
39 process.exit(code);
40}
42const argv = process.argv.slice(2);
43if (argv.includes('--help') || argv.includes('-h')) {
44 console.log(USAGE);
45 process.exit(0);
46}
47let inFile: string | null = null;
48let niterOverride: number | null = null;
49let tolerance: number | null = null;
50const wantJson = argv.includes('--json');
51for (let i = 0; i < argv.length; i++) {
52 const a = argv[i];
53 if (a === '--json') continue;
54 const valued = (name: string): string | null => {
55 if (a === `--${name}`) return argv[++i];
56 if (a.startsWith(`--${name}=`)) return a.slice(name.length + 3);
57 return null;
58 };
59 const inv = valued('in');
60 if (inv !== null) {
61 inFile = inv;
62 continue;
63 }
64 const niterv = valued('niter');
65 if (niterv !== null) {
66 niterOverride = Number(niterv);
67 if (!Number.isInteger(niterOverride) || niterOverride < 0) {
68 fail(`--niter must be an integer >= 0 (got '${niterv}')`, 2);
69 }
70 continue;
71 }
72 const tolv = valued('tolerance');
73 if (tolv !== null) {
74 tolerance = Number(tolv);
75 if (!Number.isFinite(tolerance)) fail(`--tolerance must be a number (got '${tolv}')`, 2);
76 continue;
77 }
78 fail(`unrecognized argument '${a}'\n\n${USAGE}`, 2);
79}
80if (!inFile) fail(`--in <file> is required\n\n${USAGE}`, 2);
82const attrsOf = (entity: { attrs: Record<string, { value: unknown }> }): Record<string, unknown> =>
83 Object.fromEntries(Object.entries(entity.attrs).map(([k, v]) => [k, v.value]));
85const numberAttrs = (entity: { attrs: Record<string, { value: unknown }> }): Params =>
86 Object.fromEntries(
87 Object.entries(attrsOf(entity)).map(([k, v]) => [k, Number(v)]),
88 );
90let device: GPUDevice | null = null;
91let session: ModelSession | null = null;
92let h5file: InstanceType<typeof h5wasm.File> | null = null;
94try {
95 await h5wasm.ready;
96 h5file = new h5wasm.File(inFile, 'r');
98 const rootAttrs = attrsOf(h5file);
99 const modelKey = String(rootAttrs.model);
100 const model = mModelByKey(modelKey);
101 if (!model) fail(`unknown model '${modelKey}' in ${inFile}`);
103 const specGroup = h5file.get('spec') as InstanceType<typeof h5wasm.Group>;
104 const specAttrs = attrsOf(specGroup);
105 const geometryKey = String(specAttrs.geometry);
106 const geometryModel = mGeometryByKey(geometryKey);
107 if (!geometryModel) fail(`unknown geometry '${geometryKey}' in ${inFile}`);
109 const lmax = Number(specAttrs.lmax);
110 const steps = Number(specAttrs.steps);
111 const niter = niterOverride ?? Number(specAttrs.niter);
113 const params: Params = {
114 ...defaultParams(model),
115 ...numberAttrs(specGroup.get('params') as InstanceType<typeof h5wasm.Group>),
116 };
117 const geometryParams: Params = {
118 ...defaultGeometryParams(geometryModel),
119 ...numberAttrs(specGroup.get('geometry_params') as InstanceType<typeof h5wasm.Group>),
120 };
122 const geomGroup = h5file.get('geometry') as InstanceType<typeof h5wasm.Group>;
123 const fileGeom = {
124 X: (geomGroup.get('Gx') as InstanceType<typeof h5wasm.Dataset>).value as Float32Array,
125 Y: (geomGroup.get('Gy') as InstanceType<typeof h5wasm.Dataset>).value as Float32Array,
126 Z: (geomGroup.get('Gz') as InstanceType<typeof h5wasm.Dataset>).value as Float32Array,
127 };
129 const initialGroup = h5file.get('initial') as InstanceType<typeof h5wasm.Group>;
130 const finalGroup = h5file.get('final') as InstanceType<typeof h5wasm.Group>;
131 const fileInitial: Record<string, Float32Array> = {};
132 const fileFinal: Record<string, Float32Array> = {};
133 for (const name of model.state) {
134 fileInitial[name] = (initialGroup.get(name) as InstanceType<typeof h5wasm.Dataset>)
135 .value as Float32Array;
136 fileFinal[name] = (finalGroup.get(name) as InstanceType<typeof h5wasm.Dataset>)
137 .value as Float32Array;
138 }
140 h5file.close();
141 h5file = null;
143 const runtime = await installWebGpu();
144 device = await requestShtDevice().catch((e: unknown) => {
145 throw new Error(`${errMsg(e)}\n${NO_ADAPTER_HINT}`);
146 });
147 const adapter = await describeAdapter(device);
149 session = await ModelSession.create({
150 device,
151 model,
152 params,
153 lmax,
154 geometry: geometryModel,
155 geometryParams,
156 niter,
157 });
159 const geometryError = {
160 Gx: relL2(session.geometry.X, fileGeom.X),
161 Gy: relL2(session.geometry.Y, fileGeom.Y),
162 Gz: relL2(session.geometry.Z, fileGeom.Z),
163 };
165 session.loadState(fileInitial);
166 session.step(steps);
168 const stateError: Record<string, number> = {};
169 for (const name of model.state) {
170 // Sequential: GpuModel.read() shares one staging buffer across calls.
171 const ours = await session.read(name);
172 stateError[name] = relL2(ours, fileFinal[name]);
173 }
175 const allErrors = [...Object.values(geometryError), ...Object.values(stateError)];
176 const worst = Math.max(...allErrors);
177 const pass = tolerance === null ? null : worst < tolerance;
179 if (wantJson) {
180 console.log(
181 JSON.stringify(
182 {
183 in: inFile,
184 model: model.key,
185 geometry: geometryModel.key,
186 grid: { lmax, nlm: session.sht.nlm },
187 niter,
188 steps,
189 dt: params.dt,
190 T: steps * (params.dt ?? 0),
191 backend: { adapter, runtime, precision: 'fp32' },
192 geometryError,
193 stateError,
194 tolerance,
195 pass,
196 },
197 null,
198 2,
199 ),
200 );
201 } else {
202 console.log(`ref: ${inFile}`);
203 console.log(
204 ` model ${model.label} (${model.state.join(', ')})\n` +
205 ` geometry ${geometryModel.label} ` +
206 geometryModel.params.map((p) => `${p.key}=${geometryParams[p.key]}`).join(' ') +
207 `\n grid lmax ${lmax} · nlm ${session.sht.nlm}\n` +
208 ` niter ${niter}${niterOverride !== null ? ` (file: ${specAttrs.niter})` : ''}\n` +
209 ` run ${steps} steps, dt=${params.dt} (T=${(steps * (params.dt ?? 0)).toFixed(2)})\n`,
210 );
211 console.log(` geometry check (regenerated vs file, relL2):`);
212 for (const [k, v] of Object.entries(geometryError)) console.log(` ${k} ${v.toExponential(3)}`);
213 console.log(`\n final state (this run vs file, relL2):`);
214 for (const [k, v] of Object.entries(stateError)) console.log(` ${k} ${v.toExponential(3)}`);
215 if (tolerance !== null) {
216 console.log(
217 `\n worst relL2 ${worst.toExponential(3)} vs tolerance ${tolerance.toExponential(3)}: ` +
218 (pass ? 'PASS' : 'FAIL'),
219 );
220 }
221 }
223 session.destroy();
224 device.destroy();
225 process.exit(pass === false ? 1 : 0);
226} catch (e) {
227 h5file?.close();
228 session?.destroy();
229 device?.destroy();
230 fail(errMsg(e));
231}