2 * The discretized operator of fmm3dbie's right-preconditioned combined field
3 * impedance equation (helm_rpcomb_imp), on WebGPU. A port of
4 * fmm3dbie_gpu/solver.py from the CUDA proof of concept.
5 *
6 * u = S_k[sigma] + i alpha D_k S_kappa[sigma], kappa = i|k|,
7 *
8 * with the boundary condition du/dn + i k lambda u = f. The discretized
9 * equation is
10 *
11 * z sigma + S'_k sigma + i alpha S'_kappa phi2 + i alpha (D'_k - D'_kappa) phi1
12 * + i k lambda [S_k sigma + i alpha D_k phi1 + i alpha w phi1] = f,
13 *
14 * phi1 = S_kappa sigma, phi2 = S'_kappa sigma,
15 * z = (-1)^ifinout / 2 - i alpha / 4, w = 1/2 (exterior), -1/2 (interior).
16 *
17 * Every layer potential is the sum of a far part (smooth rule on the
18 * oversampled nodes of the patches that are not near the target, the dense
19 * kernel) and a near part (precomputed weights for the near patches, the
20 * sparse kernel). The near weights are computed once, in chunks of pairs, and
21 * kept on the GPU. GMRES runs on the CPU in double precision and calls
22 * applyK for each Krylov vector.
23 */
24import { leggauss } from '../gauss.ts';
25import { oversampleMatrix, type Surface } from '../surface.ts';
26import type { CVec } from '../problem.ts';
27import { gmres, type GmresResult } from '../gmres.ts';
28import { Gpu, Pacer, UniformWriter } from './device.ts';
29import { BUILD_LEAVES, DENSE_APPLY, DENSE_WG, NEAR_WEIGHTS, NT, SPARSE_APPLY } from './shaders.ts';
31export interface OperatorOptions {
32 /** near radius factor: a patch is near a target within rfac * (patch radius) of its centroid */
33 rfac?: number;
34 rfacQ?: number;
35 nq?: number;
36 nqr?: number;
37 nqx?: number;
38 subLen?: number;
39 oscLen?: number;
40 kRMax?: number;
41 nlmax?: number;
42 nover?: number;
43 /** seconds of GPU work per dispatch */
44 pace?: number;
45 progress?: (stage: string, fraction: number) => void;
46}
48export interface NearPairs {
49 pairT: Int32Array;
50 pairP: Int32Array;
51 rowPtr: Uint32Array;
53 * per pair, 8 floats: (u_t, v_t, kind, 0, e_x, e_y, e_z, 0). kind 0: the displacement to
54 * the target is a plain difference (off-surface targets, or on-surface targets that do
55 * not map into the patch's chart); 1: an on-surface target with preimage (u_t, v_t) in
56 * the patch's chart and gap e = chart(u_t, v_t) - x_t; 2: the target's own patch.
57 */
bc338a7WebGPU port of the fmm3dbie impedance solver as an interactive web appJeremy Magland 58 pairF: Float32Array;
59 npairs: number;
60}
62/** A set of targets with everything needed to apply the operator at them. */
63interface TargetSet {
64 n: number;
65 targ: GPUBuffer; // 2 vec4 per target: [x, 0], [n, 0]
66 pairs: NearPairs;
67 pairsI: GPUBuffer;
68 pairsF: GPUBuffer;
69 rowPtr: GPUBuffer;
70 pairP: GPUBuffer;
71 W: GPUBuffer;
72 nleaves: number;
73}
75export interface SetupInfo {
76 npts: number;
77 noversampled: number;
78 npairs: number;
79 nleaves: number;
80 weightsBytes: number;
81 setupSeconds: number;
82}
84const now = () => performance.now() / 1000;
86export class Operator {
87 readonly k: number;
88 readonly alpha: number;
89 /** wall impedance; the near weights do not depend on it, so it may be changed between solves */
90 lam: number;
91 readonly ifinout: number;
92 readonly norder: number;
93 readonly npols: number;
94 readonly npts: number;
95 readonly nover: number;
96 readonly npolsO: number;
97 readonly rfac: number;
98 readonly quad: { nq: number; nqr: number; nqx: number; subLen: number; oscLen: number };
99 readonly leaf: { rfacQ: number; nlmax: number; kRMax: number };
100 readonly progress: (stage: string, fraction: number) => void;
101 readonly pace: number;
102 info!: SetupInfo;
104 private gpu: Gpu;
105 private surf: Surface;
106 private r2: Float64Array;
107 private overMat: Float64Array;
108 private gpar!: GPUBuffer;
109 private tables!: GPUBuffer;
110 private srcgeo!: GPUBuffer;
111 private srcq!: GPUBuffer;
112 private qbuf!: GPUBuffer;
113 private out!: GPUBuffer;
114 private wall!: TargetSet;
115 private nself = 0;
116 private selfItems: Float32Array[] = []; // per surface node, its polar segments (8 floats each)
117 private uniformLeaves: GPUBuffer;
118 private uniformWeights: GPUBuffer;
119 private uniformDense: GPUBuffer;
120 private uniformSparse: GPUBuffer;
121 private destroyed = false;
123 constructor(gpu: Gpu, surf: Surface, k: number, alpha: number, lam: number, ifinout = 0, opts: OperatorOptions = {}) {
124 this.gpu = gpu;
125 this.surf = surf;
126 this.k = k;
127 this.alpha = alpha;
128 this.lam = lam;
129 this.ifinout = ifinout;
130 this.norder = surf.norder;
131 this.npols = surf.npols;
132 this.npts = surf.npts;
133 if (this.norder > 10) throw new Error('order at most 10');
134 this.rfac = opts.rfac ?? (this.norder <= 6 ? 2.0 : 1.25); // fmm3dbie's get_rfacs
135 this.nover = opts.nover ?? this.norder + 2;
136 this.npolsO = (this.nover + 1) ** 2;
137 let radMax = 0;
138 for (let j = 0; j < surf.npatches; j++) radMax = Math.max(radMax, surf.rads[j]);
139 // the radial rule must resolve the basis and exp(ikr) along a ray
140 const nqr = opts.nqr ?? Math.floor(Math.min(Math.max(12, this.norder + 2 * Math.abs(k) * radMax + 6), 32));
141 this.quad = { nq: opts.nq ?? 8, nqr, nqx: opts.nqx ?? 8, subLen: opts.subLen ?? 1.0, oscLen: opts.oscLen ?? 2.0 };
142 this.leaf = { rfacQ: opts.rfacQ ?? 2.0, nlmax: opts.nlmax ?? 14, kRMax: opts.kRMax ?? 3.0 };
143 this.progress = opts.progress ?? (() => {});
144 this.pace = opts.pace ?? 0.15;
145 this.r2 = new Float64Array(surf.npatches);
146 for (let j = 0; j < surf.npatches; j++) this.r2[j] = (this.rfac * surf.rads[j]) ** 2;
147 this.overMat = oversampleMatrix(this.norder, this.nover);
148 this.uniformLeaves = gpu.uniform(64);
149 this.uniformWeights = gpu.uniform(64);
150 this.uniformDense = gpu.uniform(64);
151 this.uniformSparse = gpu.uniform(64);
152 }
154 /** Upload the geometry and compute the near weights for the surface nodes. */
155 async setup(): Promise<SetupInfo> {
156 const t0 = now();
157 const { gpu, surf, npols, npts } = this;
159 if (surf.exact) {
160 // superellipsoid face charts: 4 vec4 per patch, [n, s0], [e1, s1], [e2, t0], [0, t1]
161 gp = new Float32Array(16 * surf.npatches);
162 for (let j = 0; j < surf.npatches; j++) {
163 const par = surf.exact.params.subarray(13 * j, 13 * j + 13);
164 gp.set([par[0], par[1], par[2], par[9], par[3], par[4], par[5], par[10], par[6], par[7], par[8], par[11], 0, 0, 0, par[12]], 16 * j);
165 }
166 } else {
167 // Legendre coefficients of xyz, r_u, r_v: 9 * npols floats per patch, padded to vec4
168 gp = new Float32Array(4 * Math.ceil((9 * npts) / 4));
169 for (let j = 0; j < surf.npatches; j++) {
170 for (let c = 0; c < 9; c++) {
171 for (let n = 0; n < npols; n++) gp[j * 9 * npols + c * npols + n] = surf.srccoefs[c * npts + j * npols + n];
172 }
173 }
175 this.gpar = gpu.storage(gp);
176 // umat and the quadrature rules in one table
177 const { nq, nqr, nqx } = this.quad;
178 const gl = leggauss(nq);
179 const glr = leggauss(nqr);
180 const glx = leggauss(nqx);
181 const tab = new Float32Array(npols * npols + 2 * nq + 2 * nqr + 2 * nqx);
182 tab.set(surf.umat, 0);
183 let o = npols * npols;
184 tab.set(gl.x, o);
185 tab.set(gl.w, o + nq);
186 o += 2 * nq;
187 for (let i = 0; i < nqr; i++) {
188 tab[o + i] = 0.5 * (glr.x[i] + 1);
189 tab[o + nqr + i] = 0.5 * glr.w[i];
190 }
191 o += 2 * nqr;
192 tab.set(glx.x, o);
193 tab.set(glx.w, o + nqx);
194 this.tables = gpu.storage(tab);
195 // oversampled sources for the far field
196 const ov = surf.oversampled(this.nover);
197 const N = ov.w.length;
198 const sg = new Float32Array(8 * N);
199 for (let i = 0; i < N; i++) {
200 sg[8 * i] = ov.x[i];
201 sg[8 * i + 1] = ov.x[N + i];
202 sg[8 * i + 2] = ov.x[2 * N + i];
203 sg[8 * i + 3] = ov.w[i];
204 sg[8 * i + 4] = ov.n[i];
205 sg[8 * i + 5] = ov.n[N + i];
206 sg[8 * i + 6] = ov.n[2 * N + i];
207 }
208 this.srcgeo = gpu.storage(sg);
209 this.srcq = gpu.storage(16 * N);
210 this.qbuf = gpu.storage(2 * 8 * npts);
211 this.out = gpu.storage(6 * 8 * npts);
213 // self-patch polar segments for every node, from the first fundamental form at the node
214 this.selfItems = new Array(npts);
215 for (let i = 0; i < npts; i++) this.selfItems[i] = this.selfSegments(i);
217 const xs = surf.srcvals.subarray(0, 3 * npts);
218 const nrm = surf.srcvals.subarray(9 * npts, 12 * npts);
219 this.wall = await this.makeTargets(xs, nrm, true, 'near quadrature');
220 await gpu.sync();
221 const setupSeconds = now() - t0;
222 this.info = {
223 npts,
224 noversampled: N,
225 npairs: this.wall.pairs.npairs,
226 nleaves: this.wall.nleaves,
227 weightsBytes: this.wall.W.size,
228 setupSeconds,
229 };
230 return this.info;
231 }
233 destroy() {
234 this.destroyed = true;
235 for (const b of [this.gpar, this.tables, this.srcgeo, this.srcq, this.qbuf, this.out, this.uniformLeaves, this.uniformWeights, this.uniformDense, this.uniformSparse]) b?.destroy();
236 if (this.wall) this.freeTargets(this.wall);
237 }
239 /** Bytes of near weights for npairs pairs. */
240 weightsBytes(npairs: number): number {
241 return npairs * 6 * this.npols * 8;
242 }
244 /** Estimated bytes of near weights for the surface nodes, before setup. */
245 estimateWeightsBytes(): number {
246 const { npts } = this;
247 const xs = this.surf.srcvals.subarray(0, 3 * npts);
248 let cnt = 0;
249 const { npatches } = this.surf;
250 const cen = this.surf.centroids;
251 // sample every 16th node
252 for (let i = 0; i < npts; i += 16) {
253 for (let j = 0; j < npatches; j++) {
254 const d2 = (xs[i] - cen[j]) ** 2 + (xs[npts + i] - cen[npatches + j]) ** 2 + (xs[2 * npts + i] - cen[2 * npatches + j]) ** 2;
255 if (d2 <= this.r2[j]) cnt++;
256 }
257 }
258 return this.weightsBytes(Math.ceil((cnt * npts) / Math.ceil(npts / 16)));
259 }
261 // ---------------------------------------------------------------- near pairs
263 /** Pairs (target, patch) with |x_t - c_j|^2 <= r2_j, grouped by target, patches ascending. */
264 nearPairs(x: Float64Array, onSurface: boolean): NearPairs {
265 const n = x.length / 3;
266 const { npatches, centroids: cen, patchId, uvNodes } = this.surf;
267 const r2 = this.r2;
268 const tList: number[] = [];
269 const pList: number[] = [];
270 const fList: number[] = [];
271 const rowPtr = new Uint32Array(n + 1);
bc338a7WebGPU port of the fmm3dbie impedance solver as an interactive web appJeremy Magland 273 for (let i = 0; i < n; i++) {
274 const xi = x[i];
275 const yi = x[n + i];
276 const zi = x[2 * n + i];
277 for (let j = 0; j < npatches; j++) {
278 const d2 = (xi - cen[j]) ** 2 + (yi - cen[npatches + j]) ** 2 + (zi - cen[2 * npatches + j]) ** 2;
279 if (d2 > r2[j]) continue;
280 tList.push(i);
281 pList.push(j);
282 if (!onSurface) {
bc338a7WebGPU port of the fmm3dbie impedance solver as an interactive web appJeremy Magland 284 } else if (patchId[i] === j) {
70e5eeeAdd the cymbal problem and GIF exportJeremy Magland 285 fList.push(uvNodes[i], uvNodes[n + i], 2, 0, 0, 0, 0, 0);
286 } else if (exact) {
bc338a7WebGPU port of the fmm3dbie impedance solver as an interactive web appJeremy Magland 287 // the target mapped into this patch's face chart: c_t = x_t / (x_t . n)
70e5eeeAdd the cymbal problem and GIF exportJeremy Magland 288 const par = exact.params.subarray(13 * j, 13 * j + 13);
bc338a7WebGPU port of the fmm3dbie impedance solver as an interactive web appJeremy Magland 289 const xn = xi * par[0] + yi * par[1] + zi * par[2];
290 if (xn < 0.2) {
bc338a7WebGPU port of the fmm3dbie impedance solver as an interactive web appJeremy Magland 292 } else {
293 const cx = xi / xn;
294 const cy = yi / xn;
295 const cz = zi / xn;
296 const st = cx * par[3] + cy * par[4] + cz * par[5];
297 const tt = cx * par[6] + cy * par[7] + cz * par[8];
70e5eeeAdd the cymbal problem and GIF exportJeremy Magland 298 fList.push((2 * (st - par[9])) / (par[10] - par[9]) - 1, (2 * (tt - par[11])) / (par[12] - par[11]) - 1, 1, 0, 0, 0, 0, 0);
301 const pre = this.preimage(j, xi, yi, zi);
302 if (pre) fList.push(pre[0], pre[1], 1, 0, pre[2], pre[3], pre[4], 0);
303 else fList.push(0, 0, 0, 0, 0, 0, 0, 0);
305 }
306 rowPtr[i + 1] = tList.length;
307 }
308 return { pairT: Int32Array.from(tList), pairP: Int32Array.from(pList), rowPtr, pairF: Float32Array.from(fList), npairs: tList.length };
309 }
312 * For the interpolated geometry: the preimage (u, v) of the point x in the Legendre
313 * interpolant of patch j, extended beyond [-1, 1]^2, by Gauss-Newton from the nearest
314 * node, and the gap e = x(u, v) - x. Returns [u, v, e_x, e_y, e_z], or null when the
315 * point does not map close to the patch (then a plain difference is used).
316 */
317 private preimage(j: number, x: number, y: number, z: number): number[] | null {
318 const { npols, npts, norder } = this;
319 const n1 = norder + 1;
320 const sc = this.surf.srccoefs;
321 const sv = this.surf.srcvals;
322 const base = j * npols;
323 let best = 0;
324 let bd = Infinity;
325 for (let l = 0; l < npols; l++) {
326 const d = (sv[base + l] - x) ** 2 + (sv[npts + base + l] - y) ** 2 + (sv[2 * npts + base + l] - z) ** 2;
327 if (d < bd) {
328 bd = d;
329 best = l;
330 }
331 }
332 let u = this.surf.uvNodes[base + best];
333 let v = this.surf.uvNodes[npts + base + best];
334 const Pu = new Float64Array(n1);
335 const Pv = new Float64Array(n1);
336 const Du = new Float64Array(n1);
337 const Dv = new Float64Array(n1);
338 const leg = (t: number, P: Float64Array, Dd: Float64Array) => {
339 P[0] = 1;
340 Dd[0] = 0;
341 if (norder >= 1) {
342 P[1] = t;
343 Dd[1] = 1;
344 }
345 for (let m = 1; m < norder; m++) {
346 P[m + 1] = ((2 * m + 1) * t * P[m] - m * P[m - 1]) / (m + 1);
347 Dd[m + 1] = Dd[m - 1] + (2 * m + 1) * P[m];
348 }
349 };
350 const g = new Float64Array(9); // x, x_u, x_v of the xyz interpolant
351 const evalAt = () => {
352 leg(u, Pu, Du);
353 leg(v, Pv, Dv);
354 g.fill(0);
355 for (let c = 0; c < 3; c++) {
356 for (let b = 0; b < n1; b++) {
357 for (let a = 0; a < n1; a++) {
358 const cf = sc[c * npts + base + b * n1 + a];
359 g[c] += cf * Pu[a] * Pv[b];
360 g[3 + c] += cf * Du[a] * Pv[b];
361 g[6 + c] += cf * Pu[a] * Dv[b];
362 }
363 }
364 }
365 };
366 for (let it = 0; it < 30; it++) {
367 evalAt();
368 const r = [x - g[0], y - g[1], z - g[2]];
369 const a11 = g[3] ** 2 + g[4] ** 2 + g[5] ** 2;
370 const a12 = g[3] * g[6] + g[4] * g[7] + g[5] * g[8];
371 const a22 = g[6] ** 2 + g[7] ** 2 + g[8] ** 2;
372 const b1 = g[3] * r[0] + g[4] * r[1] + g[5] * r[2];
373 const b2 = g[6] * r[0] + g[7] * r[1] + g[8] * r[2];
374 const det = a11 * a22 - a12 * a12;
375 if (!(det > 0)) return null;
376 const du = (a22 * b1 - a12 * b2) / det;
377 const dv = (a11 * b2 - a12 * b1) / det;
378 u += du;
379 v += dv;
380 if (Math.abs(u) > 3 || Math.abs(v) > 3) return null;
381 if (Math.abs(du) + Math.abs(dv) < 1e-13) break;
382 }
383 if (Math.abs(u) > 2 || Math.abs(v) > 2) return null;
384 evalAt();
385 const e = [g[0] - x, g[1] - y, g[2] - z];
386 if (Math.hypot(e[0], e[1], e[2]) > 0.1 * this.surf.rads[j]) return null;
387 return [u, v, ...e];
388 }
bc338a7WebGPU port of the fmm3dbie impedance solver as an interactive web appJeremy Magland 390 /** The sinh-substituted segments along the four edges for the self patch of node i: 8 floats each. */
391 private selfSegments(i: number): Float32Array {
392 const { npts } = this;
393 const sv = this.surf.srcvals;
394 const u0 = this.surf.uvNodes[i];
395 const v0 = this.surf.uvNodes[npts + i];
396 const ru = [sv[3 * npts + i], sv[4 * npts + i], sv[5 * npts + i]];
397 const rv = [sv[6 * npts + i], sv[7 * npts + i], sv[8 * npts + i]];
398 const E = ru[0] ** 2 + ru[1] ** 2 + ru[2] ** 2;
399 const F = ru[0] * rv[0] + ru[1] * rv[1] + ru[2] * rv[2];
400 const G = rv[0] ** 2 + rv[1] ** 2 + rv[2] ** 2;
401 const det = Math.sqrt(Math.max(E * G - F * F, 1e-300));
402 const k = Math.abs(this.k);
403 const { subLen, oscLen } = this.quad;
404 const items: number[] = [];
405 for (let e = 0; e < 4; e++) {
406 const h = [1 - u0, -1 - u0, 1 - v0, -1 - v0][e];
407 if (Math.abs(h) < 1e-12) continue;
408 let x1: number, x2: number, xc: number, heff: number, arc: number;
409 if (e < 2) {
410 x1 = -1 - v0;
411 x2 = 1 - v0;
412 xc = (-h * F) / G;
413 heff = (Math.abs(h) * det) / G;
414 arc = Math.sqrt(G);
415 } else {
416 x1 = -1 - u0;
417 x2 = 1 - u0;
418 xc = (-h * F) / E;
419 heff = (Math.abs(h) * det) / E;
420 arc = Math.sqrt(E);
421 }
422 const xi0 = Math.asinh((x1 - xc) / heff);
423 const xi1 = Math.asinh((x2 - xc) / heff);
424 // sub-intervals of length <= subLen in xi whose extent along the edge is <= dxmax
425 const dxmax = k * arc > 0 ? oscLen / (k * arc) : Infinity;
426 let xi = xi0;
427 let cnt = 0;
428 while (xi < xi1 && cnt < 256) {
429 const xh = xc + heff * Math.sinh(xi);
430 let xn = Math.min(xi + subLen, Math.asinh((xh + dxmax - xc) / heff));
431 if (xn > xi1 || cnt === 255) xn = xi1;
432 items.push(1, e, xi, xn, h, xc, heff, 0);
433 xi = xn;
434 cnt++;
435 }
436 }
437 return Float32Array.from(items);
438 }
440 // ------------------------------------------------------------ near weights
442 /** Build the near pairs, quadtree leaves and near weights for a set of targets. */
443 private async makeTargets(x: Float64Array, nrm: Float64Array | null, onSurface: boolean, stage: string): Promise<TargetSet> {
444 const { gpu } = this;
445 const n = x.length / 3;
446 const tg = new Float32Array(8 * n);
447 for (let i = 0; i < n; i++) {
448 tg[8 * i] = x[i];
449 tg[8 * i + 1] = x[n + i];
450 tg[8 * i + 2] = x[2 * n + i];
451 if (nrm) {
452 tg[8 * i + 4] = nrm[i];
453 tg[8 * i + 5] = nrm[n + i];
454 tg[8 * i + 6] = nrm[2 * n + i];
455 }
456 }
457 const pairs = this.nearPairs(x, onSurface);
458 const wbytes = this.weightsBytes(pairs.npairs);
459 const lim = gpu.device.limits;
460 if (wbytes > lim.maxStorageBufferBindingSize || wbytes > lim.maxBufferSize) {
461 throw new Error(
462 `the near weights need ${(wbytes / 2 ** 20).toFixed(0)} MB but this GPU allows buffers of at most ` +
463 `${(Math.min(lim.maxStorageBufferBindingSize, lim.maxBufferSize) / 2 ** 20).toFixed(0)} MB; use fewer patches or a lower order`,
464 );
465 }
466 const pi = new Int32Array(2 * Math.max(1, pairs.npairs));
467 for (let p = 0; p < pairs.npairs; p++) {
468 pi[2 * p] = pairs.pairT[p];
469 pi[2 * p + 1] = pairs.pairP[p];
470 }
471 const set: TargetSet = {
472 n,
473 targ: gpu.storage(tg),
474 pairs,
475 pairsI: gpu.storage(pi),
70e5eeeAdd the cymbal problem and GIF exportJeremy Magland 476 pairsF: gpu.storage(pairs.npairs ? pairs.pairF : new Float32Array(8)),
bc338a7WebGPU port of the fmm3dbie impedance solver as an interactive web appJeremy Magland 477 rowPtr: gpu.storage(pairs.rowPtr),
478 pairP: gpu.storage(pairs.npairs ? new Uint32Array(pairs.pairP.buffer, pairs.pairP.byteOffset, pairs.npairs) : new Uint32Array(1)),
479 W: gpu.storage(Math.max(16, wbytes)),
480 nleaves: 0,
481 };
482 if (pairs.npairs > 0) set.nleaves = await this.computeWeights(set, stage);
483 return set;
484 }
486 private freeTargets(s: TargetSet) {
487 for (const b of [s.targ, s.pairsI, s.pairsF, s.rowPtr, s.pairP, s.W]) b.destroy();
488 }
490 /** Quadtree leaves and near weights for all pairs of a target set, in paced chunks. Returns the leaf count. */
491 private async computeWeights(set: TargetSet, stage: string): Promise<number> {
492 const { gpu, quad, leaf } = this;
493 const { npairs, pairF } = set.pairs;
494 const leavesPipe = await gpu.pipeline('leaves', BUILD_LEAVES, 'build_leaves');
495 const weightsPipe = await gpu.pipeline('weights', NEAR_WEIGHTS, 'near_weights');
496 const pacer = new Pacer(this.pace, 512);
497 let nleaves = 0;
498 let p0 = 0;
499 let itemsPerPair = 40;
500 const maxItemsBytes = 48 * 2 ** 20;
501 while (p0 < npairs) {
502 if (this.destroyed) throw new Error('operator destroyed');
503 const cap = Math.min(65535, npairs - p0, Math.max(256, Math.floor(maxItemsBytes / (32 * itemsPerPair))));
504 const n = pacer.next(cap);
505 const t0 = now();
506 // pass 1: leaf counts
507 const counts = gpu.storage(4 * n);
509 const ga = ex ? ex.a : 0;
510 const gpe = ex ? ex.p : 0;
511 const gmode = ex ? 0 : 1;
512 const uw = new UniformWriter(11);
513 uw.u32(n).u32(p0).u32(0).i32(leaf.nlmax).f32(leaf.rfacQ).f32(Math.abs(this.k)).f32(leaf.kRMax).f32(ga).i32(gpe).i32(gmode).i32(this.norder);
bc338a7WebGPU port of the fmm3dbie impedance solver as an interactive web appJeremy Magland 514 gpu.write(this.uniformLeaves, new Uint8Array(uw.bytes));
515 const dummyItems = gpu.storage(16);
516 const dummyOffsets = gpu.storage(16);
517 gpu.dispatch(leavesPipe, gpu.bind(leavesPipe, [this.uniformLeaves, set.pairsI, set.pairsF, set.targ, this.gpar, counts, dummyOffsets, dummyItems]), Math.ceil(n / 64));
518 const cnt = new Uint32Array(await gpu.read(counts, 4 * n));
519 // self pairs: the CPU segments
520 const ptr = new Uint32Array(n + 1);
521 for (let q = 0; q < n; q++) {
522 const p = p0 + q;
70e5eeeAdd the cymbal problem and GIF exportJeremy Magland 523 const c = pairF[8 * p + 2] === 2 ? this.selfItems[set.pairs.pairT[p]].length / 8 : cnt[q];
524 if (pairF[8 * p + 2] !== 2) nleaves += cnt[q];
bc338a7WebGPU port of the fmm3dbie impedance solver as an interactive web appJeremy Magland 525 ptr[q + 1] = ptr[q] + c;
526 }
527 const nitems = ptr[n];
528 const items = new Float32Array(8 * Math.max(1, nitems));
529 for (let q = 0; q < n; q++) {
530 const p = p0 + q;
70e5eeeAdd the cymbal problem and GIF exportJeremy Magland 531 if (pairF[8 * p + 2] === 2) items.set(this.selfItems[set.pairs.pairT[p]], 8 * ptr[q]);
533 const itemsBuf = gpu.storage(items);
534 const ptrBuf = gpu.storage(ptr);
535 // pass 2: fill the leaves
537 uw2.u32(n).u32(p0).u32(1).i32(leaf.nlmax).f32(leaf.rfacQ).f32(Math.abs(this.k)).f32(leaf.kRMax).f32(ga).i32(gpe).i32(gmode).i32(this.norder);
bc338a7WebGPU port of the fmm3dbie impedance solver as an interactive web appJeremy Magland 538 gpu.write(this.uniformLeaves, new Uint8Array(uw2.bytes));
539 gpu.dispatch(leavesPipe, gpu.bind(leavesPipe, [this.uniformLeaves, set.pairsI, set.pairsF, set.targ, this.gpar, counts, ptrBuf, itemsBuf]), Math.ceil(n / 64));
540 // the weights
542 uw3.u32(n).u32(p0).i32(this.norder).i32(this.npols).i32(quad.nq).i32(quad.nqr).i32(quad.nqx).f32(this.k).f32(ga).i32(gpe).i32(gmode);
bc338a7WebGPU port of the fmm3dbie impedance solver as an interactive web appJeremy Magland 543 gpu.write(this.uniformWeights, new Uint8Array(uw3.bytes));
544 gpu.dispatch(weightsPipe, gpu.bind(weightsPipe, [this.uniformWeights, set.pairsI, set.pairsF, set.targ, this.gpar, ptrBuf, itemsBuf, this.tables, set.W]), n);
545 await gpu.sync();
546 const dt = now() - t0;
547 pacer.record(n, dt);
548 itemsPerPair = Math.max(1, nitems / n);
549 for (const b of [counts, dummyItems, dummyOffsets, itemsBuf, ptrBuf]) b.destroy();
550 p0 += n;
551 this.progress(stage, p0 / npairs);
552 }
553 return nleaves;
554 }
556 // ------------------------------------------------------------- application
558 /** Values of q at the oversampled nodes, per patch: q_o = M q. */
559 private oversample(q: CVec, out: Float32Array, offset: number) {
560 const { npols, npolsO, overMat } = this;
561 const np = this.surf.npatches;
562 for (let j = 0; j < np; j++) {
563 for (let kk = 0; kk < npolsO; kk++) {
564 let sr = 0;
565 let si = 0;
566 for (let m = 0; m < npols; m++) {
567 const M = overMat[kk * npols + m];
568 sr += M * q[2 * (j * npols + m)];
569 si += M * q[2 * (j * npols + m) + 1];
570 }
571 const idx = 4 * (j * npolsO + kk) + offset;
572 out[idx] = sr;
573 out[idx + 1] = si;
574 }
575 }
576 }
578 /** Far and near parts of the selected kernels at a target set. sel: 6 entries in {0, 1, 2}; reqs: [kernel, density] pairs. */
579 private async applyAt(set: TargetSet, qa: CVec, qb: CVec | null, sel: number[], reqs: [number, number][]): Promise<Float32Array> {
580 const { gpu, npts, npolsO } = this;
581 const N = npolsO * this.surf.npatches;
582 const sq = new Float32Array(4 * N);
583 this.oversample(qa, sq, 0);
584 if (qb) this.oversample(qb, sq, 2);
585 gpu.write(this.srcq, sq);
586 const qd = new Float32Array(2 * 2 * npts);
587 for (let i = 0; i < 2 * npts; i++) qd[i] = qa[i];
588 if (qb) for (let i = 0; i < 2 * npts; i++) qd[2 * npts + i] = qb[i];
589 gpu.write(this.qbuf, qd);
590 const out = set.n === this.wall.n ? this.out : gpu.storage(6 * 8 * set.n);
591 gpu.clear(out);
592 const selmask = sel.reduce((acc, s, m) => acc | (s << (4 * m)), 0);
593 const densePipe = await gpu.pipeline('dense', DENSE_APPLY, 'dense_apply');
594 const group = gpu.bind(densePipe, [this.uniformDense, set.targ, this.srcgeo, this.srcq, set.rowPtr, set.pairP, out]);
595 const pacer = new Pacer(this.pace, Math.max(1, Math.floor(4e8 / (set.n * npolsO * 60))));
596 const np = this.surf.npatches;
597 let j0 = 0;
598 while (j0 < np) {
599 const n = pacer.next(np - j0);
600 const t0 = now();
601 const uw = new UniformWriter(6);
602 uw.u32(set.n).u32(j0).u32(j0 + n).u32(npolsO).f32(this.k).u32(selmask);
603 gpu.write(this.uniformDense, new Uint8Array(uw.bytes));
604 gpu.dispatch(densePipe, group, Math.ceil(set.n / DENSE_WG));
605 await gpu.sync();
606 pacer.record(n, now() - t0);
607 j0 += n;
608 }
609 if (set.pairs.npairs > 0 && reqs.length > 0) {
610 const sparsePipe = await gpu.pipeline('sparse', SPARSE_APPLY, 'sparse_apply');
611 const uw = new UniformWriter(8);
612 uw.u32(set.n).u32(this.npols).u32(npts).u32(reqs.length);
613 for (let r = 0; r < 4; r++) uw.u32(r < reqs.length ? reqs[r][0] | (reqs[r][1] << 8) : 0);
614 gpu.write(this.uniformSparse, new Uint8Array(uw.bytes));
615 gpu.dispatch(sparsePipe, gpu.bind(sparsePipe, [this.uniformSparse, set.rowPtr, set.pairP, set.W, this.qbuf, out]), Math.ceil(set.n / 128));
616 }
617 const res = new Float32Array(await gpu.read(out, 6 * 8 * set.n));
618 if (out !== this.out) out.destroy();
619 return res;
620 }
622 zid(): [number, number] {
623 return [(-1) ** this.ifinout / 2, -this.alpha / 4];
624 }
626 /** The operator without the identity term z. */
627 async applyK(sigma: CVec): Promise<CVec> {
628 const { npts, k, alpha: a, lam } = this;
629 const F = await this.applyAt(this.wall, sigma, null, [1, 1, 1, 0, 1, 0], [[0, 0], [1, 0], [2, 0], [4, 0]]);
630 const phi1 = new Float64Array(2 * npts);
631 const phi2 = new Float64Array(2 * npts);
632 const res = new Float64Array(2 * npts);
633 for (let i = 0; i < npts; i++) {
634 phi1[2 * i] = F[2 * (npts + i)];
635 phi1[2 * i + 1] = F[2 * (npts + i) + 1];
636 phi2[2 * i] = F[2 * (2 * npts + i)];
637 phi2[2 * i + 1] = F[2 * (2 * npts + i) + 1];
638 const s4r = F[2 * (4 * npts + i)];
639 const s4i = F[2 * (4 * npts + i) + 1];
640 // res = F0 + i k lam F4
641 res[2 * i] = F[2 * i] - k * lam * s4i;
642 res[2 * i + 1] = F[2 * i + 1] + k * lam * s4r;
643 }
644 const G = await this.applyAt(this.wall, phi1, phi2, [0, 0, 2, 1, 0, 1], [[2, 1], [3, 0], [5, 0]]);
645 const sgn = (-1) ** this.ifinout;
646 for (let i = 0; i < npts; i++) {
647 const g2r = G[2 * (2 * npts + i)];
648 const g2i = G[2 * (2 * npts + i) + 1];
649 const g3r = G[2 * (3 * npts + i)];
650 const g3i = G[2 * (3 * npts + i) + 1];
651 const g5r = G[2 * (5 * npts + i)];
652 const g5i = G[2 * (5 * npts + i) + 1];
653 // res += i a (G2 + G3) - k a lam G5 + sgn k a lam phi1 / 2
654 res[2 * i] += -a * (g2i + g3i) - k * a * lam * g5r + (sgn * k * a * lam * phi1[2 * i]) / 2;
655 res[2 * i + 1] += a * (g2r + g3r) - k * a * lam * g5i + (sgn * k * a * lam * phi1[2 * i + 1]) / 2;
656 }
657 return res;
658 }
660 /** phi1 = S_kappa[sigma] on the surface. */
661 async sKappa(sigma: CVec): Promise<CVec> {
662 const F = await this.applyAt(this.wall, sigma, null, [0, 1, 0, 0, 0, 0], [[1, 0]]);
663 const { npts } = this;
664 const phi1 = new Float64Array(2 * npts);
665 for (let i = 0; i < 2 * npts; i++) phi1[i] = F[2 * npts + i];
666 return phi1;
667 }
669 /** One of the six operators applied to q on the surface (principal value), for tests. */
670 async applySingle(m: number, q: CVec): Promise<CVec> {
671 const sel = [0, 0, 0, 0, 0, 0];
672 sel[m] = 1;
673 const F = await this.applyAt(this.wall, q, null, sel, [[m, 0]]);
674 const { npts } = this;
675 const out = new Float64Array(2 * npts);
676 for (let i = 0; i < 2 * npts; i++) out[i] = F[2 * m * npts + i];
677 return out;
678 }
680 async solve(rhs: CVec, epsGmres = 1e-5, numit = 200, onIter?: (it: number, rres: number) => void): Promise<GmresResult & { phi1: CVec; seconds: number }> {
681 const t0 = now();
682 const r = await gmres((x) => this.applyK(x), rhs, this.surf.wts, this.zid(), epsGmres, numit, onIter);
683 const phi1 = await this.sKappa(r.x);
684 return { ...r, phi1, seconds: now() - t0 };
685 }
687 /** u = S_k sigma + i alpha D_k phi1 at the surface nodes, as the limit from the side selected by ifinout. */
688 async evalOnSurface(sigma: CVec, phi1: CVec): Promise<CVec> {
689 const F = await this.applyAt(this.wall, sigma, phi1, [0, 0, 0, 0, 1, 2], [[4, 0], [5, 1]]);
690 return this.combine(F, this.npts, phi1);
691 }
693 /** u = S_k sigma + i alpha D_k phi1 at off-surface targets x (3, n), in chunks. */
694 async eval(x: Float64Array, sigma: CVec, phi1: CVec, stage = 'field'): Promise<CVec> {
695 const n = x.length / 3;
696 const u = new Float64Array(2 * n);
697 const chunk = 4096;
698 for (let i0 = 0; i0 < n; i0 += chunk) {
699 const m = Math.min(chunk, n - i0);
700 const xs = new Float64Array(3 * m);
701 for (let d = 0; d < 3; d++) xs.set(x.subarray(d * n + i0, d * n + i0 + m), d * m);
702 const set = await this.makeTargets(xs, null, false, stage);
703 const F = await this.applyAt(set, sigma, phi1, [0, 0, 0, 0, 1, 2], [[4, 0], [5, 1]]);
704 this.freeTargets(set);
705 u.set(this.combine(F, m, null), 2 * i0);
706 this.progress(stage, (i0 + m) / n);
707 }
708 return u;
709 }
711 private combine(F: Float32Array, n: number, phi1OnSurface: CVec | null): CVec {
712 const a = this.alpha;
713 const u = new Float64Array(2 * n);
714 const jump = phi1OnSurface ? (a * (-1) ** (this.ifinout + 1)) / 2 : 0;
715 for (let i = 0; i < n; i++) {
716 const s4r = F[2 * (4 * n + i)];
717 const s4i = F[2 * (4 * n + i) + 1];
718 const d5r = F[2 * (5 * n + i)];
719 const d5i = F[2 * (5 * n + i) + 1];
720 // u = S + i a D (+ i a (-1)^(ifinout+1) phi1 / 2 on the surface)
721 u[2 * i] = s4r - a * d5i;
722 u[2 * i + 1] = s4i + a * d5r;
723 if (phi1OnSurface) {
724 u[2 * i] += -jump * phi1OnSurface[2 * i + 1];
725 u[2 * i + 1] += jump * phi1OnSurface[2 * i];
726 }
727 }
728 return u;
729 }
731 /** The near weights of the surface nodes, for tests: (npairs, 6, npols) complex interleaved. */
732 async readWallWeights(): Promise<{ pairs: NearPairs; W: Float32Array }> {
733 const W = new Float32Array(await this.gpu.read(this.wall.W, this.weightsBytes(this.wall.pairs.npairs)));
734 return { pairs: this.wall.pairs, W };
735 }
736}
738export { NT };