1/**
2 * The medium: a .m scene file, evaluated once on the grid.
3 *
4 * A scene file is ordinary MATLAB defining one function,
5 *
6 * function [c, sig] = medium(x, y, <parameters>)
7 *
8 * over the solver's grid points. Unlike the models it is *not* compiled to
9 * WGSL: a model's step runs every frame and must lower to a fixed sequence of
10 * GPU dispatches, but a scene is evaluated exactly once and survives only as
11 * two arrays of numbers. So it runs through numbl's CPU interpreter instead,
12 * which buys the full MATLAB subset — loops, `if`, reductions, indexing,
13 * seeded randomness via `rng`/`randn`, and anything in tools/ — and f64
14 * evaluation, where the step dialect is element-wise f32.
15 *
16 * `c` is the sound speed in metres per second and `sig` the absorption rate in
17 * inverse seconds. Both are ordinary
18 * fields of position, which is what lets one function describe both the
19 * scatterer and the open boundary: the absorbing layer every scene puts around
20 * the edge (tools/sponge.m) is just the statement that the medium swallows
21 * sound out there. A scene is free to put absorption inside the domain too,
22 * which makes a lossy scatterer.
23 */
24import { parseMFile, type FunctionStmt } from 'numbl-src/numbl-core/parser/index.ts';
25import { executeCode } from 'numbl-src/numbl-core/executeCode.ts';
26import {
27 RuntimeTensor,
28 isRuntimeTensor,
29 type RuntimeValue,
30} from 'numbl-src/numbl-core/runtime/types.ts';
31import type { Grid } from '../grid.ts';
32import { C_AIR } from '../units.ts';
33import { toolFiles } from '../tools.ts';
34import { inFunction, inModel, ModelCompileError } from '../mgpu/errors.ts';
35import type { ModelParams } from '../mgpu/model.ts';
37/** The function a scene file must define. */
38export const MEDIUM_FN = 'medium';
40export interface SceneOptions {
41 grid: Grid;
42 /** Scene source (.m text). */
43 source: string;
44 /** Parameter names the .m may take beyond `x` and `y`. */
45 paramNames: string[];
46 params: ModelParams;
47}
49export class Scene {
50 /** Sound speed on the grid, npts. */
51 readonly c: Float32Array;
52 /** Absorption rate on the grid, npts. */
53 readonly sig: Float32Array;
54 readonly cmin: number;
55 readonly cmax: number;
56 /** The background speed, taken to be whatever it is in the corner of the
57 * domain — which is inside the absorbing layer, where a scene has no
58 * business putting a scatterer. What the renderer shades departures from. */
59 readonly cref: number;
60 /** The largest departure from `cref` anywhere, so the renderer's wash has a
61 * scale. Zero for a uniform medium, which draws nothing. */
62 readonly cdev: number;
64 private constructor(c: Float32Array, sig: Float32Array) {
65 this.c = c;
66 this.sig = sig;
67 let lo = Infinity;
68 let hi = 0;
69 for (const v of c) {
70 if (v < lo) lo = v;
71 if (v > hi) hi = v;
72 }
73 this.cmin = lo;
74 this.cmax = hi;
75 this.cref = c[0];
76 this.cdev = Math.max(Math.abs(hi - this.cref), Math.abs(this.cref - lo));
77 }
79 static create(opts: SceneOptions): Scene {
80 const { grid, source, paramNames, params } = opts;
81 const [c, sig] = evaluateMedium(source, paramNames, params, grid);
82 for (let i = 0; i < c.length; i++) {
83 if (!(c[i] > 0)) {
84 throw new ModelCompileError(
85 `the scene's sound speed is ${c[i]} somewhere; it must be positive ` +
86 `everywhere (the timestep is set by the fastest point, and a zero ` +
87 `or negative speed has no wave equation)`,
88 { fn: MEDIUM_FN },
89 );
90 }
91 if (!(sig[i] >= 0)) {
92 throw new ModelCompileError(
93 `the scene's absorption is ${sig[i]} somewhere; it must be zero or ` +
94 `positive (a negative one would amplify rather than absorb)`,
95 { fn: MEDIUM_FN },
96 );
97 }
98 }
99 return new Scene(c, sig);
100 }
101}
103/**
104 * Evaluate the scene file on the grid, through numbl's CPU interpreter.
105 *
106 * The .m keeps the same contract a model has: it names the arguments it wants
107 * — `x`, `y`, and any of the registry's parameters — and the host supplies
108 * them by name, so their order in the signature is the .m's own business. A
109 * one-line driver script calls `medium` with exactly the arguments its
110 * signature declares, with those names pre-bound in the driver's workspace.
111 */
112function evaluateMedium(
113 source: string,
114 paramNames: string[],
115 params: ModelParams,
116 grid: Grid,
117): [Float32Array, Float32Array] {
118 const file = `${MEDIUM_FN}.m`;
119 const ast = inModel(() => parseMFile(source, file));
120 const fn = ast.body.find(
121 (s): s is FunctionStmt =>
122 s.type === 'Function' && (s as FunctionStmt).name === MEDIUM_FN,
123 );
124 if (!fn) {
125 throw new ModelCompileError(`the scene defines no function named '${MEDIUM_FN}'`);
126 }
127 if (fn.outputs.length !== 2) {
128 throw new ModelCompileError(
129 `'${MEDIUM_FN}' must return two outputs [c, sig] — the sound speed and ` +
130 `the absorption rate — not ${fn.outputs.length}`,
131 { fn: MEDIUM_FN, start: fn.span.start, end: fn.span.end },
132 );
133 }
134 // What the grid offers a scene by name, beyond its own parameters: the
135 // coordinates themselves, in metres, and the numbers that describe them. A
136 // scene wants `L` to place its absorbing layer relative to the domain and
137 // `h` to smooth an interface over about a cell, so both are part of the
138 // contract rather than something every scene has to be told. `c0` is the
139 // speed of sound in air, which is what "the background" means here.
140 const vars: Record<string, RuntimeValue> = {
141 x: new RuntimeTensor(grid.x64, [grid.npts, 1]),
142 y: new RuntimeTensor(grid.y64, [grid.npts, 1]),
143 L: grid.L,
144 h: grid.h,
145 c0: C_AIR,
146 npts: grid.npts,
147 nx: grid.nx,
148 ny: grid.ny,
149 };
150 const known = new Set([...Object.keys(vars), ...paramNames]);
151 for (const p of fn.params) {
152 if (!known.has(p)) {
153 throw new ModelCompileError(
154 `'${MEDIUM_FN}' takes an argument '${p}' that is neither the grid ` +
155 `(${Object.keys(vars).join(', ')}) nor one of this scene's parameters` +
156 (paramNames.length ? ` (${paramNames.join(', ')})` : ''),
157 { fn: MEDIUM_FN, start: fn.span.start, end: fn.span.end },
158 );
159 }
160 }
162 for (const name of paramNames) {
163 const v = params[name];
164 // Missing parameters read as 0, as ModelPlan.setParams has it.
165 vars[name] = Number.isFinite(v) ? v : 0;
166 }
168 const driver = `[c__, sig__] = ${MEDIUM_FN}(${fn.params.join(', ')});`;
169 const result = inFunction(MEDIUM_FN, () =>
170 executeCode(
171 driver,
172 { initialVariableValues: vars, displayResults: false, implicitCwdPath: null },
173 [...toolFiles, { name: file, source }],
174 'scene-driver.m',
175 ),
176 );
178 return [
179 toGridField(result.variableValues['c__'], fn.outputs[0], grid.npts),
180 toGridField(result.variableValues['sig__'], fn.outputs[1], grid.npts),
181 ];
182}
184/** One returned field -> npts values, rounded to the solver's f32. */
185function toGridField(
186 value: RuntimeValue | undefined,
187 name: string,
188 npts: number,
189): Float32Array {
190 // A uniform field stays scalar in MATLAB; spread it over the grid.
191 if (typeof value === 'number') return new Float32Array(npts).fill(value);
192 if (value !== undefined && isRuntimeTensor(value)) {
193 if (value.imag) {
194 throw new ModelCompileError(
195 `the scene's '${name}' is complex; the medium must be real`,
196 { fn: MEDIUM_FN },
197 );
198 }
199 // A vector of npts values, either orientation. A 2-D reshape is refused
200 // rather than reordered: the tensor's column-major layout would not match
201 // the grid's x-fastest rows.
202 if (value.data.length === npts && value.shape.every((d) => d === 1 || d === npts)) {
203 return new Float32Array(value.data);
204 }
205 throw new ModelCompileError(
206 `the scene's '${name}' is ${value.shape.join(' x ')}, but the grid wants ` +
207 `one value per point (${npts} x 1)`,
208 { fn: MEDIUM_FN },
209 );
210 }
211 throw new ModelCompileError(`the scene's '${name}' is not numeric`, { fn: MEDIUM_FN });
212}