9b737b8GPU implementation of fmm3dbie's interior Helmholtz impedance solver, with a showcase room problemJeremy Magland 1"""GPU solver for the Helmholtz impedance problem with the right-preconditioned
2combined field representation of fmm3dbie (helm_rpcomb_imp):
4 u = S_k[sigma] + i alpha D_k S_kappa[sigma], kappa = i|k|,
6with the boundary condition du/dn + i k lambda u = f. The discretized equation
7(fmm3dbie, helm_rpcomb_imp.f90) is
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,
12 phi1 = S_kappa sigma, phi2 = S'_kappa sigma,
13 z = (-1)^ifinout / 2 - i alpha / 4, w = 1/2 (exterior), -1/2 (interior).
15The hypersingular operator D'_k is never applied directly: the Calderon
16identity D'_kappa S_kappa = -1/4 + S'_kappa^2 turns it into the weakly singular
17difference D'_k - D'_kappa plus a product of two S'_kappa applications.
19Every layer potential is the sum of a far part (smooth rule on the oversampled
20nodes of the patches that are not near the target, cuda_ops.dense_apply) and a
21near part (precomputed weights for the near patches, cuda_ops.sparse_apply).
22"""
24import time
26import cupy as cp
27import numpy as np
29from . import cuda_ops as ops
30from . import surface
33class 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)
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)
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))
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")
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])
89 def oversample(self, q):
90 qp = q.reshape(self.surf.npatches, self.npols)
91 return (qp @ self.over_mat.T).ravel()
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)
97 def zid(self):
98 return (-1) ** self.ifinout / 2 - 1j * self.alpha / 4
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
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]
125 def s_kappa(self, sigma):
126 """phi1 = S_kappa[sigma] on the surface."""
127 return self.apply_single(1, sigma)
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
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).
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
165def 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