1/**
2 * The instrument's body and the coupling geometry: a .m scene file, evaluated
3 * once on the air grid.
4 *
5 * A scene file is ordinary MATLAB defining one function,
6 *
7 * function [c, sig, wall, lineprof, boardprof] = medium(x, y, z, <...>)
8 *
9 * over the solver's grid points. Unlike the model it is *not* compiled to
10 * WGSL: the step runs every timestep and must lower to a fixed sequence of
11 * GPU dispatches, but a scene is evaluated exactly once and survives only as
12 * five arrays of numbers. So it runs through numbl's CPU interpreter instead,
13 * which buys the full MATLAB subset — loops, `if`, indexing, anything in
14 * tools/ — and f64 evaluation, where the step dialect is element-wise f32.
15 *
16 * The five outputs:
17 * c sound speed, m/s (normally just c0 everywhere — the walls are
18 * a mask, not a material, so they cost no timestep)
19 * sig absorption rate, 1/s (the sponge at the domain edge, plus any
20 * absorption the body's surfaces are given)
21 * wall 1 in air, 0 in the body's solid shell; what lapw masks by
22 * lineprof where the string radiates directly (a tube around its line)
23 * boardprof where the bridge force drives the air (a patch on the top plate)
24 */
25import { parseMFile, type FunctionStmt } from 'numbl-src/numbl-core/parser/index.ts';
26import { executeCode } from 'numbl-src/numbl-core/executeCode.ts';
27import {
28 RuntimeTensor,
29 isRuntimeTensor,
30 type RuntimeValue,
31} from 'numbl-src/numbl-core/runtime/types.ts';
32import type { AirGrid } from '../grid.ts';
33import { C_AIR } from '../units.ts';
34import { toolFiles } from '../tools.ts';
35import { inFunction, inModel, ModelCompileError } from '../mgpu/errors.ts';
36import type { MediumFields, ModelParams } from '../mgpu/model.ts';
38/** The function a scene file must define. */
39export const MEDIUM_FN = 'medium';
41const OUTPUTS = ['c', 'sig', 'wall', 'lineprof', 'boardprof'] as const;
43export interface SceneOptions {
44 air: AirGrid;
45 /** String length, metres — the scene builds the coupling profiles around
46 * the string, so it needs to know where the string is. */
47 Ls: number;
48 /** Scene source (.m text). */
49 source: string;
50 /** Parameter names the .m may take beyond the grid's own. */
51 paramNames: string[];
52 params: ModelParams;
53}
55export class Scene implements MediumFields {
56 readonly c: Float32Array;
57 readonly sig: Float32Array;
58 readonly wall: Float32Array;
59 readonly lineprof: Float32Array;
60 readonly boardprof: Float32Array;
61 readonly cmin: number;
62 readonly cmax: number;
63 /** The background speed, taken from a corner of the domain — inside the
64 * absorbing layer, where a scene has no business putting anything. */
65 readonly cref: number;
67 private constructor(fields: Record<(typeof OUTPUTS)[number], Float32Array>) {
68 this.c = fields.c;
69 this.sig = fields.sig;
70 this.wall = fields.wall;
71 this.lineprof = fields.lineprof;
72 this.boardprof = fields.boardprof;
73 let lo = Infinity;
74 let hi = 0;
75 for (const v of this.c) {
76 if (v < lo) lo = v;
77 if (v > hi) hi = v;
78 }
79 this.cmin = lo;
80 this.cmax = hi;
81 this.cref = this.c[0];
82 }
84 static create(opts: SceneOptions): Scene {
85 const fields = evaluateScene(opts);
86 const { c, sig, wall } = fields;
87 for (let i = 0; i < c.length; i++) {
88 if (!(c[i] > 0)) {
89 throw new ModelCompileError(
90 `the scene's sound speed is ${c[i]} somewhere; it must be positive ` +
91 `everywhere (the timestep is set by the fastest point, and a zero ` +
92 `or negative speed has no wave equation)`,
93 { fn: MEDIUM_FN },
94 );
95 }
96 if (!(sig[i] >= 0)) {
97 throw new ModelCompileError(
98 `the scene's absorption is ${sig[i]} somewhere; it must be zero or ` +
99 `positive (a negative one would amplify rather than absorb)`,
100 { fn: MEDIUM_FN },
101 );
102 }
103 // The mask multiplies Laplacian fluxes; outside [0, 1] it would add
104 // energy or invert a face. Clamp rather than refuse: a smoothed
105 // difference of indicators dips a hair below zero in f64 routinely.
106 if (wall[i] < 0) wall[i] = 0;
107 else if (wall[i] > 1) wall[i] = 1;
108 }
109 return new Scene(fields);
110 }
111}
113/**
114 * Evaluate the scene file on the grid, through numbl's CPU interpreter.
115 *
116 * The .m keeps the same contract the model has: it names the arguments it
117 * wants — the coordinates, the domain's numbers, and any of the registry's
118 * parameters — and the host supplies them by name, so their order in the
119 * signature is the .m's own business.
120 */
121function evaluateScene(
122 opts: SceneOptions,
123): Record<(typeof OUTPUTS)[number], Float32Array> {
124 const { air, Ls, source, paramNames, params } = opts;
125 const file = `${MEDIUM_FN}.m`;
126 const ast = inModel(() => parseMFile(source, file));
127 const fn = ast.body.find(
128 (s): s is FunctionStmt =>
129 s.type === 'Function' && (s as FunctionStmt).name === MEDIUM_FN,
130 );
131 if (!fn) {
132 throw new ModelCompileError(`the scene defines no function named '${MEDIUM_FN}'`);
133 }
134 if (fn.outputs.length !== OUTPUTS.length) {
135 throw new ModelCompileError(
136 `'${MEDIUM_FN}' must return ${OUTPUTS.length} outputs ` +
137 `[${OUTPUTS.join(', ')}], not ${fn.outputs.length}`,
138 { fn: MEDIUM_FN, start: fn.span.start, end: fn.span.end },
139 );
140 }
141 // What the grid offers a scene by name, beyond its own parameters: the
142 // coordinates in metres, the numbers that describe the domain, the speed
143 // of sound in air, and the string's length — the scene builds the coupling
144 // profiles around the string, so it needs to know where the string is.
145 const vars: Record<string, RuntimeValue> = {
146 x: new RuntimeTensor(air.x64, [air.npts, 1]),
147 y: new RuntimeTensor(air.y64, [air.npts, 1]),
148 z: new RuntimeTensor(air.z64, [air.npts, 1]),
149 Lx: air.Lx,
150 Ly: air.Ly,
151 Lz: air.Lz,
152 h: air.h,
153 c0: C_AIR,
154 npts: air.npts,
155 nx: air.nx,
156 ny: air.ny,
157 nz: air.nz,
158 Ls,
159 };
160 const known = new Set([...Object.keys(vars), ...paramNames]);
161 for (const p of fn.params) {
162 if (!known.has(p)) {
163 throw new ModelCompileError(
164 `'${MEDIUM_FN}' takes an argument '${p}' that is neither the grid ` +
165 `(${Object.keys(vars).join(', ')}) nor one of this scene's parameters` +
166 (paramNames.length ? ` (${paramNames.join(', ')})` : ''),
167 { fn: MEDIUM_FN, start: fn.span.start, end: fn.span.end },
168 );
169 }
170 }
172 for (const name of paramNames) {
173 const v = params[name];
174 // Missing parameters read as 0, as ModelPlan.setParams has it.
175 vars[name] = Number.isFinite(v) ? v : 0;
176 }
178 const outNames = OUTPUTS.map((o) => `${o}__`);
179 const driver = `[${outNames.join(', ')}] = ${MEDIUM_FN}(${fn.params.join(', ')});`;
180 const result = inFunction(MEDIUM_FN, () =>
181 executeCode(
182 driver,
183 { initialVariableValues: vars, displayResults: false, implicitCwdPath: null },
184 [...toolFiles, { name: file, source }],
185 'scene-driver.m',
186 ),
187 );
189 const fields = {} as Record<(typeof OUTPUTS)[number], Float32Array>;
190 OUTPUTS.forEach((name, i) => {
191 fields[name] = toGridField(
192 result.variableValues[outNames[i]],
193 fn.outputs[i],
194 air.npts,
195 );
196 });
197 return fields;
198}
200/** One returned field -> npts values, rounded to the solver's f32. */
201function toGridField(
202 value: RuntimeValue | undefined,
203 name: string,
204 npts: number,
205): Float32Array {
206 // A uniform field stays scalar in MATLAB; spread it over the grid.
207 if (typeof value === 'number') return new Float32Array(npts).fill(value);
208 if (value !== undefined && isRuntimeTensor(value)) {
209 if (value.imag) {
210 throw new ModelCompileError(
211 `the scene's '${name}' is complex; the medium must be real`,
212 { fn: MEDIUM_FN },
213 );
214 }
215 // A vector of npts values, either orientation. A reshape is refused
216 // rather than reordered: the tensor's column-major layout would not match
217 // the grid's x-fastest order.
218 if (value.data.length === npts && value.shape.every((d) => d === 1 || d === npts)) {
219 return new Float32Array(value.data);
220 }
221 throw new ModelCompileError(
222 `the scene's '${name}' is ${value.shape.join(' x ')}, but the grid wants ` +
223 `one value per point (${npts} x 1)`,
224 { fn: MEDIUM_FN },
225 );
226 }
227 throw new ModelCompileError(`the scene's '${name}' is not numeric`, { fn: MEDIUM_FN });
228}