Blame webapp/src/worker/renderWorker.ts at 4ddcc1b5d8b0781a4c34c4796d2470d8b9c9a0a4

4ddcc1b Cymbal simulator: shell-of-revolution modal model with a nonlinear energy cascade, numbl MATLAB reference, and interactive web appJeremy Magland 1 month ago1// Web Worker that renders the cymbal stroke by executing the MATLAB

2// implementation in ../matlab through numbl. The stages run as separate

3// executeCode calls so intermediates can be cached and cheap parameter

4// changes skip the expensive ones:

5// solve tier -> cymbal_basis (one shell eigensolve per nodal-diameter

6// order; the expensive stage, cached on geometry, material,

7// damping and discretization parameters)

8// radiate tier -> cymbal_strikes + cymbal_cascade + cymbal_radiate (cached

9// on the stick, the gesture, the cascade, the listener and

10// the synthesis rate)

11// output tier -> cymbal_output (always cheap, always re-run)

12// Tensors cross the boundary as numbl RuntimeValues (column-major Float64).

14import { executeCode, type RuntimeValue, type RuntimeTensor } from 'numbl';

15import { hashMFiles, type MFile } from '../mfiles';

16import { SOLVE_KEYS, RADIATE_KEYS, type CymbalParams } from '../params';

18export interface RenderRequest {

19 type: 'render';

20 id: number;

21 params: CymbalParams;

22 quality: 'preview' | 'full';

23 mfiles: MFile[];

24}

26export interface RenderResponse {

27 type: 'result';

28 id: number;

29 quality: 'preview' | 'full';

30 y: Float32Array; // final signal, [-1, 1]

31 sampleRate: number;

32 strikeT: Float32Array; // stroke times, s

33 strikeV: Float32Array; // stroke velocities, m/s

34 strikeS: Float32Array; // stroke arclength fractions

35 modeF: Float32Array; // mode frequencies, Hz

36 modeZeta: Float32Array; // damping ratio of each mode

37 modeZetaRad: Float32Array; // radiation part of the damping

38 modeSigma: Float32Array; // radiation efficiency of each mode

39 modeN: Float32Array; // nodal diameters

40 nodeR: Float32Array; // meridian nodes, m

41 nodeZ: Float32Array;

42 shapeW: Float32Array; // normal displacement per mode, nn x n col-major

43 cascadeDb: Float32Array; // energy grid, nBands x nFrames col-major, dB

44 cascadeBands: number;

45 cascadeFrames: number;

46 cascadeFLo: number;

47 cascadeFHi: number;

48 cascadeTEnd: number;

49 diag: {

50 modesKept: number;

51 maxN: number;

52 strikeCount: number;

53 meanV: number;

54 maxV: number;

55 meanTauUs: number;

56 peakForceN: number;

57 durationS: number;

58 peakPa: number;

59 peakSplDb: number;

60 directPeakPa: number;

61 listenerDistM: number;

62 peakWrmsUm: number;

63 peakStrength: number;

64 cascadeLostFrac: number;

65 plateMassKg: number;

66 warnings: string[];

67 solveMs: number; // 0 when the cached eigensolve was reused

68 radiateMs: number; // 0 when the cached radiation was reused

69 solveCached: boolean;

70 radiateCached: boolean;

71 oversample: number;

72 };

73}

75export interface ProgressMessage {

76 type: 'progress';

77 id: number;

78 stage: 'solve' | 'radiate' | 'output';

79 frac: number;

80}

82// Raised by a fault in the .m code (the user can edit it), so it is reported

83// rather than left to stall the request.

84export interface ErrorMessage {

85 type: 'error';

86 id: number;

87 message: string;

88 file?: string;

89 line?: number;

90}

92interface SolveEntry {

93 vars: Record<string, RuntimeValue>;

94 modeF: Float32Array;

95 modeZeta: Float32Array;

96 modeZetaRad: Float32Array;

97 modeSigma: Float32Array;

98 modeN: Float32Array;

99 nodeR: Float32Array;

100 nodeZ: Float32Array;

101 shapeW: Float32Array;

102 maxN: number;

103 plateMassKg: number;

104}

105

