1"""Check the GPU near-field weights against brute-force integration in NumPy.
3For a few (target, patch) pairs we compute the reference moments
4I_n = int K_m(x_t, y) P_n |J| du dv with independent quadratures: for the self
5patch, polar rules on the triangles fanned from the target to the edges, with
6a sinh substitution along each edge and very many nodes; for neighbouring patches, a uniform quadtree of depth 6 with an
78x8 Gauss-Legendre rule per leaf. The reference weights are U^T I,
8computed with the same geometry model (exact or interpolated) as the GPU code.
9"""
11import os
12import sys
13import time
15import cupy as cp
16import numpy as np
18sys.path.insert(0, os.path.join(os.path.dirname(os.path.abspath(__file__)), ".."))
20from fmm3dbie_gpu import cuda_ops as ops # noqa: E402
21from fmm3dbie_gpu import surface as S # noqa: E402
23OVER4PI = 1 / (4 * np.pi)
26def 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
50def 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
61def 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
94def 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)
105def main():
106 """Usage: python3 test_quadrature.py [norder] [k] [key=value ...]
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)
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))
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))
166if __name__ == "__main__":
167 main()