Solving acoustic scattering problems by the isogeometric boundary element method

We solve acoustic scattering problems by means of the isogeometric boundary integral equation method. In order to avoid spurious modes, we apply the combined field integral equations for either sound-hard scatterers or sound-soft scatterers. These integral equations are discretized by Galerkin's method, which especially enables the mathematically correct regularization of the hypersingular integral operator. In order to circumvent densely populated system matrices, we employ the isogeometric fast multipole method. The result is an algorithm that scales essentially linear in the number of boundary elements. Numerical experiments are performed which show the feasibility and the performance of the approach.


Introduction
Acoustic wave scattering appears in many places in engineering practice.This includes, for instance, the modelling of sonar and other methods of acoustic location, as well as outdoor noise propagation and control, especially stemming from automobiles, railways or aircrafts.Since an analytical solution of scattering problems is in general impossible, numerical approaches are called for the approximate solution.
Most acoustic scattering problems may be formulated in the frequency domain by employing the Helmholtz equation.Assume that an acoustic wave encounters an impenetrable, bounded obstacle Ω ⊂ R 3 , having a Lipschitz smooth boundary Γ := ∂Ω, and, as a consequence, gets scattered.Given the incident plane wave u inc (x) = e iκ d,x with known wavenumber κ and direction d, where d 2 = 1, the goal is to compute the scattered wave u s .The physical model behind this is as follows.The total wave u = u inc + u s satisfies the exterior Helmholtz equation The boundary condition at the scatterer's surface depends on its physical properties.If the scatterer constitutes a sound-soft obstacle, then the acoustic pressure vanishes at Γ and we have the homogeneous Dirichlet condition Whereas, if the scatterer constitutes a sound-hard obstacle, then the normal velocity vanishes at Γ and we have the homogeneous Neumann condition The behaviour towards infinity is imposed by the Sommerfeld radiation condition lim r→∞ r ∂u s ∂r − iκu s = 0, where r := x 2 . (4) It implies the asymptotic expansion as x 2 → ∞.Herein, the function u ∞ : S 1 := {x ∈ R d : x 2 = 1} → C is called the far-field pattern, which is always analytic in accordance with [6,Chapter 6].In applications, the far-field pattern is the most important quantity of interest derived from scattering problems.
To avoid the discretization of the unbounded exterior domain R 3 \Ω, one can exploit the integral equation formalism to compute the numerical solution of acoustic scattering problems.Then, one arrives at a boundary integral equation only defined on the boundary Γ .We will employ here the methodology of isogeometric analysis (IGA) to discretize this boundary integral equation.IGA has been introduced in [20] in order to incorporate simulation techniques into the design workflow of industrial development.The goal is thus to unify the CAD representation of the scatterer with the boundary element discretization of the integral equation in terms of non-uniform rational B-splines (NURBS).We refer the reader to [11,23,27] and the references therein for details of the isogeometric boundary element method.
While a reformulation of the scattering problem by means of a boundary integral equation replaces the problem posed in the unbounded domain by a problem posed on the scatterer's closed boundary, the linear operator under consideration becomes a nonlocal boundary integral operator.This results in densely populated matrices with the consequence that desirable realistic simulations would still beyond current computing capacities.Therefore, we combine the isogeometric boundary element method with the fast multipole method proposed in [17].Our particular implementation relies on interpolation as proposed in [18] and yields a black-box version of the fast multipole method which is applicable to any asymptotically smooth integral kernel.
The isogeometric boundary element method we use has been developed in [11,12] and was made accessible to the public by the software C++ library Bembel [8,9].Bembel has for example been applied successfully to engineering problems arising from electromagnetics [13,21] or from acoustics [14].It has also been used in other applications, for example, to optimize periodic structures [19], in uncertainty quantification [10,14], the coupling of FEM and BEM [15], or the partial element equivalent circuit (PEEC) method [30].In this article, we shall present results for the frequency stable solution of acoustic obstacle scattering problems by using combined field integral equations.Although we restrict ourselves to the sound-soft and sound-hard cases, the presented concepts are also suitable to treat penetrable obstacles, i.e. objects described by a different diffractive index to the free space.
The rest of this article is structured as follows.In Section 2, we introduce the frequency stable boundary integral equations which are employed to solve either sound-hard or sound-soft scattering problems.Section 3 recapitulates the basic concepts from isogemeometric analysis and introduces the discretization spaces that will be used later on.In Section 4, we discuss the discretization of the required boundary integral operators.
In particular, we address the regularization of the hypersingular operator.Moreover, we comment on the isogeometric fast multipole method for the fast assembly of the operators and the potential evaluation.The numerical experiments are presented in Section 5, while concluding remarks are stated in Section 6.

