10729f3WebGPU time loop for the cymbal2 shell model, with numbl setup and output stagesJeremy Magland 1// WGSL megakernel generator for the cymbal2 time loop. One workgroup loops
2// over CHUNK time steps inside a single dispatch, with barriers between the
3// stages of each step; every dispatch is identical, so a whole render is a
4// queue of identical dispatches with no CPU in the loop. All constants
5// (grid sizes, node indices, physics scalars, the stroke schedule thresholds)
6// are baked into the shader source; the operators live in one read-only
7// storage buffer, the evolving state and scratch in a read-write one, the
8// acceleration and force histories in a third.
9//
10// The stage structure mirrors cymbal2_sim.m line for line, with the sparse
11// radial (D1, D2) and azimuthal (T1t, T2) difference matrices applied as
12// stencils and the dense transforms (F, Gi), the per-wavenumber page solves
13// (Fi, Ai, AiB, AiC) and the beam matrices applied as stored matrices.
15export const fmt = (x) => {
16 if (!Number.isFinite(x)) throw new Error('non-finite kernel constant');
17 let s = x.toPrecision(9);
18 if (!/[.e]/.test(s)) s += '.0';
19 return s + 'f';
20};
22// Shared between the megakernel and the multi-dispatch backend: the packed
23// operator buffer, the state-buffer layout, and the baked constants.
24export function buildLayout(setup) {
25 const nrt = setup.nrt | 0;
26 const nth = setup.nth | 0;
27 const nf = setup.nfree | 0;
28 const N = nrt * nth;
29 const K = setup.K | 0;
31 const o = setup.ops;
32 const parts = [
33 ['D1', o.D1, nrt * nrt],
34 ['D2', o.D2, nrt * nrt],
35 ['F', o.F, nth * nth],
36 ['GI', o.Gi, nth * nth],
37 ['RD', o.Rd, nrt],
38 ['R2D', o.R2d, nrt],
39 ['AIB', o.AiB, nrt * nrt * nth],
40 ['AIC', o.AiC, nrt * nrt * nth],
41 ['AI', o.Ai, nrt * nrt * nth],
42 ['FI', o.Fi, nrt * nrt * nth],
43 ['ASIB', o.ASiB, nf * nf],
44 ['ASIC', o.ASiC, nf * nf],
45 ['VM', o.vM, nf],
46 ['VR', o.vR, nf],
47 ['VK', o.vK, nf],
48 ['TIPCOL', o.tipcol, nf],
49 ];
50 const ht = setup.ht;
51 const kgo = setup.tlaunch.map((tl) => {
52 let g = Math.max(0, Math.ceil(tl / ht) - 1);
53 while (g > 0 && g * ht >= tl) g -= 1;
54 while ((g + 1) * ht < tl) g += 1;
55 return g;
56 });
57 parts.push(['KGO', Float32Array.from(kgo), K]);
58 parts.push(['VSTR', Float32Array.from(setup.vstr), K]);
59 // Transposed beam matrices: each beam row thread streams its own row
60 // sequentially instead of striding across columns.
61 const tr = (m) => {
62 const t = new Float32Array(nf * nf);
63 for (let r = 0; r < nf; r++) for (let p = 0; p < nf; p++) t[p + r * nf] = m[r + p * nf];
64 return t;
65 };
66 parts.push(['ASIBT', tr(o.ASiB), nf * nf]);
67 parts.push(['ASICT', tr(o.ASiC), nf * nf]);
69 const off = {};
70 let cursor = 0;
71 for (const [name, arr, len] of parts) {
72 if (arr.length !== len) throw new Error(`ops part ${name}: ${arr.length} != ${len}`);
73 off[name] = cursor;
74 cursor += len;
75 }
76 const opsData = new Float32Array(cursor);
77 for (const [name, arr] of parts) opsData.set(arr, off[name]);
79 const SC = 0;
80 const ST_CVEL = 16;
81 const ST_D = ST_CVEL + K;
82 const ST_U = ST_D + 3 * nf;
83 const ST_UH = ST_U + 3 * N;
84 const scratchNames = [
85 'UR', 'URR', 'V1', 'UT1', 'V2', 'RHSA', 'RHSAF', 'PH3',
86 'PHI', 'PR', 'PRR', 'W1', 'PHT1', 'NL', 'NLF',
87 ];
88 const A = {};
89 let scur = ST_UH + 3 * N;
90 for (const n of scratchNames) {
91 A[n] = scur;
92 scur += N;
93 }
94 const stLen = scur;
96 const c = {
97 HT: ht,
98 HT2I: 1 / (ht * ht),
99 GAP: setup.gap,
100 KM: setup.KM,
101 ALM: setup.alM,
102 PS: setup.ps,
103 AH: setup.aH,
104 AHR: setup.aHr,
105 MISSD: setup.missDepth,
106 EHRS: setup.Eh / setup.Rs,
107 EH2: setup.Eh / 2,
108 RSI: 1 / setup.Rs,
109 SNL: (ht * ht) / setup.rhoh,
110 FSC: setup.fsc,
111 C1R: 1 / (2 * setup.hr),
112 C2R: 1 / (setup.hr * setup.hr),
113 C1T: 1 / (2 * setup.hth),
114 C2T: 1 / (setup.hth * setup.hth),
115 };
117 return {
118 nrt, nth, Nr: setup.Nr | 0, nf, N, K,
119 nsteps: setup.nsteps | 0,
120 opsData, off, stLen,
121 layout: { SC, ST_CVEL, ST_D, ST_U, ST_UH },
122 A, scratchNames, c,
123 LS0: setup.ls - 1, JS0: setup.js - 1,
124 LP0: setup.lp - 1, JP0: setup.jp - 1,
125 TIPI0: setup.tipi - 1,
126 NLR0: setup.nlr0 - 1,
127 HOLDSTEPS: setup.holdSteps | 0,
128 };
129}
131export function buildKernel(setup, opts = {}) {
132 const ws = opts.ws ?? 256;
133 const chunk = opts.chunk ?? 2048;
135 const nrt = setup.nrt | 0;
136 const nth = setup.nth | 0;
137 const Nr = setup.Nr | 0;
138 const nf = setup.nfree | 0;
139 const N = nrt * nth;
140 const K = setup.K | 0;
141 const nsteps = setup.nsteps | 0;
143 // ---- ops buffer packing ----
144 const o = setup.ops;
145 const parts = [
146 ['D1', o.D1, nrt * nrt],
147 ['D2', o.D2, nrt * nrt],
148 ['F', o.F, nth * nth],
149 ['GI', o.Gi, nth * nth],
150 ['RD', o.Rd, nrt],
151 ['R2D', o.R2d, nrt],
152 ['AIB', o.AiB, nrt * nrt * nth],
153 ['AIC', o.AiC, nrt * nrt * nth],
154 ['AI', o.Ai, nrt * nrt * nth],
155 ['FI', o.Fi, nrt * nrt * nth],
156 ['ASIB', o.ASiB, nf * nf],
157 ['ASIC', o.ASiC, nf * nf],
158 ['VM', o.vM, nf],
159 ['VR', o.vR, nf],
160 ['VK', o.vK, nf],
161 ['TIPCOL', o.tipcol, nf],
162 ];
163 // Launch-step thresholds: smallest 0-based kk with (kk+1)*ht >= tlaunch,
164 // decided in f64 exactly as MATLAB's k*ht >= tlaunch decides it.
165 const ht = setup.ht;
166 const kgo = setup.tlaunch.map((tl) => {
167 let g = Math.max(0, Math.ceil(tl / ht) - 1);
168 while (g > 0 && g * ht >= tl) g -= 1;
169 while ((g + 1) * ht < tl) g += 1;
170 return g;
171 });
172 parts.push(['KGO', Float32Array.from(kgo), K]);
173 parts.push(['VSTR', Float32Array.from(setup.vstr), K]);
175 const off = {};
176 let cursor = 0;
177 for (const [name, arr, len] of parts) {
178 if (arr.length !== len) throw new Error(`ops part ${name}: ${arr.length} != ${len}`);
179 off[name] = cursor;
180 cursor += len;
181 }
182 const opsData = new Float32Array(cursor);
183 for (const [name, arr] of parts) opsData.set(arr, off[name]);
185 // ---- state buffer layout ----
186 const SC = 0; // 16 scalars
187 const ST_CVEL = 16;
188 const ST_D = ST_CVEL + K;
189 const ST_U = ST_D + 3 * nf;
190 const ST_UH = ST_U + 3 * N;
191 const scratchNames = [
192 'UR', 'URR', 'V1', 'UT1', 'V2', 'RHSA', 'RHSAF', 'PH3',
193 'PHI', 'PR', 'PRR', 'W1', 'PHT1', 'NL', 'NLF',
194 ];
195 const A = {};
196 let scur = ST_UH + 3 * N;
197 for (const n of scratchNames) {
198 A[n] = scur;
199 scur += N;
200 }
201 const stLen = scur;
203 // ---- baked constants ----
204 const LS0 = setup.ls - 1, JS0 = setup.js - 1;
205 const LP0 = setup.lp - 1, JP0 = setup.jp - 1;
206 const TIPI0 = setup.tipi - 1;
207 const NLR0 = setup.nlr0 - 1; // 0-based first kept NL row
208 const c = {
209 HT: ht,
210 HT2I: 1 / (ht * ht),
211 GAP: setup.gap,
212 KM: setup.KM,
213 ALM: setup.alM,
214 PS: setup.ps,
215 AH: setup.aH,
216 AHR: setup.aHr,
217 MISSD: setup.missDepth,
218 EHRS: setup.Eh / setup.Rs,
219 EH2: setup.Eh / 2,
220 RSI: 1 / setup.Rs,
221 SNL: (ht * ht) / setup.rhoh,
222 FSC: setup.fsc,
223 C1R: 1 / (2 * setup.hr),
224 C2R: 1 / (setup.hr * setup.hr),
225 C1T: 1 / (2 * setup.hth),
226 C2T: 1 / (setup.hth * setup.hth),
227 };
228 const constLines = Object.entries(c)
229 .map(([k, v]) => `const ${k}: f32 = ${fmt(v)};`)
230 .join('\n');
232 const code = /* wgsl */ `
233const NRT = ${nrt}u;
234const NTH = ${nth}u;
235const N = ${N}u;
236const NF = ${nf}u;
237const WS = ${ws}u;
238const CHUNK = ${chunk}u;
239const NSTEPS = ${nsteps}u;
240const KS = ${K}u;
241const NRU = ${Nr}u;
242const LS0 = ${LS0}u; const JS0 = ${JS0}u;
243const LP0 = ${LP0}u; const JP0 = ${JP0}u;
244const TIPI = ${TIPI0}u;
245const NLR0 = ${NLR0}u;
246const HOLDSTEPS = ${setup.holdSteps | 0};
247${constLines}
248${Object.entries(off).map(([k, v]) => `const OP_${k} = ${v}u;`).join('\n')}
249const SC = ${SC}u;
250const ST_CVEL = ${ST_CVEL}u;
251const ST_D = ${ST_D}u;
252const ST_U = ${ST_U}u;
253const ST_UH = ${ST_UH}u;
254${scratchNames.map((n) => `const A_${n} = ${A[n]}u;`).join('\n')}
256@group(0) @binding(0) var<storage, read> ops: array<f32>;
257@group(0) @binding(1) var<storage, read_write> st: array<f32>;
258@group(0) @binding(2) var<storage, read_write> outb: array<f32>;
260var<workgroup> wgKbase: u32;
261var<workgroup> wgF: f32;
262var<workgroup> wgHA: f32;
263var<workgroup> wgHB: f32;
264var<workgroup> wgHC: f32;
265var<workgroup> wgShell: u32;
266var<workgroup> wgRot: u32;
268fn thm(j: u32) -> u32 { return (j + NTH - 1u) % NTH; }
269fn thp(j: u32) -> u32 { return (j + 1u) % NTH; }
271@compute @workgroup_size(${ws})
272fn main(@builtin(local_invocation_index) lid: u32) {
273 var shellA: u32 = 0u; var inC: u32 = 0u; var phase: u32 = 0u;
274 var strokeIdx: u32 = 0u; var holdLeft: i32 = 0; var vtar: f32 = 0.0;
275 var vH: f32 = 0.0; var uHm: f32 = 0.0; var uH0: f32 = 0.0; var uHp: f32 = 0.0;
276 var peakXi: f32 = 0.0; var peakDisp: f32 = 0.0;
277 if (lid == 0u) {
278 wgKbase = u32(st[SC + 0u]);
279 shellA = u32(st[SC + 1u]); inC = u32(st[SC + 2u]); phase = u32(st[SC + 3u]);
280 strokeIdx = u32(st[SC + 4u]); holdLeft = i32(st[SC + 5u]); vtar = st[SC + 6u];
281 vH = st[SC + 7u]; uHm = st[SC + 8u]; uH0 = st[SC + 9u]; uHp = st[SC + 10u];
282 peakXi = st[SC + 12u]; peakDisp = st[SC + 13u];
283 wgRot = u32(st[SC + 11u]);
284 wgShell = shellA;
285 }
286 let kbase = workgroupUniformLoad(&wgKbase);
288 for (var s: u32 = 0u; s < CHUNK; s = s + 1u) {
289 let k = kbase + s;
290 if (k >= NSTEPS) { break; }
291 let dNew = ST_D + (k % 3u) * NF;
292 let dCur = ST_D + ((k + 2u) % 3u) * NF;
293 let dOld = ST_D + ((k + 1u) % 3u) * NF;
295 // ---- stage 1: contact force + hand state machine (lane 0) ----
296 if (lid == 0u) {
297 let xi = st[ST_U + wgRot * N + LS0 + JS0 * NRT] - st[dCur + TIPI] - GAP;
298 var fM: f32 = 0.0;
299 if (xi > 0.0) {
300 fM = KM * pow(xi, ALM);
301 shellA = 1u;
302 if (inC == 0u && strokeIdx < KS && st[ST_CVEL + strokeIdx] == 0.0) {
303 st[ST_CVEL + strokeIdx] = (st[dCur + TIPI] - st[dOld + TIPI]) / HT;
304 }
305 inC = 1u;
306 if (xi > peakXi) { peakXi = xi; }
307 } else {
308 inC = 0u;
309 }
310 outb[NSTEPS + k] = fM;
311 wgF = fM;
312 wgHA = -(uHp - 2.0 * uH0 + uHm);
313 wgHB = -(HT * 0.5) * (uHp - uHm);
314 wgHC = -HT * HT * (PS * uH0 + ((1.0 - PS) * 0.5) * (uHp + uHm));
315 if (phase == 0u) {
316 if (strokeIdx < KS && f32(k) >= ops[OP_KGO + strokeIdx]) {
317 phase = 1u; vtar = ops[OP_VSTR + strokeIdx];
318 }
319 } else if (phase == 1u) {
320 vH = max(vH - AH * HT, -vtar);
321 if (inC != 0u) { phase = 3u; holdLeft = HOLDSTEPS; }
322 else if (vH <= -vtar) { phase = 2u; }
323 } else if (phase == 2u) {
324 if (inC != 0u) { phase = 3u; holdLeft = HOLDSTEPS; }
325 else if (uHp <= -MISSD) { phase = 4u; }
326 } else if (phase == 3u) {
327 holdLeft = holdLeft - 1;
328 if (holdLeft <= 0) { phase = 4u; }
329 } else if (phase == 4u) {
330 vH = min(vH + AHR * HT, 0.5 * vtar);
331 if (vH > 0.0 && uHp >= -vH * vH / (2.0 * AHR)) { phase = 5u; }
332 } else {
333 vH = vH - AHR * HT;
334 if (vH <= 0.0 || uHp >= 0.0) { vH = 0.0; phase = 0u; strokeIdx = strokeIdx + 1u; }
335 }
336 uHm = uH0; uH0 = uHp; uHp = uHp + vH * HT;
337 wgShell = shellA;
338 }
339 workgroupBarrier();
340 storageBarrier();
341 let fM1 = wgF;
342 let hA = wgHA; let hB = wgHB; let hC = wgHC;
344 // ---- stage 2: beam update ----
345 for (var r: u32 = lid; r < NF; r = r + WS) {
346 var sum: f32 = ops[OP_VM + r] * hA + ops[OP_VR + r] * hB + ops[OP_VK + r] * hC
347 + (HT * HT * fM1) * ops[OP_TIPCOL + r];
348 for (var p: u32 = 0u; p < NF; p = p + 1u) {
349 sum = sum + ops[OP_ASIB + r + p * NF] * st[dCur + p]
350 + ops[OP_ASIC + r + p * NF] * st[dOld + p];
351 }
352 st[dNew + r] = sum;
353 }
355 let sh = workgroupUniformLoad(&wgShell);
356 if (sh != 0u) {
357 let rot = workgroupUniformLoad(&wgRot);
358 let uS = ST_U + rot * N;
359 let u1S = ST_U + ((rot + 2u) % 3u) * N;
360 let unS = ST_U + ((rot + 1u) % 3u) * N;
361 let h3S = ST_UH + rot * N;
362 let h13S = ST_UH + ((rot + 2u) % 3u) * N;
363 let h3nS = ST_UH + ((rot + 1u) % 3u) * N;
365 // ---- A: radial/azimuthal derivatives of U ----
366 for (var i: u32 = lid; i < N; i = i + WS) {
367 let l = i % NRT; let j = i / NRT;
368 var ur: f32 = 0.0; var urr: f32 = 0.0;
369 if (l >= 1u && l <= NRT - 2u) {
370 let um = st[uS + (l - 1u) + j * NRT];
371 let u0 = st[uS + l + j * NRT];
372 let up = st[uS + (l + 1u) + j * NRT];
373 ur = (up - um) * C1R;
374 urr = (um - 2.0 * u0 + up) * C2R;
375 }
376 st[A_UR + i] = ur;
377 st[A_URR + i] = urr;
378 let ujm = st[uS + l + thm(j) * NRT];
379 let uj0 = st[uS + l + j * NRT];
380 let ujp = st[uS + l + thp(j) * NRT];
381 st[A_V1 + i] = ops[OP_RD + l] * ur + ops[OP_R2D + l] * ((ujm - 2.0 * uj0 + ujp) * C2T);
382 st[A_UT1 + i] = (ujp - ujm) * C1T;
383 }
384 storageBarrier();
386 // ---- B: bracket pieces and Airy right-hand side ----
387 for (var i: u32 = lid; i < N; i = i + WS) {
388 let l = i % NRT; let j = i / NRT;
389 let urt = (st[A_UR + l + thp(j) * NRT] - st[A_UR + l + thm(j) * NRT]) * C1T;
390 let v2 = ops[OP_RD + l] * urt - ops[OP_R2D + l] * st[A_UT1 + i];
391 st[A_V2 + i] = v2;
392 var rv: f32 = 0.0;
393 if (l >= 3u && l <= NRU - 2u) {
394 let urr = st[A_URR + i]; let v1 = st[A_V1 + i];
395 rv = EHRS * (urr + v1) - EH2 * 2.0 * (urr * v1 - v2 * v2);
396 }
397 st[A_RHSA + i] = rv;
398 }
399 storageBarrier();
401 // ---- C: transform the Airy rhs ----
402 for (var i: u32 = lid; i < N; i = i + WS) {
403 let l = i % NRT; let j = i / NRT;
404 var sum: f32 = 0.0;
405 for (var q: u32 = 0u; q < NTH; q = q + 1u) {
406 sum = sum + st[A_RHSA + l + q * NRT] * ops[OP_F + q + j * NTH];
407 }
408 st[A_RHSAF + i] = sum;
409 }
410 storageBarrier();
412 // ---- D: Airy solve per wavenumber column ----
413 for (var i: u32 = lid; i < N; i = i + WS) {
414 let l = i % NRT; let j = i / NRT;
415 let base = OP_FI + j * NRT * NRT;
416 var sum: f32 = 0.0;
417 for (var p: u32 = 0u; p < NRT; p = p + 1u) {
418 sum = sum + ops[base + l + p * NRT] * st[A_RHSAF + p + j * NRT];
419 }
420 st[A_PH3 + i] = sum;
421 }
422 storageBarrier();
424 // ---- E: stress function back to physical space ----
425 for (var i: u32 = lid; i < N; i = i + WS) {
426 let l = i % NRT; let j = i / NRT;
427 var sum: f32 = 0.0;
428 for (var q: u32 = 0u; q < NTH; q = q + 1u) {
429 sum = sum + st[A_PH3 + l + q * NRT] * ops[OP_GI + q + j * NTH];
430 }
431 st[A_PHI + i] = sum;
432 }
433 storageBarrier();
435 // ---- F: derivatives of the stress function ----
436 for (var i: u32 = lid; i < N; i = i + WS) {
437 let l = i % NRT; let j = i / NRT;
438 var pr: f32 = 0.0; var prr: f32 = 0.0;
439 if (l >= 1u && l <= NRT - 2u) {
440 let pm = st[A_PHI + (l - 1u) + j * NRT];
441 let p0 = st[A_PHI + l + j * NRT];
442 let pp = st[A_PHI + (l + 1u) + j * NRT];
443 pr = (pp - pm) * C1R;
444 prr = (pm - 2.0 * p0 + pp) * C2R;
445 }
446 st[A_PR + i] = pr;
447 st[A_PRR + i] = prr;
448 let pjm = st[A_PHI + l + thm(j) * NRT];
449 let pj0 = st[A_PHI + l + j * NRT];
450 let pjp = st[A_PHI + l + thp(j) * NRT];
451 st[A_W1 + i] = ops[OP_RD + l] * pr + ops[OP_R2D + l] * ((pjm - 2.0 * pj0 + pjp) * C2T);
452 st[A_PHT1 + i] = (pjp - pjm) * C1T;
453 }
454 storageBarrier();
456 // ---- G: von Karman bracket, forcing, membrane term ----
457 for (var i: u32 = lid; i < N; i = i + WS) {
458 let l = i % NRT; let j = i / NRT;
459 var nl: f32 = 0.0;
460 if (l >= NLR0 && l < NRU) {
461 let prt = (st[A_PR + l + thp(j) * NRT] - st[A_PR + l + thm(j) * NRT]) * C1T;
462 let w2 = ops[OP_RD + l] * prt - ops[OP_R2D + l] * st[A_PHT1 + i];
463 let prr = st[A_PRR + i]; let w1 = st[A_W1 + i];
464 nl = prr * st[A_V1 + i] + st[A_URR + i] * w1 - 2.0 * w2 * st[A_V2 + i]
465 - (prr + w1) * RSI;
466 if (l == LS0 && j == JS0) { nl = nl - fM1 * FSC; }
467 nl = nl * SNL;
468 }
469 st[A_NL + i] = nl;
470 }
471 storageBarrier();
473 // ---- H: transform the nonlinear load ----
474 for (var i: u32 = lid; i < N; i = i + WS) {
475 let l = i % NRT; let j = i / NRT;
476 var sum: f32 = 0.0;
477 for (var q: u32 = 0u; q < NTH; q = q + 1u) {
478 sum = sum + st[A_NL + l + q * NRT] * ops[OP_F + q + j * NTH];
479 }
480 st[A_NLF + i] = sum;
481 }
482 storageBarrier();
484 // ---- I: implicit shell update per wavenumber column ----
485 for (var i: u32 = lid; i < N; i = i + WS) {
486 let l = i % NRT; let j = i / NRT;
487 let base = j * NRT * NRT;
488 var sum: f32 = 0.0;
489 for (var p: u32 = 0u; p < NRT; p = p + 1u) {
490 let pj = p + j * NRT;
491 sum = sum + ops[OP_AIB + base + l + p * NRT] * st[h3S + pj]
492 + ops[OP_AIC + base + l + p * NRT] * st[h13S + pj]
493 + ops[OP_AI + base + l + p * NRT] * st[A_NLF + pj];
494 }
495 st[h3nS + i] = sum;
496 }
497 storageBarrier();
499 // ---- J: back to physical space; pickup acceleration ----
500 for (var i: u32 = lid; i < N; i = i + WS) {
501 let l = i % NRT; let j = i / NRT;
502 var sum: f32 = 0.0;
503 for (var q: u32 = 0u; q < NTH; q = q + 1u) {
504 sum = sum + st[h3nS + l + q * NRT] * ops[OP_GI + q + j * NTH];
505 }
506 st[unS + i] = sum;
507 if (l == LP0 && j == JP0) {
508 outb[k] = (sum - 2.0 * st[uS + i] + st[u1S + i]) * HT2I;
509 }
510 }
511 storageBarrier();
512 if (lid == 0u) {
513 wgRot = (rot + 1u) % 3u;
514 if ((k + 1u) % 512u == 0u) {
515 var pd: f32 = 0.0;
516 for (var l2: u32 = 0u; l2 < NRT; l2 = l2 + 1u) {
517 pd = max(pd, abs(st[unS + l2 + JS0 * NRT]));
518 }
519 if (pd > peakDisp) { peakDisp = pd; }
520 }
521 }
522 }
523 workgroupBarrier();
524 storageBarrier();
525 }
527 if (lid == 0u) {
528 st[SC + 0u] = f32(min(kbase + CHUNK, NSTEPS));
529 st[SC + 1u] = f32(shellA); st[SC + 2u] = f32(inC); st[SC + 3u] = f32(phase);
530 st[SC + 4u] = f32(strokeIdx); st[SC + 5u] = f32(holdLeft); st[SC + 6u] = vtar;
531 st[SC + 7u] = vH; st[SC + 8u] = uHm; st[SC + 9u] = uH0; st[SC + 10u] = uHp;
532 st[SC + 11u] = f32(wgRot); st[SC + 12u] = peakXi; st[SC + 13u] = peakDisp;
533 }
534}
535`;
537 return { code, opsData, stLen, N, chunk, ws, off, layout: { SC, ST_CVEL, ST_D, ST_U, ST_UH } };
538}