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