Comparison of Algebraic Multigrid Preconditioners for Solving Helmholtz Equations

Abstract

An algebraic multigrid (AMG) with aggregation technique to coarsen is applied to construct a better preconditioner for solving Helmholtz equations in this paper. The solution process consists of constructing the preconditioner by AMG and solving the preconditioned Helmholtz problems by Krylov subspace methods. In the setup process of AMG, we employ the double pairwise aggregation (DPA) scheme firstly proposed by Y. Notay (2006) as the coarsening method. We compare it with the smoothed aggregation algebraic multigrid and meanwhile show shifted Laplacian preconditioners. According to numerical results, we find that DPA algorithm is a good choice in AMG for Helmholtz equations in reducing time and memory. Spectral estimation of system preconditioned by the three methods and the influence of second-order and fourth-order accurate discretizations on the three techniques are also considered.

1. Introduction

In this paper, the time-harmonic wave equation in 2D homogeneous media is solved numerically. The essential equation is Helmholtz equation, which governs wave scattering and propagation phenomena arising in acoustic problems in many areas, such as geophysics, aeronautics, and optical problems. In particular, we are in search of solutions of the Helmholtz equation discretized by the finite difference method. The discrete problem becomes extremely large for very high wavenumber, because the number of gridpoints per wavelength should be sufficiently large in order to result in accepted solutions. In this case, direct methods are difficult to solve, and iterative methods are the interesting alternative. However, Krylov subspace iterative methods are not competitive without a good preconditioner. In this paper, we consider an algebraic multigrid with aggregation scheme as preconditioning to improve the convergence of Krylov subspace iterative methods.

In [1], Bayliss et al. proposed a preconditioner based on the Laplace operator for solving the discrete Helmholtz equation efficiently with CGNR. In [2], Laird and Giles proposed a preconditioner where an extra term is added to the Laplace operator for solving the discrete Helmholtz equation. Subsequently, in [3], Erlangga et al. generalized the above two kinds of preconditioners and obtained a new class of preconditioners, the so-called “shifted Laplacican” preconditioner of the form −Δ𝜙−𝛼𝑘²𝜙 with 𝛼=−(𝛽₁−𝛽₂𝑖)∈ℂ, where 𝑖=√−1 is the imaginary unit. In 2006, Erlangga et al. [4] compared multigrid with incomplete LU used to approximate the shifted Laplacian preconditioner to construct the final preconditioner for the inhomogeneous Helmholtz equation, and concluded that multigrid applied to approximate the shifted Laplacian preconditioner with Bi-CGSTAB resulted in a fast and robust method. In 2006, Erlangga et al. [5] proposed a multigrid 𝑉(1,1)-cycle with de Zeeuw’s prolongation operator, FW (full weighted) restriction, and Jacobi smoothing with the relaxation parameter 𝜔=0.5. The smallest size of the parameter 𝛽₂ in front of the imaginary Helmholtz term in the preconditioner, for which the multigrid method could be successfully employed, has been determined to 𝛽₂=0.5 and meanwhile 𝛽₁=1. In 2007, Van Gijzen et al. [6] analyzed the spectral of the discrete Helmholtz operator preconditioned with a shifted Laplacian, and proposed an optimal value for the shift by combining the results of the spectral analysis with an upper bound on the GMRES residual norm. In 2009, Umetani et al. [7] proposed a multigrid 𝑉(0,1)-cycle with AMG’s prolongation operator, FW restriction, and incomplete LU postsmoothing for a fourth-order Helmholtz discretization based on [5]. The fourth-order accurate shifted Laplacian preconditioner could be easily approximated by one 𝑉(0,1)-cycle of multigrid and enables us to choose a somewhat smaller imaginary shift parameter (𝛽₂=0.4) in the shifted Laplacian preconditioner, which improves the solvers convergence (especially for high wavenumbers on fine meshes). In 2007, Airaksinen et al. [8] proposed a preconditioner based on an algebraic multigrid approximate of the inverse of a shifted Laplacian for the Helmholtz equation. This is a generalization of the preconditioner proposed by Erlangga et al. in [5]. In 2010, Olson and Schroder [9] proposed a smoothed aggregation algebraic multigrid method for 1D and 2D scalar Helmholtz problems with exterior radiation boundary conditions discretized by 1D finite difference and 2D discontinues Galerkin. In the same year, Notay [10] proposed an aggregation-based algebraic multigrid (AGMG) method, in which double pairwise aggregation scheme firstly proposed by Notay in [11] is used.

