5778036cymbalsim: scripted plate and cymbal designer with a modal linear model and a von Karman GPU modelJeremy Magland 1// The modes of a plate up to a cutoff frequency, by spectrum slicing: the
2// eigenvalue count below a shift is the inertia of the banded factorization
3// (banded.ts), which fixes slices of about `sliceSize` modes between zero
4// and the cutoff, and each slice is solved by shift-invert Lanczos
5// (lanczos.ts). Inertia counts and slices are independent jobs, run here
6// through a JobRunner: in-process for the tests, or a pool of Web Workers
7// in the browser (pool.ts).
9import { type LinearSystem } from "../model/assemble.ts";
10import { inertia } from "./banded.ts";
11import { sliceModes } from "./lanczos.ts";
13export type Job =
14 | { kind: "inertia"; shift: number }
15 | { kind: "slice"; lo: number; hi: number; countLo: number; countHi: number; seed: number };
17export type JobResult =
18 | { kind: "inertia"; count: number }
19 | { kind: "slice"; lambda: number[]; shapes: Float64Array[]; expected: number; steps: number };
21export interface JobRunner {
22 run(job: Job): Promise<JobResult>;
23 close(): void;
24}
26/** Runs the jobs on the calling thread. */
27export function localRunner(sys: LinearSystem): JobRunner {
28 return {
29 run: async (job) => runJob(sys, job),
30 close: () => {},
31 };
32}
34export function runJob(sys: LinearSystem, job: Job): JobResult {
35 if (job.kind === "inertia") return { kind: "inertia", count: inertia(sys.band, sys.mass, sys.n, sys.bw, job.shift) };
36 const r = sliceModes(sys.band, sys.mass, sys.n, sys.bw, job.lo, job.hi, job.countLo, job.countHi, { seed: job.seed });
37 return { kind: "slice", lambda: r.lambda, shapes: r.shapes, expected: r.expected, steps: r.steps };
38}
40export interface ModeSet {
41 /** Number of unknowns of the system the shapes live on. */
42 n: number;
43 count: number;
44 /** Frequencies [Hz], ascending. */
45 freqs: Float64Array;
46 /** Mass-normalized shapes, mode-major: shapes[m * n + i]. */
47 shapes: Float32Array;
48 /** Modes the inertia counts announced but Lanczos did not deliver (zero when all is well). */
49 missing: number;
50 /** Modes below cutoff according to the inertia count (including any missing ones). */
51 expected: number;
52 /** Wall time of the computation [s]. */
53 seconds: number;
54}
56export interface ModeOptions {
76b7843Slice size follows the bandwidth (about bw/3.5); aluminum square presetJeremy Magland 57 /** Target number of modes per slice; by default defaultSliceSize(bw), which balances the factorizations against the Lanczos work. */
5778036cymbalsim: scripted plate and cymbal designer with a modal linear model and a von Karman GPU modelJeremy Magland 58 sliceSize?: number;
59 /** Progress callback with the number of finished and total slices. */
60 onProgress?: (done: number, total: number) => void;
61 /** Rigid-body modes of a free plate lie at zero; modes below this frequency are discarded [Hz]. Default 1. */
62 minFreq?: number;
63}
66 * The slice size that balances the cost of a slice's factorization, about
67 * n bw^2, against the Lanczos work per mode, whose reorthogonalization grows
68 * with the square of the slice: the optimum is near bw / 3.5, so that a
69 * finer grid, whose band is wider, uses fewer and larger slices.
70 */
71export function defaultSliceSize(bw: number): number {
72 return Math.max(20, Math.round(bw / 3.5));
73}
5778036cymbalsim: scripted plate and cymbal designer with a modal linear model and a von Karman GPU modelJeremy Magland 75interface Slice {
76 lo: number;
77 hi: number;
78 countLo: number;
79 countHi: number;
80}
82/**
83 * All modes of the system with frequencies up to fMax. The zero-frequency
84 * rigid modes of a free plate are computed and then discarded.
85 */
86export async function computeModes(sys: LinearSystem, fMax: number, runner: JobRunner, opts: ModeOptions = {}): Promise<ModeSet> {
87 const t0 = performance.now();
76b7843Slice size follows the bandwidth (about bw/3.5); aluminum square presetJeremy Magland 88 const sliceSize = opts.sliceSize ?? defaultSliceSize(sys.bw);
5778036cymbalsim: scripted plate and cymbal designer with a modal linear model and a von Karman GPU modelJeremy Magland 89 const minFreq = opts.minFreq ?? 1;
90 const w2 = (f: number) => (2 * Math.PI * f) ** 2;
91 const lamMin = -w2(1); // slightly below zero so that the rigid modes are counted in the first slice
92 const lamMax = w2(fMax);
93 const count = async (shift: number) => {
94 const r = await runner.run({ kind: "inertia", shift });
95 return (r as { count: number }).count;
96 };
97 const [nMin, nMax] = await Promise.all([count(lamMin), count(lamMax)]);
98 const expected = nMax - nMin;
99 const empty = (): ModeSet => ({ n: sys.n, count: 0, freqs: new Float64Array(0), shapes: new Float32Array(0), missing: 0, expected, seconds: (performance.now() - t0) / 1000 });
100 if (expected <= 0) return empty();
102 // Boundaries equally spaced in omega (the mode count of a plate grows about
103 // linearly with frequency), refined where a slice is too full.
104 let slices: Slice[] = [];
105 {
106 const S = Math.max(1, Math.ceil(expected / sliceSize));
107 const wMax = Math.sqrt(lamMax);
108 const bounds: number[] = [lamMin];
109 for (let i = 1; i < S; i++) bounds.push(((wMax * i) / S) ** 2);
110 bounds.push(lamMax);
111 const counts = [nMin, ...(await Promise.all(bounds.slice(1, -1).map(count))), nMax];
112 for (let i = 0; i < bounds.length - 1; i++) slices.push({ lo: bounds[i], hi: bounds[i + 1], countLo: counts[i], countHi: counts[i + 1] });
113 }
114 for (let round = 0; round < 6; round++) {
115 const full = slices.filter((s) => s.countHi - s.countLo > 1.5 * sliceSize && s.hi - s.lo > w2(0.5));
116 if (full.length === 0) break;
117 const mids = await Promise.all(full.map((s) => count(((Math.sqrt(Math.max(s.lo, 0)) + Math.sqrt(s.hi)) / 2) ** 2)));
118 const next: Slice[] = [];
119 for (const s of slices) {
120 const idx = full.indexOf(s);
121 if (idx < 0) { next.push(s); continue; }
122 const mid = ((Math.sqrt(Math.max(s.lo, 0)) + Math.sqrt(s.hi)) / 2) ** 2;
123 next.push({ lo: s.lo, hi: mid, countLo: s.countLo, countHi: mids[idx] }, { lo: mid, hi: s.hi, countLo: mids[idx], countHi: s.countHi });
124 }
125 slices = next;
126 }
127 // Merge neighbors while the union stays within the target, and drop empty slices.
128 {
129 const merged: Slice[] = [];
130 for (const s of slices) {
131 const last = merged[merged.length - 1];
132 if (last && last.countHi === s.countLo && s.countHi - last.countLo <= sliceSize) { last.hi = s.hi; last.countHi = s.countHi; }
133 else merged.push({ ...s });
134 }
135 slices = merged.filter((s) => s.countHi > s.countLo);
136 }
138 let done = 0;
139 opts.onProgress?.(0, slices.length);
140 const results = await Promise.all(
141 slices.map(async (s, i) => {
142 const r = (await runner.run({ kind: "slice", lo: s.lo, hi: s.hi, countLo: s.countLo, countHi: s.countHi, seed: 1000 + 7919 * i })) as Extract<JobResult, { kind: "slice" }>;
143 done++;
144 opts.onProgress?.(done, slices.length);
145 return r;
146 }),
147 );
148 const pairs: { lambda: number; shape: Float64Array }[] = [];
149 let missing = 0;
150 for (const r of results) {
151 missing += r.expected - r.lambda.length;
152 for (let i = 0; i < r.lambda.length; i++) pairs.push({ lambda: r.lambda[i], shape: r.shapes[i] });
153 }
154 pairs.sort((a, b) => a.lambda - b.lambda);
155 const kept = pairs.filter((p) => p.lambda > w2(minFreq));
156 const n = sys.n;
157 const shapes = new Float32Array(kept.length * n);
158 const freqs = new Float64Array(kept.length);
159 kept.forEach((p, m) => {
160 freqs[m] = Math.sqrt(p.lambda) / (2 * Math.PI);
161 shapes.set(p.shape, m * n);
162 });
163 return { n, count: kept.length, freqs, shapes, missing, expected, seconds: (performance.now() - t0) / 1000 };
164}