1793cecTime-domain nonlinear shell cymbal after Samejima (2021)Jeremy Magland 1// Web Worker that renders the cymbal2 stroke by executing the MATLAB
2// implementation in ../matlab through numbl. Unlike the sibling cymbal app
3// there is no cheap stage to cache around: the cost IS the time integration,
4// one long executeCode call, so requests carry a quality tier
5// preview -> oversample 1, tail clamped to 1.5 s (the azimuthal grid is
6// kept: hard strokes need the full wavenumber band to carry
7// the cascade, and coarsening it destabilizes the integration)
8// full -> the parameters as given
9// and the UI shows elapsed time rather than a meaningful progress fraction.
10// A cheap linear-modes stage (cymbal2_modes per wavenumber) is cached on the
11// geometry so the spectrum plot can overlay the mode table immediately.
13import { executeCode, type RuntimeValue, type RuntimeTensor } from 'numbl';
14import { hashMFiles, type MFile } from '../mfiles';
15import { MODES_KEYS, type Cymbal2Params } from '../params';
17export interface RenderRequest {
18 type: 'render';
19 id: number;
20 params: Cymbal2Params;
21 quality: 'preview' | 'full';
22 mfiles: MFile[];
23}
25export interface RenderResponse {
26 type: 'result';
27 id: number;
28 quality: 'preview' | 'full';
29 y: Float32Array; // final signal, [-1, 1]
30 sampleRate: number;
31 forceT: Float32Array; // contact force history (max-pooled), s
32 forceN: Float32Array; // and its values, N
33 modeF: Float32Array; // linear mode frequencies, Hz
34 modeM: Float32Array; // nodal diameters of each
35 diag: {
36 strokesPlayed: number;
37 meanV: number;
38 contactVel: number;
39 peakForceN: number;
40 contactMs: number;
41 contactEvents: number;
42 peakXiMm: number;
43 peakDispMm: number;
44 peakDispThickness: number;
45 peakAccel: number;
46 strikeRadiusM: number;
47 pickupRadiusM: number;
48 RsM: number;
49 nsteps: number;
50 simRate: number;
51 nDropped: number;
52 durationS: number;
53 renderMs: number;
54 oversample: number;
55 nTheta: number;
56 warnings: string[];
57 };
58}
60export interface ProgressMessage {
61 type: 'progress';
62 id: number;
63 stage: 'modes' | 'integrate' | 'output';
64 frac: number;
65}
67export interface ErrorMessage {
68 type: 'error';
69 id: number;
70 message: string;
71 file?: string;
72 line?: number;
73}
75interface ModesEntry {
76 modeF: Float32Array;
77 modeM: Float32Array;
78}
80const modesCache = new Map<string, ModesEntry>();
81const MODES_CACHE_MAX = 4;
83function lruGet<T>(cache: Map<string, T>, key: string): T | undefined {
84 const v = cache.get(key);
85 if (v !== undefined) {
86 cache.delete(key);
87 cache.set(key, v);
88 }
89 return v;
90}
92function lruSet<T>(cache: Map<string, T>, key: string, value: T, max: number): void {
93 cache.delete(key);
94 cache.set(key, value);
95 while (cache.size > max) {
96 const oldest = cache.keys().next().value as string;
97 cache.delete(oldest);
98 }
99}
101function keyFor(params: Cymbal2Params, keys: readonly (keyof Cymbal2Params)[]): string {
102 return keys.map((k) => String(params[k])).join('|');
103}
105// MATLAB source lines assigning every parameter onto a struct P.
106function pLiteral(params: Cymbal2Params): string {
107 const lines: string[] = [];
108 for (const [k, v] of Object.entries(params)) {
109 if (k === 'bit_depth') continue;
110 if (typeof v === 'number') lines.push(`P.${k} = ${v};`);
111 else if (typeof v === 'boolean') lines.push(`P.${k} = ${v ? 1 : 0};`);
112 }
113 return lines.join('\n');
114}
116function runStage(mfiles: MFile[], source: string, inputs?: Record<string, RuntimeValue>) {
117 return executeCode(
118 source,
119 { initialVariableValues: inputs, displayResults: false },
120 mfiles,
121 'stage.m',
122 );
123}
125function asScalar(v: RuntimeValue): number {
126 if (typeof v === 'number') return v;
127 const t = v as RuntimeTensor;
128 if (t && t.kind === 'tensor' && t.data.length === 1) return t.data[0];
129 if (typeof v === 'boolean') return v ? 1 : 0;
130 throw new Error('expected scalar result from numbl stage');
131}
133function toF32(v: RuntimeValue): Float32Array {
134 const t = v as RuntimeTensor;
135 if (!t || t.kind !== 'tensor') throw new Error('expected tensor result from numbl stage');
136 const out = new Float32Array(t.data.length);
137 for (let i = 0; i < t.data.length; i++) out[i] = t.data[i];
138 return out;
139}
141const SKIP_WARNING = /skipping (\S+\.m) \(syntax error at line (\d+)\)/;
143self.onmessage = (ev: MessageEvent<RenderRequest>) => {
144 const msg = ev.data;
145 if (msg.type !== 'render') return;
147 const warnings: string[] = [];
148 const origWarn = console.warn;
149 console.warn = (...args: unknown[]) => {
150 warnings.push(args.map(String).join(' '));
151 };
153 try {
154 render(msg);
155 } catch (err) {
156 const e = err as { message?: string; file?: string; line?: number };
157 let m: ErrorMessage = {
158 type: 'error',
159 id: msg.id,
160 message: e?.message ? String(e.message) : String(err),
161 file: typeof e?.file === 'string' ? e.file : undefined,
162 line: typeof e?.line === 'number' ? e.line : undefined,
163 };
164 for (const w of warnings) {
165 const hit = SKIP_WARNING.exec(w);
166 if (hit) {
167 m = { type: 'error', id: msg.id, message: 'syntax error', file: hit[1], line: Number(hit[2]) };
168 break;
169 }
170 }
171 (self as unknown as Worker).postMessage(m);
172 } finally {
173 console.warn = origWarn;
174 }
175};
177function render(msg: RenderRequest): void {
178 const { id, mfiles, quality } = msg;
179 const params = { ...msg.params };
180 if (quality === 'preview') {
181 params.oversample = 1;
182 params.tail_s = Math.min(params.tail_s, 1.5);
183 }
185 const post = (m: ProgressMessage) => (self as unknown as Worker).postMessage(m);
186 const codeKey = hashMFiles(mfiles);
188 // ---- linear modes, cheap, cached on geometry ----
189 const modesKey = codeKey + '|' + keyFor(params, MODES_KEYS);
190 let modesEntry = lruGet(modesCache, modesKey);
191 if (!modesEntry) {
192 post({ type: 'progress', id, stage: 'modes', frac: 0.05 });
193 const r = runStage(mfiles, `
194${pLiteral(params)}
195mF = [];
196mM = [];
197for mw = 0:min(10, floor(P.n_theta / 2))
198 f = cymbal2_modes(P, mw);
199 f = f(f > 0.5 & f < 4000);
200 nk = min(4, length(f));
201 mF = [mF; f(1:nk)];
202 mM = [mM; mw * ones(nk, 1)];
203end
204`);
205 modesEntry = {
206 modeF: toF32(r.variableValues.mF),
207 modeM: toF32(r.variableValues.mM),
208 };
209 lruSet(modesCache, modesKey, modesEntry, MODES_CACHE_MAX);
210 }
212 // ---- the time integration, one long call ----
213 post({ type: 'progress', id, stage: 'integrate', frac: 0.2 });
214 const t0 = performance.now();
215 const r = runStage(mfiles, `
216${pLiteral(params)}
217S = cymbal2_system(P);
218MB = cymbal2_mallet(P);
219[acc, dinfo, fhist] = cymbal2_sim(P, S, MB);
220[y, peakAccel] = cymbal2_output(P, acc);
221% max-pool the contact force history to about 2000 points for the plot
222nsim = length(fhist);
223np = min(2000, nsim);
224blk = ceil(nsim / np);
225np = ceil(nsim / blk);
226fpool = max(reshape([fhist, zeros(1, np * blk - nsim)], blk, np), [], 1);
227tpool = ((0:np - 1) * blk + blk / 2) / dinfo.simRate;
228d_strokes = dinfo.strokesPlayed;
229d_meanV = dinfo.meanV;
230d_cvel = dinfo.contactVel(1);
231d_peakF = dinfo.peakForceN;
232d_contactMs = dinfo.contactMs;
233d_events = dinfo.contactEvents;
234d_xi = dinfo.peakXiMm;
235d_disp = dinfo.peakDispM;
236d_dispH = dinfo.peakDispThickness;
237d_nsteps = dinfo.nsteps;
238d_simRate = dinfo.simRate;
239d_dropped = S.nDropped;
240d_strikeR = S.strikeRadiusM;
241d_pickupR = S.pickupRadiusM;
242d_Rs = S.Rs;
243`);
244 const renderMs = performance.now() - t0;
246 post({ type: 'progress', id, stage: 'output', frac: 0.95 });
247 const y32 = toF32(r.variableValues.y);
248 const peakAccel = asScalar(r.variableValues.peakAccel);
249 const nDropped = asScalar(r.variableValues.d_dropped);
251 const warnings: string[] = [];
252 if (nDropped > 0) {
253 warnings.push(
254 `${nDropped} azimuthal columns were truncated: the radial grid cannot resolve their inner boundary layer; raise n_radial`,
255 );
256 }
257 if (!params.normalize) {
258 let peak = 0;
259 for (let i = 0; i < y32.length; i++) peak = Math.max(peak, Math.abs(y32[i]));
260 if (peak > 1.0) {
261 warnings.push(
262 `clipping: peak acceleration ${peakAccel.toExponential(2)} m/s² exceeds full_scale_accel`,
263 );
264 }
265 }
267 const resp: RenderResponse = {
268 type: 'result',
269 id,
270 quality,
271 y: y32,
272 sampleRate: params.sample_rate,
273 forceT: toF32(r.variableValues.tpool),
274 forceN: toF32(r.variableValues.fpool),
275 modeF: modesEntry.modeF.slice(),
276 modeM: modesEntry.modeM.slice(),
277 diag: {
278 strokesPlayed: asScalar(r.variableValues.d_strokes),
279 meanV: asScalar(r.variableValues.d_meanV),
280 contactVel: asScalar(r.variableValues.d_cvel),
281 peakForceN: asScalar(r.variableValues.d_peakF),
282 contactMs: asScalar(r.variableValues.d_contactMs),
283 contactEvents: asScalar(r.variableValues.d_events),
284 peakXiMm: asScalar(r.variableValues.d_xi),
285 peakDispMm: asScalar(r.variableValues.d_disp) * 1000,
286 peakDispThickness: asScalar(r.variableValues.d_dispH),
287 peakAccel,
288 strikeRadiusM: asScalar(r.variableValues.d_strikeR),
289 pickupRadiusM: asScalar(r.variableValues.d_pickupR),
290 RsM: asScalar(r.variableValues.d_Rs),
291 nsteps: asScalar(r.variableValues.d_nsteps),
292 simRate: asScalar(r.variableValues.d_simRate),
293 nDropped,
294 durationS: y32.length / params.sample_rate,
295 renderMs,
296 oversample: params.oversample,
297 nTheta: params.n_theta,
298 warnings,
299 },
300 };
301 (self as unknown as Worker).postMessage(resp, [
302 resp.y.buffer,
303 resp.forceT.buffer,
304 resp.forceN.buffer,
305 resp.modeF.buffer,
306 resp.modeM.buffer,
307 ]);
308}