1/**
2 * The two computational grids: a rectangular box of air, and a line of
3 * string.
4 *
5 * The air grid is cell-centred and uniform, with the same spacing h along
6 * every axis; the box is twice as long in x (the string direction) as in y
7 * and z, so nx = 2*ny = 2*nz and a field is nx*ny*nz f32 values, flattened
8 * x-fastest: the point (ix, iy, iz) is element `ix + nx*(iy + ny*iz)`.
9 * Cell-centred so that no point sits exactly on the outer boundary: the
10 * stencil takes the field outside the domain to be zero, and the absorbing
11 * layer is meant to have swallowed the wave before it gets there.
12 *
13 * The string grid is node-centred: ns points from x = 0 to x = Ls inclusive,
14 * spacing hs = Ls/(ns-1), so its endpoints land exactly on the terminations
15 * that pin them. In the air's coordinates the string runs from -Ls/2 to
16 * +Ls/2 along x, at y = 0 and a height above the top plate that the scene
17 * decides.
18 */
19import { CFL } from './units.ts';
21export interface AirGrid {
22 nx: number;
23 ny: number;
24 nz: number;
25 npts: number;
26 /** Domain extents in metres, centred on the origin. */
27 Lx: number;
28 Ly: number;
29 Lz: number;
30 /** Grid spacing, the same in all three directions. */
31 h: number;
32 /** Coordinates of every grid point, npts each, x fastest — as the shaders
33 * see them (f32) and as the scene .m is evaluated at (f64). */
34 x: Float32Array;
35 y: Float32Array;
36 z: Float32Array;
37 x64: Float64Array;
38 y64: Float64Array;
39 z64: Float64Array;
40}
42export function makeAirGrid(nx: number, Lx: number): AirGrid {
43 const ny = nx / 2;
44 const nz = nx / 2;
45 if (!Number.isInteger(ny)) throw new Error(`grid nx = ${nx} must be even`);
46 const h = Lx / nx;
47 const Ly = ny * h;
48 const Lz = nz * h;
49 const npts = nx * ny * nz;
50 const x64 = new Float64Array(npts);
51 const y64 = new Float64Array(npts);
52 const z64 = new Float64Array(npts);
53 for (let iz = 0; iz < nz; iz++) {
54 const zv = -Lz / 2 + (iz + 0.5) * h;
55 for (let iy = 0; iy < ny; iy++) {
56 const yv = -Ly / 2 + (iy + 0.5) * h;
57 const base = nx * (iy + ny * iz);
58 for (let ix = 0; ix < nx; ix++) {
59 const k = base + ix;
60 x64[k] = -Lx / 2 + (ix + 0.5) * h;
61 y64[k] = yv;
62 z64[k] = zv;
63 }
64 }
65 }
66 return {
67 nx, ny, nz, npts, Lx, Ly, Lz, h,
68 x: new Float32Array(x64),
69 y: new Float32Array(y64),
70 z: new Float32Array(z64),
71 x64, y64, z64,
72 };
73}
75export interface StringGrid {
76 ns: number;
77 /** String length, metres. */
78 Ls: number;
79 /** Node spacing, Ls/(ns-1). */
80 hs: number;
81 /** Node positions, 0 to Ls inclusive. */
82 xs: Float32Array;
83 xs64: Float64Array;
84 /** 1 at interior nodes, 0 at the two pinned ends. The compiled step
85 * multiplies its update by this, which is what terminates the string. */
86 pin: Float32Array;
87}
89/**
90 * The timestep the explicit leapfrog on the air grid is stable at, seconds.
91 *
92 * Leapfrog on p_tt = c^2 lap(p) is stable while dt^2 c^2 |lap|max <= 4, and
93 * the 7-point Laplacian's extreme eigenvalue is 12/h^2, so the condition is
94 * c*dt/h <= 1/sqrt(3). `cmax` is the fastest sound speed anywhere in the
95 * medium; since the body's walls are a mask rather than a fast material (see
96 * src/mgpu/ops.ts), cmax is normally just the speed of air.
97 */
98export function stableDt(h: number, cmax: number, cfl = CFL): number {
99 return (cfl * h) / (Math.sqrt(3) * Math.max(cmax, 1e-12));
100}
102/**
103 * The coarsest string spacing the stiff-string scheme demands at timestep dt,
104 * in metres.
105 *
106 * The scheme (models/dulcimer.m) is the standard explicit leapfrog for
107 * u_tt = cs^2 u_xx - kap^2 u_xxxx - 2 sig0 u_t + 2 sig1 (u_xx)_t
108 * whose stability condition (Bilbao, Numerical Sound Synthesis, ch. 7) is
109 * hs^2 >= (a + sqrt(a^2 + 16 kap^2 dt^2)) / 2, a = cs^2 dt^2 + 4 sig1 dt.
110 */
111export function stableHs(dt: number, cs: number, kap2: number, sig1: number): number {
112 const a = cs * cs * dt * dt + 4 * sig1 * dt;
113 return Math.sqrt((a + Math.sqrt(a * a + 16 * kap2 * dt * dt)) / 2);
114}
116/**
117 * Build the string grid for length Ls at timestep dt, sized so that the
118 * scheme is stable for every value the parameter sliders can reach.
119 *
120 * The air's CFL condition fixes dt, so the string has no say in the timestep;
121 * what it gets to choose is its own spacing, and the finest stable spacing is
122 * what maximizes the string's bandwidth. Sizing for the sliders' worst case
123 * (highest fundamental, most stiffness, most frequency-dependent damping)
124 * rather than their current values means moving a slider never forces a
125 * recompile: the grid stays valid, merely a little coarser than that setting
126 * alone would need.
127 */
128export function makeStringGrid(
129 Ls: number,
130 dt: number,
131 worst: { f0: number; B: number; sig1: number },
132): StringGrid {
133 const cs = 2 * Ls * worst.f0;
134 const kap2 = (worst.B * cs * cs * Ls * Ls) / (Math.PI * Math.PI);
135 // 0.95: a little margin under the exact limit, as CFL is for the air.
136 const hsMin = stableHs(dt, cs, kap2, worst.sig1) / 0.95;
137 const ns = Math.max(8, Math.min(256, Math.floor(Ls / hsMin) + 1));
138 const hs = Ls / (ns - 1);
139 const xs64 = new Float64Array(ns);
140 for (let i = 0; i < ns; i++) xs64[i] = i * hs;
141 const pin = new Float32Array(ns).fill(1);
142 pin[0] = 0;
143 pin[ns - 1] = 0;
144 return { ns, Ls, hs, xs: new Float32Array(xs64), xs64, pin };
145}