# Optimal 25-Point Finite-Difference Subgridding Techniques for the 2D Helmholtz Equation

Authors: Tingting Wu; Zhongying Chen; Jian Chen
Journal: Mathematical Problems in Engineering
 (2016)

---

## Abstract

We present an optimal 25-point finite-difference subgridding scheme for solving the 2D Helmholtz equation with perfectly matched layer (PML). This scheme is second order in accuracy and pointwise consistent with the equation. Subgrids are used to discretize the computational domain, including the interior domain and the PML. For the transitional node in the interior domain, the finite difference equation is formulated with ghost nodes, and its weight parameters are chosen by a refined choice strategy based on minimizing the numerical dispersion. Numerical experiments are given to illustrate that the newly proposed schemes can produce highly accurate seismic modeling results with enhanced efficiency.


---

## Body

## 1. Introduction

In a Cartesian coordinate system, the 2D Helmholtz equation is(1)Δu+k2u=g,onR2,where Δ is the Laplacian, u usually represents a pressure field, k∈R, called wavenumber, is a positive number, and g denotes the Fourier transformation of the source function. The above equation is indeed the frequency domain wave equation, as it can be obtained by applying the Fourier transform with respect to time to the 2D acoustic wave equation. The Helmholtz equation has important applications in electromagnetics, acoustics, quantum mechanics, and geophysics. In geophysics, wave equation modeling in the frequency domain has many advantages over time domain modeling [1, 2]. For instance, to perform wave equation inversion and tomography, only a few frequency components are required for certain geometries. Each frequency can be computed independently, which favors parallel computing. Furthermore, frequency domain modeling allows accurate modeling of attenuation effects and allows optimal spatial discretization intervals to be chosen at each frequency. The Helmholtz equation is so important that its numerical simulation has always been under active research (cf. [2–5]).Due to the computer’s finite memory and speed limitation, the Helmholtz problem should be solved in a finite domain. Applying the PML technique (see [6, 7]) to truncate the infinite domain into a bounded rectangular domain leads to the equation:(2)∂∂xA∂u∂x+∂∂yB∂u∂y+Ck2u=g~,where(3)g~≔g,inside PML,0,outside PML,A≔eyex,B≔exey,C≔exey, in which ex≔1+iσx/ω,  ey≔1+iσy/ω,  ω is the angular frequency, and σx is a function only of x defined as(4)σx≔2πa0f0lxLPML2,inside PML,0,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. Equation (2) can be seen as a general form of the Helmholtz equation (1) with the corresponding PML equation, since in the interior domain, A=B=C=1. We call it the Helmholtz equation with PML. Moreover, (2) is also valid for variable k(x,y) (see [8, 9]).The finite difference scheme is a popular method in certain engineering fields such as geophysics as it is both easy to implement and computationally efficient (cf. [1, 2, 4, 5, 10–13]). However, it is still a difficult problem for the finite difference method 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. Using uniform structured grids always leads to unnecessary extra computations, while using unstructured grids, adaptive to velocity variations, does not always help speed up the computations, as implementation of unstructured grids consumes significant amount of computer memory. A simple and natural alternative to this problem is to use structured nonuniform grids.The subgridding algorithm is a kind of structured nonuniform grids, which usually introduce locally refined meshes into the base mesh (cf. [14–17]). As subgridding refines the mesh locally, we have major computational savings. 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. To alleviate these problems, considerable effort has been devoted (see [14–16, 18–21] and the references therein). For example, based on minimizing the numerical dispersion, Nehrbass and Lee proposed a systematic method, which optimally chooses the finite difference coefficients at the interface between the coarse 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. 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 complicated, which may not be a good choice for practical cases, and the PML technique is not considered. 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 is easier. However, we can not extend this method to the Helmholtz equation with PML, as it will not possess the property of consistency (see [22–24]). In addition, Kristek et al. proposed downsampling the wave fields in fine grid blocks for better numerical stability (see [19]). For subgridding methods in the frequency domain, there does not exist the stability problem.A consistent 25-point finite difference scheme for solving the 2D Helmholtz equation with PML was developed in [24]. This scheme was based on uniform grids and second order in accuracy. In detail, by a weighted average method, a symmetric discretization of the Laplacian term and that of the term of zero order in the Helmholtz equation were constructed, respectively. A consistent 25-point finite difference scheme for the 2D Helmholtz equation with PML was then developed by combining the approximation for the Laplacian term with that for the term of zero order. This scheme had twelve weight parameters, which could be chosen properly to improve the accuracy of the solution of the Helmholtz equation. To choose optimal coefficients for the finite difference scheme, [24] computed its normalized phase velocity and minimized the error between the normalized phase velocity and one. The improvement of the accuracy and the reduction of the numerical dispersion are significant. However, applying the optimal 25-point finite difference scheme proposed in [24] to deal with large-scale layered problems usually leads to unnecessary extra computations, as the grid size is established by the slowest velocity. To reduce the computational cost and improve the accuracy further, the main purpose of this paper is to develop a consistent 25-point finite-difference subgridding scheme for the Helmholtz equation with PML which enjoys the second order accuracy. To this end, subgrids will be used to discretize the computational domain, including the interior domain and the PML. For the transitional node in the interior domain, the weighted average method will also be used to formulate its finite difference equation, but more weight parameters will be introduced, when compared to the method proposed in [24]. In detail, the finite difference stencil for the transitional node will not be symmetric. As the finite difference stencil for the transitional node is not symmetric, the normalized phase velocity can not be obtained anymore. Therefore, we will propose an approach for choosing optimal coefficients of the finite difference stencil for the transitional node. 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 will address this issue in detail in this paper.This paper is organized as follows. In Section2, we apply the subgridding technique to discretize the computational domain. For the transitional node in the interior domain, we employ the weighted average method to formulate its finite difference scheme with ghost nodes and then prove that this scheme is consistent with the Helmholtz equation and is second order in accuracy. To choose optimal weight parameters of the scheme, a new refined choice strategy is also proposed. In Section 3, numerical experiments are given to demonstrate the efficiency of the scheme. These numerical results show that the proposed schemes can produce highly accurate seismic modeling results with enhanced efficiency. Finally, Section 4 contains the conclusions of this paper.


## 2. An Optimal 25-Point Finite-Difference Subgridding Scheme for the Helmholtz Equation with PML

In this section, subgrids are utilized to discretize the computational domain, and an optimal 25-point finite-difference subgridding scheme for the 2D Helmholtz equation with PML is then developed based on the following two aspects: firstly, as the domain is divided into different blocks, the optimal 25-point finite difference scheme proposed in [24], which is indeed a dispersion minimizing method, is directly applied in each block; secondly, for transitional nodes in the interior domain, we construct a consistent 25-point finite difference scheme with ghost nodes and propose a new refined optimization rule for choosing the scheme’s weight parameters based on minimizing the numerical dispersion.To illustrate the subgridding technique for discretizing the computation domain, we take a two-layered model as an example, which is presented in Figure1(a). In this picture, the mesh established by the red lines is the base mesh, and the green lines are added to refine the mesh locally. The domain surrounded by the blue rectangle is our area of interest (the interior area), while the rest is the perfectly matched layer. As shown in Figure 1(a), the transitional grids are categorized into ten groups. Consistent finite difference methods for them will be given one by one.Figure 1
(a) Subgridding method for a two-layered model, (b) the first kind node, (c) the second kind node, (d) the third kind node, and (e) the fifth kind node.
(a)
(b)
(c)
(d)
(e)


### 2.1. A Consistent 25-Point Finite Difference Scheme for the Second Kind Node

