65ecc03Run the time loop on the GPU in the web app when WebGPU is availableJeremy Magland 1/// <reference types="@webgpu/types" />
2// WebGPU implementation of the cymbal2 time loop (cymbal2_sim.m), ported from
3// the prototype in proof-of-concept/cymbal2-webgpu. Everything around the
4// loop stays in numbl: the caller runs the unmodified MATLAB setup stage
5// (grid, per-wavenumber inverses, mallet, stroke schedule) and hands the
6// tensors in through GpuSetup; the acceleration and contact-force histories
7// come back for the unmodified output stage.
8//
9// One time step is eight compute dispatches at full GPU width, the dispatch
10// boundaries providing global synchronization between stages:
11// AB hand state machine + Hertz contact + beam update (workgroup 0), and
12// the shell derivative stencils / Airy right-hand side (all workgroups)
13// C azimuthal analysis of the Airy right-hand side
14// D Airy solve, one small radial system per wavenumber
15// E azimuthal synthesis of the stress function
16// FG von Karman bracket + forcing (stress-function stencils inline)
17// H azimuthal analysis of the nonlinear load
18// I implicit shell update, three radial products per wavenumber
19// J azimuthal synthesis of the new displacement; pickup acceleration
20//
21// Every recorded dispatch is identical: the step index k is a chunk base
22// (state cell 0, written by the host once per chunk) plus a within-chunk
23// offset from a dynamic-offset uniform, and all slot rotation is a pure
24// function of k. Slots rotate every step whether or not the shell is active,
25// which is exact: before first contact every slot is zero, and after first
26// contact every step writes. The shell-active flag latches inside AB and is
27// read by the later stages of the same step one dispatch later; the benign
28// race inside AB itself can only skip work on a field that is still zero.
30export interface GpuSetup {
31 ops: {
32 D1: Float32Array; D2: Float32Array;
33 F: Float32Array; Gi: Float32Array;
34 Rd: Float32Array; R2d: Float32Array;
35 AiB: Float32Array; AiC: Float32Array; Ai: Float32Array; Fi: Float32Array;
36 ASiB: Float32Array; ASiC: Float32Array;
37 vM: Float32Array; vR: Float32Array; vK: Float32Array; tipcol: Float32Array;
38 };
39 tlaunch: number[];
40 vstr: number[];
41 nrt: number; Nr: number; nth: number;
42 hr: number; hth: number;
43 ht: number; Eh: number; Rs: number; rhoh: number;
44 fsc: number; nlr0: number;
45 ls: number; js: number; lp: number; jp: number;
46 tipi: number; nfree: number; ps: number;
47 KM: number; alM: number; aH: number; aHr: number;
48 holdSteps: number; missDepth: number; gap: number;
49 K: number; nsteps: number;
50}
52export interface GpuSimResult {
53 acc: Float32Array;
54 fhist: Float32Array;
55 peakXi: number;
56 peakDisp: number;
57 contactVel: number[];
58 strokesPlayed: number;
59 gpuMs: number;
60}
62const fmt = (x: number): string => {
63 if (!Number.isFinite(x)) throw new Error('non-finite kernel constant');
64 let s = x.toPrecision(9);
65 if (!/[.e]/.test(s)) s += '.0';
66 return s + 'f';
67};
69interface Layout {
70 N: number;
71 opsData: Float32Array<ArrayBuffer>;
72 stLen: number;
73 SC_CVEL: number;
74}
76function buildLayoutAndStages(su: GpuSetup, wg: number) {
77 const nrt = su.nrt | 0;
78 const nth = su.nth | 0;
79 const nf = su.nfree | 0;
80 const N = nrt * nth;
81 const K = su.K | 0;
82 const ht = su.ht;
84 // ---- packed operator buffer ----
85 const parts: [string, Float32Array, number][] = [
86 ['D1', su.ops.D1, nrt * nrt],
87 ['D2', su.ops.D2, nrt * nrt],
88 ['F', su.ops.F, nth * nth],
89 ['GI', su.ops.Gi, nth * nth],
90 ['RD', su.ops.Rd, nrt],
91 ['R2D', su.ops.R2d, nrt],
92 ['AIB', su.ops.AiB, nrt * nrt * nth],
93 ['AIC', su.ops.AiC, nrt * nrt * nth],
94 ['AI', su.ops.Ai, nrt * nrt * nth],
95 ['FI', su.ops.Fi, nrt * nrt * nth],
96 ['ASIB', su.ops.ASiB, nf * nf],
97 ['ASIC', su.ops.ASiC, nf * nf],
98 ['VM', su.ops.vM, nf],
99 ['VR', su.ops.vR, nf],
100 ['VK', su.ops.vK, nf],
101 ['TIPCOL', su.ops.tipcol, nf],
102 ];
103 // Launch-step thresholds: smallest 0-based kk with (kk+1)*ht >= tlaunch,
104 // decided in f64 exactly as MATLAB's `k*ht >= tlaunch` decides it.
105 const kgo = su.tlaunch.map((tl) => {
106 let g = Math.max(0, Math.ceil(tl / ht) - 1);
107 while (g > 0 && g * ht >= tl) g -= 1;
108 while ((g + 1) * ht < tl) g += 1;
109 return g;
110 });
111 parts.push(['KGO', Float32Array.from(kgo), K]);
112 parts.push(['VSTR', Float32Array.from(su.vstr), K]);
113 // Transposed beam matrices so each beam-row thread streams its own row.
114 const tr = (m: Float32Array) => {
115 const t = new Float32Array(nf * nf);
116 for (let r = 0; r < nf; r++) for (let p = 0; p < nf; p++) t[p + r * nf] = m[r + p * nf];
117 return t;
118 };
119 parts.push(['ASIBT', tr(su.ops.ASiB), nf * nf]);
120 parts.push(['ASICT', tr(su.ops.ASiC), nf * nf]);
122 const off: Record<string, number> = {};
123 let cursor = 0;
124 for (const [name, arr, len] of parts) {
125 if (arr.length !== len) throw new Error(`ops part ${name}: ${arr.length} != ${len}`);
126 off[name] = cursor;
127 cursor += len;
128 }
129 const opsData = new Float32Array(cursor);
130 for (const [name, arr] of parts) opsData.set(arr, off[name]);
132 // ---- state buffer: 16 scalars, cvel, 3 beam slots, 3+3 field slots, scratch ----
133 const ST_CVEL = 16;
134 const ST_D = ST_CVEL + K;
135 const ST_U = ST_D + 3 * nf;
136 const ST_UH = ST_U + 3 * N;
137 const scratchNames = ['URR', 'V1', 'V2', 'RHSA', 'RHSAF', 'PH3', 'PHI', 'NL', 'NLF'];
138 const A: Record<string, number> = {};
139 let scur = ST_UH + 3 * N;
140 for (const n of scratchNames) {
141 A[n] = scur;
142 scur += N;
143 }
144 const stLen = scur;
146 const c: Record<string, number> = {
147 HT: ht,
148 HT2I: 1 / (ht * ht),
149 GAP: su.gap,
150 KM: su.KM,
151 ALM: su.alM,
152 PS: su.ps,
153 AH: su.aH,
154 AHR: su.aHr,
155 MISSD: su.missDepth,
156 EHRS: su.Eh / su.Rs,
157 EH2: su.Eh / 2,
158 RSI: 1 / su.Rs,
159 SNL: (ht * ht) / su.rhoh,
160 FSC: su.fsc,
161 C1R: 1 / (2 * su.hr),
162 C2R: 1 / (su.hr * su.hr),
163 C1T: 1 / (2 * su.hth),
164 C2T: 1 / (su.hth * su.hth),
165 };
167 const header = `
168const NRT = ${nrt}u;
169const NTH = ${nth}u;
170const N = ${N}u;
171const NF = ${nf}u;
172const NSTEPS = ${su.nsteps}u;
173const KS = ${K}u;
174const NRU = ${su.Nr}u;
175const LS0 = ${su.ls - 1}u; const JS0 = ${su.js - 1}u;
176const LP0 = ${su.lp - 1}u; const JP0 = ${su.jp - 1}u;
177const TIPI = ${su.tipi - 1}u;
178const NLR0 = ${su.nlr0 - 1}u;
179const HOLDSTEPS = ${su.holdSteps | 0};
180${Object.entries(c).map(([k, v]) => `const ${k}: f32 = ${fmt(v)};`).join('\n')}
181${Object.entries(off).map(([k, v]) => `const OP_${k} = ${v}u;`).join('\n')}
182const ST_CVEL = ${ST_CVEL}u;
183const ST_D = ${ST_D}u;
184const ST_U = ${ST_U}u;
185const ST_UH = ${ST_UH}u;
186${scratchNames.map((n) => `const A_${n} = ${A[n]}u;`).join('\n')}
188@group(0) @binding(0) var<storage, read> ops: array<f32>;
189@group(0) @binding(1) var<storage, read_write> st: array<f32>;
190@group(0) @binding(2) var<storage, read_write> outb: array<f32>;
191@group(0) @binding(3) var<uniform> ctl: vec4<u32>;
193fn thm(j: u32) -> u32 { return (j + NTH - 1u) % NTH; }
194fn thp(j: u32) -> u32 { return (j + 1u) % NTH; }
195`;
197 const AB = `${header}
198var<workgroup> wgF: f32;
199var<workgroup> wgHA: f32;
200var<workgroup> wgHB: f32;
201var<workgroup> wgHC: f32;
202var<workgroup> shD: array<f32, ${nf}>;
203var<workgroup> shD1: array<f32, ${nf}>;
205fn u_dr(b: u32, l: u32, j: u32) -> f32 {
206 if (l < 1u || l > NRT - 2u) { return 0.0; }
207 return (st[b + (l + 1u) + j * NRT] - st[b + (l - 1u) + j * NRT]) * C1R;
208}
210@compute @workgroup_size(${wg})
211fn main(@builtin(global_invocation_id) gid: vec3<u32>,
212 @builtin(local_invocation_index) lid: u32,
213 @builtin(workgroup_id) wgid: vec3<u32>) {
214 // No early-out on k: the host records exactly nsteps steps, and the
215 // workgroupBarrier below requires uniform control flow up to it.
216 let k = u32(st[0u]) + ctl.x;
217 let uS = ST_U + (k % 3u) * N;
218 if (wgid.x == 0u) {
219 if (lid == 0u) {
220 var shellA = u32(st[1u]); var inC = u32(st[2u]); var phase = u32(st[3u]);
221 var strokeIdx = u32(st[4u]); var holdLeft = i32(st[5u]); var vtar = st[6u];
222 var vH = st[7u]; var uHm = st[8u]; var uH0 = st[9u]; var uHp = st[10u];
223 var peakXi = st[12u]; var peakDisp = st[13u];
224 let dCur = ST_D + ((k + 2u) % 3u) * NF;
225 let dOld = ST_D + ((k + 1u) % 3u) * NF;
226 if (shellA != 0u && k > 0u && (k % 512u) == 0u) {
227 var pd: f32 = 0.0;
228 for (var l: u32 = 0u; l < NRT; l = l + 1u) {
229 pd = max(pd, abs(st[uS + l + JS0 * NRT]));
230 }
231 if (pd > peakDisp) { peakDisp = pd; }
232 }
233 let xi = st[uS + LS0 + JS0 * NRT] - st[dCur + TIPI] - GAP;
234 var fM: f32 = 0.0;
235 if (xi > 0.0) {
236 fM = KM * pow(xi, ALM);
237 shellA = 1u;
238 if (inC == 0u && strokeIdx < KS && st[ST_CVEL + strokeIdx] == 0.0) {
239 st[ST_CVEL + strokeIdx] = (st[dCur + TIPI] - st[dOld + TIPI]) / HT;
240 }
241 inC = 1u;
242 if (xi > peakXi) { peakXi = xi; }
243 } else {
244 inC = 0u;
245 }
246 outb[NSTEPS + k] = fM;
247 wgF = fM;
248 wgHA = -(uHp - 2.0 * uH0 + uHm);
249 wgHB = -(HT * 0.5) * (uHp - uHm);
250 wgHC = -HT * HT * (PS * uH0 + ((1.0 - PS) * 0.5) * (uHp + uHm));
251 if (phase == 0u) {
252 if (strokeIdx < KS && f32(k) >= ops[OP_KGO + strokeIdx]) {
253 phase = 1u; vtar = ops[OP_VSTR + strokeIdx];
254 }
255 } else if (phase == 1u) {
256 vH = max(vH - AH * HT, -vtar);
257 if (inC != 0u) { phase = 3u; holdLeft = HOLDSTEPS; }
258 else if (vH <= -vtar) { phase = 2u; }
259 } else if (phase == 2u) {
260 if (inC != 0u) { phase = 3u; holdLeft = HOLDSTEPS; }
261 else if (uHp <= -MISSD) { phase = 4u; }
262 } else if (phase == 3u) {
263 holdLeft = holdLeft - 1;
264 if (holdLeft <= 0) { phase = 4u; }
265 } else if (phase == 4u) {
266 vH = min(vH + AHR * HT, 0.5 * vtar);
267 if (vH > 0.0 && uHp >= -vH * vH / (2.0 * AHR)) { phase = 5u; }
268 } else {
269 vH = vH - AHR * HT;
270 if (vH <= 0.0 || uHp >= 0.0) { vH = 0.0; phase = 0u; strokeIdx = strokeIdx + 1u; }
271 }
272 uHm = uH0; uH0 = uHp; uHp = uHp + vH * HT;
273 st[1u] = f32(shellA); st[2u] = f32(inC); st[3u] = f32(phase);
274 st[4u] = f32(strokeIdx); st[5u] = f32(holdLeft); st[6u] = vtar;
275 st[7u] = vH; st[8u] = uHm; st[9u] = uH0; st[10u] = uHp;
276 st[12u] = peakXi; st[13u] = peakDisp; st[14u] = fM;
277 }
278 if (lid < NF) {
279 shD[lid] = st[ST_D + ((k + 2u) % 3u) * NF + lid];
280 shD1[lid] = st[ST_D + ((k + 1u) % 3u) * NF + lid];
281 }
282 workgroupBarrier();
283 if (lid < NF) {
284 let r = lid;
285 var sum: f32 = ops[OP_VM + r] * wgHA + ops[OP_VR + r] * wgHB + ops[OP_VK + r] * wgHC
286 + (HT * HT * wgF) * ops[OP_TIPCOL + r];
287 for (var p: u32 = 0u; p < NF; p = p + 1u) {
288 sum = sum + ops[OP_ASIBT + p + r * NF] * shD[p]
289 + ops[OP_ASICT + p + r * NF] * shD1[p];
290 }
291 st[ST_D + (k % 3u) * NF + r] = sum;
292 }
293 }
294 let i = gid.x;
295 if (i >= N) { return; }
296 if (st[1u] == 0.0) { return; }
297 let l = i % NRT; let j = i / NRT;
298 var ur: f32 = 0.0; var urr: f32 = 0.0;
299 if (l >= 1u && l <= NRT - 2u) {
300 let um = st[uS + (l - 1u) + j * NRT];
301 let u0 = st[uS + l + j * NRT];
302 let up = st[uS + (l + 1u) + j * NRT];
303 ur = (up - um) * C1R;
304 urr = (um - 2.0 * u0 + up) * C2R;
305 }
306 let ujm = st[uS + l + thm(j) * NRT];
307 let uj0 = st[uS + l + j * NRT];
308 let ujp = st[uS + l + thp(j) * NRT];
309 let v1 = ops[OP_RD + l] * ur + ops[OP_R2D + l] * ((ujm - 2.0 * uj0 + ujp) * C2T);
310 let urt = (u_dr(uS, l, thp(j)) - u_dr(uS, l, thm(j))) * C1T;
311 let v2 = ops[OP_RD + l] * urt - ops[OP_R2D + l] * ((ujp - ujm) * C1T);
312 st[A_URR + i] = urr;
313 st[A_V1 + i] = v1;
314 st[A_V2 + i] = v2;
315 var rv: f32 = 0.0;
316 if (l >= 3u && l <= NRU - 2u) {
317 rv = EHRS * (urr + v1) - EH2 * 2.0 * (urr * v1 - v2 * v2);
318 }
319 st[A_RHSA + i] = rv;
320}
321`;
323 const gate = `
324 let k = u32(st[0u]) + ctl.x;
325 if (k >= NSTEPS) { return; }
326 let i = gid.x;
327 if (i >= N) { return; }
328 if (st[1u] == 0.0) { return; }
329 let l = i % NRT; let j = i / NRT;`;
331 const transform = (src: string, dst: string, mat: string) => `${header}
332@compute @workgroup_size(${wg})
333fn main(@builtin(global_invocation_id) gid: vec3<u32>) {
334${gate}
335 var sum: f32 = 0.0;
336 for (var q: u32 = 0u; q < NTH; q = q + 1u) {
337 sum = sum + st[${src} + l + q * NRT] * ops[${mat} + q + j * NTH];
338 }
339 st[${dst} + i] = sum;
340}
341`;
343 const D = `${header}
344@compute @workgroup_size(${wg})
345fn main(@builtin(global_invocation_id) gid: vec3<u32>) {
346${gate}
347 let base = OP_FI + j * NRT * NRT;
348 var sum: f32 = 0.0;
349 for (var p: u32 = 0u; p < NRT; p = p + 1u) {
350 sum = sum + ops[base + l + p * NRT] * st[A_RHSAF + p + j * NRT];
351 }
352 st[A_PH3 + i] = sum;
353}
354`;
356 const FG = `${header}
357fn p_dr(l: u32, j: u32) -> f32 {
358 if (l < 1u || l > NRT - 2u) { return 0.0; }
359 return (st[A_PHI + (l + 1u) + j * NRT] - st[A_PHI + (l - 1u) + j * NRT]) * C1R;
360}
361@compute @workgroup_size(${wg})
362fn main(@builtin(global_invocation_id) gid: vec3<u32>) {
363${gate}
364 var nl: f32 = 0.0;
365 if (l >= NLR0 && l < NRU) {
366 let pm = st[A_PHI + (l - 1u) + j * NRT];
367 let p0 = st[A_PHI + l + j * NRT];
368 let pp = st[A_PHI + (l + 1u) + j * NRT];
369 let pr = (pp - pm) * C1R;
370 let prr = (pm - 2.0 * p0 + pp) * C2R;
371 let pjm = st[A_PHI + l + thm(j) * NRT];
372 let pjp = st[A_PHI + l + thp(j) * NRT];
373 let w1 = ops[OP_RD + l] * pr + ops[OP_R2D + l] * ((pjm - 2.0 * p0 + pjp) * C2T);
374 let prt = (p_dr(l, thp(j)) - p_dr(l, thm(j))) * C1T;
375 let w2 = ops[OP_RD + l] * prt - ops[OP_R2D + l] * ((pjp - pjm) * C1T);
376 nl = prr * st[A_V1 + i] + st[A_URR + i] * w1 - 2.0 * w2 * st[A_V2 + i]
377 - (prr + w1) * RSI;
378 if (l == LS0 && j == JS0) { nl = nl - st[14u] * FSC; }
379 nl = nl * SNL;
380 }
381 st[A_NL + i] = nl;
382}
383`;
385 const I = `${header}
386@compute @workgroup_size(${wg})
387fn main(@builtin(global_invocation_id) gid: vec3<u32>) {
388${gate}
389 let h3S = ST_UH + (k % 3u) * N;
390 let h13S = ST_UH + ((k + 2u) % 3u) * N;
391 let h3nS = ST_UH + ((k + 1u) % 3u) * N;
392 let base = j * NRT * NRT;
393 var sum: f32 = 0.0;
394 for (var p: u32 = 0u; p < NRT; p = p + 1u) {
395 let pj = p + j * NRT;
396 sum = sum + ops[OP_AIB + base + l + p * NRT] * st[h3S + pj]
397 + ops[OP_AIC + base + l + p * NRT] * st[h13S + pj]
398 + ops[OP_AI + base + l + p * NRT] * st[A_NLF + pj];
399 }
400 st[h3nS + i] = sum;
401}
402`;
404 const J = `${header}
405@compute @workgroup_size(${wg})
406fn main(@builtin(global_invocation_id) gid: vec3<u32>) {
407${gate}
408 let h3nS = ST_UH + ((k + 1u) % 3u) * N;
409 let uS = ST_U + (k % 3u) * N;
410 let u1S = ST_U + ((k + 2u) % 3u) * N;
411 let unS = ST_U + ((k + 1u) % 3u) * N;
412 var sum: f32 = 0.0;
413 for (var q: u32 = 0u; q < NTH; q = q + 1u) {
414 sum = sum + st[h3nS + l + q * NRT] * ops[OP_GI + q + j * NTH];
415 }
416 st[unS + i] = sum;
417 if (l == LP0 && j == JP0) {
418 outb[k] = (sum - 2.0 * st[uS + i] + st[u1S + i]) * HT2I;
419 }
420}
421`;
423 const sources = [
424 AB,
425 transform('A_RHSA', 'A_RHSAF', 'OP_F'),
426 D,
427 transform('A_PH3', 'A_PHI', 'OP_GI'),
428 FG,
429 transform('A_NL', 'A_NLF', 'OP_F'),
430 I,
431 J,
432 ];
434 const layout: Layout = { N, opsData, stLen, SC_CVEL: ST_CVEL };
435 return { layout, sources, nwg: Math.ceil(N / wg) };
436}
438export async function requestGpuDevice(): Promise<GPUDevice> {
439 const gpu = (navigator as { gpu?: GPU }).gpu;
440 if (!gpu) throw new Error('WebGPU is not available in this browser');
441 const adapter = await gpu.requestAdapter();
442 if (!adapter) throw new Error('no WebGPU adapter');
443 return adapter.requestDevice();
444}
446export async function runGpuSim(
447 device: GPUDevice,
448 su: GpuSetup,
449 onProgress?: (stepsDone: number, nsteps: number) => void,
450): Promise<GpuSimResult> {
451 const wg = 64;
452 const chunk = 1024;
453 const nsteps = su.nsteps;
454 const { layout, sources, nwg } = buildLayoutAndStages(su, wg);
456 const bgl = device.createBindGroupLayout({
457 entries: [
458 { binding: 0, visibility: GPUShaderStage.COMPUTE, buffer: { type: 'read-only-storage' } },
459 { binding: 1, visibility: GPUShaderStage.COMPUTE, buffer: { type: 'storage' } },
460 { binding: 2, visibility: GPUShaderStage.COMPUTE, buffer: { type: 'storage' } },
461 { binding: 3, visibility: GPUShaderStage.COMPUTE, buffer: { type: 'uniform', hasDynamicOffset: true, minBindingSize: 16 } },
462 ],
463 });
464 const plLayout = device.createPipelineLayout({ bindGroupLayouts: [bgl] });
466 const pipelines: GPUComputePipeline[] = [];
467 for (const code of sources) {
468 const module = device.createShaderModule({ code });
469 const info = await module.getCompilationInfo();
470 const errs = info.messages.filter((m) => m.type === 'error');
471 if (errs.length) {
472 throw new Error('WGSL: ' + errs.map((m) => `${m.lineNum}:${m.linePos} ${m.message}`).join('; '));
473 }
474 pipelines.push(device.createComputePipeline({ layout: plLayout, compute: { module, entryPoint: 'main' } }));
475 }
477 const mkBuf = (size: number, usage: number) => device.createBuffer({ size, usage });
478 const opsBuf = mkBuf(layout.opsData.byteLength, GPUBufferUsage.STORAGE | GPUBufferUsage.COPY_DST);
479 device.queue.writeBuffer(opsBuf, 0, layout.opsData);
480 const stBuf = mkBuf(layout.stLen * 4, GPUBufferUsage.STORAGE | GPUBufferUsage.COPY_SRC | GPUBufferUsage.COPY_DST);
481 const outBuf = mkBuf(2 * nsteps * 4, GPUBufferUsage.STORAGE | GPUBufferUsage.COPY_SRC);
482 const uniBuf = mkBuf(chunk * 256, GPUBufferUsage.UNIFORM | GPUBufferUsage.COPY_DST);
483 {
484 const uni = new Uint32Array(chunk * 64);
485 for (let s = 0; s < chunk; s++) uni[s * 64] = s;
486 device.queue.writeBuffer(uniBuf, 0, uni);
487 }
488 const readOut = mkBuf(2 * nsteps * 4, GPUBufferUsage.MAP_READ | GPUBufferUsage.COPY_DST);
489 const readSt = mkBuf(layout.stLen * 4, GPUBufferUsage.MAP_READ | GPUBufferUsage.COPY_DST);
491 const bindGroup = device.createBindGroup({
492 layout: bgl,
493 entries: [
494 { binding: 0, resource: { buffer: opsBuf } },
495 { binding: 1, resource: { buffer: stBuf } },
496 { binding: 2, resource: { buffer: outBuf } },
497 { binding: 3, resource: { buffer: uniBuf, offset: 0, size: 16 } },
498 ],
499 });
501 const t0 = performance.now();
502 const pending: Promise<undefined>[] = [];
503 let kbase = 0;
504 while (kbase < nsteps) {
505 const n = Math.min(chunk, nsteps - kbase);
506 const kb = new Float32Array(new ArrayBuffer(4));
507 kb[0] = kbase;
508 device.queue.writeBuffer(stBuf, 0, kb);
509 const enc = device.createCommandEncoder();
510 const pass = enc.beginComputePass();
511 for (let s = 0; s < n; s++) {
512 pass.setBindGroup(0, bindGroup, [s * 256]);
513 for (const p of pipelines) {
514 pass.setPipeline(p);
515 pass.dispatchWorkgroups(nwg);
516 }
517 }
518 pass.end();
519 device.queue.submit([enc.finish()]);
520 pending.push(device.queue.onSubmittedWorkDone());
521 if (pending.length >= 2) await pending.shift();
522 kbase += n;
523 if (onProgress) onProgress(kbase, nsteps);
524 }
525 while (pending.length) await pending.shift();
526 const gpuMs = performance.now() - t0;
528 const enc = device.createCommandEncoder();
529 enc.copyBufferToBuffer(outBuf, 0, readOut, 0, 2 * nsteps * 4);
530 enc.copyBufferToBuffer(stBuf, 0, readSt, 0, layout.stLen * 4);
531 device.queue.submit([enc.finish()]);
532 await readOut.mapAsync(GPUMapMode.READ);
533 await readSt.mapAsync(GPUMapMode.READ);
534 const outData = new Float32Array(readOut.getMappedRange().slice(0));
535 const state = new Float32Array(readSt.getMappedRange().slice(0));
536 readOut.unmap();
537 readSt.unmap();
538 for (const b of [opsBuf, stBuf, outBuf, uniBuf, readOut, readSt]) b.destroy();
540 const acc = outData.subarray(0, nsteps);
541 for (let i = 0; i < acc.length; i++) {
542 if (!Number.isFinite(acc[i])) {
543 throw new Error(
544 `cymbal2: simulation diverged at t = ${(i * su.ht).toFixed(4)} s ` +
545 '(reduce strike_velocity, or raise n_theta / oversample)',
546 );
547 }
548 }
550 return {
551 acc,
552 fhist: outData.subarray(nsteps, 2 * nsteps),
553 peakXi: state[12],
554 peakDisp: state[13],
555 contactVel: Array.from(state.subarray(layout.SC_CVEL, layout.SC_CVEL + su.K)),
556 strokesPlayed: Math.min(state[4] + (state[3] > 0 ? 1 : 0), su.K),
557 gpuMs,
558 };
559}