An algebraic multigrid based shifted-Laplacian preconditioner for the Helmholtz equation

A preconditioner deﬁned by an algebraic multigrid cycle for a damped Helmholtz operator is proposed for the Helmholtz equation. This approach is well-suited for acoustic scattering problems in complicated computational domains and with varying material properties. The spectral properties of the preconditioned systems and the convergence of the GMRES method are studied with linear, quadratic, and cubic ﬁnite element discretizations. Numerical experiments are performed with two-dimensional problems describing acoustic scattering in a cross section of a car cabin and in a layered medium. Asymptotically the number of iterations grows linearly with respect to the frequency while for lower frequencies the growth is milder. The proposed preconditioner is particularly effective for low-frequency and mid-frequency problems.


Introduction
Acoustic scattering problems have applications in many disciplines.These problems can be typically modeled using wave equation and often it is sufficient to consider only time-harmonic solutions which are described by the Helmholtz equation (reduced wave equation).For numerical simulation the equations can be discretized using the finite difference method or the finite element method, for example.The solution of resulting systems of linear equations can be a computationally challenging problem.
During the past few decades, numerical methods for acoustics have been under active research.Finite element method has emerged as a generic tool for discretizing the Helmholtz equation in complex geometries.A recent review [1] offers a glance at research efforts in this field.The efficiency of these methods still often limits the feasible size of scattering problems in midfrequency and high-frequency regime.Particularly the finite element phase shift (pollution) error necessitates finer discretizations for high-frequency problems [2].The finite element method have been used successfully for interior problems like scattering in a car cabin [3] as well as for exterior problems.Since the paper [4] the research on the construction of absorbing boundary conditions and absorbing layers at the truncation boundary of the exterior domain has been active; see [1] and references therein.
The resulting systems of linear equations from the discretization of the Helmholtz equation are non-Hermitian and indefinite.Furthermore, for mid-frequency and high-frequency problems, the systems can be extremely large.These reasons make them a challenge for the current solvers.Often it is feasible to use direct methods for solving these systems for twodimensional problems, but three-dimensional problems lead to systems which cannot be solved by these methods with affordable computing effort.Hence, it is necessity to use iterative methods such as the GMRES method [5] or the BICGSTAB method [6].However, these methods require a good preconditioner for the discretized Helmholtz equations in order to have reasonably fast convergence.
Various preconditioners and iterative solution techniques have been proposed for the discrete Helmholtz equation.Several domain decomposition methods have been proposed; see [7,8,9,10,11], for example.Multigrid methods have been considered in [12,13,14].With multigrid methods, it is difficult to define a stable and sufficiently accurate coarse grid problems and smoothers for them.For problems in homogenous medium, domain imbedding/fictitious domain methods in [15,16,17] have been fairly effective.An incomplete factorization preconditioner has been considered in [18], for example, and in [19] a tensor product preconditioner is used.An alternative iterative approach for solving the Helmholtz equation has been proposed in [20] and further studied in [21].The basic idea is to find a time-periodic solution to the wave equations using a controllability method, which leads to preconditioned conjugate gradient iterations for initial data.
In this paper, we consider shifted-Laplacian preconditioners which are ob-tained from the Helmholtz operator by adding damping.The Laplace operator was proposed as a preconditioner for the Helmholtz equation in [22].A shifted-Laplacian preconditioner obtained by changing the sign of the zeroth order term in the Helmholtz operator was described in [23].As a generalization the Laplacian with a complex shift was studied in [24].Following this approach a multigrid preconditioner based on a damped Helmholtz operator was considered in [14].There, the scattering problems were posed in a rectangular domain, they were discretized using low-order finite differences, and a geometric multigrid method was used.This paper extends this approach for general shaped domains using linear, quadratic, and cubic finite element discretizations.Particularly quadratic and cubic finite elements help to reduce the number of unknowns in order to reach prescribed accuracy, as they have much smaller interpolation and phase shift errors than linear basis functions [2].Our preconditioner is based on an algebraic multigrid method which can be constructed fully algebraically when the matrix for the zeroth-order terms is also available.
There is a wide range of applications for acoustic scattering in the industry and sciences.In many applications the aim is to reduce the noise level.As an example of such a problem we consider the noise in a car cabin; see also [3].Geophysical surveys employ acoustic/elastic backscattering from different layers to reconstruct a model for the subsurface.These problems lead to very large-scale scattering problems.We consider a three layer wedge model [19,14] in our numerical experiments.Acoustic scattering simulations have also many applications in medicine, sonar, and sound preproduction, for example.This paper is organized as follows.In Section 2 we describe the Helmholtz model problem and its discretization.The iterative solution and preconditioning based on shifted-Laplacian preconditioners are discussed in Section 3. The algebraic multigrid method employed in the preconditioning is described in Section 4. Then numerical results are presented in Section 5 and finally, conclusions are given in Section 6.