We will consider using the double pairwise aggregation algorithm to coarsen in the setup process of the general algebraic multigrid method in this paper in order to improve the solution effects of Helmholtz problems.

This paper is organized as follows. In Section 2, we discuss the Helmholtz equation, and the discrete finite difference formulations of second order and fourth order. The iterative solution methods with the three different kinds of preconditionings are presented in Section 3. Numerical results are reported in Section 4, and conclusions are given in Section 5.

2. The Helmholtz Equation and Discretizations

2.1. Mathematical Problem Definition

We solve wave propagations in a two dimensional medium in a unit domain governed by the Helmholtz equation−Δ𝜙−𝑘²𝑥,𝑦𝜙=𝑔,Ω=[0,1]², (2.1) where Δ≡𝜕²/𝜕𝑥²+𝜕²/𝜕𝑦² is the Laplace operator, 𝜙(𝑥,𝑦) represents the pressure field in the frequency domain, the source term is denoted by 𝑔, and 𝑘(𝑥,𝑦) is the wavenumber, which is a constant in the homogeneous domain. Otherwise, the wavenumber 𝑘=𝜔/𝑐(𝑥,𝑦) is space dependent because of a spatially dependent speed of sound 𝑐(𝑥,𝑦) in a heterogeneous medium. With angular frequency 𝜔=2𝜋𝑓  (𝑓 is the frequency in Hertz), wavelength ℓ is defined by ℓ=𝑐/𝑓. So the wavenumber 𝑘=2𝜋/ℓ. The number of wavelengths in a domain of size 𝐿 equals 𝐿/ℓ. 𝑛_𝜔, the number of points per wavelength, is typically chosen to be 10–12 points. A dimensionless discretization step reads ℎ=ℓ/(𝑛_𝜔𝐿), and therefore 𝑘ℎ=2𝜋/(𝑛_𝜔𝐿). With domain size 𝐿=1, an accuracy requirement for second-order discretizations is that 𝑘ℎ≤𝜋/5(≈0.63) for 𝑛_𝜔=10 points per wavelength, and 𝑘ℎ≤0.53 with 𝑛_𝜔=12 points per wavelength. We select combinations of wavenumber 𝑘 and the discretization step ℎ with 𝑘ℎ=0.625, which is displayed in Table 1, in the latter numerical experiments.

Table 1 

Discretization steps related to wavenumbers.

The boundary condition can be the Dirichlet boundary condition or the first-order radiation boundary condition, and they are, respectively, as follows:𝜙=0,onΓ=𝜕Ω,Dirichletboundarycondition𝜕𝜙(),𝜕𝑛−𝑖𝑘𝜙=0,onΓ=𝜕Ω,Sommerfeldcondition (2.2)(2.3) where 𝑛 is an outward direction normal to the boundary.

2.2. Finite Difference Discretizations

The Helmholtz equations are discretized either by a second-order or by a fourth-order finite difference scheme, resulting in the linear system𝐴_ℎ𝜙_ℎ=𝑏_ℎ, (2.4) where 𝜙_ℎ and 𝑏_ℎ represent the discrete frequency-domain pressure field and the source, respectively.

The well-known 5-point stencil related to a second-order (𝑜(ℎ²)) accurate discretization for Helmholtz equation (2.1) reads𝐴_ℎ≜1ℎ²()−1−14−𝑘ℎ²−1−1. (2.5)

The fourth-order (𝑜(ℎ⁴)) accurate discretization named by HO discretization [12] for Helmholtz equation (2.1) is given with stencil𝐴^HO_ℎ≜1ℎ²−16−23−(𝑘ℎ)²−1126−23−(𝑘ℎ)²−12103−2(𝑘ℎ)²3−23−(𝑘ℎ)²−1126−23−(𝑘ℎ)²−1126, (2.6)

Compared with a second-order discretization method, a fourth-order discretization method could make the number of grid points per wavelength be smaller and discrete the equations to smaller matrices for the same level of accuracy. It is possible to lead to an algorithm that is more efficient in terms of the accuracy of the solution versus the computational cost. As in [7], the multigrid-based shifted Laplacian preconditioner with the fourth-order discretization is preferable.

