GPU: field-major set-up data, symmetric bending stiffness, shared rigid-mode sums
3 changed files+52−36
THEORY.mdmodified+1−1View file
@@ -52,7 +52,7 @@ The acceleration that radiates is the difference of consecutive *projected* velo
5252
5353 ## 5. Buffers and precision
5454
55-The GPU keeps the set-up and the state in storage buffers, at most eight in any one pipeline (the default WebGPU limit): per-element data (CST operators, area, slopes, membrane stiffness, damping factor and the 81 entries of $K_e$, 96 floats per element), element node indices, the state (q, v and the previous projected $(u, v, w)$ velocities, 13 floats per node), the element forces, per-node data (inverse mass, mass, the six rigid modes and the nodal normal, 44 floats per node, with the decimation taps at the tail), the adjacency in CSR form, the hits (static patch data and the exciter state, one float and one integer buffer), the synchronization buffer (step counters and partial sums), the radiation buffer (the $W N$ accumulator slots, the node history, the current sample's $M N$ sub-step accelerations, and the chunk's acceleration output rows when recording), the element history, the radiation coefficients (delays and interpolator taps, 36 floats per element) and the chunk's pressure output. A uniform holds the chunk's sample and frame bases. Everything on the GPU is float32. The set-up is float64 and is packed once. Diagnostics (energies, exciter state) are read back between chunks; the plate view reads $w$ from the state buffer directly in its vertex shader without a readback. The whole integrator state (state, sync, hits, radiation and element-history buffers, plus the step count) can be read back and written again, which is how a partial `.cymsim` file is continued.
55+The GPU keeps the set-up and the state in storage buffers, at most eight in any one pipeline (the default WebGPU limit): per-element data (CST operators, area, slopes, membrane stiffness, damping factor and the 45 entries of the upper triangle of $K_e$, 58 floats per element), element node indices, the state (q, v and the previous projected $(u, v, w)$ velocities, 13 floats per node), the element forces, per-node data (inverse mass, mass, the six rigid modes and the nodal normal, 43 floats per node, with the decimation taps at the tail), the adjacency in CSR form, the hits (static patch data and the exciter state, one float and one integer buffer), the synchronization buffer (step counters and partial sums), the radiation buffer (the $W N$ accumulator slots, the node history, the current sample's $M N$ sub-step accelerations, and the chunk's acceleration output rows when recording), the element history, the radiation coefficients (delays and interpolator taps, 36 floats per element) and the chunk's pressure output. The per-element and per-node data are stored field-major (field $f$ of element $e$ at $f E + e$), so that the threads of a workgroup read consecutive words; $K_e$, symmetric to round-off, is symmetrised in float64 before packing, and the product $K_e U_e$ is written out term by term in the shader. Together these halved the time of the element pass on an integrated GPU. A uniform holds the chunk's sample and frame bases. Everything on the GPU is float32. The set-up is float64 and is packed once. Diagnostics (energies, exciter state) are read back between chunks; the plate view reads $w$ from the state buffer directly in its vertex shader without a readback. The whole integrator state (state, sync, hits, radiation and element-history buffers, plus the step count) can be read back and written again, which is how a partial `.cymsim` file is continued.
5656
5757 Float32 is adequate for the leapfrog at these step sizes: the displacement increment per step is about $10^{-3}$ of the displacement, so the update loses a few bits per step in $q$ without a visible effect on the energy balance, and the run is chaotic anyway once max|w|/h exceeds a few units, so no float32 trajectory tracks a float64 one for long. What is checked is the deterministic part (section 6).
5858
src/gpu/shaders.tsmodified+38−25View file
@@ -12,8 +12,10 @@ export interface ShaderConstants {
1212 accRows: number; dt: number; sigma0: number; nl: number;
1313 }
1414
15-export const ELEM_STRIDE = 96; // Gx[3] Gy[3] A z0x z0y C11 C12 C33 cd Ke[81] pad[2]
16-export const NODE_STRIDE = 44; // Minv[5] M[5] Rm[6][5] normal[3] pad[1]; the decimation taps follow the nodes
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.
17+export const ELEM_STRIDE = 58; // fields: Gx[3] Gy[3] A z0x z0y C11 C12 C33 cd Ke[45] (upper triangle, row by row)
18+export const NODE_STRIDE = 43; // fields: Minv[5] M[5] Rm[6][5] normal[3]; the decimation taps follow the nodes
1719 export const HIT_HEADER = 12; // kPoint alpha mass gap velocity zs vs F Fmax pad[3]
1820 export const HITU_HEADER = 4; // releaseStep ncont pad[2]
1921 export const SYNC_HEADER = 4; // stepA stepB pad[2], then partials[2][G][12]
@@ -29,17 +31,31 @@ function f(x: number): string {
2931 return /[.e]/.test(s) ? s : s + ".0";
3032 }
3133
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. */
36+function 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");
46+}
47+
3248 export function shaderSource(c: ShaderConstants): string {
3349 const HF = HIT_HEADER + 5 * c.P, HU = HITU_HEADER + 3 * c.P;
3450 return /* wgsl */ `
3551 const N: u32 = ${c.N}u; const E: u32 = ${c.E}u; const K: u32 = ${c.K}u; const P: u32 = ${c.P}u;
3652 const 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;
3753 const DT: f32 = ${f(c.dt)}; const DAMP0: f32 = ${f(1 / (1 + 2 * c.sigma0 * c.dt))}; const NL: f32 = ${f(c.nl)};
38-const EF: u32 = ${ELEM_STRIDE}u; const NF: u32 = ${NODE_STRIDE}u; const HF: u32 = ${HF}u; const HU: u32 = ${HU}u; const RC: u32 = ${RADCOEF_STRIDE}u;
54+const HF: u32 = ${HF}u; const HU: u32 = ${HU}u; const RC: u32 = ${RADCOEF_STRIDE}u;
3955 const OFF_V: u32 = 5u * N; const OFF_VWP: u32 = 10u * N;
4056 const OFF_NHIST: u32 = W * N; const OFF_AHIST: u32 = (W + D) * N; const OFF_ACCOUT: u32 = (W + D + M) * N;
4157 const GE: u32 = ${Math.ceil(c.E / GATHER_ELEMS)}u;
42-const OFF_TAPS: u32 = NF * N; const OFF_PART: u32 = ${SYNC_HEADER}u;
58+const OFF_TAPS: u32 = ${NODE_STRIDE}u * N; const OFF_PART: u32 = ${SYNC_HEADER}u;
4359
4460 struct Chunk { base: i32, size: u32, frameBase: u32, every: u32, accBase: i32, accRows: u32, pad0: u32, pad1: u32 }
4561
@@ -59,14 +75,13 @@ fn elementPass(@builtin(global_invocation_id) gid: vec3<u32>) {
5975 if (id == 0u) { syncE[1] = syncE[0]; }
6076 if (id < E) {
6177 let e = id;
62- let b = e * EF;
6378 let n0 = elemU[3u * e]; let n1 = elemU[3u * e + 1u]; let n2 = elemU[3u * e + 2u];
6479 var qe: array<f32, 15>;
6580 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]; }
66- let gx = vec3<f32>(elemF[b], elemF[b + 1u], elemF[b + 2u]);
67- let gy = vec3<f32>(elemF[b + 3u], elemF[b + 4u], elemF[b + 5u]);
68- let A = elemF[b + 6u]; let z0x = elemF[b + 7u]; let z0y = elemF[b + 8u];
69- let C11 = elemF[b + 9u]; let C12 = elemF[b + 10u]; let C33 = elemF[b + 11u]; let cd = elemF[b + 12u];
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];
7085 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]);
7186 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);
7287 let ax = NL * wx + z0x; let ay = NL * wy + z0y;
@@ -86,12 +101,7 @@ fn elementPass(@builtin(global_invocation_id) gid: vec3<u32>) {
86101 }
87102 var U: array<f32, 9>;
88103 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]; }
89- for (var i = 0u; i < 9u; i++) {
90- var s = 0.0;
91- for (var j = 0u; j < 9u; j++) { s += elemF[b + 13u + 9u * i + j] * U[j]; }
92- let slot = 5u * (i / 3u) + 2u + (i % 3u);
93- out[slot] -= s;
94- }
104+${bendingProduct()}
95105 for (var i = 0u; i < 15u; i++) { elemForce[15u * e + i] = out[i]; }
96106 } else if (id < E + K) {
97107 // contact forces of hit k at this step
@@ -130,6 +140,7 @@ fn elementPass(@builtin(global_invocation_id) gid: vec3<u32>) {
130140 @group(1) @binding(8) var<uniform> chunkN: Chunk;
131141
132142 var<workgroup> red: array<f32, ${12 * WORKGROUP_NODES}>;
143+var<workgroup> csum: array<f32, 12>;
133144
134145 fn modW(x: i32) -> u32 { return u32(((x % i32(W)) + i32(W)) % i32(W)); }
135146
@@ -138,12 +149,15 @@ fn nodePass(@builtin(global_invocation_id) gid: vec3<u32>, @builtin(local_invoca
138149 let nd = gid.x; let t = lid.x;
139150 let step = bitcast<u32>(sync[1]);
140151 let parity = step & 1u; let prevp = (step + 1u) & 1u;
141- var c: array<f32, 12>;
142- for (var k = 0u; k < 12u; k++) {
152+ // twelve threads sum the previous step's partials of all workgroups; the others read the sums from shared memory
153+ if (t < 12u) {
143154 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;
155+ for (var g = 0u; g < G; g++) { s += sync[OFF_PART + (prevp * G + g) * 12u + t]; }
156+ csum[t] = s;
146157 }
158+ workgroupBarrier();
159+ var c: array<f32, 12>;
160+ for (var k = 0u; k < 12u; k++) { c[k] = csum[k]; }
147161 var part: array<f32, 12>;
148162 for (var k = 0u; k < 12u; k++) { part[k] = 0.0; }
149163 if (nd < N) {
@@ -163,19 +177,18 @@ fn nodePass(@builtin(global_invocation_id) gid: vec3<u32>, @builtin(local_invoca
163177 }
164178 }
165179 }
166- let nb = nd * NF;
167180 var acc = 0.0;
168181 for (var d = 0u; d < 5u; d++) {
169182 let idx = 5u * nd + d;
170183 var vd = state[OFF_V + idx]; var qd = state[idx];
171- for (var k = 0u; k < 6u; k++) { let r = nodeF[nb + 10u + 5u * k + d]; vd -= c[k] * r; qd -= c[6u + k] * r; }
172- let v1 = (vd + DT * fo[d] * nodeF[nb + d]) * DAMP0;
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;
173186 // normal acceleration of the previous step from consecutive projected velocities
174- if (d < 3u) { acc += nodeF[nb + 40u + d] * (vd - state[OFF_VWP + 3u * nd + d]) / DT; state[OFF_VWP + 3u * nd + d] = vd; }
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; }
175188 let q1 = qd + DT * v1;
176189 state[idx] = q1; state[OFF_V + idx] = v1;
177- let mm = nodeF[nb + 5u + d];
178- 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; }
190+ let mm = nodeF[(5u + d) * N + nd];
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; }
179192 }
180193 // decimation of the acceleration of sub-step s = step - 1 into the W accumulator slots of this node;
181194 // slot w belongs to audio sample j0 + w - W/2 + 1, kept at ring position (sample mod W)
src/gpu/simulator.tsmodified+13−10View file
@@ -49,22 +49,25 @@ export class GpuSimulator implements SimulatorBackend {
4949 this.radStateFloats = (W + D + M) * N;
5050
5151 // ---- pack the set-up
52+ // field-major: field f of element e at f * E + e, of node n at f * N + n
5253 const elemF = new Float32Array(ELEM_STRIDE * E);
5354 for (let e = 0; e < E; e++) {
54- const b = ELEM_STRIDE * e;
55- for (let i = 0; i < 3; i++) { elemF[b + i] = mesh.Gx[3 * e + i]; elemF[b + 3 + i] = mesh.Gy[3 * e + i]; }
56- elemF[b + 6] = mesh.area[e]; elemF[b + 7] = sh.z0x[e]; elemF[b + 8] = sh.z0y[e];
57- elemF[b + 9] = sh.Cm[3 * e]; elemF[b + 10] = sh.Cm[3 * e + 1]; elemF[b + 11] = sh.Cm[3 * e + 2];
58- elemF[b + 12] = 2 * setup.config.material.sigma1 * sh.rhoHA[e];
59- for (let i = 0; i < 81; i++) elemF[b + 13 + i] = sh.Ke[81 * e + i];
55+ const set = (f: number, x: number) => { elemF[f * E + e] = x; };
56+ for (let i = 0; i < 3; i++) { set(i, mesh.Gx[3 * e + i]); set(3 + i, mesh.Gy[3 * e + i]); }
57+ set(6, mesh.area[e]); set(7, sh.z0x[e]); set(8, sh.z0y[e]);
58+ set(9, sh.Cm[3 * e]); set(10, sh.Cm[3 * e + 1]); set(11, sh.Cm[3 * e + 2]);
59+ set(12, 2 * setup.config.material.sigma1 * sh.rhoHA[e]);
60+ // upper triangle of K_e, symmetrised in float64 (it is symmetric to round-off)
61+ const K = sh.Ke.subarray(81 * e, 81 * e + 81);
62+ for (let i = 0, f = 13; i < 9; i++) for (let j = i; j < 9; j++) set(f++, 0.5 * (K[9 * i + j] + K[9 * j + i]));
6063 }
6164 const elemU = new Uint32Array(3 * E); elemU.set(mesh.tri);
6265 const nodeF = new Float32Array(NODE_STRIDE * N + W * M);
6366 for (let n = 0; n < N; n++) {
64- const b = NODE_STRIDE * n;
65- for (let d = 0; d < 5; d++) { nodeF[b + d] = sh.Minv[5 * n + d]; nodeF[b + 5 + d] = sh.M[5 * n + d]; }
66- for (let k = 0; k < 6; k++) for (let d = 0; d < 5; d++) nodeF[b + 10 + 5 * k + d] = sh.Rm[(k * N + n) * 5 + d];
67- for (let d = 0; d < 3; d++) nodeF[b + 40 + d] = sh.normalNode[3 * n + d];
67+ const set = (f: number, x: number) => { nodeF[f * N + n] = x; };
68+ for (let d = 0; d < 5; d++) { set(d, sh.Minv[5 * n + d]); set(5 + d, sh.M[5 * n + d]); }
69+ for (let k = 0; k < 6; k++) for (let d = 0; d < 5; d++) set(10 + 5 * k + d, sh.Rm[(k * N + n) * 5 + d]);
70+ for (let d = 0; d < 3; d++) set(40 + d, sh.normalNode[3 * n + d]);
6871 }
6972 nodeF.set(rad.taps, NODE_STRIDE * N);
7073 const adj = new Uint32Array(N + 1 + sh.adjEntries.length); adj.set(sh.adjOffsets); adj.set(sh.adjEntries, N + 1);