2 * The external operations, as GPU dispatches — instantiated per object.
3 *
4 * Everything a compiled model does is element-wise, so these ops are the only
5 * places where a point reads anything but itself: the difference stencils on
6 * each object's own grid, the wall-masked Laplacian in the air, and the
7 * couplings that carry an object's motion into the air. Where dulcimer had a
8 * fixed set of five ops over exactly two grids, here the scene builder
9 * instantiates one op per (operation, object) pair — `dxx_s1` is the second
10 * difference on string s1's grid, `rad_m1` samples membrane m1 at each air
11 * point — each with its own geometry baked into its WGSL. The generated
12 * model .m calls them by these names.
13 *
14 * Layouts live only here and in the renderer; a .m never sees them. An air
15 * field is an npts x 1 column vector laid out x-fastest — the point
16 * (ix, iy, iz) is element `ix + nx*(iy + ny*iz)` — and an object field is n1
17 * (times n2) nodes over its local 0..L1 (x 0..L2) coordinates, first axis
18 * fastest. In world coordinates every object is centred (see src/grids.ts),
19 * which is what the coupling ops assume; where an object sits in y and z is
20 * carried by the scene-built profile fields instead, so it costs no
21 * recompile to move one.
22 *
23 * `lapw` is the wall-masked Laplacian, and it is what makes a body's shell
24 * rigid. Each of the six face terms is scaled by the mask at the neighbour it
25 * reads, so a face into the wall contributes nothing — the discrete Neumann
26 * (zero normal velocity, sound-hard) condition, in the divergence form
27 * div(w grad p). The alternative, a wall as a fast material, would cut the
28 * global timestep several-fold; a mask costs nothing. Outside the domain the
29 * field is taken to be zero (a soft outer boundary); the scene's absorbing
30 * layer is meant to have swallowed the wave before it matters.
31 */
32import { UnsupportedOnGpu, WORKGROUP_SIZE } from './wgsl.ts';
33import type { AirGrid, ObjectGrid } from '../grids.ts';
35export type OpKind = 'dxx' | 'dxxxx' | 'lap2' | 'lapw' | 'radline' | 'radrect' | 'bridge';
37/** One instantiated operation: its name in the .m, its shape contract, and
38 * its WGSL. */
39export interface OpInstance {
40 name: string;
41 kind: OpKind;
42 /** Element counts of the arguments, in order. */
43 argCounts: number[];
44 outCount: number;
45 code: string;
46}
48const lineAt = (n: number): string => `
49fn at(j: i32) -> f32 {
50 if (j < 0 || j >= ${n}) { return 0.0; }
51 return src[u32(j)];
52}
53`;
55const header = (args: number): string =>
56 [
57 `@group(0) @binding(0) var<storage, read_write> dst: array<f32>;`,
58 `@group(0) @binding(1) var<storage, read> src: array<f32>;`,
59 ...(args > 1
60 ? [`@group(0) @binding(2) var<storage, read> src2: array<f32>;`]
61 : []),
62 ].join('\n');
64/** dxx_<id>: second difference along a line grid, zero beyond the ends. */
65export function makeDxx(name: string, g: ObjectGrid): OpInstance {
66 const n = g.n1;
67 const code = `${header(1)}
68${lineAt(n)}
69@compute @workgroup_size(${WORKGROUP_SIZE})
70fn main(@builtin(global_invocation_id) gid: vec3<u32>) {
71 let i = gid.x;
72 if (i >= ${n}u) { return; }
73 let j = i32(i);
74 dst[i] = ${1 / (g.h1 * g.h1)}f * (at(j - 1) - 2.0 * at(j) + at(j + 1));
75}
76`;
77 return { name, kind: 'dxx', argCounts: [n], outCount: n, code };
78}
80/**
81 * dxxxx_<id>: fourth difference along a line grid, with *antisymmetric*
82 * ghost values (u(-k) = -u(k), and mirrored at the far end). Together with
83 * the pinned endpoints this is the simply-supported end condition (u = 0 and
84 * u'' = 0), and it makes the discrete operator exactly the square of the
85 * pinned second difference — so the modes are sines and a bar's partials
86 * land where the dispersion relation f_n = kap·lam_n/2π says, which is what
87 * the spectrum overlay draws. Zero ghosts here would quietly clamp the ends
88 * instead, and the pitch knob would lie by a factor approaching two.
89 */
90export function makeDxxxx(name: string, g: ObjectGrid): OpInstance {
91 const n = g.n1;
92 const inv4 = 1 / (g.h1 * g.h1 * g.h1 * g.h1);
93 const code = `${header(1)}
94fn at(j: i32) -> f32 {
95 var k = j;
96 var s = 1.0;
97 if (k < 0) { k = -k; s = -1.0; }
98 if (k > ${n - 1}) { k = ${2 * (n - 1)} - k; s = -s; }
99 if (k < 0 || k > ${n - 1}) { return 0.0; }
100 return s * src[u32(k)];
101}
102@compute @workgroup_size(${WORKGROUP_SIZE})
103fn main(@builtin(global_invocation_id) gid: vec3<u32>) {
104 let i = gid.x;
105 if (i >= ${n}u) { return; }
106 let j = i32(i);
107 dst[i] = ${inv4}f * (at(j - 2) - 4.0 * at(j - 1) + 6.0 * at(j) - 4.0 * at(j + 1) + at(j + 2));
108}
109`;
110 return { name, kind: 'dxxxx', argCounts: [n], outCount: n, code };
111}
113/**
114 * lap_<id>: 5-point Laplacian on a rectangular grid, zero beyond the rim.
115 * The rim itself is pinned by the object's `pin` mask, so reading zero
116 * outside is consistent with the boundary the scheme imposes.
117 */
118export function makeLap2(name: string, g: ObjectGrid): OpInstance {
119 const { n1, n2, npts, h1, h2 } = g;
120 const code = `${header(1)}
121fn at(i: i32, j: i32) -> f32 {
122 if (i < 0 || i >= ${n1} || j < 0 || j >= ${n2}) { return 0.0; }
123 return src[u32(i + ${n1} * j)];
124}
125@compute @workgroup_size(${WORKGROUP_SIZE})
126fn main(@builtin(global_invocation_id) gid: vec3<u32>) {
127 let k = gid.x;
128 if (k >= ${npts}u) { return; }
129 let i = i32(k % ${n1}u);
130 let j = i32(k / ${n1}u);
131 let c = at(i, j);
132 dst[k] = ${1 / (h1 * h1)}f * (at(i - 1, j) - 2.0 * c + at(i + 1, j))
133 + ${1 / (h2 * h2)}f * (at(i, j - 1) - 2.0 * c + at(i, j + 1));
134}
135`;
136 return { name, kind: 'lap2', argCounts: [npts], outCount: npts, code };
137}
139/** lapw: the wall-masked 7-point Laplacian in the air. One instance. */
140export function makeLapw(name: string, air: AirGrid): OpInstance {
141 const { nx, ny, nz, npts, h } = air;
142 const code = `@group(0) @binding(0) var<storage, read_write> dst: array<f32>;
143@group(0) @binding(1) var<storage, read> src: array<f32>;
144@group(0) @binding(2) var<storage, read> wall: array<f32>;
146fn idx(ix: i32, iy: i32, iz: i32) -> i32 {
147 return ix + ${nx} * (iy + ${ny} * iz);
148}
149fn inside(ix: i32, iy: i32, iz: i32) -> bool {
150 return ix >= 0 && ix < ${nx} && iy >= 0 && iy < ${ny} && iz >= 0 && iz < ${nz};
151}
152// One face's contribution: masked by the wall at the neighbour, so a face
153// into the wall carries no flux — the Neumann (rigid) condition.
154fn face(ix: i32, iy: i32, iz: i32, pc: f32) -> f32 {
155 if (!inside(ix, iy, iz)) { return -pc; } // outside the domain: p = 0, open
156 let k = u32(idx(ix, iy, iz));
157 return wall[k] * (src[k] - pc);
158}
160@compute @workgroup_size(${WORKGROUP_SIZE})
161fn main(@builtin(global_invocation_id) gid: vec3<u32>) {
162 let i = gid.x;
163 if (i >= ${npts}u) { return; }
164 let ix = i32(i % ${nx}u);
165 let iy = i32((i / ${nx}u) % ${ny}u);
166 let iz = i32(i / ${nx * ny}u);
167 let pc = src[i];
168 let s = face(ix - 1, iy, iz, pc) + face(ix + 1, iy, iz, pc)
169 + face(ix, iy - 1, iz, pc) + face(ix, iy + 1, iz, pc)
170 + face(ix, iy, iz - 1, pc) + face(ix, iy, iz + 1, pc);
171 dst[i] = ${1 / (h * h)}f * wall[i] * s;
172}
173`;
174 return { name, kind: 'lapw', argCounts: [npts, npts], outCount: npts, code };
175}
177/**
178 * rad_<id> for a line object: give each air point the object's value at that
179 * point's own x (linearly interpolated, zero beyond the ends). Multiplied by
180 * a scene-built profile, that is the object radiating from its line.
181 */
182export function makeRadLine(name: string, g: ObjectGrid, air: AirGrid): OpInstance {
183 const n = g.n1;
184 const code = `${header(1)}
185@compute @workgroup_size(${WORKGROUP_SIZE})
186fn main(@builtin(global_invocation_id) gid: vec3<u32>) {
187 let i = gid.x;
188 if (i >= ${air.npts}u) { return; }
189 let ix = i % ${air.nx}u;
190 let xw = ${-air.Lx / 2}f + (f32(ix) + 0.5) * ${air.h}f;
191 let u = (xw - ${-g.L1 / 2}f) / ${g.h1}f;
192 let j = i32(floor(u));
193 if (j < 0 || j >= ${n - 1}) { dst[i] = 0.0; return; }
194 let f = u - f32(j);
195 dst[i] = mix(src[u32(j)], src[u32(j) + 1u], f);
196}
197`;
198 return { name, kind: 'radline', argCounts: [n], outCount: air.npts, code };
199}
201/**
202 * rad_<id> for a rectangular object: give each air point the object's value
203 * at that point's own (x, y), bilinearly interpolated, zero beyond the rim.
204 * Multiplied by a scene-built sheet profile, that is the membrane radiating
205 * from its plane.
206 */
207export function makeRadRect(name: string, g: ObjectGrid, air: AirGrid): OpInstance {
208 const { n1, n2 } = g;
209 const code = `${header(1)}
210@compute @workgroup_size(${WORKGROUP_SIZE})
211fn main(@builtin(global_invocation_id) gid: vec3<u32>) {
212 let i = gid.x;
213 if (i >= ${air.npts}u) { return; }
214 let ix = i % ${air.nx}u;
215 let iy = (i / ${air.nx}u) % ${air.ny}u;
216 let xw = ${-air.Lx / 2}f + (f32(ix) + 0.5) * ${air.h}f;
217 let yw = ${-air.Ly / 2}f + (f32(iy) + 0.5) * ${air.h}f;
218 let u = (xw - ${-g.L1 / 2}f) / ${g.h1}f;
219 let v = (yw - ${-g.L2 / 2}f) / ${g.h2}f;
220 let ju = i32(floor(u));
221 let jv = i32(floor(v));
222 if (ju < 0 || ju >= ${n1 - 1} || jv < 0 || jv >= ${n2 - 1}) { dst[i] = 0.0; return; }
223 let fu = u - f32(ju);
224 let fv = v - f32(jv);
225 let a = mix(src[u32(ju + ${n1} * jv)], src[u32(ju + 1 + ${n1} * jv)], fu);
226 let b = mix(src[u32(ju + ${n1} * (jv + 1))], src[u32(ju + 1 + ${n1} * (jv + 1))], fu);
227 dst[i] = mix(a, b, fv);
228}
229`;
230 return { name, kind: 'radrect', argCounts: [g.npts], outCount: air.npts, code };
231}
233/**
234 * brg_<id>: du/dx at the last node of a line object, broadcast to every air
235 * point. The end is pinned, so this is the arriving slope — what the tension
236 * pulls the bridge with. Multiplied by a scene-built patch profile, that is
237 * the bridge force driving the top plate.
238 */
239export function makeBridge(name: string, g: ObjectGrid, air: AirGrid): OpInstance {
240 const n = g.n1;
241 const code = `${header(1)}
242@compute @workgroup_size(${WORKGROUP_SIZE})
243fn main(@builtin(global_invocation_id) gid: vec3<u32>) {
244 let i = gid.x;
245 if (i >= ${air.npts}u) { return; }
246 dst[i] = (src[${n - 1}u] - src[${n - 2}u]) * ${1 / g.h1}f;
247}
248`;
249 return { name, kind: 'bridge', argCounts: [n], outCount: air.npts, code };
250}
252/** Compiled op pipelines for one scene's set of instances. Shared by every
253 * plan compiled against that scene. */
254export class OpPlan {
255 #device: GPUDevice;
256 #instances = new Map<string, OpInstance>();
257 #layouts = new Map<number, GPUBindGroupLayout>();
258 #pipelines = new Map<string, GPUComputePipeline>();
260 constructor(device: GPUDevice, instances: OpInstance[]) {
261 this.#device = device;
262 for (const inst of instances) {
263 if (this.#instances.has(inst.name)) {
264 throw new UnsupportedOnGpu(`duplicate external op '${inst.name}'`);
265 }
266 this.#instances.set(inst.name, inst);
267 }
268 }
270 /** Every instance, for the type-rule files and the editor's highlighter. */
271 get instances(): OpInstance[] {
272 return [...this.#instances.values()];
273 }
275 has(name: string): boolean {
276 return this.#instances.has(name);
277 }
279 /** The op's shape contract, for the planner's checks. */
280 spec(name: string): { argCounts: number[]; outCount: number } {
281 const inst = this.#instances.get(name);
282 if (!inst) throw new UnsupportedOnGpu(`unknown external op '${name}'`);
283 return { argCounts: inst.argCounts, outCount: inst.outCount };
284 }
286 #layout(args: number): GPUBindGroupLayout {
287 const existing = this.#layouts.get(args);
288 if (existing) return existing;
289 const entries: GPUBindGroupLayoutEntry[] = [
290 { binding: 0, visibility: GPUShaderStage.COMPUTE, buffer: { type: 'storage' } },
291 ];
292 for (let i = 0; i < args; i++) {
293 entries.push({
294 binding: i + 1,
295 visibility: GPUShaderStage.COMPUTE,
296 buffer: { type: 'read-only-storage' },
297 });
298 }
299 const layout = this.#device.createBindGroupLayout({ label: `op-${args}`, entries });
300 this.#layouts.set(args, layout);
301 return layout;
302 }
304 /** The pipeline for one op, compiled on first use. */
305 async pipeline(name: string): Promise<{ pipeline: GPUComputePipeline; outCount: number }> {
306 const inst = this.#instances.get(name);
307 if (!inst) throw new UnsupportedOnGpu(`unknown external op '${name}'`);
308 const existing = this.#pipelines.get(name);
309 if (existing) return { pipeline: existing, outCount: inst.outCount };
310 const module = this.#device.createShaderModule({ code: inst.code, label: name });
311 let pipeline: GPUComputePipeline;
312 try {
313 pipeline = await this.#device.createComputePipelineAsync({
314 layout: this.#device.createPipelineLayout({
315 bindGroupLayouts: [this.#layout(inst.argCounts.length)],
316 }),
317 compute: { module, entryPoint: 'main' },
318 label: name,
319 });
320 } catch (e) {
321 // No error scope here either; see makePipeline in plan.ts.
322 throw new UnsupportedOnGpu(
323 `op '${name}': ${e instanceof Error ? e.message : String(e)}`,
324 );
325 }
326 this.#pipelines.set(name, pipeline);
327 return { pipeline, outCount: inst.outCount };
328 }
330 createBinding(name: string, srcs: GPUBuffer[], dst: GPUBuffer): GPUBindGroup {
331 const entries: GPUBindGroupEntry[] = [{ binding: 0, resource: { buffer: dst } }];
332 srcs.forEach((s, i) => entries.push({ binding: i + 1, resource: { buffer: s } }));
333 return this.#device.createBindGroup({ layout: this.#layout(srcs.length), entries });
334 }
336 workgroups(outCount: number): number {
337 return Math.ceil(outCount / WORKGROUP_SIZE);
338 }
339}