/** * The instrument author: writes complete instrument files for the three v1 * instruments (drum, dulcimer, chime) from one parameterized builder. * * The physics and its numbers are soundbox's, taken from the same equations * the numbl pipeline compiles there (the test suite holds the two projects * together; see scripts/test-node.ts). What this builder adds is the form: * every grid size and parameter reaches the WGSL as a pipeline override * constant, so one file covers every setting of every parameter, and the * whole precompute-vs-tuning split is carried by the parameter classes in * the file rather than by host code that knows what a drum is. * * The precompute is the reciprocal room recording: run the air alone, * kicked at the microphone, and after every step project the field onto * each mode's coupling pattern (a product of interpolated sines under the * object's radiation profile; for a string, one more pattern for its * bridge patch). Discrete reciprocity makes each projection row that * pattern's impulse response at the microphone. The render solves every * object mode in closed form (the leapfrog restricted to a sine eigenmode * is a two-term scalar recursion) and filters each kernel with a complex * one-pole per root, one GPU thread per (root, kernel) pair, tiled over * time with a reduction per tile. * * The recorded mode set is geometric: sine modes up to `mode_margin` times * what the air grid can spatially resolve across the object * (m <= mode_margin * L / (2h)). Nothing about the strike or the tuning * enters, so the bake depends on geometry and accuracy parameters alone. */ import type { InstrumentFile, MarkerSpec, ParamSpec, PassSpec, VisualSpec } from '../src/format.ts'; const WG = 64; /** Thread index linearized over a possibly two-dimensional workgroup grid * (the x dimension caps at 65535 workgroups). */ const IDX = 'let i = gid.x + gid.y * 4194240u;'; const MAIN = `@compute @workgroup_size(${WG}) fn main(@builtin(global_invocation_id) gid: vec3) { ${IDX}`; const ov = (names: [string, 'u32' | 'f32'][]): string => names.map(([n, t]) => `override ${n}: ${t};`).join('\n'); const buf = (bindings: [string, 'read' | 'read_write'][]): string => bindings .map( ([n, m], i) => `@group(0) @binding(${i}) var ${n === '__tick' ? 'tick' : n === '__events' ? 'ev' : n}: array<${n === '__tick' ? 'u32' : 'f32'}>;`, ) .join('\n'); export type ObjectKind = 'membrane' | 'string' | 'bar'; export type BodyKind = 'shell' | 'box' | 'none'; export interface BuildOptions { name: string; description: string; object: ObjectKind; body: BodyKind; /** Geometry parameter defaults overriding the built-in ones. */ defaults?: Record; } /* ------------------------------------------------------------- params --- */ const P = ( key: string, label: string, unit: string, value: number, min: number, max: number, cls: 'geometry' | 'accuracy' | 'tuning', hint?: string, ): ParamSpec => ({ key, label, unit, value, min, max, class: cls, hint }); const bodyParams: Record = { none: [], shell: [ P('shellr', 'shell radius', 'm', 0.17, 0.08, 0.3, 'geometry', 'A hair larger than the head over it.'), P('shelld', 'shell depth', 'm', 0.3, 0.05, 0.4, 'geometry', 'The cavity under the head; deeper is a lower air resonance.'), P('thick', 'wall thickness', 'm', 0.016, 0.008, 0.05, 'geometry'), P('closed', 'closed bottom', '1', 1, 0, 1, 'geometry', '1 puts a plate across the bottom; 0 leaves it open.'), P('absorb', 'surface absorption', '1/s', 800, 0, 4000, 'geometry', 'A lossy skin over the shell surfaces; 0 lets the cavity ring hard.'), ], box: [ P('boxl', 'box length', 'm', 0.66, 0.3, 0.9, 'geometry'), P('boxw', 'box width', 'm', 0.16, 0.08, 0.4, 'geometry'), P('boxd', 'box depth', 'm', 0.08, 0.03, 0.25, 'geometry', 'How far the box hangs below its top plate at z = 0.'), P('thick', 'wall thickness', 'm', 0.016, 0.008, 0.05, 'geometry'), P('holer', 'sound hole radius', 'm', 0.035, 0, 0.12, 'geometry', '0 seals the box; the hole is what lets the cavity speak.'), P('holex', 'sound hole x', 'm', -0.05, -0.3, 0.3, 'geometry'), P('absorb', 'surface absorption', '1/s', 800, 0, 4000, 'geometry'), ], }; const objectParams: Record = { membrane: { geometry: [ P('L1', 'head width', 'm', 0.32, 0.15, 0.5, 'geometry', 'Extent along x.'), P('L2', 'head depth', 'm', 0.32, 0.15, 0.5, 'geometry', 'Extent along y.'), P('oz', 'head height', 'm', 0.02, 0, 0.15, 'geometry', 'Height of the head plane above z = 0.'), ], tuning: [ P('f0', 'pitch', 'Hz', 110, 60, 300, 'tuning', 'The (1,1) fundamental; the tension follows from it and the size.'), P('t60', 'decay t60', 's', 1.5, 0.2, 6, 'tuning', 'Seconds to decay 60 dB.'), P('sig1', 'brightness decay', 'm^2/s', 0.0005, 0, 0.01, 'tuning', 'Frequency-dependent damping; the slap of the attack fades first.'), P('strikex', 'strike x', '1', 0.35, 0.1, 0.9, 'tuning', 'Where the stick lands along x, as a fraction.'), P('strikey', 'strike y', '1', 0.5, 0.1, 0.9, 'tuning'), P('striker', 'stick width', 'm', 0.04, 0.01, 0.1, 'tuning', 'Radius of the struck patch.'), P('vel', 'strike speed', 'm/s', 1.5, 0.1, 5, 'tuning', 'Linear equations: sets scale only.'), P('gair', 'radiation gain', '1', 1, 0, 2, 'tuning'), ], }, string: { geometry: [ P('L1', 'speaking length', 'm', 0.6, 0.3, 0.9, 'geometry'), P('oy', 'across', 'm', 0, -0.15, 0.15, 'geometry'), P('oz', 'height', 'm', 0.02, 0.005, 0.15, 'geometry', 'Height of the string above z = 0.'), ], tuning: [ P('f0', 'pitch', 'Hz', 294, 80, 880, 'tuning', 'The fundamental; the wave speed follows, cs = 2 L f0.'), P('B', 'inharmonicity', '1', 0.0001, 0, 0.002, 'tuning', 'Bending stiffness as the inharmonicity coefficient.'), P('t60', 'decay t60', 's', 4, 0.2, 8, 'tuning'), P('sig1', 'brightness decay', 'm^2/s', 0.005, 0, 0.05, 'tuning'), P('pluckpos', 'pluck position', '1', 0.22, 0.08, 0.92, 'tuning'), P('amp', 'pluck height', 'm', 0.002, 0.0002, 0.005, 'tuning'), P('gair', 'line radiation gain', '1', 0.3, 0, 2, 'tuning'), P('gbridge', 'bridge drive gain', '1', 1, 0, 2, 'tuning', 'The pull at the far end driving a patch above the plate.'), ], }, bar: { geometry: [ P('L1', 'length', 'm', 0.3, 0.15, 0.5, 'geometry'), P('oy', 'across', 'm', 0, -0.15, 0.15, 'geometry'), P('oz', 'height', 'm', 0.05, 0.005, 0.15, 'geometry'), ], tuning: [ P('f1', 'pitch', 'Hz', 780, 200, 2000, 'tuning', 'The first partial; higher ones land at n^2 f1.'), P('t60', 'decay t60', 's', 6, 0.5, 12, 'tuning'), P('sig1', 'brightness decay', 'm^2/s', 0.001, 0, 0.02, 'tuning'), P('strikepos', 'strike position', '1', 0.3, 0.1, 0.9, 'tuning'), P('striker', 'mallet width', 'm', 0.025, 0.01, 0.08, 'tuning'), P('vel', 'strike speed', 'm/s', 1, 0.1, 5, 'tuning'), P('gair', 'radiation gain', '1', 1, 0, 2, 'tuning'), ], }, }; /* ------------------------------------------------------------ derived --- */ function derivedOf(object: ObjectKind): Record { const air: Record = { c0: '343', cfl: '0.5', ny: 'nx / 2', nz: 'nx / 2', npts: 'nx * ny * nz', h: 'Lx / nx', Ly: 'ny * h', Lz: 'nz * h', dt: 'cfl * h / (sqrt(3) * c0)', ir_steps: 'max(1, round(ir_seconds / dt))', mic_index: 'clamp(round((mic_x + Lx/2) / h - 0.5), 0, nx - 1)' + ' + nx * (clamp(round((mic_y + Ly/2) / h - 0.5), 0, ny - 1)' + ' + ny * clamp(round((mic_z + Lz/2) / h - 0.5), 0, nz - 1))', }; if (object === 'membrane') { return { ...air, // The grid, sized for the worst case any tuning slider can reach, so // no tuning move ever invalidates it (soundbox's makeRectGrid). c_worst: '2 * f0__max / sqrt(1/(L1*L1) + 1/(L2*L2))', h_min: 'sqrt(2 * (c_worst*c_worst*dt*dt + 4*sig1__max*dt)) / 0.95', n1: 'max(5, min(256, floor(L1 / h_min) + 1))', n2: 'max(5, min(256, floor(L2 / h_min) + 1))', h1: 'L1 / (n1 - 1)', h2: 'L2 / (n2 - 1)', // The geometric mode set: what the air grid resolves, times a margin. mx: 'min(n1 - 2, ceil(mode_margin * L1 / (2 * h)))', my: 'min(n2 - 2, ceil(mode_margin * L2 / (2 * h)))', n_rows: 'mx * my', }; } const stiff: Record = object === 'string' ? { cs_worst: '2 * L1 * f0__max', kap2_worst: 'B__max * (cs_worst * cs_worst * L1 * L1) / (pi * pi)', } : { kap_worst: '2 * f1__max * (L1 * L1) / pi', kap2_worst: 'kap_worst * kap_worst', cs_worst: '0', }; return { ...air, ...stiff, a_worst: 'cs_worst*cs_worst*dt*dt + 4*sig1__max*dt', h_min: 'sqrt((a_worst + sqrt(a_worst*a_worst + 16*kap2_worst*dt*dt)) / 2) / 0.95', n1: 'max(5, min(1024, floor(L1 / h_min) + 1))', h1: 'L1 / (n1 - 1)', mx: 'n1 - 2', // line objects have few modes; record all of them n_rows: object === 'string' ? 'mx + 1' : 'mx', // + the bridge patch row }; } /* ---------------------------------------------------- medium (bodies) --- */ const bodyWgsl: Record = { none: { code: ' let shell = 0.0;', consts: [] }, shell: { code: ` // The shell: a cylindrical ring from z = -shelld to 0, optional bottom. let r = sqrt(x*x + y*y); let cyl_o = 0.5 * (1.0 - tanh((r - shellr) / sw)); let cyl_i = 0.5 * (1.0 - tanh((r - (shellr - thick)) / sw)); let cyl_z = 0.5 * (1.0 - tanh((abs(z + shelld/2.0) - shelld/2.0) / sw)); let bot = 0.5 * (1.0 - tanh((abs(z + shelld - thick/2.0) - thick/2.0) / sw)); let shell = min(1.0, max(0.0, cyl_o - cyl_i) * cyl_z + closed * (cyl_i * bot));`, consts: [['shellr', 'f32'], ['shelld', 'f32'], ['thick', 'f32'], ['closed', 'f32']], }, box: { code: ` // The box: the outer solid minus its cavity, top plate at z = 0, with a // round sound hole cut out of the top plate only. let box_ox = 0.5 * (1.0 - tanh((abs(x) - boxl/2.0) / sw)); let box_oy = 0.5 * (1.0 - tanh((abs(y) - boxw/2.0) / sw)); let box_oz = 0.5 * (1.0 - tanh((abs(z + boxd/2.0) - boxd/2.0) / sw)); let box_ix = 0.5 * (1.0 - tanh((abs(x) - (boxl/2.0 - thick)) / sw)); let box_iy = 0.5 * (1.0 - tanh((abs(y) - (boxw/2.0 - thick)) / sw)); let box_iz = 0.5 * (1.0 - tanh((abs(z + boxd/2.0) - (boxd/2.0 - thick)) / sw)); var shell = max(0.0, box_ox * box_oy * box_oz - box_ix * box_iy * box_iz); let hole_r = sqrt((x - holex)*(x - holex) + y*y); let hole = 0.5 * (1.0 - tanh((hole_r - holer) / sw)); let top = 0.5 * (1.0 - tanh((abs(z + thick/2.0) - thick/2.0) / sw)); shell = max(0.0, shell - hole * top);`, consts: [ ['boxl', 'f32'], ['boxw', 'f32'], ['boxd', 'f32'], ['thick', 'f32'], ['holer', 'f32'], ['holex', 'f32'], ], }, }; function mediumShader(object: ObjectKind, body: BodyKind): { code: string; constants: string[] } { const b = bodyWgsl[body]; const absorb = body === 'none' ? '' : ' + absorb * (4.0 * shell * (1.0 - shell))'; const profConsts: [string, 'f32'][] = object === 'membrane' ? [['oz', 'f32'], ['L1', 'f32'], ['L2', 'f32']] : [['oy', 'f32'], ['oz', 'f32'], ['L1', 'f32']]; const prof = object === 'membrane' ? ` // The head radiates from its plane: a thin sheet at z = oz over its // extent, cross-section normalized to one. let aw = 2.0 * h; let sk = exp(-((z - oz)*(z - oz)) / (aw*aw)) / (sqrt(PI) * aw); let exx = 0.5 * (1.0 + tanh((L1/2.0 - abs(x)) / (1.5*h))); let exy = 0.5 * (1.0 + tanh((L2/2.0 - abs(y)) / (1.5*h))); prof[i] = sk * exx * exy * w;` : ` // The object radiates from its line: a tube two cells wide around // (y = oy, z = oz), fading at its ends, cross-section normalized to one. let aw = 2.0 * h; let tube = exp(-((y - oy)*(y - oy) + (z - oz)*(z - oz)) / (aw*aw)) / (PI * aw * aw); let span = 0.5 * (1.0 + tanh((L1/2.0 - abs(x)) / (1.5*h))); prof[i] = tube * span * w;`; const bprof = object === 'string' ? ` // ...and its bridge force drives a patch of air just above the plate, // around the x = +L/2 end (a stand-in for a moving top plate). let pw = 1.5 * h; let pr = exp(-((x - L1/2.0)*(x - L1/2.0) + (y - oy)*(y - oy)) / (0.04*0.04)) / (PI * 0.04 * 0.04); let skin = exp(-(z - pw)*(z - pw) / (pw*pw)) / (sqrt(PI) * pw); bprof[i] = pr * skin * w;` : ''; const buffers: [string, 'read_write'][] = [ ['wall', 'read_write'], ['sig', 'read_write'], ['prof', 'read_write'], ]; if (object === 'string') buffers.push(['bprof', 'read_write']); const consts: [string, 'u32' | 'f32'][] = [ ['npts', 'u32'], ['nx', 'u32'], ['ny', 'u32'], ['h', 'f32'], ['Lx', 'f32'], ['Ly', 'f32'], ['Lz', 'f32'], ...b.consts, ...(body === 'none' ? [] : ([['absorb', 'f32']] as [string, 'f32'][])), ...profConsts, ]; const code = `// The medium: the body's rigid mask, the absorption, and where the object // couples into the air. One thread per air cell. const PI = 3.14159265358979; ${ov(consts)} ${buf(buffers)} ${MAIN} if (i >= npts) { return; } let ix = i % nx; let iy = (i / nx) % ny; let iz = i / (nx * ny); let x = -Lx/2.0 + (f32(ix) + 0.5) * h; let y = -Ly/2.0 + (f32(iy) + 0.5) * h; let z = -Lz/2.0 + (f32(iz) + 0.5) * h; let sw = 0.45 * h; ${b.code} let w = 1.0 - shell; wall[i] = w; // Open boundary at the domain edge (a quadratic absorbing ramp over the // outer 0.15 Ly of each face), plus the body's lossy surface skin. let spw = 0.15 * Ly; let dx = max(0.0, spw - (Lx/2.0 - abs(x))) / spw; let dy = max(0.0, spw - (Ly/2.0 - abs(y))) / spw; let dz = max(0.0, spw - (Lz/2.0 - abs(z))) / spw; let d = max(max(dx, dy), dz); sig[i] = 9000.0 * d * d${absorb}; ${prof} ${bprof} } `; return { code, constants: consts.map(([n]) => n) }; } /* ------------------------------------------- the rest of the precompute --- */ /** Interpolated sine tables: the pattern the rad coupling op applies is the * node sine linearly interpolated to each air column, so the projection * must use exactly that, not the continuous sine. */ function sinesShader(axis: 'x' | 'y'): { code: string; constants: string[] } { const [n, m, L, hobj, nair, Lair] = axis === 'x' ? ['n1', 'mx', 'L1', 'h1', 'nx', 'Lx'] : ['n2', 'my', 'L2', 'h2', 'ny', 'Ly']; const consts: [string, 'u32' | 'f32'][] = [ [n, 'u32'], [m, 'u32'], [nair, 'u32'], [L, 'f32'], [hobj, 'f32'], [Lair, 'f32'], ['h', 'f32'], ]; const code = `// Sine mode ${axis === 'x' ? 'm' : 'n'} of the object grid, interpolated to air // ${axis} centers: s${axis}[(mode-1) * ${nair} + i${axis}]. One thread per table entry. const PI = 3.14159265358979; ${ov(consts)} ${buf([['table', 'read_write']])} ${MAIN} if (i >= ${m} * ${nair}) { return; } let mi = i % ${m}; let ia = i / ${m}; let xw = -${Lair}/2.0 + (f32(ia) + 0.5) * h; let u = (xw + ${L}/2.0) / ${hobj}; let j = i32(floor(u)); var v = 0.0; if (j >= 0 && j < i32(${n}) - 1) { let f = u - f32(j); let mm = f32(mi + 1u); let s0 = sin(mm * PI * f32(j) / f32(${n} - 1u)); let s1 = sin(mm * PI * f32(j + 1) / f32(${n} - 1u)); v = (1.0 - f) * s0 + f * s1; } table[mi * ${nair} + ia] = v; } `; return { code, constants: consts.map(([q]) => q) }; } const impulseShader = { code: `// The reciprocal kick: one unit impulse at the microphone cell, entering // the way a source enters the first step. One thread. ${ov([['mic_index', 'u32'], ['dt', 'f32']])} ${buf([['p', 'read_write'], ['wall', 'read'], ['sig', 'read']])} @compute @workgroup_size(1) fn main() { p[mic_index] = dt * dt * wall[mic_index] / (1.0 + sig[mic_index] * dt); } `, constants: ['mic_index', 'dt'], }; const stepShader = { code: `// One leapfrog step of the air: the wall-masked Laplacian (each stencil // face scaled by the mask at the neighbour it reads, the rigid condition in // divergence form) fused with the update. One thread per air cell. ${ov([['npts', 'u32'], ['nx', 'u32'], ['ny', 'u32'], ['nz', 'u32'], ['h', 'f32'], ['dt', 'f32'], ['c0', 'f32']])} ${buf([['p', 'read'], ['pm', 'read'], ['pn', 'read_write'], ['wall', 'read'], ['sig', 'read']])} fn face(ix: i32, iy: i32, iz: i32, pc: f32) -> f32 { if (ix < 0 || ix >= i32(nx) || iy < 0 || iy >= i32(ny) || iz < 0 || iz >= i32(nz)) { return -pc; // outside the domain: p = 0, open } let k = u32(ix) + nx * (u32(iy) + ny * u32(iz)); return wall[k] * (p[k] - pc); } ${MAIN} if (i >= npts) { return; } let ix = i32(i % nx); let iy = i32((i / nx) % ny); let iz = i32(i / (nx * ny)); let pc = p[i]; let s = face(ix - 1, iy, iz, pc) + face(ix + 1, iy, iz, pc) + face(ix, iy - 1, iz, pc) + face(ix, iy + 1, iz, pc) + face(ix, iy, iz - 1, pc) + face(ix, iy, iz + 1, pc); let lp = wall[i] * s / (h * h); let sd = sig[i] * dt; pn[i] = wall[i] * ((2.0 * pc - (1.0 - sd) * pm[i] + (c0 * dt) * (c0 * dt) * lp) / (1.0 + sd)); } `, constants: ['npts', 'nx', 'ny', 'nz', 'h', 'dt', 'c0'], }; function planeShader(profName: string, outName: string): { code: string; constants: string[] } { return { code: `// Contract the field against the ${profName === 'bprof' ? 'bridge patch' : 'radiation'} profile over z: // one thread per air column. ${ov([['nx', 'u32'], ['ny', 'u32'], ['nz', 'u32']])} ${buf([[outName, 'read_write'], ['p', 'read'], [profName, 'read']])} ${MAIN} if (i >= nx * ny) { return; } var s = 0.0; for (var iz = 0u; iz < nz; iz++) { let k = i + nx * ny * iz; s += ${profName}[k] * p[k]; } ${outName}[i] = s; } `, constants: ['nx', 'ny', 'nz'], }; } const xcShader = { code: `// Contract the plane against every x sine row: tmp[m + mx * iy]. ${ov([['nx', 'u32'], ['ny', 'u32'], ['mx', 'u32']])} ${buf([['tmp', 'read_write'], ['plane', 'read'], ['sx', 'read']])} ${MAIN} if (i >= mx * ny) { return; } let m = i % mx; let iy = i / mx; var s = 0.0; for (var ix = 0u; ix < nx; ix++) { s += sx[m * nx + ix] * plane[ix + nx * iy]; } tmp[i] = s; } `, constants: ['nx', 'ny', 'mx'], }; /** The last contraction writes this step's sample of every kernel row. */ function projShader(object: ObjectKind): { code: string; constants: string[] } { if (object === 'membrane') { return { code: `// Finish each pattern's projection against its y sine row and write this // step's sample: kernels[q * ir_steps + t], q = (m-1) + mx*(n-1). ${ov([['ny', 'u32'], ['mx', 'u32'], ['my', 'u32'], ['ir_steps', 'u32']])} ${buf([['kernels', 'read_write'], ['tmp', 'read'], ['sy', 'read'], ['__tick', 'read']])} ${MAIN} if (i >= mx * my) { return; } let m = i % mx; let n = i / mx; var s = 0.0; for (var iy = 0u; iy < ny; iy++) { s += sy[n * ny + iy] * tmp[m + mx * iy]; } let t = tick[0]; if (t < ir_steps) { kernels[i * ir_steps + t] = s; } } `, constants: ['ny', 'mx', 'my', 'ir_steps'], }; } return { code: `// A line object's y factor is one everywhere: finish by summing tmp over // iy and write this step's sample of every kernel row. ${ov([['ny', 'u32'], ['mx', 'u32'], ['ir_steps', 'u32']])} ${buf([['kernels', 'read_write'], ['tmp', 'read'], ['__tick', 'read']])} ${MAIN} if (i >= mx) { return; } var s = 0.0; for (var iy = 0u; iy < ny; iy++) { s += tmp[i + mx * iy]; } let t = tick[0]; if (t < ir_steps) { kernels[i * ir_steps + t] = s; } } `, constants: ['ny', 'mx', 'ir_steps'], }; } const bridgeProjShader = { code: `// The bridge patch is a single fixed pattern: sum its plane and write the // sample into the last kernel row (row mx). One thread. ${ov([['nx', 'u32'], ['ny', 'u32'], ['mx', 'u32'], ['ir_steps', 'u32']])} ${buf([['kernels', 'read_write'], ['bplane', 'read'], ['__tick', 'read']])} @compute @workgroup_size(1) fn main() { var s = 0.0; for (var k = 0u; k < nx * ny; k++) { s += bplane[k]; } let t = tick[0]; if (t < ir_steps) { kernels[mx * ir_steps + t] = s; } } `, constants: ['nx', 'ny', 'mx', 'ir_steps'], }; /* ------------------------------------------------------------- render --- */ /** Shared complex and slot helpers for the slots shaders. A slot is one * characteristic root ready to filter one kernel row: * [cRe, cIm, zRe, zIm, row, start], six f32 each. */ const SLOT_HELPERS = ` fn cmul(ar: f32, ai: f32, br: f32, bi: f32) -> vec2f { return vec2f(ar * br - ai * bi, ar * bi + ai * br); } // c = a * gain * (z - 2 + 1/z) = a * gain * (z-1)^2 / z, the acceleration // series, with z-1 supplied stably (the naive z - 1 cancels). fn acc_scale(p: vec4f, zm1Re: f32, gain: f32) -> vec4f { let m2 = p.z * p.z + p.w * p.w; if (m2 == 0.0) { return vec4f(0.0); } let sq = cmul(zm1Re, p.w, zm1Re, p.w); let num = cmul(sq.x, sq.y, p.z, -p.w); let c = cmul(p.x, p.y, gain * num.x / m2, gain * num.y / m2); return vec4f(c.x, c.y, p.z, p.w); } // c = a * gain * z, the next-displacement series (what a string's bridge // slope is read from). fn next_scale(p: vec4f, gain: f32) -> vec4f { let c = cmul(p.x, p.y, gain * p.z, gain * p.w); return vec4f(c.x, c.y, p.z, p.w); } `; const WRITE_SLOT = ` fn write_slot(idx: u32, p: vec4f, row: f32, start: f32) { slots[idx * 6u] = p.x; slots[idx * 6u + 1u] = p.y; slots[idx * 6u + 2u] = p.z; slots[idx * 6u + 3u] = p.w; slots[idx * 6u + 4u] = row; slots[idx * 6u + 5u] = start; } `; /** * The mode's recursion solved in closed form, written to avoid the f32 * cancellations the naive discriminant carries. Restricted to eigenvalue * lam, the leapfrog is a[k+1] = cA a[k] + cB a[k-1] with * cA = (2 - x) / (1 + d), x = dt^2 w2 + 2 sig1 dt lam, * cB = (-(1 - d) + 2 sig1 dt lam) / (1 + d), d = (6.91 / t60) dt, * and the quantities the roots need reduce to * zRe = (1 - x/2)/(1+d), |z|^2 = (1 - g)/(1+d), g = d + 2 sig1 dt lam, * zIm^2 (1+d)^2 = N = dt^2 w2 - g d - x^2/4, * every term small and positive-definite where it matters, so the mode's * angle keeps full relative f32 precision (the naive cA^2 + 4 cB loses it). * N > 0 is the underdamped (audible) case; N <= 0 gives two real roots. * Excitation: a strike is a velocity (a[0] = 0, a[-1] = b), a pluck a held * displacement (a[0] = a[-1] = b). * * Leaves in scope: p0, p1 (displacement-series roots as [aRe, aIm, zRe, * zIm], zeroed when absent) and zm10, zm11 (each root's z - 1 real part, * stably). */ const rootsBlock = (excite: 'strike' | 'pluck'): string => ` let d = (6.91 / t60) * dt; let opd = 1.0 + d; let g = d + 2.0 * sig1 * dt * lam; let x = dt * dt * w2 + 2.0 * sig1 * dt * lam; let zRe = (1.0 - 0.5 * x) / opd; let zm1Re = -(0.5 * x + d) / opd; let nn = dt * dt * w2 - g * d - 0.25 * x * x; var p0 = vec4f(0.0); var p1 = vec4f(0.0); var zm10 = 0.0; var zm11 = 0.0; if (nn > 0.0) { let zIm = sqrt(nn) / opd; let mag2 = (1.0 - g) / opd; ${excite === 'strike' ? ` let aRe = 0.0; let aIm = b * mag2 / zIm;` : ` let aRe = b; let aIm = b * ((0.5 * x - g) / opd) / zIm;`} p0 = vec4f(aRe, aIm, zRe, zIm); zm10 = zm1Re; } else { let sq = sqrt(-nn) / opd; let r1 = zRe + sq; let r2 = zRe - sq; if (r1 != r2 && r1 != 0.0 && r2 != 0.0) { let det = 1.0 / r1 - 1.0 / r2; ${excite === 'strike' ? ` let a1 = b / det; let a2 = -a1;` : ` let a1 = (b - b / r2) / det; let a2 = b - a1;`} p0 = vec4f(a1, 0.0, r1, 0.0); p1 = vec4f(a2, 0.0, r2, 0.0); zm10 = zm1Re + sq; zm11 = zm1Re - sq; } }`; /** Event-row accessors for the tuning parameters a shader needs. */ const evAccess = (tuning: ParamSpec[], needed: string[]): string => needed .map((k) => { const j = tuning.findIndex((p) => p.key === k); if (j < 0) throw new Error(`event param '${k}' not in tuning list`); return ` let ${k} = ev[base + ${1 + j}u];`; }) .join('\n'); /** The per-axis excitation sums: for a strike, sum_i sin(m pi i/(n-1)) * * exp(-((i h - pos)/r)^2); a pluck replaces the Gaussian by the triangle. * One thread per (event, mode index). */ function excitationShader( object: ObjectKind, tuning: ParamSpec[], axis: 'x' | 'y', ): { code: string; constants: string[] } { const [n, m, L, hobj] = axis === 'x' ? ['n1', 'mx', 'L1', 'h1'] : ['n2', 'my', 'L2', 'h2']; let access: string; let term: string; if (object === 'string') { access = evAccess(tuning, ['pluckpos', 'amp']); term = ` let xi = f32(k) * ${hobj}; let ta = xi / (pluckpos * ${L}); let tb = (${L} - xi) / (${L} - pluckpos * ${L}); let u0 = amp * 0.5 * (ta + tb - abs(ta - tb));`; } else { const pos = object === 'bar' ? 'strikepos' : axis === 'x' ? 'strikex' : 'strikey'; access = evAccess(tuning, [pos, 'striker']); term = ` let xi = f32(k) * ${hobj}; let dx = (xi - ${pos} * ${L}) / striker; let u0 = exp(-dx * dx);`; } const consts: [string, 'u32' | 'f32'][] = [ [n, 'u32'], [m, 'u32'], ['ev_stride', 'u32'], [hobj, 'f32'], [L, 'f32'], ]; const code = `// Project this axis's excitation profile onto each recorded sine mode: // b${axis}[e * ${m} + (mode-1)]. One thread per (event, mode). const PI = 3.14159265358979; ${ov(consts)} ${buf([['bout', 'read_write'], ['__events', 'read']])} ${MAIN} if (i >= ${m} * arrayLength(&ev) / ev_stride) { return; } let mi = i % ${m}; let e = i / ${m}; let base = e * ev_stride; ${access} let mm = f32(mi + 1u); var s = 0.0; for (var k = 1u; k < ${n} - 1u; k++) { ${term} s += sin(mm * PI * f32(k) / f32(${n} - 1u)) * u0; } bout[i] = s; } `; return { code, constants: consts.map(([q]) => q) }; } /** The slots shader: every recorded mode of every event solved in closed * form and scaled into filter slots. */ function slotsShader(object: ObjectKind, tuning: ParamSpec[]): { code: string; constants: string[] } { if (object === 'membrane') { const consts: [string, 'u32' | 'f32'][] = [ ['mx', 'u32'], ['my', 'u32'], ['n_rows', 'u32'], ['n1', 'u32'], ['n2', 'u32'], ['ev_stride', 'u32'], ['h1', 'f32'], ['h2', 'f32'], ['L1', 'f32'], ['L2', 'f32'], ['dt', 'f32'], ]; const code = `// One thread per (event, recorded mode): the mode's excitation amplitude, // its recursion roots, and the acceleration-series slots that drive its // kernel row. const PI = 3.14159265358979; ${ov(consts)} ${buf([['slots', 'read_write'], ['__events', 'read'], ['bx', 'read'], ['by', 'read']])} ${SLOT_HELPERS} ${WRITE_SLOT} ${MAIN} let n_events = arrayLength(&ev) / ev_stride; if (i >= n_events * n_rows) { return; } let e = i / n_rows; let q = i % n_rows; let base = e * ev_stride; let start = ev[base]; ${evAccess(tuning, ['f0', 't60', 'sig1', 'vel', 'gair'])} let mi = q % mx; let ni = q / mx; let m = f32(mi + 1u); let n = f32(ni + 1u); let b = bx[e * mx + mi] * by[e * my + ni] * (-(dt * vel)) * 4.0 / (f32(n1 - 1u) * f32(n2 - 1u)); let s1 = sin(m * PI / (2.0 * f32(n1 - 1u))); let s2 = sin(n * PI / (2.0 * f32(n2 - 1u))); let lam = (4.0 / (h1 * h1)) * s1 * s1 + (4.0 / (h2 * h2)) * s2 * s2; let cm = 2.0 * f0 / sqrt(1.0 / (L1 * L1) + 1.0 / (L2 * L2)); let w2 = cm * cm * lam; ${rootsBlock('strike')} let ga = gair / (dt * dt); write_slot(i * 2u, acc_scale(p0, zm10, ga), f32(q), start); write_slot(i * 2u + 1u, acc_scale(p1, zm11, ga), f32(q), start); } `; return { code, constants: consts.map(([q]) => q) }; } if (object === 'bar') { const consts: [string, 'u32' | 'f32'][] = [ ['mx', 'u32'], ['n1', 'u32'], ['ev_stride', 'u32'], ['h1', 'f32'], ['L1', 'f32'], ['dt', 'f32'], ]; const code = `// One thread per (event, mode) of the bar: the tension-free stiff scheme, // partials at n^2 f1 up to the grid's own dispersion. const PI = 3.14159265358979; ${ov(consts)} ${buf([['slots', 'read_write'], ['__events', 'read'], ['bx', 'read']])} ${SLOT_HELPERS} ${WRITE_SLOT} ${MAIN} let n_events = arrayLength(&ev) / ev_stride; if (i >= n_events * mx) { return; } let e = i / mx; let mi = i % mx; let base = e * ev_stride; let start = ev[base]; ${evAccess(tuning, ['f1', 't60', 'sig1', 'vel', 'gair'])} let m = f32(mi + 1u); let b = bx[i] * (-(dt * vel)) * 2.0 / f32(n1 - 1u); let s1 = sin(m * PI / (2.0 * f32(n1 - 1u))); let lam = (4.0 / (h1 * h1)) * s1 * s1; let kap = 2.0 * f1 * (L1 * L1) / PI; let w2 = kap * kap * lam * lam; ${rootsBlock('strike')} let ga = gair / (dt * dt); write_slot(i * 2u, acc_scale(p0, zm10, ga), f32(mi), start); write_slot(i * 2u + 1u, acc_scale(p1, zm11, ga), f32(mi), start); } `; return { code, constants: consts.map(([q]) => q) }; } const consts: [string, 'u32' | 'f32'][] = [ ['mx', 'u32'], ['n1', 'u32'], ['ev_stride', 'u32'], ['h1', 'f32'], ['L1', 'f32'], ['dt', 'f32'], ]; const code = `// One thread per (event, mode) of the string: the stiff damped string's // roots, its radiation slots, and the bridge slots its arriving slope // drives (the last kernel row, row mx). const PI = 3.14159265358979; ${ov(consts)} ${buf([['slots', 'read_write'], ['__events', 'read'], ['bx', 'read']])} ${SLOT_HELPERS} ${WRITE_SLOT} ${MAIN} let n_events = arrayLength(&ev) / ev_stride; if (i >= n_events * mx) { return; } let e = i / mx; let mi = i % mx; let base = e * ev_stride; let start = ev[base]; ${evAccess(tuning, ['f0', 'B', 't60', 'sig1', 'gair', 'gbridge'])} let m = f32(mi + 1u); let b = bx[i] * 2.0 / f32(n1 - 1u); let s1 = sin(m * PI / (2.0 * f32(n1 - 1u))); let lam = (4.0 / (h1 * h1)) * s1 * s1; let cs = 2.0 * L1 * f0; let kap2 = B * (cs * L1 / PI) * (cs * L1 / PI); let w2 = cs * cs * lam + kap2 * lam * lam; ${rootsBlock('pluck')} let ga = gair / (dt * dt); write_slot(i * 4u, acc_scale(p0, zm10, ga), f32(mi), start); write_slot(i * 4u + 1u, acc_scale(p1, zm11, ga), f32(mi), start); // The bridge: the far end is pinned, so the arriving slope is // -u[n-2]/h1, and the generated step's term is // -(gbridge cs^2 / L) bprof fb; together, +gbridge cs^2/(L h1) times // sin(m pi (n1-2)/(n1-1)) times the next-displacement series. let gb = (gbridge * cs * cs / (L1 * h1)) * sin(m * PI * f32(n1 - 2u) / f32(n1 - 1u)); write_slot(i * 4u + 2u, next_scale(p0, gb), f32(mx), start); write_slot(i * 4u + 3u, next_scale(p1, gb), f32(mx), start); } `; return { code, constants: consts.map(([q]) => q) }; } const synthShader = { code: `// One thread per slot: a complex one-pole over this tile of samples, the // kernel row as input while it lasts and silence after (the closed-form // tail). State persists across tiles. ${ov([['tile', 'u32'], ['total_steps', 'u32'], ['ir_steps', 'u32'], ['n_slots', 'u32']])} ${buf([['slots', 'read'], ['state', 'read_write'], ['kernels', 'read'], ['partial', 'read_write'], ['__tick', 'read']])} ${MAIN} if (i >= n_slots) { return; } let cRe = slots[i * 6u]; let cIm = slots[i * 6u + 1u]; let zRe = slots[i * 6u + 2u]; let zIm = slots[i * 6u + 3u]; let row = u32(slots[i * 6u + 4u]); let start = u32(slots[i * 6u + 5u]); var wRe = state[i * 2u]; var wIm = state[i * 2u + 1u]; let t0 = tick[0] * tile; for (var s = 0u; s < tile; s++) { let t = t0 + s; if (t >= total_steps) { break; } var hv = 0.0; if (t >= start && t - start < ir_steps) { hv = kernels[row * ir_steps + (t - start)]; } let nRe = zRe * wRe - zIm * wIm + hv; let nIm = zRe * wIm + zIm * wRe; wRe = nRe; wIm = nIm; partial[i * tile + s] = cRe * wRe - cIm * wIm; } state[i * 2u] = wRe; state[i * 2u + 1u] = wIm; } `, constants: ['tile', 'total_steps', 'ir_steps', 'n_slots'], }; const reduceShader = { code: `// One thread per sample of the tile: sum every slot's contribution. ${ov([['tile', 'u32'], ['total_steps', 'u32'], ['n_slots', 'u32']])} ${buf([['partial', 'read'], ['trace', 'read_write'], ['__tick', 'read']])} ${MAIN} if (i >= tile) { return; } let t = tick[0] * tile + i; if (t >= total_steps) { return; } var s = 0.0; for (var q = 0u; q < n_slots; q++) { s += partial[q * tile + i]; } trace[t] = s; } `, constants: ['tile', 'total_steps', 'n_slots'], }; /* ----------------------------------------------------------- assembly --- */ export function buildInstrument(opts: BuildOptions): InstrumentFile { const { object, body } = opts; const op = objectParams[object]; const params: ParamSpec[] = [ { ...P('nx', 'air cells along x', '1', 96, 32, 256, 'accuracy', 'The air grid; ny = nz = nx/2. The simulation is honest to about c/(8h). Must be even.'), step: 2 }, P('Lx', 'domain length', 'm', 1.0, 0.6, 2.0, 'geometry', 'The air box around the instrument, twice as long in x as in y and z.'), P('ir_seconds', 'room response length', 's', 0.5, 0.1, 2.0, 'accuracy', 'How long the recorded room responses run; a room still ringing past this is truncated.'), ...(object === 'membrane' ? [P('mode_margin', 'mode set margin', '1', 2.0, 1.0, 4.0, 'accuracy', 'Head modes are recorded up to this multiple of what the air grid spatially resolves.')] : []), ...op.geometry, ...bodyParams[body], P('mic_x', 'microphone x', 'm', 0.15, -0.45, 0.45, 'geometry'), P('mic_y', 'microphone y', 'm', 0.1, -0.2, 0.2, 'geometry'), P('mic_z', 'microphone z', 'm', 0.12, -0.2, 0.2, 'geometry'), ...op.tuning, ]; for (const [k, v] of Object.entries(opts.defaults ?? {})) { const p = params.find((q) => q.key === k); if (!p) throw new Error(`default for unknown parameter '${k}'`); p.value = v; } const tuning = op.tuning; const evStride = 1 + tuning.length; const shaders: Record = { medium: mediumShader(object, body), sines_x: sinesShader('x'), impulse: impulseShader, step: stepShader, plane: planeShader('prof', 'plane'), xc: xcShader, proj: projShader(object), excite_x: excitationShader(object, tuning, 'x'), slots: slotsShader(object, tuning), synth: synthShader, reduce: reduceShader, }; if (object === 'membrane') { shaders.sines_y = sinesShader('y'); shaders.excite_y = excitationShader(object, tuning, 'y'); } if (object === 'string') { shaders.bplane = planeShader('bprof', 'bplane'); shaders.bridge_proj = bridgeProjShader; } const preBuffers = [ { name: 'wall', count: 'npts' }, { name: 'sig', count: 'npts' }, { name: 'prof', count: 'npts' }, ...(object === 'string' ? [{ name: 'bprof', count: 'npts' }] : []), { name: 'sx', count: 'mx * nx' }, ...(object === 'membrane' ? [{ name: 'sy', count: 'my * ny' }] : []), { name: 'p', count: 'npts' }, { name: 'pm', count: 'npts' }, { name: 'pn', count: 'npts' }, { name: 'plane', count: 'nx * ny' }, ...(object === 'string' ? [{ name: 'bplane', count: 'nx * ny' }] : []), { name: 'tmp', count: 'mx * ny' }, { name: 'kernels', count: 'n_rows * ir_steps' }, ]; const mediumOut = ['wall', 'sig', 'prof', ...(object === 'string' ? ['bprof'] : [])]; const preSetup: PassSpec[] = [ { shader: 'medium', buffers: mediumOut, dispatch: 'npts' }, { shader: 'sines_x', buffers: ['sx'], dispatch: 'mx * nx' }, ...(object === 'membrane' ? [{ shader: 'sines_y', buffers: ['sy'], dispatch: 'my * ny' }] : []), { shader: 'impulse', buffers: ['p', 'wall', 'sig'], dispatch: '1' }, ]; const preLoop: PassSpec[] = [ { shader: 'plane', buffers: ['plane', 'p', 'prof'], dispatch: 'nx * ny' }, { shader: 'xc', buffers: ['tmp', 'plane', 'sx'], dispatch: 'mx * ny' }, object === 'membrane' ? { shader: 'proj', buffers: ['kernels', 'tmp', 'sy', '__tick'], dispatch: 'n_rows' } : { shader: 'proj', buffers: ['kernels', 'tmp', '__tick'], dispatch: 'mx' }, ...(object === 'string' ? [ { shader: 'bplane', buffers: ['bplane', 'p', 'bprof'], dispatch: 'nx * ny' }, { shader: 'bridge_proj', buffers: ['kernels', 'bplane', '__tick'], dispatch: '1' }, ] : []), { shader: 'step', buffers: ['p', 'pm', 'pn', 'wall', 'sig'], dispatch: 'npts' }, ]; const slotsPer = object === 'string' ? 4 : 2; const slotsCount = object === 'membrane' ? 'n_events * n_rows * 2' : `n_events * mx * ${slotsPer}`; const renBuffers = [ { name: 'kernels', count: 'n_rows * ir_steps' }, { name: 'bx', count: 'n_events * mx' }, ...(object === 'membrane' ? [{ name: 'by', count: 'n_events * my' }] : []), { name: 'slots', count: 'n_slots * 6' }, { name: 'state', count: 'n_slots * 2' }, { name: 'partial', count: 'n_slots * tile' }, { name: 'trace', count: 'total_steps' }, ]; const MIC: [number, number, number] = [0.3, 0.85, 1.0]; const OUTLINE: [number, number, number] = [1.0, 0.85, 0.3]; const markers: MarkerSpec[] = [ object === 'membrane' ? { kind: 'rect', corner: ['-L1/2', '-L2/2', 'oz'], du: ['L1', '0', '0'], dv: ['0', 'L2', '0'], color: OUTLINE, label: 'head', } : { kind: 'line', from: ['-L1/2', 'oy', 'oz'], to: ['L1/2', 'oy', 'oz'], color: OUTLINE, label: object, }, { kind: 'point', at: ['mic_x', 'mic_y', 'mic_z'], color: MIC, label: 'mic' }, ]; const visual: VisualSpec = { dims: ['nx', 'ny', 'nz'], extents: ['Lx', 'Ly', 'Lz'], buffers: [ { name: 'wall', count: 'npts' }, { name: 'sig', count: 'npts' }, { name: 'prof', count: 'npts' }, ...(object === 'string' ? [{ name: 'bprof', count: 'npts' }] : []), ], setup: [{ shader: 'medium', buffers: mediumOut, dispatch: 'npts' }], volumes: [ // The body is the medium shader's own mask, so what the viewer shows // is what the solver sees, grid effects included. { buffer: 'wall', style: 'solid', invert: true, color: [0.62, 0.48, 0.34] }, { buffer: 'prof', style: 'cloud', color: [1.0, 0.45, 0.2] }, ...(object === 'string' ? [{ buffer: 'bprof', style: 'cloud' as const, color: [0.35, 0.65, 1.0] as [number, number, number] }] : []), ], markers, partials: { buffer: 'slots', stride: 6, zRe: 2, zIm: 3 }, }; const renSetup: PassSpec[] = [ { shader: 'excite_x', buffers: ['bx', '__events'], dispatch: 'n_events * mx' }, ...(object === 'membrane' ? [{ shader: 'excite_y', buffers: ['by', '__events'], dispatch: 'n_events * my' }] : []), object === 'membrane' ? { shader: 'slots', buffers: ['slots', '__events', 'bx', 'by'], dispatch: 'n_events * n_rows' } : { shader: 'slots', buffers: ['slots', '__events', 'bx'], dispatch: 'n_events * mx' }, ]; return { instrawave: 1, name: opts.name, description: opts.description, params, derived: derivedOf(object), shaders, precompute: { buffers: preBuffers, setup: preSetup, loop: { count: 'ir_steps', batch: 512, rotate: [['p', 'pm', 'pn']], body: preLoop }, outputs: ['kernels'], }, render: { buffers: renBuffers, inputs: ['kernels'], eventParams: tuning.map((p) => p.key), derived: { ev_stride: `${evStride}`, n_slots: slotsCount }, setup: renSetup, loop: { count: 'n_tiles', batch: 8, body: [ { shader: 'synth', buffers: ['slots', 'state', 'kernels', 'partial', '__tick'], dispatch: 'n_slots' }, { shader: 'reduce', buffers: ['partial', 'trace', '__tick'], dispatch: 'tile' }, ], }, trace: 'trace', }, visual, }; }