Form the Marguerre strains from the deviation, removing a noise floor on curved designs
Written as 0.5 (wx^2 - w0x^2) with the total slope, the membrane strain of a curved plate cancels two nearly equal numbers once the vibration is small compared with the profile, and in the single precision of the GPU the rounding error of the two squares is proportional to w0x^2, a constant. It acted as a broadband forcing that left the top octave of a full run at a fixed level, independent of the strike (the same 2.8 mPa for hits from 0.03 to 3 m/s on the 16 inch crash with a bell), so the high frequencies never decayed. The strains are now written so that every term carries the deviation, which is algebraically identical; the low bands of the run are unchanged to 0.1 dB, the top octave now decays with the rest, and the self-test against the float64 reference improves from 3e-6 to 6e-7 on the curved presets. A browser test runs a light hit on the curved crash and requires the top octave of the tail to lie more than 50 dB below the peak (it was 44 dB with the old expression, 66 dB now).
5 changed files+60−19
README.mdmodified+1−1View file
@@ -73,7 +73,7 @@ end
7373 - `tests/modal.test.ts`: the banded factorization and its inertia count, the tridiagonal eigensolver, all modes of a small square and a small disc against a dense solve (degenerate pairs included), the free square and disc frequencies against Gorman's and Leissa's values, and a plate on a stand.
7474 - `tests/synth.test.ts`: a single mode under the stick follows the closed-form damped sinusoid, the radiation gain of a piston equals the analytic Rayleigh sum, the attack and the tail join, the parts computed by the workers sum to the whole.
7575 - `tests/tutorial.test.ts`: every page of the tutorial renders with no design, with a plate, and with its modes, without undefined values and with balanced math delimiters; the actions run against an in-process host; the scaling pages' comparisons give the predicted ratios (2 for twice the thickness, 4 for half the size on a grid of the same cell count) to rounding.
76-- `tests/gpu.test.ts`: serves the app with Vite and drives a headless Chromium through Playwright: the in-page self-test of the WGSL kernel against the float64 reference on five presets (`?selftest=1` shows the same numbers in the browser), the linear model's modes and synthesis, the tutorial panel with its typeset mathematics and an action that computes the modes, an error in a script, short full runs on a free plate and on a stand, and a light hit simulated both ways, whose spectra must agree. It needs Playwright's Chromium (`npx playwright install chromium`) and a machine whose Chromium has WebGPU.
76+- `tests/gpu.test.ts`: serves the app with Vite and drives a headless Chromium through Playwright: the in-page self-test of the WGSL kernel against the float64 reference on five presets (`?selftest=1` shows the same numbers in the browser), the linear model's modes and synthesis, the tutorial panel with its typeset mathematics and an action that computes the modes, the loudness warning before playback, an error in a script, short full runs on a free plate and on a stand, the absence of a high-frequency noise floor in a curved design's full run, and a light hit simulated both ways, whose spectra must agree. It needs Playwright's Chromium (`npx playwright install chromium`) and a machine whose Chromium has WebGPU.
7777
7878 ## Layout
7979
THEORY.mdmodified+2−2View file
@@ -118,7 +118,7 @@ The full model keeps the nonlinear strains of Section 1 and steps the semi-discr
118118
119119 **Free plate.** A sheet hanging from a thread is treated as free, and the reaction to the stick is distributed over the plate as the inertial load of a rigid body, which is the inertia-relief equilibrium of a free body: the force of the previous step is applied as the rigid acceleration field $a_k = F\,(1/M + y_k y_c/I_x + x_k x_c/I_y)$ that balances the patch load and its moments about the center, so the plate gains no net momentum or angular momentum. The one-step lag has zero net impulse. A plate on a stand takes the reaction through the pads. This is the counterpart of dropping the rigid modes in the linear model.
120120
121-**Single precision.** With $w$ of a few millimeters and $\Delta t \approx 2\ \mu$s, a velocity below about $10^{-3}$ m/s produces an increment smaller than the resolution of $w$, and the rounding would act as a broadband noise source. The update of $w$ is therefore compensated: the rounding error of each addition is carried in a separate array and added back at the next step. The listener sum likewise forms the acceleration as the difference of successive delayed velocities at the audio rate, the average acceleration over one sample interval, rather than sampling the instantaneous acceleration, whose per-step rounding error of order $1\ \mathrm{m/s^2}$ per node would go straight into the audio.
121+**Single precision.** Three places in the kernel would otherwise turn rounding into audible noise. The first is the membrane strain of a curved design. Written directly, $\varepsilon_{xx} = u_x + \tfrac12[(w_0+w)_x^2 - w_{0,x}^2]$ subtracts two nearly equal numbers once $w$ is small compared with the profile, and the rounding error of the two squares is proportional to $w_{0,x}^2$, a constant: it does not fall as the vibration decays, and it acts as a broadband forcing that leaves a noise floor which never dies away. The strains are therefore written so that every term carries the deviation, $\varepsilon_{xx} = u_x + w_x(w_{0,x} + \tfrac12 w_x)$ and $\gamma_{xy} = u_y + v_x + w_{0,x}w_y + w_{0,y}w_x + w_x w_y$ with $w$ now the deviation alone, which is algebraically the same and has nothing large to cancel. On the 16 inch crash with a bell, whose largest profile slope is about 0.23, the direct form left the top octave at a constant level for the whole run, about 25 dB below the peak of a light hit; the rearranged form puts it more than 60 dB down and decaying. The second is the update of $w$. With $w$ of a few millimeters and $\Delta t \approx 2\ \mu$s, a velocity below about $10^{-3}$ m/s produces an increment smaller than the resolution of $w$, and the rounding would again act as a broadband noise source. The update of $w$ is therefore compensated: the rounding error of each addition is carried in a separate array and added back at the next step. The third is the listener sum, which forms the acceleration as the difference of successive delayed velocities at the audio rate, the average acceleration over one sample interval, rather than sampling the instantaneous acceleration, whose per-step rounding error of order $1\ \mathrm{m/s^2}$ per node would go straight into the audio. What remains after all three is a noise floor proportional to the signal, roughly 75 dB below the peak in the top octave for a strongly curved design and 100 dB for a flat one.
122122
123123 **The kernel.** The state is stored twice and each step reads one copy and writes the other. A step is one compute dispatch in which each workgroup of 256 threads advances a 12 by 12 tile of nodes: the threads load the tile with a two-node halo into shared memory, compute the membrane forces of the cells and the bending moments at the nodes once each, and then evaluate the update of the tile's nodes from shared memory. On grids of a few thousand nodes the cost of a step is the latency of one such workgroup rather than the arithmetic. A per-step uniform, addressed by a dynamic offset, carries the losses and the audio bookkeeping, so that a whole movie frame of steps is encoded in one pass; the retarded-time listener sum runs once per frame for all of the frame's samples from a ring buffer of the transverse velocities, and the following frame is encoded while the GPU works on the current one. A 16 inch cymbal at 8 mm runs at about 10 s per second of sound on the integrated GPU of a laptop.
124124
@@ -126,7 +126,7 @@ The full model keeps the nonlinear strains of Section 1 and steps the semi-discr
126126
127127 The tests check, without a GPU: that the linear stiffness of every configuration (uniform and tapered rectangles, a disc with a hole) has exactly the rigid modes in its null space; that the forces of the full model are the exact gradient of the discrete energy at a strongly nonlinear state, flat and bowed, uniform and tapered, free and on a stand; that the loss terms and the pad dashpots only remove energy; that the assembled $K$ equals the linearized force of the full model to $10^{-6}$ on a bowed, tapered, holed plate on a stand; that the banded factorization solves and that its inertia counts the eigenvalues; that every mode of a small square and of a small disc, degenerate pairs included, matches a dense eigensolution to $10^{-7}$ with no mode missed; that the free square and disc frequencies converge to Gorman's and Leissa's values as stated in Section 2; that the stand resists a translation with its total stiffness and a tilt with the pads' second moment, and that the lift-off rule unloads the right pad; that a single mode struck by the stick follows the closed-form damped sinusoid, that the radiation gain of a uniform piston equals the analytic Rayleigh sum, that the attack and the tail of the synthesized signal join, and that the parts computed by the workers sum to the whole.
128128
129-With a GPU, a Playwright test runs one step of the WGSL kernel from a smooth random state of about 2.5 thicknesses amplitude, with the losses and the stand with lift-off, on the flat, bowed, tapered, curved and elliptical presets, and compares the accelerations with the float64 reference: they agree to between $10^{-7}$ and $3\times10^{-6}$ relative, which is single-precision rounding. A short strike runs to completion on a free plate and on a stand. And a light hit (0.05 m/s) on the flat disc is simulated both ways: over the interval from 50 to 250 ms the magnitude spectra of the two listener signals below 6 kHz have a correlation of 0.985 and their levels agree within 13 percent. The remaining difference comes from the fitted two-parameter loss against the per-mode decay times, the pad dashpot against the loss factor, and the modal truncation at 16 kHz.
129+With a GPU, a Playwright test runs one step of the WGSL kernel from a smooth random state of about 2.5 thicknesses amplitude, with the losses and the stand with lift-off, on the flat, bowed, tapered, curved and elliptical presets, and compares the accelerations with the float64 reference: they agree to between $1.5\times10^{-7}$ and $6\times10^{-7}$ relative, which is single-precision rounding. A short strike runs to completion on a free plate and on a stand. A light hit on the curved crash is run to the end and the top octave of its tail is required to lie more than 50 dB below the peak, which catches a cancellation in the strains of the kind described in Section 5. And a light hit (0.05 m/s) on the flat disc is simulated both ways: over the interval from 50 to 250 ms the magnitude spectra of the two listener signals below 6 kHz have a correlation of 0.985 and their levels agree within 13 percent. The remaining difference comes from the fitted two-parameter loss against the per-mode decay times, the pad dashpot against the loss factor, and the modal truncation at 16 kHz.
130130
131131 ## 7. Limitations
132132
src/model/reference.tsmodified+13−10View file
@@ -125,16 +125,19 @@ export function accelerations(pl: Plate, f: Fields, out: Fields, liftoff = false
125125 }
126126 for (let gp = 0; gp < 4; gp++) {
127127 shapeDerivatives(gp, hx, hy, dNdx, dNdy);
128- let ux = 0, uy = 0, vx = 0, vy = 0, wx = 0, wy = 0, w0x = 0, w0y = 0;
128+ let ux = 0, uy = 0, vx = 0, vy = 0, dwx = 0, dwy = 0, w0x = 0, w0y = 0;
129129 for (let b = 0; b < 4; b++) {
130130 ux += dNdx[b] * U[b]; uy += dNdy[b] * U[b];
131131 vx += dNdx[b] * V[b]; vy += dNdy[b] * V[b];
132- wx += dNdx[b] * (Wn[b] + W0[b]); wy += dNdy[b] * (Wn[b] + W0[b]);
132+ dwx += dNdx[b] * Wn[b]; dwy += dNdy[b] * Wn[b];
133133 w0x += dNdx[b] * W0[b]; w0y += dNdy[b] * W0[b];
134134 }
135- const exx = ux + 0.5 * (wx * wx - w0x * w0x);
136- const eyy = vy + 0.5 * (wy * wy - w0y * w0y);
137- const gxy = uy + vx + wx * wy - w0x * w0y;
135+ const wx = w0x + dwx, wy = w0y + dwy;
136+ // See step.wgsl: every term carries the deviation, so that the single
137+ // precision of the GPU has nothing large to cancel.
138+ const exx = ux + dwx * (w0x + 0.5 * dwx);
139+ const eyy = vy + dwy * (w0y + 0.5 * dwy);
140+ const gxy = uy + vx + w0x * dwy + w0y * dwx + dwx * dwy;
138141 const Nxx = C * (exx + nu * eyy);
139142 const Nyy = C * (eyy + nu * exx);
140143 const Nxy = Gs * gxy;
@@ -255,17 +258,17 @@ export function potentialEnergy(pl: Plate, f: Fields, sr: StandRates | null = st
255258 const nn = [id(ci, cj), id(ci + 1, cj), id(ci + 1, cj + 1), id(ci, cj + 1)];
256259 for (let gp = 0; gp < 4; gp++) {
257260 shapeDerivatives(gp, hx, hy, dNdx, dNdy);
258- let ux = 0, uy = 0, vx = 0, vy = 0, wx = 0, wy = 0, w0x = 0, w0y = 0;
261+ let ux = 0, uy = 0, vx = 0, vy = 0, dwx = 0, dwy = 0, w0x = 0, w0y = 0;
259262 for (let b = 0; b < 4; b++) {
260263 const z0 = w0[nn[b]];
261264 ux += dNdx[b] * u[nn[b]]; uy += dNdy[b] * u[nn[b]];
262265 vx += dNdx[b] * v[nn[b]]; vy += dNdy[b] * v[nn[b]];
263- wx += dNdx[b] * (w[nn[b]] + z0); wy += dNdy[b] * (w[nn[b]] + z0);
266+ dwx += dNdx[b] * w[nn[b]]; dwy += dNdy[b] * w[nn[b]];
264267 w0x += dNdx[b] * z0; w0y += dNdy[b] * z0;
265268 }
266- const exx = ux + 0.5 * (wx * wx - w0x * w0x);
267- const eyy = vy + 0.5 * (wy * wy - w0y * w0y);
268- const gxy = uy + vx + wx * wy - w0x * w0y;
269+ const exx = ux + dwx * (w0x + 0.5 * dwx);
270+ const eyy = vy + dwy * (w0y + 0.5 * dwy);
271+ const gxy = uy + vx + w0x * dwy + w0y * dwx + dwx * dwy;
269272 Um += q * 0.5 * (C * (exx * exx + eyy * eyy + 2 * nu * exx * eyy) + Gs * gxy * gxy);
270273 }
271274 }
src/sim/shaders/step.wgslmodified+12−6View file
@@ -118,7 +118,7 @@ fn cellForces(ci: i32, cj: i32, lc: u32) {
118118 let Un = vec4f(U(ci, cj), U(ci + 1, cj), U(ci + 1, cj + 1), U(ci, cj + 1));
119119 let Vn = vec4f(V(ci, cj), V(ci + 1, cj), V(ci + 1, cj + 1), V(ci, cj + 1));
120120 let W0n = vec4f(W0(ci, cj), W0(ci + 1, cj), W0(ci + 1, cj + 1), W0(ci, cj + 1));
121- let Wn = vec4f(W(ci, cj), W(ci + 1, cj), W(ci + 1, cj + 1), W(ci, cj + 1)) + W0n; // total shape
121+ let Wdn = vec4f(W(ci, cj), W(ci + 1, cj), W(ci + 1, cj + 1), W(ci, cj + 1)); // deviation from the stress-free shape
122122 let q = 0.25 * P.hx * P.hy;
123123 let ihx = 1.0 / P.hx; let ihy = 1.0 / P.hy;
124124 let g3 = 0.57735026918962576;
@@ -130,12 +130,18 @@ fn cellForces(ci: i32, cj: i32, lc: u32) {
130130 let dNdy = vec4f(-(1.0 - xi), -(1.0 + xi), 1.0 + xi, 1.0 - xi) * (0.5 * ihy);
131131 let ux = dot(dNdx, Un); let uy = dot(dNdy, Un);
132132 let vx = dot(dNdx, Vn); let vy = dot(dNdy, Vn);
133- let wx = dot(dNdx, Wn); let wy = dot(dNdy, Wn);
133+ let dwx = dot(dNdx, Wdn); let dwy = dot(dNdy, Wdn);
134134 let w0x = dot(dNdx, W0n); let w0y = dot(dNdy, W0n);
135- // Marguerre strains: total slope minus the stress-free slope.
136- let exx = ux + 0.5 * (wx * wx - w0x * w0x);
137- let eyy = vy + 0.5 * (wy * wy - w0y * w0y);
138- let gxy = uy + vx + wx * wy - w0x * w0y;
135+ let wx = w0x + dwx; let wy = w0y + dwy; // total slope
136+ // Marguerre strains: the total slope minus the stress-free slope, written
137+ // so that every term carries the deviation. The direct form
138+ // 0.5 (wx^2 - w0x^2) subtracts two nearly equal numbers, and in single
139+ // precision its rounding error is proportional to w0x^2 and so does not
140+ // fall with the vibration: on a curved design it acts as a constant
141+ // broadband forcing that leaves a noise floor which never decays.
142+ let exx = ux + dwx * (w0x + 0.5 * dwx);
143+ let eyy = vy + dwy * (w0y + 0.5 * dwy);
144+ let gxy = uy + vx + w0x * dwy + w0y * dwx + dwx * dwy;
139145 let Nxx = C * (exx + P.nu * eyy);
140146 let Nyy = C * (eyy + P.nu * exx);
141147 let Nxy = G * gxy;
tests/gpu.test.tsmodified+32−0View file
@@ -253,6 +253,38 @@ test("the full simulation runs a short strike on a free plate and on a stand", a
253253 });
254254 });
255255
256+test("a curved design leaves no high-frequency noise floor in the full model", async () => {
257+ // The Marguerre strains must be formed so that every term carries the
258+ // deviation from the stress-free shape. The direct form 0.5 (wx^2 - w0x^2)
259+ // cancels two nearly equal numbers, and in the single precision of the GPU
260+ // its rounding error is proportional to w0x^2, so it does not fall with the
261+ // vibration: the top octave then sits at a constant level for the whole run
262+ // instead of decaying. A light hit is used, so that the true high-frequency
263+ // content is small and any floor dominates.
264+ await withPage("preset=crash-bell-16&mode=full&duration=0.5&spacing=0.016&strikeVelocity=0.05&fps=60", async (page, errors) => {
265+ if (!(await hasGpu(page))) { console.log("no WebGPU in this Chromium; skipping"); return; }
266+ await page.evaluate(() => (window as unknown as { cymbalsim: { runFull: () => Promise<void> } }).cymbalsim.runFull());
267+ const audio = await page.evaluate(() => { const r = (window as unknown as { cymbalsim: { app: { rec: { audio: Float32Array; samples: number } } } }).cymbalsim.app.rec; return Array.from(r.audio.subarray(0, r.samples)); });
268+ const fs = 48000, N = 4096;
269+ const peak = audio.reduce((m, v) => Math.max(m, Math.abs(v)), 0);
270+ assert.ok(peak > 1e-4, `peak ${peak} Pa`);
271+ // RMS pressure in the top octave over the last window of the run.
272+ const start = audio.length - N;
273+ const re = new Float64Array(N), im = new Float64Array(N);
274+ let wsum = 0;
275+ for (let i = 0; i < N; i++) { const w = 0.5 - 0.5 * Math.cos((2 * Math.PI * i) / N); re[i] = audio[start + i] * w; wsum += w * w; }
276+ fftInPlace(re, im);
277+ let e = 0;
278+ for (let k = Math.round((12000 / fs) * N); k < N / 2; k++) e += 2 * (re[k] * re[k] + im[k] * im[k]);
279+ const tailDb = 20 * Math.log10(Math.sqrt(e / (N * wsum)) / peak);
280+ console.log(`curved design: top octave of the tail ${tailDb.toFixed(1)} dB below the peak`);
281+ // With the cancellation the floor sits about 25 dB below the peak; without
282+ // it, below 60 dB. The bound is loose enough for any GPU's rounding.
283+ assert.ok(tailDb < -50, `top octave of the tail only ${tailDb.toFixed(1)} dB below the peak: a noise floor that does not decay`);
284+ assert.deepEqual(errors, []);
285+ });
286+});
287+
256288 test("a light hit sounds the same in the linear and the full model", async () => {
257289 await withPage("preset=flat-disc&duration=0.25&strikeVelocity=0.05&spacing=0.012&fps=60", async (page, errors) => {
258290 if (!(await hasGpu(page))) { console.log("no WebGPU in this Chromium; skipping"); return; }