The boundary discretization is also considered, and we could employ the first-order forward or backward scheme to approximate the first-order derivative in (2.3).

Next, we can obtain a linear system through substituting𝐴_ℎ𝑥=𝑏_ℎ,𝐴_ℎ∈ℂ^𝑁×𝑁, (2.7) where 𝐴_ℎ is a large, spare symmetric matrix and 𝑁 is the number of grid points. The matrix 𝐴_ℎ remains positive definite as long as 𝑘² is smaller than the minimum eigenvalue of the discrete Laplacian, is indefinite with both positive and negative eigenvalues for large values of 𝑘, and is complex-valued because of absorbing boundary conditions. However, it is non-Hermitian. Apparently, the size of the linear system (2.7) gets very large for the high frequency.

3. Iterative Solution Methods

3.1. Preconditioned Krylov Subspace Iterative Methods

Iterative methods for linear system (2.7) within the class of Krylov subspace method are based on the construction of iterants in the subspace𝐊^𝑗𝐴,𝑟₀=𝐬𝐩𝐚𝐧𝑟₀,𝐴𝑟₀,𝐴²𝑟₀,…,𝐴^𝑗−1𝑟₀, (3.1) where 𝐊^𝑗(𝐴,𝑟₀) is the 𝑗th Krylov subspace associated with 𝐴 and 𝑟₀, 𝑟₀=𝑏−𝐴𝑥₀, is the initial residual.

The Bi-CGSTAB algorithm [13] is one of the better known Krylov subspace algorithms for non-Hermitian problems, which has been used for Helmholtz problems, for example, in [5, 8]. One of the advantages of Bi-CGSTAB, compared to full GMRES, is its limited memory requirements. Bi-CGSTAB is based on the idea of computing two mutually biorthogonal bases for the Krylov subspaces based on matrix 𝐴_ℎ and its conjugate transpose 𝐴^𝐻_ℎ and is easy to implement.

Without a preconditioner, Krylov subspace methods converge very slowly, or not at all, for the problem of interest in [4]. So a preconditioner should be incorporated to improve the convergence of Krylov subspace methods. By left preconditioning, one solves a linear system premultiplied by a preconditioning matrix 𝑀_ℎ⁻¹,𝑀_ℎ⁻¹𝐴_ℎ=𝑀_ℎ⁻¹𝑏_ℎ. (3.2)

The challenge is to find a form of matrix 𝑀_ℎ, whose inverse matrix 𝑀_ℎ⁻¹ can be efficiently approximated, such that 𝑀_ℎ⁻¹𝐴_ℎ has a spectrum that is favorable for Krylov subspace iterative solution methods.

In this paper, we choose Bi-CGSTAB and GMRES with preconditionings to solve discrete Helmholtz equations. In [3], GMRES is used to solve the Helmholtz equation with the shifted Laplacian preconditioner and is compared with Bi-CGSTAB. As a result, Bi-CGSTAB is preferable since the convergence for heterogeneous high wavenumber Helmholtz problems is typically faster than that of GMRES.

3.2. Shifted Laplacian (SL) Preconditioners

In [5], a shifted Laplacian operator was proposed as a preconditioner for the Helmholtz equation, with 𝑀_ℎ defined as a discretization of𝑀_SL=−𝜕𝑥𝑥−𝜕𝑦𝑦+𝛼𝑘²𝑥,𝑦,𝛼∈ℂ,𝛼=−𝛽₁−𝛽₂𝑖, (3.3) and boundary conditions were set identically to those for the original Helmholtz equation. The readers could refer to the related contents in [4, 7, 8, 14].

