GPU: decimate once per audio sample and spread the listener sum over all elements
4 changed files+54−34
THEORY.mdmodified+3−3View file
@@ -44,15 +44,15 @@ The partial sums of two consecutive steps live in two halves of one buffer selec
4444
4545 The Python package scatters every element's acceleration into retarded-time buffers at sub-step resolution and decimates the result to the audio rate. A scatter with per-element delays needs atomics or a per-element history far too long to keep at sub-step resolution (about 300 MB at the default resolution), so the web application reverses the order: the acceleration is first decimated to the audio rate, then delayed and summed.
4646
47-The quantity that is decimated is the same as in Python: the *normal* acceleration of every node, $a_i = n_i \cdot (\ddot u_i, \ddot v_i, \ddot w_i)$ with $n_i$ the unit normal of the rest surface at the node (section 7 of cymsim's model.md), taken as the difference of consecutive *projected* velocities divided by $\Delta t$, which for the leapfrog scheme is the exact central second difference of $q$. The decimation kernel is also the same: a Kaiser-windowed sinc spanning $W = 32$ audio samples ($\beta = 8$, $-6$ dB at $0.46 f_s$, below $-80$ dB from $0.54 f_s$; `decimationKernel` in `src/model/radiation.ts` reproduces `cymsim.decimation.decimation_kernel` to $10^{-16}$). It is applied in the node pass through $W$ accumulator slots per node: at sub-step phase $k$ of audio sample $j_0$ the node adds $a_i$ times the tap $h_w(k)$ to slot $w$, which belongs to audio sample $j_0 + w - W/2 + 1$ and lives at ring position (sample mod $W$). When the sub-step count crosses a sample boundary the oldest slot is complete: its value is written to a per-node history ring (a little longer than the longest delay in audio samples), to the chunk's acceleration output when the record is on, and the slot is cleared for the newest sample. This costs $W$ multiply-adds per node per sub-step and $W N$ floats of state, and it needs no atomics because each node owns its slots. The finished node samples are what the `.cymsim` file stores (section 8).
47+The quantity that is decimated is the same as in Python: the *normal* acceleration of every node, $a_i = n_i \cdot (\ddot u_i, \ddot v_i, \ddot w_i)$ with $n_i$ the unit normal of the rest surface at the node (section 7 of cymsim's model.md), taken as the difference of consecutive *projected* velocities divided by $\Delta t$, which for the leapfrog scheme is the exact central second difference of $q$. The decimation kernel is also the same: a Kaiser-windowed sinc spanning $W = 32$ audio samples ($\beta = 8$, $-6$ dB at $0.46 f_s$, below $-80$ dB from $0.54 f_s$; `decimationKernel` in `src/model/radiation.ts` reproduces `cymsim.decimation.decimation_kernel` to $10^{-16}$). It is applied in the node pass through $W$ accumulator slots per node. The node stores $a_i$ at every sub-step phase $k$ of audio sample $j_0$, and at the last phase, $k = M - 1$, it adds $\sum_k h_w(k)\, a_i(k)$ to slot $w$ for every $w$; slot $w$ belongs to audio sample $j_0 + w - W/2 + 1$ and lives at ring position (sample mod $W$). When the sub-step count crosses a sample boundary the oldest slot is complete: its value is written to a per-node history ring (a little longer than the longest delay in audio samples), to the chunk's acceleration output when the record is on, and the slot is cleared for the newest sample. The arithmetic is the same as adding $a_i h_w(k)$ to the $W$ slots at every sub-step, but the slots are read and written once per audio sample instead of once per sub-step, which on an integrated GPU saved about 10% of the time of a sub-step. It costs $W M$ multiply-adds per node per audio sample and $(W + M) N$ floats of state, and it needs no atomics because each node owns its slots. A chunk ends one sub-step before its last phase is processed, so the stored accelerations of the current sample are part of the exported state; a state exported before they were kept continues exactly with them set to zero, since its slots already hold the earlier phases. The finished node samples are what the `.cymsim` file stores (section 8).
4848
49-The gather is one dispatch of a single workgroup per finished audio sample $j$. Each thread first forms the element means of the node sample from the node history and writes them to a per-element history ring, then reads its elements' histories at the fractional delays of the front path and of the path over the rim with a 16-tap Kaiser-windowed-sinc interpolator (`fractionalDelayTaps`), and the two sums are reduced in shared memory and written to the chunk's output buffer. The interpolator replaces the linear interpolation of the first version, which attenuated the top of the band by a few dB at a half-sample fraction. Against the Python package, whose delays are applied at sub-step resolution before the same decimation, the listener pressure now agrees to $3 \times 10^{-4}$ relative rms on the test configuration (it was a few percent), and the recorded node accelerations to about $10^{-4}$, the int16 quantisation of the record included.
49+The gather is one dispatch per finished audio sample $j$, with one thread per element. Each thread first forms its element's mean of the node sample from the node history and writes it to a per-element history ring, then reads the element's history at the fractional delays of the front path and of the path over the rim with a 16-tap Kaiser-windowed-sinc interpolator (`fractionalDelayTaps`), and the two sums are reduced in shared memory to one partial sum per workgroup, which the host adds up when it reads the chunk. (The first version ran the gather in a single workgroup, which took about 12 microseconds per sub-step on an integrated GPU, averaged over the sub-steps of a sample.) The interpolator replaces the linear interpolation of the first version, which attenuated the top of the band by a few dB at a half-sample fraction. Against the Python package, whose delays are applied at sub-step resolution before the same decimation, the listener pressure now agrees to $3 \times 10^{-4}$ relative rms on the test configuration (it was a few percent), and the recorded node accelerations to about $10^{-4}$, the int16 quantisation of the record included.
5050
5151 The acceleration that radiates is the difference of consecutive *projected* velocities. With the lagged projection, the difference of raw velocities would contain the full rigid acceleration $f_{rigid}/M$ at every step, since each step adds it and the next removes it; a plate that appears to accelerate as a whole radiates like a monopole, which is exactly the component that the front and back terms of an unbaffled plate are meant to cancel. This error was invisible in the nonlinear run and a factor 4.6 in the rms pressure of the linear reference run, where the cancellation is 25:1.
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, 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 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.
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+31−19View file
@@ -2,7 +2,8 @@
22 // contact forces of the hits), the node pass (gather, contact, lagged rigid-mode projection,
33 // leapfrog update, partial sums, exciter update, and the normal acceleration with its decimation
44 // accumulators), the gather pass (element means of a finished audio sample into the history,
5-// then the retarded-time sum at the listener) and the capture pass (displacement frames). Sizes
5+// then the retarded-time sum at the listener, one partial sum per workgroup) and the capture pass
6+// (displacement frames). Sizes
67 // and constants are baked into the source. See THEORY.md for the buffer layouts and the ordering
78 // argument behind the two step counters.
89
@@ -18,6 +19,7 @@ export const HITU_HEADER = 4; // releaseStep ncont pad[2]
1819 export const SYNC_HEADER = 4; // stepA stepB pad[2], then partials[2][G][12]
1920 export const RADCOEF_STRIDE = 36; // coefF coefB i0F i0B tapsF[16] tapsB[16]
2021 export const WORKGROUP_NODES = 256;
22+export const GATHER_ELEMS = 256; // elements per gather workgroup
2123 export const STATE_FLOATS = 13; // per node: q[5] v[5] vPrev[3]
2224 /** Chunk uniform: base size frameBase every accBase accRows pad pad (eight 32-bit words). */
2325 export const CHUNK_WORDS = 8;
@@ -35,7 +37,8 @@ const M: u32 = ${c.M}u; const W: u32 = ${c.W}u; const D: u32 = ${c.D}u; const G:
3537 const DT: f32 = ${f(c.dt)}; const DAMP0: f32 = ${f(1 / (1 + 2 * c.sigma0 * c.dt))}; const NL: f32 = ${f(c.nl)};
3638 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;
3739 const OFF_V: u32 = 5u * N; const OFF_VWP: u32 = 10u * N;
38-const OFF_NHIST: u32 = W * N; const OFF_ACCOUT: u32 = (W + D) * N;
40+const OFF_NHIST: u32 = W * N; const OFF_AHIST: u32 = (W + D) * N; const OFF_ACCOUT: u32 = (W + D + M) * N;
41+const GE: u32 = ${Math.ceil(c.E / GATHER_ELEMS)}u;
3942 const OFF_TAPS: u32 = NF * N; const OFF_PART: u32 = ${SYNC_HEADER}u;
4043
4144 struct Chunk { base: i32, size: u32, frameBase: u32, every: u32, accBase: i32, accRows: u32, pad0: u32, pad1: u32 }
@@ -189,8 +192,18 @@ fn nodePass(@builtin(global_invocation_id) gid: vec3<u32>, @builtin(local_invoca
189192 }
190193 rad[slotE * N + nd] = 0.0;
191194 }
192- let r = modW(i32(j0) + 1 - i32(W / 2u));
193- for (var w = 0u; w < W; w++) { rad[((r + w) % W) * N + nd] += acc * nodeF[OFF_TAPS + w * M + kk]; }
195+ // keep the sub-step's acceleration; once the M values of the sample are in, add them to the W slots at once
196+ rad[OFF_AHIST + kk * N + nd] = acc;
197+ if (kk == M - 1u) {
198+ let r = modW(i32(j0) + 1 - i32(W / 2u));
199+ var sl: array<f32, ${c.W}>;
200+ for (var w = 0u; w < W; w++) { sl[w] = 0.0; }
201+ for (var m = 0u; m < M; m++) {
202+ let a = rad[OFF_AHIST + m * N + nd];
203+ for (var w = 0u; w < W; w++) { sl[w] += a * nodeF[OFF_TAPS + w * M + m]; }
204+ }
205+ for (var w = 0u; w < W; w++) { rad[((r + w) % W) * N + nd] += sl[w]; }
206+ }
194207 }
195208 }
196209 for (var k = 0u; k < 12u; k++) { red[t * 12u + k] = part[k]; }
@@ -222,8 +235,8 @@ fn nodePass(@builtin(global_invocation_id) gid: vec3<u32>, @builtin(local_invoca
222235 @group(2) @binding(5) var<storage, read_write> elemHist: array<f32>;
223236 @group(2) @binding(6) var<storage, read> elemUG: array<u32>;
224237
225-var<workgroup> redF: array<f32, 256>;
226-var<workgroup> redB: array<f32, 256>;
238+var<workgroup> redF: array<f32, ${GATHER_ELEMS}>;
239+var<workgroup> redB: array<f32, ${GATHER_ELEMS}>;
227240
228241 fn delayed(e: u32, j: i32, i0: i32, tapBase: u32) -> f32 {
229242 var s = 0.0;
@@ -234,31 +247,30 @@ fn delayed(e: u32, j: i32, i0: i32, tapBase: u32) -> f32 {
234247 return s;
235248 }
236249
237-@compute @workgroup_size(256)
238-fn gatherPass(@builtin(local_invocation_id) lid: vec3<u32>) {
239- let t = lid.x;
250+@compute @workgroup_size(${GATHER_ELEMS})
251+fn gatherPass(@builtin(local_invocation_id) lid: vec3<u32>, @builtin(workgroup_id) wid: vec3<u32>) {
252+ // one thread per element; each workgroup writes its partial sums, which the host adds up
253+ let t = lid.x; let e = wid.x * ${GATHER_ELEMS}u + t;
240254 let step = bitcast<u32>(syncG[0]); // the node pass of step n has run: the finished sample belongs to s = n - 1
241255 let j = i32((step - 2u) / M) - i32(W / 2u);
242256 var fs = 0.0; var bs = 0.0;
243- if (j >= 0) {
257+ if (j >= 0 && e < E) {
244258 let hrow = OFF_NHIST + (u32(j) % D) * N;
245- for (var e = t; e < E; e += 256u) {
246- let mean = (radR[hrow + elemUG[3u * e]] + radR[hrow + elemUG[3u * e + 1u]] + radR[hrow + elemUG[3u * e + 2u]]) / 3.0;
247- elemHist[(u32(j) % D) * E + e] = mean;
248- let rc = RC * e;
249- fs += radCoef[rc] * delayed(e, j, i32(radCoef[rc + 2u]), rc + 4u);
250- bs += radCoef[rc + 1u] * delayed(e, j, i32(radCoef[rc + 3u]), rc + 4u + T);
251- }
259+ let mean = (radR[hrow + elemUG[3u * e]] + radR[hrow + elemUG[3u * e + 1u]] + radR[hrow + elemUG[3u * e + 2u]]) / 3.0;
260+ elemHist[(u32(j) % D) * E + e] = mean;
261+ let rc = RC * e;
262+ fs = radCoef[rc] * delayed(e, j, i32(radCoef[rc + 2u]), rc + 4u);
263+ bs = radCoef[rc + 1u] * delayed(e, j, i32(radCoef[rc + 3u]), rc + 4u + T);
252264 }
253265 redF[t] = fs; redB[t] = bs;
254266 workgroupBarrier();
255- for (var s = 128u; s > 0u; s = s >> 1u) {
267+ for (var s = ${GATHER_ELEMS / 2}u; s > 0u; s = s >> 1u) {
256268 if (t < s) { redF[t] += redF[t + s]; redB[t] += redB[t + s]; }
257269 workgroupBarrier();
258270 }
259271 if (t == 0u && j >= 0) {
260272 let idx = j - chunk.base;
261- if (idx >= 0 && u32(idx) < chunk.size) { pressure[u32(idx)] = redF[0]; pressure[chunk.size + u32(idx)] = redB[0]; }
273+ 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]; }
262274 }
263275 }
264276
src/gpu/simulator.tsmodified+18−11View file
@@ -6,7 +6,7 @@ import { strainEnergies } from "../model/forces.ts";
66 import { DECIMATION_SPAN, DELAY_TAPS } from "../model/radiation.ts";
77 import type { Setup } from "../model/setup.ts";
88 import type { BackendOptions, BackendState, ChunkResult, Diagnostics, SimulatorBackend } from "../model/simulator.ts";
9-import { CHUNK_WORDS, ELEM_STRIDE, HIT_HEADER, HITU_HEADER, NODE_STRIDE, RADCOEF_STRIDE, STATE_FLOATS, SYNC_HEADER, WORKGROUP_NODES, shaderSource } from "./shaders.ts";
9+import { CHUNK_WORDS, ELEM_STRIDE, GATHER_ELEMS, HIT_HEADER, HITU_HEADER, NODE_STRIDE, RADCOEF_STRIDE, STATE_FLOATS, SYNC_HEADER, WORKGROUP_NODES, shaderSource } from "./shaders.ts";
1010
1111 // evaluated lazily: under Node the WebGPU globals are installed after this module is imported
1212 const storageUsage = () => GPUBufferUsage.STORAGE | GPUBufferUsage.COPY_DST | GPUBufferUsage.COPY_SRC;
@@ -19,7 +19,7 @@ export class GpuSimulator implements SimulatorBackend {
1919 readonly setup: Setup; readonly device: GPUDevice;
2020 samplesDone = 0;
2121 private step = 0;
22- private readonly G: number; private readonly HF: number; private readonly HU: number;
22+ private readonly G: number; private readonly GE: number; private readonly HF: number; private readonly HU: number;
2323 private readonly buffers: Record<string, GPUBuffer> = {};
2424 private readonly pipelines: Record<string, GPUComputePipeline> = {};
2525 private readonly bindGroups: Record<string, GPUBindGroup> = {};
@@ -27,7 +27,7 @@ export class GpuSimulator implements SimulatorBackend {
2727 private readonly framesCapacity: number; private readonly dispCapacity: number; private readonly accRows: number;
2828 private readonly stagingPressure: GPUBuffer; private readonly stagingState: GPUBuffer; private readonly stagingHits: GPUBuffer; private readonly stagingHitsU: GPUBuffer;
2929 private readonly stagingMovie: GPUBuffer | null; private readonly stagingDisp: GPUBuffer | null; private readonly stagingAcc: GPUBuffer | null;
30- private readonly radStateFloats: number; // accumulator slots plus node history, the part of `rad` that is state
30+ private readonly radStateFloats: number; // accumulator slots, node history and the current sample's sub-step accelerations: the part of `rad` that is state
3131
3232 constructor(device: GPUDevice, setup: Setup, chunkCapacity = 8192, opts: BackendOptions | number = {}) {
3333 if (typeof opts === "number") opts = { movieEvery: opts };
@@ -43,10 +43,10 @@ export class GpuSimulator implements SimulatorBackend {
4343 this.accRows = this.recordAcceleration ? chunkCapacity + 2 : 0;
4444 const { mesh, shell: sh, schedule: sc, radiation: rad, substeps: M } = setup;
4545 const N = mesh.n, E = mesh.e, K = sc.K, P = sc.P, W = DECIMATION_SPAN, D = rad.ringLength, T = DELAY_TAPS;
46- this.G = Math.ceil(N / WORKGROUP_NODES);
46+ this.G = Math.ceil(N / WORKGROUP_NODES); this.GE = Math.ceil(E / GATHER_ELEMS);
4747 this.HF = HIT_HEADER + 5 * P; this.HU = HITU_HEADER + 3 * P;
4848 const G = this.G, HF = this.HF, HU = this.HU;
49- this.radStateFloats = (W + D) * N;
49+ this.radStateFloats = (W + D + M) * N;
5050
5151 // ---- pack the set-up
5252 const elemF = new Float32Array(ELEM_STRIDE * E);
@@ -100,11 +100,11 @@ export class GpuSimulator implements SimulatorBackend {
100100 mk("hitsF", hitsF); mk("hitsU", hitsU);
101101 mkEmpty("state", STATE_FLOATS * N); mkEmpty("elemForce", 15 * E);
102102 mkEmpty("sync", SYNC_HEADER + 2 * G * 12); mkEmpty("rad", this.radStateFloats + this.accRows * N); mkEmpty("elemHist", D * E);
103- mkEmpty("pressure", 2 * chunkCapacity);
103+ mkEmpty("pressure", 2 * chunkCapacity * this.GE); // one partial sum per gather workgroup
104104 const uniform = (name: string) => { const b = device.createBuffer({ label: name, size: 4 * CHUNK_WORDS, usage: GPUBufferUsage.UNIFORM | GPUBufferUsage.COPY_DST }); this.buffers[name] = b; return b; };
105105 uniform("chunk"); uniform("chunkD");
106106 const staging = (size: number) => device.createBuffer({ size: Math.max(16, size), usage: GPUBufferUsage.MAP_READ | GPUBufferUsage.COPY_DST });
107- this.stagingPressure = staging(8 * chunkCapacity); this.stagingState = staging(STATE_FLOATS * N * 4);
107+ this.stagingPressure = staging(8 * chunkCapacity * this.GE); this.stagingState = staging(STATE_FLOATS * N * 4);
108108 this.stagingHits = staging(HF * K * 4); this.stagingHitsU = staging(HU * K * 4);
109109 this.stagingMovie = this.framesCapacity > 0 ? (mkEmpty("movie", this.framesCapacity * N), staging(this.framesCapacity * N * 4)) : null;
110110 this.stagingDisp = this.dispCapacity > 0 ? (mkEmpty("disp", this.dispCapacity * N), staging(this.dispCapacity * N * 4)) : null;
@@ -184,18 +184,21 @@ export class GpuSimulator implements SimulatorBackend {
184184 pass.setPipeline(this.pipelines.node); pass.setBindGroup(0, this.bindGroups.empty); pass.setBindGroup(1, this.bindGroups.node); pass.dispatchWorkgroups(wgN);
185185 if (this.sampleDueAt(n) >= 0) {
186186 pass.setPipeline(this.pipelines.gather); pass.setBindGroup(0, this.bindGroups.empty); pass.setBindGroup(1, this.bindGroups.empty);
187- pass.setBindGroup(2, this.bindGroups.gather); pass.dispatchWorkgroups(1);
187+ pass.setBindGroup(2, this.bindGroups.gather); pass.dispatchWorkgroups(this.GE);
188188 }
189189 }
190190 pass.end();
191- enc.copyBufferToBuffer(this.buffers.pressure, 0, this.stagingPressure, 0, 8 * this.chunkCapacity);
191+ const GE = this.GE;
192+ enc.copyBufferToBuffer(this.buffers.pressure, 0, this.stagingPressure, 0, 8 * this.chunkCapacity * GE);
192193 const nFrames = frameSamples.length, nDisp = dispSamples.length;
193194 if (nFrames > 0) enc.copyBufferToBuffer(this.buffers.movie, 0, this.stagingMovie!, 0, nFrames * N * 4);
194195 if (nDisp > 0) enc.copyBufferToBuffer(this.buffers.disp, 0, this.stagingDisp!, 0, nDisp * N * 4);
195196 if (accRows > 0) enc.copyBufferToBuffer(this.buffers.rad, this.radStateFloats * 4, this.stagingAcc!, 0, accRows * N * 4);
196197 device.queue.submit([enc.finish()]);
197198 this.step = n1; this.samplesDone += nSamples;
198- const all = await this.readStaging(this.stagingPressure, 8 * this.chunkCapacity);
199+ const partial = await this.readStaging(this.stagingPressure, 8 * this.chunkCapacity * GE);
200+ const all = new Float32Array(2 * this.chunkCapacity);
201+ for (let i = 0; i < 2 * Math.max(count, 1); i++) { let s = 0; for (let g = 0; g < GE; g++) s += partial[i * GE + g]; all[i] = s; }
199202 const split = (data: Float32Array, k: number) => Array.from({ length: k }, (_, f) => data.slice(f * N, (f + 1) * N));
200203 const frames = nFrames > 0 ? split(await this.readStaging(this.stagingMovie!, nFrames * N * 4), nFrames) : [];
201204 const dispFrames = nDisp > 0 ? split(await this.readStaging(this.stagingDisp!, nDisp * N * 4), nDisp) : [];
@@ -264,8 +267,12 @@ export class GpuSimulator implements SimulatorBackend {
264267 if (st.kind !== "gpu") throw new Error("state was written by another backend");
265268 const { setup } = this, N = setup.mesh.n, E = setup.mesh.e, K = setup.schedule.K, D = setup.radiation.ringLength;
266269 const expect: Record<string, number> = { state: STATE_FLOATS * N, sync: SYNC_HEADER + 2 * this.G * 12, hitsF: this.HF * K, hitsU: this.HU * K, rad: this.radStateFloats, elemHist: D * E };
270+ // states written before the sub-step accelerations were kept per sample lack them; their accumulator
271+ // slots already hold the sample's earlier sub-steps, so zero accelerations continue them exactly
272+ const legacyRad = this.radStateFloats - setup.substeps * N;
267273 for (const [name, len] of Object.entries(expect)) {
268- const a = st.arrays[name];
274+ let a = st.arrays[name];
275+ if (name === "rad" && a?.length === legacyRad) { const full = new Float32Array(len); full.set(a); a = full; }
269276 if (!a || a.length !== len) throw new Error(`state array ${name} is missing or has the wrong size`);
270277 this.device.queue.writeBuffer(this.buffers[name], 0, a.buffer, a.byteOffset, a.byteLength);
271278 }
tests/gpu.check.tsmodified+2−1View file
@@ -54,7 +54,8 @@ async function runBoth(cfg: any, nSamples: number, chunk: number) {
5454 const st = await gpu.exportState();
5555 const gpu2 = new GpuSimulator(device, setup, 4096, opts);
5656 await gpu2.importState(st);
57- const [a, b] = [await gpu.runChunk(50), await gpu2.runChunk(50)];
57+ // long enough for the listener delay, so that state missing from the export would show in the pressure
58+ const [a, b] = [await gpu.runChunk(300), await gpu2.runChunk(300)];
5859 assert.equal(a.firstSample, b.firstSample);
5960 assert.ok(relErr(a.pressureFront, b.pressureFront) < 1e-5, "resumed GPU run diverges from the continued one");
6061 gpu2.dispose(); gpu.dispose();