1ee087eCymbal designer: scripted design, thickness field, mount, single strikeJeremy Magland 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// for several of the preset designs (flat and bowed rectangles, tapered
4// discs with a hole and a mount, an ellipse). Runs in the browser
5// (?selftest=1) and is driven by tests/gpu.test.ts.
7import { defaultParams, nodeWeight, type SimParams } from "./params.ts";
8import { compileDesign, evaluatePlate } from "./design.ts";
9import { PRESETS } from "./presets.ts";
10import { accelerations, addDamping, fields } from "./sim/cpuRef.ts";
11import { Simulation } from "./sim/simulation.ts";
13export interface SelfTestCase {
14 name: string;
15 /** Maximum relative error of the acceleration of each component, over the plate. */
16 accError: { u: number; v: number; w: number };
17 /** Maximum error of the compensated w update relative to the largest increment dt * w_t. */
18 positionError: number;
19 nodes: number;
20 dt: number;
21}
23export const ACC_TOLERANCE = 1e-3;
24export const POSITION_TOLERANCE = 1e-2;
26const cases: { name: string; preset: string; overrides?: Partial<SimParams> }[] = [
27 { name: "flat rectangle", preset: "flat-rectangle" },
28 { name: "bowed rectangle", preset: "bowed-rectangle" },
29 { name: "16 inch crash (taper, bell, hole, mount)", preset: "crash-16" },
30 { name: "18 inch china", preset: "china-18" },
31 { name: "ellipse with a twist", preset: "ellipse-twist" },
32 { name: "bowed rectangle, in-plane mass x4", preset: "bowed-rectangle", overrides: { inPlaneMass: 4 } },
33];
35/** A smooth random field: a few sinusoids across the plate plus a little white noise. */
36function smoothField(nx: number, ny: number, amplitude: number, noise: number, rnd: () => number): Float64Array {
37 const f = new Float64Array(nx * ny);
38 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 }));
39 for (let j = 0; j < ny; j++) {
40 for (let i = 0; i < nx; i++) {
41 let s = 0;
42 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);
43 f[j * nx + i] = amplitude * s + noise * (rnd() - 0.5);
44 }
45 }
46 return f;
47}
49export async function runSelfTest(device: GPUDevice): Promise<SelfTestCase[]> {
50 const results: SelfTestCase[] = [];
51 let seed = 2024;
52 const rnd = () => ((seed = (seed * 1664525 + 1013904223) >>> 0) / 4294967296);
53 for (const c of cases) {
54 const preset = PRESETS.find((p) => p.id === c.preset)!;
55 const params: SimParams = { ...defaultParams, ...c.overrides, spacing: 0.012 };
56 const plate = evaluatePlate(compileDesign(preset.source), params.spacing);
57 const sim = new Simulation(device, plate, params);
58 const g = sim.grid;
59 const n = g.n;
61 // State in float32 (what the GPU sees), then the same values in float64 for the reference.
62 const pos = new Float32Array(4 * n), vel = new Float32Array(3 * n);
63 const fw = smoothField(g.nx, g.ny, 2e-3, 5e-5, rnd); // about 2.5 thicknesses: strongly nonlinear
64 const fu = smoothField(g.nx, g.ny, 2e-5, 1e-6, rnd);
65 const fv = smoothField(g.nx, g.ny, 2e-5, 1e-6, rnd);
66 const fdw = smoothField(g.nx, g.ny, 0.5, 1e-2, rnd);
67 const fdu = smoothField(g.nx, g.ny, 1e-2, 1e-3, rnd);
68 const fdv = smoothField(g.nx, g.ny, 1e-2, 1e-3, rnd);
69 for (let k = 0; k < n; k++) {
70 if (nodeWeight(g, k % g.nx, Math.floor(k / g.nx)) === 0) continue;
71 pos[k] = fu[k]; pos[n + k] = fv[k]; pos[2 * n + k] = fw[k]; pos[3 * n + k] = 0;
72 vel[k] = fdu[k]; vel[n + k] = fdv[k]; vel[2 * n + k] = fdw[k];
73 }
74 const x = fields(n), v = fields(n);
75 for (let k = 0; k < n; k++) {
76 x.u[k] = pos[k]; x.v[k] = pos[n + k]; x.w[k] = pos[2 * n + k];
77 v.u[k] = vel[k]; v.v[k] = vel[n + k]; v.w[k] = vel[2 * n + k];
78 }
80 sim.writeState(pos, vel);
81 await sim.stepOnce();
82 const out = await sim.readState();
83 sim.destroy();
85 // Reference: a = elastic and mount (in-plane scaled) + losses, v' = v + dt a, w' = w + dt v'.
86 const acc = fields(n);
87 accelerations(plate, x, acc);
88 for (let k = 0; k < n; k++) { acc.u[k] /= params.inPlaneMass; acc.v[k] /= params.inPlaneMass; }
89 addDamping(plate, v, acc);
90 const dt = sim.dt;
91 const err = { u: 0, v: 0, w: 0 };
92 let scale = { u: 0, v: 0, w: 0 };
93 let posErr = 0, posScale = 0;
94 for (let k = 0; k < n; k++) {
95 if (nodeWeight(g, k % g.nx, Math.floor(k / g.nx)) === 0) continue;
96 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;
97 err.u = Math.max(err.u, Math.abs(gu - acc.u[k]));
98 err.v = Math.max(err.v, Math.abs(gv - acc.v[k]));
99 err.w = Math.max(err.w, Math.abs(gw - acc.w[k]));
100 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])) };
101 const vNew = v.w[k] + dt * acc.w[k];
102 const wNew = x.w[k] + dt * vNew;
103 posErr = Math.max(posErr, Math.abs(out.pos[2 * n + k] + out.pos[3 * n + k] - wNew));
104 posScale = Math.max(posScale, Math.abs(dt * vNew));
105 }
106 results.push({
107 name: c.name,
108 accError: { u: err.u / scale.u, v: err.v / scale.v, w: err.w / scale.w },
109 positionError: posErr / posScale,
110 nodes: n,
111 dt,
112 });
113 }
114 return results;
115}
117export function selfTestPassed(r: SelfTestCase): boolean {
118 return r.accError.u < ACC_TOLERANCE && r.accError.v < ACC_TOLERANCE && r.accError.w < ACC_TOLERANCE && r.positionError < POSITION_TOLERANCE;
119}