1// Re-render the listener signals at any position from the acceleration record, with the same
2// free-field model as cymsim's built-in listener (cymsimweb src/model/radiation.ts): a Rayleigh
3// integral over element-mean normal accelerations for the front path and a diffracted path over
4// the rim for the back term, each element delayed by its distance at a 16-tap windowed-sinc
5// fractional delay. The record already holds the audio-rate node accelerations, so this is a
6// gather over elements per output sample; cymsim's test suite validates the same reconstruction
7// against the built-in listener to 1.5e-4 relative rms. The gather runs on WebGPU when available
8// (one workgroup per output sample, the record streamed through a sliding window) and on the CPU
9// otherwise. besselI0 and fractionalDelayTaps are copied from cymsimweb.
10import type { CymsimRun } from "./io/cymsimFile.ts";
12export const DELAY_TAPS = 16;
14/** Modified Bessel function I0 by its power series (converges quickly for x < 30). */
15function besselI0(x: number): number {
16 let sum = 1, term = 1;
17 const y = (x * x) / 4;
18 for (let k = 1; k < 60; k++) { term *= y / (k * k); sum += term; if (term < 1e-17 * sum) break; }
19 return sum;
20}
22/** Windowed-sinc taps that read x(j - i0 - fr), lags i0 - T/2 + 1 ... i0 + T/2 (Kaiser, beta 6, unit sum). */
23function fractionalDelayTaps(fr: number, T = DELAY_TAPS): Float64Array {
24 const taps = new Float64Array(T), half = T / 2, i0 = besselI0(6);
25 let sum = 0;
26 for (let t = 0; t < T; t++) {
27 const k = t - half + 1, x = k - fr;
28 const sinc = x === 0 ? 1 : Math.sin(Math.PI * x) / (Math.PI * x);
29 const u = x / half, win = besselI0(6 * Math.sqrt(Math.max(1 - u * u, 0))) / i0;
30 taps[t] = sinc * win; sum += taps[t];
31 }
32 for (let t = 0; t < T; t++) taps[t] /= sum;
33 return taps;
34}
36interface Geometry {
37 E: number;
38 coefs: Float32Array<ArrayBuffer>; // [2e] front rho A / (2 pi r), [2e+1] back with the diffraction gain
39 delays: Uint32Array<ArrayBuffer>; // [2e], [2e+1]: integer delays in audio samples
40 taps: Float32Array<ArrayBuffer>; // per element 32 floats: 16 front taps then 16 back taps
41 maxLag: number; // largest lag any tap reaches back, in samples
42}
44/** Per-element coefficients, delays and interpolator taps for a listener at ``position``. */
45function buildGeometry(run: CymsimRun, position: [number, number, number]): Geometry {
46 const m = run.mesh, E = m.nElements, [Lx, Ly, Lz] = position;
47 const { c_air: c, rho_air: rho, back_diffraction_gain: g } = run.config.listener;
48 // Diffraction ring: the boundary loop, subdivided in polar coordinates to about 720 points so
49 // that the min-over-ring path length matches cymsimweb's analytic rim (which uses nRim = 720).
50 const nb = m.boundary.length, sub = Math.max(1, Math.ceil(720 / nb));
51 const rim: number[][] = [];
52 for (let k = 0; k < nb; k++) {
53 const a = m.boundary[k], b = m.boundary[(k + 1) % nb];
54 const ra = Math.hypot(m.xy[2 * a], m.xy[2 * a + 1]), rb = Math.hypot(m.xy[2 * b], m.xy[2 * b + 1]);
55 const pa = Math.atan2(m.xy[2 * a + 1], m.xy[2 * a]);
56 let pb = Math.atan2(m.xy[2 * b + 1], m.xy[2 * b]);
57 if (pb - pa > Math.PI) pb -= 2 * Math.PI;
58 if (pa - pb > Math.PI) pb += 2 * Math.PI;
59 for (let s = 0; s < sub; s++) {
60 const t = s / sub, r = ra + t * (rb - ra), phi = pa + t * (pb - pa), z = m.z0[a] + t * (m.z0[b] - m.z0[a]);
61 rim.push([r * Math.cos(phi), r * Math.sin(phi), z]);
62 }
63 }
64 const rimToL = rim.map(([x, y, z]) => Math.hypot(x - Lx, y - Ly, z - Lz));
65 const coefs = new Float32Array(2 * E), delays = new Uint32Array(2 * E), taps = new Float32Array(32 * E);
66 let maxLag = 0;
67 for (let e = 0; e < E; e++) {
68 const a = m.triangles[3 * e], b = m.triangles[3 * e + 1], d = m.triangles[3 * e + 2];
69 const cx = (m.xy[2 * a] + m.xy[2 * b] + m.xy[2 * d]) / 3;
70 const cy = (m.xy[2 * a + 1] + m.xy[2 * b + 1] + m.xy[2 * d + 1]) / 3;
71 const cz = (m.z0[a] + m.z0[b] + m.z0[d]) / 3;
72 const rF = Math.hypot(Lx - cx, Ly - cy, Lz - cz);
73 let rB = Infinity;
74 for (let k = 0; k < rim.length; k++) {
75 rB = Math.min(rB, Math.hypot(cx - rim[k][0], cy - rim[k][1], cz - rim[k][2]) + rimToL[k]);
76 }
77 coefs[2 * e] = (rho * m.area[e]) / (2 * Math.PI * rF);
78 coefs[2 * e + 1] = (g * rho * m.area[e]) / (2 * Math.PI * rB);
79 const dF = (rF / c) * run.fs, dB = (rB / c) * run.fs;
80 delays[2 * e] = Math.floor(dF); delays[2 * e + 1] = Math.floor(dB);
81 taps.set(fractionalDelayTaps(dF - delays[2 * e]), 32 * e);
82 taps.set(fractionalDelayTaps(dB - delays[2 * e + 1]), 32 * e + 16);
83 maxLag = Math.max(maxLag, delays[2 * e], delays[2 * e + 1]);
84 }
85 return { E, coefs, delays, taps, maxLag: maxLag + DELAY_TAPS / 2 + 1 };
86}
88/**
89 * Element means of the decoded record for global samples [j0, j1), sample-major (E per row);
90 * rows for samples before 0 are zero (nothing has radiated yet).
91 */
92function elementMeans(run: CymsimRun, j0: number, j1: number): Float32Array<ArrayBuffer> {
93 const acc = run.acceleration!, m = run.mesh, E = m.nElements, N = acc.nodes;
94 const out = new Float32Array((j1 - j0) * E);
95 const d0 = Math.max(0, j0);
96 if (d0 >= j1) return out;
97 const rows = acc.read(d0, j1);
98 for (let j = d0; j < j1; j++) {
99 const src = (j - d0) * N, dst = (j - j0) * E;
100 for (let e = 0; e < E; e++) {
101 const a = m.triangles[3 * e], b = m.triangles[3 * e + 1], c = m.triangles[3 * e + 2];
102 out[dst + e] = (rows[src + a] + rows[src + b] + rows[src + c]) / 3;
103 }
104 }
105 return out;
106}
108export interface RenderedListener { front: Float32Array; back: Float32Array; position: [number, number, number]; backend: "gpu" | "cpu" }
109export interface RenderOptions { onProgress?: (fraction: number, backend: "gpu" | "cpu") => void; shouldStop?: () => boolean }
111const GATHER_WGSL = /* wgsl */ `
112struct U { E: u32, D: u32, chunkLen: u32, base: u32 }
113@group(0) @binding(0) var<uniform> u: U;
114@group(0) @binding(1) var<storage, read> coefs: array<f32>;
115@group(0) @binding(2) var<storage, read> delays: array<u32>;
116@group(0) @binding(3) var<storage, read> taps: array<f32>;
117@group(0) @binding(4) var<storage, read> win: array<f32>; // (D + chunkLen) rows of E means
118@group(0) @binding(5) var<storage, read_write> outp: array<f32>; // 2 per sample: front, back
120const WG = 256u;
121var<workgroup> pf: array<f32, WG>;
122var<workgroup> pb: array<f32, WG>;
124@compute @workgroup_size(WG)
125fn gather(@builtin(workgroup_id) wg: vec3<u32>, @builtin(local_invocation_id) li: vec3<u32>) {
126 let s = wg.x; let tid = li.x;
127 var f = 0.0; var b = 0.0;
128 if (s < u.chunkLen) {
129 let jg = i32(u.base + s);
130 for (var e = tid; e < u.E; e += WG) {
131 let i0f = i32(delays[2u * e]); let i0b = i32(delays[2u * e + 1u]);
132 var af = 0.0; var ab = 0.0;
133 for (var t = 0; t < 16; t++) {
134 let lagF = i0f + t - 7;
135 if (lagF >= 0 && lagF <= jg) {
136 af += taps[e * 32u + u32(t)] * win[u32(i32(u.D + s) - lagF) * u.E + e];
137 }
138 let lagB = i0b + t - 7;
139 if (lagB >= 0 && lagB <= jg) {
140 ab += taps[e * 32u + 16u + u32(t)] * win[u32(i32(u.D + s) - lagB) * u.E + e];
141 }
142 }
143 f += coefs[2u * e] * af; b += coefs[2u * e + 1u] * ab;
144 }
145 }
146 pf[tid] = f; pb[tid] = b;
147 workgroupBarrier();
148 var step = WG / 2u;
149 while (step > 0u) {
150 if (tid < step) { pf[tid] += pf[tid + step]; pb[tid] += pb[tid + step]; }
151 workgroupBarrier();
152 step = step / 2u;
153 }
154 if (tid == 0u && s < u.chunkLen) { outp[2u * s] = pf[0]; outp[2u * s + 1u] = pb[0]; }
155}
156`;
158let devicePromise: Promise<GPUDevice | null> | null = null;
159function gpuDevice(): Promise<GPUDevice | null> {
160 devicePromise ??= (async () => {
161 try {
162 const adapter = await navigator.gpu?.requestAdapter();
163 return adapter ? await adapter.requestDevice() : null;
164 } catch { return null; }
165 })();
166 return devicePromise;
167}
169async function renderGpu(device: GPUDevice, run: CymsimRun, geom: Geometry, opts: RenderOptions): Promise<RenderedListener | null> {
170 const T = run.acceleration!.nSamples, E = geom.E, D = geom.maxLag;
171 const chunk = Math.max(256, Math.min(16384, Math.floor(48e6 / (4 * E)) - D));
172 const buf = (data: Float32Array<ArrayBuffer> | Uint32Array<ArrayBuffer>) => {
173 const b = device.createBuffer({ size: data.byteLength, usage: GPUBufferUsage.STORAGE | GPUBufferUsage.COPY_DST });
174 device.queue.writeBuffer(b, 0, data);
175 return b;
176 };
177 const coefs = buf(geom.coefs), delays = buf(geom.delays), taps = buf(geom.taps);
178 const win = device.createBuffer({ size: (D + chunk) * E * 4, usage: GPUBufferUsage.STORAGE | GPUBufferUsage.COPY_DST });
179 const outp = device.createBuffer({ size: 2 * chunk * 4, usage: GPUBufferUsage.STORAGE | GPUBufferUsage.COPY_SRC });
180 const stage = device.createBuffer({ size: 2 * chunk * 4, usage: GPUBufferUsage.MAP_READ | GPUBufferUsage.COPY_DST });
181 const uni = device.createBuffer({ size: 16, usage: GPUBufferUsage.UNIFORM | GPUBufferUsage.COPY_DST });
182 const pipe = device.createComputePipeline({
183 layout: "auto",
184 compute: { module: device.createShaderModule({ code: GATHER_WGSL }), entryPoint: "gather" },
185 });
186 const bind = device.createBindGroup({
187 layout: pipe.getBindGroupLayout(0),
188 entries: [uni, coefs, delays, taps, win, outp].map((b, i) => ({ binding: i, resource: { buffer: b } })),
189 });
190 const front = new Float32Array(T), back = new Float32Array(T);
191 try {
192 for (let j0 = 0; j0 < T; j0 += chunk) {
193 if (opts.shouldStop?.()) return null;
194 const len = Math.min(chunk, T - j0);
195 device.queue.writeBuffer(uni, 0, new Uint32Array([E, D, len, j0]));
196 device.queue.writeBuffer(win, 0, elementMeans(run, j0 - D, j0 + len));
197 const enc = device.createCommandEncoder();
198 const pass = enc.beginComputePass();
199 pass.setPipeline(pipe); pass.setBindGroup(0, bind); pass.dispatchWorkgroups(len);
200 pass.end();
201 enc.copyBufferToBuffer(outp, 0, stage, 0, 2 * len * 4);
202 device.queue.submit([enc.finish()]);
203 await stage.mapAsync(GPUMapMode.READ, 0, 2 * len * 4);
204 const rows = new Float32Array(stage.getMappedRange(0, 2 * len * 4)).slice();
205 stage.unmap();
206 for (let s = 0; s < len; s++) { front[j0 + s] = rows[2 * s]; back[j0 + s] = rows[2 * s + 1]; }
207 opts.onProgress?.((j0 + len) / T, "gpu");
208 }
209 } finally {
210 for (const b of [coefs, delays, taps, win, outp, stage, uni]) b.destroy();
211 }
212 return { front, back, position: [0, 0, 0], backend: "gpu" };
213}
215async function renderCpu(run: CymsimRun, geom: Geometry, opts: RenderOptions): Promise<RenderedListener | null> {
216 const T = run.acceleration!.nSamples, E = geom.E, D = geom.maxLag;
217 const chunk = Math.max(256, Math.min(8192, Math.floor(16e6 / (4 * E)) - D));
218 const front = new Float32Array(T), back = new Float32Array(T);
219 const { coefs, delays, taps } = geom;
220 for (let j0 = 0; j0 < T; j0 += chunk) {
221 if (opts.shouldStop?.()) return null;
222 const len = Math.min(chunk, T - j0);
223 const win = elementMeans(run, j0 - D, j0 + len);
224 for (let s = 0; s < len; s++) {
225 const jg = j0 + s;
226 let f = 0, b = 0;
227 for (let e = 0; e < E; e++) {
228 const i0f = delays[2 * e], i0b = delays[2 * e + 1];
229 let af = 0, ab = 0;
230 for (let t = 0; t < 16; t++) {
231 const lagF = i0f + t - 7;
232 if (lagF >= 0 && lagF <= jg) af += taps[32 * e + t] * win[(D + s - lagF) * E + e];
233 const lagB = i0b + t - 7;
234 if (lagB >= 0 && lagB <= jg) ab += taps[32 * e + 16 + t] * win[(D + s - lagB) * E + e];
235 }
236 f += coefs[2 * e] * af; b += coefs[2 * e + 1] * ab;
237 }
238 front[jg] = f; back[jg] = b;
239 }
240 opts.onProgress?.((j0 + len) / T, "cpu");
241 await new Promise((r) => setTimeout(r, 0)); // let the page breathe between chunks
242 }
243 return { front, back, position: [0, 0, 0], backend: "cpu" };
244}
246/** Render the baffled (front) and back signals at ``position``; null when stopped or no record. */
247export async function renderListener(run: CymsimRun, position: [number, number, number], opts: RenderOptions = {}): Promise<RenderedListener | null> {
248 if (!run.acceleration || run.acceleration.nSamples === 0) return null;
249 const geom = buildGeometry(run, position);
250 const device = await gpuDevice();
251 const out = device ? await renderGpu(device, run, geom, opts) : await renderCpu(run, geom, opts);
252 if (out) out.position = [...position];
253 return out;
254}