/ proof-of-concept / cymsimweb
Sign in
proof-of-concept / cymsimweb
2b313b8
Branches
main
cymsimweb / src / gpu / shaders.ts
305 lines · 15.0 KBCodeBlameHistory
13f4dbdexport .cymsim fileJeremy Magland 1// WGSL kernels of the time loop. Four entry points: the element pass (forces, sigma1 term,
2// contact forces of the hits), the node pass (gather, contact, lagged rigid-mode projection,
3// leapfrog update, partial sums, exciter update, and the normal acceleration with its decimation
4// accumulators), the gather pass (element means of a finished audio sample into the history,
56d33ddGPU: decimate once per audio sample and spread the listener sum over all elementsJeremy Magland 5// then the retarded-time sum at the listener, one partial sum per workgroup) and the capture pass
6// (displacement frames). Sizes
1e29264cymsimweb: WebGPU companion of cymsim with YAML/form editor, live spectrogram and playbackJeremy Magland 7// and constants are baked into the source. See THEORY.md for the buffer layouts and the ordering
8// argument behind the two step counters.
10export interface ShaderConstants {
13f4dbdexport .cymsim fileJeremy Magland 11 N: number; E: number; K: number; P: number; M: number; W: number; D: number; G: number; T: number;
12 accRows: number; dt: number; sigma0: number; nl: number;
2b313b8GPU: field-major set-up data, symmetric bending stiffness, shared rigid-mode sumsJeremy Magland 15// Per-element and per-node data are stored field-major (field f of element e at f * E + e), so that
16// neighbouring threads read neighbouring words.
17export const ELEM_STRIDE = 58; // fields: Gx[3] Gy[3] A z0x z0y C11 C12 C33 cd Ke[45] (upper triangle, row by row)
18export const NODE_STRIDE = 43; // fields: Minv[5] M[5] Rm[6][5] normal[3]; the decimation taps follow the nodes
1e29264cymsimweb: WebGPU companion of cymsim with YAML/form editor, live spectrogram and playbackJeremy Magland 19export const HIT_HEADER = 12; // kPoint alpha mass gap velocity zs vs F Fmax pad[3]
20export const HITU_HEADER = 4; // releaseStep ncont pad[2]
21export const SYNC_HEADER = 4; // stepA stepB pad[2], then partials[2][G][12]
13f4dbdexport .cymsim fileJeremy Magland 22export const RADCOEF_STRIDE = 36; // coefF coefB i0F i0B tapsF[16] tapsB[16]
56d33ddGPU: decimate once per audio sample and spread the listener sum over all elementsJeremy Magland 24export const GATHER_ELEMS = 256; // elements per gather workgroup
13f4dbdexport .cymsim fileJeremy Magland 25export const STATE_FLOATS = 13; // per node: q[5] v[5] vPrev[3]
26/** Chunk uniform: base size frameBase every accBase accRows pad pad (eight 32-bit words). */
27export const CHUNK_WORDS = 8;
29function f(x: number): string {
30 const s = x.toPrecision(9);
31 return /[.e]/.test(s) ? s : s + ".0";
2b313b8GPU: field-major set-up data, symmetric bending stiffness, shared rigid-mode sumsJeremy Magland 34/** -K_e U added to the bending forces, from the 45 stored entries of the symmetric K_e, written out
35 * in full: a loop over the packed triangle indexes private arrays dynamically, which was slower. */
36function bendingProduct(): string {
37 const lines: string[] = [], terms: string[][] = Array.from({ length: 9 }, () => []);
38 let f = 13;
39 for (let i = 0; i < 9; i++) for (let j = i; j < 9; j++) {
40 lines.push(` let k${i}${j} = elemF[${f++}u * E + e];`);
41 terms[i].push(`k${i}${j} * U[${j}]`);
42 if (j !== i) terms[j].push(`k${i}${j} * U[${i}]`);
43 }
44 for (let i = 0; i < 9; i++) lines.push(` out[${5 * Math.floor(i / 3) + 2 + (i % 3)}u] -= ${terms[i].join(" + ")};`);
45 return lines.join("\n");
1e29264cymsimweb: WebGPU companion of cymsim with YAML/form editor, live spectrogram and playbackJeremy Magland 48export function shaderSource(c: ShaderConstants): string {
49 const HF = HIT_HEADER + 5 * c.P, HU = HITU_HEADER + 3 * c.P;
50 return /* wgsl */ `
51const N: u32 = ${c.N}u; const E: u32 = ${c.E}u; const K: u32 = ${c.K}u; const P: u32 = ${c.P}u;
13f4dbdexport .cymsim fileJeremy Magland 52const M: u32 = ${c.M}u; const W: u32 = ${c.W}u; const D: u32 = ${c.D}u; const G: u32 = ${c.G}u; const T: u32 = ${c.T}u;
1e29264cymsimweb: WebGPU companion of cymsim with YAML/form editor, live spectrogram and playbackJeremy Magland 53const DT: f32 = ${f(c.dt)}; const DAMP0: f32 = ${f(1 / (1 + 2 * c.sigma0 * c.dt))}; const NL: f32 = ${f(c.nl)};
2b313b8GPU: field-major set-up data, symmetric bending stiffness, shared rigid-mode sumsJeremy Magland 54const HF: u32 = ${HF}u; const HU: u32 = ${HU}u; const RC: u32 = ${RADCOEF_STRIDE}u;
13f4dbdexport .cymsim fileJeremy Magland 55const OFF_V: u32 = 5u * N; const OFF_VWP: u32 = 10u * N;
56d33ddGPU: decimate once per audio sample and spread the listener sum over all elementsJeremy Magland 56const OFF_NHIST: u32 = W * N; const OFF_AHIST: u32 = (W + D) * N; const OFF_ACCOUT: u32 = (W + D + M) * N;
57const GE: u32 = ${Math.ceil(c.E / GATHER_ELEMS)}u;
2b313b8GPU: field-major set-up data, symmetric bending stiffness, shared rigid-mode sumsJeremy Magland 58const OFF_TAPS: u32 = ${NODE_STRIDE}u * N; const OFF_PART: u32 = ${SYNC_HEADER}u;
13f4dbdexport .cymsim fileJeremy Magland 60struct Chunk { base: i32, size: u32, frameBase: u32, every: u32, accBase: i32, accRows: u32, pad0: u32, pad1: u32 }
62// ------------------------------------------------------------------ element pass
63@group(0) @binding(0) var<storage, read> elemF: array<f32>;
64@group(0) @binding(1) var<storage, read> elemU: array<u32>;
65@group(0) @binding(2) var<storage, read> stateR: array<f32>;
66@group(0) @binding(3) var<storage, read_write> elemForce: array<f32>;
67@group(0) @binding(4) var<storage, read_write> hitsF: array<f32>;
68@group(0) @binding(5) var<storage, read> hitsUR: array<u32>;
69@group(0) @binding(6) var<storage, read_write> syncE: array<f32>;
71@compute @workgroup_size(64)
72fn elementPass(@builtin(global_invocation_id) gid: vec3<u32>) {
73 let id = gid.x;
74 let step = bitcast<u32>(syncE[0]);
75 if (id == 0u) { syncE[1] = syncE[0]; }
76 if (id < E) {
77 let e = id;
78 let n0 = elemU[3u * e]; let n1 = elemU[3u * e + 1u]; let n2 = elemU[3u * e + 2u];
79 var qe: array<f32, 15>;
80 for (var d = 0u; d < 5u; d++) { qe[d] = stateR[5u * n0 + d]; qe[5u + d] = stateR[5u * n1 + d]; qe[10u + d] = stateR[5u * n2 + d]; }
2b313b8GPU: field-major set-up data, symmetric bending stiffness, shared rigid-mode sumsJeremy Magland 81 let gx = vec3<f32>(elemF[e], elemF[1u * E + e], elemF[2u * E + e]);
82 let gy = vec3<f32>(elemF[3u * E + e], elemF[4u * E + e], elemF[5u * E + e]);
83 let A = elemF[6u * E + e]; let z0x = elemF[7u * E + e]; let z0y = elemF[8u * E + e];
84 let C11 = elemF[9u * E + e]; let C12 = elemF[10u * E + e]; let C33 = elemF[11u * E + e]; let cd = elemF[12u * E + e];
1e29264cymsimweb: WebGPU companion of cymsim with YAML/form editor, live spectrogram and playbackJeremy Magland 85 let uu = vec3<f32>(qe[0], qe[5], qe[10]); let vv = vec3<f32>(qe[1], qe[6], qe[11]); let ww = vec3<f32>(qe[2], qe[7], qe[12]);
86 let ux = dot(gx, uu); let uy = dot(gy, uu); let vx = dot(gx, vv); let vy = dot(gy, vv); let wx = dot(gx, ww); let wy = dot(gy, ww);
87 let ax = NL * wx + z0x; let ay = NL * wy + z0y;
88 let exx = ux + NL * 0.5 * wx * wx + z0x * wx;
89 let eyy = vy + NL * 0.5 * wy * wy + z0y * wy;
90 let gxy = uy + vx + NL * wx * wy + z0x * wy + z0y * wx;
91 let Nxx = C11 * exx + C12 * eyy; let Nyy = C12 * exx + C11 * eyy; let Nxy = C33 * gxy;
92 let vw = vec3<f32>(stateR[OFF_V + 5u * n0 + 2u], stateR[OFF_V + 5u * n1 + 2u], stateR[OFF_V + 5u * n2 + 2u]);
93 let vwx = dot(gx, vw); let vwy = dot(gy, vw);
94 var out: array<f32, 15>;
95 for (var i = 0u; i < 3u; i++) {
96 let g_x = gx[i]; let g_y = gy[i];
97 out[5u * i] = -A * (Nxx * g_x + Nxy * g_y);
98 out[5u * i + 1u] = -A * (Nyy * g_y + Nxy * g_x);
99 out[5u * i + 2u] = -A * (Nxx * ax * g_x + Nyy * ay * g_y + Nxy * (ay * g_x + ax * g_y)) - cd * (g_x * vwx + g_y * vwy);
100 out[5u * i + 3u] = 0.0; out[5u * i + 4u] = 0.0;
101 }
102 var U: array<f32, 9>;
103 for (var i = 0u; i < 3u; i++) { U[3u * i] = qe[5u * i + 2u]; U[3u * i + 1u] = qe[5u * i + 3u]; U[3u * i + 2u] = qe[5u * i + 4u]; }
1e29264cymsimweb: WebGPU companion of cymsim with YAML/form editor, live spectrogram and playbackJeremy Magland 105 for (var i = 0u; i < 15u; i++) { elemForce[15u * e + i] = out[i]; }
106 } else if (id < E + K) {
107 // contact forces of hit k at this step
108 let k = id - E; let hf = k * HF; let hu = k * HU;
109 var wp: array<f32, ${c.P}>;
110 var wmax = -1.0e30;
111 for (var p = 0u; p < P; p++) {
112 var w = hitsF[hf + ${HIT_HEADER}u + 3u * P + p];
113 for (var i = 0u; i < 3u; i++) { w += hitsF[hf + ${HIT_HEADER}u + 3u * p + i] * stateR[5u * hitsUR[hu + ${HITU_HEADER}u + 3u * p + i] + 2u]; }
114 wp[p] = w; wmax = max(wmax, w);
115 }
116 var zs = hitsF[hf + 5u];
117 if (step == hitsUR[hu]) { zs = wmax + hitsF[hf + 3u]; hitsF[hf + 5u] = zs; hitsF[hf + 6u] = -hitsF[hf + 4u]; }
118 let isActive = step >= hitsUR[hu];
119 let kPoint = hitsF[hf]; let alpha = hitsF[hf + 1u];
120 var F = 0.0;
121 for (var p = 0u; p < P; p++) {
122 let delta = max(wp[p] - zs, 0.0);
123 var fp = 0.0;
124 if (isActive && delta > 0.0) { fp = kPoint * pow(delta, alpha); }
125 hitsF[hf + ${HIT_HEADER}u + 4u * P + p] = fp; F += fp;
126 }
127 hitsF[hf + 7u] = F;
128 }
131// ------------------------------------------------------------------ node pass
132@group(1) @binding(0) var<storage, read_write> state: array<f32>;
133@group(1) @binding(1) var<storage, read> elemForceR: array<f32>;
134@group(1) @binding(2) var<storage, read> nodeF: array<f32>;
135@group(1) @binding(3) var<storage, read> adj: array<u32>;
136@group(1) @binding(4) var<storage, read_write> hitsFN: array<f32>;
137@group(1) @binding(5) var<storage, read_write> hitsU: array<u32>;
138@group(1) @binding(6) var<storage, read_write> sync: array<f32>;
13f4dbdexport .cymsim fileJeremy Magland 139@group(1) @binding(7) var<storage, read_write> rad: array<f32>;
140@group(1) @binding(8) var<uniform> chunkN: Chunk;
142var<workgroup> red: array<f32, ${12 * WORKGROUP_NODES}>;
2b313b8GPU: field-major set-up data, symmetric bending stiffness, shared rigid-mode sumsJeremy Magland 143var<workgroup> csum: array<f32, 12>;
13f4dbdexport .cymsim fileJeremy Magland 145fn modW(x: i32) -> u32 { return u32(((x % i32(W)) + i32(W)) % i32(W)); }
1e29264cymsimweb: WebGPU companion of cymsim with YAML/form editor, live spectrogram and playbackJeremy Magland 147@compute @workgroup_size(${WORKGROUP_NODES})
148fn nodePass(@builtin(global_invocation_id) gid: vec3<u32>, @builtin(local_invocation_id) lid: vec3<u32>, @builtin(workgroup_id) wid: vec3<u32>) {
149 let nd = gid.x; let t = lid.x;
150 let step = bitcast<u32>(sync[1]);
151 let parity = step & 1u; let prevp = (step + 1u) & 1u;
2b313b8GPU: field-major set-up data, symmetric bending stiffness, shared rigid-mode sumsJeremy Magland 152 // twelve threads sum the previous step's partials of all workgroups; the others read the sums from shared memory
153 if (t < 12u) {
2b313b8GPU: field-major set-up data, symmetric bending stiffness, shared rigid-mode sumsJeremy Magland 155 for (var g = 0u; g < G; g++) { s += sync[OFF_PART + (prevp * G + g) * 12u + t]; }
156 csum[t] = s;
159 var c: array<f32, 12>;
160 for (var k = 0u; k < 12u; k++) { c[k] = csum[k]; }
162 for (var k = 0u; k < 12u; k++) { part[k] = 0.0; }
163 if (nd < N) {
164 var fo: array<f32, 5>;
165 for (var d = 0u; d < 5u; d++) { fo[d] = 0.0; }
166 let a0 = adj[nd]; let a1 = adj[nd + 1u];
167 for (var a = a0; a < a1; a++) {
168 let ent = adj[N + 1u + a]; let e = ent / 3u; let i = ent - 3u * e;
169 for (var d = 0u; d < 5u; d++) { fo[d] += elemForceR[15u * e + 5u * i + d]; }
170 }
171 for (var k = 0u; k < K; k++) {
172 for (var p = 0u; p < P; p++) {
173 for (var i = 0u; i < 3u; i++) {
174 if (hitsU[k * HU + ${HITU_HEADER}u + 3u * p + i] == nd) {
175 fo[2] -= hitsFN[k * HF + ${HIT_HEADER}u + 4u * P + p] * hitsFN[k * HF + ${HIT_HEADER}u + 3u * p + i];
176 }
177 }
178 }
179 }
13f4dbdexport .cymsim fileJeremy Magland 180 var acc = 0.0;
182 let idx = 5u * nd + d;
183 var vd = state[OFF_V + idx]; var qd = state[idx];
2b313b8GPU: field-major set-up data, symmetric bending stiffness, shared rigid-mode sumsJeremy Magland 184 for (var k = 0u; k < 6u; k++) { let r = nodeF[(10u + 5u * k + d) * N + nd]; vd -= c[k] * r; qd -= c[6u + k] * r; }
185 let v1 = (vd + DT * fo[d] * nodeF[d * N + nd]) * DAMP0;
13f4dbdexport .cymsim fileJeremy Magland 186 // normal acceleration of the previous step from consecutive projected velocities
2b313b8GPU: field-major set-up data, symmetric bending stiffness, shared rigid-mode sumsJeremy Magland 187 if (d < 3u) { acc += nodeF[(40u + d) * N + nd] * (vd - state[OFF_VWP + 3u * nd + d]) / DT; state[OFF_VWP + 3u * nd + d] = vd; }
189 state[idx] = q1; state[OFF_V + idx] = v1;
191 for (var k = 0u; k < 6u; k++) { let r = nodeF[(10u + 5u * k + d) * N + nd] * mm; part[k] += r * v1; part[6u + k] += r * q1; }
13f4dbdexport .cymsim fileJeremy Magland 193 // decimation of the acceleration of sub-step s = step - 1 into the W accumulator slots of this node;
194 // slot w belongs to audio sample j0 + w - W/2 + 1, kept at ring position (sample mod W)
195 if (step >= 1u) {
196 let s = step - 1u; let j0 = s / M; let kk = s - j0 * M;
197 if (kk == 0u) {
198 let jE = i32(j0) - i32(W / 2u);
199 let slotE = modW(jE);
200 let val = rad[slotE * N + nd];
201 if (jE >= 0) {
202 rad[OFF_NHIST + (u32(jE) % D) * N + nd] = val;
203 let row = jE - chunkN.accBase;
204 if (chunkN.accRows > 0u && row >= 0 && u32(row) < chunkN.accRows) { rad[OFF_ACCOUT + u32(row) * N + nd] = val; }
205 }
206 rad[slotE * N + nd] = 0.0;
207 }
56d33ddGPU: decimate once per audio sample and spread the listener sum over all elementsJeremy Magland 208 // keep the sub-step's acceleration; once the M values of the sample are in, add them to the W slots at once
209 rad[OFF_AHIST + kk * N + nd] = acc;
210 if (kk == M - 1u) {
211 let r = modW(i32(j0) + 1 - i32(W / 2u));
212 var sl: array<f32, ${c.W}>;
213 for (var w = 0u; w < W; w++) { sl[w] = 0.0; }
214 for (var m = 0u; m < M; m++) {
215 let a = rad[OFF_AHIST + m * N + nd];
216 for (var w = 0u; w < W; w++) { sl[w] += a * nodeF[OFF_TAPS + w * M + m]; }
217 }
218 for (var w = 0u; w < W; w++) { rad[((r + w) % W) * N + nd] += sl[w]; }
219 }
13f4dbdexport .cymsim fileJeremy Magland 220 }
222 for (var k = 0u; k < 12u; k++) { red[t * 12u + k] = part[k]; }
223 workgroupBarrier();
224 for (var s = ${WORKGROUP_NODES / 2}u; s > 0u; s = s >> 1u) {
225 if (t < s) { for (var k = 0u; k < 12u; k++) { red[t * 12u + k] += red[(t + s) * 12u + k]; } }
226 workgroupBarrier();
227 }
228 if (t == 0u) { for (var k = 0u; k < 12u; k++) { sync[OFF_PART + (parity * G + wid.x) * 12u + k] = red[k]; } }
229 if (nd == 0u) {
230 for (var k = 0u; k < K; k++) {
231 let hf = k * HF; let hu = k * HU;
232 let F = hitsFN[hf + 7u];
233 let vs = hitsFN[hf + 6u] + DT * F / hitsFN[hf + 2u];
234 hitsFN[hf + 6u] = vs; hitsFN[hf + 5u] = hitsFN[hf + 5u] + DT * vs;
235 hitsFN[hf + 8u] = max(hitsFN[hf + 8u], F);
236 if (F > 0.0) { hitsU[hu + 1u] = hitsU[hu + 1u] + 1u; }
237 }
238 sync[0] = bitcast<f32>(step + 1u);
239 }
13f4dbdexport .cymsim fileJeremy Magland 242// ------------------------------------------------------------------ gather pass (listener), once per finished audio sample
1e29264cymsimweb: WebGPU companion of cymsim with YAML/form editor, live spectrogram and playbackJeremy Magland 243@group(2) @binding(0) var<storage, read> syncG: array<f32>;
244@group(2) @binding(1) var<storage, read> radR: array<f32>;
245@group(2) @binding(2) var<storage, read> radCoef: array<f32>;
246@group(2) @binding(3) var<storage, read_write> pressure: array<f32>;
247@group(2) @binding(4) var<uniform> chunk: Chunk;
13f4dbdexport .cymsim fileJeremy Magland 248@group(2) @binding(5) var<storage, read_write> elemHist: array<f32>;
249@group(2) @binding(6) var<storage, read> elemUG: array<u32>;
56d33ddGPU: decimate once per audio sample and spread the listener sum over all elementsJeremy Magland 251var<workgroup> redF: array<f32, ${GATHER_ELEMS}>;
252var<workgroup> redB: array<f32, ${GATHER_ELEMS}>;
13f4dbdexport .cymsim fileJeremy Magland 254fn delayed(e: u32, j: i32, i0: i32, tapBase: u32) -> f32 {
255 var s = 0.0;
256 for (var t = 0u; t < T; t++) {
257 let lag = i0 + i32(t) - i32(T / 2u) + 1; let jj = j - lag;
258 if (jj >= 0 && lag >= 0) { s += radCoef[tapBase + t] * elemHist[(u32(jj) % D) * E + e]; }
259 }
260 return s;
56d33ddGPU: decimate once per audio sample and spread the listener sum over all elementsJeremy Magland 263@compute @workgroup_size(${GATHER_ELEMS})
264fn gatherPass(@builtin(local_invocation_id) lid: vec3<u32>, @builtin(workgroup_id) wid: vec3<u32>) {
265 // one thread per element; each workgroup writes its partial sums, which the host adds up
266 let t = lid.x; let e = wid.x * ${GATHER_ELEMS}u + t;
13f4dbdexport .cymsim fileJeremy Magland 267 let step = bitcast<u32>(syncG[0]); // the node pass of step n has run: the finished sample belongs to s = n - 1
268 let j = i32((step - 2u) / M) - i32(W / 2u);
13f4dbdexport .cymsim fileJeremy Magland 271 let hrow = OFF_NHIST + (u32(j) % D) * N;
56d33ddGPU: decimate once per audio sample and spread the listener sum over all elementsJeremy Magland 272 let mean = (radR[hrow + elemUG[3u * e]] + radR[hrow + elemUG[3u * e + 1u]] + radR[hrow + elemUG[3u * e + 2u]]) / 3.0;
273 elemHist[(u32(j) % D) * E + e] = mean;
274 let rc = RC * e;
275 fs = radCoef[rc] * delayed(e, j, i32(radCoef[rc + 2u]), rc + 4u);
276 bs = radCoef[rc + 1u] * delayed(e, j, i32(radCoef[rc + 3u]), rc + 4u + T);
278 redF[t] = fs; redB[t] = bs;
279 workgroupBarrier();
56d33ddGPU: decimate once per audio sample and spread the listener sum over all elementsJeremy Magland 280 for (var s = ${GATHER_ELEMS / 2}u; s > 0u; s = s >> 1u) {
1e29264cymsimweb: WebGPU companion of cymsim with YAML/form editor, live spectrogram and playbackJeremy Magland 281 if (t < s) { redF[t] += redF[t + s]; redB[t] += redB[t + s]; }
282 workgroupBarrier();
283 }
284 if (t == 0u && j >= 0) {
285 let idx = j - chunk.base;
56d33ddGPU: decimate once per audio sample and spread the listener sum over all elementsJeremy Magland 286 if (idx >= 0 && u32(idx) < chunk.size) { pressure[u32(idx) * GE + wid.x] = redF[0]; pressure[(chunk.size + u32(idx)) * GE + wid.x] = redB[0]; }
13f4dbdexport .cymsim fileJeremy Magland 290// ------------------------------------------------------------------ capture pass (displacement frames)
3534482Linear spectrogram axis by default; capture and replay the displacement movie at a chosen slowdownJeremy Magland 291@group(3) @binding(0) var<storage, read> stateC: array<f32>;
292@group(3) @binding(1) var<storage, read_write> movie: array<f32>;
293@group(3) @binding(2) var<uniform> chunkC: Chunk;
294@group(3) @binding(3) var<storage, read> syncC: array<f32>;
296@compute @workgroup_size(256)
297fn capturePass(@builtin(global_invocation_id) gid: vec3<u32>) {
298 let nd = gid.x;
299 if (nd >= N) { return; }
300 let step = bitcast<u32>(syncC[0]);
301 let frame = (step / M) / chunkC.every - chunkC.frameBase;
302 movie[frame * N + nd] = stateC[5u * nd + 2u];
↑↓ move↵openescclose