GPU implementation of fmm3dbie's interior Helmholtz impedance solver, with a showcase room problem
18 changed files+2111−0
.gitignoreadded+3−0View file
@@ -0,0 +1,3 @@
1+results/
2+__pycache__/
3+*.pyc
README.mdadded+89−0View file
@@ -0,0 +1,89 @@
1+# fmm3dbie on the GPU: a proof of concept
2+
3+fmm3dbie solves boundary integral equations on surfaces in three dimensions with high-order quadrature and the fast multipole method. On a CPU, most of its run time for the Helmholtz equation goes into computing the near-field quadrature corrections, and the rest into FMM-accelerated matrix-vector products inside GMRES. Here we present a GPU implementation of one fmm3dbie solver, the interior Helmholtz impedance problem (`helm_rpcomb_imp_solver`). It uses the same integral equation and the same surface discretization, together with a showcase problem that runs on both codes so that their speed and accuracy can be compared directly.
4+
5+The showcase is a point source at a fixed frequency in a closed room with partially absorbing walls. The room is a slightly rounded cube, the superellipsoid |x|⁸ + |y|⁸ + |z|⁸ = 1. With time dependence e^{-iωt}, the pressure amplitude satisfies Δu + k²u = −δ(x − x₀) inside the room and ∂u/∂n + ikλu = 0 on the walls, with outward normal n and λ = −β. Here β > 0 is the normalized wall admittance, and smaller β gives more reflective walls.
6+
7+
8+
9+A single precision port of the same solver to WebGPU, which runs in a browser without CUDA, is in [fmm3dbie-webgpu-poc](https://vault1.magland.org/proof-of-concept/fmm3dbie-webgpu-poc).
10+
11+## What is being compared
12+
13+Both backends solve for the correction v = u − G, where G is the free-space field of the source. The representation is v = S_k[σ] + iα D_k S_{i|k|}[σ], and the equation is fmm3dbie's discretized impedance equation. Both use the same 6m² quadrilateral patches (m × m on each face of the cube) of order p, with (p+1)² Gauss–Legendre nodes each, so the unknowns are the same.
14+
15+- `fmm3dbie` runs fmm3dbie itself through ctypes (`reference/`), with the quadrature and FMM tolerance `--eps` (default 1e-4).
16+- `gpu` runs `fmm3dbie_gpu/`, a CuPy program with CUDA kernels compiled at run time. It does not use fmm3dbie.
17+
18+The GPU code differs from fmm3dbie in four ways, which also explain its speed and accuracy.
19+
20+1. Self-patch quadrature. fmm3dbie uses precomputed generalized Gaussian rules. We use polar coordinates on the four triangles fanned from the target to the patch edges, with a sinh substitution along each edge. This handles targets close to an edge, needs no tables, and computes every (target, patch) pair in one thread block.
21+2. Near-patch quadrature. We use a quadtree on the reference square, refined by distance and by wavelength, with an 8 × 8 Gauss–Legendre rule per leaf. This is fmm3dbie's distance-based strategy rather than its adaptive one.
22+3. Far field. Patches that are not near a target are summed directly, O(N²), rather than with an FMM. Near patches are excluded from this sum, so it has no large cancelling terms and runs in single precision, with partial sums accumulated in double precision.
23+4. Geometry. For the rounded cube we evaluate the exact surface inside the integrals. fmm3dbie uses polynomial interpolants of the node positions and tangents. The interpolants are slightly inconsistent with each other, and the double-layer kernel amplifies this at nodes next to patch edges, by up to 0.6% of the local integral in our tests. With the exact geometry the GPU solution is more accurate than fmm3dbie's on the same nodes. `--interp-geometry` switches to the interpolated model for a like-for-like comparison.
24+
25+## Quick start
26+
27+The fmm3dbie backend needs gfortran, make and OpenBLAS; the GPU backend needs an NVIDIA GPU with a CUDA 12 driver. Both need Python 3.9 or later with NumPy, and the plots need Matplotlib.
28+
29+ scripts/build_fmm3dbie.sh # clones and builds ~/src/fmm3dbie (or set FMM3DBIE_LIB)
30+ scripts/setup_gpu.sh # CuPy and CUDA libraries as pip wheels, if no CUDA toolkit
31+ source scripts/gpuenv.sh
32+
33+ python3 run.py --backend gpu --case small --selftest
34+ python3 run.py --backend fmm3dbie --case small --selftest
35+ python3 compare.py
36+ python3 plot.py results/small-gpu
37+
38+Each run writes `results/<case>-<backend>/summary.json`, with timings, iteration counts, errors and machine information, and `field.npz`. `compare.py` tabulates all runs and, for every case run with both backends, reports the speedup and the relative difference of the two wall fields. It also writes `results/comparison.md`. The runs do not have to be on the same machine: copy the result directories together before comparing. For a fair speedup, though, run fmm3dbie on the same machine as the GPU and with all of its cores (`OMP_NUM_THREADS`).
39+
40+The cases are
41+
42+| case | k | wall admittance β | room width | patches | order | points |
43+|---|---|---|---|---|---|---|
44+| small | 10 | 0.2 | 3.2 wavelengths | 96 | 6 | 4,704 |
45+| medium | 10 | 0.2 | 3.2 wavelengths | 216 | 8 | 17,496 |
46+| large | 40 | 0.05 | 12.7 wavelengths | 600 | 8 | 48,600 |
47+
48+and `--k`, `--beta`, `--m` and `--norder` override them. Run `python3 run.py --help` for the other options.
49+
50+## What each run reports
51+
52+- Timings. For the GPU, *setup* is the near quadrature for the operator and *solve* is GMRES. The one-time loading of the compiled CUDA kernels is reported separately and excluded. fmm3dbie computes its near quadrature inside the solver call, so its setup is zero and its solve includes it. Compare *setup + solve*.
53+- Self-test (`--selftest`). This solves a second problem whose solution is known, the field of a source outside the room, and reports the relative error at 50 interior points and on the walls. This is the accuracy measure. Note that the GMRES tolerance (`--eps-gmres`, default 1e-6) bounds it from below.
54+- Energy balance. The power radiated by the source, k/(4π) + Im v(x₀), must equal the power absorbed by the walls, −k Re(λ)∫|u|² dS. Their relative mismatch is a check that needs no exact solution.
55+
56+## Results
57+
58+Both backends were run with `run.py`, at the default settings, on the same workstation. That machine has an Intel Xeon w5-3435X (16 cores, 32 threads, all used by fmm3dbie) and an NVIDIA RTX 6000 Ada (48 GB). fmm3dbie was built at commit 94cb732 and run with eps = 1e-4. The GMRES tolerance was 1e-6 for both. The table is `results/comparison.md` as written by `compare.py`.
59+
60+| case | points | fmm3dbie solve | GPU setup + solve | speedup | GMRES iterations | self-test error, fmm3dbie / GPU | wall field difference |
61+|---|---|---|---|---|---|---|---|
62+| small | 4,704 | 4.7 s | 0.16 + 0.22 s | 12x | 24 | 1.5e-5 / 2.2e-6 | 3.8e-5 |
63+| medium | 17,496 | 18.4 s | 0.48 + 0.60 s | 17x | 24 | 1.5e-6 / 1.1e-6 | 1.4e-6 |
64+| large | 48,600 | 214 s | 1.6 + 9.2 s | 20x | 110 | 1.7e-5 / 1.0e-5 | 1.9e-5 |
65+
66+The two codes need the same number of GMRES iterations, which is consistent with their discretizing the same operator. Their wall fields agree to about the accuracy of each. The speedup grows with the problem size, because the fixed costs of the GPU code become less important. On the large case, 85% of the GPU time is GMRES. Most of that is presumably the O(N²) far-field sum, which is the part an FMM would reduce. The CUDA kernels are compiled or loaded once per process, taking about 2 s, and this is excluded from the GPU times. fmm3dbie's first solve in a process can be slower than its second: on the medium case the self-test solve took 81 s and the room solve 18 s. We report the room solve.
67+
68+At the default GMRES tolerance both codes' self-test errors are close to the tolerance. With a GMRES tolerance of 1e-10 the medium-case GPU error is 2.7e-7 with the single precision far field and 2.9e-10 with `--far64`. In an earlier run on a laptop, fmm3dbie reached 1.9e-8 there with eps 1e-7, at about four times the cost of eps 1e-4 (400 s against 92 s).
69+
70+## Files
71+
72+- `run.py`, `compare.py`, `plot.py`: run one backend, compare the results, and plot a field.
73+- `problem.py`: the showcase problem, shared by both backends: cases, source positions, right-hand sides and checks.
74+- `fmm3dbie_gpu/surface.py`: patch discretization and the rounded cube. It is pure NumPy.
75+- `fmm3dbie_gpu/cuda_ops.py`: the CUDA kernels: near pairs, quadtrees, near weights, far-field sum and near-field product.
76+- `fmm3dbie_gpu/solver.py`: the operator, evaluation, and GMRES, following fmm3dbie's `zgmres_guru`.
77+- `reference/`: the ctypes interface to fmm3dbie and the same rounded-cube surface built with fmm3dbie's routines. We use it because the Python bindings shipped with fmm3dbie need `numpy.distutils`, which NumPy 2 no longer has.
78+- `tests/test_quadrature.py`: compares the GPU near weights with brute-force integration. For example, `python3 tests/test_quadrature.py 8 40 m=10 tp=45` runs at order 8 and k = 40.
79+- `tests/geometry_test.py`: checks Green's identity at nodes next to patch edges with the interpolated, derivative-consistent and exact geometry models. This is the evidence for point 4 above.
80+
81+## Limitations
82+
83+The far field is a direct sum, so each product costs O(N²). At 48,600 points one product takes about 40 ms. Somewhere between 10⁵ and 2·10⁵ points a GPU FMM becomes necessary; FMM3D's Helmholtz scheme is the obvious model. The near weights would not change, since the FMM would only replace the dense sum.
84+
85+The exact geometry exists only for the rounded cube. General surfaces use the interpolated geometry, and with it fmm3dbie's sensitivity at nodes next to patch edges.
86+
87+The self-patch rule has been checked against brute-force integration for orders 6 and 8 and for k up to 40, but not for higher orders or more strongly curved patches.
88+
89+Only this one solver is ported. The two codes are compared on one family of problems, and the relative speed will differ for other geometries, orders and tolerances.
compare.pyadded+59−0View file
@@ -0,0 +1,59 @@
1+"""Put the results of run.py side by side.
2+
3+ python compare.py [results_dir]
4+
5+Prints a table of all runs found in results_dir (default results/) and, for
6+every case run with both backends, the speedup and the relative difference of
7+the wall fields. The table is also written to results_dir/comparison.md.
8+"""
9+
10+import glob
11+import json
12+import os
13+import sys
14+
15+import numpy as np
16+
17+
18+def fmt(x, spec):
19+ return "" if x is None else format(x, spec)
20+
21+
22+def main():
23+ rdir = sys.argv[1] if len(sys.argv) > 1 else os.path.join(os.path.dirname(os.path.abspath(__file__)), "results")
24+ runs = {}
25+ for f in sorted(glob.glob(os.path.join(rdir, "*", "summary.json"))):
26+ s = json.load(open(f))
27+ s["dir"] = os.path.dirname(f)
28+ runs[(s["case"], s["backend"])] = s
29+ if not runs:
30+ print(f"no results in {rdir}")
31+ return
32+
33+ lines = ["| case | k | points | backend | machine | setup (s) | solve (s) | setup + solve (s) | "
34+ "iterations | self-test error | energy mismatch |",
35+ "|---|---|---|---|---|---|---|---|---|---|---|"]
36+ order = {"small": 0, "medium": 1, "large": 2}
37+ for (case, backend), s in sorted(runs.items(), key=lambda kv: (order.get(kv[0][0], 9), kv[0][1])):
38+ m = s["machine"]
39+ where = m.get("gpu", "").split(",")[0] if backend == "gpu" else f"{m.get('cpu', '?')}, {m.get('cpu_count')} threads"
40+ lines.append(f"| {case} | {s['k']:g} | {s['npts']:,} | {backend} | {where} | {s['setup_s']:.2f} | "
41+ f"{s['solve_s']:.2f} | {s['setup_plus_solve_s']:.2f} | {s['niter']} | "
42+ f"{fmt(s.get('selftest_err'), '.1e')} | {s['energy_mismatch']:.1e} |")
43+
44+ lines += ["", "| case | speedup (setup + solve) | wall field difference |", "|---|---|---|"]
45+ for case in sorted({c for c, _ in runs}, key=lambda c: order.get(c, 9)):
46+ if (case, "gpu") in runs and (case, "fmm3dbie") in runs:
47+ g, c = runs[(case, "gpu")], runs[(case, "fmm3dbie")]
48+ ug = np.load(os.path.join(g["dir"], "field.npz"))["u_wall"]
49+ uc = np.load(os.path.join(c["dir"], "field.npz"))["u_wall"]
50+ diff = np.linalg.norm(ug - uc) / np.linalg.norm(uc) if ug.shape == uc.shape else float("nan")
51+ lines.append(f"| {case} | {c['setup_plus_solve_s'] / g['setup_plus_solve_s']:.0f}x | {diff:.1e} |")
52+ text = "\n".join(lines)
53+ print(text)
54+ with open(os.path.join(rdir, "comparison.md"), "w") as f:
55+ f.write(text + "\n")
56+
57+
58+if __name__ == "__main__":
59+ main()
docs/slices_large.pngadded+0−0View file
Binary file not shown.
fmm3dbie_gpu/__init__.pyadded+2−0View file
@@ -0,0 +1,2 @@
1+"""GPU solver for the interior Helmholtz impedance problem on quadrilateral patch
2+surfaces, following fmm3dbie's helm_rpcomb_imp formulation."""
fmm3dbie_gpu/cuda_ops.pyadded+612−0View file
@@ -0,0 +1,612 @@
1+"""CUDA kernels (compiled at run time with NVRTC through CuPy) for the
2+Helmholtz layer potentials on quadrilateral patch surfaces.
3+
4+Six kernels are handled together, in the order used by fmm3dbie's
5+helm_rpcomb_imp solver:
6+
7+ 0: S'_k 1: S_kappa 2: S'_kappa 3: D'_k - D'_kappa 4: S_k 5: D_k
8+
9+with kappa = i|k|, G_k = exp(ikr)/(4 pi r), S' the normal derivative at the
10+target, D the normal derivative at the source.
11+
12+Each layer potential at a target x is split by patch. A patch j is *near*
13+x when |x - c_j| <= rfac r_j (c_j, r_j the patch centroid and radius); this
14+test is the same device function in every kernel, so the split is
15+consistent. Far patches are summed with the smooth (oversampled Gauss-
16+Legendre) rule in dense_apply; near patches use precomputed weights
17+(near_weights) applied in sparse_apply. Because near patches never enter the
18+dense sum, the dense sum has no large cancelling terms and can be evaluated
19+in single precision.
20+
21+Near weights are computed as Legendre moments I_n = int K P_n |J| du dv,
22+converted to node weights with the vals-to-coefs matrix U: w_j = sum_n U[n, j] I_n.
23+
24+ * self patch (target on the patch): the square is split into four
25+ triangles fanned from the target to the edges. On each triangle we use
26+ polar-type coordinates (rho along the ray, x the position on the edge)
27+ and the substitution x = x_c + h_eff sinh(xi) chosen from the first
28+ fundamental form at the target, which removes the near-singular
29+ behaviour in the angular direction. Tensor Gauss-Legendre in (rho, xi),
30+ with the xi range split so that each piece is short relative to the
31+ singularity and to the wavelength.
32+ * other near patches: a quadtree on the reference square, refined while the
33+ target is within rfac_q times the radius of a leaf or the leaf is large
34+ compared with the wavelength, with a tensor Gauss-Legendre rule on each
35+ leaf (fmm3dbie's istrat=1 strategy).
36+
37+The geometry inside the integrals is either the polynomial interpolant of
38+the node data (gmode 0, as in fmm3dbie) or, for the rounded cube, the exact
39+superellipsoid parametrization (gmode 1). The exact geometry is much more
40+accurate at nodes next to patch edges: the double layer kernel is very
41+sensitive to small inconsistencies between interpolated positions and
42+normals.
43+"""
44+
45+import cupy as cp
46+import numpy as np
47+
48+_SRC = r"""
49+#include <cupy/complex.cuh>
50+typedef complex<double> cdouble;
51+
52+#define NORD1_MAX 13
53+#define NPOLS_MAX 169
54+#define NT 128
55+#define NSUB_MAX 32
56+#define OVER4PI 0.07957747154594767
57+
58+__device__ inline void legendre(double x, int norder, double* P) {
59+ P[0] = 1.0;
60+ if (norder >= 1) P[1] = x;
61+ for (int n = 1; n < norder; n++)
62+ P[n + 1] = ((2 * n + 1) * x * P[n] - n * P[n - 1]) / (n + 1);
63+}
64+
65+__device__ inline double ipow(double x, int n) {
66+ double r = 1.0;
67+ for (int i = 0; i < n; i++) r *= x;
68+ return r;
69+}
70+
71+// xyz, r_u, r_v at (u, v) of one patch.
72+// gmode 0: from the Legendre coefficients c[comp * npols + n] of the patch.
73+// gmode 1: superellipsoid face chart; gp = n(3), e1(3), e2(3), s0, s1, t0, t1.
74+__device__ inline void geometry(int gmode, const double* coefs, int norder, const double* gpar,
75+ double ga, int gpe, int patch, double u, double v, double* g) {
76+ if (gmode == 1) {
77+ const double* gp = gpar + 13 * (size_t)patch;
78+ double hs = 0.5 * (gp[10] - gp[9]), ht = 0.5 * (gp[12] - gp[11]);
79+ double s = gp[9] + (u + 1) * hs, t = gp[11] + (v + 1) * ht;
80+ double c[3], S = 0.0;
81+ for (int q = 0; q < 3; q++) { c[q] = gp[q] + s * gp[3 + q] + t * gp[6 + q]; S += ipow(c[q], gpe); }
82+ double N = pow(S, 1.0 / gpe);
83+ double dNu = 0.0, dNv = 0.0;
84+ for (int q = 0; q < 3; q++) {
85+ double gN = ipow(c[q] / N, gpe - 1);
86+ dNu += gN * gp[3 + q]; dNv += gN * gp[6 + q];
87+ }
88+ for (int q = 0; q < 3; q++) {
89+ g[q] = ga * c[q] / N;
90+ g[3 + q] = ga * hs * (gp[3 + q] / N - c[q] * dNu / (N * N));
91+ g[6 + q] = ga * ht * (gp[6 + q] / N - c[q] * dNv / (N * N));
92+ }
93+ return;
94+ }
95+ int n1 = norder + 1, npols = n1 * n1;
96+ const double* c = coefs + (size_t)patch * 9 * npols;
97+ double Lu[NORD1_MAX], Lv[NORD1_MAX];
98+ legendre(u, norder, Lu);
99+ legendre(v, norder, Lv);
100+ for (int comp = 0; comp < 9; comp++) {
101+ double s = 0.0;
102+ const double* cc = c + comp * npols;
103+ for (int b = 0; b < n1; b++) {
104+ double sb = 0.0;
105+ for (int a = 0; a < n1; a++) sb += cc[b * n1 + a] * Lu[a];
106+ s += sb * Lv[b];
107+ }
108+ g[comp] = s;
109+ }
110+}
111+
112+// The near test, shared by every kernel so that the near/far split is consistent.
113+__device__ inline bool is_near(const double* x, const double* cen, double r2) {
114+ double dx = x[0] - cen[0], dy = x[1] - cen[1], dz = x[2] - cen[2];
115+ return dx * dx + dy * dy + dz * dz <= r2;
116+}
117+
118+// The six kernels at one (target, source) pair, times a real weight.
119+template <typename T>
120+__device__ inline void kernels6(const T* xt, const T* nt, const T* y, const T* ns, T w,
121+ complex<T> zk, T zkabs, int selmask, complex<T>* K) {
122+ T dx = xt[0] - y[0], dy = xt[1] - y[1], dz = xt[2] - y[2];
123+ T r2 = dx * dx + dy * dy + dz * dz;
124+ T r = sqrt(r2);
125+ T rinv = (T)1 / r;
126+ T r2inv = rinv * rinv;
127+ T dnt = dx * nt[0] + dy * nt[1] + dz * nt[2];
128+ T dns = dx * ns[0] + dy * ns[1] + dz * ns[2];
129+ complex<T> izkr = complex<T>(0, 1) * zk * r;
130+ complex<T> Ek(0, 0);
131+ T Ec = 0;
132+ // selmask: 2 bits per kernel m at bits 2m, 2m+1
133+ if (selmask & 0xFC3) Ek = exp(izkr) * ((T)(OVER4PI) * rinv * w); // kernels 0, 3, 4, 5
134+ if (selmask & 0x0FC) Ec = exp(-zkabs * r) * ((T)(OVER4PI) * rinv * w); // kernels 1, 2, 3
135+ complex<T> fk = (T)1 - izkr;
136+ T fc = (T)1 + zkabs * r;
137+ K[4] = Ek;
138+ K[1] = complex<T>(Ec, 0);
139+ K[0] = -dnt * fk * Ek * r2inv;
140+ K[2] = complex<T>(-dnt * fc * Ec * r2inv, 0);
141+ K[5] = dns * fk * Ek * r2inv;
142+ if (selmask & 0x0C0) {
143+ T nn = ns[0] * nt[0] + ns[1] * nt[1] + ns[2] * nt[2];
144+ complex<T> gk = ((T)3 - (T)3 * izkr - zk * zk * r2);
145+ T gc = (T)3 + (T)3 * zkabs * r + zkabs * zkabs * r2;
146+ K[3] = (nn * fk * r2inv - dns * dnt * gk * r2inv * r2inv) * Ek
147+ - complex<T>((nn * fc * r2inv - dns * dnt * gc * r2inv * r2inv) * Ec, 0);
148+ } else {
149+ K[3] = complex<T>(0, 0);
150+ }
151+}
152+
153+// ---------------------------------------------------------------------------
154+// Near pairs: two passes (counts, then fill using row offsets), one thread per target.
155+extern "C" __global__ void near_pairs(
156+ int ntarg, const double* targ_x, int npatches, const double* cen, const double* r2,
157+ int* counts, const int* row_ptr, int* pair_p)
158+{
159+ int i = blockIdx.x * blockDim.x + threadIdx.x;
160+ if (i >= ntarg) return;
161+ double x[3] = {targ_x[3 * i], targ_x[3 * i + 1], targ_x[3 * i + 2]};
162+ int n = 0;
163+ for (int j = 0; j < npatches; j++) {
164+ if (is_near(x, cen + 3 * j, r2[j])) {
165+ if (pair_p != NULL) pair_p[row_ptr[i] + n] = j;
166+ n++;
167+ }
168+ }
169+ counts[i] = n;
170+}
171+
172+// ---------------------------------------------------------------------------
173+// Quadtree leaves for the (target, patch) pairs that are not self pairs.
174+// Two passes: counts only (leaves == NULL), then fill using offsets.
175+extern "C" __global__ void build_leaves(
176+ int npairs, const int* pair_t, const int* pair_p, const double* targ_x,
177+ const int* targ_patch, int gmode, const double* coefs, int norder, const double* gpar,
178+ double ga, int gpe, double rfac_q, int nlmax, double zkabs, double kR_max,
179+ int* counts, const int* offsets, double* leaves)
180+{
181+ int p = blockIdx.x * blockDim.x + threadIdx.x;
182+ if (p >= npairs) return;
183+ int it = pair_t[p], ip = pair_p[p];
184+ if (targ_patch[it] == ip) { counts[p] = 0; return; }
185+ double tx = targ_x[3 * it], ty = targ_x[3 * it + 1], tz = targ_x[3 * it + 2];
186+ double su[64], sv[64], sh[64];
187+ int sd[64];
188+ int sp = 1;
189+ su[0] = 0.0; sv[0] = 0.0; sh[0] = 1.0; sd[0] = 0;
190+ int n = 0;
191+ while (sp > 0) {
192+ sp--;
193+ double cu = su[sp], cv = sv[sp], h = sh[sp];
194+ int d = sd[sp];
195+ bool refine = false;
196+ if (d < nlmax) {
197+ // centroid and radius from a 3x3 sample grid (corners alone underestimate
198+ // the extent of strongly curved patches)
199+ double xc[9][3], cen[3] = {0, 0, 0}, g[9];
200+ for (int k = 0; k < 9; k++) {
201+ geometry(gmode, coefs, norder, gpar, ga, gpe, ip, cu + h * (k % 3 - 1), cv + h * (k / 3 - 1), g);
202+ for (int q = 0; q < 3; q++) { xc[k][q] = g[q]; cen[q] += g[q] / 9.0; }
203+ }
204+ double R2 = 0.0;
205+ for (int k = 0; k < 9; k++) {
206+ double s = 0.0;
207+ for (int q = 0; q < 3; q++) s += (xc[k][q] - cen[q]) * (xc[k][q] - cen[q]);
208+ R2 = fmax(R2, s);
209+ }
210+ double d2 = (tx - cen[0]) * (tx - cen[0]) + (ty - cen[1]) * (ty - cen[1]) + (tz - cen[2]) * (tz - cen[2]);
211+ // refine while the target is close relative to the leaf, or the leaf is
212+ // large compared with the wavelength (the leaf rule must resolve exp(ikr))
213+ refine = (d2 < rfac_q * rfac_q * R2) || (zkabs * zkabs * R2 > kR_max * kR_max);
214+ }
215+ if (refine) {
216+ double hh = 0.5 * h;
217+ const double cs[4][2] = {{-1, -1}, {1, -1}, {1, 1}, {-1, 1}};
218+ for (int k = 0; k < 4; k++) {
219+ su[sp] = cu + hh * cs[k][0]; sv[sp] = cv + hh * cs[k][1]; sh[sp] = hh; sd[sp] = d + 1; sp++;
220+ }
221+ } else {
222+ if (leaves != NULL) {
223+ size_t o = (size_t)(offsets[p] + n) * 3;
224+ leaves[o] = cu; leaves[o + 1] = cv; leaves[o + 2] = h;
225+ }
226+ n++;
227+ }
228+ }
229+ counts[p] = n;
230+}
231+
232+// ---------------------------------------------------------------------------
233+// Near-field weights: one block per (target, patch) pair.
234+struct SelfItem { int edge; double xi0, xi1; };
235+
236+extern "C" __global__ void near_weights(
237+ int npairs, const int* pair_t, const int* pair_p,
238+ const double* targ_x, const double* targ_n, const int* targ_patch, const double* targ_uv,
239+ int gmode, const double* coefs, int norder, const double* gpar, double ga, int gpe,
240+ const int* leaf_ptr, const double* leaves, int nq, const double* gl_x, const double* gl_w,
241+ int nqr, const double* glr_x, const double* glr_w, int nqx, const double* glx_x, const double* glx_w,
242+ double sub_len, double osc_len, const double* umat, cdouble zk, double zkabs, cdouble* W)
243+{
244+ __shared__ double sKr[6][NT], sKi[6][NT];
245+ __shared__ double sLu[NORD1_MAX][NT], sLv[NORD1_MAX][NT];
246+ __shared__ SelfItem sitems[4 * NSUB_MAX];
247+ __shared__ double sEdge[4][3]; // per edge: h, x_c, h_eff
248+ __shared__ int snitems;
249+
250+ int p = blockIdx.x;
251+ if (p >= npairs) return;
252+ int t = threadIdx.x;
253+ int it = pair_t[p], ip = pair_p[p];
254+ int n1 = norder + 1, npols = n1 * n1;
255+ double xt[3] = {targ_x[3 * it], targ_x[3 * it + 1], targ_x[3 * it + 2]};
256+ double nt[3] = {targ_n[3 * it], targ_n[3 * it + 1], targ_n[3 * it + 2]};
257+ bool self = (targ_patch[it] == ip);
258+ double u0 = 0, v0 = 0;
259+
260+ // per-thread moment accumulators for basis indices n = t, t + NT
261+ const int NACC = (NPOLS_MAX + NT - 1) / NT;
262+ double Ir[NACC][6], Ii[NACC][6];
263+ for (int a = 0; a < NACC; a++) for (int m = 0; m < 6; m++) { Ir[a][m] = 0; Ii[a][m] = 0; }
264+
265+ int nitems, nodes_per_item;
266+ if (self) {
267+ u0 = targ_uv[2 * it]; v0 = targ_uv[2 * it + 1];
268+ if (t == 0) {
269+ double g[9];
270+ geometry(gmode, coefs, norder, gpar, ga, gpe, ip, u0, v0, g);
271+ double E = g[3] * g[3] + g[4] * g[4] + g[5] * g[5];
272+ double F = g[3] * g[6] + g[4] * g[7] + g[5] * g[8];
273+ double G = g[6] * g[6] + g[7] * g[7] + g[8] * g[8];
274+ double det = sqrt(fmax(E * G - F * F, 1e-300));
275+ int cnt = 0;
276+ for (int e = 0; e < 4; e++) {
277+ double h, x1, x2, xc, heff;
278+ if (e == 0) { h = 1 - u0; } else if (e == 1) { h = -1 - u0; }
279+ else if (e == 2) { h = 1 - v0; } else { h = -1 - v0; }
280+ if (e < 2) { x1 = -1 - v0; x2 = 1 - v0; xc = -h * F / G; heff = fabs(h) * det / G; }
281+ else { x1 = -1 - u0; x2 = 1 - u0; xc = -h * F / E; heff = fabs(h) * det / E; }
282+ sEdge[e][0] = h; sEdge[e][1] = xc; sEdge[e][2] = heff;
283+ if (fabs(h) < 1e-12) continue;
284+ double xi0 = asinh((x1 - xc) / heff), xi1 = asinh((x2 - xc) / heff);
285+ // sub-intervals of length <= sub_len in xi (singular structure) whose
286+ // extent along the edge is <= dxmax (so that exp(ikr) is resolved)
287+ double arc = sqrt(e < 2 ? G : E);
288+ double dxmax = zkabs * arc > 0 ? osc_len / (zkabs * arc) : 1e300;
289+ double xi = xi0;
290+ while (xi < xi1 && cnt < 4 * NSUB_MAX) {
291+ double xh = xc + heff * sinh(xi);
292+ double xn = fmin(xi + sub_len, asinh((xh + dxmax - xc) / heff));
293+ if (xn > xi1 || cnt == 4 * NSUB_MAX - 1) xn = xi1;
294+ sitems[cnt].edge = e; sitems[cnt].xi0 = xi; sitems[cnt].xi1 = xn; cnt++;
295+ xi = xn;
296+ }
297+ }
298+ snitems = cnt;
299+ }
300+ __syncthreads();
301+ nitems = snitems;
302+ nodes_per_item = nqr * nqx;
303+ } else {
304+ nitems = leaf_ptr[p + 1] - leaf_ptr[p];
305+ nodes_per_item = nq * nq;
306+ }
307+
308+ int ntot = nitems * nodes_per_item;
309+ for (int base = 0; base < ntot; base += NT) {
310+ int g = base + t;
311+ double Lu[NORD1_MAX], Lv[NORD1_MAX];
312+ cdouble K[6];
313+ bool active = g < ntot;
314+ if (active) {
315+ int item = g / nodes_per_item, loc = g - item * nodes_per_item;
316+ double u, v, wq;
317+ if (self) {
318+ int e = sitems[item].edge;
319+ int ir = loc / nqx, ix = loc - ir * nqx;
320+ double h = sEdge[e][0], xc = sEdge[e][1], heff = sEdge[e][2];
321+ double rho = glr_x[ir], wr = glr_w[ir];
322+ double xi0 = sitems[item].xi0, xi1 = sitems[item].xi1;
323+ double xi = xi0 + 0.5 * (xi1 - xi0) * (glx_x[ix] + 1.0);
324+ double wxi = 0.5 * (xi1 - xi0) * glx_w[ix];
325+ double x = xc + heff * sinh(xi);
326+ double dx = heff * cosh(xi) * wxi;
327+ if (e < 2) { u = u0 + rho * h; v = v0 + rho * x; }
328+ else { u = u0 + rho * x; v = v0 + rho * h; }
329+ wq = fabs(h) * rho * wr * dx;
330+ } else {
331+ int li = leaf_ptr[p] + item;
332+ double cu = leaves[3 * li], cv = leaves[3 * li + 1], hf = leaves[3 * li + 2];
333+ int a = loc / nq, b = loc - a * nq;
334+ u = cu + hf * gl_x[a]; v = cv + hf * gl_x[b];
335+ wq = hf * hf * gl_w[a] * gl_w[b];
336+ }
337+ double gg[9];
338+ geometry(gmode, coefs, norder, gpar, ga, gpe, ip, u, v, gg);
339+ double ns[3] = {gg[4] * gg[8] - gg[5] * gg[7], gg[5] * gg[6] - gg[3] * gg[8], gg[3] * gg[7] - gg[4] * gg[6]};
340+ double jac = sqrt(ns[0] * ns[0] + ns[1] * ns[1] + ns[2] * ns[2]);
341+ ns[0] /= jac; ns[1] /= jac; ns[2] /= jac;
342+ kernels6<double>(xt, nt, gg, ns, wq * jac, zk, zkabs, 0xFFF, K);
343+ legendre(u, norder, Lu);
344+ legendre(v, norder, Lv);
345+ }
346+ for (int m = 0; m < 6; m++) { sKr[m][t] = active ? K[m].real() : 0.0; sKi[m][t] = active ? K[m].imag() : 0.0; }
347+ for (int a = 0; a < n1; a++) { sLu[a][t] = active ? Lu[a] : 0.0; sLv[a][t] = active ? Lv[a] : 0.0; }
348+ __syncthreads();
349+ for (int acc = 0; acc < NACC; acc++) {
350+ int n = t + acc * NT;
351+ if (n < npols) {
352+ int a = n % n1, b = n / n1;
353+ for (int i = 0; i < NT; i++) {
354+ double P = sLu[a][i] * sLv[b][i];
355+ for (int m = 0; m < 6; m++) { Ir[acc][m] += P * sKr[m][i]; Ii[acc][m] += P * sKi[m][i]; }
356+ }
357+ }
358+ }
359+ __syncthreads();
360+ }
361+
362+ // moments -> node weights, w_j = sum_n U[n, j] I_n. Reuse sLu/sLv as I storage.
363+ double* sIr = &sLu[0][0]; // NORD1_MAX * NT = 1664 >= 6 * NPOLS_MAX = 1014
364+ double* sIi = &sLv[0][0];
365+ for (int acc = 0; acc < NACC; acc++) {
366+ int n = t + acc * NT;
367+ if (n < npols) for (int m = 0; m < 6; m++) { sIr[m * NPOLS_MAX + n] = Ir[acc][m]; sIi[m * NPOLS_MAX + n] = Ii[acc][m]; }
368+ }
369+ __syncthreads();
370+ for (int j = t; j < npols; j += NT) {
371+ for (int m = 0; m < 6; m++) {
372+ double wr = 0, wi = 0;
373+ for (int n = 0; n < npols; n++) {
374+ double u = umat[n * npols + j];
375+ wr += u * sIr[m * NPOLS_MAX + n];
376+ wi += u * sIi[m * NPOLS_MAX + n];
377+ }
378+ W[((size_t)p * 6 + m) * npols + j] = cdouble(wr, wi);
379+ }
380+ }
381+}
382+
383+// ---------------------------------------------------------------------------
384+// Far-field pass: out[m][i] = sum over patches j not near x_i of
385+// sum_l K_m(x_i, y_l) w_l q_l over the oversampled nodes l of patch j, for the
386+// kernels selected in selmask (4 bits per kernel: 0 skip, 1 density qa, 2 density qb).
387+// Kernel evaluation in precision T; per-patch partial sums are added into double accumulators.
388+template <typename T>
389+__global__ void dense_apply(
390+ int ntarg, const double* targ_x, const T* targ_n, int npatches, int npols_o,
391+ const double* cen, const double* r2, const T* src_x, const T* src_n, const T* src_w,
392+ const complex<T>* qa, const complex<T>* qb, complex<T> zk, T zkabs, int selmask, cdouble* out)
393+{
394+ __shared__ T tile[NT][11];
395+ __shared__ double scen[4];
396+ int i = blockIdx.x * blockDim.x + threadIdx.x;
397+ double xd[3] = {0, 0, 0};
398+ T xt[3] = {0, 0, 0}, nt[3] = {0, 0, 0};
399+ if (i < ntarg) {
400+ for (int q = 0; q < 3; q++) { xd[q] = targ_x[3 * i + q]; xt[q] = (T)xd[q]; nt[q] = targ_n[3 * i + q]; }
401+ }
402+ int sel[6];
403+ int kmask = 0;
404+ for (int m = 0; m < 6; m++) { sel[m] = (selmask >> (4 * m)) & 3; if (sel[m]) kmask |= 3 << (2 * m); }
405+ cdouble acc[6];
406+ for (int m = 0; m < 6; m++) acc[m] = cdouble(0, 0);
407+ for (int j = 0; j < npatches; j++) {
408+ if (threadIdx.x < 3) scen[threadIdx.x] = cen[3 * j + threadIdx.x];
409+ if (threadIdx.x == 3) scen[3] = r2[j];
410+ __syncthreads();
411+ bool far = (i < ntarg) && !is_near(xd, scen, scen[3]);
412+ for (int base = 0; base < npols_o; base += NT) {
413+ int l = base + threadIdx.x;
414+ if (l < npols_o) {
415+ size_t s = (size_t)j * npols_o + l;
416+ for (int q = 0; q < 3; q++) { tile[threadIdx.x][q] = src_x[3 * s + q]; tile[threadIdx.x][3 + q] = src_n[3 * s + q]; }
417+ tile[threadIdx.x][6] = src_w[s];
418+ tile[threadIdx.x][7] = qa[s].real(); tile[threadIdx.x][8] = qa[s].imag();
419+ tile[threadIdx.x][9] = qb[s].real(); tile[threadIdx.x][10] = qb[s].imag();
420+ }
421+ __syncthreads();
422+ if (far) {
423+ int nl = min(NT, npols_o - base);
424+ complex<T> part[6];
425+ for (int m = 0; m < 6; m++) part[m] = complex<T>(0, 0);
426+ for (int s = 0; s < nl; s++) {
427+ complex<T> K[6];
428+ kernels6<T>(xt, nt, &tile[s][0], &tile[s][3], tile[s][6], zk, zkabs, kmask, K);
429+ complex<T> a(tile[s][7], tile[s][8]), b(tile[s][9], tile[s][10]);
430+ for (int m = 0; m < 6; m++) {
431+ if (sel[m] == 1) part[m] += K[m] * a;
432+ else if (sel[m] == 2) part[m] += K[m] * b;
433+ }
434+ }
435+ for (int m = 0; m < 6; m++) acc[m] += cdouble((double)part[m].real(), (double)part[m].imag());
436+ }
437+ __syncthreads();
438+ }
439+ }
440+ if (i < ntarg) for (int m = 0; m < 6; m++) out[(size_t)m * ntarg + i] = acc[m];
441+}
442+
443+// ---------------------------------------------------------------------------
444+// Near contributions: out[m][i] += sum_{near pairs p of i} sum_l W[p][m][l] q[d_m][patch nodes l]
445+extern "C" __global__ void sparse_apply(
446+ int ntarg, const int* row_ptr, const int* pair_p, int npols, int nreq, const int* req_m,
447+ const int* req_d, const cdouble* W, const cdouble* q, int npts, cdouble* out)
448+{
449+ int i = blockIdx.x * blockDim.x + threadIdx.x;
450+ if (i >= ntarg) return;
451+ cdouble acc[6];
452+ for (int r = 0; r < nreq; r++) acc[r] = cdouble(0, 0);
453+ for (int p = row_ptr[i]; p < row_ptr[i + 1]; p++) {
454+ int j = pair_p[p];
455+ for (int r = 0; r < nreq; r++) {
456+ int m = req_m[r];
457+ const cdouble* w = W + ((size_t)p * 6 + m) * npols;
458+ const cdouble* qq = q + (size_t)req_d[r] * npts + (size_t)j * npols;
459+ cdouble s(0, 0);
460+ for (int l = 0; l < npols; l++) s += w[l] * qq[l];
461+ acc[r] += s;
462+ }
463+ }
464+ for (int r = 0; r < nreq; r++) out[(size_t)req_m[r] * ntarg + i] += acc[r];
465+}
466+"""
467+
468+_module = None
469+
470+
471+def module():
472+ global _module
473+ if _module is None:
474+ _module = cp.RawModule(code=_SRC, options=("-std=c++17",),
475+ name_expressions=["dense_apply<double>", "dense_apply<float>"])
476+ _module.compile() # compile (or load from CuPy's cache) now rather than on first use
477+ return _module
478+
479+
480+def _i32(a):
481+ return cp.ascontiguousarray(cp.asarray(a, dtype=cp.int32))
482+
483+
484+def _f64(a):
485+ return cp.ascontiguousarray(cp.asarray(a, dtype=cp.float64))
486+
487+
488+class Geometry:
489+ """Device-side description of the patch geometry used inside the integrals."""
490+
491+ def __init__(self, surf):
492+ npols = surf.npols
493+ self.norder = surf.norder
494+ self.coefs = cp.asarray(np.ascontiguousarray(
495+ surf.srccoefs.reshape(9, surf.npatches, npols).transpose(1, 0, 2)))
496+ if getattr(surf, "exact", None) is not None:
497+ self.gmode = 1
498+ self.gpar = _f64(surf.exact["params"])
499+ self.ga = float(surf.exact["a"])
500+ self.gpe = int(surf.exact["p"])
501+ else:
502+ self.gmode, self.gpar, self.ga, self.gpe = 0, cp.zeros(13), 0.0, 0
503+
504+ def args(self):
505+ return (np.int32(self.gmode), self.coefs, np.int32(self.norder), self.gpar,
506+ np.float64(self.ga), np.int32(self.gpe))
507+
508+
509+def find_near_pairs(targ_x, cen, r2):
510+ """Pairs (target, patch) with |x_t - c_j|^2 <= r2_j, grouped by target.
511+
512+ Returns (pair_t, pair_p, row_ptr) as int32 device arrays.
513+ """
514+ ntarg = int(targ_x.shape[0])
515+ npatches = int(cen.shape[0])
516+ ker = module().get_function("near_pairs")
517+ counts = cp.zeros(ntarg, dtype=cp.int32)
518+ nb = (ntarg + 127) // 128
519+ ker((nb,), (128,), (np.int32(ntarg), targ_x, np.int32(npatches), cen, r2, counts,
520+ cp.zeros(1, cp.int32), np.uintp(0)))
521+ row_ptr = cp.zeros(ntarg + 1, dtype=cp.int32)
522+ row_ptr[1:] = cp.cumsum(counts).astype(cp.int32)
523+ npairs = int(row_ptr[-1])
524+ pair_p = cp.zeros(max(npairs, 1), dtype=cp.int32)
525+ ker((nb,), (128,), (np.int32(ntarg), targ_x, np.int32(npatches), cen, r2, counts, row_ptr, pair_p))
526+ pair_p = pair_p[:npairs]
527+ pair_t = cp.repeat(cp.arange(ntarg, dtype=cp.int32), counts.tolist()) if npairs else cp.zeros(0, cp.int32)
528+ return pair_t.astype(cp.int32), pair_p, row_ptr
529+
530+
531+def build_leaves(pair_t, pair_p, targ_x, targ_patch, geo, zk=0.0, rfac_q=1.5, nlmax=14, kR_max=3.0):
532+ npairs = int(pair_t.size)
533+ if npairs == 0:
534+ return cp.zeros(1, cp.int32), cp.zeros((1, 3), cp.float64)
535+ ker = module().get_function("build_leaves")
536+ counts = cp.zeros(npairs, dtype=cp.int32)
537+ nb = (npairs + 127) // 128
538+ args = (np.int32(npairs), pair_t, pair_p, targ_x, targ_patch) + geo.args() + (
539+ np.float64(rfac_q), np.int32(nlmax), np.float64(abs(zk)), np.float64(kR_max))
540+ ker((nb,), (128,), args + (counts, cp.zeros(1, cp.int32), np.uintp(0)))
541+ leaf_ptr = cp.zeros(npairs + 1, dtype=cp.int32)
542+ leaf_ptr[1:] = cp.cumsum(counts).astype(cp.int32)
543+ nleaves = int(leaf_ptr[-1])
544+ leaves = cp.zeros((max(nleaves, 1), 3), dtype=cp.float64)
545+ ker((nb,), (128,), args + (counts, leaf_ptr, leaves))
546+ return leaf_ptr, leaves
547+
548+
549+def near_weights(pair_t, pair_p, targ_x, targ_n, targ_patch, targ_uv, geo, leaf_ptr, leaves,
550+ umat, zk, nq=8, nqr=10, nqx=8, sub_len=2.0, osc_len=2.0):
551+ """Near weights W (npairs, 6, npols), complex128 on the device."""
552+ npairs = int(pair_t.size)
553+ npols = (geo.norder + 1) ** 2
554+ W = cp.zeros((max(npairs, 1), 6, npols), dtype=cp.complex128)
555+ if npairs == 0:
556+ return W
557+ gx, gw = np.polynomial.legendre.leggauss(nq)
558+ rx, rw = np.polynomial.legendre.leggauss(nqr)
559+ rx, rw = 0.5 * (rx + 1), 0.5 * rw # on [0, 1]
560+ xx, xw = np.polynomial.legendre.leggauss(nqx)
561+ ker = module().get_function("near_weights")
562+ ker((npairs,), (128,), (
563+ np.int32(npairs), pair_t, pair_p, targ_x, targ_n, targ_patch, targ_uv) + geo.args() + (
564+ leaf_ptr, leaves, np.int32(nq), _f64(gx), _f64(gw),
565+ np.int32(nqr), _f64(rx), _f64(rw), np.int32(nqx), _f64(xx), _f64(xw), np.float64(sub_len),
566+ np.float64(osc_len), umat, np.complex128(zk), np.float64(abs(zk)), W))
567+ return W
568+
569+
570+class FarSources:
571+ """Oversampled nodes of all patches, patch-major, in the precision of the far field."""
572+
573+ def __init__(self, x, n, w, npols_o, cen, r2, dtype):
574+ self.real = np.dtype(dtype)
575+ self.cplx = np.complex128 if self.real == np.float64 else np.complex64
576+ self.x = cp.ascontiguousarray(cp.asarray(x).astype(self.real))
577+ self.n = cp.ascontiguousarray(cp.asarray(n).astype(self.real))
578+ self.w = cp.ascontiguousarray(cp.asarray(w).astype(self.real))
579+ self.npols_o = npols_o
580+ self.cen, self.r2 = cen, r2
581+ self.npatches = int(cen.shape[0])
582+
583+
584+def dense_apply(targ_x, targ_n, far, qa, qb, zk, sel):
585+ """out (6, ntarg) complex128: far-patch sums for the kernels in sel (6 entries in {0,1,2})."""
586+ ntarg = int(targ_x.shape[0])
587+ name = "dense_apply<double>" if far.real == np.float64 else "dense_apply<float>"
588+ ker = module().get_function(name)
589+ ccast = lambda a: cp.ascontiguousarray(cp.asarray(a).astype(far.cplx, copy=False))
590+ out = cp.zeros((6, ntarg), dtype=cp.complex128)
591+ selmask = sum(int(s) << (4 * m) for m, s in enumerate(sel))
592+ nb = (ntarg + 127) // 128
593+ ker((nb,), (128,), (
594+ np.int32(ntarg), targ_x, cp.ascontiguousarray(targ_n.astype(far.real)), np.int32(far.npatches),
595+ np.int32(far.npols_o), far.cen, far.r2, far.x, far.n, far.w, ccast(qa), ccast(qb),
596+ far.cplx(zk), far.real.type(abs(zk)), np.int32(selmask), out))
597+ return out
598+
599+
600+def sparse_apply(row_ptr, pair_p, npols, W, q, reqs, out):
601+ """out (6, ntarg) += near contributions; reqs = list of (kernel m, density index d); q (ndens, npts)."""
602+ ntarg = int(row_ptr.size) - 1
603+ if len(reqs) == 0 or int(pair_p.size) == 0:
604+ return out
605+ ker = module().get_function("sparse_apply")
606+ req_m = _i32([m for m, _ in reqs])
607+ req_d = _i32([d for _, d in reqs])
608+ q = cp.ascontiguousarray(q, dtype=cp.complex128)
609+ nb = (ntarg + 127) // 128
610+ ker((nb,), (128,), (np.int32(ntarg), row_ptr, pair_p, np.int32(npols), np.int32(len(reqs)),
611+ req_m, req_d, W, q, np.int32(q.shape[1]), out))
612+ return out
fmm3dbie_gpu/solver.pyadded+211−0View file
@@ -0,0 +1,211 @@
1+"""GPU solver for the Helmholtz impedance problem with the right-preconditioned
2+combined field representation of fmm3dbie (helm_rpcomb_imp):
3+
4+ u = S_k[sigma] + i alpha D_k S_kappa[sigma], kappa = i|k|,
5+
6+with the boundary condition du/dn + i k lambda u = f. The discretized equation
7+(fmm3dbie, helm_rpcomb_imp.f90) is
8+
9+ z sigma + S'_k sigma + i alpha S'_kappa phi2 + i alpha (D'_k - D'_kappa) phi1
10+ + i k lambda [S_k sigma + i alpha D_k phi1 + i alpha w phi1] = f,
11+
12+ phi1 = S_kappa sigma, phi2 = S'_kappa sigma,
13+ z = (-1)^ifinout / 2 - i alpha / 4, w = 1/2 (exterior), -1/2 (interior).
14+
15+The hypersingular operator D'_k is never applied directly: the Calderon
16+identity D'_kappa S_kappa = -1/4 + S'_kappa^2 turns it into the weakly singular
17+difference D'_k - D'_kappa plus a product of two S'_kappa applications.
18+
19+Every layer potential is the sum of a far part (smooth rule on the oversampled
20+nodes of the patches that are not near the target, cuda_ops.dense_apply) and a
21+near part (precomputed weights for the near patches, cuda_ops.sparse_apply).
22+"""
23+
24+import time
25+
26+import cupy as cp
27+import numpy as np
28+
29+from . import cuda_ops as ops
30+from . import surface
31+
32+
33+class Operator:
34+ def __init__(self, surf, zk, alpha, zlams, ifinout=0, nover=None, rfac=None, rfac_q=2.0,
35+ nq=8, nqr=None, nqx=8, sub_len=1.0, osc_len=2.0, kR_max=3.0, nlmax=14,
36+ far_dtype=np.float32, verbose=True):
37+ t0 = time.time()
38+ self.surf = surf
39+ self.zk = complex(zk)
40+ self.alpha = complex(alpha)
41+ self.ifinout = ifinout
42+ norder, npols = surf.norder, surf.npols
43+ self.norder, self.npols, self.npts = norder, npols, surf.npts
44+ if rfac is None:
45+ rfac = 2.0 if norder <= 6 else 1.25 # fmm3dbie's get_rfacs
46+ self.rfac = rfac
47+ self.nover = norder + 2 if nover is None else nover
48+ if nqr is None: # the radial rule must resolve the basis and exp(ikr) along a ray
49+ nqr = int(min(max(12, norder + 2 * abs(zk) * surf.rads.max() + 6), 32))
50+ self.quad = dict(nq=nq, nqr=nqr, nqx=nqx, sub_len=sub_len, osc_len=osc_len)
51+ self.leafpar = dict(rfac_q=rfac_q, nlmax=nlmax, kR_max=kR_max)
52+
53+ self.zlams = cp.asarray(np.broadcast_to(np.asarray(zlams, dtype=complex), (surf.npts,)))
54+ self.geo = ops.Geometry(surf)
55+ self.umat = cp.asarray(np.ascontiguousarray(surf.umat))
56+ self.x = cp.asarray(np.ascontiguousarray(surf.srcvals[0:3].T))
57+ self.n = cp.asarray(np.ascontiguousarray(surf.srcvals[9:12].T))
58+ self.patch_id = cp.asarray(surf.patch_id.astype(np.int32))
59+ self.uv = cp.asarray(np.ascontiguousarray(surf.uv_nodes.T))
60+ self.wts = cp.asarray(surf.wts)
61+ self.cen = cp.asarray(np.ascontiguousarray(surf.centroids.T))
62+ self.r2 = cp.asarray((rfac * surf.rads) ** 2)
63+
64+ ox, on, ow = surf.oversampled(self.nover)
65+ self.npols_o = (self.nover + 1) ** 2
66+ self.far = ops.FarSources(np.ascontiguousarray(ox.T), np.ascontiguousarray(on.T), ow,
67+ self.npols_o, self.cen, self.r2, far_dtype)
68+ self.over_mat = cp.asarray(surface.oversample_matrix(norder, self.nover))
69+
70+ # near weights for on-surface targets (the nodes themselves)
71+ self.pair_t, self.pair_p, self.row_ptr = ops.find_near_pairs(self.x, self.cen, self.r2)
72+ self.W, self.nleaves = self._weights(self.pair_t, self.pair_p, self.x, self.n,
73+ self.patch_id, self.uv)
74+ cp.cuda.Device().synchronize()
75+ self.setup_time = time.time() - t0
76+ if verbose:
77+ print(f" operator setup: {self.npts} nodes, {ox.shape[1]} oversampled nodes, "
78+ f"{int(self.pair_t.size)} near pairs, {self.nleaves} leaves, "
79+ f"{'exact' if self.geo.gmode == 1 else 'interpolated'} geometry, "
80+ f"{np.dtype(far_dtype).name} far field, {self.setup_time:.2f} s")
81+
82+ def _weights(self, pair_t, pair_p, tx, tn, tpatch, tuv):
83+ leaf_ptr, leaves = ops.build_leaves(pair_t, pair_p, tx, tpatch, self.geo, zk=self.zk,
84+ **self.leafpar)
85+ W = ops.near_weights(pair_t, pair_p, tx, tn, tpatch, tuv, self.geo, leaf_ptr, leaves,
86+ self.umat, self.zk, **self.quad)
87+ return W, int(leaf_ptr[-1])
88+
89+ def oversample(self, q):
90+ qp = q.reshape(self.surf.npatches, self.npols)
91+ return (qp @ self.over_mat.T).ravel()
92+
93+ def _dense(self, tx, tn, qa, qb, sel):
94+ return ops.dense_apply(tx, tn, self.far, self.oversample(qa), self.oversample(qb),
95+ self.zk, sel)
96+
97+ def zid(self):
98+ return (-1) ** self.ifinout / 2 - 1j * self.alpha / 4
99+
100+ def apply_K(self, sigma):
101+ """The operator without the identity term z."""
102+ k, a, lam = self.zk, self.alpha, self.zlams
103+ zero = cp.zeros_like(sigma)
104+ F = self._dense(self.x, self.n, sigma, zero, [1, 1, 1, 0, 1, 0])
105+ ops.sparse_apply(self.row_ptr, self.pair_p, self.npols, self.W, sigma[None, :],
106+ [(0, 0), (1, 0), (2, 0), (4, 0)], F)
107+ phi1, phi2 = F[1], F[2]
108+ res = F[0] + 1j * k * lam * F[4]
109+ G = self._dense(self.x, self.n, phi1, phi2, [0, 0, 2, 1, 0, 1])
110+ ops.sparse_apply(self.row_ptr, self.pair_p, self.npols, self.W, cp.stack([phi1, phi2]),
111+ [(2, 1), (3, 0), (5, 0)], G)
112+ res += 1j * a * G[2] + 1j * a * G[3] - k * a * lam * G[5]
113+ res += (-1) ** self.ifinout * k * a * lam * phi1 / 2
114+ return res
115+
116+ def apply_single(self, m, q):
117+ """One of the six operators applied to q on the surface (principal value)."""
118+ q = cp.asarray(q, dtype=cp.complex128)
119+ sel = [0] * 6
120+ sel[m] = 1
121+ F = self._dense(self.x, self.n, q, cp.zeros_like(q), sel)
122+ ops.sparse_apply(self.row_ptr, self.pair_p, self.npols, self.W, q[None, :], [(m, 0)], F)
123+ return F[m]
124+
125+ def s_kappa(self, sigma):
126+ """phi1 = S_kappa[sigma] on the surface."""
127+ return self.apply_single(1, sigma)
128+
129+ def solve(self, rhs, eps_gmres=1e-6, numit=400, verbose=True):
130+ t0 = time.time()
131+ rhs = cp.asarray(rhs, dtype=cp.complex128)
132+ sigma, niter, rres = gmres(self.apply_K, rhs, self.wts, self.zid(), eps_gmres, numit)
133+ phi1 = self.s_kappa(sigma)
134+ cp.cuda.Device().synchronize()
135+ if verbose:
136+ print(f" GMRES: {niter} iterations, relative residual {rres:.2e}, {time.time() - t0:.2f} s")
137+ return sigma, phi1, niter, rres
138+
139+ def eval(self, targ_x, sigma, phi1, on_surface=False):
140+ """u = S_k sigma + i alpha D_k phi1 at targets (ntarg, 3).
141+
142+ With on_surface=True the targets are the surface nodes and the value is
143+ the limit from the side selected by ifinout.
144+ """
145+ a = self.alpha
146+ if on_surface:
147+ tx, tn = self.x, self.n
148+ row_ptr, pair_p, W = self.row_ptr, self.pair_p, self.W
149+ else:
150+ tx = cp.ascontiguousarray(cp.asarray(targ_x, dtype=cp.float64))
151+ tn = cp.zeros_like(tx)
152+ pair_t, pair_p, row_ptr = ops.find_near_pairs(tx, self.cen, self.r2)
153+ tpatch = cp.full(tx.shape[0], -1, dtype=cp.int32)
154+ tuv = cp.zeros((tx.shape[0], 2))
155+ W, _ = self._weights(pair_t, pair_p, tx, tn, tpatch, tuv)
156+ F = self._dense(tx, tn, sigma, phi1, [0, 0, 0, 0, 1, 2])
157+ ops.sparse_apply(row_ptr, pair_p, self.npols, W, cp.stack([sigma, phi1]),
158+ [(4, 0), (5, 1)], F)
159+ u = F[4] + 1j * a * F[5]
160+ if on_surface:
161+ u += 1j * a * (-1) ** (self.ifinout + 1) * phi1 / 2
162+ return u
163+
164+
165+def gmres(apply_K, rhs, wts, zid, eps, numit):
166+ """Unrestarted GMRES for (zid I + K) x = rhs with the weighted inner product
167+ <u, v> = sum conj(u) v wts, following fmm3dbie's zgmres_guru."""
168+ n = rhs.size
169+ wnorm = lambda v: float(cp.sqrt(cp.sum(cp.abs(v) ** 2 * wts)).real)
170+ rb = wnorm(rhs)
171+ V = cp.zeros((numit + 1, n), dtype=cp.complex128)
172+ H = np.zeros((numit + 1, numit), dtype=complex)
173+ cs, sn = np.zeros(numit, dtype=complex), np.zeros(numit, dtype=complex)
174+ svec = np.zeros(numit + 1, dtype=complex)
175+ svec[0] = rb
176+ V[0] = rhs / rb
177+ it_done = 0
178+ for it in range(numit):
179+ w = apply_K(V[it])
180+ # modified Gram-Schmidt in blocks (classical within one pass, repeated once)
181+ for _ in range(2):
182+ h = (cp.conj(V[:it + 1]) * (w * wts)[None, :]).sum(axis=1)
183+ w = w - h @ V[:it + 1]
184+ H[:it + 1, it] += h.get()
185+ H[it, it] += zid
186+ hn = wnorm(w)
187+ H[it + 1, it] = hn
188+ if hn > 0:
189+ V[it + 1] = w / hn
190+ for k in range(it):
191+ t1 = cs[k] * H[k, it] + sn[k] * H[k + 1, it]
192+ t2 = -np.conj(sn[k]) * H[k, it] + cs[k] * H[k + 1, it]
193+ H[k, it], H[k + 1, it] = t1, t2
194+ a, b = H[it, it], H[it + 1, it]
195+ r = np.sqrt(abs(a) ** 2 + abs(b) ** 2)
196+ if r == 0:
197+ cs[it], sn[it] = 1, 0
198+ else:
199+ cs[it] = abs(a) / r
200+ sn[it] = (a / abs(a)) * np.conj(b) / r if abs(a) > 0 else np.conj(b) / r
201+ H[it, it] = cs[it] * a + sn[it] * b
202+ H[it + 1, it] = 0
203+ svec[it + 1] = -np.conj(sn[it]) * svec[it]
204+ svec[it] = cs[it] * svec[it]
205+ it_done = it + 1
206+ if abs(svec[it + 1]) / rb < eps:
207+ break
208+ y = np.linalg.solve(np.triu(H[:it_done, :it_done]), svec[:it_done])
209+ x = cp.asarray(y) @ V[:it_done]
210+ resid = zid * x + apply_K(x) - rhs
211+ return x, it_done, wnorm(resid) / rb
fmm3dbie_gpu/surface.pyadded+182−0View file
@@ -0,0 +1,182 @@
1+"""Quadrilateral patch surfaces in the fmm3dbie convention, without the library.
2+
3+Each patch is the image of the reference square [-1, 1]^2. It carries
4+(norder+1)^2 tensor Gauss-Legendre nodes, with the u index varying fastest,
5+and an unnormalized Legendre basis P_j(u) P_i(v) with the same ordering.
6+Node data (srcvals) is xyz, dx/du, dx/dv, unit normal, and the geometry is
7+also stored as basis coefficients (srccoefs) so that it can be evaluated at
8+any point of the reference square. These are the same conventions as
9+fmm3dbie's iptype 11, and were checked against the library.
10+"""
11+
12+import numpy as np
13+
14+
15+def gauss_legendre(n):
16+ x, w = np.polynomial.legendre.leggauss(n)
17+ return x, w
18+
19+
20+def legendre_vander(x, norder):
21+ """P_0(x), ..., P_norder(x) as an array of shape (norder+1, len(x))."""
22+ x = np.asarray(x, dtype=float)
23+ P = np.zeros((norder + 1, x.size))
24+ P[0] = 1.0
25+ if norder >= 1:
26+ P[1] = x
27+ for n in range(1, norder):
28+ P[n + 1] = ((2 * n + 1) * x * P[n] - n * P[n - 1]) / (n + 1)
29+ return P
30+
31+
32+def disc_nodes(norder):
33+ """Reference nodes (2, npols) and weights (npols) of the tensor GL rule."""
34+ x, w = gauss_legendre(norder + 1)
35+ v, u = np.meshgrid(x, x, indexing="ij") # u fastest
36+ wv, wu = np.meshgrid(w, w, indexing="ij")
37+ return np.vstack([u.ravel(), v.ravel()]), (wu * wv).ravel()
38+
39+
40+def basis_2d(uv, norder):
41+ """Basis values (npols, npts): index n = i*(norder+1) + j is P_j(u) P_i(v)."""
42+ Pu = legendre_vander(uv[0], norder)
43+ Pv = legendre_vander(uv[1], norder)
44+ return (Pv[:, None, :] * Pu[None, :, :]).reshape((norder + 1) ** 2, -1)
45+
46+
47+def vals_to_coefs_matrix(norder):
48+ """U with coefs = U @ vals for values at the tensor GL nodes."""
49+ x, w = gauss_legendre(norder + 1)
50+ P = legendre_vander(x, norder)
51+ U1 = (2 * np.arange(norder + 1)[:, None] + 1) / 2 * P * w[None, :]
52+ return np.kron(U1, U1) # (i, j) x (ii, jj) with i, ii slow
53+
54+
55+def oversample_matrix(norder, nover):
56+ """M with values at the order-nover nodes = M @ values at the order-norder nodes."""
57+ uv_o, _ = disc_nodes(nover)
58+ return basis_2d(uv_o, norder).T @ vals_to_coefs_matrix(norder)
59+
60+
61+class Surface:
62+ """exact: optional dict(params=(npatches, 13), a=, p=) giving the superellipsoid face
63+ chart of every patch (n, e1, e2, s0, s1, t0, t1); when present the integrals use
64+ the exact geometry instead of the polynomial interpolant."""
65+
66+ def __init__(self, norder, srcvals, exact=None):
67+ self.exact = exact
68+ self.norder = norder
69+ self.iptype = 11
70+ self.srcvals = np.ascontiguousarray(srcvals, dtype=float) # (12, npts)
71+ self.npts = self.srcvals.shape[1]
72+ self.npols = (norder + 1) ** 2
73+ assert self.npts % self.npols == 0
74+ self.npatches = self.npts // self.npols
75+ self.uvs, self.wref = disc_nodes(norder)
76+ self.umat = vals_to_coefs_matrix(norder)
77+ vals = self.srcvals[:9].reshape(9, self.npatches, self.npols)
78+ self.srccoefs = np.einsum("nj,cpj->cpn", self.umat, vals).reshape(9, self.npts)
79+ ru, rv = self.srcvals[3:6], self.srcvals[6:9]
80+ self.jac = np.linalg.norm(np.cross(ru.T, rv.T).T, axis=0)
81+ self.wts = self.jac * np.tile(self.wref, self.npatches)
82+ # centroid and radius of each patch from a 5x5 sample grid (fmm3dbie uses the
83+ # four corners, which underestimates the extent of strongly curved patches)
84+ g5 = np.linspace(-1, 1, 5)
85+ corners = np.array(np.meshgrid(g5, g5, indexing="ij")).reshape(2, -1)
86+ Pc = basis_2d(corners, norder) # (npols, 25)
87+ xc = np.einsum("cpn,nk->cpk", self.srccoefs[:3].reshape(3, self.npatches, self.npols), Pc)
88+ self.centroids = xc.mean(axis=2) # (3, npatches)
89+ self.rads = np.linalg.norm(xc - self.centroids[:, :, None], axis=0).max(axis=1)
90+ self.patch_id = np.repeat(np.arange(self.npatches), self.npols)
91+ self.uv_nodes = np.tile(self.uvs, self.npatches) # (2, npts)
92+
93+ def eval_geometry(self, patch, uv):
94+ """xyz, r_u, r_v at reference points uv (2, n) of one patch: (9, n)."""
95+ if self.exact is not None:
96+ e = self.exact
97+ return superellipsoid_chart(e["params"][patch], e["a"], e["p"], uv)
98+ P = basis_2d(uv, self.norder)
99+ c = self.srccoefs[:, patch * self.npols:(patch + 1) * self.npols]
100+ return c @ P
101+
102+ def oversampled(self, nover):
103+ """Nodes of the order-nover rule on every patch: xyz, normal, weights."""
104+ if nover == self.norder:
105+ return self.srcvals[0:3], self.srcvals[9:12], self.wts
106+ uv_o, w_o = disc_nodes(nover)
107+ if self.exact is not None:
108+ g = np.hstack([self.eval_geometry(j, uv_o) for j in range(self.npatches)])
109+ else:
110+ P = basis_2d(uv_o, self.norder) # (npols, npols_o)
111+ g = np.einsum("cpn,nk->cpk", self.srccoefs.reshape(9, self.npatches, self.npols), P)
112+ g = g.reshape(9, -1)
113+ nrm = np.cross(g[3:6].T, g[6:9].T).T
114+ jac = np.linalg.norm(nrm, axis=0)
115+ nrm /= jac
116+ return g[0:3], nrm, jac * np.tile(w_o, self.npatches)
117+
118+
119+# ---------------------------------------------------------------------------
120+# Rounded cube (superellipsoid), the same construction as ../rounded_cube.py
121+
122+_FACES = [
123+ ((1, 0, 0), (0, 1, 0), (0, 0, 1)),
124+ ((-1, 0, 0), (0, 0, 1), (0, 1, 0)),
125+ ((0, 1, 0), (0, 0, 1), (1, 0, 0)),
126+ ((0, -1, 0), (1, 0, 0), (0, 0, 1)),
127+ ((0, 0, 1), (1, 0, 0), (0, 1, 0)),
128+ ((0, 0, -1), (0, 1, 0), (1, 0, 0)),
129+]
130+
131+
132+def graded_breaks(m, grading):
133+ t = np.linspace(-1, 1, m + 1)
134+ return (1 - grading) * t + grading * np.sin(0.5 * np.pi * t)
135+
136+
137+def superellipsoid_chart(par, a, p, uv):
138+ """xyz, r_u, r_v (9, n) of the chart par = (n, e1, e2, s0, s1, t0, t1) at uv (2, n)."""
139+ n, e1, e2 = par[0:3], par[3:6], par[6:9]
140+ s0, s1, t0, t1 = par[9:13]
141+ s = s0 + 0.5 * (uv[0] + 1) * (s1 - s0)
142+ t = t0 + 0.5 * (uv[1] + 1) * (t1 - t0)
143+ c = n[:, None] + s * e1[:, None] + t * e2[:, None]
144+ N = np.sum(c**p, axis=0) ** (1.0 / p)
145+ gradN = (c / N) ** (p - 1)
146+
147+ def deriv(dc, scale):
148+ dN = np.sum(gradN * dc[:, None], axis=0)
149+ return a * scale * (dc[:, None] / N - c * dN / N**2)
150+
151+ return np.vstack([a * c / N, deriv(e1, 0.5 * (s1 - s0)), deriv(e2, 0.5 * (t1 - t0))])
152+
153+
154+def rounded_cube(a=1.0, p=8, m=6, norder=8, grading=0.5, exact=True):
155+ """Superellipsoid |x|^p+|y|^p+|z|^p = a^p, m x m patches per face.
156+
157+ With exact=True the integrals use the analytic parametrization; otherwise the
158+ polynomial interpolant of the node data, as fmm3dbie does."""
159+ uvs, _ = disc_nodes(norder)
160+ br = graded_breaks(m, grading)
161+ blocks, params = [], []
162+ for n, e1, e2 in _FACES:
163+ for i in range(m):
164+ for j in range(m):
165+ par = np.array(list(n) + list(e1) + list(e2) + [br[i], br[i + 1], br[j], br[j + 1]],
166+ dtype=float)
167+ params.append(par)
168+ g = superellipsoid_chart(par, a, p, uvs)
169+ nrm = np.cross(g[3:6].T, g[6:9].T).T
170+ nrm /= np.linalg.norm(nrm, axis=0)
171+ blocks.append(np.vstack([g, nrm]))
172+ ex = dict(params=np.array(params), a=a, p=p) if exact else None
173+ return Surface(norder, np.hstack(blocks), exact=ex)
174+
175+
176+def superellipsoid_volume(a, p):
177+ from math import gamma
178+ return 8 * a**3 * gamma(1 + 1 / p) ** 3 / gamma(1 + 3 / p)
179+
180+
181+def inside(x, a, p, margin=0.0):
182+ return np.sum(np.abs(x / a) ** p, axis=0) < (1 - margin) ** p
plot.pyadded+98−0View file
@@ -0,0 +1,98 @@
1+"""Plot the room field saved by run.py (python plot.py results/<case>-<backend>): slices of Re(u) and level in dB,
2+plus an animation of the real pressure field Re(u exp(-i omega t)) over one period."""
3+
4+import argparse
5+import os
6+
7+import matplotlib
8+
9+matplotlib.use("Agg")
10+import matplotlib.pyplot as plt
11+import numpy as np
12+from matplotlib.animation import FuncAnimation, PillowWriter
13+
14+SLICES = [("xy", "x", "y", 2), ("xz", "x", "z", 1), ("yz", "y", "z", 0)]
15+
16+
17+def masked(a):
18+ return np.ma.masked_invalid(a)
19+
20+
21+def main():
22+ ap = argparse.ArgumentParser()
23+ ap.add_argument("out", help="a results directory written by run.py")
24+ ap.add_argument("--frames", type=int, default=32)
25+ args = ap.parse_args()
26+
27+ d = np.load(os.path.join(args.out, "field.npz"))
28+ g, x0, k = d["grid"], d["x0"], float(d["k"])
29+ ext = [g[0], g[-1], g[0], g[-1]]
30+ us = {n: d[f"u_{n}"] for n, *_ in SLICES}
31+
32+ # Color limits from the field more than half a wavelength from the source
33+ A, B = np.meshgrid(g, g, indexing="xy")
34+ far = []
35+ for n, lx, ly, fixed in SLICES:
36+ r = np.hypot(A - x0["xyz".index(lx)], B - x0["xyz".index(ly)])
37+ u = us[n]
38+ far.append(np.abs(u[np.isfinite(u) & (r > np.pi / k)]))
39+ finite = np.concatenate(far)
40+ umax = np.percentile(finite, 99)
41+ ref = np.median(finite)
42+
43+ fig, axes = plt.subplots(2, 3, figsize=(13, 8.4), constrained_layout=True)
44+ for col, (n, lx, ly, fixed) in enumerate(SLICES):
45+ u = us[n]
46+ ax = axes[0, col]
47+ im0 = ax.imshow(masked(u.real), origin="lower", extent=ext, cmap="RdBu_r",
48+ vmin=-umax, vmax=umax, interpolation="bilinear")
49+ ax.set_title(f"Re u, {n}-plane at {'xyz'[fixed]} = {x0[fixed]:.2f}")
50+ ax = axes[1, col]
51+ with np.errstate(divide="ignore", invalid="ignore"):
52+ level = 20 * np.log10(np.abs(u) / ref)
53+ im1 = ax.imshow(masked(level), origin="lower", extent=ext, cmap="Blues_r",
54+ vmin=-30, vmax=15, interpolation="bilinear")
55+ ax.set_title(f"level (dB re median), {n}-plane")
56+ for ax in axes[:, col]:
57+ ax.plot(x0["xyz".index(lx)], x0["xyz".index(ly)], marker="o", ms=5,
58+ mfc="none", mec="0.15", mew=1)
59+ ax.set_xlabel(lx)
60+ ax.set_ylabel(ly)
61+ ax.set_aspect("equal")
62+ fig.colorbar(im0, ax=axes[0, :], shrink=0.85, label="Re u")
63+ fig.colorbar(im1, ax=axes[1, :], shrink=0.85, label="dB")
64+ wl = 2 * np.pi / k
65+ fig.suptitle(f"Point source in a rounded cubical room: k = {k:g} "
66+ f"(room width = {2 * float(d['a']) / wl:.2f} wavelengths), "
67+ f"wall admittance {-float(d['lam']):g}")
68+ fig.savefig(os.path.join(args.out, "slices.png"), dpi=130)
69+ plt.close(fig)
70+
71+ # Animation: Re(u e^{-i w t}) over one period on the three planes
72+ fig, axes = plt.subplots(1, 3, figsize=(13, 4.6), constrained_layout=True)
73+ ims = []
74+ for ax, (n, lx, ly, fixed) in zip(axes, SLICES):
75+ im = ax.imshow(masked(us[n].real), origin="lower", extent=ext, cmap="RdBu_r",
76+ vmin=-umax, vmax=umax, interpolation="bilinear")
77+ ax.set_title(f"{n}-plane")
78+ ax.set_xlabel(lx)
79+ ax.set_ylabel(ly)
80+ ims.append(im)
81+ fig.colorbar(ims[-1], ax=axes, shrink=0.85, label="Re(u e^{-iωt})")
82+ title = fig.suptitle("")
83+
84+ def frame(i):
85+ ph = np.exp(-2j * np.pi * i / args.frames)
86+ for im, (n, *_) in zip(ims, SLICES):
87+ im.set_data(masked((us[n] * ph).real))
88+ title.set_text(f"pressure over one period, t/T = {i / args.frames:.2f}")
89+ return ims
90+
91+ anim = FuncAnimation(fig, frame, frames=args.frames, blit=False)
92+ anim.save(os.path.join(args.out, "animation.gif"), writer=PillowWriter(fps=12), dpi=80)
93+ plt.close(fig)
94+ print(f"wrote {args.out}/slices.png and {args.out}/animation.gif")
95+
96+
97+if __name__ == "__main__":
98+ main()
problem.pyadded+114−0View file
@@ -0,0 +1,114 @@
1+"""The showcase problem, shared by both backends.
2+
3+A point source at x0 inside a rounded cube (the superellipsoid
4+|x|^p + |y|^p + |z|^p = a^p) with impedance walls. With time dependence
5+exp(-i omega t) the pressure amplitude u satisfies
6+
7+ (Delta + k^2) u = -delta(x - x0) in the room,
8+ du/dn + i k lambda u = 0 on the walls,
9+
10+with n the outward normal and lambda = -beta, where beta > 0 is the wall
11+admittance. We write u = G + v with G the free-space field of the source, and
12+find v from the right-preconditioned combined field integral equation of
13+fmm3dbie (helm_rpcomb_imp_solver, interior mode).
14+
15+This module has no dependence on either backend.
16+"""
17+
18+import json
19+import os
20+import platform
21+import socket
22+import subprocess
23+
24+import numpy as np
25+
26+# Named cases. Both backends discretize the surface identically, so the same
27+# case gives the same unknowns on both.
28+CASES = {
29+ "small": dict(k=10.0, beta=0.2, m=4, norder=6), # 3.2 wavelengths, 4,704 points
30+ "medium": dict(k=10.0, beta=0.2, m=6, norder=8), # 3.2 wavelengths, 17,496 points
31+ "large": dict(k=40.0, beta=0.05, m=10, norder=8), # 12.7 wavelengths, 48,600 points
32+}
33+A, P = 1.0, 8 # room half-width and superellipsoid exponent
34+SOURCE = np.array([-0.35, 0.2, 0.1]) # the point source
35+EXTERIOR_SOURCE = np.array([1.7, 0.6, -0.5]) # for the self-test (outside the room)
36+ALPHA = 1.0
37+
38+
39+def green(x, x0, k):
40+ d = x - x0[:, None]
41+ r = np.linalg.norm(d, axis=0)
42+ return np.exp(1j * k * r) / (4 * np.pi * r)
43+
44+
45+def green_dn(x, nrm, x0, k):
46+ d = x - x0[:, None]
47+ r = np.linalg.norm(d, axis=0)
48+ g = np.exp(1j * k * r) / (4 * np.pi * r)
49+ return g * (1j * k - 1 / r) * np.sum(d * nrm, axis=0) / r
50+
51+
52+def inside(x, a=A, p=P, margin=0.0):
53+ return np.sum(np.abs(x / a) ** p, axis=0) < (1 - margin) ** p
54+
55+
56+def room_rhs(xs, nrm, k, lam):
57+ """Impedance data for v = u - G: -(dG/dn + i k lambda G)."""
58+ return -(green_dn(xs, nrm, SOURCE, k) + 1j * k * lam * green(xs, SOURCE, k))
59+
60+
61+def selftest_rhs(xs, nrm, k, lam):
62+ """Impedance data of the field of a source outside the room."""
63+ return green_dn(xs, nrm, EXTERIOR_SOURCE, k) + 1j * k * lam * green(xs, EXTERIOR_SOURCE, k)
64+
65+
66+def selftest_points(n=50):
67+ return np.random.default_rng(0).uniform(-0.8 * A, 0.8 * A, size=(3, n))
68+
69+
70+def slice_points(ngrid):
71+ """Points on three planes through the source; returns grid, {name: (pts, mask, shape)}."""
72+ g = np.linspace(-A, A, ngrid)
73+ out = {}
74+ for name, (i, j, fixed) in {"xy": (0, 1, 2), "xz": (0, 2, 1), "yz": (1, 2, 0)}.items():
75+ Ag, Bg = np.meshgrid(g, g, indexing="xy")
76+ pts = np.zeros((3, Ag.size))
77+ pts[i], pts[j], pts[fixed] = Ag.ravel(), Bg.ravel(), SOURCE[fixed]
78+ out[name] = (pts, inside(pts, margin=1e-3), Ag.shape)
79+ return g, out
80+
81+
82+def energy_balance(k, lam, v_at_source, wts, u_wall):
83+ """Power radiated by the source, k/(4 pi) + Im v(x0), and power absorbed by the
84+ walls, -k Re(lambda) int |u|^2 dS. They must agree."""
85+ p_src = k / (4 * np.pi) + v_at_source.imag
86+ p_wall = -k * np.real(lam) * np.sum(wts * np.abs(u_wall) ** 2)
87+ return p_src, p_wall
88+
89+
90+def machine_info():
91+ info = dict(host=socket.gethostname(), platform=platform.platform(), cpu_count=os.cpu_count())
92+ try:
93+ with open("/proc/cpuinfo") as f:
94+ for line in f:
95+ if line.startswith("model name"):
96+ info["cpu"] = line.split(":", 1)[1].strip()
97+ break
98+ except OSError:
99+ pass
100+ try:
101+ gpu = subprocess.run(["nvidia-smi", "--query-gpu=name,memory.total", "--format=csv,noheader"],
102+ capture_output=True, text=True, timeout=10)
103+ if gpu.returncode == 0:
104+ info["gpu"] = gpu.stdout.strip().splitlines()[0]
105+ except (OSError, subprocess.SubprocessError):
106+ pass
107+ return info
108+
109+
110+def save(outdir, summary, **arrays):
111+ os.makedirs(outdir, exist_ok=True)
112+ np.savez(os.path.join(outdir, "field.npz"), **arrays)
113+ with open(os.path.join(outdir, "summary.json"), "w") as f:
114+ json.dump(summary, f, indent=2)
reference/fmm3dbie_ctypes.pyadded+168−0View file
@@ -0,0 +1,168 @@
1+"""Minimal ctypes interface to the fmm3dbie shared library.
2+
3+Only the routines needed for the interior Helmholtz impedance problem are
4+wrapped. All Fortran arguments are passed by reference; integers are
5+integer*8, and arrays are Fortran (column-major) ordered.
6+"""
7+
8+import contextlib
9+import ctypes
10+import os
11+import sys
12+
13+# OpenBLAS threads compete with the library's OpenMP threads
14+os.environ.setdefault("OPENBLAS_NUM_THREADS", "1")
15+# write Fortran console output immediately so it can be redirected to a log
16+os.environ.setdefault("GFORTRAN_UNBUFFERED_PRECONNECTED", "y")
17+
18+import numpy as np
19+
20+LIBPATH = os.environ.get(
21+ "FMM3DBIE_LIB", os.path.expanduser("~/src/fmm3dbie/lib/libfmm3dbie.so")
22+)
23+_lib = ctypes.CDLL(LIBPATH)
24+
25+_i8 = ctypes.c_int64
26+_f8 = ctypes.c_double
27+
28+
29+@contextlib.contextmanager
30+def fortran_output_to(path):
31+ """Redirect the library's console output (Fortran unit 6, i.e. fd 1) to a file."""
32+ sys.stdout.flush()
33+ saved = os.dup(1)
34+ with open(path, "a") as f:
35+ os.dup2(f.fileno(), 1)
36+ try:
37+ yield
38+ finally:
39+ _libc.fflush(None)
40+ os.dup2(saved, 1)
41+ os.close(saved)
42+
43+
44+_libc = ctypes.CDLL(None)
45+
46+
47+def _ref_i8(v):
48+ return ctypes.byref(_i8(int(v)))
49+
50+
51+def _ref_f8(v):
52+ return ctypes.byref(_f8(float(v)))
53+
54+
55+def _ptr(a):
56+ return a.ctypes.data_as(ctypes.c_void_p)
57+
58+
59+def _farr(a, dtype):
60+ return np.asfortranarray(a, dtype=dtype)
61+
62+
63+class Surface:
64+ """Discretized surface in the fmm3dbie format (all patches share one order and type)."""
65+
66+ def __init__(self, norder, iptype, srcvals):
67+ self.norder = norder
68+ self.iptype = iptype
69+ self.srcvals = _farr(srcvals, np.float64) # (12, npts)
70+ self.npts = self.srcvals.shape[1]
71+ self.npols = npols_for(norder, iptype)
72+ self.npatches = self.npts // self.npols
73+ self.norders = np.full(self.npatches, norder, dtype=np.int64)
74+ self.iptypes = np.full(self.npatches, iptype, dtype=np.int64)
75+ self.ixyzs = (np.arange(self.npatches + 1, dtype=np.int64) * self.npols + 1)
76+ self.srccoefs = surf_vals_to_coefs(self, self.srcvals[0:9, :])
77+ self.wts = get_qwts(self)
78+
79+ def args(self):
80+ return (_ref_i8(self.npatches), _ptr(self.norders), _ptr(self.ixyzs),
81+ _ptr(self.iptypes), _ref_i8(self.npts), _ptr(self.srccoefs),
82+ _ptr(self.srcvals))
83+
84+
85+def npols_for(norder, iptype):
86+ if iptype == 1:
87+ return (norder + 1) * (norder + 2) // 2
88+ return (norder + 1) ** 2
89+
90+
91+def get_disc_nodes(norder, iptype):
92+ npols = npols_for(norder, iptype)
93+ uvs = np.zeros((2, npols), order="F")
94+ _lib.get_disc_nodes_(_ref_i8(norder), _ref_i8(npols), _ref_i8(iptype), _ptr(uvs))
95+ return uvs
96+
97+
98+def surf_vals_to_coefs(surf, vals):
99+ vals = _farr(vals, np.float64)
100+ nd = vals.shape[0]
101+ coefs = np.zeros_like(vals, order="F")
102+ _lib.surf_vals_to_coefs_(_ref_i8(nd), _ref_i8(surf.npatches), _ptr(surf.norders),
103+ _ptr(surf.ixyzs), _ptr(surf.iptypes), _ref_i8(surf.npts),
104+ _ptr(vals), _ptr(coefs))
105+ return coefs
106+
107+
108+def get_qwts(surf):
109+ wts = np.zeros(surf.npts)
110+ _lib.get_qwts_(_ref_i8(surf.npatches), _ptr(surf.norders), _ptr(surf.ixyzs),
111+ _ptr(surf.iptypes), _ref_i8(surf.npts), _ptr(surf.srcvals), _ptr(wts))
112+ return wts
113+
114+
115+def helm_rpcomb_imp_solver(surf, zk, alpha, zlams, rhs, ifinout, eps=1e-6,
116+ eps_gmres=1e-6, numit=400):
117+ """Solve n.grad(u) + i k lambda u = rhs with u = S_k[s] + i alpha D_k S_{i|k|}[s].
118+
119+ ifinout = 0 for the interior problem, 1 for the exterior problem.
120+ Returns (sigma, siksigma, niter, errs, rres).
121+ """
122+ zpars = np.array([zk, alpha], dtype=np.complex128)
123+ zlams = _farr(zlams, np.complex128)
124+ rhs = _farr(rhs, np.complex128)
125+ sigma = np.zeros(surf.npts, dtype=np.complex128)
126+ siksigma = np.zeros(surf.npts, dtype=np.complex128)
127+ errs = np.zeros(numit + 1)
128+ niter = _i8(0)
129+ rres = _f8(0.0)
130+ _lib.helm_rpcomb_imp_solver_(*surf.args(), _ref_f8(eps), _ptr(zpars), _ptr(zlams),
131+ _ref_i8(numit), _ref_i8(ifinout), _ptr(rhs),
132+ _ref_f8(eps_gmres), ctypes.byref(niter), _ptr(errs),
133+ ctypes.byref(rres), _ptr(sigma), _ptr(siksigma))
134+ return sigma, siksigma, niter.value, errs[: niter.value], rres.value
135+
136+
137+def helm_rpcomb_eval(surf, zk, alpha, sigma, siksigma, targs, eps=1e-6):
138+ """Evaluate u = S_k[s] + i alpha D_k S_{i|k|}[s] at off-surface targets (3, ntarg)."""
139+ targs = _farr(targs, np.float64)
140+ ntarg = targs.shape[1]
141+ zpars = np.array([zk, alpha], dtype=np.complex128)
142+ ipatch_id = np.full(ntarg, -1, dtype=np.int64)
143+ uvs_targ = np.zeros((2, ntarg), order="F")
144+ pot = np.zeros(ntarg, dtype=np.complex128)
145+ _lib.helm_rpcomb_eval_(*surf.args(), _ref_i8(3), _ref_i8(ntarg), _ptr(targs),
146+ _ptr(ipatch_id), _ptr(uvs_targ), _ref_f8(eps), _ptr(zpars),
147+ _ptr(_farr(sigma, np.complex128)),
148+ _ptr(_farr(siksigma, np.complex128)), _ptr(pot))
149+ return pot
150+
151+
152+def helm_rpcomb_eval_boundary(surf, zk, alpha, sigma, siksigma, eps=1e-6):
153+ """Interior boundary limit of u = S_k[s] + i alpha D_k S_{i|k|}[s] at the nodes.
154+
155+ The on-surface evaluation gives the principal value of D_k; with outward
156+ normals the interior limit adds -1/2 times the double-layer density.
157+ """
158+ uvs = get_disc_nodes(surf.norder, surf.iptype)
159+ ipatch_id = np.repeat(np.arange(1, surf.npatches + 1), surf.npols).astype(np.int64)
160+ uvs_targ = np.asfortranarray(np.tile(uvs, (1, surf.npatches)))
161+ zpars = np.array([zk, alpha], dtype=np.complex128)
162+ sigma = _farr(sigma, np.complex128)
163+ siksigma = _farr(siksigma, np.complex128)
164+ pot = np.zeros(surf.npts, dtype=np.complex128)
165+ _lib.helm_rpcomb_eval_(*surf.args(), _ref_i8(12), _ref_i8(surf.npts), _ptr(surf.srcvals),
166+ _ptr(ipatch_id), _ptr(uvs_targ), _ref_f8(eps), _ptr(zpars),
167+ _ptr(sigma), _ptr(siksigma), _ptr(pot))
168+ return pot - 0.5j * alpha * siksigma
reference/rounded_cube.pyadded+69−0View file
@@ -0,0 +1,69 @@
1+"""Rounded cube (superellipsoid) surface built from quadrilateral patches.
2+
3+The surface is |x|^p + |y|^p + |z|^p = a^p with p even. Each of the six cube
4+faces (s, t) in [-1, 1]^2 is mapped radially onto the surface,
5+x = a c / ||c||_p with c = n + s e1 + t e2. Larger p gives sharper edges.
6+Patch breakpoints in s and t are graded toward the face edges, where the
7+curvature is concentrated.
8+"""
9+
10+import numpy as np
11+
12+import fmm3dbie_ctypes as fb
13+
14+# (n, e1, e2) with e1 x e2 = n, so that d/ds x d/dt points outward
15+_FACES = [
16+ ((1, 0, 0), (0, 1, 0), (0, 0, 1)),
17+ ((-1, 0, 0), (0, 0, 1), (0, 1, 0)),
18+ ((0, 1, 0), (0, 0, 1), (1, 0, 0)),
19+ ((0, -1, 0), (1, 0, 0), (0, 0, 1)),
20+ ((0, 0, 1), (1, 0, 0), (0, 1, 0)),
21+ ((0, 0, -1), (0, 1, 0), (1, 0, 0)),
22+]
23+
24+
25+def graded_breaks(m, grading):
26+ """m+1 breakpoints on [-1, 1]; grading in [0, 1] blends uniform and sine-clustered."""
27+ t = np.linspace(-1, 1, m + 1)
28+ return (1 - grading) * t + grading * np.sin(0.5 * np.pi * t)
29+
30+
31+def rounded_cube(a=1.0, p=8, m=6, norder=8, grading=0.5):
32+ """Return an fmm3dbie Surface for the superellipsoid of half-width a.
33+
34+ m x m patches per face, each discretized with (norder+1)^2 Gauss-Legendre nodes.
35+ """
36+ iptype = 11
37+ uvs = fb.get_disc_nodes(norder, iptype) # (2, npols) on [-1, 1]^2
38+ br = graded_breaks(m, grading)
39+ blocks = []
40+ for n, e1, e2 in _FACES:
41+ n, e1, e2 = (np.array(v, dtype=float) for v in (n, e1, e2))
42+ for i in range(m):
43+ for j in range(m):
44+ s0, s1 = br[i], br[i + 1]
45+ t0, t1 = br[j], br[j + 1]
46+ s = s0 + 0.5 * (uvs[0] + 1) * (s1 - s0)
47+ t = t0 + 0.5 * (uvs[1] + 1) * (t1 - t0)
48+ c = n[:, None] + s * e1[:, None] + t * e2[:, None] # (3, npols)
49+ N = np.sum(c**p, axis=0) ** (1.0 / p)
50+ gradN = c ** (p - 1) / N ** (p - 1) # dN/dc
51+ xyz = a * c / N
52+
53+ def deriv(dc, scale):
54+ dN = np.sum(gradN * dc[:, None], axis=0)
55+ return a * scale * (dc[:, None] / N - c * dN / N**2)
56+
57+ du = deriv(e1, 0.5 * (s1 - s0))
58+ dv = deriv(e2, 0.5 * (t1 - t0))
59+ nrm = np.cross(du.T, dv.T).T
60+ nrm /= np.linalg.norm(nrm, axis=0)
61+ blocks.append(np.vstack([xyz, du, dv, nrm]))
62+ srcvals = np.hstack(blocks)
63+ return fb.Surface(norder, iptype, srcvals)
64+
65+
66+def superellipsoid_volume(a, p):
67+ """Exact volume of |x|^p+|y|^p+|z|^p <= a^p (used to check the discretization)."""
68+ from math import gamma
69+ return 8 * a**3 * gamma(1 + 1 / p) ** 3 / gamma(1 + 3 / p)
run.pyadded+188−0View file
@@ -0,0 +1,188 @@
1+"""Solve the showcase problem with one backend and save timings and fields.
2+
3+ python run.py --backend gpu --case small --selftest
4+ python run.py --backend fmm3dbie --case small --selftest
5+
6+Results go to results/<case>-<backend>/ (summary.json, field.npz). The two
7+backends solve the same discrete problem on the same surface nodes; compare.py
8+puts their results side by side.
9+"""
10+
11+import argparse
12+import os
13+import sys
14+import time
15+
16+import numpy as np
17+
18+HERE = os.path.dirname(os.path.abspath(__file__))
19+sys.path.insert(0, HERE)
20+
21+import problem as pb # noqa: E402
22+
23+
24+class GPUBackend:
25+ """The CuPy/CUDA solver in fmm3dbie_gpu."""
26+
27+ name = "gpu"
28+
29+ def __init__(self, case, args):
30+ import cupy as cp
31+ from fmm3dbie_gpu import cuda_ops, solver, surface
32+ self.cp = cp
33+ t0 = time.time()
34+ cuda_ops.module() # compile or load the CUDA kernels once, outside the timings
35+ self.kernel_load_s = time.time() - t0
36+ self.surf = surface.rounded_cube(a=pb.A, p=pb.P, m=case["m"], norder=case["norder"],
37+ exact=not args.interp_geometry)
38+ self.lam = -case["beta"]
39+ t0 = time.time()
40+ self.op = solver.Operator(self.surf, case["k"], pb.ALPHA, self.lam, ifinout=0,
41+ far_dtype=np.float64 if args.far64 else np.float32, verbose=False)
42+ cp.cuda.Device().synchronize()
43+ self.setup_s = time.time() - t0
44+ self.eps_gmres = args.eps_gmres
45+ self.options = dict(far_field="float64" if args.far64 else "float32",
46+ geometry="interpolated" if args.interp_geometry else "exact")
47+
48+ def solve(self, rhs):
49+ sigma, phi1, niter, rres = self.op.solve(rhs, eps_gmres=self.eps_gmres, verbose=False)
50+ return (sigma, phi1), niter, rres
51+
52+ def eval_points(self, state, pts):
53+ return self.op.eval(pts.T, *state).get()
54+
55+ def eval_wall(self, state):
56+ return self.op.eval(None, *state, on_surface=True).get()
57+
58+ def sync(self):
59+ self.cp.cuda.Device().synchronize()
60+
61+
62+class Fmm3dbieBackend:
63+ """fmm3dbie on the CPU (helm_rpcomb_imp_solver), through ctypes."""
64+
65+ name = "fmm3dbie"
66+
67+ def __init__(self, case, args):
68+ sys.path.insert(0, os.path.join(HERE, "reference"))
69+ import fmm3dbie_ctypes as fb
70+ import rounded_cube as rc
71+ self.fb = fb
72+ self.log = os.path.join(args.outdir, "fortran.log")
73+ os.makedirs(args.outdir, exist_ok=True)
74+ open(self.log, "w").close()
75+ self.kernel_load_s = 0.0
76+ with fb.fortran_output_to(self.log):
77+ self.surf = rc.rounded_cube(a=pb.A, p=pb.P, m=case["m"], norder=case["norder"])
78+ self.k, self.lam = case["k"], -case["beta"]
79+ self.setup_s = 0.0 # the near quadrature is computed inside each solver call
80+ self.eps, self.eps_gmres = args.eps, args.eps_gmres
81+ self.options = dict(eps=args.eps, threads=os.environ.get("OMP_NUM_THREADS", "all"))
82+
83+ def solve(self, rhs):
84+ with self.fb.fortran_output_to(self.log):
85+ sigma, sik, niter, errs, rres = self.fb.helm_rpcomb_imp_solver(
86+ self.surf, self.k, pb.ALPHA, self.lam * np.ones(self.surf.npts), rhs, ifinout=0,
87+ eps=self.eps, eps_gmres=self.eps_gmres)
88+ return (sigma, sik), int(niter), float(rres)
89+
90+ def eval_points(self, state, pts):
91+ with self.fb.fortran_output_to(self.log):
92+ return self.fb.helm_rpcomb_eval(self.surf, self.k, pb.ALPHA, *state, pts, eps=self.eps)
93+
94+ def eval_wall(self, state):
95+ with self.fb.fortran_output_to(self.log):
96+ return self.fb.helm_rpcomb_eval_boundary(self.surf, self.k, pb.ALPHA, *state, eps=self.eps)
97+
98+ def sync(self):
99+ pass
100+
101+
102+def timed(be, f, *a):
103+ t0 = time.time()
104+ out = f(*a)
105+ be.sync()
106+ return out, time.time() - t0
107+
108+
109+def main():
110+ ap = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter)
111+ ap.add_argument("--backend", choices=["gpu", "fmm3dbie"], required=True)
112+ ap.add_argument("--case", choices=list(pb.CASES), default="small")
113+ ap.add_argument("--k", type=float, help="override the case's wavenumber")
114+ ap.add_argument("--beta", type=float, help="override the case's wall admittance")
115+ ap.add_argument("--m", type=int, help="override the patches per face edge")
116+ ap.add_argument("--norder", type=int, help="override the order per patch")
117+ ap.add_argument("--selftest", action="store_true",
118+ help="also solve a problem with a known solution (a second solve)")
119+ ap.add_argument("--eps-gmres", type=float, default=1e-6)
120+ ap.add_argument("--eps", type=float, default=1e-4, help="fmm3dbie quadrature and FMM tolerance")
121+ ap.add_argument("--far64", action="store_true", help="gpu: double precision far field")
122+ ap.add_argument("--interp-geometry", action="store_true",
123+ help="gpu: interpolated patch geometry inside the integrals, as fmm3dbie")
124+ ap.add_argument("--ngrid", type=int, default=101, help="grid size of the field slices (0 to skip)")
125+ ap.add_argument("--outdir", help="default results/<case>-<backend>")
126+ args = ap.parse_args()
127+
128+ case = dict(pb.CASES[args.case])
129+ for key in ["k", "beta", "m", "norder"]:
130+ if getattr(args, key) is not None:
131+ case[key] = getattr(args, key)
132+ args.outdir = args.outdir or os.path.join(HERE, "results", f"{args.case}-{args.backend}")
133+ k, lam = case["k"], -case["beta"]
134+
135+ be = (GPUBackend if args.backend == "gpu" else Fmm3dbieBackend)(case, args)
136+ surf = be.surf
137+ xs, nrm = np.asarray(surf.srcvals)[0:3], np.asarray(surf.srcvals)[9:12]
138+ wts = np.asarray(surf.wts)
139+ print(f"[{args.backend}] case {args.case}: k = {k}, beta = {case['beta']}, "
140+ f"{surf.npatches} patches of order {case['norder']}, {surf.npts} points, "
141+ f"room width {2 * pb.A * k / (2 * np.pi):.2f} wavelengths")
142+ if args.backend == "gpu":
143+ print(f" kernel load {be.kernel_load_s:.2f} s, operator setup {be.setup_s:.2f} s")
144+
145+ summary = dict(backend=args.backend, case=args.case, **case, npts=int(surf.npts),
146+ eps_gmres=args.eps_gmres, options=be.options, machine=pb.machine_info(),
147+ kernel_load_s=be.kernel_load_s, setup_s=be.setup_s)
148+
149+ if args.selftest:
150+ (state, niter, rres), t = timed(be, be.solve, pb.selftest_rhs(xs, nrm, k, lam))
151+ xt = pb.selftest_points()
152+ ut = be.eval_points(state, xt)
153+ ub = be.eval_wall(state)
154+ uex, ubex = pb.green(xt, pb.EXTERIOR_SOURCE, k), pb.green(xs, pb.EXTERIOR_SOURCE, k)
155+ err = float(np.linalg.norm(ut - uex) / np.linalg.norm(uex))
156+ err_wall = float(np.linalg.norm(ub - ubex) / np.linalg.norm(ubex))
157+ print(f" self-test: {niter} iterations, {t:.2f} s; relative error {err:.2e} at 50 interior "
158+ f"points, {err_wall:.2e} on the walls")
159+ summary.update(selftest_err=err, selftest_err_wall=err_wall, selftest_solve_s=t)
160+
161+ (state, niter, rres), t_solve = timed(be, be.solve, pb.room_rhs(xs, nrm, k, lam))
162+ print(f" room problem: {niter} GMRES iterations, residual {rres:.2e}, solve {t_solve:.2f} s")
163+ v0 = complex(be.eval_points(state, pb.SOURCE[:, None])[0])
164+ u_wall = pb.green(xs, pb.SOURCE, k) + be.eval_wall(state)
165+ p_src, p_wall = pb.energy_balance(k, lam, v0, wts, u_wall)
166+ mismatch = abs(p_src - p_wall) / p_src
167+ print(f" energy balance: source {p_src:.8f}, walls {p_wall:.8f}, relative mismatch {mismatch:.2e}")
168+ summary.update(niter=niter, rres=rres, solve_s=t_solve, setup_plus_solve_s=be.setup_s + t_solve,
169+ p_src=p_src, p_wall=p_wall, energy_mismatch=mismatch)
170+
171+ arrays = dict(x0=pb.SOURCE, k=k, lam=lam, a=pb.A, p=pb.P, xs=xs, wts=wts, u_wall=u_wall)
172+ if args.ngrid > 0:
173+ g, sl = pb.slice_points(args.ngrid)
174+ t0 = time.time()
175+ for name, (pts, mask, shape) in sl.items():
176+ u = np.full(mask.size, np.nan + 1j * np.nan)
177+ u[mask] = pb.green(pts[:, mask], pb.SOURCE, k) + be.eval_points(state, pts[:, mask])
178+ arrays[f"u_{name}"] = u.reshape(shape)
179+ be.sync()
180+ summary["slices_s"] = time.time() - t0
181+ arrays["grid"] = g
182+ print(f" field on three {args.ngrid} x {args.ngrid} slices: {summary['slices_s']:.2f} s")
183+ pb.save(args.outdir, summary, **arrays)
184+ print(f" saved {args.outdir}")
185+
186+
187+if __name__ == "__main__":
188+ main()
scripts/build_fmm3dbie.shadded+28−0View file
@@ -0,0 +1,28 @@
1+#!/bin/bash
2+# Build the fmm3dbie shared library that the reference backend loads through ctypes.
3+# Usage: scripts/build_fmm3dbie.sh [install_dir] (default ~/src/fmm3dbie)
4+# Needs git, gfortran, make and an OpenBLAS (or other BLAS/LAPACK) library.
5+set -euo pipefail
6+DEST="${1:-$HOME/src/fmm3dbie}"
7+# The commit used for the results in the README; set FMM3DBIE_COMMIT= (empty) for the latest.
8+COMMIT="${FMM3DBIE_COMMIT-94cb732e1660f74864ef0f527ef713d308453d09}"
9+if [ ! -d "$DEST/.git" ]; then
10+ git clone https://github.com/fastalgorithms/fmm3dbie "$DEST"
11+fi
12+cd "$DEST"
13+if [ -n "$COMMIT" ]; then git checkout -q "$COMMIT"; fi
14+git submodule update --init
15+if [ ! -f make.inc ]; then
16+ cp make.inc.linux.gnu.openblas make.inc
17+ # Some systems ship only the versioned runtime library (libopenblas.so.0) without the
18+ # development link; link against it directly in that case.
19+ if ! ld -lopenblas -o /dev/null --unresolved-symbols=ignore-all 2>/dev/null; then
20+ for lib in /usr/lib64/libopenblas.so.0 /usr/lib/x86_64-linux-gnu/libopenblas.so.0; do
21+ if [ -f "$lib" ]; then sed -i "s|^LBLAS=.*|LBLAS=-l:libopenblas.so.0|" make.inc; break; fi
22+ done
23+ fi
24+fi
25+make lib -j"$(nproc)"
26+echo
27+echo "Built $DEST/lib/libfmm3dbie.so. If it is not at ~/src/fmm3dbie, set"
28+echo " export FMM3DBIE_LIB=$DEST/lib/libfmm3dbie.so"
scripts/gpuenv.shadded+10−0View file
@@ -0,0 +1,10 @@
1+# Source this file: points the loader at the CUDA libraries installed as pip wheels
2+# by scripts/setup_gpu.sh. Not needed when a CUDA 12 toolkit is installed system wide.
3+NV="$(python3 -c 'import nvidia; print(list(nvidia.__path__)[0])' 2>/dev/null || true)"
4+if [ -n "$NV" ]; then
5+ for d in cuda_nvrtc cuda_runtime cublas cusparse cufft curand cusolver nvjitlink; do
6+ LD_LIBRARY_PATH="$NV/$d/lib${LD_LIBRARY_PATH:+:$LD_LIBRARY_PATH}"
7+ done
8+ export LD_LIBRARY_PATH
9+fi
10+export CUPY_CACHE_DIR="${CUPY_CACHE_DIR:-$HOME/.cupy/kernel_cache}"
scripts/setup_gpu.shadded+8−0View file
@@ -0,0 +1,8 @@
1+#!/bin/bash
2+# Install CuPy and the CUDA 12 runtime libraries as pip wheels, for machines
3+# with an NVIDIA driver but no system CUDA toolkit. Then source scripts/gpuenv.sh.
4+set -euo pipefail
5+python3 -m pip install --user cupy-cuda12x numpy matplotlib \
6+ nvidia-cuda-runtime-cu12 nvidia-cuda-nvrtc-cu12 nvidia-cublas-cu12 nvidia-cusparse-cu12 \
7+ nvidia-cufft-cu12 nvidia-curand-cu12 nvidia-cusolver-cu12 nvidia-nvjitlink-cu12
8+echo "Now run: source scripts/gpuenv.sh"
tests/geometry_test.pyadded+103−0View file
@@ -0,0 +1,103 @@
1+"""Green's identity S_k[du/dn] - D_k[u] = u/2 at boundary nodes, brute force over all patches,
2+with three models of the patch geometry used inside the integrals:
3+ interp: xyz, r_u, r_v each interpolated from node values (fmm3dbie's model)
4+ consistent: xyz interpolated, r_u and r_v the derivatives of that interpolant
5+ exact: the analytic superellipsoid parametrization
6+The density (u and du/dn) is always the polynomial interpolant of node values."""
7+import os
8+import sys
9+
10+import numpy as np
11+
12+sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
13+sys.path.insert(0, os.path.join(os.path.dirname(os.path.abspath(__file__)), ".."))
14+
15+from fmm3dbie_gpu import surface as S # noqa: E402
16+from problem import green, green_dn # noqa: E402
17+from test_quadrature import kernels6 # noqa: E402
18+
19+a_, p_, m_, norder = 1.0, 8, 4, 6
20+surf = S.rounded_cube(a=a_, p=p_, m=m_, norder=norder)
21+br = S.graded_breaks(m_, 0.5)
22+npols = surf.npols
23+k = 10.0; xe = np.array([1.7, 0.6, -0.5])
24+u = green(surf.srcvals[0:3], xe, k); un = green_dn(surf.srcvals[0:3], surf.srcvals[9:12], xe, k)
25+
26+def exact_geom(patch, uv):
27+ f, rem = divmod(patch, m_ * m_); i, j = divmod(rem, m_)
28+ n, e1, e2 = (np.array(v, float) for v in S._FACES[f])
29+ s0, s1, t0, t1 = br[i], br[i + 1], br[j], br[j + 1]
30+ s = s0 + 0.5 * (uv[0] + 1) * (s1 - s0); t = t0 + 0.5 * (uv[1] + 1) * (t1 - t0)
31+ c = n[:, None] + s * e1[:, None] + t * e2[:, None]
32+ N = np.sum(c ** p_, axis=0) ** (1.0 / p_); gradN = c ** (p_ - 1) / N ** (p_ - 1)
33+ def deriv(dc, scale):
34+ dN = np.sum(gradN * dc[:, None], axis=0)
35+ return a_ * scale * (dc[:, None] / N - c * dN / N ** 2)
36+ return np.vstack([a_ * c / N, deriv(e1, 0.5 * (s1 - s0)), deriv(e2, 0.5 * (t1 - t0))])
37+
38+def dlegendre_vander(x, n):
39+ P = S.legendre_vander(x, n); D = np.zeros_like(P)
40+ for k_ in range(1, n + 1): # P_k' = sum over j = k-1, k-3, ... of (2j+1) P_j
41+ for j in range(k_ - 1, -1, -2):
42+ D[k_] += (2 * j + 1) * P[j]
43+ return D
44+
45+def consistent_geom(patch, uv):
46+ c = surf.srccoefs[0:3, patch * npols:(patch + 1) * npols]
47+ Pu, Pv = S.legendre_vander(uv[0], norder), S.legendre_vander(uv[1], norder)
48+ Du, Dv = dlegendre_vander(uv[0], norder), dlegendre_vander(uv[1], norder)
49+ B = (Pv[:, None] * Pu[None]).reshape(npols, -1)
50+ Bu = (Pv[:, None] * Du[None]).reshape(npols, -1)
51+ Bv = (Dv[:, None] * Pu[None]).reshape(npols, -1)
52+ return np.vstack([c @ B, c @ Bu, c @ Bv])
53+
54+models = {"interp": lambda P, uv: surf.eval_geometry(P, uv), "consistent": consistent_geom, "exact": exact_geom}
55+
56+def moments(geom, patch, uv, w, xt, nt):
57+ g = geom(patch, uv)
58+ ns = np.cross(g[3:6].T, g[6:9].T).T; jac = np.linalg.norm(ns, axis=0); ns = ns / jac
59+ K = kernels6(xt, nt, g[0:3], ns, k)
60+ return (K * (jac * w)[None]) @ S.basis_2d(uv, norder).T
61+
62+def self_rule(u0, v0, nx=200, nr=40):
63+ xr, wr = np.polynomial.legendre.leggauss(nr); xr, wr = 0.5 * (xr + 1), 0.5 * wr
64+ xx, wx = np.polynomial.legendre.leggauss(nx)
65+ U, W = [], []
66+ for e in range(4):
67+ h = [1 - u0, -1 - u0, 1 - v0, -1 - v0][e]
68+ x1, x2 = (-1 - v0, 1 - v0) if e < 2 else (-1 - u0, 1 - u0)
69+ for xa, xb in [(x1, 0.0), (0.0, x2)]:
70+ x = xa + 0.5 * (xb - xa) * (xx + 1); dx = 0.5 * (xb - xa) * wx
71+ R, X = np.meshgrid(xr, x, indexing="ij"); WW = np.outer(wr, dx) * abs(h) * R
72+ uv = np.vstack([(u0 + R * h).ravel(), (v0 + R * X).ravel()]) if e < 2 else \
73+ np.vstack([(u0 + R * X).ravel(), (v0 + R * h).ravel()])
74+ U.append(uv); W.append(WW.ravel())
75+ return np.hstack(U), np.hstack(W)
76+
77+def leaf_rule(depth, nq=8):
78+ gx, gw = np.polynomial.legendre.leggauss(nq); nl = 2 ** depth; h = 1.0 / nl
79+ c = -1 + h + 2 * h * np.arange(nl)
80+ Uc, Vc, A, B = np.meshgrid(c, c, gx, gx, indexing="ij")
81+ uv = np.vstack([(Uc + h * A).ravel(), (Vc + h * B).ravel()])
82+ w = (h * h * np.outer(gw, gw)[None, None]).repeat(nl, 0).repeat(nl, 1).ravel()
83+ return uv, w
84+
85+leaf6, leaf3 = leaf_rule(6), leaf_rule(3)
86+P = 13
87+for loc in [0, 3, 21, 24]:
88+ i = P * npols + loc
89+ xt, nt = surf.srcvals[0:3, i], surf.srcvals[9:12, i]
90+ out = []
91+ for name, geom in models.items():
92+ Sv = Dv = 0
93+ for j in range(surf.npatches):
94+ if j == P:
95+ uv, w = self_rule(*surf.uv_nodes[:, i])
96+ else:
97+ d = np.linalg.norm(xt - surf.centroids[:, j]) / surf.rads[j]
98+ uv, w = leaf6 if d < 3 else leaf3
99+ Wt = moments(geom, j, uv, w, xt, nt) @ surf.umat
100+ sl = slice(j * npols, (j + 1) * npols)
101+ Sv += Wt[4] @ un[sl]; Dv += Wt[5] @ u[sl]
102+ out.append(f"{name} {abs(Sv - Dv - u[i] / 2) / abs(u[i]):.1e}")
103+ print(f"node position {loc} (u,v)=({surf.uv_nodes[0, i]:+.2f},{surf.uv_nodes[1, i]:+.2f}): relative residual " + " ".join(out), flush=True)
tests/test_quadrature.pyadded+167−0View file
@@ -0,0 +1,167 @@
1+"""Check the GPU near-field weights against brute-force integration in NumPy.
2+
3+For a few (target, patch) pairs we compute the reference moments
4+I_n = int K_m(x_t, y) P_n |J| du dv with independent quadratures: for the self
5+patch, polar rules on the triangles fanned from the target to the edges, with
6+a sinh substitution along each edge and very many nodes; for neighbouring patches, a uniform quadtree of depth 6 with an
7+8x8 Gauss-Legendre rule per leaf. The reference weights are U^T I,
8+computed with the same geometry model (exact or interpolated) as the GPU code.
9+"""
10+
11+import os
12+import sys
13+import time
14+
15+import cupy as cp
16+import numpy as np
17+
18+sys.path.insert(0, os.path.join(os.path.dirname(os.path.abspath(__file__)), ".."))
19+
20+from fmm3dbie_gpu import cuda_ops as ops # noqa: E402
21+from fmm3dbie_gpu import surface as S # noqa: E402
22+
23+OVER4PI = 1 / (4 * np.pi)
24+
25+
26+def kernels6(xt, nt, y, ns, zk):
27+ """Six kernels (6, n) for target xt (3,), nt (3,), sources y (3, n), ns (3, n)."""
28+ d = xt[:, None] - y
29+ r2 = np.sum(d * d, axis=0)
30+ r = np.sqrt(r2)
31+ kabs = abs(zk)
32+ dnt = np.sum(d * nt[:, None], axis=0)
33+ dns = np.sum(d * ns, axis=0)
34+ nn = np.sum(ns * nt[:, None], axis=0)
35+ Ek = np.exp(1j * zk * r) * OVER4PI / r
36+ Ec = np.exp(-kabs * r) * OVER4PI / r
37+ fk, fc = 1 - 1j * zk * r, 1 + kabs * r
38+ gk = 3 - 3j * zk * r - zk**2 * r2
39+ gc = 3 + 3 * kabs * r + kabs**2 * r2
40+ K = np.zeros((6, y.shape[1]), dtype=complex)
41+ K[4] = Ek
42+ K[1] = Ec
43+ K[0] = -dnt * fk * Ek / r2
44+ K[2] = -dnt * fc * Ec / r2
45+ K[5] = dns * fk * Ek / r2
46+ K[3] = (nn * fk / r2 - dns * dnt * gk / r2**2) * Ek - (nn * fc / r2 - dns * dnt * gc / r2**2) * Ec
47+ return K
48+
49+
50+def moments_at(surf, patch, uv, w, xt, nt, zk):
51+ """sum_nodes K_m P_n |J| w -> (6, npols)."""
52+ g = surf.eval_geometry(patch, uv)
53+ ns = np.cross(g[3:6].T, g[6:9].T).T
54+ jac = np.linalg.norm(ns, axis=0)
55+ ns = ns / jac
56+ K = kernels6(xt, nt, g[0:3], ns, zk)
57+ P = S.basis_2d(uv, surf.norder)
58+ return (K * (jac * w)[None, :]) @ P.T
59+
60+
61+def ref_self(surf, patch, u0, v0, xt, nt, zk, nx=400, nr=80):
62+ """Polar rules on the triangles (target, edge point) with many nodes. Along each
63+ edge we split at the foot of the perpendicular and use x = |h| sinh(xi), which
64+ resolves the peak of width |h| next to the foot."""
65+ xr, wr = np.polynomial.legendre.leggauss(nr)
66+ xr, wr = 0.5 * (xr + 1), 0.5 * wr
67+ xx, wx = np.polynomial.legendre.leggauss(nx)
68+ I = np.zeros((6, surf.npols), dtype=complex)
69+ for e in range(4):
70+ h = [1 - u0, -1 - u0, 1 - v0, -1 - v0][e]
71+ if abs(h) < 1e-14:
72+ continue
73+ if e < 2:
74+ x1, x2 = -1 - v0, 1 - v0
75+ else:
76+ x1, x2 = -1 - u0, 1 - u0
77+ for xend in (x1, x2):
78+ if abs(xend) < 1e-14:
79+ continue
80+ xi1 = np.arcsinh(abs(xend) / abs(h))
81+ xi = 0.5 * xi1 * (xx + 1)
82+ x = np.sign(xend) * abs(h) * np.sinh(xi)
83+ dx = abs(h) * np.cosh(xi) * 0.5 * xi1 * wx
84+ R, X = np.meshgrid(xr, x, indexing="ij")
85+ WW = np.outer(wr, dx) * abs(h) * R
86+ if e < 2:
87+ uv = np.vstack([(u0 + R * h).ravel(), (v0 + R * X).ravel()])
88+ else:
89+ uv = np.vstack([(u0 + R * X).ravel(), (v0 + R * h).ravel()])
90+ I += moments_at(surf, patch, uv, WW.ravel(), xt, nt, zk)
91+ return I
92+
93+
94+def ref_leaf(surf, patch, xt, nt, zk, depth=6, nq=8):
95+ gx, gw = np.polynomial.legendre.leggauss(nq)
96+ nl = 2**depth
97+ h = 1.0 / nl
98+ c = -1 + h + 2 * h * np.arange(nl)
99+ U, V, A, B = np.meshgrid(c, c, gx, gx, indexing="ij")
100+ uv = np.vstack([(U + h * A).ravel(), (V + h * B).ravel()])
101+ w = (h * h * np.outer(gw, gw)[None, None]).repeat(nl, 0).repeat(nl, 1).ravel()
102+ return moments_at(surf, patch, uv, w, xt, nt, zk)
103+
104+
105+def main():
106+ """Usage: python3 test_quadrature.py [norder] [k] [key=value ...]
107+
108+ Keys: nq, nqr, nqx, sub_len, osc_len (quadrature), m (patches per face edge),
109+ tp (target patch), rfac, rfac_q, kR_max, exact (1 for the exact geometry, 0 for
110+ the interpolated one)."""
111+ norder = int(sys.argv[1]) if len(sys.argv) > 1 else 6
112+ zk = float(sys.argv[2]) if len(sys.argv) > 2 else 10.0
113+ quad = dict(nq=8, nqr=14, nqx=8, sub_len=1.0, osc_len=2.0, m=3, tp=4, rfac=2.0, rfac_q=2.0,
114+ kR_max=3.0, exact=1)
115+ for kv in sys.argv[3:]:
116+ key, val = kv.split("=")
117+ quad[key] = type(quad[key])(val)
118+ m, tp, rfac, rfac_q, kR_max, exact = (quad.pop(x) for x in ["m", "tp", "rfac", "rfac_q", "kR_max", "exact"])
119+ print("k", zk, "norder", norder, quad, "m", m, "target patch", tp, "exact geometry", bool(exact))
120+ surf = S.rounded_cube(a=1.0, p=8, m=m, norder=norder, exact=bool(exact))
121+ npols = surf.npols
122+ # targets: three nodes of patch tp (corner, centre, next to an edge); all their near pairs
123+ node_ids = [tp * npols, tp * npols + npols // 2, tp * npols + npols - 2]
124+ xt_all = surf.srcvals[0:3, node_ids]
125+ pairs = []
126+ for q, i in enumerate(node_ids):
127+ d = np.linalg.norm(surf.centroids - xt_all[:, q:q + 1], axis=0)
128+ pairs += [(i, j) for j in np.nonzero(d <= rfac * surf.rads)[0]]
129+ pair_t = cp.asarray([i for i, _ in pairs], dtype=cp.int32)
130+ pair_p = cp.asarray([j for _, j in pairs], dtype=cp.int32)
131+
132+ geo = ops.Geometry(surf)
133+ x = cp.asarray(np.ascontiguousarray(surf.srcvals[0:3].T))
134+ n = cp.asarray(np.ascontiguousarray(surf.srcvals[9:12].T))
135+ pid = cp.asarray(surf.patch_id.astype(np.int32))
136+ uv = cp.asarray(np.ascontiguousarray(surf.uv_nodes.T))
137+ umat = cp.asarray(np.ascontiguousarray(surf.umat))
138+
139+ t0 = time.time()
140+ leaf_ptr, leaves = ops.build_leaves(pair_t, pair_p, x, pid, geo, zk=zk, rfac_q=rfac_q,
141+ nlmax=14, kR_max=kR_max)
142+ W = ops.near_weights(pair_t, pair_p, x, n, pid, uv, geo, leaf_ptr, leaves, umat, zk, **quad)
143+ cp.cuda.Device().synchronize()
144+ print(f"GPU weights in {time.time() - t0:.2f} s for {len(pairs)} pairs")
145+ W = W.get()
146+ leaf_ptr = leaf_ptr.get()
147+ worst = np.zeros(6)
148+ for q, (i, j) in enumerate(pairs):
149+ xt, nt = surf.srcvals[0:3, i], surf.srcvals[9:12, i]
150+ if j == surf.patch_id[i]:
151+ I = ref_self(surf, j, *surf.uv_nodes[:, i], xt, nt, zk)
152+ kind = "self"
153+ else:
154+ I = ref_leaf(surf, j, xt, nt, zk)
155+ kind = "near"
156+ Wref = I @ surf.umat # w_j = sum_n U[n, j] I_n
157+ err = np.linalg.norm(W[q] - Wref, axis=1) / np.linalg.norm(Wref, axis=1)
158+ d = np.linalg.norm(xt - surf.centroids[:, j]) / surf.rads[j]
159+ worst = np.maximum(worst, err)
160+ if err.max() > 1e-6 or kind == "self":
161+ print(f"target {i} patch {j} ({kind}, dist/rad {d:.2f}, {leaf_ptr[q + 1] - leaf_ptr[q]} leaves): "
162+ "rel err per kernel " + " ".join(f"{e:.1e}" for e in err))
163+ print("worst per kernel:", " ".join(f"{e:.1e}" for e in worst))
164+
165+
166+if __name__ == "__main__":
167+ main()