1/**
2 * The computational grid: a square, cell-centred, uniform in both directions.
3 *
4 * Fields are flattened x-fastest — the point (ix, iy) is element `ix + nx*iy`
5 * — which is the order the stencil shader indexes in and the order the
6 * renderer reads a row of pixels in. Everything else in the project treats a
7 * field as an opaque npts x 1 column vector.
8 *
9 * Cell-centred rather than node-centred so that no grid point sits exactly on
10 * the outer boundary: the stencil takes the field outside the domain to be
11 * zero, and the absorbing layer is meant to have swallowed the wave before it
12 * gets there.
13 */
14export interface Grid {
15 n: number;
16 nx: number;
17 ny: number;
18 npts: number;
19 /** Side length of the square domain, in metres, centred on the origin. */
20 L: number;
21 /** Grid spacing, L/n, in metres. */
22 h: number;
23 /** Coordinates of every point, npts each, x fastest — as the shaders see
24 * them (f32) and as the scene .m is evaluated at (f64). */
25 x: Float32Array;
26 y: Float32Array;
27 x64: Float64Array;
28 y64: Float64Array;
29}
31export function makeGrid(n: number, L: number): Grid {
32 const h = L / n;
33 const npts = n * n;
34 const x64 = new Float64Array(npts);
35 const y64 = new Float64Array(npts);
36 for (let iy = 0; iy < n; iy++) {
37 const yv = -L / 2 + (iy + 0.5) * h;
38 for (let ix = 0; ix < n; ix++) {
39 const k = ix + n * iy;
40 x64[k] = -L / 2 + (ix + 0.5) * h;
41 y64[k] = yv;
42 }
43 }
44 return {
45 n,
46 nx: n,
47 ny: n,
48 npts,
49 L,
50 h,
51 x: new Float32Array(x64),
52 y: new Float32Array(y64),
53 x64,
54 y64,
55 };
56}
58/**
59 * The timestep the explicit leapfrog is stable at.
60 *
61 * All in SI: `h` in metres, `cmax` in metres per second, the result in
62 * seconds.
63 *
64 * Leapfrog on p_tt = c^2 L p is stable while dt^2 c^2 |L|max <= 4, and the
65 * discrete Laplacian's extreme eigenvalue is what differs between the
66 * stencils: 8/h^2 for the 5-point one, 32/(3 h^2) for the 9-point
67 * fourth-order one. That gives c*dt/h <= 1/sqrt(2) = 0.707 and
68 * c*dt/h <= sqrt(3/8) = 0.612 respectively. `cfl` is the fraction of that
69 * limit to run at, and `cmax` is the fastest sound speed anywhere in the
70 * medium — a scatterer faster than the background sets the timestep for the
71 * whole grid.
72 */
73export function stableDt(h: number, cmax: number, order: 2 | 4, cfl = 0.5): number {
74 const limit = order === 2 ? Math.SQRT1_2 : Math.sqrt(3 / 8);
75 return (cfl * limit * h) / Math.max(cmax, 1e-12);
76}