Scattering problem and finite element discretization
Under suitable assumptions on medium, acoustic scattering can be described by the wave equation where p (x, t) is pressure field, ρ (x) is the density of the material, c (x) is the speed of sound and t is time.For a time-harmonic pressure p (x, t) = u (x) e −iωt , where ω is angular velocity and i = √ −1, (1) leads to the Helmholtz equation where k (x) = ω/c (x) is the wave number.In inhomogeneous medium the wave number k varies depending on location as the sound speed c varies.
We consider three types of boundary conditions.In order to describe them, we decompose the boundary Γ = ∂Ω into three non-overlapping parts Γ d , Γ s , and Γ a such that Γ = Γ d ∪ Γ s ∪ Γ a .Some of these boundary sets can be empty.The first type of boundary is sound-soft which is described by the Dirichlet boundary condition where g describes the sound source, for example, an incident field.The second type is the impedance boundary condition where γ (x) is an absorbency coefficient in the range [0, 1] describing the amount of absorption on the boundary Γ a .The specific case γ = 0, leading to a Neumann boundary condition, corresponds to a sound-hard boundary without any absorption.
Exterior problems are truncated into a bounded domain Ω with Γ s as the truncation boundary.The boundary condition on Γ s should let outgoing waves propagate out of the domain without any reflection, as the Sommerfeld radiation condition describes.Such a perfect absorbing boundary condition is a non-local operator which is computationally difficult.Instead, it is usual to approximate it by a local operator [4].Here we use a first-order absorbing boundary condition For the weak formulation of the Helmholtz equation, we define the test function space and the solution space Now the weak form of (2) reads: Find u ∈ V g such that for all v ∈ V .
For a finite element discretization, we define a triangulation given by a set of non-overlapping triangles K h , such that Ω h = τ ∈K h τ .Here h denotes the diameter of the largest triangle and Ω h is an approximation of Ω.An example of a coarse triangulation (also called mesh) for a cross section of a car cabin is shown in Figure 1.
For the finite elements of order m a discrete test function space is where P m denotes polynomials of order m.A discrete solution space V g,h is obtained similarly by approximating g on Γ d,h instead of zero.In the following, we consider linear, quadratic, and cubic finite elements, that is, m = 1, 2, or 3. We use Lagrangian basis functions for the spaces V h and V g,h .Let the vector u contain the nodal values of u.Then using the discrete spaces instead of V and V g in ( 8) and integrating over the discrete counterparts of the domain and boundaries, we obtain a system of linear equations where A is a sparse matrix and f is a non-zero vector due to the inhomogenous Dirichlet boundary condition.The approximation properties of such finite element discretizations for the Helmholtz equation have been studied in [2].For an algebraic definition of the preconditioner described in Section 3, we define a mass matrix like M which includes the term k 2 /ρ, that is, it corresponds to the integral where u h ∈ V g,h and v h ∈ V h .Furthermore, we define the matrix K = A + M which contains the rest of the terms in the weak form.

