1e29264cymsimweb: WebGPU companion of cymsim with YAML/form editor, live spectrogram and playbackJeremy Magland 1// WGSL kernels of the time loop. Three entry points: the element pass (forces, sigma1 term,
2// contact forces of the hits, decimation accumulators of the radiated acceleration), the node
3// pass (gather, contact, lagged rigid-mode projection, leapfrog update, partial sums, exciter
4// update) and the gather pass (retarded-time sum at the listener once per audio sample). Sizes
5// and constants are baked into the source. See THEORY.md for the buffer layouts and the ordering
6// argument behind the two step counters.
8export interface ShaderConstants {
9 N: number; E: number; K: number; P: number; M: number; W: number; D: number; G: number;
10 dt: number; sigma0: number; nl: number;
11}
13export const ELEM_STRIDE = 96; // Gx[3] Gy[3] A z0x z0y C11 C12 C33 cd Ke[81] pad[2]
14export const NODE_STRIDE = 40; // Minv[5] M[5] Rm[6][5]
15export const HIT_HEADER = 12; // kPoint alpha mass gap velocity zs vs F Fmax pad[3]
16export const HITU_HEADER = 4; // releaseStep ncont pad[2]
17export const SYNC_HEADER = 4; // stepA stepB pad[2], then partials[2][G][12]
18export const WORKGROUP_NODES = 256;
20function f(x: number): string {
21 const s = x.toPrecision(9);
22 return /[.e]/.test(s) ? s : s + ".0";
23}
25export function shaderSource(c: ShaderConstants): string {
26 const HF = HIT_HEADER + 5 * c.P, HU = HITU_HEADER + 3 * c.P;
27 return /* wgsl */ `
28const N: u32 = ${c.N}u; const E: u32 = ${c.E}u; const K: u32 = ${c.K}u; const P: u32 = ${c.P}u;
29const M: u32 = ${c.M}u; const W: u32 = ${c.W}u; const D: u32 = ${c.D}u; const G: u32 = ${c.G}u;
30const DT: f32 = ${f(c.dt)}; const DAMP0: f32 = ${f(1 / (1 + 2 * c.sigma0 * c.dt))}; const NL: f32 = ${f(c.nl)};
31const EF: u32 = ${ELEM_STRIDE}u; const NF: u32 = ${NODE_STRIDE}u; const HF: u32 = ${HF}u; const HU: u32 = ${HU}u;
32const OFF_V: u32 = 5u * N; const OFF_ACC: u32 = 10u * N; const OFF_VWP: u32 = 11u * N;
33const OFF_HIST: u32 = E * W; const OFF_TAPS: u32 = EF * E; const OFF_PART: u32 = ${SYNC_HEADER}u;
35struct Chunk { base: i32, size: u32 }
37// ------------------------------------------------------------------ element pass
38@group(0) @binding(0) var<storage, read> elemF: array<f32>;
39@group(0) @binding(1) var<storage, read> elemU: array<u32>;
40@group(0) @binding(2) var<storage, read> stateR: array<f32>;
41@group(0) @binding(3) var<storage, read_write> elemForce: array<f32>;
42@group(0) @binding(4) var<storage, read_write> hitsF: array<f32>;
43@group(0) @binding(5) var<storage, read> hitsUR: array<u32>;
44@group(0) @binding(6) var<storage, read_write> syncE: array<f32>;
45@group(0) @binding(7) var<storage, read_write> rad: array<f32>;
47@compute @workgroup_size(64)
48fn elementPass(@builtin(global_invocation_id) gid: vec3<u32>) {
49 let id = gid.x;
50 let step = bitcast<u32>(syncE[0]);
51 if (id == 0u) { syncE[1] = syncE[0]; }
52 if (id < E) {
53 let e = id;
54 let b = e * EF;
55 let n0 = elemU[3u * e]; let n1 = elemU[3u * e + 1u]; let n2 = elemU[3u * e + 2u];
56 var qe: array<f32, 15>;
57 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]; }
58 let gx = vec3<f32>(elemF[b], elemF[b + 1u], elemF[b + 2u]);
59 let gy = vec3<f32>(elemF[b + 3u], elemF[b + 4u], elemF[b + 5u]);
60 let A = elemF[b + 6u]; let z0x = elemF[b + 7u]; let z0y = elemF[b + 8u];
61 let C11 = elemF[b + 9u]; let C12 = elemF[b + 10u]; let C33 = elemF[b + 11u]; let cd = elemF[b + 12u];
62 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]);
63 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);
64 let ax = NL * wx + z0x; let ay = NL * wy + z0y;
65 let exx = ux + NL * 0.5 * wx * wx + z0x * wx;
66 let eyy = vy + NL * 0.5 * wy * wy + z0y * wy;
67 let gxy = uy + vx + NL * wx * wy + z0x * wy + z0y * wx;
68 let Nxx = C11 * exx + C12 * eyy; let Nyy = C12 * exx + C11 * eyy; let Nxy = C33 * gxy;
69 let vw = vec3<f32>(stateR[OFF_V + 5u * n0 + 2u], stateR[OFF_V + 5u * n1 + 2u], stateR[OFF_V + 5u * n2 + 2u]);
70 let vwx = dot(gx, vw); let vwy = dot(gy, vw);
71 var out: array<f32, 15>;
72 for (var i = 0u; i < 3u; i++) {
73 let g_x = gx[i]; let g_y = gy[i];
74 out[5u * i] = -A * (Nxx * g_x + Nxy * g_y);
75 out[5u * i + 1u] = -A * (Nyy * g_y + Nxy * g_x);
76 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);
77 out[5u * i + 3u] = 0.0; out[5u * i + 4u] = 0.0;
78 }
79 var U: array<f32, 9>;
80 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]; }
81 for (var i = 0u; i < 9u; i++) {
82 var s = 0.0;
83 for (var j = 0u; j < 9u; j++) { s += elemF[b + 13u + 9u * i + j] * U[j]; }
84 let slot = 5u * (i / 3u) + 2u + (i % 3u);
85 out[slot] -= s;
86 }
87 for (var i = 0u; i < 15u; i++) { elemForce[15u * e + i] = out[i]; }
88 // decimation accumulators of the radiated acceleration (of the previous step)
89 if (step >= 1u) {
90 let s = step - 1u; let j0 = s / M; let kk = s - j0 * M;
91 if (kk == 0u && j0 > 0u) {
92 let j = i32(j0) - i32(W / 2u);
93 if (j >= 0) { rad[OFF_HIST + e * D + (u32(j) % D)] = rad[e * W]; }
94 for (var w = 0u; w < W - 1u; w++) { rad[e * W + w] = rad[e * W + w + 1u]; }
95 rad[e * W + W - 1u] = 0.0;
96 }
97 let a = (stateR[OFF_ACC + n0] + stateR[OFF_ACC + n1] + stateR[OFF_ACC + n2]) / 3.0;
98 for (var w = 0u; w < W; w++) { rad[e * W + w] += a * elemF[OFF_TAPS + w * M + kk]; }
99 }
100 } else if (id < E + K) {
101 // contact forces of hit k at this step
102 let k = id - E; let hf = k * HF; let hu = k * HU;
103 var wp: array<f32, ${c.P}>;
104 var wmax = -1.0e30;
105 for (var p = 0u; p < P; p++) {
106 var w = hitsF[hf + ${HIT_HEADER}u + 3u * P + p];
107 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]; }
108 wp[p] = w; wmax = max(wmax, w);
109 }
110 var zs = hitsF[hf + 5u];
111 if (step == hitsUR[hu]) { zs = wmax + hitsF[hf + 3u]; hitsF[hf + 5u] = zs; hitsF[hf + 6u] = -hitsF[hf + 4u]; }
112 let isActive = step >= hitsUR[hu];
113 let kPoint = hitsF[hf]; let alpha = hitsF[hf + 1u];
114 var F = 0.0;
115 for (var p = 0u; p < P; p++) {
116 let delta = max(wp[p] - zs, 0.0);
117 var fp = 0.0;
118 if (isActive && delta > 0.0) { fp = kPoint * pow(delta, alpha); }
119 hitsF[hf + ${HIT_HEADER}u + 4u * P + p] = fp; F += fp;
120 }
121 hitsF[hf + 7u] = F;
122 }
123}
125// ------------------------------------------------------------------ node pass
126@group(1) @binding(0) var<storage, read_write> state: array<f32>;
127@group(1) @binding(1) var<storage, read> elemForceR: array<f32>;
128@group(1) @binding(2) var<storage, read> nodeF: array<f32>;
129@group(1) @binding(3) var<storage, read> adj: array<u32>;
130@group(1) @binding(4) var<storage, read_write> hitsFN: array<f32>;
131@group(1) @binding(5) var<storage, read_write> hitsU: array<u32>;
132@group(1) @binding(6) var<storage, read_write> sync: array<f32>;
134var<workgroup> red: array<f32, ${12 * WORKGROUP_NODES}>;
136@compute @workgroup_size(${WORKGROUP_NODES})
137fn nodePass(@builtin(global_invocation_id) gid: vec3<u32>, @builtin(local_invocation_id) lid: vec3<u32>, @builtin(workgroup_id) wid: vec3<u32>) {
138 let nd = gid.x; let t = lid.x;
139 let step = bitcast<u32>(sync[1]);
140 let parity = step & 1u; let prevp = (step + 1u) & 1u;
141 var c: array<f32, 12>;
142 for (var k = 0u; k < 12u; k++) {
143 var s = 0.0;
144 for (var g = 0u; g < G; g++) { s += sync[OFF_PART + (prevp * G + g) * 12u + k]; }
145 c[k] = s;
146 }
147 var part: array<f32, 12>;
148 for (var k = 0u; k < 12u; k++) { part[k] = 0.0; }
149 if (nd < N) {
150 var fo: array<f32, 5>;
151 for (var d = 0u; d < 5u; d++) { fo[d] = 0.0; }
152 let a0 = adj[nd]; let a1 = adj[nd + 1u];
153 for (var a = a0; a < a1; a++) {
154 let ent = adj[N + 1u + a]; let e = ent / 3u; let i = ent - 3u * e;
155 for (var d = 0u; d < 5u; d++) { fo[d] += elemForceR[15u * e + 5u * i + d]; }
156 }
157 for (var k = 0u; k < K; k++) {
158 for (var p = 0u; p < P; p++) {
159 for (var i = 0u; i < 3u; i++) {
160 if (hitsU[k * HU + ${HITU_HEADER}u + 3u * p + i] == nd) {
161 fo[2] -= hitsFN[k * HF + ${HIT_HEADER}u + 4u * P + p] * hitsFN[k * HF + ${HIT_HEADER}u + 3u * p + i];
162 }
163 }
164 }
165 }
166 let nb = nd * NF;
167 for (var d = 0u; d < 5u; d++) {
168 let idx = 5u * nd + d;
169 var vd = state[OFF_V + idx]; var qd = state[idx];
170 for (var k = 0u; k < 6u; k++) { let r = nodeF[nb + 10u + 5u * k + d]; vd -= c[k] * r; qd -= c[6u + k] * r; }
171 let v1 = (vd + DT * fo[d] * nodeF[nb + d]) * DAMP0;
172 if (d == 2u) { state[OFF_ACC + nd] = (vd - state[OFF_VWP + nd]) / DT; state[OFF_VWP + nd] = vd; }
173 let q1 = qd + DT * v1;
174 state[idx] = q1; state[OFF_V + idx] = v1;
175 let mm = nodeF[nb + 5u + d];
176 for (var k = 0u; k < 6u; k++) { let r = nodeF[nb + 10u + 5u * k + d] * mm; part[k] += r * v1; part[6u + k] += r * q1; }
177 }
178 }
179 for (var k = 0u; k < 12u; k++) { red[t * 12u + k] = part[k]; }
180 workgroupBarrier();
181 for (var s = ${WORKGROUP_NODES / 2}u; s > 0u; s = s >> 1u) {
182 if (t < s) { for (var k = 0u; k < 12u; k++) { red[t * 12u + k] += red[(t + s) * 12u + k]; } }
183 workgroupBarrier();
184 }
185 if (t == 0u) { for (var k = 0u; k < 12u; k++) { sync[OFF_PART + (parity * G + wid.x) * 12u + k] = red[k]; } }
186 if (nd == 0u) {
187 for (var k = 0u; k < K; k++) {
188 let hf = k * HF; let hu = k * HU;
189 let F = hitsFN[hf + 7u];
190 let vs = hitsFN[hf + 6u] + DT * F / hitsFN[hf + 2u];
191 hitsFN[hf + 6u] = vs; hitsFN[hf + 5u] = hitsFN[hf + 5u] + DT * vs;
192 hitsFN[hf + 8u] = max(hitsFN[hf + 8u], F);
193 if (F > 0.0) { hitsU[hu + 1u] = hitsU[hu + 1u] + 1u; }
194 }
195 sync[0] = bitcast<f32>(step + 1u);
196 }
197}
199// ------------------------------------------------------------------ gather pass (listener)
200@group(2) @binding(0) var<storage, read> syncG: array<f32>;
201@group(2) @binding(1) var<storage, read> radR: array<f32>;
202@group(2) @binding(2) var<storage, read> radCoef: array<f32>;
203@group(2) @binding(3) var<storage, read_write> pressure: array<f32>;
204@group(2) @binding(4) var<uniform> chunk: Chunk;
206var<workgroup> redF: array<f32, 256>;
207var<workgroup> redB: array<f32, 256>;
209@compute @workgroup_size(256)
210fn gatherPass(@builtin(local_invocation_id) lid: vec3<u32>) {
211 let t = lid.x;
212 let step = bitcast<u32>(syncG[0]);
213 let j = i32((step - 1u) / M) - i32(W / 2u);
214 var fs = 0.0; var bs = 0.0;
215 if (j >= 0) {
216 for (var e = t; e < E; e += 256u) {
217 let rc = 6u * e;
218 let jf = j - i32(radCoef[rc + 2u]); let jb = j - i32(radCoef[rc + 4u]);
219 if (jf >= 1) {
220 let fr = radCoef[rc + 3u];
221 fs += radCoef[rc] * (radR[OFF_HIST + e * D + (u32(jf) % D)] * (1.0 - fr) + radR[OFF_HIST + e * D + (u32(jf - 1) % D)] * fr);
222 }
223 if (jb >= 1) {
224 let fr = radCoef[rc + 5u];
225 bs += radCoef[rc + 1u] * (radR[OFF_HIST + e * D + (u32(jb) % D)] * (1.0 - fr) + radR[OFF_HIST + e * D + (u32(jb - 1) % D)] * fr);
226 }
227 }
228 }
229 redF[t] = fs; redB[t] = bs;
230 workgroupBarrier();
231 for (var s = 128u; s > 0u; s = s >> 1u) {
232 if (t < s) { redF[t] += redF[t + s]; redB[t] += redB[t + s]; }
233 workgroupBarrier();
234 }
235 if (t == 0u && j >= 0) {
236 let idx = j - chunk.base;
237 if (idx >= 0 && u32(idx) < chunk.size) { pressure[u32(idx)] = redF[0]; pressure[chunk.size + u32(idx)] = redB[0]; }
238 }
239}
240`;
241}