There are different choices of the parameter pair (𝛽₁,𝛽₂) such as (𝛽₁,𝛽₂)=(0,0) (Laplace preconditioner in [1]), (𝛽₁,𝛽₂)=(−1,0) (Laird and Giles preconditioner in [2]), (𝛽₁,𝛽₂)=(0,1) (Erlangga et al. preconditioner in [3]), (𝛽₁,𝛽₂)=(1,1) (basic parameter pair choice). In [5], Erlangga et al. evaluated the influence of parameters 𝛽₁ and 𝛽₂ to the multigrid method which was applied to approximate the shifted Laplacian operator to construct the final preconditioner, and the proposed optimal parameter pair for the solver was that (𝛽₁,𝛽₂)=(1,0.5). Based on [5], Umetani et al. [7] considered using the fourth-order accurate finite difference discretization with the result that the parameter pair (𝛽₁,𝛽₂)=(1,0.4) is the best choice. In this paper, we merely consider the shifted Laplacian preconditioners, without the multigrid or incomplete LU method to approximate, which is not the same as that in [5, 7]. So we will also choose that (𝛽₁,𝛽₂)=(1,0.5) and (𝛽₁,𝛽₂)=(1,0.3) for comparison with the previously mentioned parameter pair (𝛽₁,𝛽₂) in the numerical experiments.

3.3. Smoothed Aggregation (SA) Preconditioning

In this section, we will present an overview of the smoothed aggregation setup algorithm. The function of the SA setup phase is to construct coarse operators 𝐴_𝑘 through interpolation operators𝑃_𝑘∶ℜ^𝑛𝑘+1⟶ℜ^𝑛𝑘, (3.4) where 𝑛_𝑘 and 𝑛_𝑘+1 are sizes of two successively coarser grids. The initial matrix 𝐴₀ associated with the finest level equals 𝐴. The inputs need the matrix 𝐴 and one near null-space candidate 𝑡 that is intended to form the interpolation basis. 𝑡 is an algebraically smooth mode that is slow to relax on the fine mesh. The details are described in Algorithm 1.

(1) 𝑡ℓ=𝑆𝜈1ℓ(ones(𝑛,1),0) (2) While ℓ≤𝐿 (3) [neighbour, hoods] = aggregate (𝑛ℓ,𝐴ℓ); (4) [𝑃,𝑡ℓ+1] = qr(𝑛ℓ,𝑡ℓ,hoods); (5) 𝑆ℓ= smoother (𝑛ℓ,𝐴ℓ,𝜔,neighbour); (6) 𝐼ℓℓ+1=𝑆ℓ𝑃; (7) 𝐴ℓ+1=(𝐼ℓℓ+1)𝑇𝐴ℓ𝐼ℓℓ+1; (8) ℓ=ℓ+1; (9) Endwhile

Algorithm 1 

SA_setup (𝑛_ℓ,𝐴_ℓ,ℓ).

When the current coarse level ℓ≤𝐿(𝐿 is the coarsest level), we gain two sets respectively named by “neighbour” and “hoods.” The “neighbor [𝑖]” includes the nodes which are strongly connective to the node 𝑖, and the “hoods [𝑖]” contains the nodes which are in the 𝑖th aggregation set. Next, we reset 𝑡_ℓ according to “hoods,” and QR factorization is applied to it in order to get the initial interpolation operator 𝑃 and the next 𝑡_ℓ+1. Third, we obtain the smooth operator 𝑆_ℓ according to the relative algorithm in [15]. Fourth, we get the final interpolation operator 𝐼^ℓ_ℓ+1, through the smooth operator 𝑆_ℓ is acted on the initial interpolation operator 𝑃. Finally, we can easily obtain the coarse operator 𝐴_ℓ+1. The interested readers could consult the related contents in [9]. We will employ the SA setup algorithm in AMG to approximate the original system in order to obtain the preconditioner and we name this process by SA preconditioning.

3.4. Double Pairwise Aggregation (DPA) Preconditioning

An algebraic coarsening algorithm sets up an interpolation operator 𝑃 only making use of the information available in the fine grid matrix 𝐴. The interpolation operator 𝑃 is a 𝑛×𝑛_𝑐 matrix, where 𝑛_𝑐(𝑛_𝑐≤𝑛) is the number of coarse variables. With acted on by the interpolation operator, a vector defined on the coarse variable set [1,𝑛_𝑐] could be transferred on the fine grid. The coarsen grid matrix 𝐴_𝑐 depends not only on the interpolation matrix 𝑃, but also on the restriction matrix 𝑅 and the fine grid matrix 𝐴. It is usual to take the restriction operator 𝑅 equals to the transpose of the interpolation operator 𝑃. Therefore, the coarsen grid matrix 𝐴_𝑐 is computed from the Galerkin formula𝐴_𝑐=𝑅𝐴𝑃=𝑃^𝑇𝐴𝑃. (3.5)