Iterative solution and shifted-Laplacian preconditioner
The matrix A in ( 10) is indefinite and symmetric, but not Hermitian.Hence, the generalized minimal residual (GMRES) method [5] and the BICGSTAB method [6] are suitable iterative methods for the solution of the system (10).
For these and other applicable iterative methods, see [25], for example.At each iteration, the GMRES method minimizes the norm of the residual vector on a Krylov subspace associated to the iteration.This is a desirable property leading to a monotonic reduction of residual norm over iterations, but a disadvantage is that a basis for the Krylov subspace needs to be formed and stored.Due to this the computational cost of the GMRES methods grows quadratically with iterations and the memory requirement grows linearly.
With the BICGSTAB method the computational cost grows linearly and the memory requirement is independent of number of iterations, but convergence can be erratic and slower than with the GMRES method.In the numerical experiments we use the GMRES method.
For medium-and large-scale scattering problems, the system (10) is badly conditioned, which leads to a very slow convergence of Krylov subspace methods when applied directly to the system (10).In order to improve the conditioning and the speed of convergence, we use a right preconditioner denoted by B. This leads to a preconditioned system Once ũ is solved from this system, the solution u is obtained as u = B −1 ũ.
Our aim is to find such a preconditioner B that the matrix AB −1 is well conditioned and that vectors can be multiplied by B −1 , that is, solve systems with B, with a small computational effort.These properties would lead to a fast convergence of the iterative method and to a small overall computational cost.
In 2004, Erlangga, Vuik, and Oosterlee suggested in [24] to construct a preconditioner B SL by discretizing a shifted-Laplace operator where we have added the density ρ (x) in the operator.Using the notations defined in Section 2, the preconditioner can be defined algebraically as By choosing β 1 = 1 and β 2 to be positive, B SL is the original Helmholtz operator with some additional damping.Such damping leads to good con-ditioning of AB −1 SL and it is easier to solve systems with B SL than with A [24].In [14], Erlangga, Oosterlee, and Vuik approximated the inverse of the shifted-Laplacian preconditioner B SL using one cycle of a geometric multigrid method; see [26], for example.We denote such multigrid based preconditioners by B M G .This leads to a good conditioning of AB −1 M G for lowfrequency problems, while the number of BICGSTAB iterations appeared to grow linearly with frequency for high-frequency problems.They also showed that this preconditioner is well suited for problems with a highly varying speed of sound.In this paper, we replace the geometric multigrid method with a more generic algebraic multigrid method described in Section 4.
For the GMRES method, convergence estimates can be derived based on the spectrum of a matrix and its non-normality [5,25].Similarly to [24,14], we study numerically the spectrum of the preconditioned matrices.For small problems, it is possible to compute the spectrum, while for larger problems we can only approximate it.The GMRES method forms the basis for a Krylov subspace using the Arnoldi iteration.After m iterations it has generated an m × m upper Hessenberg matrix which is usually denoted by H m .The eigenvalues of H m approximate the eigenvalues of the system matrix.

Algebraic multigrid method
The preconditioner B M G is based on an algebraic multigrid (AMG) method which approximates the multiplication by the inverse of B SL .We use an AMG method introduced by Kickinger in [27] with modifications proposed in [28].This method uses a graph based on the system matrix to construct coarse spaces.Furthermore, it eliminates the degrees of freedom associated to the Dirichlet boundaries after forming the matrices for the coarse spaces.Under these choices, the AMG method can be constructed in such a way that the coarse problems coincide with the ones obtained using a geometric multigrid method on a hierarchical linear finite element mesh.
The AMG initialization procedure is described by Algorithm 1.For linear finite elements, the initial graph G 0 is the graph defined by the sparse matrix B SL .Alternatively it can be seen as the graph defined by the triangulation.For quadratic and cubic elements, the graph is defined by a refined triangulation in which quadratic elements are divided into four triangles, and cubic elements are divided into nine triangles.The reason not to use directly the graph defined by B SL for higher-order elements is that the coarsening procedure would coarsen the graph too much, leading to a slower convergence, that is, to not so well conditioned AB −1 SL .The nodes (vertices) in the graph G 0 associated with the Dirichlet boundaries are marked.On the coarser graphs also the nodes which were marked as Dirichlet nodes on the finer graphs are marked.
The nodes onto a coarser graph G k+1 are chosen from the nodes of G k as follows.Find the node in G k which has the smallest degree, that is, the smallest number of edges associated to it.If there are several such nodes, choose the first one according to the used node numbering.This node is included onto the graph G k+1 .Eliminate this node and all its neighbors from the graph G k .There is one exception to this: if the node has a Dirichlet marked neighbor then the neighbor is not eliminated from G k .This increases the stability of the procedure by making sure that there are sufficiently many Dirichlet nodes selected to the coarse levels.Repeat this procedure until there are no nodes left in G k .Figure 1 shows an example of this coarsening strategy when all boundaries are of Dirichlet type.
After choosing the nodes on G k+1 , they are numbered following their order on G k .Then the restriction matrix R k is defined by 1 for a fine node j which is a coarse node i, for a fine node j which is a neighbor of coarse node i and has k neighboring coarse nodes, 0 otherwise, where fine and coarse refers to the graphs G k and G k+1 , respectively.The edges of the coarse graph G k+1 are formed using the restriction matrix R k .Each coarse graph node corresponds to a row in the restriction matrix and there is an edge between two nodes if and only if the corresponding rows of the restriction matrix have a non-zero element in the same column.
The prolongation matrix is the transpose of the restriction matrix, and a coarser grid system matrix is constructed by Galerkin method, that is, the fine grid matrix is multiplied by R k from the left side and by R T k from the Presmooth Restrict the residual Set x l+1 = 0 and call µ times the cycle for the next level l + 1 7.
Prolong the correction x l = x l + R T l x l+1 8.
Post-smooth x l = x l + S l (x l , f l ) 9. End if right side.The AMG cycle described by Algorithm 2 is a usual multigrid cycle given here in the general µ-cycle form.The choices µ = 1 and µ = 2 correspond to V-cycle and W-cycle, respectively.When the algorithm is used in preconditioning it is called with the approximate solution x 0 being zero.

