// Web Worker: runs the numbl script, computes linear modes, and simulates. import { growthRate, relDiff } from "./sim/analysis"; import { computeModes, isPlanar, modalDamping } from "./sim/modes"; import { numblEig, runScript } from "./sim/runScript"; import { simulateGPU } from "./sim/gpu"; import { simulate, type SimOptions } from "./sim/solver"; import { runTwin } from "./sim/twin"; import type { BeadSystem, SimResult } from "./sim/types"; const helperSources = import.meta.glob("../matlab/helpers/*.m", { query: "?raw", import: "default", eager: true }) as Record; const helpers = Object.entries(helperSources).map(([path, source]) => ({ name: path.split("/").pop()!, source })); const MAX_MODE_BEADS = 1500; // flat nets: z motion only const MAX_FULL_DOF = 1200; // curved nets: all three components of every free bead const MAX_TWIN_BYTES = 400e6; export interface RunPayload { sys: { n: number; pos0: Float64Array; pinned: Uint8Array; si: Int32Array; sj: Int32Array; nDashpots: number; nBends: number; nHinges: number; listenQuantity: string; engine: string; listenBeads: Int32Array; hitBeads: Int32Array; method: string; duration: number; }; result: SimResult; modeFreqs: Float64Array | null; modeKind: "transverse" | "full" | null; modeDecay: Float64Array | null; // small-amplitude amplitude decay rate of each mode (1/s) notes: string[]; check: { dt: number; windows: { t: number; rel: number }[] } | null; // the same run with every hit scaled by 1e-3, its output scaled back up linear: { audio: Float32Array; windows: { t: number; rel: number }[] } | null; // the same run with every hit scaled by (1 + eps): separation of the bead // positions (max over coordinates, relative to the largest displacement) twin: { eps: number; audio: Float32Array; times: Float64Array; sep: Float64Array; rate: number | null } | null; } export interface RunRequest { type: "run"; script: string; dtCheck: boolean; linear: boolean; twin: boolean; modes: boolean; // the page's "linear modes" toggle (sys.modes in the script overrides it) gpu: boolean; } export type FromWorker = | { type: "log"; text: string } | { type: "stage"; text: string } | { type: "progress"; t: number; duration: number; label: string } | { type: "result"; payload: RunPayload } | { type: "error"; message: string }; const post = (m: FromWorker, transfer: Transferable[] = []) => (self as unknown as Worker).postMessage(m, transfer); const windowsOf = (duration: number) => [0.01, 0.03, 0.1, 0.3, 1, duration].filter((t, i, a) => t <= duration && a.indexOf(t) === i); self.onmessage = async (ev: MessageEvent) => { const { script, dtCheck, linear, twin: twinWanted, modes: modesToggle, gpu } = ev.data; // the GPU path for the main run and the checks; the CPU otherwise const sim = (s: BeadSystem, o: SimOptions) => (gpu ? simulateGPU(s, o) : Promise.resolve(simulate(s, o))); try { post({ type: "stage", text: "Running script" }); const { sys, warnings } = runScript(script, helpers, text => post({ type: "log", text })); const notes = [...warnings]; let modes = null; const nFree = sys.pinned.reduce((a, p) => a + (p ? 0 : 1), 0); if (sys.modes ?? modesToggle) { const planar = isPlanar(sys); if (planar && nFree > MAX_MODE_BEADS) notes.push(`Linear modes skipped: ${nFree} free beads (limit ${MAX_MODE_BEADS} for a flat net).`); else if (!planar && 3 * nFree > MAX_FULL_DOF) { notes.push( `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).` ); } else { post({ type: "stage", text: planar ? `Computing linear modes (${nFree} free beads)` : `Computing 3D linear modes (${3 * nFree} coordinates)` }); modes = computeModes(sys, numblEig, planar ? "transverse" : "full"); const fTop = modes.freqs[modes.freqs.length - 1]; const nZero = modes.freqs.filter(f => f < 1e-6 * fTop).length; if (nZero > 0) { notes.push( `${nZero} mode${nZero > 1 ? "s are" : " is"} at (nearly) 0 Hz: the beads can move or tilt freely as a whole. ` + "Pin more beads (a cluster rather than one bead), or tie the object to pinned anchor beads with soft springs." ); } } } // positions at every frame, kept for the perturbed twin const saved: Float64Array[] = []; const savedBytes = Math.round(sys.duration * sys.frameRate + 1) * 3 * sys.n * 8; if (twinWanted && gpu) { notes.push("Perturbed twin skipped: a perturbation of one part in 1e9 is below single precision; turn off WebGPU to run it."); } const twin = twinWanted && !gpu; const keep = twin && savedBytes <= MAX_TWIN_BYTES; 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.`); post({ type: "stage", text: gpu ? "Simulating (WebGPU)" : "Simulating" }); const result = await sim(sys, { modes, onProgress: t => post({ type: "progress", t, duration: sys.duration, label: "Simulating" }), onFrame: keep ? (_f, _t, x) => saved.push(x.slice()) : undefined, }); let checkOut: RunPayload["check"] = null; if (dtCheck && !result.aborted) { post({ type: "stage", text: "Simulating again with dt/2" }); const half: BeadSystem = { ...sys, dt: result.dt / 2 }; const r2 = await sim(half, { maxAnimBytes: 0, onProgress: t => post({ type: "progress", t, duration: sys.duration, label: "Simulating with dt/2" }), }); const fs = result.sampleRate; const windows = windowsOf(sys.duration).map(t => ({ t, rel: relDiff(r2.audio, result.audio, Math.round(t * fs)) })); checkOut = { dt: r2.dt, windows }; } let linearOut: RunPayload["linear"] = null; if (linear && !result.aborted) { post({ type: "stage", text: "Simulating the small-strike reference" }); const s = 1e-3; const small: BeadSystem = { ...sys, dt: result.dt, hits: sys.hits.map(h => ({ ...h, velocity: [h.velocity[0] * s, h.velocity[1] * s, h.velocity[2] * s], force: [h.force[0] * s, h.force[1] * s, h.force[2] * s], })), }; const r3 = await sim(small, { maxAnimBytes: 0, onProgress: t => post({ type: "progress", t, duration: sys.duration, label: "Small-strike reference" }), }); const audio = r3.audio.map(v => v / s); const fs = result.sampleRate; const windows = windowsOf(sys.duration).map(t => ({ t, rel: relDiff(audio, result.audio, Math.round(t * fs)) })); linearOut = { audio, windows }; } let twinOut: RunPayload["twin"] = null; if (keep && !result.aborted) { post({ type: "stage", text: "Simulating the perturbed twin" }); const eps = 1e-9; const { audio, sep } = runTwin(sys, saved, result.dt, eps, t => post({ type: "progress", t, duration: sys.duration, label: "Perturbed twin" }) ); const times = result.frameTimes.slice(0, saved.length); twinOut = { eps, audio, times, sep, rate: growthRate(times, sep, eps) }; } const hitBeads = Int32Array.from(new Set(sys.hits.flatMap(h => Array.from(h.beads)))); const payload: RunPayload = { sys: { n: sys.n, pos0: sys.pos0, pinned: sys.pinned, si: sys.si, sj: sys.sj, nDashpots: sys.dc.filter((c, d) => c !== 0 || sys.dca[d] !== 0).length, nBends: sys.bi.length, nHinges: sys.hi.length, listenQuantity: sys.listenQuantity, engine: gpu ? "WebGPU, single precision" : "CPU, double precision", listenBeads: sys.listenBeads, hitBeads, method: sys.method, duration: sys.duration, }, result, modeFreqs: modes ? modes.freqs : null, modeKind: modes ? modes.kind : null, modeDecay: modes ? modalDamping(sys, modes) : null, notes, check: checkOut, linear: linearOut, twin: twinOut, }; const transfer: Transferable[] = [result.audio.buffer, ...result.animFrames.map(f => f.buffer)]; if (result.modeEnergy) transfer.push(result.modeEnergy.buffer); if (linearOut) transfer.push(linearOut.audio.buffer); if (twinOut) transfer.push(twinOut.audio.buffer, twinOut.sep.buffer); post({ type: "result", payload }, transfer); } catch (e) { const err = e as { message?: string; line?: number; file?: string }; const where = err.line ? ` (line ${err.line})` : ""; post({ type: "error", message: `${err.message ?? String(e)}${where}` }); } };