Boundary integral equation method
In order to solve the boundary value problem (1)-(4), we shall employ a suitable reformulation by boundary integral equations.To this end, we introduce the acoustic single layer operator the acoustic double layer operator ∂n y ρ(y) dσ y , its adjoint as well as the acoustic hypersingular operator Here, n x and n y denote the outward pointing normal vectors at the surface points x, y ∈ Γ , respectively, while G(•, •) denotes the fundamental solution for the Helmholtz equation.In three spatial dimensions, the latter is given by Although the Helmholtz problem (1)-( 4) is uniquely solvable, a respective boundary integral formulation might not if κ 2 is an eigenvalue for the Laplacian inside the scatterer Ω.In order to avoid such spurious modes, we employ combined field integral equations in the following.Then, for some real η = 0, the solution of the boundary integral equation gives rise to the scattered wave in accordance with in case of sound-soft scattering problems.In case of sound-hard obstacles, we will solve the integral equation Having solved (8), the scattered wave is computed by Notice that the boundary integral equations ( 6) and ( 8) are always uniquely solvable, independent of the wavenumber κ, compare [5,6,22].

B-splines
We shall give a brief introduction to the basic concepts of isogeometric analysis, starting with the definition of the B-spline basis, followed by the description of the scatterer by using NURBS.To this end, let K be either R or C. The original definitions (or equivalent notions) and proofs, as well as basic algorithms, can be found in most of the standard spline and isogeometric literature [7,20,25,26,29].
We define a p-open knot vector as a set where k denotes the number of control points.The associated basis functions are given by {b p j } k−1 j=0 for p = 0 as and for p > 0 via the recursive relationship cf. Figure 1.A B-spline is then defined as a function where Having the spline functions at hand, we can introduce the spline spaces which serve as fundament for the definition of the ansatz and test spaces of the boundary element method.
Definition 2 Let Ξ be a p-open knot vector containing k + p + 1 elements.We define the spline space S p (Ξ) as the space spanned by {b p j } k−1 j=0 .Finally, we should consider the relation between the spline spaces and the underlying mesh relative to a certain mesh size.
Definition 3 For a knot vector Ξ, we define the mesh size h to be the maximal distance between neighbouring knots.We call a knot vector quasi uniform, when there exists a constant θ ≥ 1 such that for all j the ratio h j • h −1 j+1 satisfies θ −1 ≤ h j • h −1 j+1 ≤ θ.B-splines on higher dimensional domains are constructed through simple tensor product relationships for p j1,..
Throughout this article, we will reserve the letter h for the mesh size (10).All knot vectors will be assumed to be quasi uniform, such that the usual spline theory is applicable [1,25,26].

Isogeometric representation of the scatterer
We assume that the boundary Γ of the scatterer is closed and Lipschitz continuous.For the remainder of this article, we assume that it is given patchwise as Γ = n j=1 Γ j , i.e. that it is induced by (smooth, nonsingular, bijective) diffeomorphisms In the spirit of isogemetric analysis, these mappings are given by NURBS mappings, i.e. by with control points c j1,j2 ∈ R 3 and weights w i1,i2 > 0. We will moreover require that, for any interface D = Γ j ∩Γ i = ∅, the NURBS mappings coincide, i.e. that, up to rotation of the reference domain, one finds F j (•, 1) ≡ F i (•, 0).