We next present a consistent 25-point finite difference method for the second kind of transitional nodes, which is second order in accuracy. Seen from Figure1(a), we find that the second kind of transitional nodes lie in the interior domain, in which the Helmholtz equation with PML (2) deteriorates to be the Helmholtz equation. With eight ghost nodes, we will construct a 25-point stencil for the second kind node. This stencil is shown in Figure 1(c), where the nodes in the bottom red line are ghost nodes.We turn to constructing a 25-point finite difference scheme for the second kind node by the weighted averaged method (see [22–24] and the references therein). For simplicity, we first introduce some notations. 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. For each m and n, we set um,n≔ux=xm,y=yn and km,n≔kx=xm,y=yn. Next, different from [24], more weight parameters will be introduced to relax the restriction for the symmetry of the finite difference stencil. In detail, let a1,a2,…,a8,  b1,b2,…,b8,  c1,c2,…,c25 be constants to be determined in the following. With these weight parameters, we then establish an approximation for the Laplacian term by the weighted average method. Let(5)Lh,xum,n+j≔um+1,n+j+um-1,n+j-2um,n+jh2,for j∈Z2≔{-1,0,1}, and(6)L2h,xum,n+j≔um+2,n+j+um-2,n+j-2um,n+j2h2,for j∈Z3≔{-2,-1,0,1,2}. Upon introducing two notations(7)L~h,xux=xm,y=yn≔a1Lh,xum,n-1+a2Lh,xum,n+a3Lh,xum,n+1,L~2h,xux=xm,y=yn≔a4L2h,xum,n-2+a5L2h,xum,n-1+a6L2h,xum,n+a7L2h,xum,n+1+a8L2h,xum,n+2,we approximate ∂2u/∂x2 as (8)∂2u∂x2x=xm,y=yn≈L~h,xux=xm,y=yn+L~2h,xux=xm,y=yn. Similar treatment with respect to ∂2u/∂y2 leads to(9)∂2u∂y2x=xm,y=yn≈L~h,yux=xm,y=yn+L~2h,yux=xm,y=yn, where(10)L~h,yux=xm,y=yn≔b1Lh,yum-1,n+b2Lh,yum,n+b3Lh,yum+1,n,L~2h,yux=xm,y=yn≔b4L2h,yum-2,n+b5L2h,yum-1,n+b6L2h,yum,n+b7L2h,yum+1,n+b8L2h,yum+2,n.By setting L~h≔L~h,x+L~h,y+L~2h,x+L~2h,y, we construct an approximation for the Laplacian term as(11)Δux=xm,y=yn≈L~hux=xm,y=yn. We then approximate k2ux=xm,y=yn by a weighted average: (12)k2ux=xm,y=yn≈I~hk2ux=xm,y=yn, where(13)I~hk2ux=xm,y=yn≔∑l=04∑j=15cj+5lk2um-3+j,n-2+l.We now establish a 25-point difference approximation for the Helmholtz equation as(14)L~hux=xm,y=yn+I~hk2ux=xm,y=yn=0.We next consider the consistency of the 25-point difference scheme (14) and give a convergence analysis. The results are presented in the following theorem.Theorem 1.
(i) If∑i=18ai=1, ∑i=18bi=1, and ∑i=125ci=1, then the 25-point finite difference approximation (14) is pointwise consistent with the Helmholtz equation (1) and is a first order scheme.
(ii) If the parameters in (14) satisfy(15)∑i=18ai=1,∑i=18bi=1,∑i=125ci=1,a3+a7+2a8-a1-2a4-a5=0,b3+b7+2b8-b1-2b4-b5=0,2∑i=15c20+i-ci+∑i=610c10+i-ci=0,∑i=04c4+5i+2c5+5i-2c1+5i-c2+5i=0,then the 25-point finite difference approximation (14) is pointwise consistent with the Helmholtz equation (1) and is a second order scheme.Proof.
Assume thatxm≤x<xm+1 and yn≤y<yn+1. It follows from the Taylor theorem that(16)L~hux=xm,y=yn=∑i=18ai∂2u∂x2+∑i=18bi∂2u∂y2+ν1h+Oh2,I~hk2ux=xm,y=yn=∑i=125cik2u+ν2h+Oh2,where(17)ν1≔a3-a1+2a8-a4+a7-a5∂3u∂x2∂y+b3-b1+2b8-b4+b7-b5∂3u∂x∂y2,ν2≔∑i=04c4+5i+2c5+5i-2c1+5i-c2+5i∂ku∂x+2∑i=15c20+i-ci+∑i=610c10+i-ci∂ku∂y. Therefore, when ∑i=18ai=1, ∑i=18bi=1, and ∑i=125ci=1 hold, we have that the left hand side of the 25-point finite difference approximation (14) is equivalent to(18)∂2u∂x2+∂2u∂y2+k2u+ηh+Oh2,where η≔ν1+ν2. Together with the concept of consistency (see [13]), we conclude the results of this theorem.We now turn to constructing a second order finite difference scheme for the second kind node in Figure1(c), which is also consistent with the Helmholtz equation. We begin with rewriting the finite difference scheme (14). Substituting (7), (10), and (13) into (14) and replacing um+i,n+j with Um+i,n+j  (i,j∈Z2) yield the 25-point finite difference equation(19)T1Um-2,n-2+T2Um-1,n-2+T3Um,n-2+T4Um+1,n-2+T5Um+2,n-2+T6Um-2,n-1+T7Um-1,n-1+T8Um,n-1+T9Um+1,n-2+T10Um+2,n-1+T11Um-2,n+T12Um-1,n+T13Um,n+T14Um+1,n+T15Um+2,n+T16Um-2,n+1+T17Um-1,n+1+T18Um,n+1+T19Um+1,n+1+T20Um+2,n+1+T21Um-2,n+2+T22Um-1,n+2+T23Um,n+2+T24Um+1,n+2+T25Um+2,n+2=0,where the coefficients are given by(20)T1≔a4+b44h2+c1km-2,n-22,T2≔b54h2+c2km-1,n-22,T3≔-2a4+b64h2+c3km,n-22,T4≔b74h2+c4km+1,n-22,T5≔a4+b84h2+c5km+2,n-22,T6≔a54h2+c6km-2,n-12,T7≔a1+b1h2+c7km-1,n-12,T8≔-2a1+b2-1/2a5h2+c8km,n-12,T9≔a1+b3h2+c9km+1,n-12,T10≔a54h2+c10km+2,n-12,T11≔a6-2b44h2+c11km-2,n2,T12≔a2-2b1-1/2b5h2+c12km-1,n2,T13≔-4a2+b2-a6+b62h2+c13km,n2,T14≔a2-2b3-1/2b7h2+c14km+1,n2,T15≔a6-2b84h2+c15km+2,n2,T16≔a74h2+c16km-2,n+12,T17≔a3+b1h2+c17km-1,n+12,T18≔-2a3+b2-1/2a7h2+c18km,n+12,T19≔a3+b3h2+c19km+1,n+12,T20≔a74h2+c20km+2,n+12,T21≔a8+b44h2+c21km-2,n+22,T22≔b54h2+c22km-1,n+22,T23≔-2a8+b64h2+c23km,n+22,T24≔b74h2+c24km+1,n+22,T25≔a8+b84h2+c25km+2,n+22.To guarantee the above scheme to be second order in accuracy, the parameters in (19) should satisfy the conditions proposed in Theorem 1. On the other hand, as the ghost nodes do not exist in the stencil for the second kind node, we let their coefficients be zero; that is, T16=0,  T17=0,  T18=0,  T19=0,  T20=0,  T21=0,  T23=0, and T25=0. This indicates that the parameters in the scheme (19) may vary according to the wavenumber k. By far, a second order finite difference scheme for the second kind node is formulated. Choice strategies for optimal parameters of the finite difference scheme will be presented in the next subsection.


### 2.2. Choice Strategies for Optimal Parameters for the Second Kind Node’s Stencil

