Update simulation, step shader, renderer and theory notes - proof-of-concept/steelsim41

proof-of-concept / steelsim41

Update simulation, step shader, renderer and theory notes

Jeremy Magland <jmagland@flatironinstitute.org> committed 20 days ago commit a46d6445ca9a parent 6201bf2 Browse files

6 changed files+128−54

THEORY.mdmodified+1−1View file

@@ -54,7 +54,7 @@ A rigid rotation is not exactly a null mode of the von Karman energy: the strain

5454

5555 ## 3. Excitation

5656

57-**Drumstick strike.** The stick tip is a free mass $M$ (default 30 g, the effective mass of a stick in a wrist stroke rather than its full mass) arriving with speed $v_0$ (default 3 m/s) at a point of the plate, given as fractions of the bounding box. The contact patch consists of the plate nodes within the tip radius (5 mm) plus the nearest one. Its contact coordinate is the mass-weighted mean transverse displacement $w_c$ of these nodes, the penetration is $\delta = w_c - x$ with $x$ the tip position, and the force is the Hertz law $F = K\,\delta^{\alpha}$ for $\delta > 0$ and zero otherwise, with $\alpha = 1.5$ for a rounded tip and $K = 5\times10^{8}\ \mathrm{N/m^{1.5}}$ for wood on metal. The force is spread over the patch in proportion to the node masses, which carry the local thickness. Because that distribution is exactly conjugate to the mean displacement, the contact has a potential energy and the coupled system is Hamiltonian apart from the losses. On a free plate the force is balanced over the whole plate by the inertia-relief field $a_k = F\,(1/M_s + y_k y_c/I_x + x_k x_c/I_y)$, where $M_s$ is the mass of the plate, $I_x$ and $I_y$ its moments of inertia about the in-plane axes, and $(x_c, y_c)$ the mass-weighted center of the patch, so that the total force and the two moments of the load vanish (Section 2). The field uses the force of the previous step, a lag of one microsecond whose net impulse over the contact is zero. On a mounted plate no such field is applied. The stick obeys $M\ddot x = F$ and is stepped together with the plate; it may bounce and strike again if the plate comes back to meet it. The recording starts at the impact, with the plate at rest in its stress-free shape. Note that a thin plate is a soft target: its driving-point impedance is $8\sqrt{D\rho h}$, about 60 N s/m for 0.8 mm steel and about the same for 1 mm bronze, so the contact lasts well under a millisecond and the peak force is a few hundred newtons. The excitation is therefore broadband, and the sound is a short bright attack followed by the ring.

57+**Drumstick strike.** The stick tip is a free mass $M$ (default 30 g, the effective mass of a stick in a wrist stroke rather than its full mass) arriving with speed $v_0$ (default 3 m/s) at a point of the plate, given as fractions of the bounding box. The contact patch consists of the plate nodes within the tip radius (5 mm) plus the nearest one. Its contact coordinate is the mass-weighted mean transverse displacement $w_c$ of these nodes, the penetration is $\delta = w_c - x$ with $x$ the tip position, and the force is the Hertz law $F = K\,\delta^{\alpha}$ for $\delta > 0$ and zero otherwise, with $\alpha = 1.5$ for a rounded tip and $K = 5\times10^{8}\ \mathrm{N/m^{1.5}}$ for wood on metal. The force is spread over the patch in proportion to the node masses, which carry the local thickness. Because that distribution is exactly conjugate to the mean displacement, the contact has a potential energy and the coupled system is Hamiltonian apart from the losses. On a free plate the force is balanced over the whole plate by the inertia-relief field $a_k = F\,(1/M_s + y_k y_c/I_x + x_k x_c/I_y)$, where $M_s$ is the mass of the plate, $I_x$ and $I_y$ its moments of inertia about the in-plane axes, and $(x_c, y_c)$ the mass-weighted center of the patch, so that the total force and the two moments of the load vanish (Section 2). The field uses the force of the previous step, a lag of one microsecond whose net impulse over the contact is zero. On a mounted plate no such field is applied. The stick obeys $M\ddot x = F$ and is stepped together with the plate, so within the hit it may chatter, releasing and re-engaging as the plate runs ahead of it and comes back, which on a thin plate happens within a fraction of a millisecond. Once the contact has been released for 5 ms the stick is withdrawn and takes no further part. Without this rule a light plate on a soft mount, which the hit sets rocking as a whole, swings back after some tens of milliseconds and strikes the stick a second time; a drummer lifts the stick well before that. The recording starts at the impact, with the plate at rest in its stress-free shape. Note that a thin plate is a soft target: its driving-point impedance is $8\sqrt{D\rho h}$, about 60 N s/m for 0.8 mm steel and about the same for 1 mm bronze, so the contact lasts well under a millisecond and the peak force is a few hundred newtons. The excitation is therefore broadband, and the sound is a short bright attack followed by the ring.

