d8f9e54Clean rebuild of the steel sheet simulatorJeremy Magland 1// Rigid-body projection: (1) mass-weighted sums of the six rigid modes for
2// displacement and velocity, kinetic energy and max |w| in one workgroup;
a2a1b68Drumstick by default, hit sequences, 170 GPa steelJeremy Magland 3// (2) subtraction of the rigid components from every node. The in-plane
4// inertia may be scaled by P.inPlaneMass; the rigid coefficients are ratios
5// and do not see the factor, the kinetic energy does.
d8f9e54Clean rebuild of the steel sheet simulatorJeremy Magland 6
7@group(0) @binding(0) var<uniform> P: Params;
8@group(0) @binding(1) var<storage, read_write> pos: array<f32>; // [u | v | w | wLo]
9@group(0) @binding(2) var<storage, read_write> vel: array<f32>; // [du | dv | dw]
10@group(0) @binding(3) var<storage, read_write> stats: array<f32>; // 0..11 coefficients, 12 KE, 13 max|w|
12const NS = 14u;
13var<workgroup> sh: array<array<f32, NS>, 256>;
15@compute @workgroup_size(256)
16fn sums(@builtin(local_invocation_id) lid: vec3u) {
17 let t = lid.x;
18 var s: array<f32, NS>;
19 for (var q = 0u; q < NS; q++) { s[q] = 0.0; }
20 for (var k = t; k < P.n; k += 256u) {
21 let i = i32(k % P.nx); let j = i32(k / P.nx);
22 let x = nodeX(i); let y = nodeY(j);
23 let m = P.rhoh * P.hx * P.hy * nodeWeight(i, j);
24 let n = P.n;
25 let p = vec3f(pos[k], pos[n + k], pos[2u * n + k]);
26 let v = vec3f(vel[k], vel[n + k], vel[2u * n + k]);
27 s[0] += m * p.x; s[1] += m * p.y; s[2] += m * p.z;
28 s[3] += m * (x * p.y - y * p.x); s[4] += m * y * p.z; s[5] += m * (-x * p.z);
29 s[6] += m * v.x; s[7] += m * v.y; s[8] += m * v.z;
30 s[9] += m * (x * v.y - y * v.x); s[10] += m * y * v.z; s[11] += m * (-x * v.z);
a2a1b68Drumstick by default, hit sequences, 170 GPa steelJeremy Magland 31 s[12] += 0.5 * m * (P.inPlaneMass * (v.x * v.x + v.y * v.y) + v.z * v.z);
33 }
34 for (var q = 0u; q < NS; q++) { sh[t][q] = s[q]; }
35 workgroupBarrier();
36 for (var stride = 128u; stride > 0u; stride >>= 1u) {
37 if (t < stride) {
38 for (var q = 0u; q < 13u; q++) { sh[t][q] += sh[t + stride][q]; }
39 sh[t][13] = max(sh[t][13], sh[t + stride][13]);
40 }
41 workgroupBarrier();
42 }
43 if (t == 0u) {
44 let inv = array<f32, 6>(1.0 / P.mass, 1.0 / P.mass, 1.0 / P.mass, 1.0 / P.Iz, 1.0 / P.Ix, 1.0 / P.Iy);
45 for (var q = 0u; q < 6u; q++) {
46 stats[q] = sh[0][q] * inv[q];
47 stats[q + 6u] = sh[0][q + 6u] * inv[q];
48 }
49 // Kinetic energy of the non-rigid motion (the six rigid modes are mass-orthogonal).
a2a1b68Drumstick by default, hit sequences, 170 GPa steelJeremy Magland 50 let ip = array<f32, 6>(P.inPlaneMass, P.inPlaneMass, 1.0, P.inPlaneMass, 1.0, 1.0);
a2a1b68Drumstick by default, hit sequences, 170 GPa steelJeremy Magland 52 for (var q = 6u; q < 12u; q++) { keRigid += 0.5 * ip[q - 6u] * sh[0][q] * sh[0][q] * inv[q - 6u]; }
d8f9e54Clean rebuild of the steel sheet simulatorJeremy Magland 53 stats[12] = max(sh[0][12] - keRigid, 0.0);
54 stats[13] = sh[0][13];
55 }
56}
58@compute @workgroup_size(256)
59fn apply(@builtin(global_invocation_id) gid: vec3u) {
60 let k = gid.x;
61 if (k >= P.n) { return; }
62 let i = i32(k % P.nx); let j = i32(k / P.nx);
63 if (nodeWeight(i, j) == 0.0) { return; } // outside the sheet
64 let x = nodeX(i); let y = nodeY(j);
65 let n = P.n;
66 pos[k] -= stats[0] - stats[3] * y;
67 pos[n + k] -= stats[1] + stats[3] * x;
68 pos[2u * n + k] -= stats[2] + stats[4] * y - stats[5] * x;
69 vel[k] -= stats[6] - stats[9] * y;
70 vel[n + k] -= stats[7] + stats[9] * x;
71 vel[2u * n + k] -= stats[8] + stats[10] * y - stats[11] * x;
72}
74// Removes the three rigid components (translation and two rotations) from a
75// per-node acceleration field so that it is self-equilibrated. One workgroup.
76@group(0) @binding(4) var<storage, read_write> field: array<f32>;
77var<workgroup> sh3: array<vec3f, 256>;
79@compute @workgroup_size(256)
80fn projectField(@builtin(local_invocation_id) lid: vec3u) {
81 let t = lid.x;
82 var s = vec3f(0.0);
83 for (var k = t; k < P.n; k += 256u) {
84 let i = i32(k % P.nx); let j = i32(k / P.nx);
85 let m = P.rhoh * P.hx * P.hy * nodeWeight(i, j);
86 let f = field[k];
87 s += vec3f(m * f, m * nodeY(j) * f, -m * nodeX(i) * f);
88 }
89 sh3[t] = s;
90 workgroupBarrier();
91 for (var stride = 128u; stride > 0u; stride >>= 1u) {
92 if (t < stride) { sh3[t] += sh3[t + stride]; }
93 workgroupBarrier();
94 }
95 let c = vec3f(sh3[0].x / P.mass, sh3[0].y / P.Ix, sh3[0].z / P.Iy);
96 for (var k = t; k < P.n; k += 256u) {
97 let i = i32(k % P.nx); let j = i32(k / P.nx);
98 if (nodeWeight(i, j) > 0.0) { field[k] -= c.x + c.y * nodeY(j) - c.z * nodeX(i); }
99 }
100}