We now go further into the problem of the choice strategies for optimal parameters for the second kind node’s stencil (see Figure1(c)). To optimize parameters for the second kind node’s stencil, we first perform classical dispersion analysis by assuming a plane-wave solution of the form U(x,y)=e-ik(xcos⁡θ+ysin⁡θ), 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.SubstitutingUm,n≔e-ik(xmcos⁡θ+ynsin⁡θ) into (19), dividing both sides by the factor e-ik(xmcos⁡θ+ynsin⁡θ), and finally applying the Euler formula eix=cos⁡x+isin⁡x lead to the following dispersion equation:(21)a4ei2khsin⁡θ+a5eikhsin⁡θ+a6+a7e-ikhsin⁡θ+a8e-i2khsin⁡θP2+b4ei2khcos⁡θ+b5eikhcos⁡θ+b6+b7e-ikhcos⁡θ+b8e-i2khcos⁡θQ2+4a1eikhsin⁡θ+a2+a3e-ikhsin⁡θP1+4b1eikhcos⁡θ+b2+b3e-ikhcos⁡θQ1+2k2h2c1ei2khcos⁡θ+c2eikhcos⁡θ+c3+c4e-ikhcos⁡θ+c5e-i2khcos⁡θei2khsin⁡θ+c6ei2khcos⁡θ+c7eikhcos⁡θ+c8+c9e-ikhcos⁡θ+c10e-i2khcos⁡θeikhsin⁡θ+c11ei2khcos⁡θ+c12eikhcos⁡θ+c13+c14e-ikhcos⁡θ+c15e-i2khcos⁡θ+c16ei2khcos⁡θ+c17eikhcos⁡θ+c18+c19e-ikhcos⁡θ+c20e-i2khcos⁡θe-ikhsin⁡θ+c21ei2khcos⁡θ+c22eikhcos⁡θ+c23+c24e-ikhcos⁡θ+c25e-i2khcos⁡θe-i2khsin⁡θ=0,where(22)P1≔cos⁡khcos⁡θ-1,Q1≔cos⁡khsin⁡θ-1,P2≔cos⁡2khcos⁡θ-1,Q2≔cos⁡2khsin⁡θ-1. As stated in Theorem 1, to obtain a second order scheme, the parameters in (19) should satisfy (15). In addition, let all the coefficients of ghost nodes in Figure 1(c) be zero; that is,(23)T16≔a74h2+c16k2=0,T17≔a3+b1h2+c17k2=0,T18≔-2a3+b2-1/2a7h2+c18k2=0,T19≔a3+b3h2+c19k2=0,T20≔a74h2+c20k2=0,T21≔a8+b44h2+c21k2=0,T23≔-2a8+b64h2+c23k2=0,T25≔a8+b84h2+c25k2=0.Combining (15) with (23) yields(24)a1=-2a4-a5-b3-2b8-k2h2c19+4c20+8c25,a2=1+a4-a6+2b3+3b8+k2h22c19+8c20+12c25,b5=12+12k2h22c17+c18+2c20+12c21+4c23-4c25,b7=12+12k2h2c18+2c19+2c20-4c21+4c23+12c25,b1=b3+k2h2-c17+c19,b2=-2b3+k2h2-c18-2c19-2c20,b4=b8+4k2h2-c21+c25,b6=-2b8-4k2h2c23+2c25,c2=2-∑i=03c7+5i-∑i=042c3+5i+3c4+5i+4c5+5i,c6=2-6c20-∑i=710ci-3∑i=1719ci-∑i=11152ci+4ci+10,c16=c20,a3=-b3-k2h2c19,a7=-4k2h2c20,a8=-b8-4k2h2c25,c1=8c20+6c24+7c25-3+∑i=02c3+4i+2c4+4i+5c15+4i+∑i=033c5+4i+4c10+4i.Inserting (24) into (21), we have that the dispersion equation for the second kind node has the imaginary part. Therefore, the numerical phase velocity or the numerical group velocity can not be obtained anymore, which is popularly used to optimize parameters for finite difference schemes [1, 2, 24]. Hence, a proper method is needed to get optimal parameters for the second kind node.We next present a method to obtain the optimal parameters of the finite difference equation (19) for the second transitional node. We first substitute (24) into (21) to eliminate the parameters a1,  a2,  a3,  a7,  a8,  b1,b2,  b4,  b5,  b6,  b7,  c1,  c2,  c6, and c16. Under the assumption that all the parameters in the finite difference equation (19) are real numbers, we then separate the dispersion equation (21) into the real part and the imaginary part. Upon introducing the notations(25)Ri,j=cos⁡khicos⁡θ+jsin⁡θ,Si,j=sin⁡khicos⁡θ+jsin⁡θ,for eachi,j∈Z3,we, respectively, let the real part and the imaginary part of the dispersion equation be zero, which lead to the following two equations:(26)41-2R0,1P1+P2R0,2a4+P2-4P1a6+2Q2-R0,2P2+43-2R0,1P1b8+P2-4P1R0,1a5+2k2h2R2,2-2R1,2+R0,2c3+R2,2-R1,2-R2,1+R1,1c7+2R2,2-3R1,2+R1,-2c4+R2,-2+3R2,2-4R1,2c5+R0,1-R2,1+2R2,2-R1,2c8+R1,-1+3R2,2-R1,2-R2,1c9+4R2,2-R1,2-R2,1+R2,-1c10+-2R2,1+R2,2+R2,0c11+2R2,2-R2,1-R1,2+R1,0c12+2P1Q2+5R2,2-2R1,2+2R2,1+R0,2c23+5R2,2-4R1,2-2R2,1+R2,0c15+3R2,2-R2,1-R1,2+R1,-1-2Q1R1,0+Q22R1,0c17+4R2,2-2R1,2-3R2,1+R0,1-2Q1+12Q2R1,0c18+5R2,2-3R1,2+R2,1+R1,1+2Q1+Q22R1,0-4P1+1Q1c19+2Q2R1,0-R2,0+3R2,2-4R2,1+R2,-2c21+8R2,2-4R1,2-6R2,1+R2,-1+R2,1+44P1-Q1+Q2R1,0-24P1+P2R0,1c20+3R2,2-2R1,2+R2,1+1c13+4R2,2-2R2,1-3R1,2+R1,0c14+23R2,2-R1,2-2R2,1c24+4R2,2-R2,1-R1,2+R1,-2c22+243-2R0,1P1-P2R0,2+22R2,2-R1,2-R2,1+Q2R2,0+R1,0-2c25=-4P1-Q2R1,0-2k2h22R1,2+R2,1-3R2,2,P2S0,2-8P1S0,1a4+P2-4P1S0,1a5+P2S0,2-8P1S0,1b8+2k2h22S2,2-3S1,2-S1,-2c4+3S2,2-4S1,2-S2,-2c5+S2,2-2S1,2+S0,2c3+S2,2-S1,2-S2,1+S1,1c7+2S2,2-S1,2-S2,1+S0,1c8+3S2,2-S1,2-S2,1-S1,-1c9+4S2,2-S1,2-S2,1-S2,-1c10+S2,2-2S2,1+S2,0c11+2S2,2-S2,1-S1,2+S1,0c12+3S2,2-2S1,2+S2,1c13+4S2,2-3S1,2-2S2,1-S1,0c14+5S2,2-4S1,2-2S2,1-S2,0c15+Q22-2Q1S1,0+3S2,2-S2,1-S1,2+S1,-1c17+4S2,2-2S1,2-3S2,1-S0,1c18+2Q1-Q22S1,0+5S2,2-S1,1-3S1,2+S2,1c19+2P2-4P1S0,1+8S2,2-4S1,2-6S2,1+S2,-1-S2,1c20+3S2,2-4S2,1+S2,-2+22S1,0-S2,0Q2c21+4S2,2-S2,1-S1,2+S1,-2c22+5S2,2-2S1,2-4S2,1-S0,2c23+6S2,2-3S1,2-4S2,1-S1,2c24+-16P1S0,1+2P2S0,2+23S2,2-2S1,2+S2,1+2Q2S2,0-2S1,0c25=-2k2h2-3S2,2+2S1,2+S2,1.Note that (26) are functions of kh and θ. To reduce the numerical dispersion and improve the accuracy of the difference scheme for the second kind node, we propose the following rule.Rule 1 (refined choice strategy).
Consider the following.
Step 1. Estimate kh, where k is the wavenumber for the second kind node and h is the step size for the fine grid.
Step 2. Let Iθ≔[0,2π) and choose θ=θm≔(m-1)2π/l∈Iθ,  m=1,2,…,l, where l is a given positive integer.
Step 3. Inserting above kh and θ to (26) yields a linear system.
Step 4. Applying the singular value decomposition method to solve the corresponding system yields the optimal weight parameters for (26). Other parameters are obtained by (24).LetG be the number of gridpoints per wavelength and there holds kh=2π/G [24]. According to G, some optimal parameters for the second kind of transitional nodes are presented in Figure 2(a). This figure only shows five parameters in (26). We find that the parameters a4,  a6,  b3, and  c3 vary smoothly, while the parameter c11 is very sensitive to the number of gridpoints per wavelength G, which may be the effect of the ghost nodes.Figure 2
: Optimal coefficients for (a) the second kind node, (b) the third kind node.
(a)
(b)


### 2.3. Optimal 25-Point Finite Difference Schemes for Other Nodes

In this subsection, we first present optimal 25-point finite difference methods for the last nine kinds of transitional nodes. According to [24], we then formulate optimal 25-point finite difference schemes for other nodes, which are not transitional nodes. All the schemes are second order in accuracy.For the first kind node (see Figure1(b)), we construct its consistent finite difference scheme by the method proposed in [21]. The formulation of a consistent finite difference scheme for the third kind node (see Figure 1(d)) is presented in the appendix. For the fourth kind node, we construct its consistent finite difference scheme with its nearest 24 coarse points. Polynomial interpolation method with fourth order in accuracy is applied to formulate a second order consistent finite difference scheme for the fifth kind of transitional node (see Figure 1(e)). The eighth kind of transitional nodes are dealt with in a similar way. For the sixth kind of transitional nodes, we construct its consistent finite difference method with its nearest 24 symmetric nodes. The finite difference schemes for the seventh and ninth kind of transitional nodes are formulated in a similar way. We refer the interested readers to equations (3.3) and (3.25) in [24] for the details. Finally, for the tenth kind of transitional nodes, we construct its finite difference scheme with its nearest 8 fine points (see equation (3.2) in [22]). All these schemes are consistent with the Helmholtz equation with PML and are second order in accuracy (see [21, 22, 24]). In addition, the parameters of the finite difference methods can be chosen with refined choice strategy.For other nodes which are not transitional nodes, we formulate the optimal 25-point finite difference scheme based on our previous work [24]. Together with the work for the transitional nodes, we have an optimal 25-point finite-difference subgridding scheme for the 2D Helmholtz equation with PML, referred to as the refined subgridding 25p for short, which is second order in accuracy and consistent with the equation.


