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, makeMaterial, nodeWeight, type PhysParams } from "./params.ts";
6import { initialShape, parseFeatures } from "./shape.ts";
7import { accelerations, addDamping, fields } from "./sim/cpuRef.ts";
8import { Simulation } from "./sim/simulation.ts";
10export interface SelfTestCase {
11 name: string;
12 /** Maximum relative error of the acceleration of each component, over the sheet. */
13 accError: { u: number; v: number; w: number };
14 /** Maximum error of the compensated w update relative to the largest increment dt * w_t. */
15 positionError: number;
16 nodes: number;
17 dt: number;
18}
20export const ACC_TOLERANCE = 1e-3;
21export const POSITION_TOLERANCE = 1e-2;
23const cases: { name: string; overrides: Partial<PhysParams> }[] = [
24 { name: "rectangle, flat", overrides: { shape: "rectangle", features: [] } },
25 { name: "rectangle, bowed across the width", overrides: { shape: "rectangle" } },
26 { name: "ellipse, dish", overrides: { shape: "ellipse", features: parseFeatures("bow:x,0.0015,0;bow:y,0.0015,0")! } },
27 { name: "rectangle, twist, dent, crease, and random field", overrides: { shape: "rectangle", features: parseFeatures("bow:x,0.002,0.5;twist:0.001;dent:0.3,0.6,0.03,0.0015;crease:0.6,0.4,25,-1.5,0.015;noise:11,0.0005,0.05")! } },
28 { name: "rectangle, bowed, in-plane mass x4", overrides: { shape: "rectangle", inPlaneMass: 4 } },
29];
31/** A smooth random field: a few sinusoids across the sheet plus a little white noise. */
32function smoothField(nx: number, ny: number, amplitude: number, noise: number, rnd: () => number): Float64Array {
33 const f = new Float64Array(nx * ny);
34 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 }));
35 for (let j = 0; j < ny; j++) {
36 for (let i = 0; i < nx; i++) {
37 let s = 0;
38 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);
39 f[j * nx + i] = amplitude * s + noise * (rnd() - 0.5);
40 }
41 }
42 return f;
43}
45export async function runSelfTest(device: GPUDevice): Promise<SelfTestCase[]> {
46 const results: SelfTestCase[] = [];
47 let seed = 2024;
48 const rnd = () => ((seed = (seed * 1664525 + 1013904223) >>> 0) / 4294967296);
49 for (const c of cases) {
50 const params: PhysParams = { ...defaultParams, ...c.overrides, spacing: 0.012 };
51 const sim = new Simulation(device, params);
52 const g = sim.grid;
53 const n = g.n;
54 const m = makeMaterial(params);
55 const w0 = Float64Array.from(initialShape(params, g));
57 // State in float32 (what the GPU sees), then the same values in float64 for the reference.
58 const pos = new Float32Array(4 * n), vel = new Float32Array(3 * n);
59 const fw = smoothField(g.nx, g.ny, 2e-3, 5e-5, rnd); // about 2.5 thicknesses: strongly nonlinear
60 const fu = smoothField(g.nx, g.ny, 2e-5, 1e-6, rnd);
61 const fv = smoothField(g.nx, g.ny, 2e-5, 1e-6, rnd);
62 const fdw = smoothField(g.nx, g.ny, 0.5, 1e-2, rnd);
63 const fdu = smoothField(g.nx, g.ny, 1e-2, 1e-3, rnd);
64 const fdv = smoothField(g.nx, g.ny, 1e-2, 1e-3, rnd);
65 for (let k = 0; k < n; k++) {
66 if (nodeWeight(g, k % g.nx, Math.floor(k / g.nx)) === 0) continue;
67 pos[k] = fu[k]; pos[n + k] = fv[k]; pos[2 * n + k] = fw[k]; pos[3 * n + k] = 0;
68 vel[k] = fdu[k]; vel[n + k] = fdv[k]; vel[2 * n + k] = fdw[k];
69 }
70 const x = fields(n), v = fields(n);
71 for (let k = 0; k < n; k++) {
72 x.u[k] = pos[k]; x.v[k] = pos[n + k]; x.w[k] = pos[2 * n + k];
73 v.u[k] = vel[k]; v.v[k] = vel[n + k]; v.w[k] = vel[2 * n + k];
74 }
76 sim.writeState(pos, vel);
77 await sim.stepOnce();
78 const out = await sim.readState();
79 sim.destroy();
81 // Reference: a = elastic (in-plane scaled) + losses, v' = v + dt a, w' = w + dt v'.
82 const acc = fields(n);
83 accelerations(g, m, x, acc, w0);
84 for (let k = 0; k < n; k++) { acc.u[k] /= params.inPlaneMass; acc.v[k] /= params.inPlaneMass; }
85 addDamping(g, params.sigma0, params.sigma1, params.sigmaIn, v, acc);
86 const dt = sim.dt;
87 const err = { u: 0, v: 0, w: 0 };
88 let scale = { u: 0, v: 0, w: 0 };
89 let posErr = 0, posScale = 0;
90 for (let k = 0; k < n; k++) {
91 if (nodeWeight(g, k % g.nx, Math.floor(k / g.nx)) === 0) continue;
92 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;
93 err.u = Math.max(err.u, Math.abs(gu - acc.u[k]));
94 err.v = Math.max(err.v, Math.abs(gv - acc.v[k]));
95 err.w = Math.max(err.w, Math.abs(gw - acc.w[k]));
96 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])) };
97 const vNew = v.w[k] + dt * acc.w[k];
98 const wNew = x.w[k] + dt * vNew;
99 posErr = Math.max(posErr, Math.abs(out.pos[2 * n + k] + out.pos[3 * n + k] - wNew));
100 posScale = Math.max(posScale, Math.abs(dt * vNew));
101 }
102 results.push({
103 name: c.name,
104 accError: { u: err.u / scale.u, v: err.v / scale.v, w: err.w / scale.w },
105 positionError: posErr / posScale,
106 nodes: n,
107 dt,
108 });
109 }
110 return results;
111}
113export function selfTestPassed(r: SelfTestCase): boolean {
114 return r.accError.u < ACC_TOLERANCE && r.accError.v < ACC_TOLERANCE && r.accError.w < ACC_TOLERANCE && r.positionError < POSITION_TOLERANCE;
115}