"""GPU solver for the Helmholtz impedance problem with the right-preconditioned combined field representation of fmm3dbie (helm_rpcomb_imp): u = S_k[sigma] + i alpha D_k S_kappa[sigma], kappa = i|k|, with the boundary condition du/dn + i k lambda u = f. The discretized equation (fmm3dbie, helm_rpcomb_imp.f90) is z sigma + S'_k sigma + i alpha S'_kappa phi2 + i alpha (D'_k - D'_kappa) phi1 + i k lambda [S_k sigma + i alpha D_k phi1 + i alpha w phi1] = f, phi1 = S_kappa sigma, phi2 = S'_kappa sigma, z = (-1)^ifinout / 2 - i alpha / 4, w = 1/2 (exterior), -1/2 (interior). The hypersingular operator D'_k is never applied directly: the Calderon identity D'_kappa S_kappa = -1/4 + S'_kappa^2 turns it into the weakly singular difference D'_k - D'_kappa plus a product of two S'_kappa applications. Every layer potential is the sum of a far part (smooth rule on the oversampled nodes of the patches that are not near the target, cuda_ops.dense_apply) and a near part (precomputed weights for the near patches, cuda_ops.sparse_apply). """ import time import cupy as cp import numpy as np from . import cuda_ops as ops from . import surface class Operator: def __init__(self, surf, zk, alpha, zlams, ifinout=0, nover=None, rfac=None, rfac_q=2.0, nq=8, nqr=None, nqx=8, sub_len=1.0, osc_len=2.0, kR_max=3.0, nlmax=14, far_dtype=np.float32, verbose=True): t0 = time.time() self.surf = surf self.zk = complex(zk) self.alpha = complex(alpha) self.ifinout = ifinout norder, npols = surf.norder, surf.npols self.norder, self.npols, self.npts = norder, npols, surf.npts if rfac is None: rfac = 2.0 if norder <= 6 else 1.25 # fmm3dbie's get_rfacs self.rfac = rfac self.nover = norder + 2 if nover is None else nover if nqr is None: # the radial rule must resolve the basis and exp(ikr) along a ray nqr = int(min(max(12, norder + 2 * abs(zk) * surf.rads.max() + 6), 32)) self.quad = dict(nq=nq, nqr=nqr, nqx=nqx, sub_len=sub_len, osc_len=osc_len) self.leafpar = dict(rfac_q=rfac_q, nlmax=nlmax, kR_max=kR_max) self.zlams = cp.asarray(np.broadcast_to(np.asarray(zlams, dtype=complex), (surf.npts,))) self.geo = ops.Geometry(surf) self.umat = cp.asarray(np.ascontiguousarray(surf.umat)) self.x = cp.asarray(np.ascontiguousarray(surf.srcvals[0:3].T)) self.n = cp.asarray(np.ascontiguousarray(surf.srcvals[9:12].T)) self.patch_id = cp.asarray(surf.patch_id.astype(np.int32)) self.uv = cp.asarray(np.ascontiguousarray(surf.uv_nodes.T)) self.wts = cp.asarray(surf.wts) self.cen = cp.asarray(np.ascontiguousarray(surf.centroids.T)) self.r2 = cp.asarray((rfac * surf.rads) ** 2) ox, on, ow = surf.oversampled(self.nover) self.npols_o = (self.nover + 1) ** 2 self.far = ops.FarSources(np.ascontiguousarray(ox.T), np.ascontiguousarray(on.T), ow, self.npols_o, self.cen, self.r2, far_dtype) self.over_mat = cp.asarray(surface.oversample_matrix(norder, self.nover)) # near weights for on-surface targets (the nodes themselves) self.pair_t, self.pair_p, self.row_ptr = ops.find_near_pairs(self.x, self.cen, self.r2) self.W, self.nleaves = self._weights(self.pair_t, self.pair_p, self.x, self.n, self.patch_id, self.uv) cp.cuda.Device().synchronize() self.setup_time = time.time() - t0 if verbose: print(f" operator setup: {self.npts} nodes, {ox.shape[1]} oversampled nodes, " f"{int(self.pair_t.size)} near pairs, {self.nleaves} leaves, " f"{'exact' if self.geo.gmode == 1 else 'interpolated'} geometry, " f"{np.dtype(far_dtype).name} far field, {self.setup_time:.2f} s") def _weights(self, pair_t, pair_p, tx, tn, tpatch, tuv): leaf_ptr, leaves = ops.build_leaves(pair_t, pair_p, tx, tpatch, self.geo, zk=self.zk, **self.leafpar) W = ops.near_weights(pair_t, pair_p, tx, tn, tpatch, tuv, self.geo, leaf_ptr, leaves, self.umat, self.zk, **self.quad) return W, int(leaf_ptr[-1]) def oversample(self, q): qp = q.reshape(self.surf.npatches, self.npols) return (qp @ self.over_mat.T).ravel() def _dense(self, tx, tn, qa, qb, sel): return ops.dense_apply(tx, tn, self.far, self.oversample(qa), self.oversample(qb), self.zk, sel) def zid(self): return (-1) ** self.ifinout / 2 - 1j * self.alpha / 4 def apply_K(self, sigma): """The operator without the identity term z.""" k, a, lam = self.zk, self.alpha, self.zlams zero = cp.zeros_like(sigma) F = self._dense(self.x, self.n, sigma, zero, [1, 1, 1, 0, 1, 0]) ops.sparse_apply(self.row_ptr, self.pair_p, self.npols, self.W, sigma[None, :], [(0, 0), (1, 0), (2, 0), (4, 0)], F) phi1, phi2 = F[1], F[2] res = F[0] + 1j * k * lam * F[4] G = self._dense(self.x, self.n, phi1, phi2, [0, 0, 2, 1, 0, 1]) ops.sparse_apply(self.row_ptr, self.pair_p, self.npols, self.W, cp.stack([phi1, phi2]), [(2, 1), (3, 0), (5, 0)], G) res += 1j * a * G[2] + 1j * a * G[3] - k * a * lam * G[5] res += (-1) ** self.ifinout * k * a * lam * phi1 / 2 return res def apply_single(self, m, q): """One of the six operators applied to q on the surface (principal value).""" q = cp.asarray(q, dtype=cp.complex128) sel = [0] * 6 sel[m] = 1 F = self._dense(self.x, self.n, q, cp.zeros_like(q), sel) ops.sparse_apply(self.row_ptr, self.pair_p, self.npols, self.W, q[None, :], [(m, 0)], F) return F[m] def s_kappa(self, sigma): """phi1 = S_kappa[sigma] on the surface.""" return self.apply_single(1, sigma) def solve(self, rhs, eps_gmres=1e-6, numit=400, verbose=True): t0 = time.time() rhs = cp.asarray(rhs, dtype=cp.complex128) sigma, niter, rres = gmres(self.apply_K, rhs, self.wts, self.zid(), eps_gmres, numit) phi1 = self.s_kappa(sigma) cp.cuda.Device().synchronize() if verbose: print(f" GMRES: {niter} iterations, relative residual {rres:.2e}, {time.time() - t0:.2f} s") return sigma, phi1, niter, rres def eval(self, targ_x, sigma, phi1, on_surface=False): """u = S_k sigma + i alpha D_k phi1 at targets (ntarg, 3). With on_surface=True the targets are the surface nodes and the value is the limit from the side selected by ifinout. """ a = self.alpha if on_surface: tx, tn = self.x, self.n row_ptr, pair_p, W = self.row_ptr, self.pair_p, self.W else: tx = cp.ascontiguousarray(cp.asarray(targ_x, dtype=cp.float64)) tn = cp.zeros_like(tx) pair_t, pair_p, row_ptr = ops.find_near_pairs(tx, self.cen, self.r2) tpatch = cp.full(tx.shape[0], -1, dtype=cp.int32) tuv = cp.zeros((tx.shape[0], 2)) W, _ = self._weights(pair_t, pair_p, tx, tn, tpatch, tuv) F = self._dense(tx, tn, sigma, phi1, [0, 0, 0, 0, 1, 2]) ops.sparse_apply(row_ptr, pair_p, self.npols, W, cp.stack([sigma, phi1]), [(4, 0), (5, 1)], F) u = F[4] + 1j * a * F[5] if on_surface: u += 1j * a * (-1) ** (self.ifinout + 1) * phi1 / 2 return u def gmres(apply_K, rhs, wts, zid, eps, numit): """Unrestarted GMRES for (zid I + K) x = rhs with the weighted inner product = sum conj(u) v wts, following fmm3dbie's zgmres_guru.""" n = rhs.size wnorm = lambda v: float(cp.sqrt(cp.sum(cp.abs(v) ** 2 * wts)).real) rb = wnorm(rhs) V = cp.zeros((numit + 1, n), dtype=cp.complex128) H = np.zeros((numit + 1, numit), dtype=complex) cs, sn = np.zeros(numit, dtype=complex), np.zeros(numit, dtype=complex) svec = np.zeros(numit + 1, dtype=complex) svec[0] = rb V[0] = rhs / rb it_done = 0 for it in range(numit): w = apply_K(V[it]) # modified Gram-Schmidt in blocks (classical within one pass, repeated once) for _ in range(2): h = (cp.conj(V[:it + 1]) * (w * wts)[None, :]).sum(axis=1) w = w - h @ V[:it + 1] H[:it + 1, it] += h.get() H[it, it] += zid hn = wnorm(w) H[it + 1, it] = hn if hn > 0: V[it + 1] = w / hn for k in range(it): t1 = cs[k] * H[k, it] + sn[k] * H[k + 1, it] t2 = -np.conj(sn[k]) * H[k, it] + cs[k] * H[k + 1, it] H[k, it], H[k + 1, it] = t1, t2 a, b = H[it, it], H[it + 1, it] r = np.sqrt(abs(a) ** 2 + abs(b) ** 2) if r == 0: cs[it], sn[it] = 1, 0 else: cs[it] = abs(a) / r sn[it] = (a / abs(a)) * np.conj(b) / r if abs(a) > 0 else np.conj(b) / r H[it, it] = cs[it] * a + sn[it] * b H[it + 1, it] = 0 svec[it + 1] = -np.conj(sn[it]) * svec[it] svec[it] = cs[it] * svec[it] it_done = it + 1 if abs(svec[it + 1]) / rb < eps: break y = np.linalg.solve(np.triu(H[:it_done, :it_done]), svec[:it_done]) x = cp.asarray(y) @ V[:it_done] resid = zid * x + apply_K(x) - rhs return x, it_done, wnorm(resid) / rb