2 * The computational grids: one rectangular box of air, and one small grid per
3 * vibrating object — a line of nodes for a string or bar, a rectangle of
4 * nodes for a membrane.
5 *
6 * The air grid is cell-centred and uniform, with the same spacing h along
7 * every axis; the box is twice as long in x as in y and z, so nx = 2*ny =
8 * 2*nz and a field is nx*ny*nz f32 values, flattened x-fastest: the point
9 * (ix, iy, iz) is element `ix + nx*(iy + ny*iz)`. Cell-centred so that no
10 * point sits exactly on the outer boundary: the stencil takes the field
11 * outside the domain to be zero, and the absorbing layer is meant to have
12 * swallowed the wave before it gets there.
13 *
14 * Object grids are node-centred, endpoints inclusive, so their edges land
15 * exactly on the terminations that pin them. Node coordinates are local —
16 * 0..L1 along the first axis (and 0..L2 along the second, for a rectangle) —
17 * because that is what the initial conditions (pluck shapes, strike bumps)
18 * are naturally written in. In the air's world coordinates every object is
19 * centred: a line runs from -L1/2 to +L1/2 along x, a rectangle spans
20 * (-L1/2..L1/2) x (-L2/2..L2/2) in the horizontal plane. Where an object
21 * sits in y and z is *not* baked in here: the coupling profiles the scene
22 * builds carry that, so moving an object is a scene re-evaluation rather
23 * than a recompile.
24 *
25 * The air's CFL condition fixes the timestep for everyone (see stableDt);
26 * each object then chooses the finest node spacing that its own scheme is
27 * stable at for anything the parameter sliders can reach, so moving a slider
28 * never forces a recompile — see makeLineGrid / makeRectGrid.
29 */
30import { CFL } from './units.ts';
32export interface AirGrid {
33 nx: number;
34 ny: number;
35 nz: number;
36 npts: number;
37 /** Domain extents in metres, centred on the origin. */
38 Lx: number;
39 Ly: number;
40 Lz: number;
41 /** Grid spacing, the same in all three directions. */
42 h: number;
43 /** Coordinates of every grid point, npts each, x fastest — as the shaders
44 * see them (f32) and as the scene .m is evaluated at (f64). */
45 x: Float32Array;
46 y: Float32Array;
47 z: Float32Array;
48 x64: Float64Array;
49 y64: Float64Array;
50 z64: Float64Array;
51}
53export function makeAirGrid(nx: number, Lx: number): AirGrid {
54 const ny = nx / 2;
55 const nz = nx / 2;
56 if (!Number.isInteger(ny)) throw new Error(`grid nx = ${nx} must be even`);
57 const h = Lx / nx;
58 const Ly = ny * h;
59 const Lz = nz * h;
60 const npts = nx * ny * nz;
61 const x64 = new Float64Array(npts);
62 const y64 = new Float64Array(npts);
63 const z64 = new Float64Array(npts);
64 for (let iz = 0; iz < nz; iz++) {
65 const zv = -Lz / 2 + (iz + 0.5) * h;
66 for (let iy = 0; iy < ny; iy++) {
67 const yv = -Ly / 2 + (iy + 0.5) * h;
68 const base = nx * (iy + ny * iz);
69 for (let ix = 0; ix < nx; ix++) {
70 const k = base + ix;
71 x64[k] = -Lx / 2 + (ix + 0.5) * h;
72 y64[k] = yv;
73 z64[k] = zv;
74 }
75 }
76 }
77 return {
78 nx, ny, nz, npts, Lx, Ly, Lz, h,
79 x: new Float32Array(x64),
80 y: new Float32Array(y64),
81 z: new Float32Array(z64),
82 x64, y64, z64,
83 };
84}
86/**
87 * The timestep the explicit leapfrog on the air grid is stable at, seconds.
88 *
89 * Leapfrog on p_tt = c^2 lap(p) is stable while dt^2 c^2 |lap|max <= 4, and
90 * the 7-point Laplacian's extreme eigenvalue is 12/h^2, so the condition is
91 * c*dt/h <= 1/sqrt(3). `cmax` is the fastest sound speed anywhere in the
92 * medium; since bodies are masks rather than fast materials (see
93 * src/mgpu/ops.ts), cmax is normally just the speed of air.
94 */
95export function stableDt(h: number, cmax: number, cfl = CFL): number {
96 return (cfl * h) / (Math.sqrt(3) * Math.max(cmax, 1e-12));
97}
99/** One vibrating object's grid: a line of nodes, or a rectangle of them. */
100export interface ObjectGrid {
101 kind: 'line' | 'rect';
102 /** Nodes along each axis; a line has n2 = 1. */
103 n1: number;
104 n2: number;
105 npts: number;
106 /** Extents in metres; a line has L2 = 0. */
107 L1: number;
108 L2: number;
109 /** Node spacing along each axis; a line has h2 = 0. */
110 h1: number;
111 h2: number;
112 /** Local node coordinates, npts each, first axis fastest: 0..L1 and
113 * 0..L2 inclusive. A line's ys is all zeros. */
114 xs: Float32Array;
115 ys: Float32Array;
116 xs64: Float64Array;
117 ys64: Float64Array;
118 /** 1 at interior nodes, 0 at the pinned boundary (the two ends of a line,
119 * the whole rim of a rectangle). The compiled step multiplies its update
120 * by this, which is what terminates the object. */
121 pin: Float32Array;
122}
124/**
125 * The coarsest node spacing the 1D stiff scheme demands at timestep dt.
126 *
127 * The scheme is the standard explicit leapfrog for
128 * u_tt = cs^2 u_xx - kap^2 u_xxxx - 2 sig0 u_t + 2 sig1 (u_xx)_t
129 * whose stability condition (Bilbao, Numerical Sound Synthesis, ch. 7) is
130 * hs^2 >= (a + sqrt(a^2 + 16 kap^2 dt^2)) / 2, a = cs^2 dt^2 + 4 sig1 dt.
131 * A bar is the cs = 0 case; an ideal string the kap = 0 case.
132 */
133export function stableH1(dt: number, cs: number, kap2: number, sig1: number): number {
134 const a = cs * cs * dt * dt + 4 * sig1 * dt;
135 return Math.sqrt((a + Math.sqrt(a * a + 16 * kap2 * dt * dt)) / 2);
136}
138/**
139 * The coarsest spacing the 2D membrane scheme demands at timestep dt,
140 * assuming (conservatively) the same spacing on both axes.
141 *
142 * For u_tt = c^2 lap(u) - 2 sig0 u_t + 2 sig1 (lap u)_t on the 5-point
143 * stencil, the extreme eigenvalue is 4/h1^2 + 4/h2^2, and stability needs
144 * (c^2 dt^2 + 4 sig1 dt)(1/h1^2 + 1/h2^2) <= 1. Requiring both spacings to
145 * be at least the h returned here (which satisfies the bound with 2/h^2) is
146 * therefore always safe.
147 */
148export function stableH2(dt: number, c: number, sig1: number): number {
149 return Math.sqrt(2 * (c * c * dt * dt + 4 * sig1 * dt));
150}
152/** Node-count caps. Objects are tiny beside the air grid; the caps exist so
e18e57eRender scenes outside the browserJeremy Magland 153 * a worst-case slider range cannot ask for an absurd grid. They are sized
154 * for the finest air grid the headless renderer is likely to use — a
155 * 256-node membrane is 65k points, still under one air slice. */
156const MAX_LINE_NODES = 1024;
157const MAX_RECT_NODES = 256;
159/**
160 * Build a line grid (string or bar) of length L1 at timestep dt, sized so
161 * the scheme is stable for every value the parameter sliders can reach.
162 *
163 * The air dictates dt, so an object has no say in the timestep; what it gets
164 * to choose is its own spacing, and the finest stable spacing is what
165 * maximizes its bandwidth. Sizing for the sliders' worst case (the stiffest,
166 * fastest, most damped settings) rather than their current values means
167 * moving a slider never forces a recompile: the grid stays valid, merely a
168 * little coarser than that setting alone would need.
169 */
170export function makeLineGrid(
171 L1: number,
172 dt: number,
173 worst: { cs: number; kap2: number; sig1: number },
174): ObjectGrid {
175 // 0.95: a little margin under the exact limit, as CFL is for the air.
176 const hMin = stableH1(dt, worst.cs, worst.kap2, worst.sig1) / 0.95;
177 const n1 = Math.max(5, Math.min(MAX_LINE_NODES, Math.floor(L1 / hMin) + 1));
178 const h1 = L1 / (n1 - 1);
179 const xs64 = new Float64Array(n1);
180 for (let i = 0; i < n1; i++) xs64[i] = i * h1;
181 const pin = new Float32Array(n1).fill(1);
182 pin[0] = 0;
183 pin[n1 - 1] = 0;
184 return {
185 kind: 'line',
186 n1, n2: 1, npts: n1,
187 L1, L2: 0, h1, h2: 0,
188 xs: new Float32Array(xs64), ys: new Float32Array(n1),
189 xs64, ys64: new Float64Array(n1),
190 pin,
191 };
192}
194/** Build a rectangular membrane grid of extent L1 x L2 at timestep dt, sized
195 * for the sliders' worst case like makeLineGrid. */
196export function makeRectGrid(
197 L1: number,
198 L2: number,
199 dt: number,
200 worst: { c: number; sig1: number },
201): ObjectGrid {
202 const hMin = stableH2(dt, worst.c, worst.sig1) / 0.95;
203 const size = (L: number): number =>
204 Math.max(5, Math.min(MAX_RECT_NODES, Math.floor(L / hMin) + 1));
205 const n1 = size(L1);
206 const n2 = size(L2);
207 const h1 = L1 / (n1 - 1);
208 const h2 = L2 / (n2 - 1);
209 const npts = n1 * n2;
210 const xs64 = new Float64Array(npts);
211 const ys64 = new Float64Array(npts);
212 const pin = new Float32Array(npts).fill(1);
213 for (let j = 0; j < n2; j++) {
214 for (let i = 0; i < n1; i++) {
215 const k = i + n1 * j;
216 xs64[k] = i * h1;
217 ys64[k] = j * h2;
218 if (i === 0 || i === n1 - 1 || j === 0 || j === n2 - 1) pin[k] = 0;
219 }
220 }
221 return {
222 kind: 'rect',
223 n1, n2, npts,
224 L1, L2, h1, h2,
225 xs: new Float32Array(xs64), ys: new Float32Array(ys64),
226 xs64, ys64,
227 pin,
228 };
229}