Coarsening by aggregation needs to define aggregates 𝐺_𝑖, which are disjoint subsets of the variable set. The number of these aggregations is the number of coarse variables 𝑛_𝑐, and the interpolation operator 𝑃 is given by𝑃_𝑖,𝑗=1,if𝑖∈𝐺_𝑗,0,otherwise1≤𝑖≤𝑛,1≤𝑗≤𝑛_𝑐. (3.6)

Note that there is no need to explicitly form interpolation operator 𝑃, and the coarse grid matrix is practically computed by𝐴_{𝑐𝑖,𝑗}=_{𝑘∈𝐺𝑖ℓ∈𝐺𝑗}𝑎_𝑘ℓ1≤𝑖,𝑗≤𝑛_𝑐. (3.7)

The double pairwise aggregation algorithm was first proposed by Notay in [11], which employed two passes of the pairwise aggregation algorithm, with the result of decreasing the number of variables by a factor slightly less than four in most cases. Subsequently, we will take the double pairwise aggregation scheme in the setup process of AMG, in order to construct a better preconditioner with Bi-CGSTAB for discrete Helmholtz equations. We also name this process by DPA preconditioning. The details of aggregation algorithms could be consulted in [10, 11].

4. Experiments

We will consider the following problem and use a uniform mesh with constant mesh size ℎ in all directions. In the experiments, we will compare the three preconditionings combined with Bi-CGSTAB or GMRES for solving discrete Helmholtz equations. The shifted Laplacian preconditioner (SLP), the smoothed aggregation preconditioning (SAP), and the double pairwise aggregation preconditioning (DPAP) are, respectively, introduced in Section 3. We let that the iterative method is terminated at the 𝑘th step if ‖𝑏−𝐴𝑥_𝑘‖₂/‖𝑏‖₂≤10⁻⁶. All the numerical results are from running relative Matlab programs in the computer with the CPU of AMD Pentium Dual Core 4000+.

4.1. The Close-off Problem and Numerical Results

We consider a problem in a rectangular homogeneous medium governed by−Δ𝜙−𝑘²𝜙=sin(𝜋𝑥)sin𝜋𝑦sin√2𝜋𝑥sin√3𝜋𝑦,𝑥=[0,1],𝑦=[0,1], (4.1)𝜙=0, at the boundaries.

Different grid resolutions are used to solve this problem with various wavenumbers in Table 1.

Firstly, we will observe the distribution of eigenvalues of the systems preconditioned by the three different preconditionings in order to estimate which preconditioning is best combined with Krylov subspace iterative methods. Since the spectral analysis of shifted Laplacian preconditioners (𝑀_SL)_ℎ with different parameter pairs (𝛽₁,𝛽₂) was considered in many papers [5–7, 14] and (𝛽₁,𝛽₂)=(1,0.5) is an almost better choice, we will only present the distribution of eigenvalues of the system preconditioned by the shifted Laplacian preconditioner (𝑀_SL)_ℎ with (𝛽₁,𝛽₂)=(1,0.5) with second-order and fourth-order accurate discretizations. We employ two-gird process to approximate 𝑉-cycle of AMG, and so 𝑀⁻¹_AMG𝐴=𝐼−𝑆(𝐼−𝑃𝐴_𝐶⁻¹𝑃^𝑇𝐴)𝑆, where 𝐼 is unit matrix and 𝑆 is smoothing matrix. Second-order and fourth-order accurate discretizations are taken into account. Next, we note preconditioners produced by the smoothed aggregation preconditioning and by the double pairwise aggregation preconditioning separately as (𝑀_SA)_ℎ and (𝑀_DPA)_ℎ. In Figures 1(a), 1(b), and 1(c) are, respectively, spectral pictures of (𝑀_SL)_ℎ⁻¹𝐴_ℎ, (𝑀_SA)_ℎ⁻¹𝐴_ℎ and (𝑀_DPA)_ℎ⁻¹𝐴_ℎ with second-order accurate discretization, and Figures 1(d), 1(e), and 1(f) are respectively spectral pictures of (𝑀_SL)_ℎ⁻¹𝐴_ℎ, (𝑀_SA)_ℎ⁻¹𝐴_ℎ, and (𝑀_DPA)_ℎ⁻¹𝐴_ℎ with fourth-order accurate discretization. Meanwhile, we compute the condition number of the preconditioned systems, and cond₂((𝑀_SL)_ℎ⁻¹𝐴_ℎ)=211.1400, cond₂((𝑀_SA)_ℎ⁻¹𝐴_ℎ)=1.3062, and cond₂((𝑀_DPA)_ℎ⁻¹𝐴_ℎ)=1.3078 separately correspond to Figures 1(a), 1(b), and 1(c).

