A multigrid-based preconditioned solver for the Helmholtz equation with a discretization by 25-point difference scheme

Available online at www.sciencedirect.com

ScienceDirect Mathematics and Computers in Simulation 117 (2015) 54–67 www.elsevier.com/locate/matcom

Original articles

A multigrid-based preconditioned solver for the Helmholtz equation with a discretization by 25-point difference scheme✩ Dongsheng Cheng a , Zhiyong Liu b,∗ , Tingting Wu c a School of Software Engineering, Shenzhen Institute of Information Technology, Shenzhen 518172, PR China b Industry center, Shenzhen Polytechnic, Shenzhen, 518055, PR China c Department of Mathematics, Shandong Normal University, Jinan, PR China

Received 4 July 2012; received in revised form 25 November 2013; accepted 12 January 2015 Available online 19 June 2015

Abstract In this paper, a preconditioned iterative method is developed to solve the Helmholtz equation with perfectly matched layer (Helmholtz-PML equation). The complex shifted-Laplacian is generalized to precondition the Helmholtz-PML equation, which is discretized by an optimal 25-point finite difference scheme that we presented in Chen et al. (2011). A spectral analysis is given for the discrete preconditioned system from the perspective of linear fractal mapping, and Bi-CGSTAB is used to solve it. The multigrid method is employed to invert the preconditioner approximately, and a new matrix-based prolongation operator is constructed in the multigrid cycle. Numerical experiments are presented to illustrate the efficiency of the multigrid-based preconditioned Bi-CGSTAB method with the new prolongation operator. Numerical results are also given to compare the performance of the new prolongation operator with that of the prolongation operator based on the algebraic multigrid (AMG) principle. c 2015 International Association for Mathematics and Computers in Simulation (IMACS). Published by Elsevier B.V. All rights ⃝ reserved. Keywords: Helmholtz equation; Preconditioner; Multigrid; Prolongation operator

1. Introduction The Helmholtz equation is ubiquitous, and it finds important applications in acoustics, electromagnetics, quantum mechanics, and geophysics. The Helmholtz equation is so important that its numerical simulation has been under active research. Due to the limited memory and speed of computers, the Helmholtz problem should be solved in a finite domain. For this end, artificial boundary condition methods are commonly used, such as the “absorbing boundary condition” (ABC, cf. [10,13]) and perfectly matched layer (PML, cf. [3,16]). The PML gradually damps the outgoing waves in absorbing ✩ This research is partially supported by the National Natural Science Foundation of China under grants 11101163, 11226105, 71101096 and 11301310, China Postdoctoral Science Foundation under grant 2013M531884, Science Research Project of Shenzhen Polytechnic under grant 601522k35010 and Opening Project of Guangdong Province Key Laboratory of Computational science at the Sun Yat-sen University under the grant 201206008. ∗ Corresponding author. E-mail addresses: [email protected] (D. Cheng), [email protected] (Z. Liu), [email protected] (T. Wu).

http://dx.doi.org/10.1016/j.matcom.2015.01.009 c 2015 International Association for Mathematics and Computers in Simulation (IMACS). Published by Elsevier B.V. All rights 0378-4754/⃝ reserved.

D. Cheng et al. / Mathematics and Computers in Simulation 117 (2015) 54–67

55

layers and gives a better performance than absorbing boundary condition, whose artificial boundary has zero width. To discretize the Helmholtz equation, we have mainly the finite difference [18,20,27,28,30] and finite element [1,5, 7,14,19] methods. With the easy implementation and less storage, finite difference methods are preferred in certain engineering fields such as geophysics. As is well known, all numerical methods require a finer mesh in order to ensure the accuracy with increasing wavenumber, which is called the “pollution effect”. During the past few decades, many finite difference methods [6,18,20,23,26] are developed to reduce the numerical dispersion caused by the pollution effect. In [20], rotated 9-point finite difference scheme was proposed to minimize the numerical dispersion, combining two discretizations on the classical Cartesian coordinated system and the 45◦ rotated system. In [26], the rotated 9-point difference scheme was extended to the 25-point scheme, which gave a better property. Though the extended 25-point scheme is a popular scheme for the Helmholtz equation, it is not a good choice for the Helmholtz-PML equation since it is not consistent. In [6], we developed an optimal 25-point finite difference scheme, which is pointwise consistent with the Helmholtz-PML equation. With a refined parameter choice strategy, the optimal 25-point scheme not only improves the accuracy, but also reduces the numerical dispersion more significantly. To reduce the numerical error, high-order schemes [27,29] have also been constructed, with the source term required to be smooth enough. After discretization, the resulting linear system from the Helmholtz equation may be extremely large, especially for large wavenumbers. We cannot afford to solve it with a direct method, and have to resort to iterative methods. However, the resulting linear system is non-Hermitian and indefinite, which leads to a very slow convergence or even divergence of the iterative method. Hence, current iterative solvers are not competitive without preconditioning, and preconditioned iterative methods arouse increasing interest. Recently, in [11], an efficient preconditioned Krylov subspace method was proposed for the Helmholtz equation. The preconditioner was based on the complex shifted Laplacian, which was a generalization of the shifted-Laplacian [2,12,15,22]. This kind of preconditioned Krylov subspace method attracted a lot of interests recently [4,33,17,21,24,31]. A spectral analysis was made in [4,33] for the system preconditioned with the complex shifted-Laplacian preconditioner. In [21] and [24], the preconditioner was generalized for three-dimensional and parallel case, respectively. In [17] and [31], the higher order finite element and finite difference methods were utilized to test the performance of the complex shifted-Laplacian preconditioner. To invert the preconditioner approximately, a multigrid method is employed, aiming to gain a more approximate inverse at a lower cost. The prolongation operator is an important component of the multigrid. For Helmholtz problems with constant wavenumber in homogeneous medium, the multigrid method with traditional linear prolongation operator performs well. However, for varying wavenumber in heterogeneous medium, the linear prolongation operator shall lead to a divergent multigrid, and hence a divergent preconditioned Krylov subspace method. To improve the performance of the multigrid, in [9,11,24], matrix-based prolongation operators were employed, with the prolongation operator originally proposed by de Zeeuw in [8]. Based on the prolongation operator of de Zeeuw’ type, we developed new matrix-based multigrid method in [4]. However, all these matrix-based prolongation operators are only restricted to 9-point finite difference stencils. That is, they are only suitable for Helmholtz problems with 9-point finite difference scheme. For the 25-point difference scheme, the corresponding prolongation operator is not available, and its construction is not so easy as that for the 9-point scheme. In [31], a prolongation operator is constructed based on algebraic multigrid (AMG) principles. In the case of 25-point difference scheme, the corresponding prolongation operator can also be constructed in this way, however, it is not efficient enough. In this paper, we consider a multigrid-based preconditioned Krylov subspace method for solving the HelmholtzPML equation with a discretization of 25-point finite difference scheme. We generalize the complex shifted-Laplacian to the complex shifted-Laplacian-PML for preconditioning the Helmholtz-PML equation. To suppress the pollution effect efficiently, the optimal 25-point finite difference scheme we developed in [6] are employed to discretize the preconditioned system. We do spectral analysis for the discrete preconditioned system from the perspective of linear fractal mapping in complex plane, which gives a clear spectral behavior. We shall present the conditions under which the spectrum of the preconditioned system will be favorable for the convergence of the Krylov subspace method. Then, the Krylov subspace method Bi-CGSTAB [32] is employed to solve the discrete preconditioned system, with the matrix-based multigrid used to invert approximately the preconditioner. For the multigrid cycle, we specially construct a matrix-based prolongation operator, according to the 25-point finite difference stencil. Numerical experiments are presented to illustrate the efficiency of the multigrid-based preconditioned Bi-CGSTAB method with the new prolongation operator. The performance of the new prolongation operator is also compared with that of the prolongation originated from algebraic multigrid (AMG) [31]. The experiments range from constant wavenumber in homogeneous medium to greatly varying one in heterogeneous medium.

