1/**
2 * Scene description -> the two MATLAB files, and everything around them.
3 *
4 * This is the studio's hinge. A SceneSpec (src/studio/types.ts) says *what*
5 * is in the scene; this module writes the model .m (init and step for every
6 * object plus the air, all advanced together) and the scene .m (bodies,
7 * absorption, and each object's coupling profiles), instantiates the
8 * external ops each object needs with its grid geometry baked in, and lists
9 * the constants, fields, state and parameters the host must offer the
10 * compiler.
11 *
12 * The generated files are not an internal detail: the Lab view shows them,
13 * and they are editable — an edited source compiles through exactly the same
14 * pipeline, the generator's output being nothing more than a starting point.
15 * That is why the generator works hard at readability: suffixed names
16 * (`u_s1`, `cs_s1`) rather than positional arrays, one object's physics kept
17 * together under its own comment, and the same idioms a hand-written model
18 * (dulcimer's) uses.
19 *
20 * Naming scheme, per object id ID:
21 * u_ID, um_ID the object's displacement and its history (state)
22 * xs_ID / ym_ID node coordinates (local, 0..L), pin_ID the pin mask
23 * n_ID, L_ID node count and length(s) — exact constants
24 * dxx_ID, lap_ID, the object's stencil and coupling ops
25 * rad_ID, brg_ID
26 * prof_ID, bprof_ID where in the air the object radiates / its bridge
27 * patch drives (scene-built, air-sized)
28 * The air keeps dulcimer's names: p, pm, c, sig, wall, lapw.
29 */
30import type { AirGrid, ObjectGrid } from '../grids.ts';
31import {
32 makeBridge, makeDxx, makeDxxxx, makeLap2, makeLapw, makeRadLine, makeRadRect,
33 type OpInstance,
34} from '../mgpu/ops.ts';
35import type { StateField } from '../mgpu/model.ts';
36import { BODY_TEMPLATES } from './bodies.ts';
37import { templateOf } from './objects.ts';
38import {
39 suffixed, type ObjectSpec, type ParamSpec, type Params, type SceneSpec,
40} from './types.ts';
42/** Everything the compiler pipeline needs, derived from one SceneSpec. */
43export interface GeneratedScene {
44 modelSource: string;
45 sceneSource: string;
46 ops: OpInstance[];
47 /** Exact scalars offered to the model .m. */
48 modelConsts: Record<string, number>;
49 /** Host-uploaded fields offered to the model .m, with initial contents.
50 * The scene-built ones (c, sig, wall, prof_*) are placeholders here —
51 * the session evaluates the scene and uploads the real thing. */
52 modelFields: Record<string, Float32Array>;
53 /** The scene-built subset of modelFields, in output order for medium(). */
54 sceneOutputs: string[];
55 /** Exact scalars offered to the scene .m. */
56 sceneConsts: Record<string, number>;
57 state: StateField[];
58 /** Model parameter names (suffixed), and their slider specs by object. */
59 modelParamNames: string[];
60 /** Scene parameter names (placement + body), and their specs. */
61 sceneParamNames: string[];
62 /** Each object's grid, by id. */
63 grids: Record<string, ObjectGrid>;
64}
66/** A parameter spec with its key suffixed for one object. */
67export const suffixedSpecs = (specs: ParamSpec[], id: string): ParamSpec[] =>
68 specs.map((s) => ({ ...s, key: suffixed(s.key, id) }));
70/** The parameter values one object contributes, suffixed, defaults filled. */
71export function suffixedParams(specs: ParamSpec[], values: Params, id: string): Params {
72 return Object.fromEntries(
73 specs.map((s) => [suffixed(s.key, id), values[s.key] ?? s.value]),
74 );
75}
77const line = (s: string): string => ` ${s}`;
79/** The generated model .m: init and step for every object plus the air. */
80function generateModel(spec: SceneSpec, grids: Record<string, ObjectGrid>): string {
81 const objs = spec.objects;
82 const stateOut = objs.flatMap((o) => [`u_${o.id}`, `um_${o.id}`]).concat(['p', 'pm']);
83 const stepOut = objs.flatMap((o) => [`un_${o.id}`, `uo_${o.id}`]).concat(['pn', 'po']);
85 // ---- init ---------------------------------------------------------------
86 const initArgs: string[] = [];
87 const initBody: string[] = [];
88 for (const o of objs) {
89 const id = o.id;
90 const t = templateOf(o);
91 initBody.push(`% -- ${t.label.toLowerCase()} ${id}: ${t.verb.toLowerCase()} as an initial condition`);
92 if (o.type === 'string') {
93 initArgs.push(`xs_${id}`, `L_${id}`, `pluckpos_${id}`, `amp_${id}`);
94 initBody.push(
95 `a_${id} = xs_${id} / (pluckpos_${id} * L_${id});`,
96 `b_${id} = (L_${id} - xs_${id}) / (L_${id} - pluckpos_${id} * L_${id});`,
97 `u_${id} = amp_${id} * 0.5 * (a_${id} + b_${id} - abs(a_${id} - b_${id}));`,
98 `um_${id} = u_${id};`,
99 );
100 } else if (o.type === 'bar') {
101 initArgs.push(`xs_${id}`, `pin_${id}`, `n_${id}`, `L_${id}`, `strikepos_${id}`, `striker_${id}`, `vel_${id}`, 'dt');
102 initBody.push(
103 `r_${id} = (xs_${id} - strikepos_${id} * L_${id}) / striker_${id};`,
104 `u_${id} = zeros(n_${id}, 1);`,
105 `um_${id} = -(dt * vel_${id}) * (pin_${id} .* exp(-(r_${id} .* r_${id})));`,
106 );
107 } else {
108 initArgs.push(
109 `xm_${id}`, `ym_${id}`, `pin_${id}`, `n_${id}`, `L1_${id}`, `L2_${id}`,
110 `strikex_${id}`, `strikey_${id}`, `striker_${id}`, `vel_${id}`, 'dt',
111 );
112 initBody.push(
113 `rr_${id} = ((xm_${id} - strikex_${id} * L1_${id}).^2 + (ym_${id} - strikey_${id} * L2_${id}).^2) / (striker_${id} * striker_${id});`,
114 `u_${id} = zeros(n_${id}, 1);`,
115 `um_${id} = -(dt * vel_${id}) * (pin_${id} .* exp(-rr_${id}));`,
116 );
117 }
118 initBody.push('');
119 }
120 initBody.push('% -- the air starts silent', 'p = zeros(npts, 1);', 'pm = zeros(npts, 1);');
121 initArgs.push('npts');
122 const initSig = `function [${stateOut.join(', ')}] = init(${[...new Set(initArgs)].join(', ')})`;
124 // ---- step ---------------------------------------------------------------
125 const stepArgs: string[] = [];
126 const stepBody: string[] = [];
127 for (const o of objs) {
128 const id = o.id;
129 const t = templateOf(o);
130 stepArgs.push(`u_${id}`, `um_${id}`, `pin_${id}`);
131 stepArgs.push(...(o.type === 'membrane' ? [`L1_${id}`, `L2_${id}`] : [`L_${id}`]));
132 stepBody.push(`% -- ${t.label.toLowerCase()} ${id} --------------------------------------------`);
133 if (o.type === 'string') {
134 stepArgs.push(`f0_${id}`, `B_${id}`, `t60_${id}`, `sig1_${id}`);
135 stepBody.push(
136 `cs_${id} = 2 * L_${id} * f0_${id}; % tension, as a wave speed`,
137 `kap2_${id} = B_${id} * (cs_${id} * L_${id} / pi)^2; % stiffness, as inharmonicity`,
138 `d_${id} = (6.91 / t60_${id}) * dt;`,
139 `lu_${id} = dxx_${id}(u_${id});`,
140 `lum_${id} = dxx_${id}(um_${id});`,
141 `l4_${id} = dxxxx_${id}(u_${id});`,
142 `un_${id} = pin_${id} .* ((2*u_${id} - (1 - d_${id})*um_${id} + (dt*dt)*((cs_${id}*cs_${id})*lu_${id} - kap2_${id}*l4_${id}) + (2*sig1_${id}*dt)*(lu_${id} - lum_${id})) ./ (1 + d_${id}));`,
143 );
144 } else if (o.type === 'bar') {
145 stepArgs.push(`f1_${id}`, `t60_${id}`, `sig1_${id}`);
146 stepBody.push(
147 `kap_${id} = 2 * f1_${id} * (L_${id} * L_${id}) / pi; % stiffness, from the first partial`,
148 `d_${id} = (6.91 / t60_${id}) * dt;`,
149 `lu_${id} = dxx_${id}(u_${id});`,
150 `lum_${id} = dxx_${id}(um_${id});`,
151 `l4_${id} = dxxxx_${id}(u_${id});`,
152 `un_${id} = pin_${id} .* ((2*u_${id} - (1 - d_${id})*um_${id} - (dt*dt)*(kap_${id}*kap_${id})*l4_${id} + (2*sig1_${id}*dt)*(lu_${id} - lum_${id})) ./ (1 + d_${id}));`,
153 );
154 } else {
155 stepArgs.push(`f0_${id}`, `t60_${id}`, `sig1_${id}`);
156 stepBody.push(
157 `cm_${id} = 2 * f0_${id} / sqrt(1/(L1_${id}*L1_${id}) + 1/(L2_${id}*L2_${id})); % tension, from the (1,1) mode`,
158 `d_${id} = (6.91 / t60_${id}) * dt;`,
159 `lu_${id} = lap_${id}(u_${id});`,
160 `lum_${id} = lap_${id}(um_${id});`,
161 `un_${id} = pin_${id} .* ((2*u_${id} - (1 - d_${id})*um_${id} + (dt*dt)*(cm_${id}*cm_${id})*lu_${id} + (2*sig1_${id}*dt)*(lu_${id} - lum_${id})) ./ (1 + d_${id}));`,
162 );
163 }
164 stepBody.push('');
165 }
167 stepBody.push(
168 '% -- the coupling: each object\'s motion, spread where its profile says --',
169 );
170 let drv = '';
171 let drvIdx = 0;
172 for (const o of objs) {
173 const id = o.id;
174 stepArgs.push(`gair_${id}`);
175 stepBody.push(
176 `acc_${id} = (un_${id} - 2*u_${id} + um_${id}) / (dt*dt);`,
177 `q_${id} = rad_${id}(acc_${id});`,
178 );
179 let term = `gair_${id} * (prof_${id} .* q_${id})`;
180 if (o.type === 'string') {
181 stepArgs.push(`gbridge_${id}`);
182 stepBody.push(`fb_${id} = brg_${id}(un_${id});`);
183 term += ` - (gbridge_${id} * (cs_${id}*cs_${id}) / L_${id}) * (bprof_${id} .* fb_${id})`;
184 }
185 const next = `drv${++drvIdx}`;
186 stepBody.push(`${next} = ${drv ? `${drv} + ` : ''}${term};`);
187 drv = next;
188 }
189 stepBody.push(
190 '',
191 '% -- the air ------------------------------------------------------------',
192 'lp = lapw(p, wall);',
193 'sd = sig * dt;',
194 `pn = wall .* ((2*p - (1 - sd) .* pm + (c*dt).^2 .* lp + (dt*dt) * ${drv}) ./ (1 + sd));`,
195 '',
196 "% This step's fields become the next step's history. Lines of their own,",
197 '% so each plans as the copy it is.',
198 ...objs.map((o) => `uo_${o.id} = u_${o.id};`),
199 'po = p;',
200 );
201 stepArgs.push('p', 'pm', 'c', 'sig', 'wall');
202 for (const o of objs) {
203 stepArgs.push(`prof_${o.id}`);
204 if (o.type === 'string') stepArgs.push(`bprof_${o.id}`);
205 }
206 stepArgs.push('dt');
207 const stepSig = `function [${stepOut.join(', ')}] = step(${[...new Set(stepArgs)].join(', ')})`;
209 const objectSummary = objs
210 .map((o) => {
211 const g = grids[o.id];
212 const t = templateOf(o);
213 return `% ${o.id}: ${t.label.toLowerCase()}, ${
214 g.kind === 'line' ? `${g.n1} nodes over ${g.L1} m` : `${g.n1} x ${g.n2} nodes over ${g.L1} x ${g.L2} m`
215 }`;
216 })
217 .join('\n');
219 return `% Generated by soundbox from the scene you built — and editable: this file
220% is the model, compiled to WebGPU kernels, and Run edits recompiles whatever
221% you make of it. (Rebuilding the scene from the Build tab regenerates it.)
222%
223% Objects in this scene:
224${objectSummary}
225%
226% Every object is advanced by the same explicit leapfrog family: displacement
227% u and its history um, one update per timestep, pinned at the boundary by
228% its pin mask. The air carries pressure p on the box grid,
229%
230% p_tt + 2 sig p_t = c^2 lapw(p, wall) + s
231%
232% where lapw is the wall-masked Laplacian (the rigid body, at no timestep
233% cost) and s collects every object's contribution: its acceleration, spread
234% along the scene-built profile of where it sits (and, for a string, its
235% bridge force driving a patch above the plate). The air does not push back;
236% at these amplitudes that back-reaction is far below what the models already
237% neglect.
239${initSig}
240${initBody.map(line).join('\n')}
241end
243${stepSig}
244${stepBody.map(line).join('\n')}
245end
246`;
247}
249/** The generated scene .m: bodies, absorption, coupling profiles. */
250function generateScene(spec: SceneSpec): { source: string; outputs: string[]; paramNames: string[] } {
251 const body = BODY_TEMPLATES[spec.body.type];
252 const outputs = ['c', 'sig', 'wall'];
253 const args = ['x', 'y', 'z', 'h', 'npts', 'Lx', 'Ly', 'Lz', 'c0'];
254 const paramNames: string[] = [];
256 for (const b of body.paramSpecs) {
257 args.push(b.key);
258 paramNames.push(b.key);
259 }
261 const lines: string[] = [...body.matlab, 'wall = 1 - shell;', ''];
262 lines.push(
263 '% Uniform air. The walls being a mask is what keeps this uniform — and',
264 '% the timestep as large as air allows.',
265 'c = c0;',
266 '',
267 "% Open boundary at the domain edge" +
268 (body.paramSpecs.some((p) => p.key === 'absorb')
269 ? ", plus the body's surface absorption."
270 : '.'),
271 body.paramSpecs.some((p) => p.key === 'absorb')
272 ? 'sig = sponge3(x, y, z, Lx, Ly, Lz, 0.15*Ly, 9000) + absorb * (4 * shell .* (1 - shell));'
273 : 'sig = sponge3(x, y, z, Lx, Ly, Lz, 0.15*Ly, 9000);',
274 '',
275 );
277 for (const o of spec.objects) {
278 const id = o.id;
279 const t = templateOf(o);
280 // The object's baked size(s) come in as exact constants; its placement
281 // as scene parameters.
282 args.push(...(o.type === 'membrane' ? [`L1_${id}`, `L2_${id}`] : [`L_${id}`]));
283 for (const p of t.placeSpecs) {
284 const k = suffixed(p.key, id);
285 args.push(k);
286 paramNames.push(k);
287 }
288 if (o.type === 'membrane') {
289 outputs.push(`prof_${id}`);
290 lines.push(
291 `% ${t.label.toLowerCase()} ${id} radiates from its plane: a thin sheet at z = oz_${id}`,
292 '% over its extent, normalized so its cross-section integrates to one.',
293 `aw_${id} = 2*h;`,
294 `sk_${id} = exp(-((z - oz_${id}).^2) / (aw_${id}*aw_${id})) / (sqrt(pi) * aw_${id});`,
295 `exx_${id} = 0.5 * (1 + tanh((L1_${id}/2 - abs(x)) / (1.5*h)));`,
296 `exy_${id} = 0.5 * (1 + tanh((L2_${id}/2 - abs(y)) / (1.5*h)));`,
297 `prof_${id} = sk_${id} .* exx_${id} .* exy_${id} .* wall;`,
298 '',
299 );
300 } else {
301 outputs.push(`prof_${id}`);
302 lines.push(
303 `% ${t.label.toLowerCase()} ${id} radiates from its line: a tube two cells wide`,
304 `% around (y = oy_${id}, z = oz_${id}), fading out at its ends, normalized so`,
305 '% its cross-section integrates to one.',
306 `aw_${id} = 2*h;`,
307 `tube_${id} = exp(-((y - oy_${id}).^2 + (z - oz_${id}).^2) / (aw_${id}*aw_${id})) / (pi * aw_${id} * aw_${id});`,
308 `span_${id} = 0.5 * (1 + tanh((L_${id}/2 - abs(x)) / (1.5*h)));`,
309 `prof_${id} = tube_${id} .* span_${id} .* wall;`,
310 );
311 if (o.type === 'string') {
312 outputs.push(`bprof_${id}`);
313 lines.push(
314 `% ...and its bridge force drives a patch of air just above the plate,`,
315 `% around its x = +L/2 end. This stands in for a top plate moving; the`,
316 '% plate\'s own modes are not modelled.',
317 `pw_${id} = 1.5*h;`,
318 `pr_${id} = exp(-((x - L_${id}/2).^2 + (y - oy_${id}).^2) / (0.04*0.04)) / (pi * 0.04 * 0.04);`,
319 `skin_${id} = exp(-(z - pw_${id}).^2 / (pw_${id}*pw_${id})) / (sqrt(pi) * pw_${id});`,
320 `bprof_${id} = pr_${id} .* skin_${id} .* wall;`,
321 );
322 }
323 lines.push('');
324 }
325 }
327 const source = `% Generated by soundbox — the scene: ${body.label.toLowerCase()}, and where each
328% object couples into the air. Evaluated once on the CPU (numbl's
329% interpreter, full MATLAB subset available), so editing it costs no
330% recompile; the model .m is the compiled half.
331%
332% The body is a *mask*: \`wall\` is 0 inside solid and 1 in air, and the
333% solver's masked Laplacian drops any flux into a masked cell — the rigid
334% (Neumann) condition, at no timestep cost. What a rigid mask cannot do is
335% absorb, so \`absorb\` paints a thin lossy skin over the shell's surfaces
336% (4*shell*(1-shell) peaks exactly at the surface).
337function [${outputs.join(', ')}] = medium(${args.join(', ')})
338${lines.map(line).join('\n')}
339end
340`;
341 return { source, outputs, paramNames };
342}
344/** Derive everything the pipeline needs from a scene description. */
345export function generate(
346 spec: SceneSpec,
347 air: AirGrid,
348 dt: number,
349): GeneratedScene {
350 if (spec.objects.length === 0) {
351 throw new Error('a scene needs at least one object (add a string, bar or membrane)');
352 }
353 const seen = new Set<string>();
354 for (const o of spec.objects) {
355 if (seen.has(o.id)) throw new Error(`duplicate object id '${o.id}'`);
356 seen.add(o.id);
357 }
359 const grids: Record<string, ObjectGrid> = {};
360 for (const o of spec.objects) grids[o.id] = templateOf(o).makeGrid(o.size, dt);
362 // Ops, consts, fields, state.
363 const ops: OpInstance[] = [makeLapw('lapw', air)];
364 const modelConsts: Record<string, number> = {
365 npts: air.npts, nx: air.nx, ny: air.ny, nz: air.nz,
366 };
367 const sceneConsts: Record<string, number> = {};
368 const modelFields: Record<string, Float32Array> = {
369 x: air.x, y: air.y, z: air.z,
370 };
371 const state: StateField[] = [];
372 const modelParamNames: string[] = [];
374 for (const o of spec.objects) {
375 const id = o.id;
376 const g = grids[id];
377 modelConsts[`n_${id}`] = g.npts;
378 if (g.kind === 'line') {
379 modelConsts[`L_${id}`] = g.L1;
380 sceneConsts[`L_${id}`] = g.L1;
381 modelFields[`xs_${id}`] = g.xs;
382 ops.push(makeDxx(`dxx_${id}`, g), makeDxxxx(`dxxxx_${id}`, g), makeRadLine(`rad_${id}`, g, air));
383 if (o.type === 'string') ops.push(makeBridge(`brg_${id}`, g, air));
384 } else {
385 modelConsts[`L1_${id}`] = g.L1;
386 modelConsts[`L2_${id}`] = g.L2;
387 sceneConsts[`L1_${id}`] = g.L1;
388 sceneConsts[`L2_${id}`] = g.L2;
389 modelFields[`xm_${id}`] = g.xs;
390 modelFields[`ym_${id}`] = g.ys;
391 ops.push(makeLap2(`lap_${id}`, g), makeRadRect(`rad_${id}`, g, air));
392 }
393 modelFields[`pin_${id}`] = g.pin;
394 state.push({ name: `u_${id}`, count: g.npts }, { name: `um_${id}`, count: g.npts });
395 for (const p of templateOf(o).paramSpecs) modelParamNames.push(suffixed(p.key, id));
396 }
397 state.push({ name: 'p', count: air.npts }, { name: 'pm', count: air.npts });
399 const scene = generateScene(spec);
400 // Scene-built fields are uploaded after evaluation; placeholders size them.
401 for (const name of scene.outputs) modelFields[name] = new Float32Array(air.npts);
403 return {
404 modelSource: generateModel(spec, grids),
405 sceneSource: scene.source,
406 ops,
407 modelConsts,
408 modelFields,
409 sceneOutputs: scene.outputs,
410 sceneConsts,
411 state,
412 modelParamNames,
413 sceneParamNames: scene.paramNames,
414 grids,
415 };
416}