## 2.1. A Consistent 25-Point Finite Difference Scheme for the Second Kind Node

We next present a consistent 25-point finite difference method for the second kind of transitional nodes, which is second order in accuracy. Seen from Figure1(a), we find that the second kind of transitional nodes lie in the interior domain, in which the Helmholtz equation with PML (2) deteriorates to be the Helmholtz equation. With eight ghost nodes, we will construct a 25-point stencil for the second kind node. This stencil is shown in Figure 1(c), where the nodes in the bottom red line are ghost nodes.We turn to constructing a 25-point finite difference scheme for the second kind node by the weighted averaged method (see [22–24] and the references therein). For simplicity, we first introduce some notations. 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. For each m and n, we set um,n≔ux=xm,y=yn and km,n≔kx=xm,y=yn. Next, different from [24], more weight parameters will be introduced to relax the restriction for the symmetry of the finite difference stencil. In detail, let a1,a2,…,a8,  b1,b2,…,b8,  c1,c2,…,c25 be constants to be determined in the following. With these weight parameters, we then establish an approximation for the Laplacian term by the weighted average method. Let(5)Lh,xum,n+j≔um+1,n+j+um-1,n+j-2um,n+jh2,for j∈Z2≔{-1,0,1}, and(6)L2h,xum,n+j≔um+2,n+j+um-2,n+j-2um,n+j2h2,for j∈Z3≔{-2,-1,0,1,2}. Upon introducing two notations(7)L~h,xux=xm,y=yn≔a1Lh,xum,n-1+a2Lh,xum,n+a3Lh,xum,n+1,L~2h,xux=xm,y=yn≔a4L2h,xum,n-2+a5L2h,xum,n-1+a6L2h,xum,n+a7L2h,xum,n+1+a8L2h,xum,n+2,we approximate ∂2u/∂x2 as (8)∂2u∂x2x=xm,y=yn≈L~h,xux=xm,y=yn+L~2h,xux=xm,y=yn. Similar treatment with respect to ∂2u/∂y2 leads to(9)∂2u∂y2x=xm,y=yn≈L~h,yux=xm,y=yn+L~2h,yux=xm,y=yn, where(10)L~h,yux=xm,y=yn≔b1Lh,yum-1,n+b2Lh,yum,n+b3Lh,yum+1,n,L~2h,yux=xm,y=yn≔b4L2h,yum-2,n+b5L2h,yum-1,n+b6L2h,yum,n+b7L2h,yum+1,n+b8L2h,yum+2,n.By setting L~h≔L~h,x+L~h,y+L~2h,x+L~2h,y, we construct an approximation for the Laplacian term as(11)Δux=xm,y=yn≈L~hux=xm,y=yn. We then approximate k2ux=xm,y=yn by a weighted average: (12)k2ux=xm,y=yn≈I~hk2ux=xm,y=yn, where(13)I~hk2ux=xm,y=yn≔∑l=04∑j=15cj+5lk2um-3+j,n-2+l.We now establish a 25-point difference approximation for the Helmholtz equation as(14)L~hux=xm,y=yn+I~hk2ux=xm,y=yn=0.We next consider the consistency of the 25-point difference scheme (14) and give a convergence analysis. The results are presented in the following theorem.Theorem 1.
(i) If∑i=18ai=1, ∑i=18bi=1, and ∑i=125ci=1, then the 25-point finite difference approximation (14) is pointwise consistent with the Helmholtz equation (1) and is a first order scheme.
(ii) If the parameters in (14) satisfy(15)∑i=18ai=1,∑i=18bi=1,∑i=125ci=1,a3+a7+2a8-a1-2a4-a5=0,b3+b7+2b8-b1-2b4-b5=0,2∑i=15c20+i-ci+∑i=610c10+i-ci=0,∑i=04c4+5i+2c5+5i-2c1+5i-c2+5i=0,then the 25-point finite difference approximation (14) is pointwise consistent with the Helmholtz equation (1) and is a second order scheme.Proof.
Assume thatxm≤x<xm+1 and yn≤y<yn+1. It follows from the Taylor theorem that(16)L~hux=xm,y=yn=∑i=18ai∂2u∂x2+∑i=18bi∂2u∂y2+ν1h+Oh2,I~hk2ux=xm,y=yn=∑i=125cik2u+ν2h+Oh2,where(17)ν1≔a3-a1+2a8-a4+a7-a5∂3u∂x2∂y+b3-b1+2b8-b4+b7-b5∂3u∂x∂y2,ν2≔∑i=04c4+5i+2c5+5i-2c1+5i-c2+5i∂ku∂x+2∑i=15c20+i-ci+∑i=610c10+i-ci∂ku∂y. Therefore, when ∑i=18ai=1, ∑i=18bi=1, and ∑i=125ci=1 hold, we have that the left hand side of the 25-point finite difference approximation (14) is equivalent to(18)∂2u∂x2+∂2u∂y2+k2u+ηh+Oh2,where η≔ν1+ν2. Together with the concept of consistency (see [13]), we conclude the results of this theorem.We now turn to constructing a second order finite difference scheme for the second kind node in Figure1(c), which is also consistent with the Helmholtz equation. We begin with rewriting the finite difference scheme (14). Substituting (7), (10), and (13) into (14) and replacing um+i,n+j with Um+i,n+j  (i,j∈Z2) yield the 25-point finite difference equation(19)T1Um-2,n-2+T2Um-1,n-2+T3Um,n-2+T4Um+1,n-2+T5Um+2,n-2+T6Um-2,n-1+T7Um-1,n-1+T8Um,n-1+T9Um+1,n-2+T10Um+2,n-1+T11Um-2,n+T12Um-1,n+T13Um,n+T14Um+1,n+T15Um+2,n+T16Um-2,n+1+T17Um-1,n+1+T18Um,n+1+T19Um+1,n+1+T20Um+2,n+1+T21Um-2,n+2+T22Um-1,n+2+T23Um,n+2+T24Um+1,n+2+T25Um+2,n+2=0,where the coefficients are given by(20)T1≔a4+b44h2+c1km-2,n-22,T2≔b54h2+c2km-1,n-22,T3≔-2a4+b64h2+c3km,n-22,T4≔b74h2+c4km+1,n-22,T5≔a4+b84h2+c5km+2,n-22,T6≔a54h2+c6km-2,n-12,T7≔a1+b1h2+c7km-1,n-12,T8≔-2a1+b2-1/2a5h2+c8km,n-12,T9≔a1+b3h2+c9km+1,n-12,T10≔a54h2+c10km+2,n-12,T11≔a6-2b44h2+c11km-2,n2,T12≔a2-2b1-1/2b5h2+c12km-1,n2,T13≔-4a2+b2-a6+b62h2+c13km,n2,T14≔a2-2b3-1/2b7h2+c14km+1,n2,T15≔a6-2b84h2+c15km+2,n2,T16≔a74h2+c16km-2,n+12,T17≔a3+b1h2+c17km-1,n+12,T18≔-2a3+b2-1/2a7h2+c18km,n+12,T19≔a3+b3h2+c19km+1,n+12,T20≔a74h2+c20km+2,n+12,T21≔a8+b44h2+c21km-2,n+22,T22≔b54h2+c22km-1,n+22,T23≔-2a8+b64h2+c23km,n+22,T24≔b74h2+c24km+1,n+22,T25≔a8+b84h2+c25km+2,n+22.To guarantee the above scheme to be second order in accuracy, the parameters in (19) should satisfy the conditions proposed in Theorem 1. On the other hand, as the ghost nodes do not exist in the stencil for the second kind node, we let their coefficients be zero; that is, T16=0,  T17=0,  T18=0,  T19=0,  T20=0,  T21=0,  T23=0, and T25=0. This indicates that the parameters in the scheme (19) may vary according to the wavenumber k. By far, a second order finite difference scheme for the second kind node is formulated. Choice strategies for optimal parameters of the finite difference scheme will be presented in the next subsection.


## 2.2. Choice Strategies for Optimal Parameters for the Second Kind Node’s Stencil

