d8f9e54Clean rebuild of the steel sheet simulatorJeremy 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
4// batch encodes one compute dispatch per time step (with a per-step uniform
5// addressed by dynamic offset), a rigid-body projection every RB_INTERVAL
a2a1b68Drumstick by default, hit sequences, 170 GPa steelJeremy Magland 6// steps, one retarded-time listener sum over the batch's audio samples, and
7// one movie-frame extraction at the end. Audio, frame, and statistics are
8// read back after every batch. One batch is kept encoded ahead of the one
9// being read back, so the GPU does not wait for the readback or for the
10// encoding; a decision that depends on a readback (the release of the held
11// sheet) therefore takes effect one batch later.
d8f9e54Clean rebuild of the steel sheet simulatorJeremy Magland 12
13import {
14 type PhysParams,
15 type Grid,
16 type Material,
17 makeGrid,
18 makeMaterial,
19 initialShape,
20 strikePoint,
21 pushPoint,
22 contactPatch,
24 type Hit,
26 nodeY,
27 nodeWeight,
28 nodeWeights,
30 HIST_MARGIN,
32 RB_INTERVAL_HELD,
33 RHO_AIR,
35 SETTLE_WINDOW,
37} from "../params.ts";
38import { maxAngularFrequency } from "./cpuRef.ts";
39import commonWgsl from "./shaders/common.wgsl?raw";
40import stepWgsl from "./shaders/step.wgsl?raw";
41import rigidWgsl from "./shaders/rigid.wgsl?raw";
42import audioWgsl from "./shaders/audio.wgsl?raw";
43import extractWgsl from "./shaders/extract.wgsl?raw";
45export type Phase = "push" | "hold" | "recording" | "done";
47export interface StickState {
48 x: number;
49 v: number;
50 force: number;
51 maxForce: number;
52}
54export interface SimStatus {
55 phase: Phase;
56 /** Simulation time relative to the release or the impact [s]. */
57 t: number;
58 totalSteps: number;
59 ke: number;
60 wmax: number;
61 samples: number;
62 frames: number;
63 peakPressure: number;
64 delayClamped: boolean;
65 /** Drumstick state (strike excitation only). */
66 stick: StickState | null;
a2a1b68Drumstick by default, hit sequences, 170 GPa steelJeremy Magland 67 /** Hits delivered so far, of the scheduled total (strike excitation only). */
68 hits: { done: number; total: number } | null;
282f706Single excitation panel; rigid-body projection optional while ringingJeremy Magland 71const PARAMS_LEN = 36; // floats in the Params uniform
d8f9e54Clean rebuild of the steel sheet simulatorJeremy Magland 72const STEP_LEN = 12; // floats in the StepParams uniform
73const STATS_LEN = 16;
74const STICK_LEN = 8;
a2a1b68Drumstick by default, hit sequences, 170 GPa steelJeremy Magland 75/** Nodes per workgroup along each axis in step.wgsl. */
76const TILE = 12;
78/** A batch that has been encoded and submitted but not read back. */
79interface Pending {
80 staging: GPUBuffer;
81 steps: number;
82 nAudio: number;
83 recording: boolean;
84 /** Preparation clock at the end of the batch [s]. */
85 prepTimeEnd: number;
86 /** Which of the two stick-state slots the batch's last step wrote. */
87 stickParity: number;
88}
d8f9e54Clean rebuild of the steel sheet simulatorJeremy Magland 89
90export class Simulation {
91 readonly grid: Grid;
92 readonly mat: Material;
93 readonly dt: number;
94 /** Time steps per audio sample. */
95 readonly stepsPerSample: number;
96 readonly samplesPerBatch: number;
97 readonly stepsPerBatch: number;
98 readonly omegaMax: number;
a2a1b68Drumstick by default, hit sequences, 170 GPa steelJeremy Magland 99 /** Ring buffer length of the listener sum [samples]. */
100 readonly histLen: number;
d8f9e54Clean rebuild of the steel sheet simulatorJeremy Magland 101
102 /** Listener pressure [Pa] at fs, from the release or the impact. */
103 audio: Float32Array = new Float32Array(0);
104 samples = 0;
105 /** Movie frames of w [m]; frame k is at t = k / fps. */
106 frames: Float32Array<ArrayBuffer>[] = [];
108 phase: Phase = "push";
109 totalSteps = 0;
110 ke = 0;
111 wmax = 0;
112 peakPressure = 0;
113 delayClamped = false;
a2a1b68Drumstick by default, hit sequences, 170 GPa steelJeremy Magland 115 /** Steps enqueued since the start of the current clock (preparation or recording). */
a2a1b68Drumstick by default, hit sequences, 170 GPa steelJeremy Magland 117 /** Preparation time enqueued [s]. */
119 private released = false;
120 private lastPrepFrame: Float32Array<ArrayBuffer> | null = null;
121 private stick: StickState | null = null;
a2a1b68Drumstick by default, hit sequences, 170 GPa steelJeremy Magland 122 /** Drumstick hit schedule (strike excitation only) and the index of the next hit to arm. */
123 readonly hits: Hit[];
124 private hitSteps: number[];
125 private nextHit = 1;
126 private inFlight: Pending | null = null;
127 /** The held sheet has settled; release once the batch in flight has been read back. */
128 private releaseRequested = false;
129 /** Kinetic energy after each batch of the hold, for the plateau criterion. */
130 private keHistory: { t: number; ke: number }[] = [];
d8f9e54Clean rebuild of the steel sheet simulatorJeremy Magland 131
132 private device: GPUDevice;
133 private params: PhysParams;
134 private paramsBuf: GPUBuffer;
135 private stepBuf: GPUBuffer;
136 private stepStride: number;
137 private stepData: ArrayBuffer;
138 private pos: GPUBuffer[];
139 private vel: GPUBuffer[];
140 private hist: GPUBuffer;
141 private react: GPUBuffer;
142 /** Stress-free initial shape (f32 per node), also used by the renderer. */
143 readonly w0: GPUBuffer;
144 /** Node weights (f32 per node, 0 outside the sheet), also used by the renderer. */
145 readonly nw: GPUBuffer;
146 private lw: GPUBuffer;
147 private audioOut: GPUBuffer;
148 private frameBuf: GPUBuffer;
149 private statsBuf: GPUBuffer;
150 private staging: GPUBuffer[];
151 private stagingIdx = 0;
152 private stagingSize: number;
153 /** Index of the buffer pair holding the current state. */
154 private cur = 0;
155 /** Plain f32 array of w for the live view. */
156 readonly displayW: GPUBuffer;
158 private stepPipe: GPUComputePipeline;
159 private sumsPipe: GPUComputePipeline;
160 private applyPipe: GPUComputePipeline;
161 private projectPipe: GPUComputePipeline;
162 private audioPipe: GPUComputePipeline;
163 private extractPipe: GPUComputePipeline;
164 private stepBG: GPUBindGroup[];
165 private rbBG: GPUBindGroup[];
166 private projectBG: GPUBindGroup;
167 private audioBG: GPUBindGroup;
168 private extractBG: GPUBindGroup[];
169 private displayBG: GPUBindGroup[];
171 constructor(device: GPUDevice, params: PhysParams) {
172 this.device = device;
173 this.params = { ...params };
174 const g = makeGrid(params);
175 const m = makeMaterial(params);
176 this.grid = g;
177 this.mat = m;
179 // Stable time step, rounded so that an integer number of steps make one audio sample.
a2a1b68Drumstick by default, hit sequences, 170 GPa steelJeremy Magland 180 this.omegaMax = maxAngularFrequency(g, m, params.inPlaneMass);
181 const dtMax = (CFL * 2) / this.omegaMax;
d8f9e54Clean rebuild of the steel sheet simulatorJeremy Magland 182 this.stepsPerSample = Math.ceil(1 / (params.fs * dtMax));
183 this.dt = 1 / (params.fs * this.stepsPerSample);
184 this.samplesPerBatch = Math.max(1, Math.round(params.fs / params.fps));
185 this.stepsPerBatch = this.samplesPerBatch * this.stepsPerSample;
a2a1b68Drumstick by default, hit sequences, 170 GPa steelJeremy Magland 186 this.histLen = this.samplesPerBatch + HIST_MARGIN;
d8f9e54Clean rebuild of the steel sheet simulatorJeremy Magland 187
188 const n = g.n;
189 const f32 = (size: number, usage = GPUBufferUsage.STORAGE | GPUBufferUsage.COPY_SRC | GPUBufferUsage.COPY_DST) =>
190 device.createBuffer({ size: Math.max(16, size * 4), usage });
192 this.pos = [f32(4 * n), f32(4 * n)]; // [u | v | w | wLo]
193 this.vel = [f32(3 * n), f32(3 * n)]; // [du | dv | dw]
a2a1b68Drumstick by default, hit sequences, 170 GPa steelJeremy Magland 194 this.hist = f32(this.histLen * n);
d8f9e54Clean rebuild of the steel sheet simulatorJeremy Magland 195 this.react = f32(n + STICK_LEN);
196 this.w0 = f32(n);
197 device.queue.writeBuffer(this.w0, 0, initialShape(params, g));
198 this.nw = f32(n);
199 device.queue.writeBuffer(this.nw, 0, nodeWeights(g));
200 this.lw = f32(2 * n);
201 this.audioOut = f32(this.samplesPerBatch);
202 this.frameBuf = f32(n);
203 this.displayW = f32(n);
204 this.statsBuf = f32(STATS_LEN);
205 this.stagingSize = (this.samplesPerBatch + n + STATS_LEN + STICK_LEN) * 4;
206 this.staging = [0, 1].map(() =>
207 device.createBuffer({ size: this.stagingSize, usage: GPUBufferUsage.MAP_READ | GPUBufferUsage.COPY_DST }),
208 );
210 // Mass and moments of inertia of the sheet.
211 let mass = 0, Ix = 0, Iy = 0;
212 for (let j = 0; j < g.ny; j++) {
213 for (let i = 0; i < g.nx; i++) {
214 const mk = m.rhoh * g.hx * g.hy * nodeWeight(g, i, j);
215 mass += mk;
216 Ix += mk * nodeY(g, j) ** 2;
217 Iy += mk * nodeX(g, i) ** 2;
218 }
219 }
a2a1b68Drumstick by default, hit sequences, 170 GPa steelJeremy Magland 220 // Drumstick contact patch, hit schedule (in steps, strictly increasing),
221 // and initial stick state (position 0, moving toward the sheet at the speed of the first hit).
d8f9e54Clean rebuild of the steel sheet simulatorJeremy Magland 222 const cp = contactPatch(params, g);
a2a1b68Drumstick by default, hit sequences, 170 GPa steelJeremy Magland 223 this.hits = hitSchedule(params);
224 this.hitSteps = [];
225 for (const h of this.hits) this.hitSteps.push(Math.max(this.hitSteps.length === 0 ? 0 : this.hitSteps[this.hitSteps.length - 1] + 1, Math.round(h.t / this.dt)));
d8f9e54Clean rebuild of the steel sheet simulatorJeremy Magland 226 const stick = new Float32Array(STICK_LEN);
a2a1b68Drumstick by default, hit sequences, 170 GPa steelJeremy Magland 227 stick[1] = -this.hits[0].v; stick[5] = -this.hits[0].v;
d8f9e54Clean rebuild of the steel sheet simulatorJeremy Magland 228 device.queue.writeBuffer(this.react, n * 4, stick);
230 const [px, py] = pushPoint(params);
231 const [sx, sy] = strikePoint(params);
232 const pf = new Float32Array(PARAMS_LEN);
233 const pu = new Uint32Array(pf.buffer);
a2a1b68Drumstick by default, hit sequences, 170 GPa steelJeremy Magland 234 pu[0] = g.nx; pu[1] = g.ny; pu[2] = n; pu[3] = this.histLen;
d8f9e54Clean rebuild of the steel sheet simulatorJeremy Magland 235 pf[4] = g.hx; pf[5] = g.hy; pf[6] = this.dt; pf[7] = m.nu;
236 pf[8] = m.D; pf[9] = m.C; pf[10] = m.G; pf[11] = m.rhoh;
237 pf[12] = px; pf[13] = py; pf[14] = params.fingerRadius ** 2; pf[15] = mass;
238 pf[16] = Ix; pf[17] = Iy; pf[18] = Ix + Iy; pf[19] = params.fs;
239 pf[20] = sx; pf[21] = sy; pf[22] = params.contactRadius ** 2; pf[23] = cp.mass;
240 pu[24] = cp.nearest; pf[25] = params.contactStiffness; pf[26] = params.contactExponent; pf[27] = params.stickMass;
241 pu[28] = cp.box[0]; pu[29] = cp.box[1]; pu[30] = cp.box[2]; pu[31] = cp.box[3];
282f706Single excitation panel; rigid-body projection optional while ringingJeremy Magland 242 pf[32] = cp.center[0]; pf[33] = cp.center[1];
a2a1b68Drumstick by default, hit sequences, 170 GPa steelJeremy Magland 243 pf[34] = params.inPlaneMass; pf[35] = 1 / params.inPlaneMass;
d8f9e54Clean rebuild of the steel sheet simulatorJeremy Magland 244 this.paramsBuf = device.createBuffer({ size: pf.byteLength, usage: GPUBufferUsage.UNIFORM | GPUBufferUsage.COPY_DST });
245 device.queue.writeBuffer(this.paramsBuf, 0, pf);
247 // Per-step parameters, one struct per step at a dynamic offset.
248 this.stepStride = Math.max(STEP_LEN * 4, device.limits.minUniformBufferOffsetAlignment);
249 this.stepData = new ArrayBuffer(this.stepStride * this.stepsPerBatch);
250 this.stepBuf = device.createBuffer({ size: this.stepData.byteLength, usage: GPUBufferUsage.UNIFORM | GPUBufferUsage.COPY_DST });
252 // Pipelines.
253 const shader = (code: string) => device.createShaderModule({ code: commonWgsl + code });
254 const entry = (binding: number, type: GPUBufferBindingType, hasDynamicOffset = false): GPUBindGroupLayoutEntry => ({
255 binding,
256 visibility: GPUShaderStage.COMPUTE,
257 buffer: { type, hasDynamicOffset },
258 });
259 const pipeline = (module: GPUShaderModule, entryPoint: string, layout: GPUBindGroupLayout) =>
260 device.createComputePipeline({ layout: device.createPipelineLayout({ bindGroupLayouts: [layout] }), compute: { module, entryPoint } });
262 const stepLayout = device.createBindGroupLayout({
263 entries: [
264 entry(0, "uniform"),
265 entry(1, "uniform", true),
266 entry(2, "read-only-storage"),
267 entry(3, "read-only-storage"),
268 entry(4, "storage"),
269 entry(5, "storage"),
270 entry(6, "storage"),
271 entry(7, "storage"),
272 entry(8, "read-only-storage"),
273 entry(9, "read-only-storage"),
274 ],
275 });
276 this.stepPipe = pipeline(shader(stepWgsl), "main", stepLayout);
278 const rigidModule = shader(rigidWgsl);
279 const rbLayout = device.createBindGroupLayout({
280 entries: [entry(0, "uniform"), entry(1, "storage"), entry(2, "storage"), entry(3, "storage"), entry(9, "read-only-storage")],
281 });
282 this.sumsPipe = pipeline(rigidModule, "sums", rbLayout);
283 this.applyPipe = pipeline(rigidModule, "apply", rbLayout);
284 const projectLayout = device.createBindGroupLayout({
285 entries: [entry(0, "uniform"), entry(4, "storage"), entry(9, "read-only-storage")],
286 });
287 this.projectPipe = pipeline(rigidModule, "projectField", projectLayout);
289 const audioLayout = device.createBindGroupLayout({
290 entries: [entry(0, "uniform"), entry(1, "uniform", true), entry(2, "read-only-storage"), entry(3, "read-only-storage"), entry(4, "storage")],
291 });
292 this.audioPipe = pipeline(shader(audioWgsl), "main", audioLayout);
294 const extractLayout = device.createBindGroupLayout({
295 entries: [entry(0, "uniform"), entry(1, "read-only-storage"), entry(2, "storage")],
296 });
297 this.extractPipe = pipeline(shader(extractWgsl), "main", extractLayout);
299 // Bind groups. The state is double-buffered: bind group p reads pair p and writes pair 1 - p.
300 const sp = { buffer: this.stepBuf, size: STEP_LEN * 4 };
301 this.stepBG = [0, 1].map((p) =>
302 device.createBindGroup({
303 layout: stepLayout,
304 entries: [
305 { binding: 0, resource: { buffer: this.paramsBuf } },
306 { binding: 1, resource: sp },
307 { binding: 2, resource: { buffer: this.pos[p] } },
308 { binding: 3, resource: { buffer: this.vel[p] } },
309 { binding: 4, resource: { buffer: this.pos[1 - p] } },
310 { binding: 5, resource: { buffer: this.vel[1 - p] } },
311 { binding: 6, resource: { buffer: this.hist } },
312 { binding: 7, resource: { buffer: this.react } },
313 { binding: 8, resource: { buffer: this.w0 } },
314 { binding: 9, resource: { buffer: this.nw } },
315 ],
316 }),
317 );
318 this.rbBG = [0, 1].map((p) =>
319 device.createBindGroup({
320 layout: rbLayout,
321 entries: [
322 { binding: 0, resource: { buffer: this.paramsBuf } },
323 { binding: 1, resource: { buffer: this.pos[p] } },
324 { binding: 2, resource: { buffer: this.vel[p] } },
325 { binding: 3, resource: { buffer: this.statsBuf } },
326 { binding: 9, resource: { buffer: this.nw } },
327 ],
328 }),
329 );
330 this.projectBG = device.createBindGroup({
331 layout: projectLayout,
332 entries: [
333 { binding: 0, resource: { buffer: this.paramsBuf } },
334 { binding: 4, resource: { buffer: this.react } },
335 { binding: 9, resource: { buffer: this.nw } },
336 ],
337 });
338 this.audioBG = device.createBindGroup({
339 layout: audioLayout,
340 entries: [
341 { binding: 0, resource: { buffer: this.paramsBuf } },
342 { binding: 1, resource: sp },
343 { binding: 2, resource: { buffer: this.hist } },
344 { binding: 3, resource: { buffer: this.lw } },
345 { binding: 4, resource: { buffer: this.audioOut } },
346 ],
347 });
348 const extractBG = (p: number, out: GPUBuffer) =>
349 device.createBindGroup({
350 layout: extractLayout,
351 entries: [
352 { binding: 0, resource: { buffer: this.paramsBuf } },
353 { binding: 1, resource: { buffer: this.pos[p] } },
354 { binding: 2, resource: { buffer: out } },
355 ],
356 });
357 this.extractBG = [0, 1].map((p) => extractBG(p, this.frameBuf));
358 this.displayBG = [0, 1].map((p) => extractBG(p, this.displayW));
360 this.setListener(params.listenerDist, params.listenerAzimuth, params.listenerElevation);
361 if (params.excitation === "strike") {
362 // No preparation: the sheet is at rest in its natural shape when the stick arrives.
363 this.phase = "recording";
364 this.released = true;
365 this.frames.push(new Float32Array(n));
366 }
367 }
369 /** Recomputes the per-node weights and delays of the listener sum. Takes effect on later samples. */
370 setListener(dist: number, azimuthDeg: number, elevationDeg: number): void {
371 const g = this.grid;
372 const az = (azimuthDeg * Math.PI) / 180;
373 const el = (elevationDeg * Math.PI) / 180;
374 const L = [dist * Math.sin(az) * Math.cos(el), dist * Math.sin(el), dist * Math.cos(az) * Math.cos(el)];
375 const data = new Float32Array(2 * g.n);
376 let dmin = Infinity;
377 for (let j = 0; j < g.ny; j++) {
378 for (let i = 0; i < g.nx; i++) {
379 const k = j * g.nx + i;
380 const rx = L[0] - nodeX(g, i), ry = L[1] - nodeY(g, j), rz = L[2];
381 const r = Math.sqrt(rx * rx + ry * ry + rz * rz);
382 const dA = g.hx * g.hy * nodeWeight(g, i, j);
383 data[2 * k] = (RHO_AIR * dA) / (2 * Math.PI * r);
384 data[2 * k + 1] = (r / SPEED_OF_SOUND) * this.params.fs;
385 dmin = Math.min(dmin, data[2 * k + 1]);
386 }
387 }
388 // Only the spread of delays across the sheet is kept; the common delay is dropped.
389 this.delayClamped = false;
390 for (let k = 0; k < g.n; k++) {
391 let d = data[2 * k + 1] - dmin;
a2a1b68Drumstick by default, hit sequences, 170 GPa steelJeremy Magland 392 if (d > HIST_MARGIN - 3) {
393 d = HIST_MARGIN - 3;
395 }
396 data[2 * k + 1] = d;
397 }
398 this.device.queue.writeBuffer(this.lw, 0, data);
399 this.params.listenerDist = dist;
400 this.params.listenerAzimuth = azimuthDeg;
401 this.params.listenerElevation = elevationDeg;
402 }
404 setDamping(sigma0: number, sigma1: number, sigmaIn: number): void {
405 this.params.sigma0 = sigma0;
406 this.params.sigma1 = sigma1;
407 this.params.sigmaIn = sigmaIn;
408 }
410 setDuration(duration: number): void {
411 this.params.duration = duration;
412 if (this.canContinue()) this.phase = "recording";
413 }
415 /** True when a finished or stopped run can be extended. */
416 canContinue(): boolean {
417 return this.phase === "done" && this.released && this.samples < this.targetSamples();
418 }
420 /** Ends the run; the state, audio, and movie are kept. */
421 stop(): void {
422 this.phase = "done";
423 }
425 targetSamples(): number {
426 return Math.round(this.params.duration * this.params.fs);
427 }
429 status(): SimStatus {
a2a1b68Drumstick by default, hit sequences, 170 GPa steelJeremy Magland 430 // Report the time of the state that has been read back, not the one enqueued.
431 const ahead = this.inFlight ? this.inFlight.steps * this.dt : 0;
432 const t = this.phase === "push" || this.phase === "hold" ? this.prepTime - ahead - this.prepDuration() : this.stepCount * this.dt - ahead;
434 phase: this.phase,
435 t,
436 totalSteps: this.totalSteps,
437 ke: this.ke,
438 wmax: this.wmax,
439 samples: this.samples,
440 frames: this.frames.length,
441 peakPressure: this.peakPressure,
442 delayClamped: this.delayClamped,
443 stick: this.stick,
a2a1b68Drumstick by default, hit sequences, 170 GPa steelJeremy Magland 444 hits: this.params.excitation === "strike" ? { done: this.hitSteps.filter((h) => h * this.dt <= t).length, total: this.hits.length } : null,
446 }
448 private prepDuration(): number {
449 return this.params.rampTime + this.params.holdMax;
450 }
452 /** Enqueues an extraction of the current w field into `displayW` (for the live view). */
453 encodeDisplay(encoder: GPUCommandEncoder): void {
454 const pass = encoder.beginComputePass();
455 pass.setPipeline(this.extractPipe);
456 pass.setBindGroup(0, this.displayBG[this.cur]);
457 pass.dispatchWorkgroups(Math.ceil(this.grid.n / 256));
458 pass.end();
459 }
a2a1b68Drumstick by default, hit sequences, 170 GPa steelJeremy Magland 461 private stepWorkgroups(): [number, number] {
462 return [Math.ceil(this.grid.nx / TILE), Math.ceil(this.grid.ny / TILE)];
463 }
d8f9e54Clean rebuild of the steel sheet simulatorJeremy Magland 465 /** The recorded listener pressure so far [Pa]. */
466 recordedAudio(): Float32Array {
467 return this.audio.subarray(0, this.samples);
468 }
470 private appendAudio(raw: Float32Array): void {
471 const need = this.samples + raw.length;
472 if (this.audio.length < need) {
473 const a = new Float32Array(Math.max(need, 2 * this.audio.length, this.params.fs));
474 a.set(this.audio.subarray(0, this.samples));
475 this.audio = a;
476 }
477 this.audio.set(raw, this.samples);
478 for (let k = 0; k < raw.length; k++) this.peakPressure = Math.max(this.peakPressure, Math.abs(raw[k]));
479 this.samples = need;
480 }
482 /** Writes the per-step uniform for step s of a batch. */
a2a1b68Drumstick by default, hit sequences, 170 GPa steelJeremy Magland 483 private writeStepParams(s: number, fields: { flags: number; slot: number; outIdx: number; pushW: number; pushV: number; sigma0: number; releaseGain: number; strikeV?: number }): void {
d8f9e54Clean rebuild of the steel sheet simulatorJeremy Magland 484 const o = (s * this.stepStride) / 4;
485 const u32 = new Uint32Array(this.stepData);
486 const i32 = new Int32Array(this.stepData);
487 const f32 = new Float32Array(this.stepData);
488 u32[o] = this.totalSteps + s;
489 u32[o + 1] = fields.flags;
490 i32[o + 2] = fields.slot;
491 i32[o + 3] = fields.outIdx;
492 f32[o + 4] = fields.pushW;
493 f32[o + 5] = fields.pushV;
494 f32[o + 6] = fields.sigma0;
495 f32[o + 7] = this.params.sigma1;
496 f32[o + 8] = fields.releaseGain;
497 f32[o + 9] = this.params.sigmaIn;
498 u32[o + 10] = (this.totalSteps + s) & 1;
a2a1b68Drumstick by default, hit sequences, 170 GPa steelJeremy Magland 499 f32[o + 11] = fields.strikeV ?? 0;
503 * Runs one batch (one movie frame of time). Resolves when its results have
504 * been read back. The following batch is encoded before this one is awaited,
505 * unless its kind depends on this one's result.
506 */
d8f9e54Clean rebuild of the steel sheet simulatorJeremy Magland 507 async runBatch(): Promise<void> {
508 if (this.phase === "done") return;
a2a1b68Drumstick by default, hit sequences, 170 GPa steelJeremy Magland 509 const cur = this.inFlight ?? this.encodeBatch();
510 this.inFlight = null;
511 if (!cur) {
513 return;
514 }
a2a1b68Drumstick by default, hit sequences, 170 GPa steelJeremy Magland 515 if (this.canEncodeAhead()) this.inFlight = this.encodeBatch();
516 await this.collect(cur);
517 }
519 /** True when the batch after the one in flight is certainly needed and of a known kind. */
520 private canEncodeAhead(): boolean {
521 if (this.phase === "recording") return this.targetSamples() * this.stepsPerSample - this.stepCount > 0;
522 return !this.releaseRequested && this.prepTime < this.prepDuration();
523 }
525 /** Encodes and submits one batch; null when no steps remain. */
526 private encodeBatch(): Pending | null {
527 const p = this.params;
528 const recording = this.phase === "recording";
529 const stepsLeft = recording ? this.targetSamples() * this.stepsPerSample - this.stepCount : this.stepsPerBatch;
530 const steps = Math.min(this.stepsPerBatch, Math.max(0, stepsLeft));
531 if (steps === 0) return null;
d8f9e54Clean rebuild of the steel sheet simulatorJeremy Magland 532
533 // Per-step parameters: the push ramp and settle damping before the release,
a2a1b68Drumstick by default, hit sequences, 170 GPa steelJeremy Magland 534 // the release ramp or the drumstick (re-armed at each scheduled hit) and
535 // the audio bookkeeping afterwards.
d8f9e54Clean rebuild of the steel sheet simulatorJeremy Magland 536 const audioSteps: number[] = [];
537 for (let s = 0; s < steps; s++) {
538 const local = this.stepCount + s;
a2a1b68Drumstick by default, hit sequences, 170 GPa steelJeremy Magland 539 let flags = 0, pushW = 0, pushV = 0, sigma0 = p.sigma0, slot = -1, outIdx = -1, releaseGain = 0, strikeV = 0;
541 const t = this.prepTime + (s + 1) * this.dt; // state after this step
542 const tau = Math.min(1, t / p.rampTime);
543 const smooth = tau * tau * (3 - 2 * tau);
544 const dsmooth = tau < 1 ? (6 * tau * (1 - tau)) / p.rampTime : 0;
545 flags = 1 | (tau >= 1 ? 2 : 0); // record the reaction once the ramp is complete
546 pushW = p.pushAmount * smooth;
547 pushV = p.pushAmount * dsmooth;
548 sigma0 = p.settleSigma0;
549 } else {
550 const t = local * this.dt;
551 if (p.excitation === "strike") {
552 flags = 8;
a2a1b68Drumstick by default, hit sequences, 170 GPa steelJeremy Magland 553 if (this.nextHit < this.hitSteps.length && this.hitSteps[this.nextHit] === local) {
554 flags |= 16;
555 strikeV = this.hits[this.nextHit].v;
556 this.nextHit++;
557 }
d8f9e54Clean rebuild of the steel sheet simulatorJeremy Magland 558 } else if (p.releaseTime > 0 && t < p.releaseTime) {
559 const tau = t / p.releaseTime;
560 releaseGain = 1 - tau * tau * (3 - 2 * tau);
561 flags = 4;
562 }
563 if (local % this.stepsPerSample === 0) {
a2a1b68Drumstick by default, hit sequences, 170 GPa steelJeremy Magland 564 slot = (local / this.stepsPerSample) % this.histLen;
566 audioSteps.push(s);
567 }
568 }
a2a1b68Drumstick by default, hit sequences, 170 GPa steelJeremy Magland 569 this.writeStepParams(s, { flags, slot, outIdx, pushW, pushV, sigma0, releaseGain, strikeV });
571 this.device.queue.writeBuffer(this.stepBuf, 0, this.stepData, 0, steps * this.stepStride);
573 // Encode the batch in chunks so that the GPU starts while JS is still encoding.
574 const CHUNK = 1024;
575 let enc = this.device.createCommandEncoder();
576 let pass = enc.beginComputePass();
a2a1b68Drumstick by default, hit sequences, 170 GPa steelJeremy Magland 577 const [wgX, wgY] = this.stepWorkgroups();
d8f9e54Clean rebuild of the steel sheet simulatorJeremy Magland 578 const wgN = Math.ceil(this.grid.n / 256);
579 let cur = this.cur;
282f706Single excitation panel; rigid-body projection optional while ringingJeremy Magland 580 // While held, the projection defines the inertia-relief holding and runs
581 // every RB_INTERVAL_HELD steps. While ringing it is optional; the sums
582 // still run once per batch for the statistics.
583 const project = !recording || p.rigidProjection;
d8f9e54Clean rebuild of the steel sheet simulatorJeremy Magland 584 const rbInterval = recording ? RB_INTERVAL : RB_INTERVAL_HELD;
585 let stepBound = false;
586 for (let s = 0; s < steps; s++) {
587 if (s > 0 && s % CHUNK === 0) {
588 pass.end();
589 this.device.queue.submit([enc.finish()]);
590 enc = this.device.createCommandEncoder();
591 pass = enc.beginComputePass();
592 stepBound = false;
593 }
594 if (!stepBound) {
595 pass.setPipeline(this.stepPipe);
596 stepBound = true;
597 }
598 pass.setBindGroup(0, this.stepBG[cur], [s * this.stepStride]);
599 pass.dispatchWorkgroups(wgX, wgY);
600 cur = 1 - cur;
601 const global = this.totalSteps + s + 1;
282f706Single excitation panel; rigid-body projection optional while ringingJeremy Magland 602 if ((project && global % rbInterval === 0) || s === steps - 1) {
d8f9e54Clean rebuild of the steel sheet simulatorJeremy Magland 603 pass.setPipeline(this.sumsPipe);
604 pass.setBindGroup(0, this.rbBG[cur]);
605 pass.dispatchWorkgroups(1);
282f706Single excitation panel; rigid-body projection optional while ringingJeremy Magland 606 if (project) {
607 pass.setPipeline(this.applyPipe);
608 pass.setBindGroup(0, this.rbBG[cur]);
609 pass.dispatchWorkgroups(wgN);
610 }
612 }
614 // The listener sum for every audio sample of the batch, one workgroup each;
615 // the ring buffer is long enough that no slot of the batch has been reused.
616 const nAudio = audioSteps.length;
617 if (nAudio > 0) {
618 pass.setPipeline(this.audioPipe);
619 pass.setBindGroup(0, this.audioBG, [audioSteps[0] * this.stepStride]);
620 pass.dispatchWorkgroups(nAudio);
622 pass.setPipeline(this.extractPipe);
623 pass.setBindGroup(0, this.extractBG[cur]);
624 pass.dispatchWorkgroups(wgN);
625 pass.end();
626 const staging = this.staging[this.stagingIdx];
627 this.stagingIdx = 1 - this.stagingIdx;
628 const n = this.grid.n;
629 enc.copyBufferToBuffer(this.audioOut, 0, staging, 0, Math.max(4, nAudio * 4));
630 enc.copyBufferToBuffer(this.frameBuf, 0, staging, this.samplesPerBatch * 4, n * 4);
631 enc.copyBufferToBuffer(this.statsBuf, 0, staging, (this.samplesPerBatch + n) * 4, STATS_LEN * 4);
632 enc.copyBufferToBuffer(this.react, n * 4, staging, (this.samplesPerBatch + n + STATS_LEN) * 4, STICK_LEN * 4);
633 this.device.queue.submit([enc.finish()]);
634 this.cur = cur;
635 this.totalSteps += steps;
636 this.stepCount += steps;
637 if (!recording) this.prepTime += steps * this.dt;
a2a1b68Drumstick by default, hit sequences, 170 GPa steelJeremy Magland 638 return { staging, steps, nAudio, recording, prepTimeEnd: this.prepTime, stickParity: this.totalSteps & 1 };
639 }
d8f9e54Clean rebuild of the steel sheet simulatorJeremy Magland 640
a2a1b68Drumstick by default, hit sequences, 170 GPa steelJeremy Magland 641 /** Reads back a batch and applies its results and the phase transitions they imply. */
642 private async collect(b: Pending): Promise<void> {
643 const p = this.params;
644 const n = this.grid.n;
645 await b.staging.mapAsync(GPUMapMode.READ);
646 const view = new Float32Array(b.staging.getMappedRange());
d8f9e54Clean rebuild of the steel sheet simulatorJeremy Magland 647 const stats = view.subarray(this.samplesPerBatch + n, this.samplesPerBatch + n + STATS_LEN);
648 this.ke = stats[12];
649 this.wmax = stats[13];
650 if (p.excitation === "strike") {
651 // The slot written by the last step is the one the next step reads.
a2a1b68Drumstick by default, hit sequences, 170 GPa steelJeremy Magland 652 const base = this.samplesPerBatch + n + STATS_LEN + 4 * b.stickParity;
d8f9e54Clean rebuild of the steel sheet simulatorJeremy Magland 653 this.stick = { x: view[base], v: view[base + 1], force: view[base + 2], maxForce: view[base + 3] };
654 }
655 const frame = view.slice(this.samplesPerBatch, this.samplesPerBatch + n);
657 this.appendAudio(view.slice(0, b.nAudio));
659 } else {
660 this.lastPrepFrame = frame;
661 }
664 if (!b.recording) {
665 if (this.phase === "push" && b.prepTimeEnd >= p.rampTime) this.phase = "hold";
666 if (this.phase === "hold" && !this.releaseRequested && this.settled(b.prepTimeEnd)) this.releaseRequested = true;
667 if (b.prepTimeEnd >= this.prepDuration()) this.releaseRequested = true;
668 // Release only once every preparation batch has been read back, so that
669 // frame 0 is the last held state.
670 if (this.releaseRequested && !this.inFlight) this.release();
d8f9e54Clean rebuild of the steel sheet simulatorJeremy Magland 671 } else if (this.samples >= this.targetSamples()) {
672 this.phase = "done";
673 }
674 }
677 * Settle criterion of the held sheet at hold time t: the non-rigid kinetic
678 * energy is below settleKE, or it has stopped falling (it has not halved
679 * over the last window). The second case catches the rounding floor, which
680 * can lie above any fixed threshold.
681 */
682 private settled(t: number): boolean {
683 const p = this.params;
684 this.keHistory.push({ t, ke: this.ke });
685 if (this.ke < p.settleKE) return true;
686 const window = Math.max(SETTLE_WINDOW, 3 / p.settleSigma0);
687 let ref: { t: number; ke: number } | null = null;
688 for (const h of this.keHistory) if (h.t <= t - window) ref = h;
689 return ref !== null && this.ke >= SETTLE_RATIO * ref.ke;
690 }
693 this.phase = "recording";
694 this.released = true;
a2a1b68Drumstick by default, hit sequences, 170 GPa steelJeremy Magland 695 this.releaseRequested = false;
696 this.keHistory = [];
698 this.totalSteps = 0;
699 // The finger force and the distributed holding reaction relax together:
700 // remove the rigid-body components of the recorded reaction.
701 const enc = this.device.createCommandEncoder();
702 const pass = enc.beginComputePass();
703 pass.setPipeline(this.projectPipe);
704 pass.setBindGroup(0, this.projectBG);
705 pass.dispatchWorkgroups(1);
706 pass.end();
707 this.device.queue.submit([enc.finish()]);
708 // Frame 0 is the held, settled shape.
709 if (this.lastPrepFrame) this.frames.push(this.lastPrepFrame);
710 }
712 // ---------- Testing support ----------
714 /**
715 * Advances the state by exactly one free step (no constraint, no drumstick,
716 * no rigid-body projection, no audio) with the current loss parameters.
717 */
718 async stepOnce(): Promise<void> {
719 this.writeStepParams(0, { flags: 0, slot: -1, outIdx: -1, pushW: 0, pushV: 0, sigma0: this.params.sigma0, releaseGain: 0 });
720 this.device.queue.writeBuffer(this.stepBuf, 0, this.stepData, 0, this.stepStride);
721 const enc = this.device.createCommandEncoder();
722 const pass = enc.beginComputePass();
723 pass.setPipeline(this.stepPipe);
724 pass.setBindGroup(0, this.stepBG[this.cur], [0]);
a2a1b68Drumstick by default, hit sequences, 170 GPa steelJeremy Magland 725 pass.dispatchWorkgroups(...this.stepWorkgroups());
727 this.device.queue.submit([enc.finish()]);
728 this.cur = 1 - this.cur;
729 this.totalSteps++;
730 await this.device.queue.onSubmittedWorkDone();
731 }
733 /** Overwrites the state: pos as [u | v | w | wLo], vel as [du | dv | dw]. */
734 writeState(pos: Float32Array<ArrayBuffer>, vel: Float32Array<ArrayBuffer>): void {
735 const n = this.grid.n;
736 if (pos.length !== 4 * n || vel.length !== 3 * n) throw new Error("state arrays have the wrong length");
737 this.device.queue.writeBuffer(this.pos[this.cur], 0, pos);
738 this.device.queue.writeBuffer(this.vel[this.cur], 0, vel);
739 }
741 /** Reads back the full state as [u | v | w | wLo] and [du | dv | dw]. */
742 async readState(): Promise<{ pos: Float32Array; vel: Float32Array }> {
743 const n = this.grid.n;
744 const stg = this.device.createBuffer({ size: 7 * n * 4, usage: GPUBufferUsage.MAP_READ | GPUBufferUsage.COPY_DST });
745 const enc = this.device.createCommandEncoder();
746 enc.copyBufferToBuffer(this.pos[this.cur], 0, stg, 0, 4 * n * 4);
747 enc.copyBufferToBuffer(this.vel[this.cur], 0, stg, 4 * n * 4, 3 * n * 4);
748 this.device.queue.submit([enc.finish()]);
749 await stg.mapAsync(GPUMapMode.READ);
750 const all = new Float32Array(stg.getMappedRange()).slice();
751 stg.unmap();
752 stg.destroy();
753 return { pos: all.subarray(0, 4 * n), vel: all.subarray(4 * n) };
754 }
756 destroy(): void {
d8f9e54Clean rebuild of the steel sheet simulatorJeremy Magland 758 for (const b of [...this.pos, ...this.vel, this.hist, this.react, this.w0, this.nw, this.lw, this.audioOut, this.frameBuf, this.displayW, this.statsBuf, this.paramsBuf, this.stepBuf, ...this.staging]) b.destroy();
759 }
761 get fs(): number {
762 return this.params.fs;
763 }
764 get fps(): number {
765 return this.params.fps;
766 }
767 get physParams(): PhysParams {
768 return this.params;
769 }
770}