- Email: [email protected]
A dispersion minimizing subgridding finite difference scheme for the Helmholtz equation with PML
Journal of Computational and Applied Mathematics 267 (2014) 82–95
Contents lists available at ScienceDirect
Journal of Computational and Applied Mathematics journal homepage: www.elsevier.com/locate/cam
A dispersion minimizing subgridding finite difference scheme for the Helmholtz equation with PML✩ Tingting Wu a,∗ , Zhongying Chen b a
School of Mathematical Sciences, Shandong Normal University, Jinan 250014, PR China
b
Guangdong Province Key Laboratory of Computational Science, Sun Yat-sen University, Guangzhou 510275, PR China
article
abstract
info
Article history: Received 25 June 2013 Received in revised form 1 December 2013 Keywords: Helmholtz equation PML Subgridding Numerical dispersion
In this paper, we present a dispersion minimizing subgridding finite difference scheme for solving the Helmholtz equation with perfectly matched layer (PML) in the two dimensional domain, which is a second order scheme and pointwise consistent with the equation. Subgrids are used to discretize the computational domain, and a refined choice strategy based on minimizing the numerical dispersion is proposed for choosing weight parameters for transitional nodes. Numerical experiments are given to illustrate that the newly proposed schemes can produce highly accurate seismic modeling results with enhanced efficiency, compared to uniform grids. © 2014 Elsevier B.V. All rights reserved.
1. Introduction In a Cartesian coordinate system, the 2D Helmholtz equation is
∆u + k2 u = 0,
on R2 ,
(1.1)
where ∆ := ∂ /∂ x + ∂ 2 /∂ y2 is the Laplacian, k is the wavenumber defined as k := 2π f /v in which f and v represent the frequency and the velocity respectively, and u is the unknown indicating the pressure of the wave field. The above equation appears directly or indirectly in almost all wave-related problems arisen from many science, engineering, and industry applications. Therefore, solving the Helmholtz equation, in various forms, has always been the hot topic in wave computation (cf. [1,2] and the references therein). To compute the solution of the above problem, due to finite memory and speed limitation of the computer, we need first to formulate the problem as a finite domain problem. Applying the PML technique to truncate the infinite domain into a bounded rectangular domain and taking the same notations as in [3] lead to the equation 2
2
∂ ∂u ∂ ∂u A + B + Ck2 u = 0, ∂x ∂x ∂y ∂y
(1.2)
where A :=
ey ex
,
B :=
ex ey
,
C := ex ey ,
✩ This research is partially supported by the Natural Science Foundation of China under grants 11301310, 11226303, 10771224 and 11071264, Guangdong Provincial Government of China through the ‘‘Computational Science Innovative Research Team’’ program 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] (T. Wu), [email protected] (Z. Chen).
http://dx.doi.org/10.1016/j.cam.2014.01.031 0377-0427/© 2014 Elsevier B.V. All rights reserved.
T. Wu, Z. Chen / Journal of Computational and Applied Mathematics 267 (2014) 82–95
83
σ
in which ex := 1 − i σωx , ey := 1 − i ωy , ω := 2π f is the angular frequency, σx is a function only of x defined as
σx :=
2 π a0 f 0 0,
lx LPML
2
,
inside PML,
(1.3)
outside PML,
with the dominant frequency f0 of the source, the thickness LPML of PML, the distance lx from the point (x, y) inside PML and the interface between the interior region and PML region, and a constant a0 , and σy is defined similarly. Eq. (1.2) can be seen as a general form of the Helmholtz equation (1.1) with the corresponding PML equation, since in the interior domain, A = B = C = 1. We call it the Helmholtz equation with PML. Moreover, Eq. (1.2) is also valid for variable k(x, y) (see, [4,5]). The finite difference method is a popular and powerful computational technique for numerical seismic wave propagation modeling (cf. [6,7]). It is both easy to implement and computationally efficient. Furthermore, we can easily extend finite difference methods to the three dimensional cases, and improve the accuracy of the solution of the Helmholtz equation by choosing optimal parameters in finite difference formulas [3,8–12]. However, one of the problems with the finite difference method is their inability to accurately handle different grid spacings. For many cases, there may exist small isolated regions with a finer resolution requirement than the rest of the problem space. These locations occur when there is a high material property contrast. When using uniform structured grids, we must use the grid size established by the slowest velocity. As a result, heterogeneous models are always oversampled, which leads to unnecessary extra computations. Using unstructured grids, adaptive to velocity variations, is an attractive method to discretize earth models but it does not always help speed up the computations. For realistic size models, construction of unstructured grids takes computing time and representation of these grids consumes significant amount of computer memory. Efficient implementation of accurate numerical methods for general irregular grids is also a challenging task. A simple and natural alternative to this problem is to use structured nonuniform grids. There are mainly two groups of nonuniform grid finite difference algorithms for wave modeling (see [13] and the references therein): the first group uses discontinuous grids, while the second group uses continuous grids but varying grid spacings. Among the first group, subgridding algorithms are popular in the computation, which usually introduce locally refined meshes into the base mesh (cf. [14–16]). As subgridding consists in refining the mesh locally, it can lead to major computational savings. In computation, conventional uniform grid finite difference schemes can be directly applied in each block of the grids. Along the boundaries of blocks, interpolation schemes are often used to facilitate finite difference calculations. However, common problems associated with the subgridding algorithm in the time domain are spurious reflections, unconditional instabilities and accuracy problems. Spurious reflections and accuracy problems occur mainly because of aliasing between fine and coarse representations due to suboptimal application of interpolation/decimation, and numerical impedance mismatch due to the different discrete impedances at the coarse and fine grids. The quality of subgridding can be measured by observing the level of reflections from the dense to coarse mesh interface. For lowreflection levels, the accuracy of the subgridding is the highest. To alleviate these problems, considerable effort has been devoted (see [17,18,14–16,19] and the references therein). For example, to improve the accuracy and reduce the reflections, algorithms for the interpolation/decimation between the fine and coarse interfaces have been modified, which may depend on polynomial approximations (see, [19]). Based on minimizing the numerical dispersion, J.W. Nehrbass et al. proposed a systematic method that shows how to optimally choose the finite difference coefficients at the interface between the course grid and the fine grid to approximate the Helmholtz equation (see, [16]). This method maintains the second order accuracy of the original finite difference method, while the method based on linear interpolation/decimation is only first order accurate. Moreover, this theory is general and can be applied to any dimension and any desired stencil in a straightforward manner. However, the computation of the coefficients is too complicated, which may not be a good choice for practical cases, and the PML technique is not considered. Y. Li et al. did some improvements for the method proposed in [16] based on the idea of the rotated finite difference method (cf. [14]). This scheme is second order in accuracy and computation for its coefficients are implemented easily. However, we cannot extend this method to the Helmholtz equation with PML, as it will not possess the property of consistency (see, [3,9,10]). In addition, Kristek et al. proposed to downsample the wave fields in fine grid blocks for better numerical stability (see, [18]). For subgridding methods in the frequency domain, there does not exist the stability problem. We next focus on nonuniform grid finite difference algorithms with continuous grids. The readers are referred to the paper [20] and [21] for the detailed construction of continuous grids. Different from the discontinuous grids, no interpolations are required. However, finite difference operators for each grid node where the nodal distribution pattern changes must be designed. In general, the classic Taylor series expansion approach (see, [20]) or using polynomial interpolators (see, [21]) is applied to compute the finite difference operators for nonuniform grids. From the construction of continuous grids, we know that continuous grids still lead to unnecessary extra computations for some cases. In practice, the above mentioned two groups of methods are often used together for better flexibility and efficiency (cf. [22]). In this paper, subgrids will be used to discretize the computational domain, including the interior domain and the PML. Then, with the ghost nodes we will propose a new subgridding finite difference scheme, which not only maintains the second accuracy of the original method but also is consistent with the Helmholtz equation with PML. Furthermore, to reduce reflections from the fine to coarse mesh interface, a new method for choosing optimal coefficients for the proposed scheme will be
84
T. Wu, Z. Chen / Journal of Computational and Applied Mathematics 267 (2014) 82–95
a
b
x
y
Fig. 1. (a) Subgridding method for a two-layered model, and (b) the first kind node.
given. This method is based on minimizing the numerical dispersion and easily implemented. The only restriction is that the material properties within the transition region must be homogeneous. We shall address this issue in detail in this paper. This paper is organized as follows. In Section 2, we use subgrids to discretize the computational domain, including the interior domain and the PML. The optimal 9-point finite difference scheme proposed in [3] is directly applied in each block. For the transitional node in the interior domain, the finite difference equation is formulated with the ghost nodes, and it is proved that it is consistent with the Helmholtz equation and is a second order scheme. Then, a new refined choice strategy for choosing optimal parameters of the scheme based on minimizing the numerical dispersion is proposed. In Section 3, numerical experiments are given to demonstrate the efficiency of the scheme. These numerical results show that the newly proposed schemes can produce highly accurate seismic modeling results with enhanced efficiency, compared to uniform grids. Finally, Section 4 contains the conclusions of this paper. 2. A dispersion minimizing subgridding finite difference scheme for the 2D Helmholtz equation with PML In this section, we firstly use subgrids to discretize the computational domain. Then, we propose a dispersion minimizing subgridding finite difference scheme for 2D Helmholtz equation with PML based on the following two aspects: firstly, as the domain is divided into different blocks, the optimal 9-point finite difference scheme proposed in [3], which is indeed a dispersion minimizing method, is directly applied in each block; secondly, for transitional nodes in the interior domain, we propose a consistent finite difference method with ghost nodes, and present a new refined optimization rule for choosing the parameters of the finite difference scheme based on minimizing the numerical dispersion. To construct finite difference schemes, we firstly apply the subgridding technique to discretize the computational domain. In detail, the mesh for the computational domain is established by the following steps: first, a coarse mesh is established and placed throughout the computational domain; second, sensitive regions are defined where a finer resolution is desired; third, in these sensitive regions, the resolution is doubled so that every coarse cell is divided into four finer cells; finally, for the consistency of the finite difference method, the mesh in some parts of the perfectly layer is also refined. The region in which adjacent cells do not have the same spatial resolution is referred to as the transition region. If still finer resolutions are required, the process is again repeated. As an example, in Fig. 1(a), we present a two-layered model, in which the top layer needs to be discretized with fine grids, compared with the bottom layer. For the two-layered model, the mesh established by the red lines is the base mesh, and the green lines are added to refine the mesh locally. In addition, the domain surrounded by the blue rectangle is our interested area (the interior area), while the rest is the perfectly matched layer. As is shown in Fig. 1(a), the transitional grids can be categorized into five groups. We will present consistent finite difference methods for them one by one, and propose choice strategies for optimal parameters for their corresponding finite difference stencils. 2.1. A consistent finite difference scheme for the first kind node We next present a consistent finite difference method for the first kind of transitional nodes, which is also second order in accuracy. Seen from Fig. 1(a), we find that the first kind of transitional nodes lie in the interior domain, in which the Helmholtz equation with PML (1.2) deteriorates to be the Helmholtz equation, as A = B = C = 1 holds for the interior area. With three ghost nodes, we firstly construct a symmetric 9-point stencil for the first kind node. This stencil is shown in Fig. 1(b), where the nodes in the bottom red line are ghost nodes. Then, we formulate a finite difference scheme for the first kind node based on the averaged idea (see, [3,9,10] and the references therein). For simplicity, we introduce some notations here. Let h and H respectively represent the step size for the fine grid and the coarse grid, and define xm := x0 + (m − 1)h and yn := y0 + (n − 1)h. Let um,n := u|x=xm ,y=yn represent the pressure of the wave field at the location (xm , yn ), and similarly let km,n := k|x=xm ,y=yn . Next, different from the paper [3],
T. Wu, Z. Chen / Journal of Computational and Applied Mathematics 267 (2014) 82–95
85
we introduce more parameters to relax the restriction for the symmetry of the finite difference stencil. In detail, let
Lh,x u|(m,n+j) :=
um+1,n+j + um−1,n+j − 2um,n+j h2
,
(2.1)
for j ∈ Z2 := {−1, 0, 1}, and define
L˜ h,x u|x=xm ,y=yn := a1 Lh,x u|(m,n−1) + a2 Lh,x u|(m,n) + a3 Lh,x u|(m,n+1) ,
(2.2) ∂2u
where a1 , a2 , a3 are constants to be determined in the following. Then we approximate ∂ x2 as
∂ 2 u ≈ L˜ h,x u|x=xm ,y=yn . ∂ x2 x=xm ,y=yn 2 We deal with the approximation of ∂∂ y2u in a similar way, that is,
∂ 2 u ≈ L˜ h,y u|x=xm ,y=yn , ∂ y2 x=xm ,y=yn where
L˜ h,y u|x=xm ,y=yn := b1 Lh,y u|(m−1,n) + b2 Lh,y u|(m,n) + b3 Lh,y u|(m+1,n) ,
(2.3)
in which b1 , b2 , b3 are constants to be determined in the following. Let L˜ h := L˜ h,x + L˜ h,y . Moreover, k2 u|x=xm ,y=yn is approximated by a weighted average: k2 u|x=xm ,y=yn ≈ I˜h k2 u |x=xm ,y=yn ,
where
I˜h k2 u |x=xm ,y=yn := c1 k2m−1,n−1 um−1,n−1 + c2 k2m,n−1 um,n−1 + c3 k2m+1,n−1 um+1,n−1
+ c4 k2m−1,n um−1,n + c5 k2m,n um,n + c6 k2m+1,n um+1,n + c7 k2m−1,n+1 um−1,n+1 + c8 k2m,n+1 um,n+1 + c9 k2m+1,n+1 um+1,n+1 ,
(2.4)
9
in which i=1 ci = 1. These yield the 9-point difference approximation with ghost nodes for the Helmholtz equation as
L˜ h u|x=xm ,y=yn + I˜h k2 u |x=xm ,y=yn = 0.
(2.5)
To give a consistent finite difference scheme for the first kind node, we next present the convergence analysis for the 9-point difference scheme (2.5). Proposition 2.1. (i) If a1 + a2 + a3 = 1, b1 + b2 + b3 = 1 and i=1 ci = 1, then the 9-point finite difference approximation (2.5) is pointwise consistent with the Helmholtz equation (1.1) and is a first order scheme. (ii) If a1 + a2 + a3 = 1, b1 + b2 + b3 = 1, 9i=1 ci = 1, a1 = a3 , b1 = b3 , −c1 + c3 − c4 + c6 − c7 + c9 = 0, and
9
−
3
9
ci + i=7 ci = 0, then the 9-point finite difference approximation (2.5) is pointwise consistent with the Helmholtz equation (1.1) and is a second order scheme. i =1
Proof. Assume that xm ≤ x < xm+1 and yn ≤ y < yn+1 . It follows from the Taylor theorem that
L˜ h u|x=xm ,y=yn = (a1 + a2 + a3 )
∂ 2u ∂ 2u + (b1 + b2 + b3 ) 2 + ν1 h + O (h2 ), 2 ∂x ∂y
(2.6)
and
2
I˜h k u |x=xm ,y=yn =
9
ci
k2 u + ν2 h + O (h2 ),
i=1
where
ν1 := (a3 − a1 )
∂ 3u ∂ 3u + ( b − b ) , 3 1 ∂ x2 ∂ y ∂ x∂ y2
3 9 ∂(ku) ∂(ku) ν2 := (−c1 + c3 − c4 + c6 − c7 + c9 ) + − ci + ci . ∂x ∂y i=1 i =7
(2.7)
86
T. Wu, Z. Chen / Journal of Computational and Applied Mathematics 267 (2014) 82–95
Therefore, under the assumptions of part (i) in this proposition, we have that the left hand side of the 9-point finite difference approximation (2.5) is equivalent to
∂ 2u ∂ 2u + 2 + k2 u + ηh + O (h2 ), ∂ x2 ∂y
(2.8)
where η := ν1 + ν2 . From (2.8) and the concept of consistency, we conclude the results of this proposition.
Now, we show how to construct a second-order finite difference scheme for the first kind node in Fig. 1(a), which is also consistent with the Helmholtz equation. Substituting the Eqs (2.2)–(2.4) into Eq. (2.5) and replacing um+i,n+j with Um+i,n+j (i, j ∈ Z2 ) yield the 9-point finite difference equation T1 Um−1,n−1 + T2 Um,n−1 + T3 Um+1,n−1 + T4 Um−1,n + T5 Um,n + T6 Um+1,n
+ T7 Um−1,n+1 + T8 Um,n+1 + T9 Um+1,n+1 = 0,
(2.9)
in which the coefficients are given by T1 := T4 := T7 :=
a1 + b 1 h2
+ c1 k2m−1,n−1 ,
a2 − 2b1 h2 a3 + b 1 h2
T2 :=
−2a1 + b2 h2
+ c4 k2m−1,n ,
T5 := −
+ c7 k2m−1,n+1 ,
T8 :=
2a2 + 2b2 h2
−2a3 + b2 h2
+ c2 k2m,n−1 , + c5 k2m,n ,
T3 := T6 :=
+ c8 k2m,n+1 ,
a1 + b3
h2 a2 − 2b3
+ c6 k2m+1,n ,
h2
T9 :=
a3 + b3 h2
+ c3 k2m+1,n−1 ,
+ c9 k2m+1,n+1 .
To guarantee the above scheme to be second order in accuracy, the parameters in Eq. (2.9) should satisfy the conditions proposed in Proposition 2.1. Moreover, as the ghost nodes do not exist in the stencil, we should let the parameters of the ghost nodes to be zero, that is, T7 = 0, T8 = 0, T9 = 0, which indicates that the parameters may vary according to the wavenumber k. Choice strategies for optimal parameters of the finite difference scheme will be present in the next subsection. 2.2. Choice strategies for optimal parameters for the first kind node’s stencil We now go further into the problem of the choice strategies for optimal parameters for the first kind node’s stencil. To optimize the 9-point scheme (2.9) for the first transitional node (see, Fig. 1(b)), we first perform classical dispersion analysis by assuming a plane-wave solution of the form U (x, y) = e−ik(x cos θ+y sin θ ) , where θ is the propagation angle from the y-axis. The following analysis is based on the assumption that the wavenumber k is a positive constant. Moreover, let v be the velocity of propagation, λ be the wavelength, and G be the number of gridpoints per wavelength, that is, G = λh . Since and k = ωv , we have kh = 2Gπ . Substituting Um,n := e−ik(xm cos θ+yn sin θ) into the Eq. (2.9), dividing both sides by the factor e−ik(xm cos θ+yn sin θ ) and finally applying the Euler formula eix = cos x + i sin x lead to the following dispersion equation
λ=
2πv
ω
2 [cos(kh cos θ ) − 1] a1 eikh sin θ + a2 + a3 e−ikh sin θ + [cos(kh sin θ ) − 1] b1 eikh cos θ + b2 + b3 e−ikh cos θ
+ k2 h2 c1 eikh(cos θ+sin θ) + c2 eikh sin θ + c3 e−ikh(cos θ −sin θ) + c4 eikh cos θ + c5 + c6 e−ikh cos θ + c7 e−ikh(− cos θ +sin θ) + c8 e−ikh sin θ + c9 e−ikh(cos θ+sin θ) = 0.
(2.10)
As stated in Proposition 2.1, to obtain a second order scheme, the parameters in the Eq. (2.9) should satisfy a1 + a2 + a3 = 1,
b1 + b2 + b3 = 1,
9
ci = 1,
(2.11)
i=1
a1 = a3 ,
b1 = b3 ,
−c1 + c3 − c4 + c6 − c7 + c9 = 0,
−
3 i =1
ci +
9
ci = 0.
(2.12)
i=7
Also, all the parameters of ghost nodes in Fig. 1(b) should be zero, that is T7 = T8 = T9 =
a3 + b 1
+ c7 k2 = 0,
h2 −2a3 + b2 h2 a3 + b 3 h2
+ c8 k2 = 0,
+ c9 k2 = 0.
(2.13) (2.14) (2.15)
T. Wu, Z. Chen / Journal of Computational and Applied Mathematics 267 (2014) 82–95
87
Combining the above Eqs (2.11)–(2.15) yields
a3 = a1 , a2 = 1 − 2a1 , b3 = b1 , b2 = 1 − 2b1 , c8 k2 h2 = 2(a1 + b1 ) − 1, c9 k2 h2 = −(a1 + b1 ), c1 = c3 + c5 + 2c6 + 2c8 + 4c9 − 1, c2 = −2c3 − c5 − 2c6 − c8 − 2c9 + 1, c4 = −c5 − c6 − 2c8 − 4c9 + 1, c7 = c9 .
(2.16)
From Eq. (2.16), we see that for the first kind node, the symmetry of the finite difference stencil is destroyed. Hence, different from the paper [3], the dispersion equation (2.10) has the imaginary part. Therefore, we cannot obtain the numerical phase velocity or the numerical group velocity any more. Hence, a proper method is needed to get optimal parameters for the first kind node. We next present a new way to obtain the optimal parameters of the finite difference equation (2.9) for the first transitional node. Firstly, we substitute the Eq. (2.16) into (2.10) to eliminate the parameters a2 , a3 , b2 , b3 , c1 , c2 , c4 , c7 , c8 , c9 . Then, under the assumption that all the parameters in the finite difference equation (2.9) are real numbers, we separate the dispersion equation (2.10) into the real part and the imaginary part. The real part of the dispersion equation (2.10) 2 [cos (kh cos θ ) cos (kh sin θ ) − cos (kh sin θ) − 2 cos (kh cos θ) + 2] (a1 + b1 )
+ k2 h2 2 cos (kh sin θ ) [cos (kh cos θ ) − 1] c3 + [cos (kh(cos θ + sin θ )) − cos (kh sin θ) − cos (kh cos θ ) + 1]c5 + 2 [cos (kh(cos θ + sin θ )) − cos (kh sin θ )] c6 = −4 cos (kh cos θ ) − 2 cos (kh sin θ ) + 4 + 2 cos (kh(cos θ + sin θ )) − k2 h2 [cos (kh(cos θ + sin θ )) − cos (kh sin θ) − cos (kh cos θ)] ,
(2.17)
and the corresponding imaginary part is
2 sin (kh sin θ) [cos (kh cos θ ) − 1] (a1 + b1 ) + k2 h2 2 sin (kh sin θ) [cos (kh cos θ) − 1]c3
+ [sin (kh (cos θ + sin θ)) − sin (kh sin θ ) − sin (kh cos θ)](c5 + 2c6 ) = [sin (kh (cos θ + sin θ )) − sin (kh sin θ ) − sin (kh cos θ)](2 + k2 h2 ).
(2.18)
Note that Eqs. (2.17) and (2.18) are functions of kh and θ . To reduce the numerical dispersion and improve the accuracy of the difference scheme for the first kind node, we propose the following rule. Rule 2.2 (Refined Choice Strategy). Step 1. Estimate kh, where k is the wavenumber for the first kind node and h is the step size for the fine grid. (m−1)2π ∈ Iθ , m = 1, 2, . . . , l, where l is a given positive integer. Step 2. Let Iθ := [0, 2π ), and choose θ = θm := l Step 3. Inserting above kh and θ to Eqs. (2.17) and (2.18) yields a linear system. Step 4. Applying the singular value decomposition method to solve the corresponding system yields the optimal weight parameters for Eqs. (2.17) and (2.18). Other parameters are obtained by Eq. (2.16). Different from the paper [3], here we let the number of gridpoints per wavelength G be a fixed number, which seems to be a proper choice, as in the above analysis we assume the wavenumber to be constant and we can obtain other parameters only by Eq. (2.16). In addition, the angle θ belongs to the interval [0, 2π ), which arises from that the symmetry of the finite difference stencil is destroyed. As an example, in Fig. 2, we present some optimal parameters for the first kind of transitional nodes. In this figure, we only show the parameters which appear in the Eqs. (2.17) and (2.18), as other parameters can be easily obtained by Eq. (2.16). We find that the parameters a1 , b1 , c3 vary smoothly, while the parameters c5 and c6 are very sensitive to the number of gridpoints per wavelength G, which may be the effect of the ghost nodes. 2.3. Optimal 9-point finite difference schemes for other nodes In this subsection, we firstly present optimal 9-point finite difference methods for the last four kinds of transitional nodes, which are also second order in accuracy. Then, based on the paper [3], we formulate optimal 9-point finite difference scheme for other nodes, which are not transitional nodes. For the second kind node, which also lie in the interior domain, we construct its consistent finite difference scheme with its nearest 8 coarse points, which is different from the first kind node. For the third kind of transitional nodes, we construct its consistent finite difference method with its nearest 8 symmetric nodes. The finite difference schemes for the fourth and fifth kind of transitional nodes are formulated similarly. We refer the interested readers to the Eqs (3.2) and (3.25) in [3] for the details. All these schemes are consistent with the Helmholtz equation with PML and are second order in accuracy (see, [3]). In addition, the parameters of the finite difference methods can be chosen with refined choice strategy (see, Rule 3.6 and Rule 3.8 in [3]).
88
T. Wu, Z. Chen / Journal of Computational and Applied Mathematics 267 (2014) 82–95
Fig. 2. Optimal coefficients for the first kind node.
a
b
Fig. 3. (a) A near-surface low-velocity layer model. (b) A low-velocity interlayer model.
For other nodes which are not transitional nodes, we formulate the optimal 9-point finite difference scheme for the Helmholtz equation with PML based on our previous work [3]. Finally, together with the work for the transitional nodes, we have a dispersion minimizing subgridding finite difference scheme for the 2D Helmholtz equation with PML, referred to as the refined optimal subgridding 9-point finite difference scheme for short, which is second order in accuracy and consistent with the equation. 3. Numerical experiments In this section, we present three numerical experiments to illustrate the efficiency of the schemes described in the last section. All the experiments are performed with Matlab 7v on an Intel Xeon (2-core) with 2.40 GHz and 2.32 Gb RAM. 3.1. A near-surface low-velocity layer model To validate the efficiency of the refined optimal subgridding 9-point finite difference scheme (the refined subgridding 9p) over the refined optimal 9-point finite difference scheme based on uniform grids proposed in [3] (the refined uniform 9p) for the case of variational wavenumbers, which often results from heterogeneous medium, we consider a near-surface low-velocity layer model shown in Fig. 3(a). The domain is [0 m, 1600 m] × [0 m, 1600 m] and there are four velocities: from the top, v := 1200 m/s, 2000 m/s, 2400 m/s and 2800 m/s. A point source δ(x − xs , y − ys ) R(f , f0 ) is located at the point (xs , ys ) = (0 m, 800 m), where R(f , f0 ) is the Ricker wavelet in the frequency domain defined by
R(f , f0 ) :=
+∞ −∞
(1 − 2π 2 f02 t 2 ) exp(−π 2 f02 t 2 )e−i2π ft dt ,
T. Wu, Z. Chen / Journal of Computational and Applied Mathematics 267 (2014) 82–95
89
whose dominant frequency is f0 = 25 Hz. In addition, all the receivers are placed on the top of the ground, that is, they are in the line of y = 0. In the following computation, the time sampling is 8 ms and the highest frequency we compute is 60 Hz. For the refined uniform 9p, the uniform step size hu := ∆x = ∆y is used for both variables x and y. For the refined subgridding 9p, H and h denote the coarse and fine step size, respectively. In this experiment, we let H = 8 m and h = 4 m. In detail, we only use fine grids for the top layer whose velocity is v = 1200 m/s. We next generate common-shot-point records (shot profiles) for the wave field by the refined uniform 9p and the refined subgridding 9p. The monofrequency wave field (real part) for f = 30 Hz obtained by the refined subgridding 9p with (H , h) = (8 m, 4 m) is present in Fig. 4(a). It is easy to see that PML’s absorption of the outgoing waves is efficient, and almost no boundary reflections exist. Additionally, the downward incident waves and transmissive waves are clear. In the first velocity layer, the incident waves are interfered with the reflected waves returned from the bottom of the first velocity layer, and this phenomena obeys Snell’s law. In Fig. 4(b)–(d), we show the common-shot-point records generated by the refined uniform 9p with uniform grid size hu = 4 m, the refined uniform 9p with hu = 8 m and the refined subgridding 9p with (H , h) = (8 m, 4 m), respectively. In these three pictures, from the top, the first wave represents the direct wave, the second one is the wave reflected by the bottom of the first layer (the low-velocity layer), the third one denotes the wave reflected by the bottom of the second layer, and the fourth wave is the wave reflected by the bottom of the third layer (see, the reference number in Fig. 4(d)). Comparing these three figures, we find that Fig. 4(d) is almost as clear as Fig. 4(b), while there are some nonphysical oscillations in Fig. 4(c), which is caused by the numerical dispersion. To see clearly, Fig. 4(e) and Fig. 4(f) present local enlargement for Fig. 4(c) and Fig. 4(d), respectively. From the last two pictures in Fig. 4, some nonphysical oscillations in the background, such as in the time interval [0.04 s, 0.40 s], can be found in the former, but these do not exist in the latter, which verify that the refined subgridding 9p cannot only improve the accuracy but also suppress nonphysical oscillations. Furthermore, when comparing the cost of Fig. 4(d) with that of Fig. 4(b), there is a 56.44% decrease in memory and a 76.63% decrease in the computational time. Therefore, the efficiency of the refined subgridding 9p is confirmed. 3.2. A low-velocity interlayer model To validate the refined subgridding 9p’s availability for the case of variational wavenumbers, we also consider a lowvelocity interlayer model shown in Fig. 3(b). The domain is [0 m, 1600 m] × [0 m, 1600 m] and there are four velocities: from the top, v := 2000 m/s, 1200 m/s, 2400 m/s and 2800 m/s. A point source, which is the same as the one used in the above example, is located at the point (xs , ys ) = (0 m, 800 m). Also, all the receivers are placed on the top of the ground, the time sampling is 8 ms and the highest frequency we compute is 60 Hz. For the refined subgridding 9p, we let the coarse grid size H = 8 m and the fine grid size h = 4 m. In detail, we only use fine grids for the interlayer whose velocity is v = 1200 m/s. In Fig. 5, we generate common-shot-point records for the wave field by the refined uniform 9p and the refined subgridding 9p. The monofrequency wave field (real part) for f = 35 Hz obtained by the refined subgridding 9p with (H , h) = (8 m, 4 m) is present in Fig. 5(a). It is easy to see that PML’s absorption of the outgoing waves is efficient, and almost no boundary reflections exist. Additionally, the downward incident waves and transmissive waves are clear. In the second velocity layer, the incident waves are interfered with the reflected waves returned from the bottom of the second velocity layer, and this phenomena obeys Snell’s law. In Fig. 5(b)–(d), we show the common-shot-point records generated by the refined uniform 9p with uniform grid size hu = 4 m, the refined uniform 9p with hu = 8 m and the refined subgridding 9p with (H , h) = (8 m, 4 m), respectively. In these three pictures, from the top, the first wave represents the direct wave, the second one is the wave reflected by the bottom of the first layer, the third one denotes the wave reflected by the bottom of the second layer (the low-velocity layer), and the fourth wave is the wave reflected by the bottom of the third layer (see, the reference number in Fig. 5(d)). Comparing these three figures, we find that Fig. 5(d) is almost as clear as Fig. 5(b), while there are some nonphysical oscillations in Fig. 5(c), which is caused by the numerical dispersion. To see clearly, Fig. 5(e) and Fig. 5(f) present local enlargement for Fig. 5(c) and Fig. 5(d), respectively. From the last two pictures in Fig. 5, some nonphysical oscillations in the background, such as in the time interval [0.88 s, 1.20 s], can be found in the former, but these do not exist in the latter, which verify that the refined subgridding 9p cannot only improve the accuracy but also suppress nonphysical oscillations. Furthermore, when comparing the cost of Fig. 5(d) with that of Fig. 5(b), there is a 54.31% decrease in memory and a 77.25% decrease in the computational time. Therefore, the efficiency of the refined subgridding 9p is confirmed. 3.3. A two-layered model To compare the refined subgridding 9p with the refined uniform 9p in a more quantitative way, we consider a two-layered model [23]. For the observation point x = (x, y) and source point xs = (xs , ys ) with ys < 0, we consider the solution G(x, xs ) of the Helmholtz equation in a two-layered medium in R2 which satisfies
∆G(x, xs ) + k2 (x)G(x, xs ) = −δ(x − xs ),
(3.1)
with continuity conditions G(x, xs )|y=0+ = G(x, xs )|y=0− ,
(3.2)
90
T. Wu, Z. Chen / Journal of Computational and Applied Mathematics 267 (2014) 82–95
a
b
c
d
e
f
Fig. 4. (a) The real part of the solution for f = 30 Hz obtained by the refined subgridding 9p with (H , h) = (8 m, 4 m). (b) The common-shot-point records for the wave field by the refined uniform 9p with hu = 4 m. (c) The common-shot-point records for the wave field by the refined uniform 9p with hu = 8 m. (d) The common-shot-point records for the wave field by the refined subgridding 9p with (H , h) = (8 m, 4 m). (e) Local enlargement for (c). (f) Local enlargement for (d).
∂ G(x, xs ) ∂ G(x, xs ) = , ∂ y y=0+ ∂ y y=0−
(3.3)
where the wavenumber k(x) =
k1 , k2 ,
for y > 0, for y < 0,
in which k1 and k2 are given positive numbers satisfying k1 ̸= k2 . Without loss of generality, we assume that k1 < k2 holds for this problem.
T. Wu, Z. Chen / Journal of Computational and Applied Mathematics 267 (2014) 82–95
a
b
c
d
e
f
91
Fig. 5. (a) The real part of the solution for f = 35 Hz obtained by the refined subgridding 9p with (H , h) = (8 m, 4 m). (b) The common-shot-point records for the wave field by the refined uniform 9p with hu = 4 m. (c) The common-shot-point records for the wave field by the refined uniform 9p with hu = 8 m. (d) The common-shot-point records for the wave field by the refined subgridding 9p with (H , h) = (8 m, 4 m). (e) Local enlargement for (c). (f) Local enlargement for (d).
We next present the analytical formula of the solution for the above two-layered problem. It follows from the Fourier transform that the solution is given by [23]
• y < 0, G(x, xs ) = Φ (x, xs ) +
i 4π
+∞
g1 (x, xs , ξ )dξ , −∞
(3.4)
92
T. Wu, Z. Chen / Journal of Computational and Applied Mathematics 267 (2014) 82–95
Table 1 The performance of the refined uniform 9p. The pair (·, ·) denotes the error for the layer with k(x) = 0.1571 and the layer with k(x) = 0.0628, respectively. hu
10 m
5m
2.5 m
1.25 m
n Errors Time
8281 (19.5045, 10.0745) 4.5708 s
25 921 (5.4717, 2.4063) 50.7783 s
68 121 (1.8233, 0.6844) 326.9001 s
231 361 (0.6093,0.2453) 3471.7 s
Table 2 The performance of the refined subgridding 9p. The pair (·, ·) denotes the error for the layer with k(x) = 0.1571 and the layer with k(x) = 0.0628, respectively.
(H , h)
(10 m, 5 m)
(5 m, 2.5 m)
(2.5 m, 1.25 m)
(1.25 m, 0.625 m)
n Errors Time
10 286 (11.8947, 4.8360) 7.0356 s
32 631 (3.6300, 1.2230) 80.4185 s
87 541 (1.1214, 0.3889) 553.7567 s
300 282 (0.4225, 0.1705) 5851.8 s
• y > 0, G(x, xs ) =
i 2π
+∞
g2 (x, xs , ξ )dξ ,
(3.5)
−∞
where
Φ (x, xs ) =
i 4
(1)
H0 (k2 |x − xs |) ,
(3.6)
(1)
in which H0 is the Hankel function of the first kind with order zero, i2 = −1,
2 − k2 − √2 2 ξ ξ 2 − k21 2 ξ −k2 (y+ys ) iξ (x−xs ) e , e 2 2 2 2 2 2 i ξ − k2 ξ − k1 + ξ − k2 √2 2 k22 − ξ 2 − i ξ 2 − k21 e−i k2 −ξ (y+ys ) eiξ (x−xs ) , g1 (x, xs , ξ ) = k22 − ξ 2 i ξ 2 − k21 + k22 − ξ 2 √2 2 k22 − ξ 2 − k21 − ξ 2 e−i k2 −ξ (y+ys ) eiξ (x−xs ) , k22 − ξ 2 k21 − ξ 2 + k22 − ξ 2 and
−√ξ 2 −k2 y+√ξ 2 −k2 y 1 2 s e eiξ (x−xs ) , i ξ 2 − k21 + ξ 2 − k22 √ √ − ξ 2 −k21 y−i k22 −ξ 2 ys e g2 (x, xs , ξ ) = eiξ (x−xs ) , 2 2 i ξ 2 − k1 + k2 − ξ 2 √ √ i k21 −ξ 2 y− k22 −ξ 2 ys e eiξ (x−xs ) , 2 2 2 2 k1 − ξ + k2 − ξ
for |ξ | > k2 ,
for k1 < |ξ | < k2 ,
for |ξ | < k1 ,
for |ξ | > k2 ,
for k1 < |ξ | < k2 ,
for |ξ | < k1 .
Then, we consider a two-layered problem (3.1) with the source point fixed at (250 m, −52 m), k1 = 0.0628 and k2 = 0.1571. The bounded domain [0 m, 500 m] × [−52 m, 448 m] is our interested area. The corresponding solution is given visa (3.4)–(3.6). Thus, numerical integrals have to be dealt with to obtain the solution. In this experiment, we apply the Gauss-type quadratures [24] for weakly singular integrals appearing in (3.4) and (3.5). This kind of quadrature rules possess polynomial order of accuracy, and enable us to obtain the solution with high accuracy. In the following, we take the solution obtained by Gauss-type quadratures as the analytical solution of the two-layered problem. We next compare the performance of the refined uniform 9p and the refined subgridding 9p, including errors, computational times, number of unknowns and the errors’ decay exponents (briefly, decay exp.). The error is measured
T. Wu, Z. Chen / Journal of Computational and Applied Mathematics 267 (2014) 82–95
93
Table 3 The decay exponent for the refined subgridding 9p. n
10 286
32 631
87 541
300 282
Errors Decay exp.
11.8947
3.6300 1.0280
1.1214 1.1903
0.4225 0.7920
a
b
c
d
Fig. 6. (a) The real part of the analytical solution. (b) The real part of the solution generated by the refined uniform 9p with hu = 1.25 m. (c) The real part of the solution generated by the refined subgridding 9p with (H , h) = (2.5 m, 1.25 m). (d) The real part of solutions in the line of x = 125 m generated by different methods.
in C-norm, which is defined as: for any complex vector z = [z1 , z2 , . . . , zM ],
∥z∥C := max |zj |, 1≤j≤M
where |zj | is the complex modulus of zj . Let n represent the number of all unknowns in the two dimensional computational domain in one test, that is, n is the sum of the number of the unknowns in the interior area and that of the unknowns in the ln(en /em ) PML. Each decay exponent is computed by ln(m , where en denotes the numerical solution’s error for one test in which the /n) number of all unknowns is n. Also, for the refined uniform 9p, the uniform step size hu := ∆x = ∆y is used for both variables x and y. For the refined subgridding 9p, H and h denote the coarse and fine step size, respectively. In this experiment, the coarse grids are refined only once in the layer with bigger wavenumber k(x) = 0.1571. Therefore, h = H /2 holds. In Tables 1 and 2, we present the number of all unknowns in the two dimensional computational domain, errors and computational time for the refined uniform 9p and subgridding 9p, respectively. Seen from Tables 1 and 2, for the two methods, the error of the layer with bigger wavenumber is almost two times larger than that of the layer with smaller wavenumber, which arises from the varying of the wavenumber in the two layers. In addition, for the refined subgridding 9p with a pair of given step sizes (H , h), its corresponding error is almost half of that obtained by the refined uniform 9p
94
T. Wu, Z. Chen / Journal of Computational and Applied Mathematics 267 (2014) 82–95
with hu = H, and is almost twice as big as that obtained by the refined uniform 9p with hu = h. However, we find that there is a sharp decrease in both the number of unknowns and the computational time for the refined subgridding 9p with (H , h), when compared with the refined uniform 9p with hu = h. In Table 3, we present the decay exponent for the refined subgridding 9p. We find that the decay exponent is one. As n denotes the number of all unknowns in the two dimensional computational domain, the decay exponent being one indicates that the overall global convergence for the refined subgridding 9p is second order. For further comparison, in Fig. 6(a)–(c) we present the real part of the solution for the two-layered problem in the domain [0 m, 500 m] × [−52 m, 448 m] by the analytical method, the refined uniform 9p with hu = 1.25 m and the refined subgridding 9p with (H , h) = (2.5 m, 1.25 m), respectively. In each of these three pictures, we see that the downward incident waves and transmissive waves are clear. Moreover, in the first layer, the incident waves are interfered with the reflected waves returned from the bottom of the first layer, and this phenomena obeys Snell’s law. It is easy to see that both Fig. 6(b) and Fig. 6(c) are similar with Fig. 6(a), but Fig. 6(c) needs less memory and computational time than those needed by Fig. 6(b) (see, Tables 1 and 2), which verifies the efficiency of the refined subgridding 9p. Also, almost no boundary reflections exist in Fig. 6(b) and Fig. 6(c), which confirms the efficiency of the PML’s absorption. Finally, Fig. 6(d) presents the real part of the two-layered problem’s solution in the line of x = 125 m respectively by the analytical method, the refined uniform 9p with hu = 2.5 m, the refined uniform 9p with hu = 1.25 m, the refined subgridding 9p with (H , h) = (2.5 m, 1.25 m) and the refined subgridding 9p with (H , h) = (1.25 m, 0.625 m). From this figure, we find that the solution obtained by the refined subgridding 9p with (H , h) = (1.25 m, 0.625 m) is a good approximation for the analytical solution. Furthermore, comparing the two solutions obtained by the refined subgridding 9p with different step sizes, we see that the solution obtained by the refined subgridding 9p converges to the analytical solution. In addition, for the results obtained by the refined subgridding 9p, there are nonphysical oscillations near the line y = 0 m, which is the interface between the coarse grids and the fine grids, and this phenomena indicates the validity of our refined choice strategy for the transitional nodes. Now, we come to our conclusion that, when balancing the accuracy against the computational cost, the refined subgridding 9p possesses much improvement over the refined uniform 9p, especially for large-scale layered problems. 4. Conclusions We summarize and comment on the numerical results. Firstly, we applied subgrids to discretize the computational domain, and presented a consistent 9-point subgridding difference scheme for the Helmholtz equation with PML, which is second order in accuracy. Then, for the transitional nodes, we proposed refined choice strategies for choosing optimal parameters based on minimizing the numerical dispersion. Finally, numerical experiments were presented to confirm that the refined subgridding 9p is a good choice for the Helmholtz equation with varying wavenumbers, especially for large-scale layered problems. References [1] R.G. Pratt, C. Shin, G.J. Hicks, Gauss–Newton and full Newton methods in frequency-space seismic waveform inversion, Geophys. J. Int. 133 (1998) 341–362. [2] R.G. Pratt, M.H. Worthington, Inverse theory applied to multi-source cross-hole tompgraphy. Part 1: acoustic wave-equation method, Geophys. Prospect. 38 (1990) 287–310. [3] Z. Chen, D. Cheng, W. Feng, T. Wu, An optimal 9-point finite difference scheme for the Helmholtz equation with PML, Int. J. Numer. Anal. Model. 10 (2013) 389–410. [4] J.W. Kang, L.F. Kallivokas, Mixed unsplit-field perfectly matched layers for transient simulations of scalar waves in heterogeneous domains, Comput. Geosci. 14 (2010) 623–648. [5] J.W. Kang, L.F. Kallivokas, The inverse medium problem in heterogeneous PML-truncated domains using scalar probing waves, Comput. Methods Appl. Mech. Engrg. 200 (2011) 265–283. [6] Y.A. Erlangga, C.W. Oosterlee, C. Vuik, A novel multigrid based preconditioner for heterogeneous Helmholtz problems, SIAM J. Sci. Comput. 27 (2006) 1471–1492. [7] 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. [8] Z. Chen, D. Cheng, W. Feng, T. Wu, H. Yang, A new multigrid based preconditioned Krylov subspace method for the Helmholtz equation with PML, J. Math. Anal. Appl. 383 (2011) 522–540. [9] Z. Chen, D. Cheng, T. Wu, A dispersion minimizing finite difference scheme and preconditioned solver for the 3D Helmholtz equation, J. Comput. Phys. 231 (2012) 8152–8175. [10] 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. [11] C.-H. Jo, C. Shin, J.H. Suh, An optimal 9-point, finite-difference, frequency-space, 2-D scalar wave extrapolator, Geophysics 61 (1996) 529–537. [12] C. Shin, H. Sohn, A frequency-space 2-D scalar wave extrapolator using extended 25-point finite-difference operator, Geophysics 63 (1998) 289–296. [13] C. Chu, P.L. Stoffa, Nonuniform grid implicit spatial finite difference method for acoustic wave modeling in tilted transversely isotropic media, J. Appl. Geophys. 76 (2012) 44–49. [14] Y. Li, L. Sun, W. Hong, Helmholtz equation sub-grid method for multi-transition region, National Microwave Meeting, Qingdao China (2011) pp. 1150–1152. [15] P. Moczo, E. Bystrický, J. Kristek, J.M. Carcione, M. Bouchon, Hybrid modeling of P-SV seismic motion at inhomogeneous viscoelastic topographic structures, Bull. Seismol. Soc. Am. 87 (1997) 1305–1323. [16] J.W. Nehrbass, R. Lee, Optimal finite-difference sub-gridding techniques applied to the Helmholtz equation, IEEE Trans. Microw. Theory Tech. 48 (2000) 976–984. [17] B. Donderici, F.L. Teixeira, Improved FDTD subgridding algorithms via digital filtering and domain overriding, IEEE Trans. Antennas and Propagation 53 (2005) 2938–2951.
T. Wu, Z. Chen / Journal of Computational and Applied Mathematics 267 (2014) 82–95
95
[18] J. Kristek, P. Moczo, M. Galis, Stable discontinuous staggered grid in the finite difference modelling of seismic motion, Geophys. J. Int. 183 (2010) 1401–1407. [19] D.T. Prescott, N.V. Shuley, A method for incorporating different sized cells into the finite-difference time-domain analysis technique, IEEE Microw. Guided Wave Lett. 2 (1992) 434–436. [20] A. Pitarka, 3D elastic finite-difference modeling of seismic motion using staggered grids with nonuniform spacing, Bull. Seismol. Soc. Am. 89 (1999) 54–68. [21] S.A.M. Oliveira, A fourth-order finite-difference method for the acoustic wave equation on irregular grids, Geophysics 68 (2003) 672–676. [22] C. Jastram, E. Tessmer, Elastic modelling on a grid with vertically varying spacing, Geophys. Prospect. 42 (1994) 357–370. [23] J. Coyle, Locating the support of objects contained in a two-layered background medium in two dimensions, Inverse Problems 16 (2000) 275–292. [24] H. Kaneko, Y. Xu, Gauss-type quadratures for weakly singular integrals and their application to Fredholm integral equations of the second kind, Math. Comput. 62 (1994) 739–753.
Contents lists available at ScienceDirect
Journal of Computational and Applied Mathematics journal homepage: www.elsevier.com/locate/cam
A dispersion minimizing subgridding finite difference scheme for the Helmholtz equation with PML✩ Tingting Wu a,∗ , Zhongying Chen b a
School of Mathematical Sciences, Shandong Normal University, Jinan 250014, PR China
b
Guangdong Province Key Laboratory of Computational Science, Sun Yat-sen University, Guangzhou 510275, PR China
article
abstract
info
Article history: Received 25 June 2013 Received in revised form 1 December 2013 Keywords: Helmholtz equation PML Subgridding Numerical dispersion
In this paper, we present a dispersion minimizing subgridding finite difference scheme for solving the Helmholtz equation with perfectly matched layer (PML) in the two dimensional domain, which is a second order scheme and pointwise consistent with the equation. Subgrids are used to discretize the computational domain, and a refined choice strategy based on minimizing the numerical dispersion is proposed for choosing weight parameters for transitional nodes. Numerical experiments are given to illustrate that the newly proposed schemes can produce highly accurate seismic modeling results with enhanced efficiency, compared to uniform grids. © 2014 Elsevier B.V. All rights reserved.
1. Introduction In a Cartesian coordinate system, the 2D Helmholtz equation is
∆u + k2 u = 0,
on R2 ,
(1.1)
where ∆ := ∂ /∂ x + ∂ 2 /∂ y2 is the Laplacian, k is the wavenumber defined as k := 2π f /v in which f and v represent the frequency and the velocity respectively, and u is the unknown indicating the pressure of the wave field. The above equation appears directly or indirectly in almost all wave-related problems arisen from many science, engineering, and industry applications. Therefore, solving the Helmholtz equation, in various forms, has always been the hot topic in wave computation (cf. [1,2] and the references therein). To compute the solution of the above problem, due to finite memory and speed limitation of the computer, we need first to formulate the problem as a finite domain problem. Applying the PML technique to truncate the infinite domain into a bounded rectangular domain and taking the same notations as in [3] lead to the equation 2
2
∂ ∂u ∂ ∂u A + B + Ck2 u = 0, ∂x ∂x ∂y ∂y
(1.2)
where A :=
ey ex
,
B :=
ex ey
,
C := ex ey ,
✩ This research is partially supported by the Natural Science Foundation of China under grants 11301310, 11226303, 10771224 and 11071264, Guangdong Provincial Government of China through the ‘‘Computational Science Innovative Research Team’’ program 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] (T. Wu), [email protected] (Z. Chen).
http://dx.doi.org/10.1016/j.cam.2014.01.031 0377-0427/© 2014 Elsevier B.V. All rights reserved.
T. Wu, Z. Chen / Journal of Computational and Applied Mathematics 267 (2014) 82–95
83
σ
in which ex := 1 − i σωx , ey := 1 − i ωy , ω := 2π f is the angular frequency, σx is a function only of x defined as
σx :=
2 π a0 f 0 0,
lx LPML
2
,
inside PML,
(1.3)
outside PML,
with the dominant frequency f0 of the source, the thickness LPML of PML, the distance lx from the point (x, y) inside PML and the interface between the interior region and PML region, and a constant a0 , and σy is defined similarly. Eq. (1.2) can be seen as a general form of the Helmholtz equation (1.1) with the corresponding PML equation, since in the interior domain, A = B = C = 1. We call it the Helmholtz equation with PML. Moreover, Eq. (1.2) is also valid for variable k(x, y) (see, [4,5]). The finite difference method is a popular and powerful computational technique for numerical seismic wave propagation modeling (cf. [6,7]). It is both easy to implement and computationally efficient. Furthermore, we can easily extend finite difference methods to the three dimensional cases, and improve the accuracy of the solution of the Helmholtz equation by choosing optimal parameters in finite difference formulas [3,8–12]. However, one of the problems with the finite difference method is their inability to accurately handle different grid spacings. For many cases, there may exist small isolated regions with a finer resolution requirement than the rest of the problem space. These locations occur when there is a high material property contrast. When using uniform structured grids, we must use the grid size established by the slowest velocity. As a result, heterogeneous models are always oversampled, which leads to unnecessary extra computations. Using unstructured grids, adaptive to velocity variations, is an attractive method to discretize earth models but it does not always help speed up the computations. For realistic size models, construction of unstructured grids takes computing time and representation of these grids consumes significant amount of computer memory. Efficient implementation of accurate numerical methods for general irregular grids is also a challenging task. A simple and natural alternative to this problem is to use structured nonuniform grids. There are mainly two groups of nonuniform grid finite difference algorithms for wave modeling (see [13] and the references therein): the first group uses discontinuous grids, while the second group uses continuous grids but varying grid spacings. Among the first group, subgridding algorithms are popular in the computation, which usually introduce locally refined meshes into the base mesh (cf. [14–16]). As subgridding consists in refining the mesh locally, it can lead to major computational savings. In computation, conventional uniform grid finite difference schemes can be directly applied in each block of the grids. Along the boundaries of blocks, interpolation schemes are often used to facilitate finite difference calculations. However, common problems associated with the subgridding algorithm in the time domain are spurious reflections, unconditional instabilities and accuracy problems. Spurious reflections and accuracy problems occur mainly because of aliasing between fine and coarse representations due to suboptimal application of interpolation/decimation, and numerical impedance mismatch due to the different discrete impedances at the coarse and fine grids. The quality of subgridding can be measured by observing the level of reflections from the dense to coarse mesh interface. For lowreflection levels, the accuracy of the subgridding is the highest. To alleviate these problems, considerable effort has been devoted (see [17,18,14–16,19] and the references therein). For example, to improve the accuracy and reduce the reflections, algorithms for the interpolation/decimation between the fine and coarse interfaces have been modified, which may depend on polynomial approximations (see, [19]). Based on minimizing the numerical dispersion, J.W. Nehrbass et al. proposed a systematic method that shows how to optimally choose the finite difference coefficients at the interface between the course grid and the fine grid to approximate the Helmholtz equation (see, [16]). This method maintains the second order accuracy of the original finite difference method, while the method based on linear interpolation/decimation is only first order accurate. Moreover, this theory is general and can be applied to any dimension and any desired stencil in a straightforward manner. However, the computation of the coefficients is too complicated, which may not be a good choice for practical cases, and the PML technique is not considered. Y. Li et al. did some improvements for the method proposed in [16] based on the idea of the rotated finite difference method (cf. [14]). This scheme is second order in accuracy and computation for its coefficients are implemented easily. However, we cannot extend this method to the Helmholtz equation with PML, as it will not possess the property of consistency (see, [3,9,10]). In addition, Kristek et al. proposed to downsample the wave fields in fine grid blocks for better numerical stability (see, [18]). For subgridding methods in the frequency domain, there does not exist the stability problem. We next focus on nonuniform grid finite difference algorithms with continuous grids. The readers are referred to the paper [20] and [21] for the detailed construction of continuous grids. Different from the discontinuous grids, no interpolations are required. However, finite difference operators for each grid node where the nodal distribution pattern changes must be designed. In general, the classic Taylor series expansion approach (see, [20]) or using polynomial interpolators (see, [21]) is applied to compute the finite difference operators for nonuniform grids. From the construction of continuous grids, we know that continuous grids still lead to unnecessary extra computations for some cases. In practice, the above mentioned two groups of methods are often used together for better flexibility and efficiency (cf. [22]). In this paper, subgrids will be used to discretize the computational domain, including the interior domain and the PML. Then, with the ghost nodes we will propose a new subgridding finite difference scheme, which not only maintains the second accuracy of the original method but also is consistent with the Helmholtz equation with PML. Furthermore, to reduce reflections from the fine to coarse mesh interface, a new method for choosing optimal coefficients for the proposed scheme will be
84
T. Wu, Z. Chen / Journal of Computational and Applied Mathematics 267 (2014) 82–95
a
b
x
y
Fig. 1. (a) Subgridding method for a two-layered model, and (b) the first kind node.
given. This method is based on minimizing the numerical dispersion and easily implemented. The only restriction is that the material properties within the transition region must be homogeneous. We shall address this issue in detail in this paper. This paper is organized as follows. In Section 2, we use subgrids to discretize the computational domain, including the interior domain and the PML. The optimal 9-point finite difference scheme proposed in [3] is directly applied in each block. For the transitional node in the interior domain, the finite difference equation is formulated with the ghost nodes, and it is proved that it is consistent with the Helmholtz equation and is a second order scheme. Then, a new refined choice strategy for choosing optimal parameters of the scheme based on minimizing the numerical dispersion is proposed. In Section 3, numerical experiments are given to demonstrate the efficiency of the scheme. These numerical results show that the newly proposed schemes can produce highly accurate seismic modeling results with enhanced efficiency, compared to uniform grids. Finally, Section 4 contains the conclusions of this paper. 2. A dispersion minimizing subgridding finite difference scheme for the 2D Helmholtz equation with PML In this section, we firstly use subgrids to discretize the computational domain. Then, we propose a dispersion minimizing subgridding finite difference scheme for 2D Helmholtz equation with PML based on the following two aspects: firstly, as the domain is divided into different blocks, the optimal 9-point finite difference scheme proposed in [3], which is indeed a dispersion minimizing method, is directly applied in each block; secondly, for transitional nodes in the interior domain, we propose a consistent finite difference method with ghost nodes, and present a new refined optimization rule for choosing the parameters of the finite difference scheme based on minimizing the numerical dispersion. To construct finite difference schemes, we firstly apply the subgridding technique to discretize the computational domain. In detail, the mesh for the computational domain is established by the following steps: first, a coarse mesh is established and placed throughout the computational domain; second, sensitive regions are defined where a finer resolution is desired; third, in these sensitive regions, the resolution is doubled so that every coarse cell is divided into four finer cells; finally, for the consistency of the finite difference method, the mesh in some parts of the perfectly layer is also refined. The region in which adjacent cells do not have the same spatial resolution is referred to as the transition region. If still finer resolutions are required, the process is again repeated. As an example, in Fig. 1(a), we present a two-layered model, in which the top layer needs to be discretized with fine grids, compared with the bottom layer. For the two-layered model, the mesh established by the red lines is the base mesh, and the green lines are added to refine the mesh locally. In addition, the domain surrounded by the blue rectangle is our interested area (the interior area), while the rest is the perfectly matched layer. As is shown in Fig. 1(a), the transitional grids can be categorized into five groups. We will present consistent finite difference methods for them one by one, and propose choice strategies for optimal parameters for their corresponding finite difference stencils. 2.1. A consistent finite difference scheme for the first kind node We next present a consistent finite difference method for the first kind of transitional nodes, which is also second order in accuracy. Seen from Fig. 1(a), we find that the first kind of transitional nodes lie in the interior domain, in which the Helmholtz equation with PML (1.2) deteriorates to be the Helmholtz equation, as A = B = C = 1 holds for the interior area. With three ghost nodes, we firstly construct a symmetric 9-point stencil for the first kind node. This stencil is shown in Fig. 1(b), where the nodes in the bottom red line are ghost nodes. Then, we formulate a finite difference scheme for the first kind node based on the averaged idea (see, [3,9,10] and the references therein). For simplicity, we introduce some notations here. Let h and H respectively represent the step size for the fine grid and the coarse grid, and define xm := x0 + (m − 1)h and yn := y0 + (n − 1)h. Let um,n := u|x=xm ,y=yn represent the pressure of the wave field at the location (xm , yn ), and similarly let km,n := k|x=xm ,y=yn . Next, different from the paper [3],
T. Wu, Z. Chen / Journal of Computational and Applied Mathematics 267 (2014) 82–95
85
we introduce more parameters to relax the restriction for the symmetry of the finite difference stencil. In detail, let
Lh,x u|(m,n+j) :=
um+1,n+j + um−1,n+j − 2um,n+j h2
,
(2.1)
for j ∈ Z2 := {−1, 0, 1}, and define
L˜ h,x u|x=xm ,y=yn := a1 Lh,x u|(m,n−1) + a2 Lh,x u|(m,n) + a3 Lh,x u|(m,n+1) ,
(2.2) ∂2u
where a1 , a2 , a3 are constants to be determined in the following. Then we approximate ∂ x2 as
∂ 2 u ≈ L˜ h,x u|x=xm ,y=yn . ∂ x2 x=xm ,y=yn 2 We deal with the approximation of ∂∂ y2u in a similar way, that is,
∂ 2 u ≈ L˜ h,y u|x=xm ,y=yn , ∂ y2 x=xm ,y=yn where
L˜ h,y u|x=xm ,y=yn := b1 Lh,y u|(m−1,n) + b2 Lh,y u|(m,n) + b3 Lh,y u|(m+1,n) ,
(2.3)
in which b1 , b2 , b3 are constants to be determined in the following. Let L˜ h := L˜ h,x + L˜ h,y . Moreover, k2 u|x=xm ,y=yn is approximated by a weighted average: k2 u|x=xm ,y=yn ≈ I˜h k2 u |x=xm ,y=yn ,
where
I˜h k2 u |x=xm ,y=yn := c1 k2m−1,n−1 um−1,n−1 + c2 k2m,n−1 um,n−1 + c3 k2m+1,n−1 um+1,n−1
+ c4 k2m−1,n um−1,n + c5 k2m,n um,n + c6 k2m+1,n um+1,n + c7 k2m−1,n+1 um−1,n+1 + c8 k2m,n+1 um,n+1 + c9 k2m+1,n+1 um+1,n+1 ,
(2.4)
9
in which i=1 ci = 1. These yield the 9-point difference approximation with ghost nodes for the Helmholtz equation as
L˜ h u|x=xm ,y=yn + I˜h k2 u |x=xm ,y=yn = 0.
(2.5)
To give a consistent finite difference scheme for the first kind node, we next present the convergence analysis for the 9-point difference scheme (2.5). Proposition 2.1. (i) If a1 + a2 + a3 = 1, b1 + b2 + b3 = 1 and i=1 ci = 1, then the 9-point finite difference approximation (2.5) is pointwise consistent with the Helmholtz equation (1.1) and is a first order scheme. (ii) If a1 + a2 + a3 = 1, b1 + b2 + b3 = 1, 9i=1 ci = 1, a1 = a3 , b1 = b3 , −c1 + c3 − c4 + c6 − c7 + c9 = 0, and
9
−
3
9
ci + i=7 ci = 0, then the 9-point finite difference approximation (2.5) is pointwise consistent with the Helmholtz equation (1.1) and is a second order scheme. i =1
Proof. Assume that xm ≤ x < xm+1 and yn ≤ y < yn+1 . It follows from the Taylor theorem that
L˜ h u|x=xm ,y=yn = (a1 + a2 + a3 )
∂ 2u ∂ 2u + (b1 + b2 + b3 ) 2 + ν1 h + O (h2 ), 2 ∂x ∂y
(2.6)
and
2
I˜h k u |x=xm ,y=yn =
9
ci
k2 u + ν2 h + O (h2 ),
i=1
where
ν1 := (a3 − a1 )
∂ 3u ∂ 3u + ( b − b ) , 3 1 ∂ x2 ∂ y ∂ x∂ y2
3 9 ∂(ku) ∂(ku) ν2 := (−c1 + c3 − c4 + c6 − c7 + c9 ) + − ci + ci . ∂x ∂y i=1 i =7
(2.7)
86
T. Wu, Z. Chen / Journal of Computational and Applied Mathematics 267 (2014) 82–95
Therefore, under the assumptions of part (i) in this proposition, we have that the left hand side of the 9-point finite difference approximation (2.5) is equivalent to
∂ 2u ∂ 2u + 2 + k2 u + ηh + O (h2 ), ∂ x2 ∂y
(2.8)
where η := ν1 + ν2 . From (2.8) and the concept of consistency, we conclude the results of this proposition.
Now, we show how to construct a second-order finite difference scheme for the first kind node in Fig. 1(a), which is also consistent with the Helmholtz equation. Substituting the Eqs (2.2)–(2.4) into Eq. (2.5) and replacing um+i,n+j with Um+i,n+j (i, j ∈ Z2 ) yield the 9-point finite difference equation T1 Um−1,n−1 + T2 Um,n−1 + T3 Um+1,n−1 + T4 Um−1,n + T5 Um,n + T6 Um+1,n
+ T7 Um−1,n+1 + T8 Um,n+1 + T9 Um+1,n+1 = 0,
(2.9)
in which the coefficients are given by T1 := T4 := T7 :=
a1 + b 1 h2
+ c1 k2m−1,n−1 ,
a2 − 2b1 h2 a3 + b 1 h2
T2 :=
−2a1 + b2 h2
+ c4 k2m−1,n ,
T5 := −
+ c7 k2m−1,n+1 ,
T8 :=
2a2 + 2b2 h2
−2a3 + b2 h2
+ c2 k2m,n−1 , + c5 k2m,n ,
T3 := T6 :=
+ c8 k2m,n+1 ,
a1 + b3
h2 a2 − 2b3
+ c6 k2m+1,n ,
h2
T9 :=
a3 + b3 h2
+ c3 k2m+1,n−1 ,
+ c9 k2m+1,n+1 .
To guarantee the above scheme to be second order in accuracy, the parameters in Eq. (2.9) should satisfy the conditions proposed in Proposition 2.1. Moreover, as the ghost nodes do not exist in the stencil, we should let the parameters of the ghost nodes to be zero, that is, T7 = 0, T8 = 0, T9 = 0, which indicates that the parameters may vary according to the wavenumber k. Choice strategies for optimal parameters of the finite difference scheme will be present in the next subsection. 2.2. Choice strategies for optimal parameters for the first kind node’s stencil We now go further into the problem of the choice strategies for optimal parameters for the first kind node’s stencil. To optimize the 9-point scheme (2.9) for the first transitional node (see, Fig. 1(b)), we first perform classical dispersion analysis by assuming a plane-wave solution of the form U (x, y) = e−ik(x cos θ+y sin θ ) , where θ is the propagation angle from the y-axis. The following analysis is based on the assumption that the wavenumber k is a positive constant. Moreover, let v be the velocity of propagation, λ be the wavelength, and G be the number of gridpoints per wavelength, that is, G = λh . Since and k = ωv , we have kh = 2Gπ . Substituting Um,n := e−ik(xm cos θ+yn sin θ) into the Eq. (2.9), dividing both sides by the factor e−ik(xm cos θ+yn sin θ ) and finally applying the Euler formula eix = cos x + i sin x lead to the following dispersion equation
λ=
2πv
ω
2 [cos(kh cos θ ) − 1] a1 eikh sin θ + a2 + a3 e−ikh sin θ + [cos(kh sin θ ) − 1] b1 eikh cos θ + b2 + b3 e−ikh cos θ
+ k2 h2 c1 eikh(cos θ+sin θ) + c2 eikh sin θ + c3 e−ikh(cos θ −sin θ) + c4 eikh cos θ + c5 + c6 e−ikh cos θ + c7 e−ikh(− cos θ +sin θ) + c8 e−ikh sin θ + c9 e−ikh(cos θ+sin θ) = 0.
(2.10)
As stated in Proposition 2.1, to obtain a second order scheme, the parameters in the Eq. (2.9) should satisfy a1 + a2 + a3 = 1,
b1 + b2 + b3 = 1,
9
ci = 1,
(2.11)
i=1
a1 = a3 ,
b1 = b3 ,
−c1 + c3 − c4 + c6 − c7 + c9 = 0,
−
3 i =1
ci +
9
ci = 0.
(2.12)
i=7
Also, all the parameters of ghost nodes in Fig. 1(b) should be zero, that is T7 = T8 = T9 =
a3 + b 1
+ c7 k2 = 0,
h2 −2a3 + b2 h2 a3 + b 3 h2
+ c8 k2 = 0,
+ c9 k2 = 0.
(2.13) (2.14) (2.15)
T. Wu, Z. Chen / Journal of Computational and Applied Mathematics 267 (2014) 82–95
87
Combining the above Eqs (2.11)–(2.15) yields
a3 = a1 , a2 = 1 − 2a1 , b3 = b1 , b2 = 1 − 2b1 , c8 k2 h2 = 2(a1 + b1 ) − 1, c9 k2 h2 = −(a1 + b1 ), c1 = c3 + c5 + 2c6 + 2c8 + 4c9 − 1, c2 = −2c3 − c5 − 2c6 − c8 − 2c9 + 1, c4 = −c5 − c6 − 2c8 − 4c9 + 1, c7 = c9 .
(2.16)
From Eq. (2.16), we see that for the first kind node, the symmetry of the finite difference stencil is destroyed. Hence, different from the paper [3], the dispersion equation (2.10) has the imaginary part. Therefore, we cannot obtain the numerical phase velocity or the numerical group velocity any more. Hence, a proper method is needed to get optimal parameters for the first kind node. We next present a new way to obtain the optimal parameters of the finite difference equation (2.9) for the first transitional node. Firstly, we substitute the Eq. (2.16) into (2.10) to eliminate the parameters a2 , a3 , b2 , b3 , c1 , c2 , c4 , c7 , c8 , c9 . Then, under the assumption that all the parameters in the finite difference equation (2.9) are real numbers, we separate the dispersion equation (2.10) into the real part and the imaginary part. The real part of the dispersion equation (2.10) 2 [cos (kh cos θ ) cos (kh sin θ ) − cos (kh sin θ) − 2 cos (kh cos θ) + 2] (a1 + b1 )
+ k2 h2 2 cos (kh sin θ ) [cos (kh cos θ ) − 1] c3 + [cos (kh(cos θ + sin θ )) − cos (kh sin θ) − cos (kh cos θ ) + 1]c5 + 2 [cos (kh(cos θ + sin θ )) − cos (kh sin θ )] c6 = −4 cos (kh cos θ ) − 2 cos (kh sin θ ) + 4 + 2 cos (kh(cos θ + sin θ )) − k2 h2 [cos (kh(cos θ + sin θ )) − cos (kh sin θ) − cos (kh cos θ)] ,
(2.17)
and the corresponding imaginary part is
2 sin (kh sin θ) [cos (kh cos θ ) − 1] (a1 + b1 ) + k2 h2 2 sin (kh sin θ) [cos (kh cos θ) − 1]c3
+ [sin (kh (cos θ + sin θ)) − sin (kh sin θ ) − sin (kh cos θ)](c5 + 2c6 ) = [sin (kh (cos θ + sin θ )) − sin (kh sin θ ) − sin (kh cos θ)](2 + k2 h2 ).
(2.18)
Note that Eqs. (2.17) and (2.18) are functions of kh and θ . To reduce the numerical dispersion and improve the accuracy of the difference scheme for the first kind node, we propose the following rule. Rule 2.2 (Refined Choice Strategy). Step 1. Estimate kh, where k is the wavenumber for the first kind node and h is the step size for the fine grid. (m−1)2π ∈ Iθ , m = 1, 2, . . . , l, where l is a given positive integer. Step 2. Let Iθ := [0, 2π ), and choose θ = θm := l Step 3. Inserting above kh and θ to Eqs. (2.17) and (2.18) yields a linear system. Step 4. Applying the singular value decomposition method to solve the corresponding system yields the optimal weight parameters for Eqs. (2.17) and (2.18). Other parameters are obtained by Eq. (2.16). Different from the paper [3], here we let the number of gridpoints per wavelength G be a fixed number, which seems to be a proper choice, as in the above analysis we assume the wavenumber to be constant and we can obtain other parameters only by Eq. (2.16). In addition, the angle θ belongs to the interval [0, 2π ), which arises from that the symmetry of the finite difference stencil is destroyed. As an example, in Fig. 2, we present some optimal parameters for the first kind of transitional nodes. In this figure, we only show the parameters which appear in the Eqs. (2.17) and (2.18), as other parameters can be easily obtained by Eq. (2.16). We find that the parameters a1 , b1 , c3 vary smoothly, while the parameters c5 and c6 are very sensitive to the number of gridpoints per wavelength G, which may be the effect of the ghost nodes. 2.3. Optimal 9-point finite difference schemes for other nodes In this subsection, we firstly present optimal 9-point finite difference methods for the last four kinds of transitional nodes, which are also second order in accuracy. Then, based on the paper [3], we formulate optimal 9-point finite difference scheme for other nodes, which are not transitional nodes. For the second kind node, which also lie in the interior domain, we construct its consistent finite difference scheme with its nearest 8 coarse points, which is different from the first kind node. For the third kind of transitional nodes, we construct its consistent finite difference method with its nearest 8 symmetric nodes. The finite difference schemes for the fourth and fifth kind of transitional nodes are formulated similarly. We refer the interested readers to the Eqs (3.2) and (3.25) in [3] for the details. All these schemes are consistent with the Helmholtz equation with PML and are second order in accuracy (see, [3]). In addition, the parameters of the finite difference methods can be chosen with refined choice strategy (see, Rule 3.6 and Rule 3.8 in [3]).
88
T. Wu, Z. Chen / Journal of Computational and Applied Mathematics 267 (2014) 82–95
Fig. 2. Optimal coefficients for the first kind node.
a
b
Fig. 3. (a) A near-surface low-velocity layer model. (b) A low-velocity interlayer model.
For other nodes which are not transitional nodes, we formulate the optimal 9-point finite difference scheme for the Helmholtz equation with PML based on our previous work [3]. Finally, together with the work for the transitional nodes, we have a dispersion minimizing subgridding finite difference scheme for the 2D Helmholtz equation with PML, referred to as the refined optimal subgridding 9-point finite difference scheme for short, which is second order in accuracy and consistent with the equation. 3. Numerical experiments In this section, we present three numerical experiments to illustrate the efficiency of the schemes described in the last section. All the experiments are performed with Matlab 7v on an Intel Xeon (2-core) with 2.40 GHz and 2.32 Gb RAM. 3.1. A near-surface low-velocity layer model To validate the efficiency of the refined optimal subgridding 9-point finite difference scheme (the refined subgridding 9p) over the refined optimal 9-point finite difference scheme based on uniform grids proposed in [3] (the refined uniform 9p) for the case of variational wavenumbers, which often results from heterogeneous medium, we consider a near-surface low-velocity layer model shown in Fig. 3(a). The domain is [0 m, 1600 m] × [0 m, 1600 m] and there are four velocities: from the top, v := 1200 m/s, 2000 m/s, 2400 m/s and 2800 m/s. A point source δ(x − xs , y − ys ) R(f , f0 ) is located at the point (xs , ys ) = (0 m, 800 m), where R(f , f0 ) is the Ricker wavelet in the frequency domain defined by
R(f , f0 ) :=
+∞ −∞
(1 − 2π 2 f02 t 2 ) exp(−π 2 f02 t 2 )e−i2π ft dt ,
T. Wu, Z. Chen / Journal of Computational and Applied Mathematics 267 (2014) 82–95
89
whose dominant frequency is f0 = 25 Hz. In addition, all the receivers are placed on the top of the ground, that is, they are in the line of y = 0. In the following computation, the time sampling is 8 ms and the highest frequency we compute is 60 Hz. For the refined uniform 9p, the uniform step size hu := ∆x = ∆y is used for both variables x and y. For the refined subgridding 9p, H and h denote the coarse and fine step size, respectively. In this experiment, we let H = 8 m and h = 4 m. In detail, we only use fine grids for the top layer whose velocity is v = 1200 m/s. We next generate common-shot-point records (shot profiles) for the wave field by the refined uniform 9p and the refined subgridding 9p. The monofrequency wave field (real part) for f = 30 Hz obtained by the refined subgridding 9p with (H , h) = (8 m, 4 m) is present in Fig. 4(a). It is easy to see that PML’s absorption of the outgoing waves is efficient, and almost no boundary reflections exist. Additionally, the downward incident waves and transmissive waves are clear. In the first velocity layer, the incident waves are interfered with the reflected waves returned from the bottom of the first velocity layer, and this phenomena obeys Snell’s law. In Fig. 4(b)–(d), we show the common-shot-point records generated by the refined uniform 9p with uniform grid size hu = 4 m, the refined uniform 9p with hu = 8 m and the refined subgridding 9p with (H , h) = (8 m, 4 m), respectively. In these three pictures, from the top, the first wave represents the direct wave, the second one is the wave reflected by the bottom of the first layer (the low-velocity layer), the third one denotes the wave reflected by the bottom of the second layer, and the fourth wave is the wave reflected by the bottom of the third layer (see, the reference number in Fig. 4(d)). Comparing these three figures, we find that Fig. 4(d) is almost as clear as Fig. 4(b), while there are some nonphysical oscillations in Fig. 4(c), which is caused by the numerical dispersion. To see clearly, Fig. 4(e) and Fig. 4(f) present local enlargement for Fig. 4(c) and Fig. 4(d), respectively. From the last two pictures in Fig. 4, some nonphysical oscillations in the background, such as in the time interval [0.04 s, 0.40 s], can be found in the former, but these do not exist in the latter, which verify that the refined subgridding 9p cannot only improve the accuracy but also suppress nonphysical oscillations. Furthermore, when comparing the cost of Fig. 4(d) with that of Fig. 4(b), there is a 56.44% decrease in memory and a 76.63% decrease in the computational time. Therefore, the efficiency of the refined subgridding 9p is confirmed. 3.2. A low-velocity interlayer model To validate the refined subgridding 9p’s availability for the case of variational wavenumbers, we also consider a lowvelocity interlayer model shown in Fig. 3(b). The domain is [0 m, 1600 m] × [0 m, 1600 m] and there are four velocities: from the top, v := 2000 m/s, 1200 m/s, 2400 m/s and 2800 m/s. A point source, which is the same as the one used in the above example, is located at the point (xs , ys ) = (0 m, 800 m). Also, all the receivers are placed on the top of the ground, the time sampling is 8 ms and the highest frequency we compute is 60 Hz. For the refined subgridding 9p, we let the coarse grid size H = 8 m and the fine grid size h = 4 m. In detail, we only use fine grids for the interlayer whose velocity is v = 1200 m/s. In Fig. 5, we generate common-shot-point records for the wave field by the refined uniform 9p and the refined subgridding 9p. The monofrequency wave field (real part) for f = 35 Hz obtained by the refined subgridding 9p with (H , h) = (8 m, 4 m) is present in Fig. 5(a). It is easy to see that PML’s absorption of the outgoing waves is efficient, and almost no boundary reflections exist. Additionally, the downward incident waves and transmissive waves are clear. In the second velocity layer, the incident waves are interfered with the reflected waves returned from the bottom of the second velocity layer, and this phenomena obeys Snell’s law. In Fig. 5(b)–(d), we show the common-shot-point records generated by the refined uniform 9p with uniform grid size hu = 4 m, the refined uniform 9p with hu = 8 m and the refined subgridding 9p with (H , h) = (8 m, 4 m), respectively. In these three pictures, from the top, the first wave represents the direct wave, the second one is the wave reflected by the bottom of the first layer, the third one denotes the wave reflected by the bottom of the second layer (the low-velocity layer), and the fourth wave is the wave reflected by the bottom of the third layer (see, the reference number in Fig. 5(d)). Comparing these three figures, we find that Fig. 5(d) is almost as clear as Fig. 5(b), while there are some nonphysical oscillations in Fig. 5(c), which is caused by the numerical dispersion. To see clearly, Fig. 5(e) and Fig. 5(f) present local enlargement for Fig. 5(c) and Fig. 5(d), respectively. From the last two pictures in Fig. 5, some nonphysical oscillations in the background, such as in the time interval [0.88 s, 1.20 s], can be found in the former, but these do not exist in the latter, which verify that the refined subgridding 9p cannot only improve the accuracy but also suppress nonphysical oscillations. Furthermore, when comparing the cost of Fig. 5(d) with that of Fig. 5(b), there is a 54.31% decrease in memory and a 77.25% decrease in the computational time. Therefore, the efficiency of the refined subgridding 9p is confirmed. 3.3. A two-layered model To compare the refined subgridding 9p with the refined uniform 9p in a more quantitative way, we consider a two-layered model [23]. For the observation point x = (x, y) and source point xs = (xs , ys ) with ys < 0, we consider the solution G(x, xs ) of the Helmholtz equation in a two-layered medium in R2 which satisfies
∆G(x, xs ) + k2 (x)G(x, xs ) = −δ(x − xs ),
(3.1)
with continuity conditions G(x, xs )|y=0+ = G(x, xs )|y=0− ,
(3.2)
90
T. Wu, Z. Chen / Journal of Computational and Applied Mathematics 267 (2014) 82–95
a
b
c
d
e
f
Fig. 4. (a) The real part of the solution for f = 30 Hz obtained by the refined subgridding 9p with (H , h) = (8 m, 4 m). (b) The common-shot-point records for the wave field by the refined uniform 9p with hu = 4 m. (c) The common-shot-point records for the wave field by the refined uniform 9p with hu = 8 m. (d) The common-shot-point records for the wave field by the refined subgridding 9p with (H , h) = (8 m, 4 m). (e) Local enlargement for (c). (f) Local enlargement for (d).
∂ G(x, xs ) ∂ G(x, xs ) = , ∂ y y=0+ ∂ y y=0−
(3.3)
where the wavenumber k(x) =
k1 , k2 ,
for y > 0, for y < 0,
in which k1 and k2 are given positive numbers satisfying k1 ̸= k2 . Without loss of generality, we assume that k1 < k2 holds for this problem.
T. Wu, Z. Chen / Journal of Computational and Applied Mathematics 267 (2014) 82–95
a
b
c
d
e
f
91
Fig. 5. (a) The real part of the solution for f = 35 Hz obtained by the refined subgridding 9p with (H , h) = (8 m, 4 m). (b) The common-shot-point records for the wave field by the refined uniform 9p with hu = 4 m. (c) The common-shot-point records for the wave field by the refined uniform 9p with hu = 8 m. (d) The common-shot-point records for the wave field by the refined subgridding 9p with (H , h) = (8 m, 4 m). (e) Local enlargement for (c). (f) Local enlargement for (d).
We next present the analytical formula of the solution for the above two-layered problem. It follows from the Fourier transform that the solution is given by [23]
• y < 0, G(x, xs ) = Φ (x, xs ) +
i 4π
+∞
g1 (x, xs , ξ )dξ , −∞
(3.4)
92
T. Wu, Z. Chen / Journal of Computational and Applied Mathematics 267 (2014) 82–95
Table 1 The performance of the refined uniform 9p. The pair (·, ·) denotes the error for the layer with k(x) = 0.1571 and the layer with k(x) = 0.0628, respectively. hu
10 m
5m
2.5 m
1.25 m
n Errors Time
8281 (19.5045, 10.0745) 4.5708 s
25 921 (5.4717, 2.4063) 50.7783 s
68 121 (1.8233, 0.6844) 326.9001 s
231 361 (0.6093,0.2453) 3471.7 s
Table 2 The performance of the refined subgridding 9p. The pair (·, ·) denotes the error for the layer with k(x) = 0.1571 and the layer with k(x) = 0.0628, respectively.
(H , h)
(10 m, 5 m)
(5 m, 2.5 m)
(2.5 m, 1.25 m)
(1.25 m, 0.625 m)
n Errors Time
10 286 (11.8947, 4.8360) 7.0356 s
32 631 (3.6300, 1.2230) 80.4185 s
87 541 (1.1214, 0.3889) 553.7567 s
300 282 (0.4225, 0.1705) 5851.8 s
• y > 0, G(x, xs ) =
i 2π
+∞
g2 (x, xs , ξ )dξ ,
(3.5)
−∞
where
Φ (x, xs ) =
i 4
(1)
H0 (k2 |x − xs |) ,
(3.6)
(1)
in which H0 is the Hankel function of the first kind with order zero, i2 = −1,
2 − k2 − √2 2 ξ ξ 2 − k21 2 ξ −k2 (y+ys ) iξ (x−xs ) e , e 2 2 2 2 2 2 i ξ − k2 ξ − k1 + ξ − k2 √2 2 k22 − ξ 2 − i ξ 2 − k21 e−i k2 −ξ (y+ys ) eiξ (x−xs ) , g1 (x, xs , ξ ) = k22 − ξ 2 i ξ 2 − k21 + k22 − ξ 2 √2 2 k22 − ξ 2 − k21 − ξ 2 e−i k2 −ξ (y+ys ) eiξ (x−xs ) , k22 − ξ 2 k21 − ξ 2 + k22 − ξ 2 and
−√ξ 2 −k2 y+√ξ 2 −k2 y 1 2 s e eiξ (x−xs ) , i ξ 2 − k21 + ξ 2 − k22 √ √ − ξ 2 −k21 y−i k22 −ξ 2 ys e g2 (x, xs , ξ ) = eiξ (x−xs ) , 2 2 i ξ 2 − k1 + k2 − ξ 2 √ √ i k21 −ξ 2 y− k22 −ξ 2 ys e eiξ (x−xs ) , 2 2 2 2 k1 − ξ + k2 − ξ
for |ξ | > k2 ,
for k1 < |ξ | < k2 ,
for |ξ | < k1 ,
for |ξ | > k2 ,
for k1 < |ξ | < k2 ,
for |ξ | < k1 .
Then, we consider a two-layered problem (3.1) with the source point fixed at (250 m, −52 m), k1 = 0.0628 and k2 = 0.1571. The bounded domain [0 m, 500 m] × [−52 m, 448 m] is our interested area. The corresponding solution is given visa (3.4)–(3.6). Thus, numerical integrals have to be dealt with to obtain the solution. In this experiment, we apply the Gauss-type quadratures [24] for weakly singular integrals appearing in (3.4) and (3.5). This kind of quadrature rules possess polynomial order of accuracy, and enable us to obtain the solution with high accuracy. In the following, we take the solution obtained by Gauss-type quadratures as the analytical solution of the two-layered problem. We next compare the performance of the refined uniform 9p and the refined subgridding 9p, including errors, computational times, number of unknowns and the errors’ decay exponents (briefly, decay exp.). The error is measured
T. Wu, Z. Chen / Journal of Computational and Applied Mathematics 267 (2014) 82–95
93
Table 3 The decay exponent for the refined subgridding 9p. n
10 286
32 631
87 541
300 282
Errors Decay exp.
11.8947
3.6300 1.0280
1.1214 1.1903
0.4225 0.7920
a
b
c
d
Fig. 6. (a) The real part of the analytical solution. (b) The real part of the solution generated by the refined uniform 9p with hu = 1.25 m. (c) The real part of the solution generated by the refined subgridding 9p with (H , h) = (2.5 m, 1.25 m). (d) The real part of solutions in the line of x = 125 m generated by different methods.
in C-norm, which is defined as: for any complex vector z = [z1 , z2 , . . . , zM ],
∥z∥C := max |zj |, 1≤j≤M
where |zj | is the complex modulus of zj . Let n represent the number of all unknowns in the two dimensional computational domain in one test, that is, n is the sum of the number of the unknowns in the interior area and that of the unknowns in the ln(en /em ) PML. Each decay exponent is computed by ln(m , where en denotes the numerical solution’s error for one test in which the /n) number of all unknowns is n. Also, for the refined uniform 9p, the uniform step size hu := ∆x = ∆y is used for both variables x and y. For the refined subgridding 9p, H and h denote the coarse and fine step size, respectively. In this experiment, the coarse grids are refined only once in the layer with bigger wavenumber k(x) = 0.1571. Therefore, h = H /2 holds. In Tables 1 and 2, we present the number of all unknowns in the two dimensional computational domain, errors and computational time for the refined uniform 9p and subgridding 9p, respectively. Seen from Tables 1 and 2, for the two methods, the error of the layer with bigger wavenumber is almost two times larger than that of the layer with smaller wavenumber, which arises from the varying of the wavenumber in the two layers. In addition, for the refined subgridding 9p with a pair of given step sizes (H , h), its corresponding error is almost half of that obtained by the refined uniform 9p
94
T. Wu, Z. Chen / Journal of Computational and Applied Mathematics 267 (2014) 82–95
with hu = H, and is almost twice as big as that obtained by the refined uniform 9p with hu = h. However, we find that there is a sharp decrease in both the number of unknowns and the computational time for the refined subgridding 9p with (H , h), when compared with the refined uniform 9p with hu = h. In Table 3, we present the decay exponent for the refined subgridding 9p. We find that the decay exponent is one. As n denotes the number of all unknowns in the two dimensional computational domain, the decay exponent being one indicates that the overall global convergence for the refined subgridding 9p is second order. For further comparison, in Fig. 6(a)–(c) we present the real part of the solution for the two-layered problem in the domain [0 m, 500 m] × [−52 m, 448 m] by the analytical method, the refined uniform 9p with hu = 1.25 m and the refined subgridding 9p with (H , h) = (2.5 m, 1.25 m), respectively. In each of these three pictures, we see that the downward incident waves and transmissive waves are clear. Moreover, in the first layer, the incident waves are interfered with the reflected waves returned from the bottom of the first layer, and this phenomena obeys Snell’s law. It is easy to see that both Fig. 6(b) and Fig. 6(c) are similar with Fig. 6(a), but Fig. 6(c) needs less memory and computational time than those needed by Fig. 6(b) (see, Tables 1 and 2), which verifies the efficiency of the refined subgridding 9p. Also, almost no boundary reflections exist in Fig. 6(b) and Fig. 6(c), which confirms the efficiency of the PML’s absorption. Finally, Fig. 6(d) presents the real part of the two-layered problem’s solution in the line of x = 125 m respectively by the analytical method, the refined uniform 9p with hu = 2.5 m, the refined uniform 9p with hu = 1.25 m, the refined subgridding 9p with (H , h) = (2.5 m, 1.25 m) and the refined subgridding 9p with (H , h) = (1.25 m, 0.625 m). From this figure, we find that the solution obtained by the refined subgridding 9p with (H , h) = (1.25 m, 0.625 m) is a good approximation for the analytical solution. Furthermore, comparing the two solutions obtained by the refined subgridding 9p with different step sizes, we see that the solution obtained by the refined subgridding 9p converges to the analytical solution. In addition, for the results obtained by the refined subgridding 9p, there are nonphysical oscillations near the line y = 0 m, which is the interface between the coarse grids and the fine grids, and this phenomena indicates the validity of our refined choice strategy for the transitional nodes. Now, we come to our conclusion that, when balancing the accuracy against the computational cost, the refined subgridding 9p possesses much improvement over the refined uniform 9p, especially for large-scale layered problems. 4. Conclusions We summarize and comment on the numerical results. Firstly, we applied subgrids to discretize the computational domain, and presented a consistent 9-point subgridding difference scheme for the Helmholtz equation with PML, which is second order in accuracy. Then, for the transitional nodes, we proposed refined choice strategies for choosing optimal parameters based on minimizing the numerical dispersion. Finally, numerical experiments were presented to confirm that the refined subgridding 9p is a good choice for the Helmholtz equation with varying wavenumbers, especially for large-scale layered problems. References [1] R.G. Pratt, C. Shin, G.J. Hicks, Gauss–Newton and full Newton methods in frequency-space seismic waveform inversion, Geophys. J. Int. 133 (1998) 341–362. [2] R.G. Pratt, M.H. Worthington, Inverse theory applied to multi-source cross-hole tompgraphy. Part 1: acoustic wave-equation method, Geophys. Prospect. 38 (1990) 287–310. [3] Z. Chen, D. Cheng, W. Feng, T. Wu, An optimal 9-point finite difference scheme for the Helmholtz equation with PML, Int. J. Numer. Anal. Model. 10 (2013) 389–410. [4] J.W. Kang, L.F. Kallivokas, Mixed unsplit-field perfectly matched layers for transient simulations of scalar waves in heterogeneous domains, Comput. Geosci. 14 (2010) 623–648. [5] J.W. Kang, L.F. Kallivokas, The inverse medium problem in heterogeneous PML-truncated domains using scalar probing waves, Comput. Methods Appl. Mech. Engrg. 200 (2011) 265–283. [6] Y.A. Erlangga, C.W. Oosterlee, C. Vuik, A novel multigrid based preconditioner for heterogeneous Helmholtz problems, SIAM J. Sci. Comput. 27 (2006) 1471–1492. [7] 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. [8] Z. Chen, D. Cheng, W. Feng, T. Wu, H. Yang, A new multigrid based preconditioned Krylov subspace method for the Helmholtz equation with PML, J. Math. Anal. Appl. 383 (2011) 522–540. [9] Z. Chen, D. Cheng, T. Wu, A dispersion minimizing finite difference scheme and preconditioned solver for the 3D Helmholtz equation, J. Comput. Phys. 231 (2012) 8152–8175. [10] 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. [11] C.-H. Jo, C. Shin, J.H. Suh, An optimal 9-point, finite-difference, frequency-space, 2-D scalar wave extrapolator, Geophysics 61 (1996) 529–537. [12] C. Shin, H. Sohn, A frequency-space 2-D scalar wave extrapolator using extended 25-point finite-difference operator, Geophysics 63 (1998) 289–296. [13] C. Chu, P.L. Stoffa, Nonuniform grid implicit spatial finite difference method for acoustic wave modeling in tilted transversely isotropic media, J. Appl. Geophys. 76 (2012) 44–49. [14] Y. Li, L. Sun, W. Hong, Helmholtz equation sub-grid method for multi-transition region, National Microwave Meeting, Qingdao China (2011) pp. 1150–1152. [15] P. Moczo, E. Bystrický, J. Kristek, J.M. Carcione, M. Bouchon, Hybrid modeling of P-SV seismic motion at inhomogeneous viscoelastic topographic structures, Bull. Seismol. Soc. Am. 87 (1997) 1305–1323. [16] J.W. Nehrbass, R. Lee, Optimal finite-difference sub-gridding techniques applied to the Helmholtz equation, IEEE Trans. Microw. Theory Tech. 48 (2000) 976–984. [17] B. Donderici, F.L. Teixeira, Improved FDTD subgridding algorithms via digital filtering and domain overriding, IEEE Trans. Antennas and Propagation 53 (2005) 2938–2951.
T. Wu, Z. Chen / Journal of Computational and Applied Mathematics 267 (2014) 82–95
95
[18] J. Kristek, P. Moczo, M. Galis, Stable discontinuous staggered grid in the finite difference modelling of seismic motion, Geophys. J. Int. 183 (2010) 1401–1407. [19] D.T. Prescott, N.V. Shuley, A method for incorporating different sized cells into the finite-difference time-domain analysis technique, IEEE Microw. Guided Wave Lett. 2 (1992) 434–436. [20] A. Pitarka, 3D elastic finite-difference modeling of seismic motion using staggered grids with nonuniform spacing, Bull. Seismol. Soc. Am. 89 (1999) 54–68. [21] S.A.M. Oliveira, A fourth-order finite-difference method for the acoustic wave equation on irregular grids, Geophysics 68 (2003) 672–676. [22] C. Jastram, E. Tessmer, Elastic modelling on a grid with vertically varying spacing, Geophys. Prospect. 42 (1994) 357–370. [23] J. Coyle, Locating the support of objects contained in a two-layered background medium in two dimensions, Inverse Problems 16 (2000) 275–292. [24] H. Kaneko, Y. Xu, Gauss-type quadratures for weakly singular integrals and their application to Fredholm integral equations of the second kind, Math. Comput. 62 (1994) 739–753.