5858

5959 ## 4. Spatial discretization

6060

src/main.tsmodified+4−3View file

@@ -232,7 +232,7 @@ function buildViewToolbar(): void {

232232 const hint = document.createElement("span");

233233 hint.className = "sub";

234234 hint.style.margin = "0";

235- hint.textContent = "3D: drag to orbit, wheel to zoom";

235+ hint.textContent = "3D: drag to orbit, wheel to zoom, double-click to reframe";

236236 tb.append(seg, labeled("Color", colorSel), labeled("Exaggeration ×", exagInput), labeled("Color range (mm, 0 = auto)", rangeInput), labeled("Show profile", shapeCheck), hint);

237237 }

238238

@@ -289,7 +289,8 @@ function buildMovieToolbar(): void {

289289 levelSel.appendChild(o);

290290 }

291291 gainSlider.type = "range";

292- gainSlider.min = "-40"; gainSlider.max = "60"; gainSlider.step = "1"; gainSlider.value = "20";

292+ // A hard strike heard from 0.3 m peaks near 40 Pa, so absolute level needs a large cut to stay under full scale.

293+ gainSlider.min = "-60"; gainSlider.max = "60"; gainSlider.step = "1"; gainSlider.value = "-34";

293294 gainSlider.style.width = "120px";

294295 gainSlider.style.flex = "none";

