1// Web Worker: runs the numbl script, computes linear modes, and simulates.
3import { growthRate, relDiff } from "./sim/analysis";
4import { computeModes, isPlanar, modalDamping } from "./sim/modes";
5import { numblEig, runScript } from "./sim/runScript";
6import { simulate } from "./sim/solver";
7import { runTwin } from "./sim/twin";
8import type { BeadSystem, SimResult } from "./sim/types";
10const helperSources = import.meta.glob("../matlab/helpers/*.m", { query: "?raw", import: "default", eager: true }) as Record<string, string>;
11const helpers = Object.entries(helperSources).map(([path, source]) => ({ name: path.split("/").pop()!, source }));
13const MAX_MODE_BEADS = 1500; // flat nets: z motion only
14const MAX_FULL_DOF = 1200; // curved nets: all three components of every free bead
15const MAX_TWIN_BYTES = 400e6;
17export interface RunPayload {
18 sys: {
19 n: number;
20 pos0: Float64Array;
21 pinned: Uint8Array;
22 si: Int32Array;
23 sj: Int32Array;
24 nDashpots: number;
25 nBends: number;
26 nHinges: number;
27 listenQuantity: string;
28 listenBeads: Int32Array;
29 hitBeads: Int32Array;
30 method: string;
31 duration: number;
32 };
33 result: SimResult;
34 modeFreqs: Float64Array | null;
35 modeKind: "transverse" | "full" | null;
36 modeDecay: Float64Array | null; // small-amplitude amplitude decay rate of each mode (1/s)
37 notes: string[];
38 check: { dt: number; windows: { t: number; rel: number }[] } | null;
39 // the same run with every hit scaled by 1e-3, its output scaled back up
40 linear: { audio: Float32Array; windows: { t: number; rel: number }[] } | null;
41 // the same run with every hit scaled by (1 + eps): separation of the bead
42 // positions (max over coordinates, relative to the largest displacement)
43 twin: { eps: number; audio: Float32Array; times: Float64Array; sep: Float64Array; rate: number | null } | null;
44}
46export interface RunRequest {
47 type: "run";
48 script: string;
49 dtCheck: boolean;
50 linear: boolean;
51 twin: boolean;
52}
54export type FromWorker =
55 | { type: "log"; text: string }
56 | { type: "stage"; text: string }
57 | { type: "progress"; t: number; duration: number; label: string }
58 | { type: "result"; payload: RunPayload }
59 | { type: "error"; message: string };
61const post = (m: FromWorker, transfer: Transferable[] = []) => (self as unknown as Worker).postMessage(m, transfer);
63const windowsOf = (duration: number) =>
64 [0.01, 0.03, 0.1, 0.3, 1, duration].filter((t, i, a) => t <= duration && a.indexOf(t) === i);
66self.onmessage = (ev: MessageEvent<RunRequest>) => {
67 const { script, dtCheck, linear, twin } = ev.data;
68 try {
69 post({ type: "stage", text: "Running script" });
70 const { sys, warnings } = runScript(script, helpers, text => post({ type: "log", text }));
71 const notes = [...warnings];
73 let modes = null;
74 const nFree = sys.pinned.reduce((a, p) => a + (p ? 0 : 1), 0);
75 if (sys.modes) {
76 const planar = isPlanar(sys);
77 if (planar && nFree > MAX_MODE_BEADS) notes.push(`Linear modes skipped: ${nFree} free beads (limit ${MAX_MODE_BEADS} for a flat net).`);
78 else if (!planar && 3 * nFree > MAX_FULL_DOF) {
79 notes.push(
80 `Linear modes skipped: the net is curved, so its modes need all ${3 * nFree} coordinates of the free beads (limit ${MAX_FULL_DOF}, about ${MAX_FULL_DOF / 3} free beads).`
81 );
82 } else {
83 post({ type: "stage", text: planar ? `Computing linear modes (${nFree} free beads)` : `Computing 3D linear modes (${3 * nFree} coordinates)` });
84 modes = computeModes(sys, numblEig, planar ? "transverse" : "full");
85 const fTop = modes.freqs[modes.freqs.length - 1];
86 const nZero = modes.freqs.filter(f => f < 1e-6 * fTop).length;
87 if (nZero > 0) {
88 notes.push(
89 `${nZero} mode${nZero > 1 ? "s are" : " is"} at (nearly) 0 Hz: the beads can move or tilt freely as a whole. ` +
90 "Pin more beads (a cluster rather than one bead), or tie the object to pinned anchor beads with soft springs."
91 );
92 }
93 }
94 }
96 // positions at every frame, kept for the perturbed twin
97 const saved: Float64Array[] = [];
98 const savedBytes = Math.round(sys.duration * sys.frameRate + 1) * 3 * sys.n * 8;
99 const keep = twin && savedBytes <= MAX_TWIN_BYTES;
100 if (twin && !keep) notes.push(`Perturbed twin skipped: it would store ${(savedBytes / 1e6).toFixed(0)} MB of positions (limit ${MAX_TWIN_BYTES / 1e6} MB). Lower sys.frameRate or sys.duration.`);
102 post({ type: "stage", text: "Simulating" });
103 const result = simulate(sys, {
104 modes,
105 onProgress: t => post({ type: "progress", t, duration: sys.duration, label: "Simulating" }),
106 onFrame: keep ? (_f, _t, x) => saved.push(x.slice()) : undefined,
107 });
109 let checkOut: RunPayload["check"] = null;
110 if (dtCheck && !result.aborted) {
111 post({ type: "stage", text: "Simulating again with dt/2" });
112 const half: BeadSystem = { ...sys, dt: result.dt / 2 };
113 const r2 = simulate(half, {
114 maxAnimBytes: 0,
115 onProgress: t => post({ type: "progress", t, duration: sys.duration, label: "Simulating with dt/2" }),
116 });
117 const fs = result.sampleRate;
118 const windows = windowsOf(sys.duration).map(t => ({ t, rel: relDiff(r2.audio, result.audio, Math.round(t * fs)) }));
119 checkOut = { dt: r2.dt, windows };
120 }
122 let linearOut: RunPayload["linear"] = null;
123 if (linear && !result.aborted) {
124 post({ type: "stage", text: "Simulating the small-strike reference" });
125 const s = 1e-3;
126 const small: BeadSystem = {
127 ...sys,
128 dt: result.dt,
129 hits: sys.hits.map(h => ({
130 ...h,
131 velocity: [h.velocity[0] * s, h.velocity[1] * s, h.velocity[2] * s],
132 force: [h.force[0] * s, h.force[1] * s, h.force[2] * s],
133 })),
134 };
135 const r3 = simulate(small, {
136 maxAnimBytes: 0,
137 onProgress: t => post({ type: "progress", t, duration: sys.duration, label: "Small-strike reference" }),
138 });
139 const audio = r3.audio.map(v => v / s);
140 const fs = result.sampleRate;
141 const windows = windowsOf(sys.duration).map(t => ({ t, rel: relDiff(audio, result.audio, Math.round(t * fs)) }));
142 linearOut = { audio, windows };
143 }
145 let twinOut: RunPayload["twin"] = null;
146 if (keep && !result.aborted) {
147 post({ type: "stage", text: "Simulating the perturbed twin" });
148 const eps = 1e-9;
149 const { audio, sep } = runTwin(sys, saved, result.dt, eps, t =>
150 post({ type: "progress", t, duration: sys.duration, label: "Perturbed twin" })
151 );
152 const times = result.frameTimes.slice(0, saved.length);
153 twinOut = { eps, audio, times, sep, rate: growthRate(times, sep, eps) };
154 }
156 const hitBeads = Int32Array.from(new Set(sys.hits.flatMap(h => Array.from(h.beads))));
157 const payload: RunPayload = {
158 sys: {
159 n: sys.n,
160 pos0: sys.pos0,
161 pinned: sys.pinned,
162 si: sys.si,
163 sj: sys.sj,
164 nDashpots: sys.dc.filter((c, d) => c !== 0 || sys.dca[d] !== 0).length,
165 nBends: sys.bi.length,
166 nHinges: sys.hi.length,
167 listenQuantity: sys.listenQuantity,
168 listenBeads: sys.listenBeads,
169 hitBeads,
170 method: sys.method,
171 duration: sys.duration,
172 },
173 result,
174 modeFreqs: modes ? modes.freqs : null,
175 modeKind: modes ? modes.kind : null,
176 modeDecay: modes ? modalDamping(sys, modes) : null,
177 notes,
178 check: checkOut,
179 linear: linearOut,
180 twin: twinOut,
181 };
182 const transfer: Transferable[] = [result.audio.buffer, ...result.animFrames.map(f => f.buffer)];
183 if (result.modeEnergy) transfer.push(result.modeEnergy.buffer);
184 if (linearOut) transfer.push(linearOut.audio.buffer);
185 if (twinOut) transfer.push(twinOut.audio.buffer, twinOut.sep.buffer);
186 post({ type: "result", payload }, transfer);
187 } catch (e) {
188 const err = e as { message?: string; line?: number; file?: string };
189 const where = err.line ? ` (line ${err.line})` : "";
190 post({ type: "error", message: `${err.message ?? String(e)}${where}` });
191 }
192};