2 * The computational grid: a cube, cell-centred, uniform in all three
3 * directions.
4 *
5 * Fields are flattened x-fastest then y — the point (ix, iy, iz) is element
6 * `ix + n*iy + n*n*iz` — which is the order the stencil shader indexes in and
7 * the order the renderer reads a voxel in. Everything else treats a field as
8 * an opaque npts-long array of f32.
9 *
10 * Cell-centred rather than node-centred so that no grid point sits exactly on
11 * the outer boundary: the stencil takes the field outside the domain to be
12 * zero, and the absorbing layer is meant to have swallowed the wave before it
13 * gets there.
14 *
15 * `n` is required to be a multiple of 16 because the update shader dispatches
16 * 8x8x4 workgroups and does not want a partial one at the far corner. That is
17 * a convenience, not a physical constraint.
18 */
19export interface Grid {
20 n: number;
21 npts: number;
22 /** Side length of the cubic domain, in metres, centred on the origin. */
23 L: number;
24 /** Grid spacing, L/n, in metres. */
25 h: number;
26}
28export function makeGrid(n: number, L: number): Grid {
29 return { n, npts: n * n * n, L, h: L / n };
30}
32/** The coordinate of cell index `i` along one axis, in metres. */
33export const coord = (i: number, g: Grid): number => -g.L / 2 + (i + 0.5) * g.h;
35/**
36 * The timestep the explicit leapfrog is stable at, in seconds.
37 *
38 * Leapfrog on p_tt = c^2 lap(p) is stable while dt^2 c^2 |lap|max <= 4, and
39 * the 7-point Laplacian's extreme eigenvalue is 12/h^2, so the condition is
40 * c*dt/h <= 1/sqrt(3) = 0.577. That is stricter than the 0.707 of the flat
41 * case, which is the usual price of a dimension. `cfl` is the fraction of the
42 * limit to run at, and `cmax` is the fastest sound speed anywhere in the
43 * medium: a scatterer faster than the background sets the timestep for the
44 * whole grid.
45 */
46export function stableDt(h: number, cmax: number, cfl = 0.5): number {
47 return (cfl * h) / (Math.sqrt(3) * Math.max(cmax, 1e-12));
48}