5778036cymbalsim: scripted plate and cymbal designer with a modal linear model and a von Karman GPU modelJeremy Magland 1// WebGPU time-stepping driver of the full (geometrically nonlinear)
2// simulation. The run advances in batches of one movie frame (1/fps s). A
3// batch encodes one compute dispatch per time step (with a per-step uniform
4// addressed by dynamic offset), one retarded-time listener sum over the
5// batch's audio samples, and one extraction of the transverse displacement
6// at the end. Audio, frame, and stick state are read back after every batch;
7// one batch is kept encoded ahead of the one being read back so that the
8// GPU never waits for JavaScript. The recording starts at the impact, with
9// the plate at rest in its stress-free shape.
11import {
12 type Plate,
13 type SimParams,
14 CFL,
15 HIST_MARGIN,
16 massProperties,
17 strikePatch,
18 strikePoint,
19 thicknessRange,
20 nodeWeights,
21 listenerTerms,
22} from "../params.ts";
23import { type FullLoss, maxAngularFrequency } from "../model/reference.ts";
24import commonWgsl from "./shaders/common.wgsl?raw";
25import stepWgsl from "./shaders/step.wgsl?raw";
26import audioWgsl from "./shaders/audio.wgsl?raw";
27import extractWgsl from "./shaders/extract.wgsl?raw";
29export type Phase = "recording" | "done";
31export interface StickState {
32 x: number;
33 v: number;
34 force: number;
35 maxForce: number;
36 withdrawn: boolean;
37}
39export interface SimStatus {
40 phase: Phase;
41 /** Simulation time since the impact of the state read back so far [s]. */
42 t: number;
43 totalSteps: number;
44 samples: number;
45 frames: number;
46 peakPressure: number;
47 /** Largest transverse displacement seen in the frames [m]. */
48 wmax: number;
49 delayClamped: boolean;
50 stick: StickState;
51}
5778036cymbalsim: scripted plate and cymbal designer with a modal linear model and a von Karman GPU modelJeremy Magland 54const STEP_LEN = 12;
55const STICK_LEN = 16;
56/** Nodes per workgroup along each axis in step.wgsl. */
57export const TILE = 12;
59interface Pending {
60 staging: GPUBuffer;
61 steps: number;
62 nAudio: number;
63 stickParity: number;
64}
66export class Simulation {
67 readonly plate: Plate;
68 readonly dt: number;
69 readonly stepsPerSample: number;
70 readonly samplesPerBatch: number;
71 readonly stepsPerBatch: number;
72 readonly omegaMax: number;
73 readonly histLen: number;
74 readonly mass: number;
76 audio: Float32Array = new Float32Array(0);
77 samples = 0;
78 /** Movie frames of w [m]; frame k is at t = k / fps. */
79 frames: Float32Array[] = [];
80 phase: Phase = "recording";
81 totalSteps = 0;
82 peakPressure = 0;
83 wmax = 0;
84 delayClamped = false;
86 private stick: StickState;
87 private inFlight: Pending | null = null;
88 private device: GPUDevice;
89 private params: SimParams;
90 private loss: FullLoss;
91 private paramsBuf: GPUBuffer;
92 private stepBuf: GPUBuffer;
93 private stepStride: number;
94 private stepData: ArrayBuffer;
95 private pos: GPUBuffer[];
96 private vel: GPUBuffer[];
97 private hist: GPUBuffer;
98 private stickBuf: GPUBuffer;
99 private nd: GPUBuffer;
100 private lw: GPUBuffer;
101 private audioOut: GPUBuffer;
102 private frameBuf: GPUBuffer;
103 private staging: GPUBuffer[];
104 private stagingIdx = 0;
105 private cur = 0;
106 private stepPipe: GPUComputePipeline;
107 private audioPipe: GPUComputePipeline;
108 private extractPipe: GPUComputePipeline;
109 private stepBG: GPUBindGroup[];
110 private audioBG: GPUBindGroup;
111 private extractBG: GPUBindGroup[];
113 constructor(device: GPUDevice, plate: Plate, params: SimParams, loss: FullLoss) {
114 this.device = device;
115 this.plate = plate;
116 this.params = { ...params };
117 this.loss = loss;
118 const g = plate.grid, m = plate.mat;
120 // Stable time step, rounded so that an integer number of steps make one audio sample.
121 this.omegaMax = maxAngularFrequency(g, m, thicknessRange(plate).max);
122 const dtMax = (CFL * 2) / this.omegaMax;
123 this.stepsPerSample = Math.ceil(1 / (params.fs * dtMax));
124 this.dt = 1 / (params.fs * this.stepsPerSample);
125 this.samplesPerBatch = Math.max(1, Math.round(params.fs / params.fps));
126 this.stepsPerBatch = this.samplesPerBatch * this.stepsPerSample;
127 this.histLen = this.samplesPerBatch + HIST_MARGIN;
129 const n = g.n;
130 const f32 = (size: number, usage = GPUBufferUsage.STORAGE | GPUBufferUsage.COPY_SRC | GPUBufferUsage.COPY_DST) =>
131 device.createBuffer({ size: Math.max(16, size * 4), usage });
132 this.pos = [f32(4 * n), f32(4 * n)];
133 this.vel = [f32(3 * n), f32(3 * n)];
134 this.hist = f32(this.histLen * n);
135 this.stickBuf = f32(STICK_LEN);
136 const weights = nodeWeights(g);
137 const packed = new Float32Array(4 * n);
138 for (let k = 0; k < n; k++) {
139 packed[4 * k] = plate.w0[k]; packed[4 * k + 1] = weights[k]; packed[4 * k + 2] = plate.h[k];
140 }
141 this.nd = f32(4 * n);
142 device.queue.writeBuffer(this.nd, 0, packed);
143 this.lw = f32(2 * n);
144 this.audioOut = f32(this.samplesPerBatch);
145 this.frameBuf = f32(n);
146 const stagingSize = (this.samplesPerBatch + n + STICK_LEN) * 4;
147 this.staging = [0, 1].map(() => device.createBuffer({ size: stagingSize, usage: GPUBufferUsage.MAP_READ | GPUBufferUsage.COPY_DST }));
149 const { mass, Ix, Iy } = massProperties(plate);
150 this.mass = mass;
151 const cp = strikePatch(plate, params);
152 let i0 = g.nx, i1 = 0, j0 = g.ny, j1 = 0;
153 for (const k of cp.nodes) {
154 const i = k % g.nx, j = Math.floor(k / g.nx);
155 i0 = Math.min(i0, i); i1 = Math.max(i1, i); j0 = Math.min(j0, j); j1 = Math.max(j1, j);
156 }
157 const stick = new Float32Array(STICK_LEN);
158 stick[1] = -params.strikeVelocity; stick[9] = -params.strikeVelocity;
159 device.queue.writeBuffer(this.stickBuf, 0, stick);
160 this.stick = { x: 0, v: -params.strikeVelocity, force: 0, maxForce: 0, withdrawn: false };
162 const sp = strikePoint(plate, params);
163 const pf = new Float32Array(PARAMS_LEN);
164 const pu = new Uint32Array(pf.buffer);
165 pu[0] = g.nx; pu[1] = g.ny; pu[2] = n; pu[3] = this.histLen;
166 pf[4] = g.hx; pf[5] = g.hy; pf[6] = this.dt; pf[7] = m.nu;
167 pf[8] = m.E / (12 * (1 - m.nu * m.nu)); pf[9] = m.E / (1 - m.nu * m.nu); pf[10] = m.E / (2 * (1 + m.nu)); pf[11] = m.rho;
168 pf[12] = mass; pf[13] = Ix; pf[14] = Iy; pf[15] = params.fs;
169 pf[16] = sp.x; pf[17] = sp.y; pf[18] = params.contactRadius ** 2; pf[19] = cp.mass;
170 pu[20] = cp.nearest; pf[21] = params.contactStiffness; pf[22] = params.contactExponent; pf[23] = params.stickMass;
171 pu[24] = i0; pu[25] = i1; pu[26] = j0; pu[27] = j1;
172 pf[28] = cp.cx; pf[29] = cp.cy;
173 this.paramsBuf = device.createBuffer({ size: pf.byteLength, usage: GPUBufferUsage.UNIFORM | GPUBufferUsage.COPY_DST });
174 device.queue.writeBuffer(this.paramsBuf, 0, pf);
176 this.stepStride = Math.max(STEP_LEN * 4, device.limits.minUniformBufferOffsetAlignment);
177 this.stepData = new ArrayBuffer(this.stepStride * this.stepsPerBatch);
178 this.stepBuf = device.createBuffer({ size: this.stepData.byteLength, usage: GPUBufferUsage.UNIFORM | GPUBufferUsage.COPY_DST });
180 const shader = (code: string) => device.createShaderModule({ code: commonWgsl + code });
181 const entry = (binding: number, type: GPUBufferBindingType, hasDynamicOffset = false): GPUBindGroupLayoutEntry => ({
182 binding, visibility: GPUShaderStage.COMPUTE, buffer: { type, hasDynamicOffset },
183 });
184 const pipeline = (module: GPUShaderModule, layout: GPUBindGroupLayout) =>
185 device.createComputePipeline({ layout: device.createPipelineLayout({ bindGroupLayouts: [layout] }), compute: { module, entryPoint: "main" } });
186 const stepLayout = device.createBindGroupLayout({
2481a63Remove the standJeremy Magland 187 entries: [entry(0, "uniform"), entry(1, "uniform", true), entry(2, "read-only-storage"), entry(3, "read-only-storage"), entry(4, "storage"), entry(5, "storage"), entry(6, "storage"), entry(7, "storage"), entry(8, "read-only-storage")],
5778036cymbalsim: scripted plate and cymbal designer with a modal linear model and a von Karman GPU modelJeremy Magland 188 });
189 this.stepPipe = pipeline(shader(stepWgsl), stepLayout);
190 const audioLayout = device.createBindGroupLayout({
191 entries: [entry(0, "uniform"), entry(1, "uniform", true), entry(2, "read-only-storage"), entry(3, "read-only-storage"), entry(4, "storage"), entry(8, "read-only-storage")],
192 });
193 this.audioPipe = pipeline(shader(audioWgsl), audioLayout);
194 const extractLayout = device.createBindGroupLayout({ entries: [entry(0, "uniform"), entry(1, "read-only-storage"), entry(2, "storage"), entry(8, "read-only-storage")] });
195 this.extractPipe = pipeline(shader(extractWgsl), extractLayout);
197 const spb = { buffer: this.stepBuf, size: STEP_LEN * 4 };
198 this.stepBG = [0, 1].map((p) =>
199 device.createBindGroup({
200 layout: stepLayout,
201 entries: [
202 { binding: 0, resource: { buffer: this.paramsBuf } },
203 { binding: 1, resource: spb },
204 { binding: 2, resource: { buffer: this.pos[p] } },
205 { binding: 3, resource: { buffer: this.vel[p] } },
206 { binding: 4, resource: { buffer: this.pos[1 - p] } },
207 { binding: 5, resource: { buffer: this.vel[1 - p] } },
208 { binding: 6, resource: { buffer: this.hist } },
209 { binding: 7, resource: { buffer: this.stickBuf } },
210 { binding: 8, resource: { buffer: this.nd } },
211 ],
212 }),
213 );
214 this.audioBG = device.createBindGroup({
215 layout: audioLayout,
216 entries: [
217 { binding: 0, resource: { buffer: this.paramsBuf } },
218 { binding: 1, resource: spb },
219 { binding: 2, resource: { buffer: this.hist } },
220 { binding: 3, resource: { buffer: this.lw } },
221 { binding: 4, resource: { buffer: this.audioOut } },
222 { binding: 8, resource: { buffer: this.nd } },
223 ],
224 });
225 this.extractBG = [0, 1].map((p) =>
226 device.createBindGroup({
227 layout: extractLayout,
228 entries: [
229 { binding: 0, resource: { buffer: this.paramsBuf } },
230 { binding: 1, resource: { buffer: this.pos[p] } },
231 { binding: 2, resource: { buffer: this.frameBuf } },
232 { binding: 8, resource: { buffer: this.nd } },
233 ],
234 }),
235 );
236 this.setListener(params);
237 this.frames.push(new Float32Array(n));
238 }
240 /** Recomputes the per-node weights and delays of the listener sum; takes effect on later samples. */
241 setListener(p: SimParams): void {
242 const g = this.plate.grid;
243 this.params.listenerDist = p.listenerDist; this.params.listenerAngle = p.listenerAngle; this.params.listenerAzimuth = p.listenerAzimuth;
244 const { weight, delay } = listenerTerms(this.plate, this.params);
245 const data = new Float32Array(2 * g.n);
246 let dmin = Infinity;
247 for (let k = 0; k < g.n; k++) if (weight[k] > 0) dmin = Math.min(dmin, delay[k]);
248 this.delayClamped = false;
249 for (let k = 0; k < g.n; k++) {
250 data[2 * k] = weight[k];
251 let d = weight[k] > 0 ? (delay[k] - dmin) * this.params.fs : 0;
252 if (d > HIST_MARGIN - 3) { d = HIST_MARGIN - 3; this.delayClamped = true; }
253 data[2 * k + 1] = d;
254 }
255 this.device.queue.writeBuffer(this.lw, 0, data);
256 }
258 setDuration(duration: number): void {
259 this.params.duration = duration;
260 if (this.canContinue()) this.phase = "recording";
261 }
262 canContinue(): boolean {
263 return this.phase === "done" && this.samples < this.targetSamples();
264 }
265 stop(): void {
266 this.phase = "done";
267 }
268 targetSamples(): number {
269 return Math.round(this.params.duration * this.params.fs);
270 }
271 status(): SimStatus {
272 const ahead = this.inFlight ? this.inFlight.steps * this.dt : 0;
273 return { phase: this.phase, t: this.totalSteps * this.dt - ahead, totalSteps: this.totalSteps, samples: this.samples, frames: this.frames.length, peakPressure: this.peakPressure, wmax: this.wmax, delayClamped: this.delayClamped, stick: this.stick };
274 }
275 recordedAudio(): Float32Array {
276 return this.audio.subarray(0, this.samples);
277 }
278 private stepWorkgroups(): [number, number] {
279 return [Math.ceil(this.plate.grid.nx / TILE), Math.ceil(this.plate.grid.ny / TILE)];
280 }
281 private appendAudio(raw: Float32Array): void {
282 const need = this.samples + raw.length;
283 if (this.audio.length < need) {
284 const a = new Float32Array(Math.max(need, 2 * this.audio.length, this.params.fs));
285 a.set(this.audio.subarray(0, this.samples));
286 this.audio = a;
287 }
288 this.audio.set(raw, this.samples);
289 for (let k = 0; k < raw.length; k++) this.peakPressure = Math.max(this.peakPressure, Math.abs(raw[k]));
290 this.samples = need;
291 }
292 private writeStepParams(s: number, fields: { flags: number; slot: number; outIdx: number }): void {
293 const o = (s * this.stepStride) / 4;
294 const u32 = new Uint32Array(this.stepData), i32 = new Int32Array(this.stepData), f32 = new Float32Array(this.stepData);
295 u32[o] = this.totalSteps + s;
296 u32[o + 1] = fields.flags;
297 i32[o + 2] = fields.slot;
298 i32[o + 3] = fields.outIdx;
299 f32[o + 4] = this.loss.sigma0;
300 f32[o + 5] = this.loss.sigma1;
301 f32[o + 6] = this.loss.sigmaIn;
302 u32[o + 7] = (this.totalSteps + s) & 1;
303 }
305 /** Runs one batch (one movie frame of time); resolves when its results have been read back. */
306 async runBatch(): Promise<void> {
307 if (this.phase === "done") return;
308 const cur = this.inFlight ?? this.encodeBatch();
309 this.inFlight = null;
310 if (!cur) { this.phase = "done"; return; }
311 if (this.stepsLeft() > 0) this.inFlight = this.encodeBatch();
312 await this.collect(cur);
313 }
314 private stepsLeft(): number {
315 return this.targetSamples() * this.stepsPerSample - this.totalSteps;
316 }
317 private encodeBatch(): Pending | null {
318 const steps = Math.min(this.stepsPerBatch, Math.max(0, this.stepsLeft()));
319 if (steps === 0) return null;
320 const audioSteps: number[] = [];
321 for (let s = 0; s < steps; s++) {
322 const local = this.totalSteps + s;
323 let slot = -1, outIdx = -1;
324 if (local % this.stepsPerSample === 0) {
325 slot = (local / this.stepsPerSample) % this.histLen;
326 outIdx = audioSteps.length;
327 audioSteps.push(s);
328 }
329 this.writeStepParams(s, { flags: 1, slot, outIdx });
330 }
331 this.device.queue.writeBuffer(this.stepBuf, 0, this.stepData, 0, steps * this.stepStride);
332 const CHUNK = 1024;
333 let enc = this.device.createCommandEncoder();
334 let pass = enc.beginComputePass();
335 const [wgX, wgY] = this.stepWorkgroups();
336 const n = this.plate.grid.n;
337 let cur = this.cur;
338 pass.setPipeline(this.stepPipe);
339 for (let s = 0; s < steps; s++) {
340 if (s > 0 && s % CHUNK === 0) {
341 pass.end();
342 this.device.queue.submit([enc.finish()]);
343 enc = this.device.createCommandEncoder();
344 pass = enc.beginComputePass();
345 pass.setPipeline(this.stepPipe);
346 }
347 pass.setBindGroup(0, this.stepBG[cur], [s * this.stepStride]);
348 pass.dispatchWorkgroups(wgX, wgY);
349 cur = 1 - cur;
350 }
351 const nAudio = audioSteps.length;
352 if (nAudio > 0) {
353 pass.setPipeline(this.audioPipe);
354 pass.setBindGroup(0, this.audioBG, [audioSteps[0] * this.stepStride]);
355 pass.dispatchWorkgroups(nAudio);
356 }
357 pass.setPipeline(this.extractPipe);
358 pass.setBindGroup(0, this.extractBG[cur]);
359 pass.dispatchWorkgroups(Math.ceil(n / 256));
360 pass.end();
361 const staging = this.staging[this.stagingIdx];
362 this.stagingIdx = 1 - this.stagingIdx;
363 enc.copyBufferToBuffer(this.audioOut, 0, staging, 0, Math.max(4, nAudio * 4));
364 enc.copyBufferToBuffer(this.frameBuf, 0, staging, this.samplesPerBatch * 4, n * 4);
365 enc.copyBufferToBuffer(this.stickBuf, 0, staging, (this.samplesPerBatch + n) * 4, STICK_LEN * 4);
366 this.device.queue.submit([enc.finish()]);
367 this.cur = cur;
368 this.totalSteps += steps;
369 return { staging, steps, nAudio, stickParity: this.totalSteps & 1 };
370 }
371 private async collect(b: Pending): Promise<void> {
372 const n = this.plate.grid.n;
373 await b.staging.mapAsync(GPUMapMode.READ);
374 const view = new Float32Array(b.staging.getMappedRange());
375 const base = this.samplesPerBatch + n + 8 * b.stickParity;
376 this.stick = { x: view[base], v: view[base + 1], force: view[base + 2], maxForce: view[base + 3], withdrawn: view[base + 5] !== 0 };
377 this.appendAudio(view.slice(0, b.nAudio));
378 const frame = view.slice(this.samplesPerBatch, this.samplesPerBatch + n);
379 for (let k = 0; k < n; k++) this.wmax = Math.max(this.wmax, Math.abs(frame[k]));
380 this.frames.push(frame);
381 b.staging.unmap();
382 if (this.samples >= this.targetSamples()) this.phase = "done";
383 }
385 // ---------- Testing support ----------
2481a63Remove the standJeremy Magland 387 /** Advances the state by exactly one free step (no drumstick, no audio) with the losses. */
5778036cymbalsim: scripted plate and cymbal designer with a modal linear model and a von Karman GPU modelJeremy Magland 388 async stepOnce(): Promise<void> {
389 this.writeStepParams(0, { flags: 0, slot: -1, outIdx: -1 });
390 this.device.queue.writeBuffer(this.stepBuf, 0, this.stepData, 0, this.stepStride);
391 const enc = this.device.createCommandEncoder();
392 const pass = enc.beginComputePass();
393 pass.setPipeline(this.stepPipe);
394 pass.setBindGroup(0, this.stepBG[this.cur], [0]);
395 pass.dispatchWorkgroups(...this.stepWorkgroups());
396 pass.end();
397 this.device.queue.submit([enc.finish()]);
398 this.cur = 1 - this.cur;
399 this.totalSteps++;
400 await this.device.queue.onSubmittedWorkDone();
401 }
2481a63Remove the standJeremy Magland 402 /** Overwrites the state: pos as [u | v | w | wLo], vel as [du | dv | dw]. */
403 writeState(pos: Float32Array<ArrayBuffer>, vel: Float32Array<ArrayBuffer>): void {
5778036cymbalsim: scripted plate and cymbal designer with a modal linear model and a von Karman GPU modelJeremy Magland 404 const n = this.plate.grid.n;
405 if (pos.length !== 4 * n || vel.length !== 3 * n) throw new Error("state arrays have the wrong length");
406 this.device.queue.writeBuffer(this.pos[this.cur], 0, pos);
407 this.device.queue.writeBuffer(this.vel[this.cur], 0, vel);
408 }
2481a63Remove the standJeremy Magland 409 /** Reads back the full state as [u | v | w | wLo] and [du | dv | dw]. */
410 async readState(): Promise<{ pos: Float32Array; vel: Float32Array }> {
411 const n = this.plate.grid.n;
412 const stg = this.device.createBuffer({ size: 7 * n * 4, usage: GPUBufferUsage.MAP_READ | GPUBufferUsage.COPY_DST });
5778036cymbalsim: scripted plate and cymbal designer with a modal linear model and a von Karman GPU modelJeremy Magland 413 const enc = this.device.createCommandEncoder();
414 enc.copyBufferToBuffer(this.pos[this.cur], 0, stg, 0, 4 * n * 4);
415 enc.copyBufferToBuffer(this.vel[this.cur], 0, stg, 4 * n * 4, 3 * n * 4);
416 this.device.queue.submit([enc.finish()]);
417 await stg.mapAsync(GPUMapMode.READ);
418 const all = new Float32Array(stg.getMappedRange()).slice();
419 stg.unmap();
420 stg.destroy();
2481a63Remove the standJeremy Magland 421 return { pos: all.subarray(0, 4 * n), vel: all.subarray(4 * n) };
5778036cymbalsim: scripted plate and cymbal designer with a modal linear model and a von Karman GPU modelJeremy Magland 422 }
424 destroy(): void {
425 this.inFlight = null;
2481a63Remove the standJeremy Magland 426 for (const b of [...this.pos, ...this.vel, this.hist, this.stickBuf, this.nd, this.lw, this.audioOut, this.frameBuf, this.paramsBuf, this.stepBuf, ...this.staging]) b.destroy();
5778036cymbalsim: scripted plate and cymbal designer with a modal linear model and a von Karman GPU modelJeremy Magland 427 }
428 get fs(): number { return this.params.fs; }
429 get fps(): number { return this.params.fps; }
430 get simParams(): SimParams { return this.params; }
431}