1/**
2 * The room recording: every pattern-to-mic impulse response, from one run.
3 *
4 * The fast render path (src/studio/modal.ts) needs, per coupling pattern —
5 * a membrane mode's prof .* rad(phi_mn), a string or bar mode's
6 * prof .* rad(phi_m), a string's bridge patch bprof — the air's response at
7 * the microphone to an impulse of that pattern. Recording them one run per
8 * pattern would cost thousands of runs; discrete reciprocity collapses them
9 * to one. The air update's transfer matrix is symmetric — the masked
10 * Laplacian in divergence form has w_i·w_j/h² off-diagonal entries and
11 * everything else is diagonal — so the response at the mic to pattern g
12 * equals the projection onto g of the field one impulse *at the mic*
13 * produces. (This needs the uniform sound speed every generated scene has;
14 * scripts/modal-poc.ts checks the symmetry to machine precision.) So: run
15 * the air alone, kicked at the microphone, and after every step project the
16 * field onto every pattern at once.
17 *
18 * The air runs through the same pipeline as everything else — an air-only .m
19 * compiled by the same compiler against the same lapw op, so the fast path's
20 * numerics are the real solver's, not a reimplementation. The impulse is the
21 * init: kicking the update with source s at step 0 produces
22 * p[1] = dt²·wall·s/(1 + sig·dt), which init writes directly.
23 *
24 * The projections ride each step's submission, the way the Recorder's mic
25 * tap does, organized as *pattern groups*: each group contracts the field
26 * against one 3D profile over z (the sheet, the tube, the bridge patch),
27 * then against per-axis row tables sx and sy — rad of a product of sines
28 * factorizes per axis, a line object's y factor is a row of ones, and a
29 * fixed pattern like the bridge is a single (ones, ones) pair. Kernels
30 * accumulate on the GPU in chunks and come back to the CPU between chunks.
31 */
32import type { AirGrid } from '../grids.ts';
33import { GpuModel } from '../mgpu/model.ts';
34import { makeLapw, OpPlan } from '../mgpu/ops.ts';
35import { WORKGROUP_SIZE } from '../mgpu/wgsl.ts';
37/** The air-only model the recording runs: the generated step's air lines,
38 * with the reciprocal impulse as the initial condition. */
39const roomModel = (): string => `% Generated by soundbox for the room recording: the air alone, kicked at
40% the microphone. The step is the generated model's air update, verbatim.
41function [p, pm] = init(npts, wall, sig, imp, dt)
42 p = (dt*dt) * ((wall .* imp) ./ (1 + sig*dt));
43 pm = zeros(npts, 1);
44end
46function [pn, po] = step(p, pm, c, sig, wall, dt)
47 lp = lapw(p, wall);
48 sd = sig * dt;
49 pn = wall .* ((2*p - (1 - sd) .* pm + (c*dt).^2 .* lp) ./ (1 + sd));
50 po = p;
51end
52`;
54/** Thrown when a recording is cancelled through `RoomOptions.cancelled`. */
55export class RoomCancelled extends Error {
56 constructor() {
57 super('room recording cancelled');
58 }
59}
61/** One family of patterns sharing a 3D profile: each pair (xi, yi) is the
62 * pattern prof(x,y,z) · sx[xi](ix) · sy[yi](iy), and produces one kernel
63 * row. */
64export interface RoomPatternGroup {
65 label: string;
66 /** The 3D profile, air-sized; contracted over z first. */
67 prof: Float32Array;
68 /** Row tables: sx is sxRows × nx, sy is syRows × ny. */
69 sx: Float32Array;
70 sxRows: number;
71 sy: Float32Array;
72 syRows: number;
73 /** The (sx row, sy row) of each kernel row this group records. */
74 pairs: { xi: number; yi: number }[];
75}
77export interface RoomOptions {
78 device: GPUDevice;
79 air: AirGrid;
80 /** The scene's evaluated medium. */
81 fields: { c: Float32Array; sig: Float32Array; wall: Float32Array };
82 dt: number;
83 /** Flat air index of the microphone cell. */
84 micIndex: number;
85 groups: RoomPatternGroup[];
86 /** Kernel length, steps. */
87 steps: number;
88 onProgress?: (done: number, total: number) => void;
89 /** Polled between batches; returning true abandons the recording (the
90 * buffers are torn down and RoomCancelled is thrown). */
91 cancelled?: () => boolean;
92}
94export interface Room {
95 /** One row of `steps` samples per pattern, row-major: the groups'
96 * `pairs`, concatenated in group order. */
97 kernels: Float32Array;
98 steps: number;
99 rows: number;
100}
102const planeShader = (air: AirGrid): string => `
103@group(0) @binding(0) var<storage, read_write> plane: array<f32>;
104@group(0) @binding(1) var<storage, read> field: array<f32>;
105@group(0) @binding(2) var<storage, read> prof: array<f32>;
106@compute @workgroup_size(${WORKGROUP_SIZE})
107fn main(@builtin(global_invocation_id) gid: vec3<u32>) {
108 let c = gid.x;
109 if (c >= ${air.nx * air.ny}u) { return; }
110 var s = 0.0;
111 for (var iz = 0u; iz < ${air.nz}u; iz++) {
112 let i = c + ${air.nx * air.ny}u * iz;
113 s += prof[i] * field[i];
114 }
115 plane[c] = s;
116}
117`;
119const xcShader = (air: AirGrid, rows: number): string => `
120@group(0) @binding(0) var<storage, read_write> tmp: array<f32>;
121@group(0) @binding(1) var<storage, read> plane: array<f32>;
122@group(0) @binding(2) var<storage, read> sx: array<f32>;
123@compute @workgroup_size(${WORKGROUP_SIZE})
124fn main(@builtin(global_invocation_id) gid: vec3<u32>) {
125 let k = gid.x;
126 if (k >= ${rows * air.ny}u) { return; }
127 let m = k % ${rows}u;
128 let iy = k / ${rows}u;
129 var s = 0.0;
130 for (var ix = 0u; ix < ${air.nx}u; ix++) {
131 s += sx[m * ${air.nx}u + ix] * plane[ix + ${air.nx}u * iy];
132 }
133 tmp[k] = s;
134}
135`;
137const projShader = (
138 air: AirGrid,
139 sxRows: number,
140 nPairs: number,
141 rowOffset: number,
142 totalRows: number,
143 chunk: number,
144): string => `
145@group(0) @binding(0) var<storage, read_write> chunkBuf: array<f32>;
146@group(0) @binding(1) var<storage, read> tmp: array<f32>;
147@group(0) @binding(2) var<storage, read> sy: array<f32>;
148@group(0) @binding(3) var<storage, read> pairs: array<u32>; // xi | yi<<16
149@group(0) @binding(4) var<storage, read> head: array<u32>;
150@compute @workgroup_size(${WORKGROUP_SIZE})
151fn main(@builtin(global_invocation_id) gid: vec3<u32>) {
152 let q = gid.x;
153 if (q >= ${nPairs}u) { return; }
154 let mn = pairs[q];
155 let xi = mn & 0xffffu;
156 let yi = mn >> 16u;
157 var s = 0.0;
158 for (var iy = 0u; iy < ${air.ny}u; iy++) {
159 s += sy[yi * ${air.ny}u + iy] * tmp[xi + ${sxRows}u * iy];
160 }
161 let i = head[0];
162 if (i < ${chunk}u) {
163 chunkBuf[i * ${totalRows}u + ${rowOffset}u + q] = s;
164 }
165}
166`;
168const bumpShader = `
169@group(0) @binding(0) var<storage, read_write> head: array<u32>;
170@compute @workgroup_size(1)
171fn main() {
172 head[0] = head[0] + 1u;
173}
174`;
176/** Run the reciprocal recording and return every pattern's kernel. */
177export async function recordRoom(opts: RoomOptions): Promise<Room> {
178 const { device, air, fields, dt, micIndex, groups, steps } = opts;
179 const rows = groups.reduce((a, g) => a + g.pairs.length, 0);
180 if (rows === 0) return { kernels: new Float32Array(0), steps, rows };
181 for (const g of groups) {
182 for (const p of g.pairs) {
183 if (p.xi < 0 || p.xi >= g.sxRows || p.yi < 0 || p.yi >= g.syRows) {
184 throw new Error(`room group '${g.label}': pair (${p.xi},${p.yi}) outside its tables`);
185 }
186 }
187 }
189 // The air-only model, through the ordinary pipeline.
190 const imp = new Float32Array(air.npts);
191 imp[micIndex] = 1;
192 const gpu = await GpuModel.create({
193 device,
194 ops: new OpPlan(device, [makeLapw('lapw', air)]),
195 consts: { npts: air.npts, nx: air.nx, ny: air.ny, nz: air.nz },
196 fields: { c: fields.c, sig: fields.sig, wall: fields.wall, imp },
197 source: roomModel(),
198 paramNames: [],
199 state: [
200 { name: 'p', count: air.npts },
201 { name: 'pm', count: air.npts },
202 ],
203 });
204 gpu.setDt(dt);
205 gpu.setParams({});
207 // Kernels come back in chunks of whole samples; each sample is one f32 per
208 // pattern row. Cap a chunk at 64 MB.
209 const chunk = Math.max(1, Math.min(steps, Math.floor((64 << 20) / (4 * rows))));
210 const owned: GPUBuffer[] = [];
211 const mk = (label: string, size: number): GPUBuffer => {
212 const b = device.createBuffer({
213 label: `room-${label}`,
214 size,
215 usage: GPUBufferUsage.STORAGE | GPUBufferUsage.COPY_SRC | GPUBufferUsage.COPY_DST,
216 });
217 owned.push(b);
218 return b;
219 };
220 try {
221 const headBuf = mk('head', 4);
222 const chunkBuf = mk('chunk', 4 * rows * chunk);
223 const readback = device.createBuffer({
224 label: 'room-readback',
225 size: 4 * rows * chunk,
226 usage: GPUBufferUsage.MAP_READ | GPUBufferUsage.COPY_DST,
227 });
228 owned.push(readback);
229 device.queue.writeBuffer(headBuf, 0, new Uint32Array([0]));
231 const pipe = async (code: string, label: string): Promise<GPUComputePipeline> =>
232 device.createComputePipelineAsync({
233 label: `room-${label}`,
234 layout: 'auto',
235 compute: {
236 module: device.createShaderModule({ code, label: `room-${label}` }),
237 entryPoint: 'main',
238 },
239 });
240 const bind = (p: GPUComputePipeline, buffers: GPUBuffer[]): GPUBindGroup =>
241 device.createBindGroup({
242 layout: p.getBindGroupLayout(0),
243 entries: buffers.map((buffer, i) => ({ binding: i, resource: { buffer } })),
244 });
245 const upload = (label: string, data: Float32Array | Uint32Array): GPUBuffer => {
246 const b = mk(label, data.byteLength);
247 device.queue.writeBuffer(b, 0, data as Float32Array<ArrayBuffer>);
248 return b;
249 };
251 const pField = gpu.stateBuffer('p')!;
252 const planePipe = await pipe(planeShader(air), 'plane');
253 const bumpPipe = await pipe(bumpShader, 'bump');
254 const bumpBind = bind(bumpPipe, [headBuf]);
256 interface GroupRun {
257 xcPipe: GPUComputePipeline;
258 projPipe: GPUComputePipeline;
259 planeBind: GPUBindGroup;
260 xcBind: GPUBindGroup;
261 projBind: GPUBindGroup;
262 xcThreads: number;
263 nPairs: number;
264 }
265 const runs: GroupRun[] = [];
266 let rowOffset = 0;
267 for (const g of groups) {
268 const planeBuf = mk(`plane-${g.label}`, 4 * air.nx * air.ny);
269 const tmpBuf = mk(`tmp-${g.label}`, 4 * g.sxRows * air.ny);
270 const profBuf = upload(`prof-${g.label}`, g.prof);
271 const sxBuf = upload(`sx-${g.label}`, g.sx);
272 const syBuf = upload(`sy-${g.label}`, g.sy);
273 const pairsBuf = upload(
274 `pairs-${g.label}`,
275 new Uint32Array(g.pairs.map((p) => p.xi | (p.yi << 16))),
276 );
277 const xcPipe = await pipe(xcShader(air, g.sxRows), `xc-${g.label}`);
278 const projPipe = await pipe(
279 projShader(air, g.sxRows, g.pairs.length, rowOffset, rows, chunk),
280 `proj-${g.label}`,
281 );
282 runs.push({
283 xcPipe,
284 projPipe,
285 planeBind: bind(planePipe, [planeBuf, pField, profBuf]),
286 xcBind: bind(xcPipe, [tmpBuf, planeBuf, sxBuf]),
287 projBind: bind(projPipe, [chunkBuf, tmpBuf, syBuf, pairsBuf, headBuf]),
288 xcThreads: g.sxRows * air.ny,
289 nPairs: g.pairs.length,
290 });
291 rowOffset += g.pairs.length;
292 }
294 const encodeProjection = (enc: GPUCommandEncoder): void => {
295 const pass = enc.beginComputePass({ label: 'room-projection' });
296 const go = (p: GPUComputePipeline, b: GPUBindGroup, threads: number): void => {
297 pass.setPipeline(p);
298 pass.setBindGroup(0, b);
299 pass.dispatchWorkgroups(Math.ceil(threads / WORKGROUP_SIZE));
300 };
301 for (const r of runs) {
302 go(planePipe, r.planeBind, air.nx * air.ny);
303 go(r.xcPipe, r.xcBind, r.xcThreads);
304 go(r.projPipe, r.projBind, r.nPairs);
305 }
306 go(bumpPipe, bumpBind, 1);
307 pass.end();
308 };
310 // Sample 0 is the projection of the post-init field (the kicked air, one
311 // step in); every later sample rides its step's own submission.
312 const kernels = new Float32Array(rows * steps);
313 gpu.init();
314 {
315 const enc = device.createCommandEncoder({ label: 'room-first-sample' });
316 encodeProjection(enc);
317 device.queue.submit([enc.finish()]);
318 }
319 let got = 0; // samples already on the CPU
320 let pending = 1; // samples sitting in the chunk buffer
321 const flush = async (): Promise<void> => {
322 const enc = device.createCommandEncoder({ label: 'room-flush' });
323 enc.copyBufferToBuffer(chunkBuf, 0, readback, 0, 4 * rows * pending);
324 device.queue.submit([enc.finish()]);
325 await readback.mapAsync(GPUMapMode.READ, 0, 4 * rows * pending);
326 const data = new Float32Array(readback.getMappedRange(0, 4 * rows * pending));
327 // Chunk samples are row-minor; kernels are row-major.
328 for (let s = 0; s < pending; s++)
329 for (let q = 0; q < rows; q++) kernels[q * steps + got + s] = data[s * rows + q];
330 readback.unmap();
331 device.queue.writeBuffer(headBuf, 0, new Uint32Array([0]));
332 got += pending;
333 pending = 0;
334 };
335 // Steps per submission, as the plain renderer batches; the sync between
336 // batches keeps the queue from running unboundedly ahead, and is what
337 // paces the progress reports and the cancellation checks.
338 const batch = 2048;
339 while (got + pending < steps) {
340 if (opts.cancelled?.()) throw new RoomCancelled();
341 const n = Math.min(batch, chunk - pending, steps - got - pending);
342 gpu.step(n, encodeProjection);
343 pending += n;
344 await device.queue.onSubmittedWorkDone();
345 opts.onProgress?.(got + pending, steps);
346 if (pending === chunk) await flush();
347 }
348 if (pending > 0) await flush();
349 return { kernels, steps, rows };
350 } finally {
351 for (const b of owned) b.destroy();
352 gpu.destroy();
353 }
354}