(a)

(b)

(c)

(d)

(e)

(f)

Figure 1 

Spectral pictures.

Next, we will present numerical results in Tables 2, 3, and 4 and explain the meaning of signs in them. Tables 2, 3, and 4, respectively, show the computational performance of GMRES with SLPs, Bi-CGSTAB with SAP, and Bi-CGSTAB with DPAP for solving the closed-off problem with finite difference discretizations, in terms of the number of iterations and computational time, to reach the specified convergence. We also test Bi-CGSTAB with SLPs, but it leads to more iterations and solve time. So we do not present the corresponding results. Meanwhile, we find that Bi-CGSTAB with SLPs when 𝛽₁=1 and 𝛽₂ in the interval (0,1) is better than GMRES, and the fourth-order accurate discretization could make the iterative solution rate more quick. Nevertheless, we make the focus only on the second-order accurate discretization for shifted Laplacian preconditioners, so the second-order accurate discretization is used in Table 2. “Iter” denotes the number of iterations, and “—” means that the corresponding experimental result is not given. “𝑇_setup (s),” and “𝑇_solve (s)” respectively, denote the setup time of preconditioner and the time of iterative solution by the second, and “𝑇_𝑂” and “𝐹_𝑂,” respectively, mean the second-order and fourth-order finite difference discretizations shown in Section 2. The values of “presmooth” and “postsmooth” respectively, denote the number of pre- and postsmoothing iterations in the solve process of AMG, “smoother” denotes the chosen smoothing method, and “setup_Alg” denotes the coarsening algorithm used in the setup process of AMG.

Table 2 

Computational performance of GMRES with SLPs.

Table 3 

Computational performance of Bi-CGSTAB with SAP.

Table 4 

Computational performance of Bi-CGSTAB with DPAP.

4.2. Analysis of Numerical Results

We will analyze the numerical results from the following aspects.

Firstly, from Figure 1, by comparing pictures (a), (b), (c), we find that the eigenvalues of (b) and (c) are more farther from zero point and more close with one, but the eigenvalues of (b) contain zero point. So we could boldly deduce that DPAP and SAP are more beneficial for Krylov subspace iterative method. Perhaps DPAP is better than SAP since the eigenvalues of (b) contain zero point. Next, we contrast (a), (b), (c), respectively, with (d), (e), (f). We notice that the higher accurate discretization is possible to improve the solution efficiency of SLP, since the fourth-order discretization makes eigenvalues more close to one in (d), is probably good for SAP for avoiding the zero point eigenvalue, and is likely bad for DPAP because of appearing the zero point eigenvalue. As we all know, the larger the condition number of matrix is, the worse the numerical stability of solution of the corresponding matrix equation is. So the system preconditioned by the shifted Laplacian operator has worse numerical stability, which means that the small change of the right hand will bring large change of solution, since cond₂((𝑀_SL)_ℎ⁻¹𝐴_ℎ)=211.1400 is much larger. It is perhaps the reason for using multigrid or LU to approximate the shifted Laplacian operator to obtain final preconditioner, such as [4, 5, 7, 8].

Secondly, from Table 2, we could easily see that the number and time of iterative solution are both going up largely with the increasing wavenumber. The shifted Laplacian preconditioner degenerates into the Laplace preconditioner when (𝛽₁,𝛽₂)=(0,0). Next, we find out that the number of iterations is smallest when (𝛽₁,𝛽₂)=(1,0.3), but it is not best in terms of the time of iterative solution. In general, the best choice is that (𝛽₁,𝛽₂)=(1,0.5) by considering both aspects of the number and time of iterative solution.