56

D. Cheng et al. / Mathematics and Computers in Simulation 117 (2015) 54–67

This paper is organized as follows. In Section 2, we generalize the complex shifted-Laplacian preconditioner to the complex shifted-Laplacian-PML preconditioner, and give some spectral analysis for the discrete preconditioned system. In Section 3, we develop a new prolongation operator for the matrix-based multigrid, which is used to invert approximately the preconditioner in the preconditioned Bi-CGSTAB method. In Section 4, numerical experiments are presented to confirm the efficiency of the multigrid-based preconditioned Bi-CGSTAB method with the new prolongation operator. Finally, in Section 5, some conclusions are drawn. 2. Spectral analysis for the discrete Helmholtz-PML equation preconditioned with complex shifted-LaplacianPML In this section, we present the complex shifted-Laplacian-PML for the Helmholtz-PML equation, and further study the spectral distribution of the discrete preconditioned system, from the perspective of linear fractal mapping. 2.1. The complex shifted-Laplacian-PML preconditioner We consider the 2D Helmholtz equation − ∆u − (1 − αi)k 2 u = g

in R2 ,

(2.1)

where ∆ is the Laplacian, unknown u usually √ represents a pressure field, α is the nonnegative real number indicating the fraction of damping in the medium, i = −1 is the imaginary unit, k := 2π f /v is the wavenumber with f , v indicating the frequency and wave speed respectively, finally the right side g represents the source term. The wavenumber k, depending on speed v, is a constant in the homogeneous medium, and varies in the heterogeneous medium. For α = 0, (2.1) is the common Helmholtz equation, while for α > 0 it is a damped one, with the maximum value of α being 0.05 (5%) in geophysical applications. To solve (2.1) in a finite domain, artificial boundary conditions should be employed. Here, we consider the boundary condition of perfect match layer (PML), which has an astonishing property of generating almost no reflection. With the PML, (2.1) can be rewritten as     ∂ e y ∂u ∂ ex ∂u Au := − − − (1 − αi)ex e y k 2 u = g, (2.2) ∂ x ex ∂ x ∂ y ey ∂ y where the 1D functions ex (x) and e y (y) are damping functions implemented in the PML, and they equal to one in the interior area. The thickness of PML is usually chosen to be a half-wavelength which is defined by ι := v/ f . The number of wavelengths in a square domain of size H is H/ι, and the dimensionless wavenumber [19] is 2π f H/v. The wavenumber refers to dimensionless wavenumber (also denoted by k) in the remainder of this text. To make iterative methods competitive for solving the linear system stemming from the Helmholtz equation, it is often preconditioned with an efficient preconditioner. In [11], the complex shifted-Laplacian was proposed as a preconditioner, which is a generalization of shifted-Laplacian preconditioners [2,12,15]. For the Helmholtz-PML equation, the preconditioner we consider is based on the operator     ∂ e y ∂u ∂ ex ∂u M := − − − (1 − βi)ex e y k 2 , (2.3) ∂ x ex ∂ x ∂ y ey ∂ y with β > α ≥ 0. We call the operator (2.3) as the complex shifted-Laplacian-PML, which is reduced to the complex shifted-Laplacian when ex = e y ≡ 1. The complex shifted-Laplacian-PML preconditioner contributes to a nice spectral distribution of the discrete preconditioned system, which shall be presented later. To precondition Eq. (2.2) with the preconditioner (2.3), we have the preconditioned system in the operator form AM−1 v = g,

Mu = v.

(2.4)

We now discretize the operator equation (2.4) by using the optimal 25-point finite difference scheme we proposed in [6]. For the Helmholtz equation, some efficient difference schemes [18,20,23,26] are proposed to reduce the numerical dispersion. For 2D Helmholtz problems, the extended 25-point difference scheme in [26] is a popular method, which has a better performance than the rotated 9-point scheme [20]. However, the extended 25-point scheme is

D. Cheng et al. / Mathematics and Computers in Simulation 117 (2015) 54–67

57

not consistent with Helmholtz-PML equation, due to the PML. In [6], we developed the optimal 25-point difference scheme, which is consistent with the Helmholtz-PML equation. Based on minimizing the numerical dispersion, we then proposed the refined strategy for choosing optimal parameters. The optimal 25-point scheme gives a good performance with respect to the wavenumber k, and numerical results show that the pollution effect of large wavenumbers can be suppressed to a certain degree. After discretization, we obtain the discrete preconditioned linear system AM−1 v = g,

(2.5)

where A, M ∈ C N ×N , u, g ∈ C N , C is the complex number set and N is the number of unknowns. Matrices A and M are the discrete Helmholtz operator and preconditioner respectively, and they have the form A := L − z 1 D,

z 1 := 1 − αi

M := L − z 2 D,

z 2 := 1 − βi

(2.6)

and

where matrix L is the discretization of high order term − ∂∂x

(2.7) 

e y ∂u ex ∂ x



−

∂ ∂y



ex ∂u ey ∂ y



and matrix D is the discretization

of zeroth order term −(1 − βi)ex e y Moreover, A, M have complex entries, and are sparse, nonsymmetric, and diagonal-distributional. It is seen that M differs from A only in the damping factor, that is, β > α. The choice of β is a compromise between the convergence rate of the multigrid method and a favorable spectrum of the preconditioned system (2.5) (cf. [4]). k 2 u.

