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 {
57 /** Target number of modes per slice; 80 by default. */
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}
65interface Slice {
66 lo: number;
67 hi: number;
68 countLo: number;
69 countHi: number;
70}
72/**
73 * All modes of the system with frequencies up to fMax. The zero-frequency
74 * rigid modes of a free plate are computed and then discarded.
75 */
76export async function computeModes(sys: LinearSystem, fMax: number, runner: JobRunner, opts: ModeOptions = {}): Promise<ModeSet> {
77 const t0 = performance.now();
78 const sliceSize = opts.sliceSize ?? 80;
79 const minFreq = opts.minFreq ?? 1;
80 const w2 = (f: number) => (2 * Math.PI * f) ** 2;
81 const lamMin = -w2(1); // slightly below zero so that the rigid modes are counted in the first slice
82 const lamMax = w2(fMax);
83 const count = async (shift: number) => {
84 const r = await runner.run({ kind: "inertia", shift });
85 return (r as { count: number }).count;
86 };
87 const [nMin, nMax] = await Promise.all([count(lamMin), count(lamMax)]);
88 const expected = nMax - nMin;
89 const empty = (): ModeSet => ({ n: sys.n, count: 0, freqs: new Float64Array(0), shapes: new Float32Array(0), missing: 0, expected, seconds: (performance.now() - t0) / 1000 });
90 if (expected <= 0) return empty();
92 // Boundaries equally spaced in omega (the mode count of a plate grows about
93 // linearly with frequency), refined where a slice is too full.
94 let slices: Slice[] = [];
95 {
96 const S = Math.max(1, Math.ceil(expected / sliceSize));
97 const wMax = Math.sqrt(lamMax);
98 const bounds: number[] = [lamMin];
99 for (let i = 1; i < S; i++) bounds.push(((wMax * i) / S) ** 2);
100 bounds.push(lamMax);
101 const counts = [nMin, ...(await Promise.all(bounds.slice(1, -1).map(count))), nMax];
102 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] });
103 }
104 for (let round = 0; round < 6; round++) {
105 const full = slices.filter((s) => s.countHi - s.countLo > 1.5 * sliceSize && s.hi - s.lo > w2(0.5));
106 if (full.length === 0) break;
107 const mids = await Promise.all(full.map((s) => count(((Math.sqrt(Math.max(s.lo, 0)) + Math.sqrt(s.hi)) / 2) ** 2)));
108 const next: Slice[] = [];
109 for (const s of slices) {
110 const idx = full.indexOf(s);
111 if (idx < 0) { next.push(s); continue; }
112 const mid = ((Math.sqrt(Math.max(s.lo, 0)) + Math.sqrt(s.hi)) / 2) ** 2;
113 next.push({ lo: s.lo, hi: mid, countLo: s.countLo, countHi: mids[idx] }, { lo: mid, hi: s.hi, countLo: mids[idx], countHi: s.countHi });
114 }
115 slices = next;
116 }
117 // Merge neighbors while the union stays within the target, and drop empty slices.
118 {
119 const merged: Slice[] = [];
120 for (const s of slices) {
121 const last = merged[merged.length - 1];
122 if (last && last.countHi === s.countLo && s.countHi - last.countLo <= sliceSize) { last.hi = s.hi; last.countHi = s.countHi; }
123 else merged.push({ ...s });
124 }
125 slices = merged.filter((s) => s.countHi > s.countLo);
126 }
128 let done = 0;
129 opts.onProgress?.(0, slices.length);
130 const results = await Promise.all(
131 slices.map(async (s, i) => {
132 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" }>;
133 done++;
134 opts.onProgress?.(done, slices.length);
135 return r;
136 }),
137 );
138 const pairs: { lambda: number; shape: Float64Array }[] = [];
139 let missing = 0;
140 for (const r of results) {
141 missing += r.expected - r.lambda.length;
142 for (let i = 0; i < r.lambda.length; i++) pairs.push({ lambda: r.lambda[i], shape: r.shapes[i] });
143 }
144 pairs.sort((a, b) => a.lambda - b.lambda);
145 const kept = pairs.filter((p) => p.lambda > w2(minFreq));
146 const n = sys.n;
147 const shapes = new Float32Array(kept.length * n);
148 const freqs = new Float64Array(kept.length);
149 kept.forEach((p, m) => {
150 freqs[m] = Math.sqrt(p.lambda) / (2 * Math.PI);
151 shapes.set(p.shape, m * n);
152 });
153 return { n, count: kept.length, freqs, shapes, missing, expected, seconds: (performance.now() - t0) / 1000 };
154}