Model problems with homogenous medium
We use two different model geometries with homogenous medium: the unit square and a cross section of a car cabin.Figure 2 shows typical solutions for these geometries.The same problem in the unit square was considered also in [14].It has a point source at the middle and the absorbing boundary condition ( 5) is posed on all boundaries.The car cabin problem with a nonconvex geometry resembles more real-world applications.The height of the car cabin is 1.5 m and its width is 3 m.The noise source is modeled using the Dirichlet boundary condition (3) with g = 1 on the wall behind pedals and on other boundaries the impedance boundary condition (4) with γ = 0.2 is used.The meshes for the car cabin problem were generated using Netgen [29] by refining a coarse mesh depicted in Figure 1.

Eigenvalue study
We study the eigenvalues for both problems by computing both the full spectrum and Arnoldi approximations discussed in the end of Section 2. First we consider the Helmholtz problem with k = 20 in the unit square domain discretized on a 31 × 31 structured mesh.Figures 3 and 4 demonstrate the influence of β 2 , that is, the amount of damping to the spectrum of AB −1 SL .We use the values β 2 = 0.5 and β 2 = 1.0.In these and all following results we have β 1 = 1. Figure 3 shows that with the Neumann boundary condition the real parts of the eigenvalues are between zero and one with both β 2 s, while the density of eigenvalues near zero is higher with β 2 = 1.0.The differences are more pronounced with the absorbing boundary condition.With a smaller β 2 , the matrix AB −1 SL is closer to identity and this is seen as tighter clustering around one in Figure 4.The spectrum of AB −1 SL is very similar with quadratic and cubic elements to one with linear elements shown in these figures.
Next we study the quality of the AMG approximations of the inverse of the discrete shifted-Laplacians.For this we use one W-cycle based on one presmoothing and postsmoothing iteration performed by the underrelaxed Jacobi with the relaxation parameter chosen according to Table 2. around one than with quadratic elements which indicates that the AMG approximation of the inverse is better with linear elements.M G with the Neumann and absorbing boundary conditions together with their Arnoldi approximations.These plots show that with quadratic elements some of the eigenvalues have near zero or even negative real parts, while for linear and cubic elements the real parts are clearly positive.This shows up as higher number of iterations with quadratic elements in Section 5.1.2.
For the car cabin problem with k = 15, similar plots of eigenvalues are given in Figure 7.The figures show that the spectrums are fairly similar for the unit square and car cabin problems.This suggest that the quality of the AMG preconditioner is not particularly sensitive to the geometry.