We now go further into the problem of the choice strategies for optimal parameters for the second kind node’s stencil (see Figure1(c)). To optimize parameters for the second kind node’s stencil, we first perform classical dispersion analysis by assuming a plane-wave solution of the form U(x,y)=e-ik(xcos⁡θ+ysin⁡θ), 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.SubstitutingUm,n≔e-ik(xmcos⁡θ+ynsin⁡θ) into (19), dividing both sides by the factor e-ik(xmcos⁡θ+ynsin⁡θ), and finally applying the Euler formula eix=cos⁡x+isin⁡x lead to the following dispersion equation:(21)a4ei2khsin⁡θ+a5eikhsin⁡θ+a6+a7e-ikhsin⁡θ+a8e-i2khsin⁡θP2+b4ei2khcos⁡θ+b5eikhcos⁡θ+b6+b7e-ikhcos⁡θ+b8e-i2khcos⁡θQ2+4a1eikhsin⁡θ+a2+a3e-ikhsin⁡θP1+4b1eikhcos⁡θ+b2+b3e-ikhcos⁡θQ1+2k2h2c1ei2khcos⁡θ+c2eikhcos⁡θ+c3+c4e-ikhcos⁡θ+c5e-i2khcos⁡θei2khsin⁡θ+c6ei2khcos⁡θ+c7eikhcos⁡θ+c8+c9e-ikhcos⁡θ+c10e-i2khcos⁡θeikhsin⁡θ+c11ei2khcos⁡θ+c12eikhcos⁡θ+c13+c14e-ikhcos⁡θ+c15e-i2khcos⁡θ+c16ei2khcos⁡θ+c17eikhcos⁡θ+c18+c19e-ikhcos⁡θ+c20e-i2khcos⁡θe-ikhsin⁡θ+c21ei2khcos⁡θ+c22eikhcos⁡θ+c23+c24e-ikhcos⁡θ+c25e-i2khcos⁡θe-i2khsin⁡θ=0,where(22)P1≔cos⁡khcos⁡θ-1,Q1≔cos⁡khsin⁡θ-1,P2≔cos⁡2khcos⁡θ-1,Q2≔cos⁡2khsin⁡θ-1. As stated in Theorem 1, to obtain a second order scheme, the parameters in (19) should satisfy (15). In addition, let all the coefficients of ghost nodes in Figure 1(c) be zero; that is,(23)T16≔a74h2+c16k2=0,T17≔a3+b1h2+c17k2=0,T18≔-2a3+b2-1/2a7h2+c18k2=0,T19≔a3+b3h2+c19k2=0,T20≔a74h2+c20k2=0,T21≔a8+b44h2+c21k2=0,T23≔-2a8+b64h2+c23k2=0,T25≔a8+b84h2+c25k2=0.Combining (15) with (23) yields(24)a1=-2a4-a5-b3-2b8-k2h2c19+4c20+8c25,a2=1+a4-a6+2b3+3b8+k2h22c19+8c20+12c25,b5=12+12k2h22c17+c18+2c20+12c21+4c23-4c25,b7=12+12k2h2c18+2c19+2c20-4c21+4c23+12c25,b1=b3+k2h2-c17+c19,b2=-2b3+k2h2-c18-2c19-2c20,b4=b8+4k2h2-c21+c25,b6=-2b8-4k2h2c23+2c25,c2=2-∑i=03c7+5i-∑i=042c3+5i+3c4+5i+4c5+5i,c6=2-6c20-∑i=710ci-3∑i=1719ci-∑i=11152ci+4ci+10,c16=c20,a3=-b3-k2h2c19,a7=-4k2h2c20,a8=-b8-4k2h2c25,c1=8c20+6c24+7c25-3+∑i=02c3+4i+2c4+4i+5c15+4i+∑i=033c5+4i+4c10+4i.Inserting (24) into (21), we have that the dispersion equation for the second kind node has the imaginary part. Therefore, the numerical phase velocity or the numerical group velocity can not be obtained anymore, which is popularly used to optimize parameters for finite difference schemes [1, 2, 24]. Hence, a proper method is needed to get optimal parameters for the second kind node.We next present a method to obtain the optimal parameters of the finite difference equation (19) for the second transitional node. We first substitute (24) into (21) to eliminate the parameters a1,  a2,  a3,  a7,  a8,  b1,b2,  b4,  b5,  b6,  b7,  c1,  c2,  c6, and c16. Under the assumption that all the parameters in the finite difference equation (19) are real numbers, we then separate the dispersion equation (21) into the real part and the imaginary part. Upon introducing the notations(25)Ri,j=cos⁡khicos⁡θ+jsin⁡θ,Si,j=sin⁡khicos⁡θ+jsin⁡θ,for eachi,j∈Z3,we, respectively, let the real part and the imaginary part of the dispersion equation be zero, which lead to the following two equations:(26)41-2R0,1P1+P2R0,2a4+P2-4P1a6+2Q2-R0,2P2+43-2R0,1P1b8+P2-4P1R0,1a5+2k2h2R2,2-2R1,2+R0,2c3+R2,2-R1,2-R2,1+R1,1c7+2R2,2-3R1,2+R1,-2c4+R2,-2+3R2,2-4R1,2c5+R0,1-R2,1+2R2,2-R1,2c8+R1,-1+3R2,2-R1,2-R2,1c9+4R2,2-R1,2-R2,1+R2,-1c10+-2R2,1+R2,2+R2,0c11+2R2,2-R2,1-R1,2+R1,0c12+2P1Q2+5R2,2-2R1,2+2R2,1+R0,2c23+5R2,2-4R1,2-2R2,1+R2,0c15+3R2,2-R2,1-R1,2+R1,-1-2Q1R1,0+Q22R1,0c17+4R2,2-2R1,2-3R2,1+R0,1-2Q1+12Q2R1,0c18+5R2,2-3R1,2+R2,1+R1,1+2Q1+Q22R1,0-4P1+1Q1c19+2Q2R1,0-R2,0+3R2,2-4R2,1+R2,-2c21+8R2,2-4R1,2-6R2,1+R2,-1+R2,1+44P1-Q1+Q2R1,0-24P1+P2R0,1c20+3R2,2-2R1,2+R2,1+1c13+4R2,2-2R2,1-3R1,2+R1,0c14+23R2,2-R1,2-2R2,1c24+4R2,2-R2,1-R1,2+R1,-2c22+243-2R0,1P1-P2R0,2+22R2,2-R1,2-R2,1+Q2R2,0+R1,0-2c25=-4P1-Q2R1,0-2k2h22R1,2+R2,1-3R2,2,P2S0,2-8P1S0,1a4+P2-4P1S0,1a5+P2S0,2-8P1S0,1b8+2k2h22S2,2-3S1,2-S1,-2c4+3S2,2-4S1,2-S2,-2c5+S2,2-2S1,2+S0,2c3+S2,2-S1,2-S2,1+S1,1c7+2S2,2-S1,2-S2,1+S0,1c8+3S2,2-S1,2-S2,1-S1,-1c9+4S2,2-S1,2-S2,1-S2,-1c10+S2,2-2S2,1+S2,0c11+2S2,2-S2,1-S1,2+S1,0c12+3S2,2-2S1,2+S2,1c13+4S2,2-3S1,2-2S2,1-S1,0c14+5S2,2-4S1,2-2S2,1-S2,0c15+Q22-2Q1S1,0+3S2,2-S2,1-S1,2+S1,-1c17+4S2,2-2S1,2-3S2,1-S0,1c18+2Q1-Q22S1,0+5S2,2-S1,1-3S1,2+S2,1c19+2P2-4P1S0,1+8S2,2-4S1,2-6S2,1+S2,-1-S2,1c20+3S2,2-4S2,1+S2,-2+22S1,0-S2,0Q2c21+4S2,2-S2,1-S1,2+S1,-2c22+5S2,2-2S1,2-4S2,1-S0,2c23+6S2,2-3S1,2-4S2,1-S1,2c24+-16P1S0,1+2P2S0,2+23S2,2-2S1,2+S2,1+2Q2S2,0-2S1,0c25=-2k2h2-3S2,2+2S1,2+S2,1.Note that (26) are functions of kh and θ. To reduce the numerical dispersion and improve the accuracy of the difference scheme for the second kind node, we propose the following rule.Rule 1 (refined choice strategy).
Consider the following.
Step 1. Estimate kh, where k is the wavenumber for the second kind node and h is the step size for the fine grid.
Step 2. Let Iθ≔[0,2π) and choose θ=θm≔(m-1)2π/l∈Iθ,  m=1,2,…,l, where l is a given positive integer.
Step 3. Inserting above kh and θ to (26) yields a linear system.
Step 4. Applying the singular value decomposition method to solve the corresponding system yields the optimal weight parameters for (26). Other parameters are obtained by (24).LetG be the number of gridpoints per wavelength and there holds kh=2π/G [24]. According to G, some optimal parameters for the second kind of transitional nodes are presented in Figure 2(a). This figure only shows five parameters in (26). We find that the parameters a4,  a6,  b3, and  c3 vary smoothly, while the parameter c11 is very sensitive to the number of gridpoints per wavelength G, which may be the effect of the ghost nodes.Figure 2
: Optimal coefficients for (a) the second kind node, (b) the third kind node.
(a)
(b)