106interface RadiateEntry {

107 p: RuntimeValue;

108 strikeT: Float32Array;

109 strikeV: Float32Array;

110 strikeS: Float32Array;

111 cascadeDb: Float32Array;

112 cascadeBands: number;

113 cascadeFrames: number;

114 cascadeTEnd: number;

115 strikeCount: number;

116 meanV: number;

117 maxV: number;

118 meanTauUs: number;

119 peakForceN: number;

120 durationS: number;

121 directPeakPa: number;

122 listenerDistM: number;

123 listenerInside: boolean;

124 peakWrmsUm: number;

125 peakStrength: number;

126 cascadeLostFrac: number;

127 oversample: number;

128}

129

130// Small LRU caches. The solve cache does not depend on the oversample rate,

131// so a preview render right after a full one reuses the same eigensolves.

132const solveCache = new Map<string, SolveEntry>();

133const radiateCache = new Map<string, RadiateEntry>();

134const SOLVE_CACHE_MAX = 2;

135const RADIATE_CACHE_MAX = 3;

136

137function lruGet<T>(cache: Map<string, T>, key: string): T | undefined {

138 const v = cache.get(key);

139 if (v !== undefined) {

140 cache.delete(key);

141 cache.set(key, v);

142 }

143 return v;

144}

145

146function lruSet<T>(cache: Map<string, T>, key: string, value: T, max: number): void {

147 cache.delete(key);

148 cache.set(key, value);

149 while (cache.size > max) {

150 const oldest = cache.keys().next().value as string;

151 cache.delete(oldest);

152 }

153}

154

155function keyFor(params: CymbalParams, keys: readonly (keyof CymbalParams)[]): string {

156 return keys.map((k) => String(params[k])).join('|');

157}

158

159// MATLAB source lines assigning every parameter onto a struct P.

160function pLiteral(params: CymbalParams): string {

161 const lines: string[] = [];

162 for (const [k, v] of Object.entries(params)) {

163 if (typeof v === 'number') lines.push(`P.${k} = ${v};`);

164 else if (typeof v === 'boolean') lines.push(`P.${k} = ${v ? 1 : 0};`);

165 }

166 return lines.join('\n');

167}

168

169function runStage(mfiles: MFile[], source: string, inputs?: Record<string, RuntimeValue>) {

170 return executeCode(

171 source,

172 { initialVariableValues: inputs, displayResults: false },

173 mfiles,

174 'stage.m',

175 );

176}

177

178function asTensor(v: RuntimeValue): RuntimeTensor {

179 const t = v as RuntimeTensor;

180 if (!t || t.kind !== 'tensor') throw new Error('expected tensor result from numbl stage');

181 return t;

182}

183

184function asScalar(v: RuntimeValue): number {

185 if (typeof v === 'number') return v;

186 const t = v as RuntimeTensor;

187 if (t && t.kind === 'tensor' && t.data.length === 1) return t.data[0];

188 if (typeof v === 'boolean') return v ? 1 : 0;

189 throw new Error('expected scalar result from numbl stage');

190}

191

192function toF32(v: RuntimeValue): Float32Array {

193 const t = asTensor(v);

194 const out = new Float32Array(t.data.length);

195 for (let i = 0; i < t.data.length; i++) out[i] = t.data[i];

196 return out;

197}

198

199// The basis fields handed from the solve stage into the radiate stage, in the

200// order the reassembly script expects.

201const BASIS_FIELDS = [

202 'f', 'omega', 'zeta', 'zetaRad', 'zetaMount', 'sigma', 'ndia', 'ldip',

203 'NU', 'NW', 'NT',

204 'rnode', 'znode', 'snode', 'drdsnode', 'dzdsnode', 'thknode', 'csnode', 'ctnode',

205 'surfR', 'surfZ', 'surfDr', 'surfDz', 'surfLen', 'SW',

206 'S', 'area', 'mass', 'thkMean',

207];

208

209const CASCADE_BANDS = 48;

210const CASCADE_FRAMES = 240;

211

212const SKIP_WARNING = /skipping (\S+\.m) $syntax error at line (\d+)$/;

213