2.2. Spectral analysis of the discrete preconditioned system We now give a spectral analysis of the discrete preconditioned system (2.5) from the perspective of linear fractal mapping in complex variable function. For the shifted-Laplacian preconditioner, [33] analyzed the spectrum of the discrete preconditioned Helmholtz equation, with the homogeneous Neumann condition, homogeneous Dirichlet condition and Sommerfeld condition. In [4], we presented a spectral analysis from the perspective of the linear fractional mapping, which makes the spectral behavior easy to understand. Here, we continue the analysis method and give a further study for the spectral distribution of the discrete preconditioned system. For any matrix A, B ∈ C N ×N , we denote the spectrum of A by σ (A), and the spectrum of A with respect to B by σ (A, B). When B = I is the identity matrix, σ (A) = σ (A, I). It is clear that σ (AB−1 ) = σ (B−1 A) = σ (A, B) if B is nonsingular. To analyze the spectrum of (2.5), we first recall a lemma in [4]. Lemma 2.1. Let L, D ∈ C N ×N , A := L − z 1 D and M := L − z 2 D with z 1 , z 2 ∈ C. If D and M are nonsingular, −1 −1 1 µ ∈ C and λ := µ−z µ−z 2 with µ ̸= z 2 then λ ∈ σ (A, M) if and only if µ ∈ σ (L, D), and M A, D L share the same eigenvectors. Similar results have also been described in [33], and a different presentation in Lemma 2.1 helps to develop a linear fractional mapping. Denote the real and imaginary parts of µ by Re(µ) and Im(µ) respectively, and define λ = λ(µ) :=

µ − z1 , µ − z2

(2.8)

which is a linear fractional mapping. Then, we have the following proposition, which is a generalization of Lemma 2.2 in [4]. Proposition 2.2. Let the linear fractional mapping λ : C → C be defined by (2.8). Then the straight ¯  line γ µ+  γ¯ µ+d  γ (z 1 −z 2 )  γ z 1 +γ¯ z¯ 2 +d = 0 with γ ∈ C and d ∈ R is mapped to the circle O(c, R) with c := γ z 2 +γ¯ z¯ 2 +d and R :=  γ z 2 +γ¯ z¯ 2 +d , and the half-plane γ µ + γ¯ µ¯ + d > 0 and γ µ + γ¯ µ¯ + d < 0 are mapped inside and outside O(c, R), respectively. Proof. In the complex plane C = {µ : µ = x + yi, x, y ∈ R}, a circle can be represented by a(x 2 + y 2 ) + bx + cy + d = 0

(2.9)

58

D. Cheng et al. / Mathematics and Computers in Simulation 117 (2015) 54–67

a

b

Fig. 1. The linear fractional mapping with lines mapped into circles. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) b c with parameters a, b, c, d ∈ R satisfying a ≥ 0 and b2 + c2 > 4ad. The center of this circle is ( 2a , 2a ) and the radius √ 2 2 b +c −4ad is R := . Noting that 2a

µ + µ¯ µ − µ¯ , y= , 2 2i (2.9) can be rewritten as x=

x 2 + y 2 = µµ, ¯

aµµ¯ + γ µ¯ + γ¯ µ + d = 0 with γ := 21 (b + ci), |γ |2 > ad. The center and the radius can be represented by c = − γa and R = respectively. When a = 0, (2.10) is reduced to a straight line γ µ + γ¯ µ¯ + d = 0,

(2.10) √ |γ |2 −ad a

(2.11)

which can be considered as a circle with R = ∞ on the complex plane. From (2.8), we have z2λ − z1 z¯ 2 λ¯ − z¯ 1 , µ¯ = . λ−1 λ¯ − 1 Substitution of (2.12) into (2.11) yields µ=

(2.12)

a ′ λλ¯ + γ ′ λ¯ + γ¯ ′ λ + d ′ = 0,

(2.13)

where a ′ := γ z 2 + γ¯ z¯ 2 + d, γ ′ := −(γ z 1 + γ¯ z¯ 2 + d) and d ′ := (γ z 1 + γ¯ z¯ 1 + d).√Thus, (2.13) represents a circle,   ′  γ (z 1 −z 2 )  |γ ′ |2 −a ′ d ′ γ¯ z¯ 2 +d and the radius R := = denoted by O(c, R), with the center c := − γa ′ = γγ zz 12 +  ′ +γ¯ z¯ 2 +d γ z 2 +γ¯ z¯ 2 +d . It can a be easily obtained that γ µ + γ¯ µ¯ + d > 0 is mapped inside O(c, R), and γ µ + γ¯ µ¯ + d < 0 is mapped outside O(c, R).  In the remainder of the text, we let z 2 = 1 − βi (β > α ≥ 0). Fig. 1 presents the linear fractional mapping (2.8) for z 1 = 1 (α = 0) and z 2 = 1 − 0.5i (β = 0.5). According to Proposition 2.2, we can obtain that the line l1 : µ − µ¯ = 0 (real axis, plotted in red in Fig. 1(a)) is mapped into the circle l1′ (plotted in red in Fig. 1(b)) with the center c = 0.5 and radius R = 0.5. The interval [0.8, 1.2] (dotted line in Fig. 1(a)) is mapped into a part of l1′ (dotted line in Fig. 1(a)) near the imaginary axis. The upper half-plane (the shadow region, denoted by Ω + ) and the lower half-plane (Ω − ) are mapped inside and outside of the circle l1′ respectively. Moreover, the points µ0 = 1 + 0.5i, µ1 = ∞, µ2 = 0.5, µ3 = 1, µ4 = 1.5 in Fig. 2(a) are mapped to λ0 = 0.5 (the center of l1′ ), λ1 = 1, λ2 = 0.5 + 0.5i, λ3 = 0, λ4 = 0.5 − 0.5i in Fig. 2(b), respectively. We also can see that the line segment, which is a part of l2 : γ µ + γ¯ µ¯ = 0 (γ = 3 + 10i, plotted in blue in Fig. 1(a)), is mapped into a part of the circle l1′ with c = 0.6875 + 0.0938i and radius R = 0.3263 (plotted in blue in Fig. 1(b)). According to above results, we can obtain that if σ (L, D) has a nonnegative imaginary part, then σ (A, M) would be enclosed by certain circle on the right half-plane, which is favorable for the convergence of Krylov subspace, such as GMRES [25] and Bi-CGSTAB method [32]. With the optimal 25-point difference scheme and the PML, matrix L is

D. Cheng et al. / Mathematics and Computers in Simulation 117 (2015) 54–67

a

b

59

c