295296 const restart = () => {

@@ -622,7 +623,7 @@ function updateInfo(st: SimStatus): void {

622623 ["Max |w|", `${(st.wmax * 1e3).toFixed(3)} mm`],

623624 ["Audio", `${(st.samples / sim.fs).toFixed(2)} s, ${st.frames} frames`],

624625 ["Peak pressure", st.peakPressure > 0 ? `${st.peakPressure.toFixed(3)} Pa (${spl.toFixed(1)} dB SPL)` : "–"],

625- ["Stick", `${(st.stick.x * 1e3).toFixed(2)} mm, ${st.stick.v.toFixed(2)} m/s`],

626+ ["Stick", st.stick.withdrawn ? "withdrawn after the hit" : `${(st.stick.x * 1e3).toFixed(2)} mm, ${st.stick.v.toFixed(2)} m/s`],

626627 ["Contact force", `${st.stick.force.toFixed(1)} N (max ${st.stick.maxForce.toFixed(1)} N)`],

627628 ];

628629 let html = rows.map(([a, b]) => `<div class="row"><span>${a}</span><span>${b}</span></div>`).join("");

src/params.tsmodified+1−1View file

@@ -33,7 +33,7 @@ export const defaultParams: SimParams = {

3333 contactStiffness: 5e8,

3434 contactExponent: 1.5,

3535 contactRadius: 5e-3,

36- listenerDist: 1.0,

36+ listenerDist: 0.3,

3737 listenerAzimuth: 0,

3838 listenerElevation: 0,

3939 spacing: 0.008,

src/render/renderer.tsmodified+86−32View file

@@ -3,7 +3,7 @@

33 // or by the thickness. Both views paint the mount and the strike point on the

44 // surface as rings; the 3D view also stands a post at the mount, so that it

55 // stays visible when the ring is edge-on or lost on a steep bell, and stands a

6-// microphone on the line to the listener.

6+// microphone where the listener is, framing the view to reach both.

77

88 import { type Grid, cellMask, insideCell } from "../params.ts";

99 import { lookAt, multiply, perspective } from "./mat4.ts";

@@ -22,26 +22,42 @@ export interface Marker {

2222 const MOUNT_VERTS = 24 * 6 + 2 * 24 * 3;

2323 /** Vertices of the microphone: a 16 by 24 surface of revolution, plus a fan closing its back (see shaders.wgsl). */

2424 const MIC_VERTS = 16 * 24 * 6 + 16 * 3;

25-/** Its widest radius, as a fraction of the longer side of the plate's bounding box. */

26-const MIC_RADIUS_FRACTION = 0.038;

25+/** Its widest radius and the length of its body [m]: a large-diaphragm microphone, the same whatever the plate. */

26+const MIC_RADIUS = 0.028;

27+const MIC_ASPECT = 4.2; // must match shaders.wgsl

28+const MIC_LENGTH = MIC_RADIUS * MIC_ASPECT;

29+

30+/** Vertical field of view of the 3D view [rad], the framing margin around the scene, and the orbit limits. */

31+const FOV = 0.7;

32+const FIT_MARGIN = 1.1;

33+const DIST_MIN = 0.05;

34+const DIST_MAX = 200;

35+

36+type Vec3 = [number, number, number];

37+const dot = (a: Vec3, b: Vec3) => a[0] * b[0] + a[1] * b[1] + a[2] * b[2];

38+const cross = (a: Vec3, b: Vec3): Vec3 => [a[1] * b[2] - a[2] * b[1], a[2] * b[0] - a[0] * b[2], a[0] * b[1] - a[1] * b[0]];

39+const unit = (a: Vec3): Vec3 => { const n = Math.hypot(a[0], a[1], a[2]) || 1; return [a[0] / n, a[1] / n, a[2] / n]; };

2740

2841 /**

29- * How the microphone that stands for the listener is drawn, given the true

30- * listener distance and the longer side of the plate's bounding box.

31- *

32- * A listener a meter from a 16 inch cymbal is four times as far away as the

33- * cymbal is wide, so no framing can hold both at true scale; and the plate

34- * already fills most of the view, which leaves only a thin annulus around it

35- * for the microphone to stand in. The drawn distance therefore follows the

36- * true one through a saturating map that stays inside that annulus, and the

37- * microphone shrinks a little as the listener recedes, for the same reason a

38- * distant object looks small. The direction is drawn faithfully; the distance

39- * itself is the field in the Listener group, and it is that field, not this

40- * drawing, that sets the 1/r spreading and the spread of the delays.

42+ * The smallest orbit radius, about the plate's center along `dir`, at which

43+ * every one of `points` still falls inside the frustum, each point taken with

44+ * the radius it needs around it. The view is much wider than it is tall, so

45+ * fitting what is actually on screen rather than a bounding sphere makes a

46+ * large difference to how much of it the plate fills.

4147 */

42-export function drawnListener(dist: number, maxL: number): { at: number; radius: number } {

43- const f = dist / (dist + maxL); // 0 at the plate, approaching 1 far away

44- return { at: maxL * (0.26 + 0.36 * f), radius: MIC_RADIUS_FRACTION * maxL * (1.25 - 0.5 * f) };

48+function fitDistance(points: [Vec3, number][], dir: Vec3, aspect: number): number {

49+ const fwd: Vec3 = [-dir[0], -dir[1], -dir[2]];

50+ const right = unit(cross(fwd, [0, 1, 0]));

51+ const up = cross(right, fwd);

52+ const tanH = Math.tan(FOV / 2);

53+ let need = 0;

54+ for (const [p, r] of points) {

55+ const depth = dot(p, dir);

56+ const v = (Math.abs(dot(p, up)) + r) / tanH;

57+ const h = (Math.abs(dot(p, right)) + r) / (tanH * aspect);

58+ need = Math.max(need, depth + Math.max(v, h));

59+ }

60+ return Math.max(DIST_MIN, Math.min(DIST_MAX, need * FIT_MARGIN));

4561 }

4662

4763 export class Renderer {

@@ -55,6 +71,14 @@ export class Renderer {

5571 colorMode: ColorMode = "field";

5672 /** Orbit about the plate's center: yaw about the vertical, pitch above the plane of the plate. */

5773 camera = { yaw: 0.7, pitch: 0.5, dist: 0.95 };

74+ /**

75+ * Until the wheel is used, the orbit radius is chosen to frame the whole

76+ * scene, the microphone included. A listener a meter from a 16 inch cymbal

77+ * is drawn a meter from its center, so the cymbal is small in the view; that

78+ * is what a microphone a meter up really looks like. A double-click returns

79+ * to this framing, and to the opening angle, after moving the view.

80+ */

81+ private userZoomed = false;

5882 /** Markers drawn in both views; read at every render. */

5983 strike: Marker = { x: 0, y: 0, r: 0 };

6084 mount: Marker | null = null;

@@ -195,9 +219,34 @@ export class Renderer {

195219 this.canvas.addEventListener("wheel", (e) => {

196220 if (this.mode !== "3d") return;

197221 e.preventDefault();

198- this.camera.dist *= Math.exp(e.deltaY * 0.001);

199- this.camera.dist = Math.max(0.2, Math.min(5, this.camera.dist));

222+ this.userZoomed = true;

223+ this.camera.dist = Math.max(DIST_MIN, Math.min(DIST_MAX, this.camera.dist * Math.exp(e.deltaY * 0.001)));

200224 }, { passive: false });

225+ this.canvas.addEventListener("dblclick", () => {

226+ this.camera.yaw = 0.7;

227+ this.camera.pitch = 0.5;

228+ this.userZoomed = false;

229+ });

230+ }

231+

232+ /**

233+ * The points the 3D framing has to hold: the corners of the box around the

234+ * plate at its displayed height, and the two ends of the microphone, each

235+ * with the radius it needs around it.

236+ */

237+ private fitPoints(L: [number, number, number] | null, Ln: number, zmax: number): [Vec3, number][] {

238+ const g = this.grid!;

239+ const hx = 0.5 * (g.nx - 1) * g.hx, hy = 0.5 * (g.ny - 1) * g.hy, hz = Math.max(zmax * this.exaggeration, 1e-4);

240+ // The plate's frame in the world (see `world` in shaders.wgsl): x, height, -y.

241+ const pts: [Vec3, number][] = [];

242+ for (const sx of [-1, 1]) for (const sy of [-1, 1]) for (const sz of [-1, 1]) pts.push([[sx * hx, sz * hz, -sy * hy], 0]);

243+ if (L && Ln > 0) {

244+ for (const d of [Ln, Ln + MIC_LENGTH]) {

245+ const s = d / Ln;

246+ pts.push([[L[0] * s, L[2] * s, -L[1] * s], MIC_RADIUS]);

247+ }

248+ }

249+ return pts;

201250 }

202251

203252 private resize(): void {

@@ -221,12 +270,22 @@ export class Renderer {

221270 const range = this.colorRange > 0 ? this.colorRange : Math.max(auto, 1e-6);

222271 const u = new Float32Array(44);

223272 const ui = new Uint32Array(u.buffer);

224- const eye = [

225- this.camera.dist * Math.cos(this.camera.pitch) * Math.sin(this.camera.yaw),

226- this.camera.dist * Math.sin(this.camera.pitch),

227- this.camera.dist * Math.cos(this.camera.pitch) * Math.cos(this.camera.yaw),

273+ // The microphone stands where the listener is, at the listener distance from the center of the plate,

274+ // so the framing has to reach out to it as well as around the plate.

275+ const L = this.listener;

276+ const Ln = L ? Math.hypot(L[0], L[1], L[2]) : 0;

277+ const drawListener = !!L && Ln > 0;

278+ const dir: Vec3 = [

279+ Math.cos(this.camera.pitch) * Math.sin(this.camera.yaw),

280+ Math.sin(this.camera.pitch),

281+ Math.cos(this.camera.pitch) * Math.cos(this.camera.yaw),

228282 ];

229- const mvp = multiply(perspective(0.7, aspect, 0.02, 20), lookAt(eye, [0, 0, 0], [0, 1, 0]));

283+ const sceneRadius = this.fitPoints(drawListener && this.mode === "3d" ? L : null, Ln, auto);

284+ if (!this.userZoomed && this.mode === "3d") this.camera.dist = fitDistance(sceneRadius, dir, aspect);

285+ const eye: Vec3 = [this.camera.dist * dir[0], this.camera.dist * dir[1], this.camera.dist * dir[2]];

286+ const reach = Math.max(...sceneRadius.map(([p]) => Math.hypot(p[0], p[1], p[2])));

287+ const near = 0.02 * this.camera.dist, far = 4 * this.camera.dist + 4 * reach;

288+ const mvp = multiply(perspective(FOV, aspect, near, far), lookAt(eye, [0, 0, 0], [0, 1, 0]));

230289 u.set(mvp, 0);

231290 ui[16] = this.grid.nx; ui[17] = this.grid.ny;

232291 u[18] = this.grid.hx; u[19] = this.grid.hy;

@@ -239,14 +298,9 @@ export class Renderer {

239298 ui[34] = this.colorMode === "thickness" ? 1 : 0;

240299 ui[35] = this.mount ? 1 : 0;

241300 u[36] = this.hRange[0]; u[37] = this.hRange[1];

242- // Only the direction of the listener is drawn; its distance is far outside the framing of the plate.

243- const L = this.listener;

244- const Ln = L ? Math.hypot(L[0], L[1], L[2]) : 0;

245- const drawListener = !!L && Ln > 0;

246301 if (drawListener) { u[38] = L![0] / Ln; u[39] = L![1] / Ln; u[40] = L![2] / Ln; }

247- const mic = drawnListener(Ln, Math.max((this.grid.nx - 1) * this.grid.hx, (this.grid.ny - 1) * this.grid.hy));

248- u[41] = mic.at;

249- u[42] = mic.radius;

302+ u[41] = Ln;

303+ u[42] = MIC_RADIUS;

250304 ui[43] = drawListener ? 1 : 0;

251305 this.device.queue.writeBuffer(this.uniform, 0, u);

252306

src/sim/shaders/step.wgslmodified+29−12View file

@@ -31,7 +31,9 @@

3131 @group(0) @binding(4) var<storage, read_write> posOut: array<f32>;

3232 @group(0) @binding(5) var<storage, read_write> velOut: array<f32>;

3333 @group(0) @binding(6) var<storage, read_write> hist: array<f32>; // [slot * n + node]: transverse velocity at the audio sample instants

34-// Two 4-float slots of drumstick state (position, velocity, force, max force).

34+// Two 8-float slots of drumstick state: position, velocity, force, max force,

35+// time since the contact released, withdrawn flag, and two spare.

36+const STICK_RELEASE: f32 = 5e-3; // [s] a release longer than this ends the hit

3537 @group(0) @binding(7) var<storage, read_write> stick: array<f32>;

3638

3739 const TILE = 12; // nodes per workgroup along each axis (must match simulation.ts)

@@ -268,9 +270,14 @@ fn main(@builtin(local_invocation_id) lid: vec3u, @builtin(workgroup_id) wid: ve

268270 // plate nodes within the tip radius, plus the nearest one). The contact

269271 // coordinate is the mass-weighted mean displacement of the patch, and the

270272 // force is spread over the patch in proportion to the node masses, so the

271- // pair is conjugate and the contact has a potential energy.

273+ // pair is conjugate and the contact has a potential energy. The hit may

274+ // chatter (release and re-engage within a fraction of a millisecond) and

275+ // the plate may rebound onto the stick within a few milliseconds; both

276+ // count as the hit. Once the contact has been released for STICK_RELEASE

277+ // the stick is withdrawn, so a plate that swings back later as a whole

278+ // (a light plate rocking on its mount) cannot be struck a second time.

272279 if ((SP.flags & 1u) != 0u) {

273- let sIn = 4u * SP.parity;

280+ let sIn = 8u * SP.parity;

274281 if (P.mountOn == 0u) {

275282 // Free plate: inertia relief. The force of the previous step is applied

276283 // to the whole plate as the rigid acceleration field that balances the

@@ -282,24 +289,34 @@ fn main(@builtin(local_invocation_id) lid: vec3u, @builtin(workgroup_id) wid: ve

282289 }

283290 if (inContactPatch(i, j)) {

284291 let sx = stick[sIn]; let sv = stick[sIn + 1u]; let sMax = stick[sIn + 3u];

285- var wc = 0.0;

286- for (var pj = i32(P.patchJ0); pj <= i32(P.patchJ1); pj++) {

287- for (var pi = i32(P.patchI0); pi <= i32(P.patchI1); pi++) {

288- if (inContactPatch(pi, pj)) { wc += nodeMass(pi, pj) * posIn[2u * n + id(pi, pj)]; }

292+ var released = stick[sIn + 4u]; // time since the contact last released [s]

293+ var withdrawn = stick[sIn + 5u];

294+ var Fc = 0.0;

295+ if (withdrawn == 0.0) {

296+ var wc = 0.0;

297+ for (var pj = i32(P.patchJ0); pj <= i32(P.patchJ1); pj++) {

298+ for (var pi = i32(P.patchI0); pi <= i32(P.patchI1); pi++) {

299+ if (inContactPatch(pi, pj)) { wc += nodeMass(pi, pj) * posIn[2u * n + id(pi, pj)]; }

300+ }

301+ }

302+ wc /= P.patchMass;

303+ let pen = wc - sx; // the stick sits on the +z side and moves toward -z

304+ if (pen > 0.0) { Fc = P.K * pow(pen, P.alpha); released = 0.0; }

305+ else if (sMax > 0.0) {

306+ released += P.dt;

307+ if (released > STICK_RELEASE) { withdrawn = 1.0; }

289308 }

290309 }

291- wc /= P.patchMass;

292- let pen = wc - sx; // the stick sits on the +z side and moves toward -z

293- var Fc = 0.0;

294- if (pen > 0.0) { Fc = P.K * pow(pen, P.alpha); }

295310 acc.z -= Fc / P.patchMass;

296311 if (k == P.strikeNode) {

297- let sOut = 4u * (1u - SP.parity);

312+ let sOut = 8u * (1u - SP.parity);

298313 let svn = sv + P.dt * Fc / P.stickMass;

299314 stick[sOut] = sx + P.dt * svn;

300315 stick[sOut + 1u] = svn;

301316 stick[sOut + 2u] = Fc;

302317 stick[sOut + 3u] = max(sMax, Fc);

318+ stick[sOut + 4u] = released;

319+ stick[sOut + 5u] = withdrawn;

303320 }

304321 }

305322 }

src/sim/simulation.tsmodified+7−5View file

@@ -44,6 +44,8 @@ export interface StickState {

4444 v: number;

4545 force: number;

4646 maxForce: number;

47+ /** True once the contact has been released for longer than the shader's STICK_RELEASE; the stick then no longer touches the plate. */

48+ withdrawn: boolean;

4749 }

4850

4951 export interface SimStatus {

@@ -63,7 +65,7 @@ export interface SimStatus {

6365 const PARAMS_LEN = 40; // floats in the Params uniform

6466 const STEP_LEN = 12; // floats in the StepParams uniform

6567 const STATS_LEN = 16;

66-const STICK_LEN = 8;

68+const STICK_LEN = 16; // two 8-float slots of stick state

6769 /** Nodes per workgroup along each axis in step.wgsl. */

6870 const TILE = 12;

6971

@@ -206,9 +208,9 @@ export class Simulation {

206208 const cp = contactPatch(params, plate);

207209 const mp = mountPatch(plate);

208210 const stick = new Float32Array(STICK_LEN);

209- stick[1] = -params.strikeVelocity; stick[5] = -params.strikeVelocity;

211+ stick[1] = -params.strikeVelocity; stick[9] = -params.strikeVelocity;

210212 device.queue.writeBuffer(this.stickBuf, 0, stick);

211- this.stick = { x: 0, v: -params.strikeVelocity, force: 0, maxForce: 0 };

213+ this.stick = { x: 0, v: -params.strikeVelocity, force: 0, maxForce: 0, withdrawn: false };

212214

213215 const [sx, sy] = strikePoint(params, plate);

214216 const pf = new Float32Array(PARAMS_LEN);

@@ -563,8 +565,8 @@ export class Simulation {

563565 this.ke = stats[12];

564566 this.wmax = stats[13];

565567 // The slot written by the last step is the one the next step reads.

566- const base = this.samplesPerBatch + n + STATS_LEN + 4 * b.stickParity;

567- this.stick = { x: view[base], v: view[base + 1], force: view[base + 2], maxForce: view[base + 3] };

568+ const base = this.samplesPerBatch + n + STATS_LEN + 8 * b.stickParity;

569+ this.stick = { x: view[base], v: view[base + 1], force: view[base + 2], maxForce: view[base + 3], withdrawn: view[base + 5] !== 0 };

568570 this.appendAudio(view.slice(0, b.nAudio));

569571 this.frames.push(view.slice(this.samplesPerBatch, this.samplesPerBatch + n));

570572 b.staging.unmap();