214self.onmessage = (ev: MessageEvent<RenderRequest>) => {

215 const msg = ev.data;

216 if (msg.type !== 'render') return;

217

218 // numbl reports a file it could not parse through console.warn and then

219 // skips it, so the failure surfaces at the call site as "Undefined

220 // function". Capture the warning to report the actual syntax error instead.

221 const warnings: string[] = [];

222 const origWarn = console.warn;

223 console.warn = (...args: unknown[]) => {

224 warnings.push(args.map(String).join(' '));

225 };

226

227 try {

228 render(msg);

229 } catch (err) {

230 const e = err as { message?: string; file?: string; line?: number };

231 let m: ErrorMessage = {

232 type: 'error',

233 id: msg.id,

234 message: e?.message ? String(e.message) : String(err),

235 file: typeof e?.file === 'string' ? e.file : undefined,

236 line: typeof e?.line === 'number' ? e.line : undefined,

237 };

238 for (const w of warnings) {

239 const hit = SKIP_WARNING.exec(w);

240 if (hit) {

241 m = { type: 'error', id: msg.id, message: 'syntax error', file: hit[1], line: Number(hit[2]) };

242 break;

243 }

244 }

245 (self as unknown as Worker).postMessage(m);

246 } finally {

247 console.warn = origWarn;

248 }

249};

250