Performance analysis
We mainly use the car cabin problem to study the performance of the iterative solver while we also report some results for the unit square problem.
Figure 2 shows usual time-harmonic scattering patterns for these problems.Table 1 describes the different meshes used with the car cabin problem.
The preconditioner B −1 M G is defined by one algebraic multigrid cycle described in Section 3. Our aim is to choose the parameters defining the preconditioner in such a way that the overall performance is optimal.We use the value β 1 = 1 while for β 2 we examine the values 0.5 and 1.0.There are several choices related to the AMG method.Our smoother is the underrelaxed Jacobi iteration with the relaxation parameter ω chosen according to Table 2.These relaxation parameters minimize the overall solution time in numerical tests.We use the W-cycle in the AMG method as it leads to faster solution times than the V-cycle in experiments.
In Table 3, the number of iterations are reported for the unit square problem.We have chosen the wave numbers k and the mesh step sizes h to be the same as in [14].Their discretization was performed using a finite difference method and they employed a tuned geometric multigrid method instead of the AMG method used here.Furthermore, they used the BICGSTAB method with a slightly more strict stopping criterion.The number of iterations required here are higher.When excluding quadratic elements we need up to 1.5 times more iterations with k = 40 while the difference grows with k.Nevertheless the results in here and in [14] suggest that for higher frequencies the number of iterations roughly doubles when the frequency is doubled.The AMG preconditioner leads to particularly good results with cubic finite elements.
The convergence results for the car cabin problem are presented in Tables 4-6.In these tables, the wave number doubles from a row to the next and the Table 2 The optimal choice of the Jacobi relaxation parameter ω for different finite elements and different values of β 2 found by extensive numerical experiments.mesh step size h is halved from a column to the next.This leads to constant khs on diagonals.Along them, we can again observe that for higher frequencies the number of iterations roughly doubles when the wave number is doubled.The lower triangles of the tables correspond to discretizations which do not have sufficiently high number of nodes per wavelength to capture the oscillatory behavior of solutions.This shows up as unusually high number of iterations.Based on these results, the value β 2 = 0.5 leads to much faster convergence on higher wave numbers.
With a large number of iterations, the time spent in forming the basis vectors in the GMRES method takes a larger part of the CPU time.This effect is seen in Table 7 which reports the time spent in different parts of the solver for three problems.These results were obtained with a computer which has Intel Core Duo U2500 processor.

Wedge problem with inhomogenous medium
The wedge problem is defined by three layers with different speed of sound c in the rectangle 600 × 1000 m 2 , as shown in Figure 8.This model problem was considered in [14,19].Meshes for different frequencies were constructed with Comsol Multiphysics mesh generator in such a way that the mesh step size h was approximately one tenth of one wavelength, that is, h ≈ λ/10.The wedge model has a point source at the middle of the top boundary.The first-order absorbing boundary condition is posed on all boundaries.A coarse mesh for the frequency f = 5 Hz and the solutions for the frequencies f = 30 Hz and f = 50 Hz are shown in Figure 8.
The performance results for different frequencies are presented in Table 8.
According to these iteration counts, the convergence with β 2 = 0.5 is about 15 % faster than with β 2 = 1.0 with linear elements and over 40 % faster with cubic elements.We need again more iterations when compared to the results in [14], but our finite element discretizations have three advantages over the finite differences used in there: we can accurately model the interface between layers, we can use coarser meshes with cubic elements, and we can use coarser mesh where the speed of sound is higher.Figure 9 plots the eigenvalues of the preconditioned system for the low frequency f = 5 Hz with β 2 = 1.0 and β 2 = 0.5.

Conclusions
We have studied a preconditioner based on an algebraic multigrid (AMG) approximation of the inverse of a shifted-Laplacian for the Helmholtz equation.This is a generalization of the preconditioner proposed by Erlangga, Oosterlee, and Vuik in [14].They used finite difference discretizations on rectangular domains and a geometrical multigrid.With our finite element discretizations and the AMG method we can solve problems in complicated domains and use higher-order finite elements.A big advantage of the AMG method is that the solver does not need hierarchical meshes nor operators discretized on different meshes.When the matrix for the zeroth-order term in a discretized Helmholtz equation or the mass matrix for a constant wave number problem is also available, the preconditioner can be constructed fully algebraically.Thus, in this case the preconditioned iteration can be seen as a "black box solver".
The numerical results demonstrated the capability to solve efficiently problems in complicated domains and varying wave numbers using the proposed preconditioner.Furthermore, the preconditioner was shown to be effective with linear, quadratic, and cubic finite elements.The proposed approach is especially well-suited for low-frequency and mid-frequency problems while for high-frequency problems the number of iterations roughly doubles when the frequency is doubled.The same behavior was also observed in [14].