Ansatz and test spaces
The mappings of (12) give rise to the transformations which can be utilized to define discrete spaces patchwise, by mapping the space of tensor product B-splines as in (11) with Ξ p,m := 0, . . ., 0 p+1 times , 1/2 m , . . ., (2 m − 1)/2 m , 1, . . ., 1 p+1 times to the geometry.Here, the variable m denotes the level of uniform refinement.For the purposes of discretizing V, K, and K , the global function space on Γ defined by for some g ∈ S p,p (Ξ p,m , Ξ p,m ) , as commonly done in the isogeometric literature, see e.g.[3,4], is sufficient.Note that the spline space S 2 p,m (Γ ) is of dimension n•(2 m +p) 2 , where n denotes the number of patches involved in the description of the geometry.For the purposes of discretizing W, we also require the space for some g ∈ S p,p (Ξ p,m , Ξ p,m ) , see, e.g., also [3,4].Note that that S 0 p,m (Γ ) ⊂ S 2 p,m (Γ ) consists of globally continuous B-splines whereas S 2 p,m (Γ ) is discontinuous across patch boundaries.

Galerkin method
With the boundary integral equations and a collection of spline spaces available, we are now in the position to discretize ( 6) and (8).We consider a Galerkin discretization in the L 2 (Γ )-duality product with the spline spaces S 2 p,m (Γ ) and S 0 p,m (Γ ) as ansatz and test spaces.Thus, the discrete variational formulation for (6) reads with the Galerkin approximation t h ≈ ∂u/∂n.Choosing a basis S 2 p,m (Γ ) = span{ψ 2,1 , . . ., ψ 2,N } leads to the system of linear equations with , and t being the coefficient vector of t h .The discrete variational formulation for (8) reads with the Galerkin approximation g h ≈ u| Γ .Choosing a basis S 0 p,m (Γ ) = span{ψ 0,1 , . . ., ψ 0,M } leads to the linear system of equations with , and g being the coefficient vector of g h .It is well known that the matrices V 2 , K 0 , K 2 , and W 0 are dense, which makes the assembly and storage of these matrices as well as the solution of the corresponding linear systems of equations computationally prohibitively expensive for higher resolution of the ansatz spaces, i.e., large M or N .This is why we shall apply the multipole method presented in Subsection 4.4.

Reformulation on the reference domain
Due to the isogeometric representations of the geometry, the bilinear forms for the computation of the matrix entries can entirely pulled back to the reference domain [18].To this end, let A with be one of the operators V, K, or K and µ, ν : Γ → C be functions of sufficient regularity.Defining the surface measure of a mapping F j for x = (x, y) ∈ [0, 1] 2 as the bilinear forms for the matrix entries can be recast as with the pull-back of the kernel function and the ansatz and test functions Applying a similar reasoning to the right-hand side yields Due to the additional derivative, the hypersingular operator W requires a special treatment which we will elaborate next.