251function render(msg: RenderRequest): void {

252 const { id, params, mfiles } = msg;

253 let quality = msg.quality;

254

255 // The .m sources are part of the cache identity: editing the code must

256 // invalidate everything computed from the previous version.

257 const codeKey = hashMFiles(mfiles);

258 const solveKey = codeKey + '|' + keyFor(params, SOLVE_KEYS);

259

260 const fullOversample = params.oversample;

261 const radiateParamsKey = solveKey + '||' + keyFor(params, RADIATE_KEYS);

262 if (quality === 'preview' && radiateCache.has(radiateParamsKey + `|os=${fullOversample}`)) {

263 quality = 'full'; // full-rate radiation already cached

264 }

265 const oversample = quality === 'preview' ? Math.min(1, params.oversample) : params.oversample;

266 const fsHi = params.sample_rate * oversample;

267

268 const post = (m: ProgressMessage) => (self as unknown as Worker).postMessage(m);

269

270 // ---- solve stage: one shell eigensolve per nodal-diameter order ----

271 let solveMs = 0;

272 let solveCached = true;

273 let solveEntry = lruGet(solveCache, solveKey);

274 if (!solveEntry) {

275 solveCached = false;

276 post({ type: 'progress', id, stage: 'solve', frac: 0.12 });

277 const t0 = performance.now();

278 const r = runStage(mfiles, `

279${pLiteral(params)}

280basis = cymbal_basis(P);

281${BASIS_FIELDS.map((f) => `b_${f} = basis.${f};`).join('\n')}

282maxN = max(basis.ndia);

283`);

284 solveMs = performance.now() - t0;

285

286 const vars: Record<string, RuntimeValue> = {};

287 for (const f of BASIS_FIELDS) vars[`b_${f}`] = r.variableValues[`b_${f}`];

288 solveEntry = {

289 vars,

290 modeF: toF32(r.variableValues.b_f),

291 modeZeta: toF32(r.variableValues.b_zeta),

292 modeZetaRad: toF32(r.variableValues.b_zetaRad),

293 modeSigma: toF32(r.variableValues.b_sigma),

294 modeN: toF32(r.variableValues.b_ndia),

295 nodeR: toF32(r.variableValues.b_rnode),

296 nodeZ: toF32(r.variableValues.b_znode),

297 shapeW: toF32(r.variableValues.b_NW),

298 maxN: asScalar(r.variableValues.maxN),

299 plateMassKg: asScalar(r.variableValues.b_mass),

300 };

301 lruSet(solveCache, solveKey, solveEntry, SOLVE_CACHE_MAX);

302 }

303

304 // ---- radiate stage: strikes + cascade + modal synthesis ----

305 const radiateKey = radiateParamsKey + `|os=${oversample}`;

306 let radiateMs = 0;

307 let radiateCached = true;

308 let radEntry = lruGet(radiateCache, radiateKey);

309 if (!radEntry) {

310 radiateCached = false;

311 post({ type: 'progress', id, stage: 'radiate', frac: 0.5 });

312 const t0 = performance.now();

313 const r = runStage(mfiles, `

314${pLiteral(params)}

315${BASIS_FIELDS.map((f) => `basis.${f} = b_${f};`).join('\n')}

316basis.n = length(basis.f);

317[tstr, vstr, jstr, taustr, fstr, sstr, sinfo] = cymbal_strikes(P, basis);

318[env, tframe, DV, kf, cinfo] = cymbal_cascade(P, basis, tstr, jstr, taustr, sstr, sinfo.durationS);

319[p, radinfo] = cymbal_radiate(P, basis, env, tframe, DV, kf, tstr, taustr, fstr, sstr, ${fsHi});

320% collapse the modal energies onto a log-frequency by time grid for the plot

321nb = ${CASCADE_BANDS};

322nt = ${CASCADE_FRAMES};

323fLo = max(1, min(basis.f));

324fHi = max(fLo * 1.2, max(basis.f));

325bandOf = min(nb, max(1, 1 + floor(nb * (log(basis.f) - log(fLo)) / (log(fHi) - log(fLo) + 1e-12))));

326nf = size(env, 2);

327grid = zeros(nb, nt);

328for q = 1:nt

329 c0 = 1 + floor((q - 1) * nf / nt);

330 c1 = max(c0, floor(q * nf / nt));

331 e = sum(env(:, c0:c1).^2, 2) / (c1 - c0 + 1);

332 grid(:, q) = accumarray(bandOf, e, [nb, 1]);

333end

334cascadeDb = 10 * log10(grid + 1e-30);

335cascadeTEnd = tframe(length(tframe));

336strikeCount = sinfo.count;

337meanV = sinfo.meanV;

338maxV = sinfo.maxV;

339meanTauUs = sinfo.meanTauUs;

340peakForceN = sinfo.peakForceN;

341durationS = sinfo.durationS;

342directPeakPa = radinfo.directPeakPa;

343listenerDistM = radinfo.listenerDistM;

344listenerInside = radinfo.listenerInside;

345peakWrmsUm = cinfo.peakWrmsUm;

346peakStrength = cinfo.peakStrength;

347cascadeLostFrac = cinfo.cascadeLostJ / max(1e-30, cinfo.injectedJ);

348`, solveEntry.vars);

349 radiateMs = performance.now() - t0;

350 radEntry = {

351 p: r.variableValues.p,

352 strikeT: toF32(r.variableValues.tstr),

353 strikeV: toF32(r.variableValues.vstr),

354 strikeS: toF32(r.variableValues.sstr),

355 cascadeDb: toF32(r.variableValues.cascadeDb),

356 cascadeBands: CASCADE_BANDS,

357 cascadeFrames: CASCADE_FRAMES,

358 cascadeTEnd: asScalar(r.variableValues.cascadeTEnd),

359 strikeCount: asScalar(r.variableValues.strikeCount),

360 meanV: asScalar(r.variableValues.meanV),

361 maxV: asScalar(r.variableValues.maxV),

362 meanTauUs: asScalar(r.variableValues.meanTauUs),

363 peakForceN: asScalar(r.variableValues.peakForceN),

364 durationS: asScalar(r.variableValues.durationS),

365 directPeakPa: asScalar(r.variableValues.directPeakPa),

366 listenerDistM: asScalar(r.variableValues.listenerDistM),

367 listenerInside: asScalar(r.variableValues.listenerInside) !== 0,

368 peakWrmsUm: asScalar(r.variableValues.peakWrmsUm),

369 peakStrength: asScalar(r.variableValues.peakStrength),

370 cascadeLostFrac: asScalar(r.variableValues.cascadeLostFrac),

371 oversample,

372 };

373 lruSet(radiateCache, radiateKey, radEntry, RADIATE_CACHE_MAX);

374 }

375

376 // ---- output stage ----

377 post({ type: 'progress', id, stage: 'output', frac: 0.9 });

378 const rOut = runStage(mfiles, `

379${pLiteral(params)}

380[y, peakPa] = cymbal_output(P, p, ${radEntry.oversample});

381`, { p: radEntry.p });

382

383 const y32 = toF32(rOut.variableValues.y);

384 const peakPa = asScalar(rOut.variableValues.peakPa);

385

386 const warnings: string[] = [];

387 if (radEntry.listenerInside) {

388 warnings.push(

389 'the listener is within 2 cm of the plate; the dipole field is evaluated at a clamped minimum distance',

390 );

391 }

392 if (!params.normalize) {

393 let peak = 0;

394 for (let i = 0; i < y32.length; i++) peak = Math.max(peak, Math.abs(y32[i]));

395 if (peak > 1.0) {

396 warnings.push(

397 `clipping: peak is ${peakPa.toFixed(2)} Pa but full_scale_pa is ${params.full_scale_pa.toFixed(2)}`,

398 );

399 }

400 }

401 if (params.n_azimuth < 2.5 * solveEntry.maxN) {

402 warnings.push(

403 `only ${params.n_azimuth} azimuth samples for ${solveEntry.maxN} nodal diameters; the listener transfer of the high-order modes is aliased`,

404 );

405 }

406 if (radEntry.cascadeLostFrac > 0.25) {

407 warnings.push(

408 `${(100 * radEntry.cascadeLostFrac).toFixed(0)}% of the struck energy left the top of the retained band; raise max_mode_hz or lower cascade_leak`,

409 );

410 }

411 const topMode = solveEntry.modeF[solveEntry.modeF.length - 1];

412 if (topMode < 0.6 * params.max_mode_hz) {

413 warnings.push(

414 `the highest retained mode is only ${topMode.toFixed(0)} Hz; raise max nodal diameters to fill the band`,

415 );

416 }

417

418 const resp: RenderResponse = {

419 type: 'result',

420 id,

421 quality,

422 y: y32,

423 sampleRate: params.sample_rate,

424 strikeT: radEntry.strikeT.slice(),

425 strikeV: radEntry.strikeV.slice(),

426 strikeS: radEntry.strikeS.slice(),

427 modeF: solveEntry.modeF.slice(),

428 modeZeta: solveEntry.modeZeta.slice(),

429 modeZetaRad: solveEntry.modeZetaRad.slice(),

430 modeSigma: solveEntry.modeSigma.slice(),

431 modeN: solveEntry.modeN.slice(),

432 nodeR: solveEntry.nodeR.slice(),

433 nodeZ: solveEntry.nodeZ.slice(),

434 shapeW: solveEntry.shapeW.slice(),

435 cascadeDb: radEntry.cascadeDb.slice(),

436 cascadeBands: radEntry.cascadeBands,

437 cascadeFrames: radEntry.cascadeFrames,

438 cascadeFLo: solveEntry.modeF[0],

439 cascadeFHi: solveEntry.modeF[solveEntry.modeF.length - 1],

440 cascadeTEnd: radEntry.cascadeTEnd,

441 diag: {

442 modesKept: solveEntry.modeF.length,

443 maxN: solveEntry.maxN,

444 strikeCount: radEntry.strikeCount,

445 meanV: radEntry.meanV,

446 maxV: radEntry.maxV,

447 meanTauUs: radEntry.meanTauUs,

448 peakForceN: radEntry.peakForceN,

449 durationS: radEntry.durationS,

450 peakPa,

451 peakSplDb: 20 * Math.log10(peakPa / 20e-6 + 1e-30),

452 directPeakPa: radEntry.directPeakPa,

453 listenerDistM: radEntry.listenerDistM,

454 peakWrmsUm: radEntry.peakWrmsUm,

455 peakStrength: radEntry.peakStrength,

456 cascadeLostFrac: radEntry.cascadeLostFrac,

457 plateMassKg: solveEntry.plateMassKg,

458 warnings,

459 solveMs,

460 radiateMs,

461 solveCached,

462 radiateCached,

463 oversample: radEntry.oversample,

464 },

465 };

466 (self as unknown as Worker).postMessage(resp, [

467 resp.y.buffer,

468 resp.strikeT.buffer,

469 resp.strikeV.buffer,

470 resp.strikeS.buffer,

471 resp.modeF.buffer,

472 resp.modeZeta.buffer,

473 resp.modeZetaRad.buffer,

474 resp.modeSigma.buffer,

475 resp.modeN.buffer,

476 resp.nodeR.buffer,

477 resp.nodeZ.buffer,

478 resp.shapeW.buffer,

479 resp.cascadeDb.buffer,

480 ]);

481}