Fig. 2. The spectral distribution for matrices (a) A, (b) D−1 L, (c) AM−1 . (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

neither Hermitian nor symmetric, hence we cannot prove theoretically that σ (L, D) has a nonnegative imaginary part. However, for different wavenumber and discretization step, numerical results show that there holds the sufficient condition. Fig. 2 presents the spectral distribution of A, D−1 L and AM−1 respectively, for z 1 = 1, z 2 = 1 − 0.5i, k = 50, 1 and h = 40 . Fig. 2(a) is the original spectrum for A, located on both the left and right half-planes, which means that A is indefinite. Fig. 2(b) shows the spectrum for σ (L, D), which has a nonnegative imaginary part. According to above results, we can expect a clustered spectrum which is enclosed by the circle O(0.5, 0.5), as presented in Fig. 2(c). As can be seen, the pink part (nearly a straight line) in Fig. 2(b) is mapped into the pink counterpart (nearly a circle) in Fig. 2(c). As is known, it is very favorable for the convergence of Krylov subspace methods when the spectrum of the preconditioned system is clustered enough in the right half-plane, which gives a very small condition number. However, it should be pointed out that Krylov subspace methods are still efficient for the system with its spectrum far away from the origin in the right half-plane [33], even though it is not so clustered as being enclosed by certain circle. The following proposition indicates under what conditions the spectrum σ (A, M) shall be far away from the origin in the right half-plane. Proposition 2.3. Let the linear fractional mapping λ : C → C be defined by (2.8) with z 1 := 1 − αi, z 2 := 1 − βi β−α ¯ and β > α ≥ 0. Then the circle O(c, R) with c = (1, 2βd−α−β 2(1−l) ) and R = 2(1−q) is mapped to the line λ + λ = 2q (0 ≤ q < 1). The outside of O(c, R) is mapped to the half-plane λ + λ¯ > 2q, and the inside of O(c, R) is mapped to the half-plane λ + λ¯ < 2q. Proof. According to the proof of Proposition 2.2, on the complex plane, a circle can be represented by aµµ¯ + γ µ¯ + γ¯ µ + d = 0

(2.14)

with the center and the radius can be represented by c = − γa and R = β−α with c = (1, 2βq−α−β 2(1−q) ) and R = 2(1−q) , we have a = 1, γ = −1 − substituting (2.12) into (2.14) with above a, γ , d yields

√ |γ |2 −ad respectively. For the circle O(c, R) a 2βq−α−β β(βq−α) 2(1−q) i, and d = 1 − 1−q . Similarly,

λ + λ¯ = 2q.

(2.15)

Moreover, the outside of O(c, R) is mapped to half-plane on the right of line λ + λ¯ = 2q, that is, λ + λ¯ > 2q. Correspondingly, the inside of O(c, R) is mapped to the half-plane λ + λ¯ < 2q.  According to Proposition 2.2, for z 1 := 1 − αi and z 2 := 1 − βi, (2.8) maps the real axis (i.e., µ − µ¯ = 0) β−α into the circle O(c, R) with c := β+α 2β , R := 2β . Denote the distance between the origin and O(c, R) by ϵ := infµ∈O(c,R) |µ|, then we have ϵ = βα . Let q = ϵ = βα in Proposition 2.3, we have the following proposition. Proposition 2.4. Let the linear fractional mapping λ : C → C be defined by (2.8) with z 1 := 1 − αi and z 2 := 1 − βi β > α ≥ 0. The circle O(c, R) with c = (1, − β2 ) and R = β2 is mapped to the line λ + λ¯ = 2α β . The outside and 2α 2α ¯ ¯ inside of O(c, R) are mapped to the half-planes λ + λ > and λ + λ < , respectively. β

β

60

D. Cheng et al. / Mathematics and Computers in Simulation 117 (2015) 54–67

a

b

Fig. 3. The linear fractal mapping with the circle mapped into the line. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

a

b

c

Fig. 4. (a) The spectrum for D−1 L; (b) a snapshot of (a) and the circle O(c, R) with c = (1, −0.25) and R = 0.25; (c) the spectrum for AM−1 with line λ + λ¯ = 0.2 plotted. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 3 shows the linear fractal mapping with z 1 = 1 − 0.05i (α = 0.05), z 2 = 1 − 0.5i (β = 0.5). In Fig. 3(a), the real axis l1 (plotted in red) is mapped to the circle l1′ : O(0.55, 0.45) in Fig. 3(b). The circle l2 : O(c, R) with c = (1, −0.25), R = 0.25 (plotted in blue) is mapped into the line l2′ : λ + λ¯ = 2q with q = βα = 0.1. The points µ0 = 1 − 0.25i, µ1 = 1, µ2 = 0.75 − 0.25i, µ3 = 1 − 0.5i, µ4 = 1.25 − 0.25i are mapped to λ0 = −0.8, λ1 = 0.1, λ2 = 0.1 + 0.9i, λ3 = ∞, λ4 = 0.1 − 0.9i, respectively. In addition, the outside (denoted by shadow region) and inside of l2 are mapped to the half-plane λ + λ¯ > 0.2 and λ + λ¯ < 0.2 respectively. According to Proposition 2.4, if σ (L, D) is far from the small circle O(c, R) with c = (1, β2 ) and R = β2 , then σ (A, M) would be on the right of the vertical line λ + λ¯ = 2α β and far from it. We can see this in Fig. 4 with z 1 =

1 1 − 0.05i, z 2 = 1 − 0.5i, k = 50, and h = 40 . Fig. 4(a) and (c) show the spectrum for D−1 L and AM−1 . Fig. 4(b) presents a snapshot of σ (L, D) and the circle O(c, R) with c = (1, −0.25) and R = 0.25. We can observe that σ (L, D) is located far from the circle O(c, R) (plotted in red). With the linear fractal mapping (2.8), σ (L, D) and O(c, R) in Fig. 4(b) are mapped to σ (A, M) and λ + λ¯ − 0.2 = 0 respectively in Fig. 4(c). The outside (shadow region) and inside of the circle are mapped to λ + λ¯ > 0.2 (shadow region) and λ + λ¯ < 0.2, respectively. Since σ (L, D) is located outside of the circle, σ (A, M) lies in the half-plane λ + λ¯ > 0.2. We can also see that the center (1, −β) and radius β2 of the circle are only with respect to β, and independence of α. When there is no damping, namely, α = 0, the image of circle O(c, R) is reduced to the imaginary axis. With the optimal 25-point finite difference scheme, for different k and h, we can obtain numerically that σ (L, D) is located outside of the circle O(c, R), hence σ (A, M) lies in the right-half plane, which favors the convergence of the Krylov subspace method. Finally, we point out that the choice of β is a compromise between the convergence rate of the multigrid method and a favorable spectrum of the preconditioned system. In this text, we choose β = 0.5, which is a proper choice in practice. If β is too small, a good spectrum of AM−1 can be obtained, but it is difficult to invert approximately the

D. Cheng et al. / Mathematics and Computers in Simulation 117 (2015) 54–67

a

b

61

c

1 , and (a) β = 0.1; (b) β = 0.5; (c) β = 1.0. Fig. 5. The spectrum AM−1 for k = 50, h = 40

preconditioner M with the multigrid method, since M has a poor property similar to A. However, if β is too large, we would obtain an unfavorable spectrum near the origin, which leads to a poor performance of the preconditioning. Fig. 5 presents the spectral distribution of AM−1 for β = 0.1, β = 0.5 and β = 1.0, respectively. Both the cases in Fig. 5(a) and (c) would weaken the efficiency of the Krylov subspace method. Factually, in the pioneering work [11] for the complex shifted-Laplacian preconditioner, β = 0.5 was suggested since it contributed to a fast convergence of the preconditioned Krylov subspace method. 3. A new prolongation operator for the multigrid-based preconditioned Bi-CGSTAB method In this section, we propose a new matrix-based prolongation operator for the multigrid cycle, which is used to approximate the inverse of the preconditioner M in the Bi-CGSTAB method for solving the discrete preconditioned system (2.5). By the spectral analysis in Section 2, we obtain that the spectrum of the discrete preconditioned system is favorable for the convergence of the Krylov subspace method. We choose the Bi-CGSTAB method [32] to solve the discrete preconditioned system (2.5), which is called as the preconditioned Bi-CGSTAB method. Then, we have to invert approximately the preconditioner M in the form of u := M−1 v, where v is the given vector with u to be determined. To avoid inverting M directly, we usually choose to solve the additional linear systems Mu = v. In order to obtain a good approximate solution at a low cost, the multigrid method is employed to solve the additional linear system. For the multigrid method, the prolongation operator is an important component, and has much effect on the convergence. The traditional linear prolongation operator is commonly used, however, it does not work well for solving Mu = f with varying wavenumbers, and leads to a poor performance of the preconditioned Bi-CGSTAB method. To remedy this problem, matrix-based prolongation operators have been designed, which give a good performance (see [4,9,11,24]). The matrix-based prolongation operator is originated from the one proposed by de Zeeuw in [8], and is constructed based on the finite difference stencil, namely, the algebraic information of the linear system. However, previous matrixbased prolongation operators are all based on 9-point finite difference schemes. For the optimal 25-point difference scheme employed in this paper, we have to design a new matrix-based prolongation operator accordingly. Next, we shall present the construction of the new prolongation operator, based on [8] and the property of the optimal 25-point difference scheme. For a clear description, we first begin with two graphs. Fig. 6(a) presents the numbering in optimal 25-point finite difference stencil, with m j,l ( j, l = −2, −1, 0, 1, 2) representing the component of the stencil. Fig. 6(b) shows one coarse and four fine grid cells with capital letters indicating coarse grid points (plotted in red) and small letters indicating fine grid points (plotted in blue). The fine grid denoted by Ωh is split into four disjunct sub-grid Ωh = {Ωh,(0,0) , Ωh,(1,0) , Ωh,(0,1) , Ωh,(1,1) }. Here, Ωh,(0,0) consists of the fine grid point which is also the coarse grid point. Ωh,(1,0) and Ωh,(0,1) consist of the fine grid point located between two coarse grid points along the x axis and y axis, respectively. Ωh,(1,1) consists of the fine grid point which is not next to any coarse grid point. For example in Fig. 6, A, B, C, D ∈ Ωh,(0,0) , q, q ′ ∈ Ωh,(1,0) , p, p ′ ∈ Ωh,(0,1) , and t ∈ Ωh,(1,1) . The symbols eh , e H represent the grid functions (corrections) for fine and coarse

62

D. Cheng et al. / Mathematics and Computers in Simulation 117 (2015) 54–67

a

b

Fig. 6. (a) 25-point stencil with numbering. (b) Coarse and fine grid points, which are designated by capital and small letters respectively. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

a

b

c

d

Fig. 7. Different 9-point stencils employed in the construction of the optimal 25-point finite difference scheme.

grid respectively. If e H is given, eh can be obtained from the prolongation of e H by  eh (A) := e2h (A), A ∈ Ωh,(0,0) ,      eh ( p) := W A ( p)e2h (A) + W B ( p)e2h (B), p ∈ Ωh,(0,1) , eh (q) := W A (q)e2h (A) + W D (q)e2h (D), q ∈ Ωh,(1,0) ,    eh (t) := W A (t)e2h (A) + W B (t)e2h (B)    + WC (t)e2h (C) + W D (t)e2h (D), t ∈ Ωh,(1,1) ,

(3.1)

where eh (·), e H (·) denote the components of corrections on fine and coarse grids respectively, and W A (·), W B (·), WC (·), W D (·) represent prolongation weights to be determined. When prolongation weights are determined on all the fine grid points, the prolongation operator (matrix), denoted by P, is obtained, and eh = Pe H . Now, we shall describe the construction of a new matrix-based prolongation operator P for the optimal 25-point finite difference scheme. Since M is nonsymmetric, we split it into a symmetric and an anti-symmetric part. That is, let M = Ms + Mu with Ms :=

1 (M + MT ), 2

Mu :=

1 (M − MT ). 2

Let the 25-point stencil of M be m

−2,−2

m −2,−1  M :=  m −2,0 m −2,1 m −2,2

m −1,−2 m −1,−1 m −1,0 m −1,1 m −1,2

m 0,−2 m 0,−1 m 0,0 m 0,1 m 0,2

m 1,−2 m 1,−1 m 1,0 m 1,1 m 1,2

m 2,−2  m 2,−1   m 2,0  . m 2,1  m 2,2

(3.2)

Then, we can get corresponding stencils of Ms and Mu , which are denoted by Ms and Mu with m s, j,l , m u, j,l ( j, l = −2, −1, 0, 1, 2) representing their components, respectively. For the construction of the optimal 25-point finite difference scheme in [6], it employs four 9-point stencils, which are presented in Fig. 7. For the Laplacian-PML, Fig. 7(a), (b), (d) are used to discretize it along the x axis, and Fig. 7(a), (c), (d) are used to discretize it along the y axis. To determine prolongation weights on p ∈ Ωh,(0,1) and

D. Cheng et al. / Mathematics and Computers in Simulation 117 (2015) 54–67

63

q ∈ Ωh,(1,0) , with the stencil Ms of the symmetric part, we define dn1 := max(|m s,−1,−1 + m s,0,−1 + m s,1,−1 |, |m s,−1,−1 |, |m s,1,−1 |), dn2 := max(|m s,−1,−2 + m s,0,−2 + m s,1,−2 |, |m s,−1,−2 |, |m s,1,−2 |), dn3 := max(|m s,−2,−2 + m s,0,−2 + m s,2,−2 |, |m s,−2,−2 |, |m s,2,−2 |), ds1 := max(|m s,−1,1 + m s,0,−1 + m s,1,1 |, |m s,−1,1 |, |m s,1,1 |), ds2 := max(|m s,−1,2 + m s,0,2 + m s,1,2 |, |m s,−1,2 |, |m s,1,2 |), ds3 := max(|m s,−2,2 + m s,0,2 + m s,2,2 |, |m s,−2,2 |, |m s,2,2 |), and dw1 := max(|m s,−1,−1 + m s,−1,0 + m s,−1,1 |, |m s,−1,−1 |, |m s,−1,1 |), dw2 := max(|m s,−2,−1 + m s,−2,0 + m s,−2,1 |, |m s,−2,−1 |, |m s,−2,1 |), dw3 := max(|m s,−2,−2 + m s,−2,0 + m s,−2,2 |, |m s,−2,0 |, |m s,−2,2 |), de1 := max(|m s,1,−1 + m s,1,0 + m s,1,1 |, |m s,1,−1 |, |m s,1,1 |), de2 := max(|m s,2,−1 + m s,2,0 + m s,2,1 |, |m s,2,−1 |, |m s,2,1 |), de3 := max(|m s,2,−2 + m s,2,0 + m s,2,2 |, |m s,2,0 |, |m s,2,2 |). Let dn := max(dn1 , dn2 , dn3 ),

ds := max(ds1 , ds2 , ds3 ),

dw := max(dw1 , dw2 , dw3 ),

de := max(de1 , de2 , de3 ),

and 2       σ := min 1, 1 − m s, j,l /m s,0,0  . j,l=−2

For the stencil Mu of the anti-symmetric part, we define     (m u, j,−1 − m u, j,1 ) + (m u, j,−2 − m u, j,2 ) + c1 :=  j=−1,0,1

   c2 :=  (m u,−1,l − m u,1,l ) + l=−1,0,1

j=−1,0,1



(m u,−2,l − m u,2,l ) +

l=−1,0,1

Then, prolongation weights on p and q can be derived by  1 1 d − d  1 c1 n s   W A ( p) = σ  + +  , 2 2 dn + ds 2 dw + de + dn + ds

 j=−2,0,2



  (m u, j,−2 − m u, j,2 ).

  (m u,−2,l − m u,2,l ).

l=−2,0,2

W D ( p) = σ (1 − W A ( p)),

(3.3)

W B (q) = σ (1 − W A (q)).

(3.4)

and  1 1 d − d  1 c2 w e   W A (q) = σ  + +  , 2 2 dw + de 2 dw + de + dn + ds

We can similarly determine prolongation weights on p ′ ∈ Ωh,(0,1) and q ′ ∈ Ωh,(1,0) . With weights obtained above, we now determine prolongation weights on the gridpoint t ∈ Ωh,(1,1) . Applying the multigrid to solve Mu = v, after ˜ with the residual being r˜ := v − Mu. ˜ After several relaxation sweeps, we get an approximate solution, denoted by u, coarse grid correction, we get a correction solution uˆ := u˜ + eh = u˜ + Pe H , and the residual is rˆ := v − Muˆ = r˜ − MPe H . As in [8], in order to prevent huge jumps in the norm of the residual after prolongation, we require MPe H |t = 0,

64

D. Cheng et al. / Mathematics and Computers in Simulation 117 (2015) 54–67

which is equivalent to 2 2  

m j,l (t)eh,( j,l) (t) = 0.

(3.5)

j=−2 l=−2

Here, if ( j, l) = (0, 0), eh,( j,l) (t) denotes the component of eh restricted on t, otherwise it denotes the component restricted on the neighbor of t. We rewrite (3.5) as 1 1  

2  

m j,l (t)eh,( j,l) (t) +

j=−1 l=−1

1  

m j,l (t)eh,( j,l) (t) +

j=−2,2 l=−2

m j,l (t)eh,( j,l) (t) = 0.

(3.6)

j=−1 l=−2,2

Substituting  eh,(−1,0) (t) := W A ( p)e2h (A) + W B ( p)e2h (B),    eh,(1,0) (t) := WC ( p ′ )e2h (C) + W D ( p ′ )e2h (D), eh,(0,−1) (t) := W A (q)e2h (A) + W D (q)e2h (D),    eh,(0,1) (t) := W B (q ′ )e2h (B) + WC (q ′ )e2h (C),

(3.7)

into the first part of (3.6) and neglecting remaining parts, we can determine weights W A (t), W B (t), WC (t) and W D (t) in eh,(0,0) (t) := W A (t)e2h (A) + W B (t)e2h (B) + WC (t)e2h (C) + W D (t)e2h (D). However, the second part of (3.6) should not be neglected in order to get a good prolongation at r . Then, we assume that eh,( j,−2) (t) ≈ eh,( j,−1) (t), eh,(−2,l) (t) ≈ eh,(−1,l) (t),

eh,( j,2) (t) ≈ eh,( j,1) (t), j = −2, −1, 0, 1, 2; eh,(2,l) (t) ≈ eh,(1,l) (t), l = −2, −1, 0, 1, 2;

(3.8)

which is a proper approximation. Substituting (3.7) and (3.8) into (3.6) and choosing for equality, we then have 2 2  

m ′j,l (t)eh,( j,l) (t) = 0,

(3.9)

j=−2 l=−2

which is an approximation of (3.5) with m ′h,(−1,−1) (t) =

−1  −1 

m h,( j,l) (t),

m ′h,(1,−1) (t) =

j=−2 l=−2

m ′h,(−1,1) (t) =

−1  2 

m h,( j,l) (t),

±2  l=±1

m h,(0,l) (t),

m h,( j,l) (t),

j=1 l=−2

m ′h,(1,1) (t) =

j=−2 l=1

m ′h,(0,±1) (t) =

2  −1 

2  2 

m h,( j,l) (t),

j=1 l=1

m ′h,(±1,0) (t) =

±2 

m h,( j,0) (t),

m ′h,(0,0) (t) = m h,(0,0) (t).

j=±1

From (3.9), we can compute weights W A (t), W B (t), WC (t) and W D (t) as follows:   W A (t) := − m ′h,(−1,−1) (t) + W A ( p)m ′h,(−1,0) (t) + W A (q)m ′h,(0,−1) (t) /m ′h,(0,0) (t),   W B (t) := − m ′h,(−1,1) (t) + W B ( p)m ′h,(−1,0) (t) + W B (q ′ )m ′h,(0,1) (t) /m ′h,(0,0) (t),   WC (t) := − m ′h,(1,1) (t) + WC ( p ′ )m ′h,(1,0) (t) + WC (q ′ )m ′h,(0,1) (t) /m ′h,(0,0) (t),   W D (t) := − m ′h,(1,−1) (t) + W D ( p ′ )m ′h,(1,0) (t) + W D (q)m ′h,(0,−1) (t) /m ′h,(0,0) (t). We remark that the construction of the prolongation weights for gridpoints in Ωh,(0,1) and Ωh,(1,0) is based on the symmetric and anti-symmetric part of M. To obtain good prolongation weights for gridpoints in Ωh,(1,1) , a new 9-point

D. Cheng et al. / Mathematics and Computers in Simulation 117 (2015) 54–67

a

65

b

Fig. 8. (a) The real part of numerical solution at k = 50 without damping. (b) The real part of numerical solution at k = 100 without damping. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

stencil is used to approximate the original 25-point stencil. The resulting weights would be efficient in the multigridbased preconditioned Bi-CGSTAB method, which shall be presented later in next section. When the Helmholtz-PML with 25-point stencil is reduced to the Laplacian with classical 5-point stencil, the prolongation operator above is just the bilinear prolongation operator, which contributes to an efficient multigrid method. 4. Numerical experiments In this section, numerical experiments related to geophysical applications are presented to show the efficiency and robustness of the multigrid-based preconditioned Bi-CGSTAB method for the discrete preconditioned system (2.5), which is discretized from the Helmholtz-PML equation preconditioned with a complex shifted Laplacian-PML by an optimal 25-point difference scheme. We chose β1 = 1 and β2 = 0.5 in (2.5) as basic parameters for computation, which is a proper choice. To approximate the inverse of the preconditioner M, the full multigrid V-cycle (FMG) is chosen since it possesses both the robustness of the W-cycle and the efficiency of the V-cycle, and only one iteration with two relaxation sweeps is enough. For multigrid components, the pointwise Jacobi relaxation with underrelaxation (ω-JAC) is used as a smoother, and weighting restriction operator and the newly proposed prolongation operator are adopted as intergrid transfer operators. All the experiments here are performed serially on an Intel Xeon (8-core) with 3.33 GHz and 96 Gb RAM, with Matlab and Visual C++. In the computation, a zero initial guess has been chosen, and the preconditioned Bi-CGSTAB iteration terminates when the Euclidean norm of the relative residual error is reduced to the order of 10−6 . The experiments range from constant wavenumber in homogeneous medium to greatly varying wavenumber in heterogeneous medium. For the constant wavenumber, we consider a computational domain Ω = (0, 1) × (0, 1), with a point source is located at the center of the domain. The number of gridpoints per wavelength is chosen to be G = 8, and an π accuracy requirement for second order discretizations is that kh ≤ 2π G = 4 . The thickness of the PML is set to be 20, that is, the PML possesses 20 gridpoints in each direction. Fig. 8(a) and (b) present the numerical solution corresponding to k = 60 and k = 120, respectively. As can be seen, there are no physical reflections at the boundaries due to the PML. Table 1 gives the number of preconditioned Bi-CGSTAB iterations and the CPU time needed (in parentheses) in seconds, with and without damping for different wavenumbers respectively. The largest wavenumber for test is k = 800, and the number of unknowns (without PML) is 10242 = 1040,400 (about 1 million). In Table 1, we also compare the performance of the new prolongation operator with that of the prolongation in [31], where the prolongation operator is chosen based on the algebraic multigrid (AMG) principles. (I) and (II) represent the preconditioned BI-CGSTAB method based on the multigrid with the new prolongation operator in this paper, and the prolongation operator in [31], respectively. As can be seen, the method (I) gives a better performance. It can also be observed that adding some damping in the Helmholtz problem will contribute to a faster convergence of the preconditioned BI-CGSTAB method, which can be expected from the spectral analysis in Section 2. Next, we shall test the multigrid-based preconditioned BI-CGSTAB method for solving the Helmholtz-PML equation with varying wavenumbers in geophysical applications. We consider the salt dome model which mimics the subsurface geology under the sea. The computational domain is defined to be a rectangle of dimension 4540×9660 m2 .

66

D. Cheng et al. / Mathematics and Computers in Simulation 117 (2015) 54–67

b

a

Fig. 9. (a) The salt dome model with velocity profile indicated. (b) The real part of numerical solution at f = 30 Hz. Table 1 Number of Bi-CGSTAB iterations and CPU time (in parentheses) for different wavenumbers k with and without damping. k = 50

k = 100

k = 200

k = 400

k = 800

Grid

65 × 65

129 × 129

256 × 256

511 × 511

1020 × 1020

(I) α = 0.00 α = 0.01 (II) α = 0.00 α = 0.01

15 (0.65) 14 (0.62) 18 (0.77) 16 (0.70)

29 (4.72) 26 (4.21) 35 (5.67) 31 (5.10)

60 (38.7) 45 (29.0) 74 (47.2) 56 (36.3)

126 (307.6) 75 (182.9) 158 (386.2) 95 (232.4)

263 (2402) 98 (913.4) 327 (3018) 124 (1143)

Table 2 Bi-CGSTAB convergence for the salt dome problem with and without damping. Iterations and CPU time in seconds (in parentheses) are shown. f (Hz)

10

15

20

25

30

Grid

227 × 483

284 × 604

378 × 805

454 × 966

568 × 1208

(I) α = 0.00 α = 0.01 (II) α = 0.00 α = 0.01

43 (44.0) 35 (35.7) 49 (49.2) 40 (40.6)

56 (84.1) 40 (60.4) 68 (101.8) 47 (69.7)

81 (226.8) 48 (133.9) 100 (279.3) 56 (157.2)

98 (378.2) 55 (217.6) 120 (475.2) 68 (268.7)

130 (750.4) 70 (405.7) 159 (917.4) 87 (502.5)

A point source is located at point (2270 m, 4830 m) under the surface, with the upper surface assigned to be y = 0. Fig. 9(a) presents the domain, the salt dome, and the variation of speed in the medium. The values for the speed of sound are irregularly structured throughout the domain. The lowest velocity is 1524 m/s and the highest velocity is 4480 m/s in the salt dome (the red part in Fig. 8(a)). The real part of the numerical solution at f = 30 Hz without damping is plotted in Fig. 8(b). Table 2 presents the preconditioned Bi-CGSTAB convergence for the case with and without damping for different frequencies, which vary from 10 Hz to 30 Hz. The numbers of Bi-CGSTAB iterations and CPU time in minutes (in parentheses) are shown. As is observed, method (I) still gains faster convergence speed than method (II). For f = 30 Hz, 568×1208 gridpoints are used, with the total unknowns being 686,144. The number of iterations for methods (I) and (II) without damping are 130 and 159, respectively. As can be seen, for the case with varying wavenumbers, the multigrid-based preconditioned Bi-CGSTAB method with the new prolongation operator is more efficient and robust. 5. Conclusions In this paper, we develop a preconditioned iterative method to solve the Helmholtz-PML equation. The complex shifted-Laplacian is generalized to precondition the Helmholtz-PML equation. An optimal 25-point finite difference scheme is used to discretize the preconditioned system, and a spectral analysis is given for the preconditioned system

D. Cheng et al. / Mathematics and Computers in Simulation 117 (2015) 54–67

67

from the perspective of linear fractional mapping. Then, we adopt the Krylov subspace method Bi-CGSTAB to solve the discrete preconditioned system. To invert the preconditioner approximately, the multigrid method is employed, and a new prolongation operator is constructed for the multigrid cycle. Numerical experiments, including both cases of constant wavenumbers and varying wavenumbers, illustrate the efficiency and robustness of the multigrid-based preconditioned Bi-CGSTAB method with the new prolongation operator. References [1] I. Babuˇska, S. Sauter, Is the pollution effect of the FEM avoidable for the Helmholtz equation considering high wave numbers?, SIAM Rev. 42 (2000) 451–484. [2] A. Bayliss, C. Goldstein, E. Turkel, An iterative method for Helmholtz equation, J. Comput. Phys. 49 (1983) 631–644. [3] J. B´erenger, A perfectly matched layer for the absorption of electromagnetic waves, J. Comput. Phys. 114 (1994) 185–200. [4] Z. Chen, D. Cheng, W. Feng, T. Wu, H. Yang, A multigrid-based preconditioned Krylov subspace method for the Helmholtz equation with PML, J. Math. Anal. Appl. 383 (2011) 522–540. [5] Z. Chen, H. Wu, An adaptive finite element method with perfectly matched absorbing layers for the wave scattering by periodic structures, SIAM J. Numer. Anal. 41 (2003) 799–826. [6] Z. Chen, T. Wu, H. Yang, An optimal 25-point finite difference scheme for the Helmholtz equation with PML, J. Comput. Appl. Math. 236 (2011) 1240–1258. [7] A. Deraemaeker, I. Babuˇska, P. Bouillard, Dispersion and pollution of the FEM sollution for the Helmholtz equation in one, two and three dimensions, Int. J. Numer. Methods Eng. 46 (1999) 471–499. [8] P. de Zeeuw, Matrix-dependent prolongations and restrictions in a blackbox multigrid solver, J. Comput. Appl. Math. 33 (1990) 1–27. [9] I. Duff, S. Gratton, X. Pinel, X. Vasseur, Multigrid based preconditioners for the numerical solution of two-dimensional heterogeneous problems in geophysics, Int. J. Comput. Math. 84 (2007) 1167–1181. [10] B. Engquist, A. Majda, Radiation boundary conditions for acoustic and elastic wave calculations, Comm. Pure Appl. Math. 32 (1979) 314–358. [11] Y. Erlangga, C. Oosterlee, C. Vuik, A novel multigrid based preconditioner for heterogeneous Helmholtz problems, SIAM J. Sci. Comput. 27 (2006) 1471–1492. [12] Y. Erlangga, C. Vuik, C. Oosterlee, On a class of preconditioners for the Helmholtz equation, Appl. Numer. Math. 50 (2004) 629–651. [13] X. Feng, Absorbing boundary conditions for electromagnetic wave propagation, Math. Comp. 68 (1999) 145–168. [14] X. Feng, H. Wu, Discontinuous galerkin methods for the Helmholtz equation with large wave number, SIAM J. Numer. Anal. 47 (2009) 2872–2896. [15] J. Gozani, A. Nachshon, E. Turkel, Conjuate gradient compiled with multigrid for an indefinite problems, in: R. Vichnevestsky, R.S. Tepelman (Eds.), Advances in Computer Methods for Partial Differential Equation, vol. V, IMACS, New Brunswick, NJ, 1984, pp. 425–427. [16] F. Hastings, J. Schneider, S. Broschat, Application of the perfectly matched layer (PML) absorbing boundary condition to elastic wave propagation, J. Acoust. Soc. Am. 100 (1996) 3061–3069. [17] P. Hooghiemstra, D. Van der Heul, C. Vuik, Application of the shifted-Laplace preconditioner for iterative solution of a higher order finite element discretisation of the vector wave equation: First experiences, Appl. Numer. Math. 60 (2010) 1157–1170. [18] B. Hustedt, S. Operto, J. Virieux, Mixed-grid and staggered-grid finite difference methods for frequency-domain acoustic wave modelling, Geophys. J. Int. 157 (2004) 1269–1296. [19] F. Ihlenburg, I. Babuˇska, Dispersion analysis and error estimation of galerkin finite element methods for the Helmholtz equation, Int. J. Numer. Methods Eng. 38 (1995) 4207–4235. [20] C. Jo, C. Shin, J. Suh, An optimal 9-point, finite-difference, frequency-space, 2-D scalar wave extrapolator, Geophysics 61 (1996) 529–537. [21] H. Knibbe, C. Oosterlee, C. Vuik, GPU implementation of a Helmholtz Krylov solver preconditioned by a shifted Laplace multigrid method, J. Comput. Appl. Math. 236 (2011) 281–293. [22] A. Laird, M. Giles, Preconditioned Iterative Solution of the 2D Helmholtz equation, Report 02/12, Oxford Computer Laboratory, Oxford, UK, 2002. [23] S. Operto, J. Virieux, P. Amestoy, J. L’Excellent, L. Giraud, H. Ali, 3D finite-difference frequency-domain modeling of visco-acoustic wave propagation using a massively parallel direct solver: A feasibility study, Geophysics 72 (2007) 195–211. [24] C. Riyanti, A. Kononov, Y. Erlangga, C. Vuik, C. Oosterlee, R. Plessix, W. Mulder, A parallel multigrid-based preconditioner for the 3D heterogeneous high-frequency Helmholtz equation, J. Comput. Phys. 224 (2007) 431–448. [25] Y. Saad, W. Schultz, GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear system, SIAM J. Sci. Stat. Comput. 7 (1986) 1–176. [26] C. Shin, H. Sohn, A frequency-space 2-D scalar wave extrapolator using extended 25-point finite-difference operator, Geophysics 63 (1998) 289–296. [27] I. Singer, E. Turkel, High-order finite difference methods for the Helmholtz equation, Comput. Methods Appl. Mech. Engrg. 163 (1998) 343–358. [28] I. Singer, E. Turkel, A perfectly matched layer for the Helmholtz equation in a semi-infinite strip, J. Comput. Phys. 201 (2004) 439–465. [29] G. Sutmann, Compact finite difference schemes of sixth order for the Helmholtz equation, J. Comput. Appl. Math. 203 (2007) 15–31. ˇ [30] I. Stekl, R. Pratt, Accurate viscoelastic modeling by frequency-domain finite differences using rotated operators, Geophysics 63 (1998) 1779–1794. [31] N. Umetani, S. MacLachlan, C. Oosterlee, A multigrid-based shifted Laplacian preconditioner for a fourth-order Helmholtz discretization, Numer. Linear Algebra Appl. 16 (2009) 603–626. [32] H. Van Der Vorst, Bi-CGSTAB: A fast and smoothly converging variant of Bi-CG for the solution of nonsymmetric linear systems, SIAM J. Sci. Stat. Comput. 13 (1992) 631–644. [33] M. Van Gijzen, Y. Erlangga, C. Vuik, Spectral analysis of the discrete Helmholtz operator with a shifted Laplacian, SIAM J. Sci. Comput. 29 (2007) 1942–1958.