f8a8a6cNonlinear steel sheet gong simulator with WebGPUJeremy Magland 1// One explicit time step (symplectic Euler) of the von Karman plate.
2//
3// State layout (structure of arrays, n = nx * ny nodes):
4// pos: [u | v | w | wLo] wLo carries the rounding error of the w update
5// vel: [du | dv | dw]
c63c76aAdd elliptical sheet shape on a masked gridJeremy Magland 6// Each 16x16 workgroup loads its tile plus a 2-node halo of w and of the node
7// weights into shared memory once, so the stencils below never touch global
8// memory. Nodes outside the sheet (weight 0) are skipped: a second difference
9// exists only where the node and both neighbours belong to the sheet, and a
10// cell contributes only when its four corners do. That reproduces the free
11// edge of the rectangle and gives a staircase free edge for an ellipse.
f8a8a6cNonlinear steel sheet gong simulator with WebGPUJeremy Magland 12
13@group(0) @binding(0) var<uniform> P: Params;
14@group(0) @binding(1) var<uniform> SP: StepParams;
15@group(0) @binding(2) var<storage, read> posIn: array<f32>;
16@group(0) @binding(3) var<storage, read> velIn: array<f32>;
17@group(0) @binding(4) var<storage, read_write> posOut: array<f32>;
18@group(0) @binding(5) var<storage, read_write> velOut: array<f32>;
bfde5e4High-frequency ringing: in-plane loss, band-limited listener sum, damping defaultsJeremy Magland 19@group(0) @binding(6) var<storage, read_write> hist: array<f32>; // [slot * n + node] velocity of w at the audio sample instants
10998acDrumstick strike excitation (Hertz contact with a free mass)Jeremy Magland 20// Constraint reaction (acceleration of w) per node, followed by two 4-float
21// slots of drumstick state (position, velocity, force, max force).
22@group(0) @binding(7) var<storage, read_write> react: array<f32>;
8dcc11fInitial-shape imperfection (Marguerre shallow-shell strains)Jeremy Magland 23@group(0) @binding(8) var<storage, read> w0: array<f32>; // stress-free initial shape
f8a8a6cNonlinear steel sheet gong simulator with WebGPUJeremy Magland 24
25const TILE = 16;
26const HALO = 2;
27const TW = 20; // TILE + 2 * HALO
28var<workgroup> tw: array<f32, 400>;
c63c76aAdd elliptical sheet shape on a masked gridJeremy Magland 29// Node weights of the same tile in quarters (0..4), packed four to a word to stay
30// within the default 16 KB of workgroup storage.
31var<workgroup> tn: array<u32, 100>;
f8a8a6cNonlinear steel sheet gong simulator with WebGPUJeremy Magland 32// Membrane force contributions of the 17x17 cells around the tile, 4 nodes x 3 components each.
33const CT = 17;
34var<workgroup> cellF: array<f32, 3468>;
36// Tile origin in global node coordinates (set once per invocation).
37var<private> ox: i32;
38var<private> oy: i32;
39var<private> ihx2: f32;
40var<private> ihy2: f32;
41var<private> ihxy: f32;
43fn id(i: i32, j: i32) -> u32 { return u32(j) * P.nx + u32(i); }
44fn lidx(i: i32, j: i32) -> u32 { return u32(j - oy) * u32(TW) + u32(i - ox); }
45// Global indices; the tile covers [ox, ox + TW) x [oy, oy + TW).
46fn W(i: i32, j: i32) -> f32 { return tw[lidx(i, j)]; }
c63c76aAdd elliptical sheet shape on a masked gridJeremy Magland 47fn NWq(i: i32, j: i32) -> u32 { let l = lidx(i, j); return (tn[l >> 2u] >> ((l & 3u) * 8u)) & 0xFFu; }
48fn NW(i: i32, j: i32) -> f32 { return 0.25 * f32(NWq(i, j)); }
49fn INS(i: i32, j: i32) -> bool { return NWq(i, j) != 0u; }
50fn insCell(ci: i32, cj: i32) -> bool { return INS(ci, cj) && INS(ci + 1, cj) && INS(ci + 1, cj + 1) && INS(ci, cj + 1); }
f8a8a6cNonlinear steel sheet gong simulator with WebGPUJeremy Magland 51fn U(i: i32, j: i32) -> f32 { return posIn[id(i, j)]; }
52fn V(i: i32, j: i32) -> f32 { return posIn[P.n + id(i, j)]; }
53fn DW(i: i32, j: i32) -> f32 { return velIn[2u * P.n + id(i, j)]; }
8dcc11fInitial-shape imperfection (Marguerre shallow-shell strains)Jeremy Magland 54fn W0(i: i32, j: i32) -> f32 { return w0[id(i, j)]; }
f8a8a6cNonlinear steel sheet gong simulator with WebGPUJeremy Magland 55
ebf81f1Add an optional k^4 loss term (sigma2), off by defaultJeremy Magland 56// Graph Laplacian of the transverse velocity at node (i, j), over neighbours in the sheet.
57fn lapDW(i: i32, j: i32) -> f32 {
58 let c = DW(i, j);
59 var lap = 0.0;
60 if (INS(i - 1, j)) { lap += (DW(i - 1, j) - c) * ihx2; }
61 if (INS(i + 1, j)) { lap += (DW(i + 1, j) - c) * ihx2; }
62 if (INS(i, j - 1)) { lap += (DW(i, j - 1) - c) * ihy2; }
63 if (INS(i, j + 1)) { lap += (DW(i, j + 1) - c) * ihy2; }
64 return lap;
65}
f8a8a6cNonlinear steel sheet gong simulator with WebGPUJeremy Magland 67// Weighted bending "moments" at nodes; zero where the second difference is undefined.
68fn Ma(i: i32, j: i32) -> f32 {
c63c76aAdd elliptical sheet shape on a masked gridJeremy Magland 69 if (!(INS(i - 1, j) && INS(i, j) && INS(i + 1, j))) { return 0.0; }
f8a8a6cNonlinear steel sheet gong simulator with WebGPUJeremy Magland 70 let a = (W(i + 1, j) - 2.0 * W(i, j) + W(i - 1, j)) * ihx2;
71 var b = 0.0;
c63c76aAdd elliptical sheet shape on a masked gridJeremy Magland 72 if (INS(i, j - 1) && INS(i, j + 1)) { b = (W(i, j + 1) - 2.0 * W(i, j) + W(i, j - 1)) * ihy2; }
73 return NW(i, j) * (a + P.nu * b);
75fn Mb(i: i32, j: i32) -> f32 {
c63c76aAdd elliptical sheet shape on a masked gridJeremy Magland 76 if (!(INS(i, j - 1) && INS(i, j) && INS(i, j + 1))) { return 0.0; }
f8a8a6cNonlinear steel sheet gong simulator with WebGPUJeremy Magland 77 let b = (W(i, j + 1) - 2.0 * W(i, j) + W(i, j - 1)) * ihy2;
78 var a = 0.0;
c63c76aAdd elliptical sheet shape on a masked gridJeremy Magland 79 if (INS(i - 1, j) && INS(i + 1, j)) { a = (W(i + 1, j) - 2.0 * W(i, j) + W(i - 1, j)) * ihx2; }
80 return NW(i, j) * (b + P.nu * a);
c63c76aAdd elliptical sheet shape on a masked gridJeremy Magland 82// Mixed second difference at a cell center; zero outside the sheet.
f8a8a6cNonlinear steel sheet gong simulator with WebGPUJeremy Magland 83fn Cc(ci: i32, cj: i32) -> f32 {
c63c76aAdd elliptical sheet shape on a masked gridJeremy Magland 84 if (!insCell(ci, cj)) { return 0.0; }
f8a8a6cNonlinear steel sheet gong simulator with WebGPUJeremy Magland 85 return (W(ci + 1, cj + 1) - W(ci + 1, cj) - W(ci, cj + 1) + W(ci, cj)) * ihxy;
86}
88// Membrane forces of the bilinear cell (ci, cj) on its four nodes (2x2 Gauss),
89// written to cellF at local cell index lc.
90fn cellForces(ci: i32, cj: i32, lc: u32) {
91 let base = lc * 12u;
f8a8a6cNonlinear steel sheet gong simulator with WebGPUJeremy Magland 93 for (var q = 0u; q < 12u; q++) { cellF[base + q] = 0.0; }
94 return;
95 }
96 let Un = vec4f(U(ci, cj), U(ci + 1, cj), U(ci + 1, cj + 1), U(ci, cj + 1));
97 let Vn = vec4f(V(ci, cj), V(ci + 1, cj), V(ci + 1, cj + 1), V(ci, cj + 1));
8dcc11fInitial-shape imperfection (Marguerre shallow-shell strains)Jeremy Magland 98 let W0n = vec4f(W0(ci, cj), W0(ci + 1, cj), W0(ci + 1, cj + 1), W0(ci, cj + 1));
99 let Wn = vec4f(W(ci, cj), W(ci + 1, cj), W(ci + 1, cj + 1), W(ci, cj + 1)) + W0n; // total shape
f8a8a6cNonlinear steel sheet gong simulator with WebGPUJeremy Magland 100 let q = 0.25 * P.hx * P.hy;
101 let ihx = 1.0 / P.hx; let ihy = 1.0 / P.hy;
102 let g3 = 0.57735026918962576;
103 var Fu = vec4f(0.0); var Fv = vec4f(0.0); var Fw = vec4f(0.0);
104 for (var gp = 0; gp < 4; gp++) {
105 var xi = -g3; if ((gp & 1) != 0) { xi = g3; }
106 var eta = -g3; if ((gp & 2) != 0) { eta = g3; }
107 let dNdx = vec4f(-(1.0 - eta), 1.0 - eta, 1.0 + eta, -(1.0 + eta)) * (0.5 * ihx);
108 let dNdy = vec4f(-(1.0 - xi), -(1.0 + xi), 1.0 + xi, 1.0 - xi) * (0.5 * ihy);
109 let ux = dot(dNdx, Un); let uy = dot(dNdy, Un);
110 let vx = dot(dNdx, Vn); let vy = dot(dNdy, Vn);
111 let wx = dot(dNdx, Wn); let wy = dot(dNdy, Wn);
8dcc11fInitial-shape imperfection (Marguerre shallow-shell strains)Jeremy Magland 112 let w0x = dot(dNdx, W0n); let w0y = dot(dNdy, W0n);
113 // Marguerre strains: total slope minus the stress-free slope.
114 let exx = ux + 0.5 * (wx * wx - w0x * w0x);
115 let eyy = vy + 0.5 * (wy * wy - w0y * w0y);
116 let gxy = uy + vx + wx * wy - w0x * w0y;
f8a8a6cNonlinear steel sheet gong simulator with WebGPUJeremy Magland 117 let Nxx = P.C * (exx + P.nu * eyy);
118 let Nyy = P.C * (eyy + P.nu * exx);
119 let Nxy = P.G * gxy;
120 Fu -= q * (Nxx * dNdx + Nxy * dNdy);
121 Fv -= q * (Nyy * dNdy + Nxy * dNdx);
122 Fw -= q * ((Nxx * wx + Nxy * wy) * dNdx + (Nxy * wx + Nyy * wy) * dNdy);
123 }
124 for (var a = 0u; a < 4u; a++) {
125 cellF[base + a * 3u] = Fu[a];
126 cellF[base + a * 3u + 1u] = Fv[a];
127 cellF[base + a * 3u + 2u] = Fw[a];
128 }
129}
131// Stored membrane force on local node `a` of cell (ci, cj).
132fn cellForce(ci: i32, cj: i32, a: u32) -> vec3f {
133 let lc = u32(cj - oy - 1) * u32(CT) + u32(ci - ox - 1);
134 let base = lc * 12u + a * 3u;
135 return vec3f(cellF[base], cellF[base + 1u], cellF[base + 2u]);
136}
138@compute @workgroup_size(16, 16)
139fn main(@builtin(global_invocation_id) gid: vec3u, @builtin(local_invocation_id) lid: vec3u, @builtin(workgroup_id) wid: vec3u) {
140 let n = P.n;
141 ox = i32(wid.x) * TILE - HALO;
142 oy = i32(wid.y) * TILE - HALO;
143 ihx2 = 1.0 / (P.hx * P.hx); ihy2 = 1.0 / (P.hy * P.hy); ihxy = 1.0 / (P.hx * P.hy);
c63c76aAdd elliptical sheet shape on a masked gridJeremy Magland 144 // Cooperative load of the w and weight tiles (400 entries each, 256 threads).
f8a8a6cNonlinear steel sheet gong simulator with WebGPUJeremy Magland 145 let lt = lid.y * 16u + lid.x;
146 for (var e = lt; e < 400u; e += 256u) {
147 let gi = ox + i32(e % u32(TW)); let gj = oy + i32(e / u32(TW));
148 var ww = 0.0;
149 if (gi >= 0 && gi < i32(P.nx) && gj >= 0 && gj < i32(P.ny)) { ww = posIn[2u * n + id(gi, gj)]; }
150 tw[e] = ww;
151 }
153 var word = 0u;
154 for (var q = 0u; q < 4u; q++) {
155 let e = lt * 4u + q;
156 let gi = ox + i32(e % u32(TW)); let gj = oy + i32(e / u32(TW));
157 var wq = 0u;
158 if (gi >= 0 && gi < i32(P.nx) && gj >= 0 && gj < i32(P.ny)) { wq = u32(round(4.0 * nw[id(gi, gj)])); }
159 word |= wq << (q * 8u);
160 }
161 tn[lt] = word;
162 }
164 // Membrane forces of the 17x17 cells adjacent to the tile's nodes.
165 for (var e = lt; e < 289u; e += 256u) {
166 let lci = i32(e % u32(CT)); let lcj = i32(e / u32(CT));
167 cellForces(ox + 1 + lci, oy + 1 + lcj, e);
168 }
169 workgroupBarrier();
171 let i = i32(gid.x); let j = i32(gid.y);
172 if (i >= i32(P.nx) || j >= i32(P.ny)) { return; }
173 let k = id(i, j);
175 if (weight == 0.0) {
176 // Outside the sheet: the node stays at rest.
177 posOut[k] = 0.0; posOut[n + k] = 0.0; posOut[2u * n + k] = 0.0; posOut[3u * n + k] = 0.0;
178 velOut[k] = 0.0; velOut[n + k] = 0.0; velOut[2u * n + k] = 0.0;
179 if ((SP.flags & 2u) != 0u) { react[k] = 0.0; }
180 return;
181 }
f8a8a6cNonlinear steel sheet gong simulator with WebGPUJeremy Magland 182
183 let bend = (Ma(i - 1, j) - 2.0 * Ma(i, j) + Ma(i + 1, j)) * ihx2
184 + (Mb(i, j - 1) - 2.0 * Mb(i, j) + Mb(i, j + 1)) * ihy2
185 + 2.0 * (1.0 - P.nu) * (Cc(i, j) - Cc(i - 1, j) - Cc(i, j - 1) + Cc(i - 1, j - 1)) * ihxy;
186 var F = vec3f(0.0, 0.0, -P.D * P.hx * P.hy * bend);
187 F += cellForce(i, j, 0u);
188 F += cellForce(i - 1, j, 1u);
189 F += cellForce(i - 1, j - 1, 2u);
190 F += cellForce(i, j - 1, 3u);
c63c76aAdd elliptical sheet shape on a masked gridJeremy Magland 191 let mass = P.rhoh * P.hx * P.hy * weight;
f8a8a6cNonlinear steel sheet gong simulator with WebGPUJeremy Magland 192 var acc = F * (1.0 / mass);
ebf81f1Add an optional k^4 loss term (sigma2), off by defaultJeremy Magland 194 // Loss: frequency-independent on all components; on w also the graph Laplacian of
195 // the velocity (sigma1, ~k^2) and, when enabled, its biharmonic (sigma2, ~k^4).
f8a8a6cNonlinear steel sheet gong simulator with WebGPUJeremy Magland 196 let vel = vec3f(velIn[k], velIn[n + k], velIn[2u * n + k]);
ebf81f1Add an optional k^4 loss term (sigma2), off by defaultJeremy Magland 197 let lap = lapDW(i, j);
f8a8a6cNonlinear steel sheet gong simulator with WebGPUJeremy Magland 198 acc -= 2.0 * SP.sigma0 * vel;
199 acc.z += 2.0 * SP.sigma1 * lap;
ebf81f1Add an optional k^4 loss term (sigma2), off by defaultJeremy Magland 200 if (SP.sigma2 > 0.0) {
201 var bih = 0.0;
202 if (INS(i - 1, j)) { bih += (lapDW(i - 1, j) - lap) * ihx2; }
203 if (INS(i + 1, j)) { bih += (lapDW(i + 1, j) - lap) * ihx2; }
204 if (INS(i, j - 1)) { bih += (lapDW(i, j - 1) - lap) * ihy2; }
205 if (INS(i, j + 1)) { bih += (lapDW(i, j + 1) - lap) * ihy2; }
206 acc.z -= 2.0 * SP.sigma2 * bih;
207 }
bfde5e4High-frequency ringing: in-plane loss, band-limited listener sum, damping defaultsJeremy Magland 208 // In-plane (stretching) waves lie at several kHz; they get their own loss.
209 acc.x -= 2.0 * SP.sigmaIn * vel.x;
210 acc.y -= 2.0 * SP.sigmaIn * vel.y;
f8a8a6cNonlinear steel sheet gong simulator with WebGPUJeremy Magland 211
adc495dFinite finger release: ramp the recorded reaction force to zero over a release timeJeremy Magland 212 // Finger letting go: the reaction recorded at the end of the hold, ramped to zero.
213 if ((SP.flags & 4u) != 0u) { acc.z += react[k] * SP.releaseGain; }
10998acDrumstick strike excitation (Hertz contact with a free mass)Jeremy Magland 215 // Drumstick: Hertz contact between a free mass and the contact patch.
216 if ((SP.flags & 8u) != 0u) {
217 let sIn = n + 4u * SP.parity;
218 let sx = react[sIn]; let sv = react[sIn + 1u]; let sMax = react[sIn + 3u];
219 let wc = posIn[2u * n + P.strikeNode];
220 let pen = wc - sx; // the stick sits on the +z side and moves toward -z
221 var F = 0.0;
222 if (pen > 0.0) { F = P.K * pow(pen, P.alpha); }
223 let dxs = nodeX(i) - P.strikeX; let dys = nodeY(j) - P.strikeY;
224 if (dxs * dxs + dys * dys <= P.strikeR2 || k == P.strikeNode) { acc.z -= F / P.patchMass; }
225 if (k == P.strikeNode) {
226 let sOut = n + 4u * (1u - SP.parity);
227 let svn = sv + P.dt * F / P.stickMass;
228 react[sOut] = sx + P.dt * svn;
229 react[sOut + 1u] = svn;
230 react[sOut + 2u] = F;
231 react[sOut + 3u] = max(sMax, F);
232 }
233 }
f8a8a6cNonlinear steel sheet gong simulator with WebGPUJeremy Magland 235 var vNew = vel + P.dt * acc;
bfde5e4High-frequency ringing: in-plane loss, band-limited listener sum, damping defaultsJeremy Magland 236 if (SP.audioSlot >= 0) { hist[u32(SP.audioSlot) * n + k] = vNew.z; }
f8a8a6cNonlinear steel sheet gong simulator with WebGPUJeremy Magland 237 let pOld = vec3f(posIn[k], posIn[n + k], posIn[2u * n + k]);
238 var pNew = pOld + P.dt * vNew;
239 // Compensated update of w: the increment dt * v can be far below the
240 // float32 resolution of w, so the rounding error is carried separately.
241 let y = P.dt * vNew.z + posIn[3u * n + k];
242 let t = pOld.z + y;
243 var lo = (pOld.z - t) + y;
244 pNew.z = t;
245 if ((SP.flags & 1u) != 0u) {
246 let dxp = nodeX(i) - P.pushX; let dyp = nodeY(j) - P.pushY;
adc495dFinite finger release: ramp the recorded reaction force to zero over a release timeJeremy Magland 247 let inPatch = dxp * dxp + dyp * dyp <= P.pushR2;
248 if ((SP.flags & 2u) != 0u) {
249 // Acceleration the constraint has to supply to hold this node.
250 var r = 0.0;
251 if (inPatch) { r = (SP.pushV - vNew.z) / P.dt; }
252 react[k] = r;
253 }
254 if (inPatch) { pNew.z = SP.pushW; vNew.z = SP.pushV; lo = 0.0; }
256 posOut[k] = pNew.x; posOut[n + k] = pNew.y; posOut[2u * n + k] = pNew.z; posOut[3u * n + k] = lo;
257 velOut[k] = vNew.x; velOut[n + k] = vNew.y; velOut[2u * n + k] = vNew.z;
258}