2 * The scene: bodies, absorption, and every object's coupling geometry, as a
3 * .m file evaluated once on the air grid.
4 *
5 * A scene file is ordinary MATLAB defining one function,
6 *
7 * function [c, sig, wall, prof_s1, ...] = 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 * 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 first three outputs are fixed:
17 * c sound speed, m/s (normally just c0 everywhere — bodies are a
18 * 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 a body's surfaces are given)
21 * wall 1 in air, 0 in a body's solid shell; what lapw masks by
22 * After them come whatever coupling profiles this scene's objects need
23 * (prof_s1, brg_s1, ...), named by the caller — the studio's code generator
24 * decides the list, and the same list is offered to the compiled model as
25 * uploadable fields.
26 */
27import { parseMFile, type FunctionStmt } from 'numbl-src/numbl-core/parser/index.ts';
28import { executeCode } from 'numbl-src/numbl-core/executeCode.ts';
29import {
30 RuntimeTensor,
31 isRuntimeTensor,
32 type RuntimeValue,
33} from 'numbl-src/numbl-core/runtime/types.ts';
34import type { AirGrid } from '../grids.ts';
35import { C_AIR } from '../units.ts';
36import { toolFiles } from '../tools.ts';
37import { inFunction, inModel, ModelCompileError } from '../mgpu/errors.ts';
38import type { ModelParams } from '../mgpu/model.ts';
40/** The function a scene file must define. */
41export const MEDIUM_FN = 'medium';
43export interface SceneOptions {
44 air: AirGrid;
45 /** Scene source (.m text). */
46 source: string;
47 /** Output names, in order: 'c', 'sig', 'wall', then the profiles. */
48 outputs: string[];
49 /** Exact scalars offered by name beyond the grid's own (object lengths and
50 * the like). */
51 consts: Record<string, number>;
52 /** Parameter names the .m may take beyond the grid and the consts. */
53 paramNames: string[];
54 params: ModelParams;
55}
57export class Scene {
58 /** Every output field, air-sized, by name. */
59 readonly fields: Record<string, Float32Array>;
60 readonly cmin: number;
61 readonly cmax: number;
63 private constructor(fields: Record<string, Float32Array>) {
64 this.fields = fields;
65 let lo = Infinity;
66 let hi = 0;
67 for (const v of fields.c) {
68 if (v < lo) lo = v;
69 if (v > hi) hi = v;
70 }
71 this.cmin = lo;
72 this.cmax = hi;
73 }
75 get c(): Float32Array {
76 return this.fields.c;
77 }
78 get sig(): Float32Array {
79 return this.fields.sig;
80 }
81 get wall(): Float32Array {
82 return this.fields.wall;
83 }
85 static create(opts: SceneOptions): Scene {
86 const fields = evaluateScene(opts);
87 const { c, sig, wall } = fields;
88 for (let i = 0; i < c.length; i++) {
89 if (!(c[i] > 0)) {
90 throw new ModelCompileError(
91 `the scene's sound speed is ${c[i]} somewhere; it must be positive ` +
92 `everywhere (the timestep is set by the fastest point, and a zero ` +
93 `or negative speed has no wave equation)`,
94 { fn: MEDIUM_FN },
95 );
96 }
97 if (!(sig[i] >= 0)) {
98 throw new ModelCompileError(
99 `the scene's absorption is ${sig[i]} somewhere; it must be zero or ` +
100 `positive (a negative one would amplify rather than absorb)`,
101 { fn: MEDIUM_FN },
102 );
103 }
104 // The mask multiplies Laplacian fluxes; outside [0, 1] it would add
105 // energy or invert a face. Clamp rather than refuse: a smoothed
106 // difference of indicators dips a hair below zero in f64 routinely.
107 if (wall[i] < 0) wall[i] = 0;
108 else if (wall[i] > 1) wall[i] = 1;
109 }
110 return new Scene(fields);
111 }
112}
114/**
115 * Evaluate the scene file on the grid, through numbl's CPU interpreter.
116 *
117 * The .m keeps the same contract the model has: it names the arguments it
118 * wants — the coordinates, the domain's numbers, and any of the offered
119 * parameters — and the host supplies them by name, so their order in the
120 * signature is the .m's own business.
121 */
122function evaluateScene(opts: SceneOptions): Record<string, Float32Array> {
123 const { air, source, outputs, paramNames, params } = opts;
124 const file = `${MEDIUM_FN}.m`;
125 const ast = inModel(() => parseMFile(source, file));
126 const fn = ast.body.find(
127 (s): s is FunctionStmt =>
128 s.type === 'Function' && (s as FunctionStmt).name === MEDIUM_FN,
129 );
130 if (!fn) {
131 throw new ModelCompileError(`the scene defines no function named '${MEDIUM_FN}'`);
132 }
133 if (fn.outputs.length !== outputs.length) {
134 throw new ModelCompileError(
135 `'${MEDIUM_FN}' must return ${outputs.length} outputs ` +
136 `[${outputs.join(', ')}], not ${fn.outputs.length}`,
137 { fn: MEDIUM_FN, start: fn.span.start, end: fn.span.end },
138 );
139 }
140 // What the grid offers a scene by name, beyond its own parameters: the
141 // coordinates in metres, the numbers that describe the domain, the speed
142 // of sound in air, and whatever consts the caller adds (object lengths —
143 // the scene builds the coupling profiles around each object, so it needs
144 // to know how big each one 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 ...opts.consts,
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: Record<string, Float32Array> = {};
190 outputs.forEach((name, i) => {
191 fields[name] = toGridField(
192 result.variableValues[outNames[i]],
193 name,
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}