Figure 1 .
Figure1.A mesh for a cross-section of a car cabin.Gray nodes are selected to next coarse level.

Figure 2 .
Figure 2. In the left plot a solution for the unit square problem.In the right plot the solution for the car cabin problem with the wave number k = 18.3 which corresponds to the frequency f ≈ 1 kHz.

Figure 5 plotsFigure 3 .Figure 4 .
Figure 3.The eigenvalues of AB −1 SL with β 2 = 1.0 in the left plot and β 2 = 0.5 in the right plot for the unit square with the Neumann boundary condition discretized using linear finite elements.

Figure 6
Figure 6 depicts the eigenvalues of AB −1M G with the Neumann and absorbing boundary conditions together with their Arnoldi approximations.These plots show that with quadratic elements some of the eigenvalues have near zero or even negative real parts, while for linear and cubic elements the real parts are clearly positive.This shows up as higher number of iterations with quadratic elements in Section 5.1.2.

Figure 5 .
Figure 5.The eigenvalues of B SL B −1 M G for the unit square problem with the absorbing boundary condition.The left and right plots are computed using linear and quadratic elements, respectively.

Figure 6 .
Figure 6.For the unit square problem the eigenvalues of AB −1M G are marked with • and their Arnoldi approximations with ×.The left and right plots are based on β 2 = 1.0 and β 2 = 0.5, respectively.The discretizations have been performed using linear (plots on the top), quadratic (plots in the middle), and cubic (plots on the bottom) finite elements.

Figure 7 .
Figure 7.For the car cabin problem discretized with linear elements the eigenvalues of AB −1 M G marked with • and their Arnoldi approximations marked with ×.The left and right plots are based on β 2 = 1.0 and β 2 = 0.5, respectively.

Figure 8 . 5 β 2 = 1
Figure 8.The plot on the left shows the mesh for the frequency f = 5 Hz and the definition of the problem.On the middle and right the solutions for f = 30 Hz and f = 50 Hz, respectively, are shown.Table 8 For the wedge problem the number of GMRES iterations required to reduce the norm of the residual by the factor 10 −6 as a function of the frequency.element type f = 5 Hz f = 30 Hz f = 50 Hz linear, β 2 = 1.0 24 97 148 linear, β 2 = 0.5 17 83 124 cubic, β 2 = 1.0 24 105 171 cubic, β 2 = 0.5 19 62 97

Figure 9 .
Figure 9.The eigenvalues of AB −1 M G marked with • and their Arnoldi approximations marked with × for the wedge problem with the frequency f = 5 Hz discretized using linear elements.The left and right plots are based on β 2 = 1.0 and β 2 = 0.5, respectively.
Algorithm 1 AMG initialization Input: Matrix B 0 , initial graph G 0 , the maximum size of the coarsest system n c 1. k = 0 2. Do while the size of B k is greater than n c 3. Select the set of coarse nodes from the graph G k Eliminate the rows and columns of B k marked in the graph G k 13.Factorize B k Algorithm 2 Recursive algorithm for the AMG cycle.Input: Matrix B l , approximate solution x l , right-hand side vector f l Output: Improved approximate solution x l 1.If on the coarsest level, that is, l

Table 3
The number of iterations for the unit square problem for different element types when the iterations are terminated once the norm of the residual is reduced by the factor 10 −6 .

Table 4
For the car cabin problem discretized with linear elements, the number of GM-RES iterations required to reduce the norm of the residual by the factor 10 −6 as a function of the number of refinements n r and the wave number k.

Table 5
For the car cabin problem discretized with quadratic elements, the number of GM-RES iterations required to reduce the norm of the residual by the factor 10 −6 as a function of the number of refinements n r and the wave number k.

Table 6
For the car cabin problem discretized with cubic elements, the number of GMRES iterations required to reduce the norm of the residual by the factor 10 −6 as a function of the number of refinements n r and the wave number k.

Table 7
The time spent in the solver with the car cabin problem discretized using quadratic elements on the n r = 6 mesh.The AMG uses W-cycle and β 2 = 1.0.