9b737b8GPU implementation of fmm3dbie's interior Helmholtz impedance solver, with a showcase room problemJeremy Magland 1"""Minimal ctypes interface to the fmm3dbie shared library.
3Only the routines needed for the interior Helmholtz impedance problem are
4wrapped. All Fortran arguments are passed by reference; integers are
5integer*8, and arrays are Fortran (column-major) ordered.
6"""
8import contextlib
9import ctypes
10import os
11import sys
13# OpenBLAS threads compete with the library's OpenMP threads
14os.environ.setdefault("OPENBLAS_NUM_THREADS", "1")
15# write Fortran console output immediately so it can be redirected to a log
16os.environ.setdefault("GFORTRAN_UNBUFFERED_PRECONNECTED", "y")
18import numpy as np
20LIBPATH = os.environ.get(
21 "FMM3DBIE_LIB", os.path.expanduser("~/src/fmm3dbie/lib/libfmm3dbie.so")
22)
23_lib = ctypes.CDLL(LIBPATH)
25_i8 = ctypes.c_int64
26_f8 = ctypes.c_double
30def 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)
44_libc = ctypes.CDLL(None)
47def _ref_i8(v):
48 return ctypes.byref(_i8(int(v)))
51def _ref_f8(v):
52 return ctypes.byref(_f8(float(v)))
55def _ptr(a):
56 return a.ctypes.data_as(ctypes.c_void_p)
59def _farr(a, dtype):
60 return np.asfortranarray(a, dtype=dtype)
63class Surface:
64 """Discretized surface in the fmm3dbie format (all patches share one order and type)."""
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)
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))
85def npols_for(norder, iptype):
86 if iptype == 1:
87 return (norder + 1) * (norder + 2) // 2
88 return (norder + 1) ** 2
91def 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
98def 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
108def 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
115def 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].
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
137def 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
bb9c969Add the cymbal problem: a listener near a rigid cymbal in free spaceJeremy Magland 152def helm_rpcomb_eval_boundary(surf, zk, alpha, sigma, siksigma, eps=1e-6, ifinout=0):
153 """Boundary limit of u = S_k[s] + i alpha D_k S_{i|k|}[s] at the nodes, from the
154 inside (ifinout = 0) or the outside (ifinout = 1).
9b737b8GPU implementation of fmm3dbie's interior Helmholtz impedance solver, with a showcase room problemJeremy Magland 155
156 The on-surface evaluation gives the principal value of D_k; with outward
bb9c969Add the cymbal problem: a listener near a rigid cymbal in free spaceJeremy Magland 157 normals the interior limit adds -1/2 times the double-layer density and the
158 exterior limit +1/2 times it.
9b737b8GPU implementation of fmm3dbie's interior Helmholtz impedance solver, with a showcase room problemJeremy Magland 159 """
160 uvs = get_disc_nodes(surf.norder, surf.iptype)
161 ipatch_id = np.repeat(np.arange(1, surf.npatches + 1), surf.npols).astype(np.int64)
162 uvs_targ = np.asfortranarray(np.tile(uvs, (1, surf.npatches)))
163 zpars = np.array([zk, alpha], dtype=np.complex128)
164 sigma = _farr(sigma, np.complex128)
165 siksigma = _farr(siksigma, np.complex128)
166 pot = np.zeros(surf.npts, dtype=np.complex128)
167 _lib.helm_rpcomb_eval_(*surf.args(), _ref_i8(12), _ref_i8(surf.npts), _ptr(surf.srcvals),
168 _ptr(ipatch_id), _ptr(uvs_targ), _ref_f8(eps), _ptr(zpars),
169 _ptr(sigma), _ptr(siksigma), _ptr(pot))
bb9c969Add the cymbal problem: a listener near a rigid cymbal in free spaceJeremy Magland 170 return pot + (-1) ** (ifinout + 1) * 0.5j * alpha * siksigma