## 2.3. Optimal 25-Point Finite Difference Schemes for Other Nodes

In this subsection, we first present optimal 25-point finite difference methods for the last nine kinds of transitional nodes. According to [24], we then formulate optimal 25-point finite difference schemes for other nodes, which are not transitional nodes. All the schemes are second order in accuracy.For the first kind node (see Figure1(b)), we construct its consistent finite difference scheme by the method proposed in [21]. The formulation of a consistent finite difference scheme for the third kind node (see Figure 1(d)) is presented in the appendix. For the fourth kind node, we construct its consistent finite difference scheme with its nearest 24 coarse points. Polynomial interpolation method with fourth order in accuracy is applied to formulate a second order consistent finite difference scheme for the fifth kind of transitional node (see Figure 1(e)). The eighth kind of transitional nodes are dealt with in a similar way. For the sixth kind of transitional nodes, we construct its consistent finite difference method with its nearest 24 symmetric nodes. The finite difference schemes for the seventh and ninth kind of transitional nodes are formulated in a similar way. We refer the interested readers to equations (3.3) and (3.25) in [24] for the details. Finally, for the tenth kind of transitional nodes, we construct its finite difference scheme with its nearest 8 fine points (see equation (3.2) in [22]). All these schemes are consistent with the Helmholtz equation with PML and are second order in accuracy (see [21, 22, 24]). In addition, the parameters of the finite difference methods can be chosen with refined choice strategy.For other nodes which are not transitional nodes, we formulate the optimal 25-point finite difference scheme based on our previous work [24]. Together with the work for the transitional nodes, we have an optimal 25-point finite-difference subgridding scheme for the 2D Helmholtz equation with PML, referred to as the refined subgridding 25p for short, which is second order in accuracy and consistent with the equation.


## 3. Numerical Experiments

In this section, we present two 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 3.40 GHz and 4.00 Gb RAM.


### 3.1. A Near-Surface Low-Velocity Layer Model

To validate the efficiency of the refined subgridding 25p over the optimal 25-point finite difference scheme (the refined uniform 25p) proposed in [24], we consider a near-surface low-velocity layer model shown in Figure 3(a). The domain is [0m,1600m]×[0m,1600m] and there are four velocities: from the top, v≔1200 m/s, 2800 m/s, 3000 m/s, and 3500 m/s. A point source δ(x-xs,y-ys)R^(f,f0) is located at the point (xs,ys)=(800m,20m), where R^(f,f0) is the Ricker wavelet in the frequency domain defined by (27)R^f,f0≔∫-∞+∞1-2π2f02t2exp⁡-π2f02t2ei2πftdt, whose dominant frequency is f0=25Hz. 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 subgridding 25p, H and h denote the coarse and fine step size, respectively. In this experiment, we let H=10m and h=5m. Specifically, we only use fine grids for the top layer whose velocity is v=1200m/s.Figure 3
(a) A near-surface low-velocity layer model, (b) the real part of the solution forf=20Hz obtained by the refined subgridding 25p with (H,h)=(10m,5m), (c) the common-shot-point records for the wavefield by the refined uniform 25p with uniform grid size hu=5m, (d) the common-shot-point records for the wavefield by the refined subgridding 25p with (H,h)=(10m,5m), (e) synthetic seismograms generated by the refined uniform 25p with hu=5m for the point (300 m, 0 m), and (f) synthetic seismograms generated by the refined subgridding 25p with (H,h)=(10m,5m) for the point (300 m, 0 m).
(a)
(b)
(c)
(d)
(e)
(f)In Figure3(b), the monofrequency wavefield (real part) for f=20Hz obtained by the refined subgridding 25p with (H,h)=(10m,5m) is presented. 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 phenomenon obeys Snell’s law. In Figures 3(c) and 3(d), the common-shot-point records by the refined uniform 25p with hu=5m and the refined subgridding 25p with (H,h)=(10m,5m) are shown, respectively. In the two 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 Figure 3(d)). It is found that Figure 3(d) is almost as clear as Figure 3(c). However, when comparing the cost of Figure 3(d) with that of Figure 3(c), there is a 53.5208% decrease in memory and a 71.0987% decrease in the computational time. Specifically, the computational time for Figure 3(c) is about 55.9861 hours, while that for Figure 3(d) is only 16.1807 hours. Therefore, the efficiency of the refined subgridding 25p is confirmed. Finally, in Figures 3(e) and 3(f) we present synthetic seismograms generated by the refined uniform 25p with hu=5m and that generated by the refined subgridding 25p with (H,h)=(10m,5m) for the point (300 m, 0 m), respectively. These two figures look like the same, although there are some nonphysical oscillations in Figure 3(f), which may be caused by the reflections from the dense to coarse mesh interface. However, there is a sharp decrease in the computational time for obtaining Figure 3(f).


### 3.2. A Two-Layered Model

