2 * The showcase problem: a point source at x0 inside a rounded cube (the
3 * superellipsoid |x|^p + |y|^p + |z|^p = a^p) with impedance walls. With time
4 * dependence exp(-i omega t) the pressure amplitude u satisfies
5 *
6 * (Delta + k^2) u = -delta(x - x0) in the room,
7 * du/dn + i k lambda u = 0 on the walls,
8 *
9 * with n the outward normal and lambda = -beta, where beta > 0 is the wall
10 * admittance. We write u = G + v with G the free-space field of the source and
11 * find v from fmm3dbie's right-preconditioned combined field equation.
12 *
13 * Complex vectors are interleaved Float64Arrays: [re0, im0, re1, im1, ...].
14 */
16export type CVec = Float64Array;
18export const ALPHA = 1.0;
20/** Free-space Green's function exp(ikr)/(4 pi r) at points x (3, n) for a source x0. */
21export function green(x: Float64Array, x0: ArrayLike<number>, k: number): CVec {
22 const n = x.length / 3;
23 const out = new Float64Array(2 * n);
24 for (let i = 0; i < n; i++) {
25 const r = Math.hypot(x[i] - x0[0], x[n + i] - x0[1], x[2 * n + i] - x0[2]);
26 const s = 1 / (4 * Math.PI * r);
27 out[2 * i] = s * Math.cos(k * r);
28 out[2 * i + 1] = s * Math.sin(k * r);
29 }
30 return out;
31}
33/** Normal derivative of the Green's function at x with unit normals nrm (3, n). */
34export function greenDn(x: Float64Array, nrm: Float64Array, x0: ArrayLike<number>, k: number): CVec {
35 const n = x.length / 3;
36 const out = new Float64Array(2 * n);
37 for (let i = 0; i < n; i++) {
38 const dx = x[i] - x0[0];
39 const dy = x[n + i] - x0[1];
40 const dz = x[2 * n + i] - x0[2];
41 const r = Math.hypot(dx, dy, dz);
42 const s = 1 / (4 * Math.PI * r);
43 const gr = s * Math.cos(k * r);
44 const gi = s * Math.sin(k * r);
45 const dn = (dx * nrm[i] + dy * nrm[n + i] + dz * nrm[2 * n + i]) / r;
46 // g * (ik - 1/r) * dn
47 const fr = -1 / r;
48 const fi = k;
49 out[2 * i] = (gr * fr - gi * fi) * dn;
50 out[2 * i + 1] = (gr * fi + gi * fr) * dn;
51 }
52 return out;
53}
55export function inside(x: ArrayLike<number>, a: number, p: number, margin = 0): boolean {
56 return Math.pow(Math.abs(x[0] / a), p) + Math.pow(Math.abs(x[1] / a), p) + Math.pow(Math.abs(x[2] / a), p) < (1 - margin) ** p;
57}
59/** Impedance data for v = u - G: -(dG/dn + i k lambda G) at the wall nodes. */
60export function roomRhs(xs: Float64Array, nrm: Float64Array, x0: ArrayLike<number>, k: number, lam: number): CVec {
61 const g = green(xs, x0, k);
62 const gn = greenDn(xs, nrm, x0, k);
63 const n = g.length / 2;
64 const out = new Float64Array(2 * n);
65 for (let i = 0; i < n; i++) {
66 // i k lam g = (-k lam g_im, k lam g_re)
67 out[2 * i] = -(gn[2 * i] - k * lam * g[2 * i + 1]);
68 out[2 * i + 1] = -(gn[2 * i + 1] + k * lam * g[2 * i]);
69 }
70 return out;
71}
73/** Impedance data of the field of a source outside the room (the self-test). */
74export function selftestRhs(xs: Float64Array, nrm: Float64Array, xe: ArrayLike<number>, k: number, lam: number): CVec {
75 const r = roomRhs(xs, nrm, xe, k, lam);
76 for (let i = 0; i < r.length; i++) r[i] = -r[i];
77 return r;
78}
80/** n deterministic pseudo-random points in the cube of half-width 0.8 a. */
81export function selftestPoints(n: number, a: number): Float64Array {
82 const out = new Float64Array(3 * n);
83 let s = 12345;
84 const rnd = () => {
85 s = (Math.imul(s, 1664525) + 1013904223) >>> 0;
86 return s / 4294967296;
87 };
88 for (let i = 0; i < 3 * n; i++) out[i] = 0.8 * a * (2 * rnd() - 1);
89 return out;
90}
92export interface Slice {
93 name: 'xy' | 'xz' | 'yz';
94 axes: [number, number];
95 fixed: number;
70e5eeeAdd the cymbal problem and GIF exportJeremy Magland 96 /** the value of the fixed coordinate */
97 value: number;
98 /** [lo, hi] along the two axes */
99 extent: [[number, number], [number, number]];
bc338a7WebGPU port of the fmm3dbie impedance solver as an interactive web appJeremy Magland 100 pts: Float64Array; // (3, ngrid^2), row-major over the grid with the first axis varying fastest
70e5eeeAdd the cymbal problem and GIF exportJeremy Magland 101 /** 1 where the field is computed (inside the room, or outside the cymbal) */
bc338a7WebGPU port of the fmm3dbie impedance solver as an interactive web appJeremy Magland 102 mask: Uint8Array;
103}
105/** Points on the three planes through x0; each grid runs from -a to a. */
106export function slicePoints(ngrid: number, x0: ArrayLike<number>, a: number, p: number): { grid: Float64Array; slices: Slice[] } {
107 const grid = new Float64Array(ngrid);
108 for (let i = 0; i < ngrid; i++) grid[i] = -a + (2 * a * i) / (ngrid - 1);
109 const defs: [Slice['name'], [number, number], number][] = [
110 ['xy', [0, 1], 2],
111 ['xz', [0, 2], 1],
112 ['yz', [1, 2], 0],
113 ];
114 const slices: Slice[] = [];
115 for (const [name, axes, fixed] of defs) {
116 const N = ngrid * ngrid;
117 const pts = new Float64Array(3 * N);
118 const mask = new Uint8Array(N);
119 for (let row = 0; row < ngrid; row++) {
120 for (let col = 0; col < ngrid; col++) {
121 const q = row * ngrid + col;
122 const x = [0, 0, 0];
123 x[axes[0]] = grid[col];
124 x[axes[1]] = grid[row];
125 x[fixed] = x0[fixed];
126 pts[q] = x[0];
127 pts[N + q] = x[1];
128 pts[2 * N + q] = x[2];
129 mask[q] = inside(x, a, p, 1e-3) ? 1 : 0;
130 }
131 }
70e5eeeAdd the cymbal problem and GIF exportJeremy Magland 132 slices.push({ name, axes, fixed, value: x0[fixed], extent: [[-a, a], [-a, a]], pts, mask });
134 return { grid, slices };
135}
137/** Power radiated by the source, k/(4 pi) + Im v(x0), and absorbed by the walls, -k Re(lambda) int |u|^2 dS. */
138export function energyBalance(k: number, lam: number, vAtSource: [number, number], wts: Float64Array, uWall: CVec) {
139 const pSrc = k / (4 * Math.PI) + vAtSource[1];
140 let s = 0;
141 for (let i = 0; i < wts.length; i++) s += wts[i] * (uWall[2 * i] ** 2 + uWall[2 * i + 1] ** 2);
142 const pWall = -k * lam * s;
143 return { pSrc, pWall, mismatch: Math.abs(pSrc - pWall) / Math.abs(pSrc) };
144}
146/** Relative 2-norm difference of two complex vectors. */
147export function relError(u: CVec, uex: CVec): number {
148 let num = 0;
149 let den = 0;
150 for (let i = 0; i < u.length; i++) {
151 num += (u[i] - uex[i]) ** 2;
152 den += uex[i] ** 2;
153 }
154 return Math.sqrt(num / den);
155}
157/** Select the columns of x (3, n) where mask is set. */
158export function maskedPoints(x: Float64Array, mask: Uint8Array): { pts: Float64Array; index: Int32Array } {
159 const n = mask.length;
160 let cnt = 0;
161 for (let i = 0; i < n; i++) cnt += mask[i];
162 const pts = new Float64Array(3 * cnt);
163 const index = new Int32Array(cnt);
164 let q = 0;
165 for (let i = 0; i < n; i++) {
166 if (!mask[i]) continue;
167 pts[q] = x[i];
168 pts[cnt + q] = x[n + i];
169 pts[2 * cnt + q] = x[2 * n + i];
170 index[q++] = i;
171 }
172 return { pts, index };
173}