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 { name: "rectangle, bowed, in-plane mass x4", overrides: { shape: "rectangle", bowShape: "width", inPlaneMass: 4 } },
27];
29/** A smooth random field: a few sinusoids across the sheet plus a little white noise. */
30function smoothField(nx: number, ny: number, amplitude: number, noise: number, rnd: () => number): Float64Array {
31 const f = new Float64Array(nx * ny);
32 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 }));
33 for (let j = 0; j < ny; j++) {
34 for (let i = 0; i < nx; i++) {
35 let s = 0;
36 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);
37 f[j * nx + i] = amplitude * s + noise * (rnd() - 0.5);
38 }
39 }
40 return f;
41}
43export async function runSelfTest(device: GPUDevice): Promise<SelfTestCase[]> {
44 const results: SelfTestCase[] = [];
45 let seed = 2024;
46 const rnd = () => ((seed = (seed * 1664525 + 1013904223) >>> 0) / 4294967296);
47 for (const c of cases) {
48 const params: PhysParams = { ...defaultParams, ...c.overrides, spacing: 0.012 };
49 const sim = new Simulation(device, params);
50 const g = sim.grid;
51 const n = g.n;
52 const m = makeMaterial(params);
53 const w0 = Float64Array.from(initialShape(params, g));
55 // State in float32 (what the GPU sees), then the same values in float64 for the reference.
56 const pos = new Float32Array(4 * n), vel = new Float32Array(3 * n);
57 const fw = smoothField(g.nx, g.ny, 2e-3, 5e-5, rnd); // about 2.5 thicknesses: strongly nonlinear
58 const fu = smoothField(g.nx, g.ny, 2e-5, 1e-6, rnd);
59 const fv = smoothField(g.nx, g.ny, 2e-5, 1e-6, rnd);
60 const fdw = smoothField(g.nx, g.ny, 0.5, 1e-2, rnd);
61 const fdu = smoothField(g.nx, g.ny, 1e-2, 1e-3, rnd);
62 const fdv = smoothField(g.nx, g.ny, 1e-2, 1e-3, rnd);
63 for (let k = 0; k < n; k++) {
64 if (nodeWeight(g, k % g.nx, Math.floor(k / g.nx)) === 0) continue;
65 pos[k] = fu[k]; pos[n + k] = fv[k]; pos[2 * n + k] = fw[k]; pos[3 * n + k] = 0;
66 vel[k] = fdu[k]; vel[n + k] = fdv[k]; vel[2 * n + k] = fdw[k];
67 }
68 const x = fields(n), v = fields(n);
69 for (let k = 0; k < n; k++) {
70 x.u[k] = pos[k]; x.v[k] = pos[n + k]; x.w[k] = pos[2 * n + k];
71 v.u[k] = vel[k]; v.v[k] = vel[n + k]; v.w[k] = vel[2 * n + k];
72 }
74 sim.writeState(pos, vel);
75 await sim.stepOnce();
76 const out = await sim.readState();
77 sim.destroy();
79 // Reference: a = elastic (in-plane scaled) + losses, v' = v + dt a, w' = w + dt v'.
80 const acc = fields(n);
81 accelerations(g, m, x, acc, w0);
82 for (let k = 0; k < n; k++) { acc.u[k] /= params.inPlaneMass; acc.v[k] /= params.inPlaneMass; }
83 addDamping(g, params.sigma0, params.sigma1, params.sigmaIn, v, acc);
84 const dt = sim.dt;
85 const err = { u: 0, v: 0, w: 0 };
86 let scale = { u: 0, v: 0, w: 0 };
87 let posErr = 0, posScale = 0;
88 for (let k = 0; k < n; k++) {
89 if (nodeWeight(g, k % g.nx, Math.floor(k / g.nx)) === 0) continue;
90 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;
91 err.u = Math.max(err.u, Math.abs(gu - acc.u[k]));
92 err.v = Math.max(err.v, Math.abs(gv - acc.v[k]));
93 err.w = Math.max(err.w, Math.abs(gw - acc.w[k]));
94 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])) };
95 const vNew = v.w[k] + dt * acc.w[k];
96 const wNew = x.w[k] + dt * vNew;
97 posErr = Math.max(posErr, Math.abs(out.pos[2 * n + k] + out.pos[3 * n + k] - wNew));
98 posScale = Math.max(posScale, Math.abs(dt * vNew));
99 }
100 results.push({
101 name: c.name,
102 accError: { u: err.u / scale.u, v: err.v / scale.v, w: err.w / scale.w },
103 positionError: posErr / posScale,
104 nodes: n,
105 dt,
106 });
107 }
108 return results;
109}
111export function selfTestPassed(r: SelfTestCase): boolean {
112 return r.accError.u < ACC_TOLERANCE && r.accError.v < ACC_TOLERANCE && r.accError.w < ACC_TOLERANCE && r.positionError < POSITION_TOLERANCE;
113}