To compare the refined subgridding 25p with the refined uniform 25p in a more quantitative way, we consider a two-layered model [25].For the observation pointx=(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(28)ΔGx,xs+k2xGx,xs=-δx-xs,with continuity conditions(29)Gx,xsy=0+=Gx,xsy=0-,∂Gx,xs∂yy=0+=∂Gx,xs∂yy=0-,with the wavenumber(30)kx=k1,fory>0,k2,fory<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. We remark that the condition that k1<k2 indicates the velocity in the layer with y<0 is slower than that in the layer with y>0, and the solution of the above problem characterizes the propagation of the wave for a fixed frequency in a two-layered model, which is generated by a point source.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 [25]:(i)
y<0,(31)Gx,xs=Φx,xs+i4π∫-∞+∞g1x,xs,ξdξ,(ii)
y>0,(32)Gx,xs=i2π∫-∞+∞g2x,xs,ξdξ, where(33)Φx,xs=i4H01k2x-xs,in which H0(1) is the Hankel function of the first kind with order zero, i2=-1, and(34)
g
1
x
,
x
s
,
ξ
=
ξ
2
-
k
2
2
-
ξ
2
-
k
1
2
i
ξ
2
-
k
2
2
ξ
2
-
k
1
2
+
ξ
2
-
k
2
2
e
ξ
2
-
k
2
2
y
+
y
s
e
i
ξ
x
-
x
s
,
for
ξ
>
k
2
,
k
2
2
-
ξ
2
-
i
ξ
2
-
k
1
2
k
2
2
-
ξ
2
i
ξ
2
-
k
1
2
+
k
2
2
-
ξ
2
e
-
i
k
2
2
-
ξ
2
y
+
y
s
e
i
ξ
x
-
x
s
,
for
k
1
<
ξ
<
k
2
,
k
2
2
-
ξ
2
-
k
1
2
-
ξ
2
k
2
2
-
ξ
2
k
1
2
-
ξ
2
+
k
2
2
-
ξ
2
e
-
i
k
2
2
-
ξ
2
y
+
y
s
e
i
ξ
x
-
x
s
,
for
ξ
<
k
1
,
g
2
x
,
x
s
,
ξ
=
e
-
ξ
2
-
k
1
2
y
+
ξ
2
-
k
2
2
y
s
i
ξ
2
-
k
1
2
+
ξ
2
-
k
2
2
e
i
ξ
x
-
x
s
,
for
ξ
>
k
2
,
e
-
ξ
2
-
k
1
2
y
-
i
k
2
2
-
ξ
2
y
s
i
ξ
2
-
k
1
2
+
k
2
2
-
ξ
2
e
i
ξ
x
-
x
s
,
for
k
1
<
ξ
<
k
2
,
e
i
k
1
2
-
ξ
2
y
-
k
2
2
-
ξ
2
y
s
k
1
2
-
ξ
2
+
k
2
2
-
ξ
2
e
i
ξ
x
-
x
s
,
for
ξ
<
k
1
.In practice, the first term and the second term on the left hand side of (31) are the downward incident waves and the reflected waves, respectively. Equation (31) indicates that the incident waves will be interfered with the reflected waves in the layer with y<0. Furthermore, (32) represents the transmitted waves, which presents that there are only transmitted waves in the layer with y>0. Visual explanations will be given in Figure 4.Figure 4
(a) The real part of the analytical solution, (b) the real part of the solution generated by the refined uniform 25p withhu=3m, (c) the real part of the solution generated by the refined subgridding 25p with (H,h)=(6m,3m), and (d) the real part of solutions in the line of x=120m generated by different methods.
(a)
(b)
(c)
(d)We then consider a two-layered problem (28) with the source point fixed at (300m,-88m), k1=0.0838 and k2=0.2094. The bounded domain [0m,600m]×[-112m,488m] is our area of interest. The corresponding solution is given via (31)~(33). Thus, numerical integrals have to be dealt with to obtain the solution. In this experiment, we apply the Gauss-type quadratures [26] for weakly singular integrals appearing in (31) and (32). This kind of quadrature rules possesses polynomial order of accuracy and enables 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 25p and the refined subgridding 25p including errors, computational times, number of unknowns, and the errors’ decay exponents (briefly, decay exp.). The error is measured inC-norm, which is defined as follows: for any complex vector z=z1,z2,…,zM, (35)zC≔max1≤j≤M⁡zj, where |zj| is the complex modulus of zj. The C-norm measures the error related to the amplitude of the waves. In addition, 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 PML. In addition, let n~ represent the number of all unknowns which lie in the interior domain in the layer with k(x)=0.0838. Each decay exponent is computed by ln⁡(en~/em~)/ln⁡(m~/n~), where en~ denotes the numerical solution’s error in the layer with k(x)=0.0838 for one test. For the refined subgridding 25p, 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.2094. Therefore, h=H/2 holds.In Tables1 and 2, we present the number of all unknowns in the two-dimensional computational domain, the errors and computational time for the refined uniform 25p and subgridding 25p, respectively. We find that the error of the refined subgridding 25p with (H,h) is much smaller than that by the refined uniform 25p with hu=H, while it is a little bigger than that by the refined uniform 25p with hu=h. However, there is a sharp decrease in both the number of unknowns and the computational time for the refined subgridding 25p with (H,h), when compared with the refined uniform 25p with hu=h. It is also found that, for the two methods, the error for the layer with k(x)=0.2094 decreases slowly when we reduce the step sizes. This is because in each test, the point source function produces a solution which is nonsmooth close to the source and the interface between two layers (see Figure 4(a)). In fact, the solution can be subdivided into a near- and far-field contribution. As in the top layer the incident waves are interfered with the reflected waves returned from the bottom of the top layer (see Figure 4(a)), we take the solution for the layer with k(x)=0.0838 as the far-field contribution.Table 1
The performance of the refined uniform 25p. The pair(⋅,⋅) denotes the error for the layer with k(x)=0.2094 and the layer with k(x)=0.0838, respectively.
h
u
12 m
6 m
3 m
1.5 m
n
8281
25921
68121
231361
Errors
(0.2915, 0.0888)
(0.0282, 0.0119)
(0.0205, 0.0075)
(0.0185, 0.0025)
Time
39.7491 s
396.5701 s
2.5210
e
+
003 s
2.5780
e
+
004 sTable 2
The performance of the refined subgridding 25p. The pair(⋅,⋅) denotes the error for the layer with k(x)=0.2094 and the layer with k(x)=0.0838, respectively.
(
H
,
h
)
(12 m, 6 m)
(6 m, 3 m)
(3 m, 1.5 m)
n
11432
36602
100100
Errors
(0.0370, 0.0261)
(0.0208, 0.0077)
(0.0185, 0.0029)
Time
73.3205 s
796.7439 s
5.6267
e
+
003 sWe next consider the decay exponent of the refined subgridding 25p in the far-field contribution. Table3 presents the decay exponent for the refined subgridding 25p. We find that the decay exponent is one. As n~ denotes the number of the unknowns which lie in the layer with k(x)=0.0838, the decay exponent being one indicates that the refined subgridding 25p is second order in accuracy.Table 3
The decay exponent for the refined subgridding 25p on the layer withk(x)=0.0838.
n
~
2393
8984
33965
Errors
0.0261
0.0077
0.0029
Decay exp.
0.9194
0.7395For further comparison, in Figures4(a)–4(c) we present the real part of the solution for the two-layered problem in the domain [0m,600m]×[-112m,488m] by the analytical method, the refined uniform 25p with hu=3m, and the refined subgridding 25p with (H,h)=(6m,3m), 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 phenomenon obeys Snell’s law. It is easy to see that both Figures 4(b) and 4(c) are similar to Figure 4(a), but Figure 4(c) needs less memory and computational time than those needed by Figure 4(b) (see Tables 1 and 2), which verifies the efficiency of the refined subgridding 25p. In addition, almost no boundary reflections exist in Figures 4(b) and 4(c), which confirms the efficiency of the PML’s absorption. Finally, Figure 4(d) presents the real part of the two-layered problem’s solution in the line of x=120m, respectively, by the analytical method, the refined uniform 25p with hu=12m, the refined subgridding 25p with (H,h)=(12m,6m), and the refined subgridding 25p with (H,h)=(3m,1.5m). From this figure, it is obvious that the solution by the refined uniform 25p with hu=12m has a wavelength which is different from the analytical solution. This is called the “dispersion,’’ which leads to “the pollution error’’ that has a global character (see [3]). We find that the solution by the refined subgridding 25p with (H,h)=(12m,6m) seems better. Furthermore, comparing the two solutions obtained by the refined subgridding 25p with different step sizes, we see that the solution obtained by the refined subgridding 25p converges to the analytical solution.Now, we come to our conclusion that when balancing the accuracy against the computational cost, the refined subgridding 25p possesses much improvement over the refined uniform 25p, especially for large-scale layered problems.


## 3.1. A Near-Surface Low-Velocity Layer Model

To validate the efficiency of the refined subgridding 25p over the optimal 25-point finite difference scheme (the refined uniform 25p) proposed in [24], we consider a near-surface low-velocity layer model shown in Figure 3(a). The domain is [0m,1600m]×[0m,1600m] and there are four velocities: from the top, v≔1200 m/s, 2800 m/s, 3000 m/s, and 3500 m/s. A point source δ(x-xs,y-ys)R^(f,f0) is located at the point (xs,ys)=(800m,20m), where R^(f,f0) is the Ricker wavelet in the frequency domain defined by (27)R^f,f0≔∫-∞+∞1-2π2f02t2exp⁡-π2f02t2ei2πftdt, whose dominant frequency is f0=25Hz. 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 subgridding 25p, H and h denote the coarse and fine step size, respectively. In this experiment, we let H=10m and h=5m. Specifically, we only use fine grids for the top layer whose velocity is v=1200m/s.Figure 3
(a) A near-surface low-velocity layer model, (b) the real part of the solution forf=20Hz obtained by the refined subgridding 25p with (H,h)=(10m,5m), (c) the common-shot-point records for the wavefield by the refined uniform 25p with uniform grid size hu=5m, (d) the common-shot-point records for the wavefield by the refined subgridding 25p with (H,h)=(10m,5m), (e) synthetic seismograms generated by the refined uniform 25p with hu=5m for the point (300 m, 0 m), and (f) synthetic seismograms generated by the refined subgridding 25p with (H,h)=(10m,5m) for the point (300 m, 0 m).
(a)
(b)
(c)
(d)
(e)
(f)In Figure3(b), the monofrequency wavefield (real part) for f=20Hz obtained by the refined subgridding 25p with (H,h)=(10m,5m) is presented. 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 phenomenon obeys Snell’s law. In Figures 3(c) and 3(d), the common-shot-point records by the refined uniform 25p with hu=5m and the refined subgridding 25p with (H,h)=(10m,5m) are shown, respectively. In the two 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 Figure 3(d)). It is found that Figure 3(d) is almost as clear as Figure 3(c). However, when comparing the cost of Figure 3(d) with that of Figure 3(c), there is a 53.5208% decrease in memory and a 71.0987% decrease in the computational time. Specifically, the computational time for Figure 3(c) is about 55.9861 hours, while that for Figure 3(d) is only 16.1807 hours. Therefore, the efficiency of the refined subgridding 25p is confirmed. Finally, in Figures 3(e) and 3(f) we present synthetic seismograms generated by the refined uniform 25p with hu=5m and that generated by the refined subgridding 25p with (H,h)=(10m,5m) for the point (300 m, 0 m), respectively. These two figures look like the same, although there are some nonphysical oscillations in Figure 3(f), which may be caused by the reflections from the dense to coarse mesh interface. However, there is a sharp decrease in the computational time for obtaining Figure 3(f).


## 3.2. A Two-Layered Model

