1// GPU-versus-CPU self-test: one time step of the WGSL kernel on a smooth,
2// strongly nonlinear random state is compared with the float64 reference.
3// Runs in the browser (?selftest=1) and is driven by tests/gpu.test.ts.
5import { defaultParams, initialShape, makeMaterial, nodeWeight, type PhysParams } from "./params.ts";
6import { accelerations, addDamping, fields } from "./sim/cpuRef.ts";
7import { Simulation } from "./sim/simulation.ts";
9export interface SelfTestCase {
10 name: string;
11 /** Maximum relative error of the acceleration of each component, over the sheet. */
12 accError: { u: number; v: number; w: number };
13 /** Maximum error of the compensated w update relative to the largest increment dt * w_t. */
14 positionError: number;
15 nodes: number;
16 dt: number;
17}
19export const ACC_TOLERANCE = 1e-3;
20export const POSITION_TOLERANCE = 1e-2;
22const cases: { name: string; overrides: Partial<PhysParams> }[] = [
23 { name: "rectangle, flat", overrides: { shape: "rectangle", bowShape: "none" } },
24 { name: "rectangle, bowed across the width", overrides: { shape: "rectangle", bowShape: "width" } },
25 { name: "ellipse, dish", overrides: { shape: "ellipse", bowShape: "dish" } },
26];
28/** A smooth random field: a few sinusoids across the sheet plus a little white noise. */
29function smoothField(nx: number, ny: number, amplitude: number, noise: number, rnd: () => number): Float64Array {
30 const f = new Float64Array(nx * ny);
31 const terms = Array.from({ length: 6 }, () => ({ kx: 1 + Math.floor(4 * rnd()), ky: 1 + Math.floor(4 * rnd()), px: 2 * Math.PI * rnd(), py: 2 * Math.PI * rnd(), a: rnd() - 0.5 }));
32 for (let j = 0; j < ny; j++) {
33 for (let i = 0; i < nx; i++) {
34 let s = 0;
35 for (const t of terms) s += t.a * Math.cos((t.kx * Math.PI * i) / (nx - 1) + t.px) * Math.cos((t.ky * Math.PI * j) / (ny - 1) + t.py);
36 f[j * nx + i] = amplitude * s + noise * (rnd() - 0.5);
37 }
38 }
39 return f;
40}
42export async function runSelfTest(device: GPUDevice): Promise<SelfTestCase[]> {
43 const results: SelfTestCase[] = [];
44 let seed = 2024;
45 const rnd = () => ((seed = (seed * 1664525 + 1013904223) >>> 0) / 4294967296);
46 for (const c of cases) {
47 const params: PhysParams = { ...defaultParams, ...c.overrides, spacing: 0.012 };
48 const sim = new Simulation(device, params);
49 const g = sim.grid;
50 const n = g.n;
51 const m = makeMaterial(params);
52 const w0 = Float64Array.from(initialShape(params, g));
54 // State in float32 (what the GPU sees), then the same values in float64 for the reference.
55 const pos = new Float32Array(4 * n), vel = new Float32Array(3 * n);
56 const fw = smoothField(g.nx, g.ny, 2e-3, 5e-5, rnd); // about 2.5 thicknesses: strongly nonlinear
57 const fu = smoothField(g.nx, g.ny, 2e-5, 1e-6, rnd);
58 const fv = smoothField(g.nx, g.ny, 2e-5, 1e-6, rnd);
59 const fdw = smoothField(g.nx, g.ny, 0.5, 1e-2, rnd);
60 const fdu = smoothField(g.nx, g.ny, 1e-2, 1e-3, rnd);
61 const fdv = smoothField(g.nx, g.ny, 1e-2, 1e-3, rnd);
62 for (let k = 0; k < n; k++) {
63 if (nodeWeight(g, k % g.nx, Math.floor(k / g.nx)) === 0) continue;
64 pos[k] = fu[k]; pos[n + k] = fv[k]; pos[2 * n + k] = fw[k]; pos[3 * n + k] = 0;
65 vel[k] = fdu[k]; vel[n + k] = fdv[k]; vel[2 * n + k] = fdw[k];
66 }
67 const x = fields(n), v = fields(n);
68 for (let k = 0; k < n; k++) {
69 x.u[k] = pos[k]; x.v[k] = pos[n + k]; x.w[k] = pos[2 * n + k];
70 v.u[k] = vel[k]; v.v[k] = vel[n + k]; v.w[k] = vel[2 * n + k];
71 }
73 sim.writeState(pos, vel);
74 await sim.stepOnce();
75 const out = await sim.readState();
76 sim.destroy();
78 // Reference: a = elastic + losses, v' = v + dt a, w' = w + dt v'.
79 const acc = fields(n);
80 accelerations(g, m, x, acc, w0);
81 addDamping(g, params.sigma0, params.sigma1, params.sigmaIn, v, acc);
82 const dt = sim.dt;
83 const err = { u: 0, v: 0, w: 0 };
84 let scale = { u: 0, v: 0, w: 0 };
85 let posErr = 0, posScale = 0;
86 for (let k = 0; k < n; k++) {
87 if (nodeWeight(g, k % g.nx, Math.floor(k / g.nx)) === 0) continue;
88 const gu = (out.vel[k] - vel[k]) / dt, gv = (out.vel[n + k] - vel[n + k]) / dt, gw = (out.vel[2 * n + k] - vel[2 * n + k]) / dt;
89 err.u = Math.max(err.u, Math.abs(gu - acc.u[k]));
90 err.v = Math.max(err.v, Math.abs(gv - acc.v[k]));
91 err.w = Math.max(err.w, Math.abs(gw - acc.w[k]));
92 scale = { u: Math.max(scale.u, Math.abs(acc.u[k])), v: Math.max(scale.v, Math.abs(acc.v[k])), w: Math.max(scale.w, Math.abs(acc.w[k])) };
93 const vNew = v.w[k] + dt * acc.w[k];
94 const wNew = x.w[k] + dt * vNew;
95 posErr = Math.max(posErr, Math.abs(out.pos[2 * n + k] + out.pos[3 * n + k] - wNew));
96 posScale = Math.max(posScale, Math.abs(dt * vNew));
97 }
98 results.push({
99 name: c.name,
100 accError: { u: err.u / scale.u, v: err.v / scale.v, w: err.w / scale.w },
101 positionError: posErr / posScale,
102 nodes: n,
103 dt,
104 });
105 }
106 return results;
107}
109export function selfTestPassed(r: SelfTestCase): boolean {
110 return r.accError.u < ACC_TOLERANCE && r.accError.v < ACC_TOLERANCE && r.accError.w < ACC_TOLERANCE && r.positionError < POSITION_TOLERANCE;
111}