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
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.
13import {
14 type PhysParams,
15 type Grid,
16 type Material,
17 makeGrid,
18 makeMaterial,
19 strikePoint,
20 pushPoint,
21 contactPatch,
22 hitSchedule,
23 type Hit,
24 nodeX,
25 nodeY,
26 nodeWeight,
27 nodeWeights,
28 CFL,
29 HIST_MARGIN,
30 RB_INTERVAL,
31 RB_INTERVAL_HELD,
32 RHO_AIR,
33 SETTLE_RATIO,
34 SETTLE_WINDOW,
35 SPEED_OF_SOUND,
36} from "../params.ts";
37import { initialShape } from "../shape.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;
67 /** Hits delivered so far, of the scheduled total (strike excitation only). */
68 hits: { done: number; total: number } | null;
69}
71const PARAMS_LEN = 36; // floats in the Params uniform
72const STEP_LEN = 12; // floats in the StepParams uniform
73const STATS_LEN = 16;
74const STICK_LEN = 8;
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}
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;
99 /** Ring buffer length of the listener sum [samples]. */
100 readonly histLen: number;
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;
115 /** Steps enqueued since the start of the current clock (preparation or recording). */
116 private stepCount = 0;
117 /** Preparation time enqueued [s]. */
118 private prepTime = 0;
119 private released = false;
120 private lastPrepFrame: Float32Array<ArrayBuffer> | null = null;
121 private stick: StickState | null = null;
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 }[] = [];
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 /** Largest magnitude of the initial shape [m]. */
145 readonly w0Max: number;
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.
182 this.omegaMax = maxAngularFrequency(g, m, params.inPlaneMass);
183 const dtMax = (CFL * 2) / this.omegaMax;
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;
188 this.histLen = this.samplesPerBatch + HIST_MARGIN;
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]
196 this.hist = f32(this.histLen * n);
197 this.react = f32(n + STICK_LEN);
198 this.w0 = f32(n);
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);
202 this.nw = f32(n);
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 }
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).
226 const cp = contactPatch(params, g);
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)));
230 const stick = new Float32Array(STICK_LEN);
231 stick[1] = -this.hits[0].v; stick[5] = -this.hits[0].v;
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);
238 pu[0] = g.nx; pu[1] = g.ny; pu[2] = n; pu[3] = this.histLen;
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];
246 pf[32] = cp.center[0]; pf[33] = cp.center[1];
247 pf[34] = params.inPlaneMass; pf[35] = 1 / params.inPlaneMass;
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;
396 if (d > HIST_MARGIN - 3) {
397 d = HIST_MARGIN - 3;
398 this.delayClamped = true;
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 {
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;
437 return {
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,
448 hits: this.params.excitation === "strike" ? { done: this.hitSteps.filter((h) => h * this.dt <= t).length, total: this.hits.length } : null,
449 };
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 }
465 private stepWorkgroups(): [number, number] {
466 return [Math.ceil(this.grid.nx / TILE), Math.ceil(this.grid.ny / TILE)];
467 }
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. */
487 private writeStepParams(s: number, fields: { flags: number; slot: number; outIdx: number; pushW: number; pushV: number; sigma0: number; releaseGain: number; strikeV?: number }): void {
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;
503 f32[o + 11] = fields.strikeV ?? 0;
504 }
506 /**
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 */
511 async runBatch(): Promise<void> {
512 if (this.phase === "done") return;
513 const cur = this.inFlight ?? this.encodeBatch();
514 this.inFlight = null;
515 if (!cur) {
516 this.phase = "done";
517 return;
518 }
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;
537 // Per-step parameters: the push ramp and settle damping before the release,
538 // the release ramp or the drumstick (re-armed at each scheduled hit) and
539 // the audio bookkeeping afterwards.
540 const audioSteps: number[] = [];
541 for (let s = 0; s < steps; s++) {
542 const local = this.stepCount + s;
543 let flags = 0, pushW = 0, pushV = 0, sigma0 = p.sigma0, slot = -1, outIdx = -1, releaseGain = 0, strikeV = 0;
544 if (!recording) {
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;
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 }
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) {
568 slot = (local / this.stepsPerSample) % this.histLen;
569 outIdx = audioSteps.length;
570 audioSteps.push(s);
571 }
572 }
573 this.writeStepParams(s, { flags, slot, outIdx, pushW, pushV, sigma0, releaseGain, strikeV });
574 }
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();
581 const [wgX, wgY] = this.stepWorkgroups();
582 const wgN = Math.ceil(this.grid.n / 256);
583 let cur = this.cur;
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;
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;
606 if ((project && global % rbInterval === 0) || s === steps - 1) {
607 pass.setPipeline(this.sumsPipe);
608 pass.setBindGroup(0, this.rbBG[cur]);
609 pass.dispatchWorkgroups(1);
610 if (project) {
611 pass.setPipeline(this.applyPipe);
612 pass.setBindGroup(0, this.rbBG[cur]);
613 pass.dispatchWorkgroups(wgN);
614 }
615 stepBound = false;
616 }
617 }
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);
625 }
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;
642 return { staging, steps, nAudio, recording, prepTimeEnd: this.prepTime, stickParity: this.totalSteps & 1 };
643 }
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());
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.
656 const base = this.samplesPerBatch + n + STATS_LEN + 4 * b.stickParity;
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);
660 if (b.recording) {
661 this.appendAudio(view.slice(0, b.nAudio));
662 this.frames.push(frame);
663 } else {
664 this.lastPrepFrame = frame;
665 }
666 b.staging.unmap();
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();
675 } else if (this.samples >= this.targetSamples()) {
676 this.phase = "done";
677 }
678 }
680 /**
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 }
696 private release(): void {
697 this.phase = "recording";
698 this.released = true;
699 this.releaseRequested = false;
700 this.keHistory = [];
701 this.stepCount = 0;
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]);
729 pass.dispatchWorkgroups(...this.stepWorkgroups());
730 pass.end();
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 {
761 this.inFlight = null;
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}