Regularization of the Helmholtz hypersingular operator
The hypersingular operator W from ( 5) does not have a well defined integral operator representation as in (16).Instead, it is common knowledge that the operator can be replaced by a regularized one in case of a Galerkin discretization.Namely, for the computation of the matrix entries, the representation i, j = 1, . . ., M , can be used, cf.e.g.[24].Therein, curl Γ ψ 0,i denotes the surface curl which maps a scalar valued function on the surface into a vector field in the tangential space of Γ .On any given patch Γ j , the isogeometric representations of the boundary of the scatterer allow for its explicit representation for all x = F j (x) ∈ Γ j , x ∈ [0, 1] 2 , see [12] for example for the precise derivation.This amounts to the following expression of the hypersingular operator in closed form where the pull-back of the kernel k j,i is given by and K j,i denotes the first fundamental tensor of differential geometry, Compared to the Laplace case, see [12], we note the occurrence of a second term in the regularized representation (19).However, this additional term behaves similar to the single layer operator and thus poses no further challenges for implementation.
For the numerical evaluation of the first term in (19), recall that an ansatz function ψ 0,j | Γi on the patch Γ i is given by ψ 0,j = ι −1 i ( ψ) for some ψ ∈ S p,p (Ξ p,m , Ξ p,m ), see (13).There therefore holds Thus, for the purposes of implementation, one only has to provide S p,p (Ξ p,m , Ξ p,m ) and the directional derivatives of S p,p (Ξ p,m , Ξ p,m ), which are given by where Ξ p,m denotes the truncation of Ξ p,m , i.e., the knot vector Ξ p,m without its first and last knot.These spline spaces are readily available.

Fast multipole method
The black-box fast multipole method, cp.[16], relies on a degenerate kernel approximation of the integral kernel under consideration.Such an approximation is available in the kernel's far-field, which means that the supports of the trial and test functions have to be sufficiently distant from each other -they are admissible.
One arrives at an efficient algorithm, if one subdivides the set of trial functions hierarchically into socalled clusters.Then, the kernel interaction of two clusters is computed by using the degenerate kernel approximation if the clusters are admissible.This means a huge matrix block in the system matrix is replaced by a low-rank matrix.If the clusters are not admissible, then one subdivides them and considers the interactions of the respective children.That way, the assembly of the Galerkin matrix can be perfomed in essentially linear complexity.
For the realization of the multipole method in the present context of isogeomtric boundary element methods, we refer the reader to [11,12].A particular advantage of the referred compression method is that the isogeometric setting allows to perform the compression of the system matrix in the reference domain rather than the computational domain.This means that we consider the pull-back of the kernel (17) instead of the kernel in free space, as originally proposed in [18].Thus, the rank of the low-rank blocks in the number of onedimensional interpolation points p decreases from O(p 3 ) Fig. 2: Torus represented by 16 patches and illustration of its dimensions.
to O(p 2 ).The compressed matrix is finally represented in the H 2 -matrix format as usual, see [2].
For the potential evaluation, i.e., for evaluating ( 7) and ( 9), we exploit a similar approximation of the kernel function.However, this time we perform the lowrank approximation in physical space, that is, we employ a degenerate kernel approximation for the kernel k.Rather than clustering elements as before, we directly cluster evaluation and quadrature points and realize the potential evaluations by means of matrix-vector multiplications.The rank of the low-rank blocks is in this case O(p 3 ).In particular, we may employ a matrix-free version, as all blocks are only required once.The advantage of this approach becomes immanent if the number of potential evaluation points increases proportionally to the number of boundary elements.In this case, the cost of the proposed potential evaluations scales essentially linearly instead of quadratically.

Setup
The numerical experiments are performed by using the publicly available C++ library Bembel, see [8,9].To this end, the previously not available operators (double layer, adjoint double layer, and hypersingular operator) were implemented.Each of the matrices in the combined field integral equations ( 14) and ( 15) was computed separately in compressed form as H 2 -matrix by using the fast multipole method on the reference domain from [9,12].The compression parameters for the fast multipole method were set to the default values (η = 1.6, nine interpolation points per direction), see [9,12] for more details.The product of the matrix sums with vectors was implemented using lazy evaluation and the arising systems of linear equations 14 and 15 were solved up to relative machine precision by means of a restarted GMRES method with a restart after 30 iterations.Finally, all computations were performed in parallel by using the built-in OpenMP-parallelization of Bembel on a compute server with 1.3 terrabyte RAM and four Intel(R) Xeon(R) E7-4850 v2 CPU with twelve 2.30GHz cores each and hyperthreading disabled.

