1// Web Worker that renders the cymbal2 stroke. The MATLAB implementation in
2// ../matlab remains the source of truth; what varies is who runs the time
3// loop. When the browser exposes WebGPU, the worker executes the unmodified
4// setup stage (grid, per-wavenumber inverses, mallet, stroke schedule)
5// through numbl, hands the operators to a WGSL port of cymbal2_sim.m
6// (./gpuSim.ts), and feeds the returned acceleration back through the
7// unmodified cymbal2_output.m; this is 20x or more faster than interpreting
8// the loop. Without WebGPU, or when cymbal2_sim.m has been edited in the app
9// (the WGSL mirrors the shipped step, so edits there must fall back), the
10// original interpreter path runs instead. Requests carry a quality tier
11// preview -> oversample 1, tail clamped to 1.5 s (the azimuthal grid is
12// kept: hard strokes need the full wavenumber band to carry
13// the cascade, and coarsening it destabilizes the integration)
14// full -> the parameters as given
15// and the UI shows elapsed time, plus integration progress on the GPU path.
16// A cheap linear-modes stage (cymbal2_modes per wavenumber) is cached on the
17// geometry so the spectrum plot can overlay the mode table immediately.
19import { executeCode, RTV, type RuntimeValue, type RuntimeTensor } from 'numbl';
20import { hashMFiles, DEFAULT_MFILES, type MFile } from '../mfiles';
21import { MODES_KEYS, type Cymbal2Params } from '../params';
22import { requestGpuDevice, runGpuSim, type GpuSetup } from './gpuSim';
24export interface RenderRequest {
25 type: 'render';
26 id: number;
27 params: Cymbal2Params;
28 quality: 'preview' | 'full';
29 mfiles: MFile[];
30}
32export interface RenderResponse {
33 type: 'result';
34 id: number;
35 quality: 'preview' | 'full';
36 y: Float32Array; // final signal, [-1, 1]
37 sampleRate: number;
38 forceT: Float32Array; // contact force history (max-pooled), s
39 forceN: Float32Array; // and its values, N
40 modeF: Float32Array; // linear mode frequencies, Hz
41 modeM: Float32Array; // nodal diameters of each
42 diag: {
43 strokesPlayed: number;
44 meanV: number;
45 contactVel: number;
46 peakForceN: number;
47 contactMs: number;
48 contactEvents: number;
49 peakXiMm: number;
50 peakDispMm: number;
51 peakDispThickness: number;
52 peakAccel: number;
53 strikeRadiusM: number;
54 pickupRadiusM: number;
55 RsM: number;
56 nsteps: number;
57 simRate: number;
58 nDropped: number;
59 durationS: number;
60 renderMs: number;
61 oversample: number;
62 nTheta: number;
63 backend: 'webgpu' | 'cpu';
64 warnings: string[];
65 };
66}
68export interface ProgressMessage {
69 type: 'progress';
70 id: number;
71 stage: 'modes' | 'integrate' | 'output';
72 frac: number;
73}
75export interface ErrorMessage {
76 type: 'error';
77 id: number;
78 message: string;
79 file?: string;
80 line?: number;
81}
83interface ModesEntry {
84 modeF: Float32Array;
85 modeM: Float32Array;
86}
88const modesCache = new Map<string, ModesEntry>();
89const MODES_CACHE_MAX = 4;
91function lruGet<T>(cache: Map<string, T>, key: string): T | undefined {
92 const v = cache.get(key);
93 if (v !== undefined) {
94 cache.delete(key);
95 cache.set(key, v);
96 }
97 return v;
98}
100function lruSet<T>(cache: Map<string, T>, key: string, value: T, max: number): void {
101 cache.delete(key);
102 cache.set(key, value);
103 while (cache.size > max) {
104 const oldest = cache.keys().next().value as string;
105 cache.delete(oldest);
106 }
107}
109function keyFor(params: Cymbal2Params, keys: readonly (keyof Cymbal2Params)[]): string {
110 return keys.map((k) => String(params[k])).join('|');
111}
113// MATLAB source lines assigning every parameter onto a struct P.
114function pLiteral(params: Cymbal2Params): string {
115 const lines: string[] = [];
116 for (const [k, v] of Object.entries(params)) {
117 if (k === 'bit_depth') continue;
118 if (typeof v === 'number') lines.push(`P.${k} = ${v};`);
119 else if (typeof v === 'boolean') lines.push(`P.${k} = ${v ? 1 : 0};`);
120 }
121 return lines.join('\n');
122}
124function runStage(mfiles: MFile[], source: string, inputs?: Record<string, RuntimeValue>) {
125 return executeCode(
126 source,
127 { initialVariableValues: inputs, displayResults: false },
128 mfiles,
129 'stage.m',
130 );
131}
133function asScalar(v: RuntimeValue): number {
134 if (typeof v === 'number') return v;
135 const t = v as RuntimeTensor;
136 if (t && t.kind === 'tensor' && t.data.length === 1) return t.data[0];
137 if (typeof v === 'boolean') return v ? 1 : 0;
138 throw new Error('expected scalar result from numbl stage');
139}
141function toF32(v: RuntimeValue): Float32Array {
142 const t = v as RuntimeTensor;
143 if (!t || t.kind !== 'tensor') throw new Error('expected tensor result from numbl stage');
144 const out = new Float32Array(t.data.length);
145 for (let i = 0; i < t.data.length; i++) out[i] = t.data[i];
146 return out;
147}
149// Max-pool the contact force history to about 2000 points for the plot.
150function poolForce(fhist: Float32Array, simRate: number): { tpool: Float32Array; fpool: Float32Array } {
151 const nsim = fhist.length;
152 const np0 = Math.min(2000, nsim);
153 const blk = Math.ceil(nsim / np0);
154 const np = Math.ceil(nsim / blk);
155 const fpool = new Float32Array(np);
156 const tpool = new Float32Array(np);
157 for (let b = 0; b < np; b++) {
158 let m = 0;
159 const end = Math.min(nsim, (b + 1) * blk);
160 for (let i = b * blk; i < end; i++) m = Math.max(m, fhist[i]);
161 fpool[b] = m;
162 tpool[b] = (b * blk + blk / 2) / simRate;
163 }
164 return { tpool, fpool };
165}
167const SKIP_WARNING = /skipping (\S+\.m) \(syntax error at line (\d+)\)/;
169self.onmessage = (ev: MessageEvent<RenderRequest>) => {
170 const msg = ev.data;
171 if (msg.type !== 'render') return;
173 const warnings: string[] = [];
174 const origWarn = console.warn;
175 console.warn = (...args: unknown[]) => {
176 warnings.push(args.map(String).join(' '));
177 };
179 render(msg)
180 .catch((err) => {
181 const e = err as { message?: string; file?: string; line?: number };
182 let m: ErrorMessage = {
183 type: 'error',
184 id: msg.id,
185 message: e?.message ? String(e.message) : String(err),
186 file: typeof e?.file === 'string' ? e.file : undefined,
187 line: typeof e?.line === 'number' ? e.line : undefined,
188 };
189 for (const w of warnings) {
190 const hit = SKIP_WARNING.exec(w);
191 if (hit) {
192 m = { type: 'error', id: msg.id, message: 'syntax error', file: hit[1], line: Number(hit[2]) };
193 break;
194 }
195 }
196 (self as unknown as Worker).postMessage(m);
197 })
198 .finally(() => {
199 console.warn = origWarn;
200 });
201};
203// Everything the response assembly needs, produced by either backend.
204interface SimOutcome {
205 y32: Float32Array;
206 peakAccel: number;
207 nDropped: number;
208 tpool: Float32Array;
209 fpool: Float32Array;
210 backend: 'webgpu' | 'cpu';
211 d: {
212 strokesPlayed: number;
213 meanV: number;
214 contactVel: number;
215 peakForceN: number;
216 contactMs: number;
217 contactEvents: number;
218 peakXiMm: number;
219 peakDispM: number;
220 peakDispThickness: number;
221 nsteps: number;
222 simRate: number;
223 strikeRadiusM: number;
224 pickupRadiusM: number;
225 RsM: number;
226 };
227}
229// The stroke-schedule block from cymbal2_sim.m, verbatim, so the GPU render
230// integrates exactly the schedule the interpreter would.
231const SCHEDULE_SRC = `
232K = max(1, round(P.strike_count));
233tstr = P.pre_delay_ms * 1e-3 + (0:K - 1) / max(1e-3, P.strike_rate_hz);
234vstr = zeros(1, K);
235state = P.seed * 2654435761;
236for i = 1:K
237 state = mod(state * 1664525 + 1013904223, 4294967296);
238 x = state / 4294967296;
239 vstr(i) = max(0.05, P.strike_velocity * (1 + P.velocity_spread * (2 * x - 1)));
240end
241gap = P.strike_gap_m;
242if gap <= 0
243 gap = 1.05 * max(vstr)^2 / (2 * P.strike_accel) + 0.005;
244end
245travel = gap ./ vstr + vstr / (2 * P.strike_accel);
246tlaunch = tstr - travel;
247off = -min(0, min(tlaunch));
248tlaunch = tlaunch + off;
249duration = off + (K - 1) / max(1e-3, P.strike_rate_hz) + P.tail_s;
250nsteps = floor(duration * P.sample_rate * P.oversample);
251`;
253async function renderOnGpu(
254 mfiles: MFile[],
255 params: Cymbal2Params,
256 onProgress: (frac: number) => void,
257): Promise<SimOutcome> {
258 // A WebGPU device first: if there is none, fail before spending seconds on
259 // the setup stage.
260 const device = await requestGpuDevice();
261 try {
262 const r = runStage(mfiles, `
263${pLiteral(params)}
264S = cymbal2_system(P);
265MB = cymbal2_mallet(P);
266G = S.G;
267o_D1 = G.D1; o_D2 = G.D2;
268o_F = G.F; o_Gi = G.Gi;
269o_Rd = G.Rd; o_R2d = G.R2d;
270o_AiB = S.AiB; o_AiC = S.AiC; o_Ai = S.Ai; o_Fi = S.Fi;
271o_ASiB = MB.ASiB; o_ASiC = MB.ASiC;
272o_vM = MB.vM; o_vR = MB.vR; o_vK = MB.vK; o_tipcol = MB.tipcol;
273${SCHEDULE_SRC}
274o_tlaunch = tlaunch; o_vstr = vstr;
275c_nrt = G.nrt; c_Nr = G.Nr; c_nth = G.nth; c_hr = G.hr; c_hth = G.hth;
276c_ht = S.ht; c_Eh = S.Eh; c_Rs = S.Rs; c_rhoh = S.rhoh;
277c_fsc = S.forceScale; c_nlr0 = S.eomRow0; c_nDropped = S.nDropped;
278c_ls = S.ls; c_js = S.js; c_lp = S.lp; c_jp = S.jp;
279c_strikeR = S.strikeRadiusM; c_pickupR = S.pickupRadiusM;
280c_tipi = MB.tipi; c_nfree = MB.nfree; c_ps = MB.ps;
281c_holdSteps = max(1, round(P.strike_hold_ms * 1e-3 / S.ht));
282c_missDepth = gap + 0.06;
283c_gap = gap; c_K = K; c_nsteps = nsteps;
284c_simRate = P.sample_rate * P.oversample;
285`);
286 const v = r.variableValues;
287 const s = (n: string) => asScalar(v[n]);
288 const t = (n: string) => toF32(v[n]);
289 const setup: GpuSetup = {
290 ops: {
291 D1: t('o_D1'), D2: t('o_D2'),
292 F: t('o_F'), Gi: t('o_Gi'),
293 Rd: t('o_Rd'), R2d: t('o_R2d'),
294 AiB: t('o_AiB'), AiC: t('o_AiC'), Ai: t('o_Ai'), Fi: t('o_Fi'),
295 ASiB: t('o_ASiB'), ASiC: t('o_ASiC'),
296 vM: t('o_vM'), vR: t('o_vR'), vK: t('o_vK'), tipcol: t('o_tipcol'),
297 },
298 tlaunch: Array.from((v.o_tlaunch as RuntimeTensor).data),
299 vstr: Array.from((v.o_vstr as RuntimeTensor).data),
300 nrt: s('c_nrt'), Nr: s('c_Nr'), nth: s('c_nth'),
301 hr: s('c_hr'), hth: s('c_hth'),
302 ht: s('c_ht'), Eh: s('c_Eh'), Rs: s('c_Rs'), rhoh: s('c_rhoh'),
303 fsc: s('c_fsc'), nlr0: s('c_nlr0'),
304 ls: s('c_ls'), js: s('c_js'), lp: s('c_lp'), jp: s('c_jp'),
305 tipi: s('c_tipi'), nfree: s('c_nfree'), ps: s('c_ps'),
306 KM: params.head_stiffness, alM: params.head_alpha,
307 aH: params.strike_accel, aHr: params.retract_accel,
308 holdSteps: s('c_holdSteps'), missDepth: s('c_missDepth'),
309 gap: s('c_gap'), K: s('c_K'), nsteps: s('c_nsteps'),
310 };
311 const simRate = s('c_simRate');
313 const sim = await runGpuSim(device, setup, (done, n) => onProgress(done / n));
315 // The output stage through the unmodified cymbal2_output.m.
316 const accT = RTV.tensor(Float64Array.from(sim.acc), [1, sim.acc.length]);
317 const ro = runStage(mfiles, `
318${pLiteral(params)}
319[y, peakAccel] = cymbal2_output(P, acc);
320`, { acc: accT });
322 // Contact bookkeeping from the force history, as cymbal2_sim does it.
323 let peakF = 0, contactSteps = 0, events = 0, prev = 0;
324 for (let i = 0; i < sim.fhist.length; i++) {
325 const f = sim.fhist[i];
326 if (f > peakF) peakF = f;
327 if (f > 0) { contactSteps++; if (prev === 0) events++; }
328 prev = f;
329 }
330 const { tpool, fpool } = poolForce(sim.fhist, simRate);
331 const meanV = setup.vstr.reduce((a, b) => a + b, 0) / setup.vstr.length;
333 return {
334 y32: toF32(ro.variableValues.y),
335 peakAccel: asScalar(ro.variableValues.peakAccel),
336 nDropped: s('c_nDropped'),
337 tpool,
338 fpool,
339 backend: 'webgpu',
340 d: {
341 strokesPlayed: sim.strokesPlayed,
342 meanV,
343 contactVel: sim.contactVel[0] ?? 0,
344 peakForceN: peakF,
345 contactMs: (1000 * contactSteps) / simRate,
346 contactEvents: events,
347 peakXiMm: sim.peakXi * 1000,
348 peakDispM: sim.peakDisp,
349 peakDispThickness: sim.peakDisp / params.thickness_m,
350 nsteps: setup.nsteps,
351 simRate,
352 strikeRadiusM: s('c_strikeR'),
353 pickupRadiusM: s('c_pickupR'),
354 RsM: setup.Rs,
355 },
356 };
357 } finally {
358 device.destroy();
359 }
360}
362function renderOnCpu(mfiles: MFile[], params: Cymbal2Params): SimOutcome {
363 const r = runStage(mfiles, `
364${pLiteral(params)}
365S = cymbal2_system(P);
366MB = cymbal2_mallet(P);
367[acc, dinfo, fhist] = cymbal2_sim(P, S, MB);
368[y, peakAccel] = cymbal2_output(P, acc);
369% max-pool the contact force history to about 2000 points for the plot
370nsim = length(fhist);
371np = min(2000, nsim);
372blk = ceil(nsim / np);
373np = ceil(nsim / blk);
374fpool = max(reshape([fhist, zeros(1, np * blk - nsim)], blk, np), [], 1);
375tpool = ((0:np - 1) * blk + blk / 2) / dinfo.simRate;
376d_strokes = dinfo.strokesPlayed;
377d_meanV = dinfo.meanV;
378d_cvel = dinfo.contactVel(1);
379d_peakF = dinfo.peakForceN;
380d_contactMs = dinfo.contactMs;
381d_events = dinfo.contactEvents;
382d_xi = dinfo.peakXiMm;
383d_disp = dinfo.peakDispM;
384d_dispH = dinfo.peakDispThickness;
385d_nsteps = dinfo.nsteps;
386d_simRate = dinfo.simRate;
387d_dropped = S.nDropped;
388d_strikeR = S.strikeRadiusM;
389d_pickupR = S.pickupRadiusM;
390d_Rs = S.Rs;
391`);
392 const v = r.variableValues;
393 return {
394 y32: toF32(v.y),
395 peakAccel: asScalar(v.peakAccel),
396 nDropped: asScalar(v.d_dropped),
397 tpool: toF32(v.tpool),
398 fpool: toF32(v.fpool),
399 backend: 'cpu',
400 d: {
401 strokesPlayed: asScalar(v.d_strokes),
402 meanV: asScalar(v.d_meanV),
403 contactVel: asScalar(v.d_cvel),
404 peakForceN: asScalar(v.d_peakF),
405 contactMs: asScalar(v.d_contactMs),
406 contactEvents: asScalar(v.d_events),
407 peakXiMm: asScalar(v.d_xi),
408 peakDispM: asScalar(v.d_disp),
409 peakDispThickness: asScalar(v.d_dispH),
410 nsteps: asScalar(v.d_nsteps),
411 simRate: asScalar(v.d_simRate),
412 strikeRadiusM: asScalar(v.d_strikeR),
413 pickupRadiusM: asScalar(v.d_pickupR),
414 RsM: asScalar(v.d_Rs),
415 },
416 };
417}
419async function render(msg: RenderRequest): Promise<void> {
420 const { id, mfiles, quality } = msg;
421 const params = { ...msg.params };
422 if (quality === 'preview') {
423 params.oversample = 1;
424 params.tail_s = Math.min(params.tail_s, 1.5);
425 }
427 const post = (m: ProgressMessage) => (self as unknown as Worker).postMessage(m);
428 const codeKey = hashMFiles(mfiles);
430 // ---- linear modes, cheap, cached on geometry ----
431 const modesKey = codeKey + '|' + keyFor(params, MODES_KEYS);
432 let modesEntry = lruGet(modesCache, modesKey);
433 if (!modesEntry) {
434 post({ type: 'progress', id, stage: 'modes', frac: 0.05 });
435 const r = runStage(mfiles, `
436${pLiteral(params)}
437mF = [];
438mM = [];
439for mw = 0:min(10, floor(P.n_theta / 2))
440 f = cymbal2_modes(P, mw);
441 f = f(f > 0.5 & f < 4000);
442 nk = min(4, length(f));
443 mF = [mF; f(1:nk)];
444 mM = [mM; mw * ones(nk, 1)];
445end
446`);
447 modesEntry = {
448 modeF: toF32(r.variableValues.mF),
449 modeM: toF32(r.variableValues.mM),
450 };
451 lruSet(modesCache, modesKey, modesEntry, MODES_CACHE_MAX);
452 }
454 // ---- the time integration ----
455 post({ type: 'progress', id, stage: 'integrate', frac: 0.1 });
456 const t0 = performance.now();
457 const warnings: string[] = [];
459 // The WGSL mirrors the shipped cymbal2_sim.m, so an edited step falls back
460 // to the interpreter, where the edit is live. Edits to the setup and output
461 // stages flow through numbl on both paths.
462 const simSource = mfiles.find((f) => f.name === 'cymbal2_sim.m')?.source;
463 const simShipped = DEFAULT_MFILES.find((f) => f.name === 'cymbal2_sim.m')?.source;
464 const simEdited = simSource !== simShipped;
466 let out: SimOutcome | null = null;
467 if (!simEdited) {
468 try {
469 out = await renderOnGpu(mfiles, params, (frac) =>
470 post({ type: 'progress', id, stage: 'integrate', frac: 0.1 + 0.85 * frac }),
471 );
472 } catch (err) {
473 const detail = err instanceof Error ? err.message : String(err);
474 if (/diverged/.test(detail)) throw err;
475 // console.warn is intercepted during a render; use info so the reason
476 // is visible in the devtools console.
477 console.info(`cymbal2: WebGPU unavailable (${detail}); rendering on the CPU`);
478 }
479 }
480 if (!out) {
481 out = renderOnCpu(mfiles, params);
482 if (simEdited) {
483 // Not a fault: just say why the render took the slow path.
484 warnings.push('cymbal2_sim.m differs from the shipped source, so the render used the interpreter (the GPU path mirrors the shipped step)');
485 }
486 }
487 const renderMs = performance.now() - t0;
489 post({ type: 'progress', id, stage: 'output', frac: 0.97 });
490 const { y32, peakAccel, nDropped } = out;
492 if (nDropped > 0) {
493 warnings.push(
494 `${nDropped} azimuthal columns were truncated: the radial grid cannot resolve their inner boundary layer; raise n_radial`,
495 );
496 }
497 if (!params.normalize) {
498 let peak = 0;
499 for (let i = 0; i < y32.length; i++) peak = Math.max(peak, Math.abs(y32[i]));
500 if (peak > 1.0) {
501 warnings.push(
502 `clipping: peak acceleration ${peakAccel.toExponential(2)} m/s² exceeds full_scale_accel`,
503 );
504 }
505 }
507 const resp: RenderResponse = {
508 type: 'result',
509 id,
510 quality,
511 y: y32,
512 sampleRate: params.sample_rate,
513 forceT: out.tpool,
514 forceN: out.fpool,
515 modeF: modesEntry.modeF.slice(),
516 modeM: modesEntry.modeM.slice(),
517 diag: {
518 strokesPlayed: out.d.strokesPlayed,
519 meanV: out.d.meanV,
520 contactVel: out.d.contactVel,
521 peakForceN: out.d.peakForceN,
522 contactMs: out.d.contactMs,
523 contactEvents: out.d.contactEvents,
524 peakXiMm: out.d.peakXiMm,
525 peakDispMm: out.d.peakDispM * 1000,
526 peakDispThickness: out.d.peakDispThickness,
527 peakAccel,
528 strikeRadiusM: out.d.strikeRadiusM,
529 pickupRadiusM: out.d.pickupRadiusM,
530 RsM: out.d.RsM,
531 nsteps: out.d.nsteps,
532 simRate: out.d.simRate,
533 nDropped,
534 durationS: y32.length / params.sample_rate,
535 renderMs,
536 oversample: params.oversample,
537 nTheta: params.n_theta,
538 backend: out.backend,
539 warnings,
540 },
541 };
542 (self as unknown as Worker).postMessage(resp, [
543 resp.y.buffer,
544 resp.forceT.buffer,
545 resp.forceN.buffer,
546 resp.modeF.buffer,
547 resp.modeM.buffer,
548 ]);
549}