f8a8a6cNonlinear steel sheet gong simulator with WebGPUJeremy Magland 1// WebGPU time-stepping driver for the von Karman plate.
2//
3// The simulation advances in batches of one movie frame (1/fps seconds). A batch
4// encodes one compute dispatch per time step (with a per-step uniform addressed
5// by dynamic offset), a rigid-body projection every RB_INTERVAL steps, one
6// retarded-time audio sum per audio sample, and one movie-frame extraction at
7// the end. Audio, frame, and statistics are read back after every batch.
9import {
10 type PhysParams,
11 type Grid,
12 type Material,
13 makeGrid,
14 makeMaterial,
18 nodeY,
19 trapWeight,
20 HIST_LEN,
21 RB_INTERVAL,
22 RB_INTERVAL_HELD,
23 RHO_AIR,
24 SPEED_OF_SOUND,
25} from "../params.ts";
26import { maxAngularFrequency } from "./cpuRef.ts";
27import commonWgsl from "./shaders/common.wgsl?raw";
28import stepWgsl from "./shaders/step.wgsl?raw";
29import rigidWgsl from "./shaders/rigid.wgsl?raw";
30import audioWgsl from "./shaders/audio.wgsl?raw";
31import extractWgsl from "./shaders/extract.wgsl?raw";
33export type Phase = "push" | "hold" | "recording" | "done";
35export interface SimStatus {
36 phase: Phase;
37 /** Simulation time relative to the release [s]. */
38 t: number;
39 totalSteps: number;
40 ke: number;
41 wmax: number;
42 samples: number;
43 frames: number;
44 peakPressure: number;
45 delayClamped: boolean;
10998acDrumstick strike excitation (Hertz contact with a free mass)Jeremy Magland 46 /** Drumstick state (strike excitation only). */
47 stick: { x: number; v: number; force: number; maxForce: number } | null;
50const STATS_LEN = 16;
10998acDrumstick strike excitation (Hertz contact with a free mass)Jeremy Magland 51const STICK_LEN = 8;
f8a8a6cNonlinear steel sheet gong simulator with WebGPUJeremy Magland 52
53export class Simulation {
54 readonly grid: Grid;
55 readonly mat: Material;
56 readonly dt: number;
57 /** Time steps per audio sample. */
58 readonly stepsPerSample: number;
59 readonly samplesPerBatch: number;
60 readonly stepsPerBatch: number;
61 readonly omegaMax: number;
63 /** Raw listener pressure [Pa] at fs, from the release. */
64 audio: Float32Array = new Float32Array(0);
65 samples = 0;
66 /** Movie frames of w [m]; frame k is at t = k / fps. */
67 frames: Float32Array<ArrayBuffer>[] = [];
69 phase: Phase = "push";
70 /** Steps since the start of the current phase clock (prep or recording). */
71 private stepCount = 0;
72 private prepTime = 0;
73 totalSteps = 0;
74 ke = 0;
75 wmax = 0;
76 peakPressure = 0;
77 delayClamped = false;
78 private hpState = { y: 0, x: 0 };
80 private device: GPUDevice;
81 private params: PhysParams;
82 private paramsBuf: GPUBuffer;
83 private stepBuf: GPUBuffer;
84 private stepStride: number;
85 private stepData: ArrayBuffer;
86 private pos: GPUBuffer[];
87 private vel: GPUBuffer[];
88 private hist: GPUBuffer;
adc495dFinite finger release: ramp the recorded reaction force to zero over a release timeJeremy Magland 89 private react: GPUBuffer;
8dcc11fInitial-shape imperfection (Marguerre shallow-shell strains)Jeremy Magland 90 /** Stress-free initial shape (f32 per node), also used by the renderer. */
91 readonly w0: GPUBuffer;
93 private audioOut: GPUBuffer;
94 private frameBuf: GPUBuffer;
95 private statsBuf: GPUBuffer;
96 private staging: GPUBuffer[];
97 private stagingIdx = 0;
98 /** Index of the buffer pair holding the current state. */
99 private cur = 0;
100 /** The `displayW` buffer is a plain f32 array of w for the renderer. */
101 readonly displayW: GPUBuffer;
103 private stepPipe: GPUComputePipeline;
104 private sumsPipe: GPUComputePipeline;
105 private applyPipe: GPUComputePipeline;
b9a5a7eSelf-equilibrated release force and loudness normalizationJeremy Magland 106 private projectPipe: GPUComputePipeline;
107 private projectBG: GPUBindGroup;
f8a8a6cNonlinear steel sheet gong simulator with WebGPUJeremy Magland 108 private audioPipe: GPUComputePipeline;
109 private extractPipe: GPUComputePipeline;
110 private stepBG: GPUBindGroup[];
111 private rbBG: GPUBindGroup[];
112 private audioBG: GPUBindGroup;
113 private extractBG: GPUBindGroup[];
114 private displayBG: GPUBindGroup[];
115 private stagingSize: number;
117 constructor(device: GPUDevice, params: PhysParams) {
118 this.device = device;
119 this.params = { ...params };
120 const g = makeGrid(params);
121 const m = makeMaterial(params);
122 this.grid = g;
123 this.mat = m;
125 // Stable time step, aligned so that an integer number of steps make one audio sample.
126 this.omegaMax = maxAngularFrequency(g, m);
127 const dtMax = (0.85 * 2) / this.omegaMax;
128 this.stepsPerSample = Math.ceil(1 / (params.fs * dtMax));
129 this.dt = 1 / (params.fs * this.stepsPerSample);
130 this.samplesPerBatch = Math.max(1, Math.round(params.fs / params.fps));
131 this.stepsPerBatch = this.samplesPerBatch * this.stepsPerSample;
133 const n = g.n;
134 const f32 = (size: number, usage = GPUBufferUsage.STORAGE | GPUBufferUsage.COPY_SRC | GPUBufferUsage.COPY_DST) =>
135 device.createBuffer({ size: Math.max(16, size * 4), usage });
137 this.pos = [f32(4 * n), f32(4 * n)]; // [u | v | w | wLo]
138 this.vel = [f32(3 * n), f32(3 * n)]; // [du | dv | dw]
139 this.hist = f32(HIST_LEN * n);
10998acDrumstick strike excitation (Hertz contact with a free mass)Jeremy Magland 140 this.react = f32(n + STICK_LEN);
8dcc11fInitial-shape imperfection (Marguerre shallow-shell strains)Jeremy Magland 141 this.w0 = f32(n);
142 device.queue.writeBuffer(this.w0, 0, initialShape(params, g));
144 this.audioOut = f32(this.samplesPerBatch);
145 this.frameBuf = f32(n);
146 this.displayW = f32(n);
147 this.statsBuf = f32(STATS_LEN);
10998acDrumstick strike excitation (Hertz contact with a free mass)Jeremy Magland 148 this.stagingSize = (this.samplesPerBatch + n + STATS_LEN + STICK_LEN) * 4;
f8a8a6cNonlinear steel sheet gong simulator with WebGPUJeremy Magland 149 this.staging = [0, 1].map(() =>
150 device.createBuffer({ size: this.stagingSize, usage: GPUBufferUsage.MAP_READ | GPUBufferUsage.COPY_DST }),
151 );
153 // Constant parameters.
154 let mass = 0, Ix = 0, Iy = 0;
155 for (let j = 0; j < g.ny; j++) {
156 for (let i = 0; i < g.nx; i++) {
157 const mk = m.rhoh * g.hx * g.hy * trapWeight(g, i, j);
158 mass += mk;
159 Ix += mk * nodeY(g, j) ** 2;
160 Iy += mk * nodeX(g, i) ** 2;
161 }
162 }
10998acDrumstick strike excitation (Hertz contact with a free mass)Jeremy Magland 163 // Drumstick contact patch.
164 const [sx, sy] = strikePoint(params);
165 let strikeNode = 0, best = Infinity, patchMass = 0;
166 for (let j = 0; j < g.ny; j++) {
167 for (let i = 0; i < g.nx; i++) {
168 const d2 = (nodeX(g, i) - sx) ** 2 + (nodeY(g, j) - sy) ** 2;
169 if (d2 < best) { best = d2; strikeNode = j * g.nx + i; }
170 }
171 }
172 for (let j = 0; j < g.ny; j++) {
173 for (let i = 0; i < g.nx; i++) {
174 const d2 = (nodeX(g, i) - sx) ** 2 + (nodeY(g, j) - sy) ** 2;
175 if (d2 <= params.contactRadius ** 2 || j * g.nx + i === strikeNode) patchMass += m.rhoh * g.hx * g.hy * trapWeight(g, i, j);
176 }
177 }
178 const stick = new Float32Array(STICK_LEN);
179 stick[0] = 0; stick[1] = -params.strikeVelocity; stick[4] = 0; stick[5] = -params.strikeVelocity;
180 device.queue.writeBuffer(this.react, n * 4, stick);
182 const pf = new Float32Array(28);
f8a8a6cNonlinear steel sheet gong simulator with WebGPUJeremy Magland 183 const pu = new Uint32Array(pf.buffer);
184 pu[0] = g.nx; pu[1] = g.ny; pu[2] = n; pu[3] = HIST_LEN;
185 pf[4] = g.hx; pf[5] = g.hy; pf[6] = this.dt; pf[7] = m.nu;
186 pf[8] = m.D; pf[9] = m.C; pf[10] = m.G; pf[11] = m.rhoh;
187 pf[12] = params.Lx / 2; pf[13] = -params.Ly / 2; pf[14] = params.fingerRadius ** 2; pf[15] = mass;
bfde5e4High-frequency ringing: in-plane loss, band-limited listener sum, damping defaultsJeremy Magland 188 pf[16] = Ix; pf[17] = Iy; pf[18] = Ix + Iy; pf[19] = params.fs;
10998acDrumstick strike excitation (Hertz contact with a free mass)Jeremy Magland 189 pf[20] = sx; pf[21] = sy; pf[22] = params.contactRadius ** 2; pf[23] = patchMass;
190 pu[24] = strikeNode; pf[25] = params.contactStiffness; pf[26] = params.contactExponent; pf[27] = params.stickMass;
f8a8a6cNonlinear steel sheet gong simulator with WebGPUJeremy Magland 191 this.paramsBuf = device.createBuffer({ size: pf.byteLength, usage: GPUBufferUsage.UNIFORM | GPUBufferUsage.COPY_DST });
192 device.queue.writeBuffer(this.paramsBuf, 0, pf);
adc495dFinite finger release: ramp the recorded reaction force to zero over a release timeJeremy Magland 194 // Per-step parameters, one 48-byte struct per step at a dynamic offset.
195 this.stepStride = Math.max(48, device.limits.minUniformBufferOffsetAlignment);
f8a8a6cNonlinear steel sheet gong simulator with WebGPUJeremy Magland 196 this.stepData = new ArrayBuffer(this.stepStride * this.stepsPerBatch);
197 this.stepBuf = device.createBuffer({ size: this.stepData.byteLength, usage: GPUBufferUsage.UNIFORM | GPUBufferUsage.COPY_DST });
199 const mk = (code: string, entry: string, layout: GPUPipelineLayout | "auto" = "auto") =>
200 device.createComputePipeline({
201 layout,
202 compute: { module: device.createShaderModule({ code: commonWgsl + code }), entryPoint: entry },
203 });
204 const entry = (binding: number, type: GPUBufferBindingType, hasDynamicOffset = false): GPUBindGroupLayoutEntry => ({
205 binding,
206 visibility: GPUShaderStage.COMPUTE,
207 buffer: { type, hasDynamicOffset },
208 });
209 const stepLayout = device.createBindGroupLayout({
210 entries: [
211 entry(0, "uniform"),
212 entry(1, "uniform", true),
213 entry(2, "read-only-storage"),
214 entry(3, "read-only-storage"),
215 entry(4, "storage"),
216 entry(5, "storage"),
217 entry(6, "storage"),
adc495dFinite finger release: ramp the recorded reaction force to zero over a release timeJeremy Magland 218 entry(7, "storage"),
8dcc11fInitial-shape imperfection (Marguerre shallow-shell strains)Jeremy Magland 219 entry(8, "read-only-storage"),
221 });
222 this.stepPipe = mk(stepWgsl, "main", device.createPipelineLayout({ bindGroupLayouts: [stepLayout] }));
223 const audioLayout = device.createBindGroupLayout({
224 entries: [entry(0, "uniform"), entry(1, "uniform", true), entry(2, "read-only-storage"), entry(3, "read-only-storage"), entry(4, "storage")],
225 });
226 const rigidModule = device.createShaderModule({ code: commonWgsl + rigidWgsl });
227 const rbLayout = device.createBindGroupLayout({
228 entries: [
229 { binding: 0, visibility: GPUShaderStage.COMPUTE, buffer: { type: "uniform" } },
230 { binding: 1, visibility: GPUShaderStage.COMPUTE, buffer: { type: "storage" } },
231 { binding: 2, visibility: GPUShaderStage.COMPUTE, buffer: { type: "storage" } },
232 { binding: 3, visibility: GPUShaderStage.COMPUTE, buffer: { type: "storage" } },
233 ],
234 });
235 const rbPipeLayout = device.createPipelineLayout({ bindGroupLayouts: [rbLayout] });
236 this.sumsPipe = device.createComputePipeline({ layout: rbPipeLayout, compute: { module: rigidModule, entryPoint: "sums" } });
237 this.applyPipe = device.createComputePipeline({ layout: rbPipeLayout, compute: { module: rigidModule, entryPoint: "apply" } });
b9a5a7eSelf-equilibrated release force and loudness normalizationJeremy Magland 238 const projectLayout = device.createBindGroupLayout({
239 entries: [
240 { binding: 0, visibility: GPUShaderStage.COMPUTE, buffer: { type: "uniform" } },
241 { binding: 4, visibility: GPUShaderStage.COMPUTE, buffer: { type: "storage" } },
242 ],
243 });
244 this.projectPipe = device.createComputePipeline({
245 layout: device.createPipelineLayout({ bindGroupLayouts: [projectLayout] }),
246 compute: { module: rigidModule, entryPoint: "projectField" },
247 });
f8a8a6cNonlinear steel sheet gong simulator with WebGPUJeremy Magland 248 this.audioPipe = mk(audioWgsl, "main", device.createPipelineLayout({ bindGroupLayouts: [audioLayout] }));
249 this.extractPipe = mk(extractWgsl, "main");
adc495dFinite finger release: ramp the recorded reaction force to zero over a release timeJeremy Magland 251 const sp = { buffer: this.stepBuf, size: 48 };
f8a8a6cNonlinear steel sheet gong simulator with WebGPUJeremy Magland 252 this.stepBG = [0, 1].map((p) =>
253 device.createBindGroup({
254 layout: stepLayout,
255 entries: [
256 { binding: 0, resource: { buffer: this.paramsBuf } },
257 { binding: 1, resource: sp },
258 { binding: 2, resource: { buffer: this.pos[p] } },
259 { binding: 3, resource: { buffer: this.vel[p] } },
260 { binding: 4, resource: { buffer: this.pos[1 - p] } },
261 { binding: 5, resource: { buffer: this.vel[1 - p] } },
262 { binding: 6, resource: { buffer: this.hist } },
adc495dFinite finger release: ramp the recorded reaction force to zero over a release timeJeremy Magland 263 { binding: 7, resource: { buffer: this.react } },
8dcc11fInitial-shape imperfection (Marguerre shallow-shell strains)Jeremy Magland 264 { binding: 8, resource: { buffer: this.w0 } },
266 }),
267 );
268 this.rbBG = [0, 1].map((p) =>
269 device.createBindGroup({
270 layout: rbLayout,
271 entries: [
272 { binding: 0, resource: { buffer: this.paramsBuf } },
273 { binding: 1, resource: { buffer: this.pos[p] } },
274 { binding: 2, resource: { buffer: this.vel[p] } },
275 { binding: 3, resource: { buffer: this.statsBuf } },
276 ],
277 }),
278 );
b9a5a7eSelf-equilibrated release force and loudness normalizationJeremy Magland 279 this.projectBG = device.createBindGroup({
280 layout: projectLayout,
281 entries: [
282 { binding: 0, resource: { buffer: this.paramsBuf } },
283 { binding: 4, resource: { buffer: this.react } },
284 ],
285 });
f8a8a6cNonlinear steel sheet gong simulator with WebGPUJeremy Magland 286 this.audioBG = device.createBindGroup({
287 layout: audioLayout,
288 entries: [
289 { binding: 0, resource: { buffer: this.paramsBuf } },
290 { binding: 1, resource: sp },
291 { binding: 2, resource: { buffer: this.hist } },
292 { binding: 3, resource: { buffer: this.lw } },
293 { binding: 4, resource: { buffer: this.audioOut } },
294 ],
295 });
296 const extractBG = (p: number, out: GPUBuffer) =>
297 device.createBindGroup({
298 layout: this.extractPipe.getBindGroupLayout(0),
299 entries: [
300 { binding: 0, resource: { buffer: this.paramsBuf } },
301 { binding: 1, resource: { buffer: this.pos[p] } },
302 { binding: 2, resource: { buffer: out } },
303 ],
304 });
305 this.extractBG = [0, 1].map((p) => extractBG(p, this.frameBuf));
306 this.displayBG = [0, 1].map((p) => extractBG(p, this.displayW));
308 this.setListener(params.listenerDist, params.listenerAzimuth, params.listenerElevation, params.radiation);
10998acDrumstick strike excitation (Hertz contact with a free mass)Jeremy Magland 309 if (params.excitation === "strike") {
310 // No push phase: the sheet is at rest in its natural shape when the stick arrives.
311 this.phase = "recording";
176c8fbStop button: end a run and unlock the parameter fieldsJeremy Magland 312 this.released = true;
10998acDrumstick strike excitation (Hertz contact with a free mass)Jeremy Magland 313 this.frames.push(new Float32Array(n));
314 }
317 /** Recomputes the per-node weights and delays of the listener sum. Takes effect on later samples. */
318 setListener(dist: number, azimuthDeg: number, elevationDeg: number, radiation: "rayleigh" | "dipole"): void {
319 const g = this.grid;
320 const az = (azimuthDeg * Math.PI) / 180;
321 const el = (elevationDeg * Math.PI) / 180;
322 const L = [dist * Math.sin(az) * Math.cos(el), dist * Math.sin(el), dist * Math.cos(az) * Math.cos(el)];
323 const data = new Float32Array(2 * g.n);
324 let dmin = Infinity;
325 for (let j = 0; j < g.ny; j++) {
326 for (let i = 0; i < g.nx; i++) {
327 const k = j * g.nx + i;
328 const rx = L[0] - nodeX(g, i), ry = L[1] - nodeY(g, j), rz = L[2];
329 const r = Math.sqrt(rx * rx + ry * ry + rz * rz);
330 const dA = g.hx * g.hy * trapWeight(g, i, j);
331 let wgt = (RHO_AIR * dA) / (2 * Math.PI * r);
332 if (radiation === "dipole") wgt *= Math.abs(rz) / r;
333 data[2 * k] = wgt;
334 data[2 * k + 1] = (r / SPEED_OF_SOUND) * this.params.fs;
335 dmin = Math.min(dmin, data[2 * k + 1]);
336 }
337 }
338 this.delayClamped = false;
339 for (let k = 0; k < g.n; k++) {
340 let d = data[2 * k + 1] - dmin;
bfde5e4High-frequency ringing: in-plane loss, band-limited listener sum, damping defaultsJeremy Magland 341 if (d > HIST_LEN - 3) {
342 d = HIST_LEN - 3;
344 }
345 data[2 * k + 1] = d;
346 }
347 this.device.queue.writeBuffer(this.lw, 0, data);
348 this.params.listenerDist = dist;
349 this.params.listenerAzimuth = azimuthDeg;
350 this.params.listenerElevation = elevationDeg;
351 this.params.radiation = radiation;
352 }
bfde5e4High-frequency ringing: in-plane loss, band-limited listener sum, damping defaultsJeremy Magland 354 setDamping(sigma0: number, sigma1: number, sigmaIn: number): void {
f8a8a6cNonlinear steel sheet gong simulator with WebGPUJeremy Magland 355 this.params.sigma0 = sigma0;
356 this.params.sigma1 = sigma1;
bfde5e4High-frequency ringing: in-plane loss, band-limited listener sum, damping defaultsJeremy Magland 357 this.params.sigmaIn = sigmaIn;
360 setDuration(duration: number): void {
361 this.params.duration = duration;
176c8fbStop button: end a run and unlock the parameter fieldsJeremy Magland 362 if (this.phase === "done" && this.released && this.samples < this.targetSamples()) this.phase = "recording";
176c8fbStop button: end a run and unlock the parameter fieldsJeremy Magland 365 /** True when a finished or stopped run can be extended. */
366 canContinue(): boolean {
367 return this.phase === "done" && this.released && this.samples < this.targetSamples();
368 }
370 /** Ends the run; the state, audio, and movie are kept. */
371 stop(): void {
372 this.phase = "done";
373 }
374 private released = false;
377 return Math.round(this.params.duration * this.params.fs);
378 }
380 status(): SimStatus {
381 const t = this.phase === "push" || this.phase === "hold" ? this.prepTime - this.prepDuration() : this.stepCount * this.dt;
382 return {
383 phase: this.phase,
384 t,
385 totalSteps: this.totalSteps,
386 ke: this.ke,
387 wmax: this.wmax,
388 samples: this.samples,
389 frames: this.frames.length,
390 peakPressure: this.peakPressure,
391 delayClamped: this.delayClamped,
10998acDrumstick strike excitation (Hertz contact with a free mass)Jeremy Magland 392 stick: this.stick,
394 }
396 private prepDuration(): number {
397 return this.params.rampTime + this.params.holdMax;
398 }
400 /** Enqueues an extraction of the current w field into `displayW` (for the live view). */
401 encodeDisplay(encoder: GPUCommandEncoder): void {
402 const pass = encoder.beginComputePass();
403 pass.setPipeline(this.extractPipe);
404 pass.setBindGroup(0, this.displayBG[this.cur]);
405 pass.dispatchWorkgroups(Math.ceil(this.grid.n / 256));
406 pass.end();
407 }
409 /** Listener audio with the radiation post-processing applied [Pa]. */
410 processedAudio(): Float32Array {
411 return this.processed.subarray(0, this.samples);
412 }
413 private processed: Float32Array = new Float32Array(0);
415 private appendAudio(raw: Float32Array): void {
416 const need = this.samples + raw.length;
417 if (this.audio.length < need) {
418 const grow = Math.max(need, 2 * this.audio.length, this.params.fs);
419 const a = new Float32Array(grow);
420 a.set(this.audio.subarray(0, this.samples));
421 this.audio = a;
422 const b = new Float32Array(grow);
423 b.set(this.processed.subarray(0, this.samples));
424 this.processed = b;
425 }
426 this.audio.set(raw, this.samples);
427 // Dipole option: first-order high-pass at the baffle cutoff c / (2 sqrt(Lx Ly)).
428 const dipole = this.params.radiation === "dipole";
429 const fc = SPEED_OF_SOUND / (2 * Math.sqrt(this.params.Lx * this.params.Ly));
430 const alpha = 1 / (1 + (2 * Math.PI * fc) / this.params.fs);
431 for (let k = 0; k < raw.length; k++) {
432 const x = raw[k];
433 let y = x;
434 if (dipole) {
435 y = alpha * (this.hpState.y + x - this.hpState.x);
436 }
437 this.hpState = { y, x };
438 this.processed[this.samples + k] = y;
439 this.peakPressure = Math.max(this.peakPressure, Math.abs(y));
440 }
441 this.samples = need;
442 }
444 /** Runs one batch (one movie frame of time). Resolves when its results have been read back. */
445 async runBatch(): Promise<void> {
446 if (this.phase === "done") return;
447 const p = this.params;
448 const recording = this.phase === "recording";
449 const nSteps = this.stepsPerBatch;
450 const stepsLeft = recording ? (this.targetSamples() - this.samples) * this.stepsPerSample : nSteps;
451 const steps = Math.min(nSteps, Math.max(0, stepsLeft));
452 if (steps === 0) {
453 this.phase = "done";
454 return;
455 }
457 // Fill the per-step parameters.
458 const u32 = new Uint32Array(this.stepData);
459 const i32 = new Int32Array(this.stepData);
460 const f32 = new Float32Array(this.stepData);
461 const stride = this.stepStride / 4;
462 const audioStepsInBatch: number[] = [];
463 for (let s = 0; s < steps; s++) {
464 const o = s * stride;
465 const local = this.stepCount + s;
466 u32[o] = this.totalSteps + s;
adc495dFinite finger release: ramp the recorded reaction force to zero over a release timeJeremy Magland 467 let flags = 0, pushW = 0, pushV = 0, sigma0 = p.sigma0, slot = -1, outIdx = -1, releaseGain = 0;
10998acDrumstick strike excitation (Hertz contact with a free mass)Jeremy Magland 468 const parity = (this.totalSteps + s) & 1;
470 const t = this.prepTime + (s + 1) * this.dt; // state after this step
471 const tau = Math.min(1, t / p.rampTime);
472 const smooth = tau * tau * (3 - 2 * tau);
473 const dsmooth = tau < 1 ? (6 * tau * (1 - tau)) / p.rampTime : 0;
adc495dFinite finger release: ramp the recorded reaction force to zero over a release timeJeremy Magland 474 flags = 1 | (tau >= 1 ? 2 : 0); // record the reaction once the ramp is complete
f8a8a6cNonlinear steel sheet gong simulator with WebGPUJeremy Magland 475 pushW = p.pushAmount * smooth;
476 pushV = p.pushAmount * dsmooth;
477 sigma0 = p.settleSigma0;
adc495dFinite finger release: ramp the recorded reaction force to zero over a release timeJeremy Magland 478 } else {
479 const t = local * this.dt;
10998acDrumstick strike excitation (Hertz contact with a free mass)Jeremy Magland 480 if (p.excitation === "strike") {
481 flags = 8;
482 } else if (p.releaseTime > 0 && t < p.releaseTime) {
adc495dFinite finger release: ramp the recorded reaction force to zero over a release timeJeremy Magland 483 const tau = t / p.releaseTime;
484 releaseGain = 1 - tau * tau * (3 - 2 * tau);
485 flags = 4;
486 }
487 if (local % this.stepsPerSample === 0) {
488 const sampleIdx = local / this.stepsPerSample;
489 slot = sampleIdx % HIST_LEN;
490 outIdx = audioStepsInBatch.length;
491 audioStepsInBatch.push(s);
492 }
494 u32[o + 1] = flags;
495 i32[o + 2] = slot;
496 i32[o + 3] = outIdx;
497 f32[o + 4] = pushW;
498 f32[o + 5] = pushV;
499 f32[o + 6] = sigma0;
500 f32[o + 7] = p.sigma1;
adc495dFinite finger release: ramp the recorded reaction force to zero over a release timeJeremy Magland 501 f32[o + 8] = releaseGain;
bfde5e4High-frequency ringing: in-plane loss, band-limited listener sum, damping defaultsJeremy Magland 502 f32[o + 9] = p.sigmaIn;
10998acDrumstick strike excitation (Hertz contact with a free mass)Jeremy Magland 503 u32[o + 10] = parity;
505 this.device.queue.writeBuffer(this.stepBuf, 0, this.stepData, 0, steps * this.stepStride);
507 // Encode the batch in chunks so that the GPU starts while JS is still encoding.
508 const CHUNK = 1024;
509 let enc = this.device.createCommandEncoder();
510 let pass = enc.beginComputePass();
511 const wgX = Math.ceil(this.grid.nx / 16), wgY = Math.ceil(this.grid.ny / 16);
512 const wgN = Math.ceil(this.grid.n / 256);
513 let cur = this.cur;
514 let audioPtr = 0;
515 const rbInterval = recording ? RB_INTERVAL : RB_INTERVAL_HELD;
516 let stepBound = false;
517 for (let s = 0; s < steps; s++) {
518 if (s > 0 && s % CHUNK === 0) {
519 pass.end();
520 this.device.queue.submit([enc.finish()]);
521 enc = this.device.createCommandEncoder();
522 pass = enc.beginComputePass();
523 stepBound = false;
524 }
525 if (!stepBound) {
526 pass.setPipeline(this.stepPipe);
527 stepBound = true;
528 }
529 pass.setBindGroup(0, this.stepBG[cur], [s * this.stepStride]);
530 pass.dispatchWorkgroups(wgX, wgY);
531 cur = 1 - cur;
532 const global = this.totalSteps + s + 1;
533 if (global % rbInterval === 0 || s === steps - 1) {
534 pass.setPipeline(this.sumsPipe);
535 pass.setBindGroup(0, this.rbBG[cur]);
536 pass.dispatchWorkgroups(1);
537 pass.setPipeline(this.applyPipe);
538 pass.setBindGroup(0, this.rbBG[cur]);
539 pass.dispatchWorkgroups(wgN);
540 stepBound = false;
541 }
542 if (audioPtr < audioStepsInBatch.length && audioStepsInBatch[audioPtr] === s) {
543 pass.setPipeline(this.audioPipe);
544 pass.setBindGroup(0, this.audioBG, [s * this.stepStride]);
545 pass.dispatchWorkgroups(1);
546 audioPtr++;
547 stepBound = false;
548 }
549 }
550 pass.setPipeline(this.extractPipe);
551 pass.setBindGroup(0, this.extractBG[cur]);
552 pass.dispatchWorkgroups(wgN);
553 pass.end();
554 const staging = this.staging[this.stagingIdx];
555 this.stagingIdx = 1 - this.stagingIdx;
556 const nAudio = audioStepsInBatch.length;
557 enc.copyBufferToBuffer(this.audioOut, 0, staging, 0, Math.max(4, nAudio * 4));
558 enc.copyBufferToBuffer(this.frameBuf, 0, staging, this.samplesPerBatch * 4, this.grid.n * 4);
559 enc.copyBufferToBuffer(this.statsBuf, 0, staging, (this.samplesPerBatch + this.grid.n) * 4, STATS_LEN * 4);
10998acDrumstick strike excitation (Hertz contact with a free mass)Jeremy Magland 560 enc.copyBufferToBuffer(this.react, this.grid.n * 4, staging, (this.samplesPerBatch + this.grid.n + STATS_LEN) * 4, STICK_LEN * 4);
f8a8a6cNonlinear steel sheet gong simulator with WebGPUJeremy Magland 561 this.device.queue.submit([enc.finish()]);
562 this.cur = cur;
563 this.totalSteps += steps;
564 this.stepCount += steps;
565 if (!recording) this.prepTime += steps * this.dt;
567 await staging.mapAsync(GPUMapMode.READ);
568 const view = new Float32Array(staging.getMappedRange());
569 const stats = view.subarray(this.samplesPerBatch + this.grid.n, this.samplesPerBatch + this.grid.n + STATS_LEN);
570 this.ke = stats[12];
571 this.wmax = stats[13];
10998acDrumstick strike excitation (Hertz contact with a free mass)Jeremy Magland 572 if (p.excitation === "strike") {
573 // The slot written by the last step is the one the next step reads.
574 const base = this.samplesPerBatch + this.grid.n + STATS_LEN + 4 * (this.totalSteps & 1);
575 this.stick = { x: view[base], v: view[base + 1], force: view[base + 2], maxForce: view[base + 3] };
576 }
578 this.appendAudio(view.slice(0, nAudio));
579 this.frames.push(view.slice(this.samplesPerBatch, this.samplesPerBatch + this.grid.n));
580 } else {
581 this.lastPrepFrame = view.slice(this.samplesPerBatch, this.samplesPerBatch + this.grid.n);
582 }
583 staging.unmap();
585 if (!recording) {
586 if (this.phase === "push" && this.prepTime >= p.rampTime) this.phase = "hold";
587 const settled = this.phase === "hold" && this.ke < p.settleKE;
588 if (settled || this.prepTime >= this.prepDuration()) this.release();
589 } else if (this.samples >= this.targetSamples()) {
590 this.phase = "done";
591 }
592 }
593 private lastPrepFrame: Float32Array<ArrayBuffer> | null = null;
10998acDrumstick strike excitation (Hertz contact with a free mass)Jeremy Magland 594 private stick: SimStatus["stick"] = null;
f8a8a6cNonlinear steel sheet gong simulator with WebGPUJeremy Magland 595
596 private release(): void {
597 this.phase = "recording";
176c8fbStop button: end a run and unlock the parameter fieldsJeremy Magland 598 this.released = true;
600 this.totalSteps = 0;
b9a5a7eSelf-equilibrated release force and loudness normalizationJeremy Magland 601 // The finger force and the distributed holding reaction relax together:
602 // remove the rigid-body components of the recorded reaction.
603 const enc = this.device.createCommandEncoder();
604 const pass = enc.beginComputePass();
605 pass.setPipeline(this.projectPipe);
606 pass.setBindGroup(0, this.projectBG);
607 pass.dispatchWorkgroups(1);
608 pass.end();
609 this.device.queue.submit([enc.finish()]);
f8a8a6cNonlinear steel sheet gong simulator with WebGPUJeremy Magland 610 // Frame 0 is the held, settled shape.
611 if (this.lastPrepFrame) this.frames.push(this.lastPrepFrame);
612 }
614 /** Reads back the full state (debugging and tests). */
615 async readState(): Promise<{ pos: Float32Array; vel: Float32Array }> {
616 const n = this.grid.n;
617 const stg = this.device.createBuffer({ size: 7 * n * 4, usage: GPUBufferUsage.MAP_READ | GPUBufferUsage.COPY_DST });
618 const enc = this.device.createCommandEncoder();
619 enc.copyBufferToBuffer(this.pos[this.cur], 0, stg, 0, 4 * n * 4);
620 enc.copyBufferToBuffer(this.vel[this.cur], 0, stg, 4 * n * 4, 3 * n * 4);
621 this.device.queue.submit([enc.finish()]);
622 await stg.mapAsync(GPUMapMode.READ);
623 const all = new Float32Array(stg.getMappedRange()).slice();
624 stg.unmap();
625 stg.destroy();
626 // Returned as [u | v | w | wLo] and [du | dv | dw].
627 return { pos: all.subarray(0, 4 * n), vel: all.subarray(4 * n) };
628 }
630 destroy(): void {
8dcc11fInitial-shape imperfection (Marguerre shallow-shell strains)Jeremy Magland 631 for (const b of [...this.pos, ...this.vel, this.hist, this.react, this.w0, this.lw, this.audioOut, this.frameBuf, this.displayW, this.statsBuf, this.paramsBuf, this.stepBuf, ...this.staging]) b.destroy();
634 get fs(): number {
635 return this.params.fs;
636 }
637 get fps(): number {
638 return this.params.fps;
639 }
640 get physParams(): PhysParams {
641 return this.params;
642 }
643}