Convergence benchmark
In order to study convergence rates, we consider a torus with major radius two and minor radius 0.5 that is represented by 16 patches, see Figure 2 for an illustration.On this geometry, we aim at computing the scattered wave of a plane incident wave in x direction with wavenumber 2.5.The scattered wave is then measured on 100 points distributed on a sphere with radius 5 around the origin.We refer to Figure 3 for an illustration of the Dirichlet data of the total wave (top plot) in case of a sound-hard torus and the Neumann data of the total wave (bottom plot) in case of a sound-soft torus.
The optimal convergence rates for the potential evaluation in case of splines of degree p are O h 2p+2 for the boundary integral equation ( 6) which corresponds to sound-soft obstacles and O h 2p+1 for the bound- ary integral equation ( 8) which corresponds to soundhard obstacles.Since the obstacle under consideration is smooth, we should achieve these convergence rates.Note that these rates are twice as high as for the collocation method and are known as the superconvergence of the Galerkin formulation, compare [28] for example.Figure 4 validates that we indeed reach these theoretical achievable convergence rates when compared to solutions obtained from an indirect formulation using a single layer or adjoint double layer ansatz, respectively.
Figure 5 illustrates the scaling of the runtimes of the computations.Instead of a quadratic scaling of the runtimes, which we would have in the case of a traditional boundary element method, one figures out that the multipole-accelerated isogeometric boundary element method scales essentially linearly as expected.This enables large-scale calculations as we will consider in the next example.

Computational benchmark
As a computational benchmark, we consider a turbine with ten blades that is parametrized by 120 patches as illustrated in Figure 6.Thereof, it can be figured out that the turbine has a diameter of 5. Again, we compute the scattered wave of a plane incident wave in x direction, but with wavenumber 1.0.
We choose cubic B-splines and three refinement levels to discretize the Cauchy data u and ∂u/∂n on the surface geometry.This results in 14'520 degrees of freedom in case of a sound-soft turbine and 12'000 degrees of freedom in case of a sound-hard turbine, respectively.The overall solution time for assembly and solution of the underlying systems of linear equations requires only about a few hours.
We compute next the scattered wave in a cylinder on up to 3'664'832 points, compare Figure 7 for an illustration.To demonstrate the efficiency of the fast potential evaluation, we compare the scaling of the multipole-Fig.6: Turbine geometry with 120 patches.accelerated potential evaluations with the traditional potential evaluations.Figure 8 illustrates that -after a certain warm-up phase for only a few potential points -the H 2 -matrix accelerated potential evaluation is indeed superior to the conventional one when increasing the number of evaluation points.Consequently, the calculation of the scattered wave also in free space becomes feasible and very efficient.

Conclusion
We have solved acoustic scattering problems for soundsoft and sound-hard scatterers by means of frequency stable combined field integral equations and the isogeometric boundary integral equation method.The major advantage of this approach is that no mesh for the unbounded exterior domain is required.Our discretization is based on Galerkin's method together with an appropriate regularization of the hypersingular operator.The method becomes computationally efficient by the use of a black-box fast multipole method tailored to isogeometric surfaces.A similar approach is employed for the postprocessing step to efficiently evaluate the solution in free space.We have presented convergence benchmarks that impressively demonstrate the high ac-curacy of the isogeometric boundary element method.
In addition, we have considered a complex computational benchmark on a complex geometry, which corroborates the feasibility of the approach in the engineering practice.Fig. 7: Scattered wave evaluated in 3'664'832 points for the sound-hard case.

Fig. 3 :
Fig. 3: The Dirichlet data u (top) of the total wave in case of a sound-hard torus and the Neumann data ∂u/∂n (bottom) of the total wave in case of a soundsoft torus.

2 p = 3 p = 4 Fig. 5 :
Fig. 5: Scaling of the combined field integral equations for various polynomial degrees.The dashed lines illustrate log-linear scaling.

Fig. 8 :
Fig. 8: Computation time of conventional and H 2 -matrix accelerated potential evaluation for various number of points.