d8f9e54Clean rebuild of the steel sheet simulatorJeremy 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]
a2a1b68Drumstick by default, hit sequences, 170 GPa steelJeremy Magland 6// A workgroup of 16x16 threads advances a 12x12 tile of nodes in three
7// phases separated by barriers. (1) Every thread loads one entry of the tile
8// plus its 2-node halo (16x16 values of w, u, v, w0, w_t and the node
9// weights) into shared memory. (2) The threads compute, once each and again
10// into shared memory, the membrane forces of the 13x13 cells around the tile
11// and the weighted bending moments and mixed differences that the bending
12// stencil gathers. (3) Each of the 144 tile nodes evaluates its update from
13// shared memory alone. On grids of a few thousand nodes the kernel is bound
14// by the latency of one workgroup rather than by throughput, so the shared
15// phases, which replace a tenfold per-node recomputation of the stencil
16// ingredients, are what make a step cheap.
17//
18// Nodes outside the sheet (weight 0) are skipped: a second difference exists
19// only where the node and both neighbors belong to the sheet, and a cell
d8f9e54Clean rebuild of the steel sheet simulatorJeremy Magland 20// contributes only when its four corners do. That reproduces the free edge of
21// the rectangle and gives a staircase free edge for an ellipse.
23@group(0) @binding(0) var<uniform> P: Params;
24@group(0) @binding(1) var<uniform> SP: StepParams;
25@group(0) @binding(2) var<storage, read> posIn: array<f32>;
26@group(0) @binding(3) var<storage, read> velIn: array<f32>;
27@group(0) @binding(4) var<storage, read_write> posOut: array<f32>;
28@group(0) @binding(5) var<storage, read_write> velOut: array<f32>;
29@group(0) @binding(6) var<storage, read_write> hist: array<f32>; // [slot * n + node]: transverse velocity at the audio sample instants
30// Constraint reaction (acceleration of w) per node, followed by two 4-float
31// slots of drumstick state (position, velocity, force, max force).
32@group(0) @binding(7) var<storage, read_write> react: array<f32>;
33@group(0) @binding(8) var<storage, read> w0: array<f32>; // stress-free initial shape
a2a1b68Drumstick by default, hit sequences, 170 GPa steelJeremy Magland 35const TILE = 12; // nodes per workgroup along each axis (must match simulation.ts)
a2a1b68Drumstick by default, hit sequences, 170 GPa steelJeremy Magland 37const TW = 16; // TILE + 2 * HALO: the loaded tile, one thread per entry
38const CT = 13; // TILE + 1: cells adjacent to the tile's nodes
39const MT = 14; // TILE + 2: nodes at which the bending moments are needed
40// Tile with halo, local index (j - oy) * TW + (i - ox).
41var<workgroup> tw: array<f32, 256>; // w
42var<workgroup> tu: array<f32, 256>; // u
43var<workgroup> tv: array<f32, 256>; // v
44var<workgroup> tw0: array<f32, 256>; // stress-free shape
45var<workgroup> tdw: array<f32, 256>; // w_t
46var<workgroup> tn: array<u32, 64>; // node weights in quarters (0..4), four to a word
47// Membrane force contributions of the CT x CT cells, 4 nodes x 3 components each.
48var<workgroup> cellF: array<f32, 2028>;
49// Weighted bending moments at the MT x MT inner nodes and mixed differences at the CT x CT cells.
50var<workgroup> ma: array<f32, 196>;
51var<workgroup> mb: array<f32, 196>;
52var<workgroup> cc: array<f32, 169>;
d8f9e54Clean rebuild of the steel sheet simulatorJeremy Magland 53
54// Tile origin in global node coordinates (set once per invocation).
55var<private> ox: i32;
56var<private> oy: i32;
57var<private> ihx2: f32;
58var<private> ihy2: f32;
59var<private> ihxy: f32;
61fn id(i: i32, j: i32) -> u32 { return u32(j) * P.nx + u32(i); }
62fn lidx(i: i32, j: i32) -> u32 { return u32(j - oy) * u32(TW) + u32(i - ox); }
a2a1b68Drumstick by default, hit sequences, 170 GPa steelJeremy Magland 63fn midx(i: i32, j: i32) -> u32 { return u32(j - oy - 1) * u32(MT) + u32(i - ox - 1); }
64fn cidx(ci: i32, cj: i32) -> u32 { return u32(cj - oy - 1) * u32(CT) + u32(ci - ox - 1); }
d8f9e54Clean rebuild of the steel sheet simulatorJeremy Magland 65// Tile lookups with global indices; the tile covers [ox, ox + TW) x [oy, oy + TW).
66fn W(i: i32, j: i32) -> f32 { return tw[lidx(i, j)]; }
a2a1b68Drumstick by default, hit sequences, 170 GPa steelJeremy Magland 67fn U(i: i32, j: i32) -> f32 { return tu[lidx(i, j)]; }
68fn V(i: i32, j: i32) -> f32 { return tv[lidx(i, j)]; }
69fn W0(i: i32, j: i32) -> f32 { return tw0[lidx(i, j)]; }
70fn DW(i: i32, j: i32) -> f32 { return tdw[lidx(i, j)]; }
d8f9e54Clean rebuild of the steel sheet simulatorJeremy Magland 71fn NWq(i: i32, j: i32) -> u32 { let l = lidx(i, j); return (tn[l >> 2u] >> ((l & 3u) * 8u)) & 0xFFu; }
72fn NW(i: i32, j: i32) -> f32 { return 0.25 * f32(NWq(i, j)); }
73fn INS(i: i32, j: i32) -> bool { return NWq(i, j) != 0u; }
74fn insCell(ci: i32, cj: i32) -> bool { return INS(ci, cj) && INS(ci + 1, cj) && INS(ci + 1, cj + 1) && INS(ci, cj + 1); }
76// Graph Laplacian of the transverse velocity at node (i, j), over neighbors in the sheet.
77fn lapDW(i: i32, j: i32) -> f32 {
78 let c = DW(i, j);
79 var lap = 0.0;
80 if (INS(i - 1, j)) { lap += (DW(i - 1, j) - c) * ihx2; }
81 if (INS(i + 1, j)) { lap += (DW(i + 1, j) - c) * ihx2; }
82 if (INS(i, j - 1)) { lap += (DW(i, j - 1) - c) * ihy2; }
83 if (INS(i, j + 1)) { lap += (DW(i, j + 1) - c) * ihy2; }
84 return lap;
85}
87// Weighted bending moments at nodes; zero where the second difference is undefined.
88fn Ma(i: i32, j: i32) -> f32 {
89 if (!(INS(i - 1, j) && INS(i, j) && INS(i + 1, j))) { return 0.0; }
90 let a = (W(i + 1, j) - 2.0 * W(i, j) + W(i - 1, j)) * ihx2;
91 var b = 0.0;
92 if (INS(i, j - 1) && INS(i, j + 1)) { b = (W(i, j + 1) - 2.0 * W(i, j) + W(i, j - 1)) * ihy2; }
93 return NW(i, j) * (a + P.nu * b);
94}
95fn Mb(i: i32, j: i32) -> f32 {
96 if (!(INS(i, j - 1) && INS(i, j) && INS(i, j + 1))) { return 0.0; }
97 let b = (W(i, j + 1) - 2.0 * W(i, j) + W(i, j - 1)) * ihy2;
98 var a = 0.0;
99 if (INS(i - 1, j) && INS(i + 1, j)) { a = (W(i + 1, j) - 2.0 * W(i, j) + W(i - 1, j)) * ihx2; }
100 return NW(i, j) * (b + P.nu * a);
101}
102// Mixed second difference at a cell center; zero outside the sheet.
103fn Cc(ci: i32, cj: i32) -> f32 {
104 if (!insCell(ci, cj)) { return 0.0; }
105 return (W(ci + 1, cj + 1) - W(ci + 1, cj) - W(ci, cj + 1) + W(ci, cj)) * ihxy;
106}
108// Membrane forces of the bilinear cell (ci, cj) on its four nodes (2x2 Gauss),
109// written to cellF at local cell index lc.
110fn cellForces(ci: i32, cj: i32, lc: u32) {
111 let base = lc * 12u;
112 if (!insCell(ci, cj)) {
113 for (var q = 0u; q < 12u; q++) { cellF[base + q] = 0.0; }
114 return;
115 }
116 let Un = vec4f(U(ci, cj), U(ci + 1, cj), U(ci + 1, cj + 1), U(ci, cj + 1));
117 let Vn = vec4f(V(ci, cj), V(ci + 1, cj), V(ci + 1, cj + 1), V(ci, cj + 1));
118 let W0n = vec4f(W0(ci, cj), W0(ci + 1, cj), W0(ci + 1, cj + 1), W0(ci, cj + 1));
119 let Wn = vec4f(W(ci, cj), W(ci + 1, cj), W(ci + 1, cj + 1), W(ci, cj + 1)) + W0n; // total shape
120 let q = 0.25 * P.hx * P.hy;
121 let ihx = 1.0 / P.hx; let ihy = 1.0 / P.hy;
122 let g3 = 0.57735026918962576;
123 var Fu = vec4f(0.0); var Fv = vec4f(0.0); var Fw = vec4f(0.0);
124 for (var gp = 0; gp < 4; gp++) {
125 var xi = -g3; if ((gp & 1) != 0) { xi = g3; }
126 var eta = -g3; if ((gp & 2) != 0) { eta = g3; }
127 let dNdx = vec4f(-(1.0 - eta), 1.0 - eta, 1.0 + eta, -(1.0 + eta)) * (0.5 * ihx);
128 let dNdy = vec4f(-(1.0 - xi), -(1.0 + xi), 1.0 + xi, 1.0 - xi) * (0.5 * ihy);
129 let ux = dot(dNdx, Un); let uy = dot(dNdy, Un);
130 let vx = dot(dNdx, Vn); let vy = dot(dNdy, Vn);
131 let wx = dot(dNdx, Wn); let wy = dot(dNdy, Wn);
132 let w0x = dot(dNdx, W0n); let w0y = dot(dNdy, W0n);
133 // Marguerre strains: total slope minus the stress-free slope.
134 let exx = ux + 0.5 * (wx * wx - w0x * w0x);
135 let eyy = vy + 0.5 * (wy * wy - w0y * w0y);
136 let gxy = uy + vx + wx * wy - w0x * w0y;
137 let Nxx = P.C * (exx + P.nu * eyy);
138 let Nyy = P.C * (eyy + P.nu * exx);
139 let Nxy = P.G * gxy;
140 Fu -= q * (Nxx * dNdx + Nxy * dNdy);
141 Fv -= q * (Nyy * dNdy + Nxy * dNdx);
142 Fw -= q * ((Nxx * wx + Nxy * wy) * dNdx + (Nxy * wx + Nyy * wy) * dNdy);
143 }
144 for (var a = 0u; a < 4u; a++) {
145 cellF[base + a * 3u] = Fu[a];
146 cellF[base + a * 3u + 1u] = Fv[a];
147 cellF[base + a * 3u + 2u] = Fw[a];
148 }
149}
151// True for the sheet nodes of the drumstick contact patch (same rule as contactPatch in params.ts).
152fn inContactPatch(i: i32, j: i32) -> bool {
153 if (nodeWeight(i, j) == 0.0) { return false; }
154 let dxs = nodeX(i) - P.strikeX; let dys = nodeY(j) - P.strikeY;
155 return dxs * dxs + dys * dys <= P.strikeR2 || id(i, j) == P.strikeNode;
156}
158// Stored membrane force on local node `a` of cell (ci, cj).
159fn cellForce(ci: i32, cj: i32, a: u32) -> vec3f {
a2a1b68Drumstick by default, hit sequences, 170 GPa steelJeremy Magland 160 let base = cidx(ci, cj) * 12u + a * 3u;
d8f9e54Clean rebuild of the steel sheet simulatorJeremy Magland 161 return vec3f(cellF[base], cellF[base + 1u], cellF[base + 2u]);
162}
164@compute @workgroup_size(16, 16)
a2a1b68Drumstick by default, hit sequences, 170 GPa steelJeremy Magland 165fn main(@builtin(local_invocation_id) lid: vec3u, @builtin(workgroup_id) wid: vec3u) {
167 ox = i32(wid.x) * TILE - HALO;
168 oy = i32(wid.y) * TILE - HALO;
169 ihx2 = 1.0 / (P.hx * P.hx); ihy2 = 1.0 / (P.hy * P.hy); ihxy = 1.0 / (P.hx * P.hy);
170 let lt = lid.y * 16u + lid.x;
a2a1b68Drumstick by default, hit sequences, 170 GPa steelJeremy Magland 171
172 // Phase 1: one entry of the tile with halo per thread.
173 {
174 let gi = ox + i32(lid.x); let gj = oy + i32(lid.y);
175 var ww = 0.0; var uu = 0.0; var vv = 0.0; var z0 = 0.0; var dw = 0.0;
176 if (gi >= 0 && gi < i32(P.nx) && gj >= 0 && gj < i32(P.ny)) {
177 let k = id(gi, gj);
178 ww = posIn[2u * n + k]; uu = posIn[k]; vv = posIn[n + k]; z0 = w0[k]; dw = velIn[2u * n + k];
179 }
180 tw[lt] = ww; tu[lt] = uu; tv[lt] = vv; tw0[lt] = z0; tdw[lt] = dw;
181 if (lt < 64u) {
182 var word = 0u;
183 for (var q = 0u; q < 4u; q++) {
184 let e = lt * 4u + q;
185 let ei = ox + i32(e % u32(TW)); let ej = oy + i32(e / u32(TW));
186 var wq = 0u;
187 if (ei >= 0 && ei < i32(P.nx) && ej >= 0 && ej < i32(P.ny)) { wq = u32(round(4.0 * nw[id(ei, ej)])); }
188 word |= wq << (q * 8u);
189 }
190 tn[lt] = word;
192 }
193 workgroupBarrier();
a2a1b68Drumstick by default, hit sequences, 170 GPa steelJeremy Magland 194
195 // Phase 2: membrane forces of the CT x CT cells, bending moments at the MT x MT
196 // inner nodes, mixed differences at the cells; one item per thread each.
197 if (lt < 169u) {
198 let ci = ox + 1 + i32(lt % u32(CT)); let cj = oy + 1 + i32(lt / u32(CT));
199 cellForces(ci, cj, lt);
200 cc[lt] = Cc(ci, cj);
201 }
202 if (lt < 196u) {
203 let mi = ox + 1 + i32(lt % u32(MT)); let mj = oy + 1 + i32(lt / u32(MT));
204 ma[lt] = Ma(mi, mj);
205 mb[lt] = Mb(mi, mj);
207 workgroupBarrier();
a2a1b68Drumstick by default, hit sequences, 170 GPa steelJeremy Magland 209 // Phase 3: the update of the TILE x TILE nodes (no barriers past this point).
210 if (lt >= 144u) { return; }
211 let i = ox + HALO + i32(lt % u32(TILE)); let j = oy + HALO + i32(lt / u32(TILE));
d8f9e54Clean rebuild of the steel sheet simulatorJeremy Magland 212 if (i >= i32(P.nx) || j >= i32(P.ny)) { return; }
213 let k = id(i, j);
214 let weight = NW(i, j);
215 if (weight == 0.0) {
216 // Outside the sheet: the node stays at rest.
217 posOut[k] = 0.0; posOut[n + k] = 0.0; posOut[2u * n + k] = 0.0; posOut[3u * n + k] = 0.0;
218 velOut[k] = 0.0; velOut[n + k] = 0.0; velOut[2u * n + k] = 0.0;
219 if ((SP.flags & 2u) != 0u) { react[k] = 0.0; }
220 return;
221 }
223 // Elastic forces: bending from the discrete Kirchhoff energy, membrane from the cells.
a2a1b68Drumstick by default, hit sequences, 170 GPa steelJeremy Magland 224 let bend = (ma[midx(i - 1, j)] - 2.0 * ma[midx(i, j)] + ma[midx(i + 1, j)]) * ihx2
225 + (mb[midx(i, j - 1)] - 2.0 * mb[midx(i, j)] + mb[midx(i, j + 1)]) * ihy2
226 + 2.0 * (1.0 - P.nu) * (cc[cidx(i, j)] - cc[cidx(i - 1, j)] - cc[cidx(i, j - 1)] + cc[cidx(i - 1, j - 1)]) * ihxy;
d8f9e54Clean rebuild of the steel sheet simulatorJeremy Magland 227 var F = vec3f(0.0, 0.0, -P.D * P.hx * P.hy * bend);
228 F += cellForce(i, j, 0u);
229 F += cellForce(i - 1, j, 1u);
230 F += cellForce(i - 1, j - 1, 2u);
231 F += cellForce(i, j - 1, 3u);
232 let mass = P.rhoh * P.hx * P.hy * weight;
233 var acc = F * (1.0 / mass);
a2a1b68Drumstick by default, hit sequences, 170 GPa steelJeremy Magland 234 // Optional scaling of the in-plane inertia (lengthens the stable time step).
235 acc.x *= P.invInPlaneMass; acc.y *= P.invInPlaneMass;
d8f9e54Clean rebuild of the steel sheet simulatorJeremy Magland 236
237 // Losses: frequency-independent on all components, the graph Laplacian of the
238 // transverse velocity (~k^2), and an extra frequency-independent in-plane loss.
239 let vel = vec3f(velIn[k], velIn[n + k], velIn[2u * n + k]);
240 acc -= 2.0 * SP.sigma0 * vel;
241 acc.z += 2.0 * SP.sigma1 * lapDW(i, j);
242 acc.x -= 2.0 * SP.sigmaIn * vel.x;
243 acc.y -= 2.0 * SP.sigmaIn * vel.y;
245 // Finger letting go: the reaction recorded at the end of the hold, ramped to zero.
246 if ((SP.flags & 4u) != 0u) { acc.z += react[k] * SP.releaseGain; }
248 // Drumstick: Hertz contact between a free mass and the contact patch (the
249 // sheet nodes within the tip radius, plus the nearest one). The contact
250 // coordinate is the mass-weighted mean displacement of the patch, and the
251 // force is spread over the patch in proportion to the node masses, so the
252 // pair is conjugate and the contact has a potential energy.
253 if ((SP.flags & 8u) != 0u) {
282f706Single excitation panel; rigid-body projection optional while ringingJeremy Magland 254 let sIn = n + 4u * SP.parity;
255 // Inertia relief: the force of the previous step, applied to the whole sheet
256 // as the rigid acceleration field that balances the patch load and its
257 // moments about the center, so the sheet gains no net momentum (the thread
258 // would take it up in reality). The one-step lag integrates to zero.
259 let Fprev = react[sIn + 2u];
260 acc.z += Fprev * (1.0 / P.mass + nodeY(j) * P.patchCY / P.Ix + nodeX(i) * P.patchCX / P.Iy);
262 let sx = react[sIn]; let sv = react[sIn + 1u]; let sMax = react[sIn + 3u];
263 var wc = 0.0;
264 for (var pj = i32(P.patchJ0); pj <= i32(P.patchJ1); pj++) {
265 for (var pi = i32(P.patchI0); pi <= i32(P.patchI1); pi++) {
266 if (inContactPatch(pi, pj)) { wc += P.rhoh * P.hx * P.hy * nodeWeight(pi, pj) * posIn[2u * n + id(pi, pj)]; }
267 }
268 }
269 wc /= P.patchMass;
270 let pen = wc - sx; // the stick sits on the +z side and moves toward -z
271 var Fc = 0.0;
272 if (pen > 0.0) { Fc = P.K * pow(pen, P.alpha); }
273 acc.z -= Fc / P.patchMass;
274 if (k == P.strikeNode) {
275 let sOut = n + 4u * (1u - SP.parity);
a2a1b68Drumstick by default, hit sequences, 170 GPa steelJeremy Magland 276 if ((SP.flags & 16u) != 0u) {
277 // Next hit of a series: a fresh stick arrives at the patch, wherever
278 // it is at this instant, with the scheduled speed; the previous
279 // stick's rebound is discarded.
280 react[sOut] = wc;
281 react[sOut + 1u] = -SP.strikeV;
282 react[sOut + 2u] = 0.0;
283 react[sOut + 3u] = sMax;
284 } else {
285 let svn = sv + P.dt * Fc / P.stickMass;
286 react[sOut] = sx + P.dt * svn;
287 react[sOut + 1u] = svn;
288 react[sOut + 2u] = Fc;
289 react[sOut + 3u] = max(sMax, Fc);
290 }
292 }
293 }
295 var vNew = vel + P.dt * acc;
296 if (SP.audioSlot >= 0) { hist[u32(SP.audioSlot) * n + k] = vNew.z; }
297 let pOld = vec3f(posIn[k], posIn[n + k], posIn[2u * n + k]);
298 var pNew = pOld + P.dt * vNew;
299 // Compensated update of w: the increment dt * v can be far below the
300 // float32 resolution of w, so the rounding error is carried separately.
301 let y = P.dt * vNew.z + posIn[3u * n + k];
302 let t = pOld.z + y;
303 var lo = (pOld.z - t) + y;
304 pNew.z = t;
305 if ((SP.flags & 1u) != 0u) {
306 let dxp = nodeX(i) - P.pushX; let dyp = nodeY(j) - P.pushY;
307 let inPatch = dxp * dxp + dyp * dyp <= P.pushR2;
308 if ((SP.flags & 2u) != 0u) {
309 // Acceleration the constraint has to supply to hold this node.
310 var r = 0.0;
311 if (inPatch) { r = (SP.pushV - vNew.z) / P.dt; }
312 react[k] = r;
313 }
314 if (inPatch) { pNew.z = SP.pushW; vNew.z = SP.pushV; lo = 0.0; }
315 }
316 posOut[k] = pNew.x; posOut[n + k] = pNew.y; posOut[2u * n + k] = pNew.z; posOut[3u * n + k] = lo;
317 velOut[k] = vNew.x; velOut[n + k] = vNew.y; velOut[2u * n + k] = vNew.z;
318}