To compare the refined subgridding 25p with the refined uniform 25p in a more quantitative way, we consider a two-layered model [25].For the observation pointx=(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(28)ΔGx,xs+k2xGx,xs=-δx-xs,with continuity conditions(29)Gx,xsy=0+=Gx,xsy=0-,∂Gx,xs∂yy=0+=∂Gx,xs∂yy=0-,with the wavenumber(30)kx=k1,fory>0,k2,fory<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. We remark that the condition that k1<k2 indicates the velocity in the layer with y<0 is slower than that in the layer with y>0, and the solution of the above problem characterizes the propagation of the wave for a fixed frequency in a two-layered model, which is generated by a point source.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 [25]:(i)
y<0,(31)Gx,xs=Φx,xs+i4π∫-∞+∞g1x,xs,ξdξ,(ii)
y>0,(32)Gx,xs=i2π∫-∞+∞g2x,xs,ξdξ, where(33)Φx,xs=i4H01k2x-xs,in which H0(1) is the Hankel function of the first kind with order zero, i2=-1, and(34)
g
1
x
,
x
s
,
ξ
=
ξ
2
-
k
2
2
-
ξ
2
-
k
1
2
i
ξ
2
-
k
2
2
ξ
2
-
k
1
2
+
ξ
2
-
k
2
2
e
ξ
2
-
k
2
2
y
+
y
s
e
i
ξ
x
-
x
s
,
for
ξ
>
k
2
,
k
2
2
-
ξ
2
-
i
ξ
2
-
k
1
2
k
2
2
-
ξ
2
i
ξ
2
-
k
1
2
+
k
2
2
-
ξ
2
e
-
i
k
2
2
-
ξ
2
y
+
y
s
e
i
ξ
x
-
x
s
,
for
k
1
<
ξ
<
k
2
,
k
2
2
-
ξ
2
-
k
1
2
-
ξ
2
k
2
2
-
ξ
2
k
1
2
-
ξ
2
+
k
2
2
-
ξ
2
e
-
i
k
2
2
-
ξ
2
y
+
y
s
e
i
ξ
x
-
x
s
,
for
ξ
<
k
1
,
g
2
x
,
x
s
,
ξ
=
e
-
ξ
2
-
k
1
2
y
+
ξ
2
-
k
2
2
y
s
i
ξ
2
-
k
1
2
+
ξ
2
-
k
2
2
e
i
ξ
x
-
x
s
,
for
ξ
>
k
2
,
e
-
ξ
2
-
k
1
2
y
-
i
k
2
2
-
ξ
2
y
s
i
ξ
2
-
k
1
2
+
k
2
2
-
ξ
2
e
i
ξ
x
-
x
s
,
for
k
1
<
ξ
<
k
2
,
e
i
k
1
2
-
ξ
2
y
-
k
2
2
-
ξ
2
y
s
k
1
2
-
ξ
2
+
k
2
2
-
ξ
2
e
i
ξ
x
-
x
s
,
for
ξ
<
k
1
.In practice, the first term and the second term on the left hand side of (31) are the downward incident waves and the reflected waves, respectively. Equation (31) indicates that the incident waves will be interfered with the reflected waves in the layer with y<0. Furthermore, (32) represents the transmitted waves, which presents that there are only transmitted waves in the layer with y>0. Visual explanations will be given in Figure 4.Figure 4
(a) The real part of the analytical solution, (b) the real part of the solution generated by the refined uniform 25p withhu=3m, (c) the real part of the solution generated by the refined subgridding 25p with (H,h)=(6m,3m), and (d) the real part of solutions in the line of x=120m generated by different methods.
(a)
(b)
(c)
(d)We then consider a two-layered problem (28) with the source point fixed at (300m,-88m), k1=0.0838 and k2=0.2094. The bounded domain [0m,600m]×[-112m,488m] is our area of interest. The corresponding solution is given via (31)~(33). Thus, numerical integrals have to be dealt with to obtain the solution. In this experiment, we apply the Gauss-type quadratures [26] for weakly singular integrals appearing in (31) and (32). This kind of quadrature rules possesses polynomial order of accuracy and enables 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 25p and the refined subgridding 25p including errors, computational times, number of unknowns, and the errors’ decay exponents (briefly, decay exp.). The error is measured inC-norm, which is defined as follows: for any complex vector z=z1,z2,…,zM, (35)zC≔max1≤j≤M⁡zj, where |zj| is the complex modulus of zj. The C-norm measures the error related to the amplitude of the waves. In addition, 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 PML. In addition, let n~ represent the number of all unknowns which lie in the interior domain in the layer with k(x)=0.0838. Each decay exponent is computed by ln⁡(en~/em~)/ln⁡(m~/n~), where en~ denotes the numerical solution’s error in the layer with k(x)=0.0838 for one test. For the refined subgridding 25p, 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.2094. Therefore, h=H/2 holds.In Tables1 and 2, we present the number of all unknowns in the two-dimensional computational domain, the errors and computational time for the refined uniform 25p and subgridding 25p, respectively. We find that the error of the refined subgridding 25p with (H,h) is much smaller than that by the refined uniform 25p with hu=H, while it is a little bigger than that by the refined uniform 25p with hu=h. However, there is a sharp decrease in both the number of unknowns and the computational time for the refined subgridding 25p with (H,h), when compared with the refined uniform 25p with hu=h. It is also found that, for the two methods, the error for the layer with k(x)=0.2094 decreases slowly when we reduce the step sizes. This is because in each test, the point source function produces a solution which is nonsmooth close to the source and the interface between two layers (see Figure 4(a)). In fact, the solution can be subdivided into a near- and far-field contribution. As in the top layer the incident waves are interfered with the reflected waves returned from the bottom of the top layer (see Figure 4(a)), we take the solution for the layer with k(x)=0.0838 as the far-field contribution.Table 1
The performance of the refined uniform 25p. The pair(⋅,⋅) denotes the error for the layer with k(x)=0.2094 and the layer with k(x)=0.0838, respectively.
h
u
12 m
6 m
3 m
1.5 m
n
8281
25921
68121
231361
Errors
(0.2915, 0.0888)
(0.0282, 0.0119)
(0.0205, 0.0075)
(0.0185, 0.0025)
Time
39.7491 s
396.5701 s
2.5210
e
+
003 s
2.5780
e
+
004 sTable 2
The performance of the refined subgridding 25p. The pair(⋅,⋅) denotes the error for the layer with k(x)=0.2094 and the layer with k(x)=0.0838, respectively.
(
H
,
h
)
(12 m, 6 m)
(6 m, 3 m)
(3 m, 1.5 m)
n
11432
36602
100100
Errors
(0.0370, 0.0261)
(0.0208, 0.0077)
(0.0185, 0.0029)
Time
73.3205 s
796.7439 s
5.6267
e
+
003 sWe next consider the decay exponent of the refined subgridding 25p in the far-field contribution. Table3 presents the decay exponent for the refined subgridding 25p. We find that the decay exponent is one. As n~ denotes the number of the unknowns which lie in the layer with k(x)=0.0838, the decay exponent being one indicates that the refined subgridding 25p is second order in accuracy.Table 3
The decay exponent for the refined subgridding 25p on the layer withk(x)=0.0838.
n
~
2393
8984
33965
Errors
0.0261
0.0077
0.0029
Decay exp.
0.9194
0.7395For further comparison, in Figures4(a)–4(c) we present the real part of the solution for the two-layered problem in the domain [0m,600m]×[-112m,488m] by the analytical method, the refined uniform 25p with hu=3m, and the refined subgridding 25p with (H,h)=(6m,3m), 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 phenomenon obeys Snell’s law. It is easy to see that both Figures 4(b) and 4(c) are similar to Figure 4(a), but Figure 4(c) needs less memory and computational time than those needed by Figure 4(b) (see Tables 1 and 2), which verifies the efficiency of the refined subgridding 25p. In addition, almost no boundary reflections exist in Figures 4(b) and 4(c), which confirms the efficiency of the PML’s absorption. Finally, Figure 4(d) presents the real part of the two-layered problem’s solution in the line of x=120m, respectively, by the analytical method, the refined uniform 25p with hu=12m, the refined subgridding 25p with (H,h)=(12m,6m), and the refined subgridding 25p with (H,h)=(3m,1.5m). From this figure, it is obvious that the solution by the refined uniform 25p with hu=12m has a wavelength which is different from the analytical solution. This is called the “dispersion,’’ which leads to “the pollution error’’ that has a global character (see [3]). We find that the solution by the refined subgridding 25p with (H,h)=(12m,6m) seems better. Furthermore, comparing the two solutions obtained by the refined subgridding 25p with different step sizes, we see that the solution obtained by the refined subgridding 25p converges to the analytical solution.Now, we come to our conclusion that when balancing the accuracy against the computational cost, the refined subgridding 25p possesses much improvement over the refined uniform 25p, especially for large-scale layered problems.


## 4. Conclusions

In this paper, we proposed an optimal 25-point finite-difference subgridding scheme to solve the 2D Helmholtz equation with PML. Firstly, subgrids were applied to discretize the computational domain, and a consistent 25-point finite-difference subgridding scheme for the Helmholtz equation with PML was presented, which is second order in accuracy. For the transitional nodes, new refined choice strategies for choosing optimal parameters were then proposed, based on minimizing the numerical dispersion. Finally, numerical experiments were presented to confirm that the refined subgridding 25p is a good choice for the Helmholtz equation with varying wavenumbers, especially for large-scale layered problems.



---

*Source: 1719846-2016-03-20.xml*