Thirdly, from Table 3, we find that the fourth-order accurate discretization only affects the iterative solution a little, and it is little better than the second-order accurate discretization when the wavenumber 𝑘=10,20,40. However, only the low wavenumbers 𝑘 are computed in Table 3 because of the too long time of iterative solution or the inadequate memory for higher wavenumbers. From Table 4, we discover that the fourth-order accurate discretization does not evidently improve the efficiency of iterative solution and that the second-order accurate discretization looks more stable in the aspect of the number of iterations. The fourth-order accurate discretization improves the efficiency of iterative solution both on the sides of setup-time and solve time when the wavenumber 𝑘=10,20 but costs more time and iterations for the higher wavenumber 𝑘=30,…,90. So we could conclude that the second-order accurate discretization combined with DPAP is preferable.

Finally, to compare results of Tables 2, 3, and 4, according to the number of iterations, we can easily find that Bi-CGSTAB with DPAP and second-order accurate discretization is preferable. However, presmooth =0 is chosen in Table 3, and, when we take presmooth =1 for Bi-CGSTAB with SAP, the iterative steps will be smaller but still not better than Bi-CGSTAB with DPAP, with leading to the longer solve time for the added smoothing process. Following, we contrast solve time and setup time with second-order accurate discretization in Tables 2, 3, and 4. Then, we see that GMRES with SLP with (𝛽₁,𝛽₂)=(−1,0) need the least solve-time for wavenumber 𝑘=10, Bi-CGSTAB with SAP has the least time of iterative solution for wavenumber 𝑘=20, and Bi-CGSTAB with DPAP has the quickest iterative rate for wavenumber 𝑘=30,40. SAP has much more setup time than DPAP for the same wavenumber, and setup time of SAP increases more rapidly than that of DPAP. Though setup time will lead to the whole solution time longer, if there are some matrix equations with the same left matrix and the different right hands, DPAP will excel for higher wavenumber because of needing the setup process once. The results of preconditioning for wavenumber 𝑘=50,…,90 are only shown in Table 4, and the reason is that the larger the wavenumber 𝑘 is, the longer the solution time is for higher wavenumber. So we do not present the results for the higher wavenumber in Tables 2 and 3 due to the limited experiment condition. Next, we observe that, with the increasing wavenumber 𝑘, the number of iterative solution changes between linearly and squarely in Table 2, but it fluctuates between 4 and 7 in Table 3 and between 2 and 4 in Table 4. In other words, the number of iterative solution is not related to the wavenumber 𝑘 when Bi-CGSTAB with DPAP and SAP is employed. In contrast with Bi-CGSTAB with SAP in Table 3, Bi-CGSTAB with DPAP in Table 4 makes the time of iterative solution change more slowly while the wavenumber 𝑘 is increasing.

5. Conclusions

As is known to all, Helmholtz equations are difficult to solve, especially with high wavenumbers. The reason is that the high wavenumber is able to make the corresponding matrices of Helmholtz equations highly indefinite. The original Krylov subspace iterative method is not efficient to solve these problems. Preconditioning is needed to apply into Krylov subspace iterative methods in order to improve the efficiency of solution. We consider using the algebraic multigrid method to construct the preconditioner and look for a concrete AMG method which is fit for Helmholtz problems. According to numerical results and analysis in Section 4, the double pairwise aggregation is applied in the setup process of AMG and finally better results are gained. Meanwhile, the number of iterative solution is independent on the wavenumber, with the result that Helmholtz equations with high wavenumbers can be solved more efficiently in terms of time and memory. The methods that are Bi-CGSTAB with DPAP and with SAP get fewer number of iterative solution however, SAP prolongs the setup time in contrast with DPAP. The two kinds of preconditionings are obviously more efficient than the shifted Laplacian preconditioners. In view of the above, we can draw the conclusion that the double pair aggregation algorithm in the setup process of AMG is a good choice to solve Helmholtz equations especially with high wavenumbers. Though only the 2 D close-off problem with Dirichlet boundary condition is experimented in Section 4, this conclusion is also possible to generalize to Helmholtz equations with other boundary conditions.

Acknowledgments

This research is supported by NSFC (60973015, 61170311), Sichuan Province Sci. and Tech. Research Project (2011JY0002, 12ZC1802), the tianyuan foundation. Chinese Universities Specialized Research Fund for the Doctoral Program (20110185110020).