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.
1. Introduction
In a Cartesian coordinate system, the 2D Helmholtz equation iswhere is the Laplacian, usually represents a pressure field, , called wavenumber, is a positive number, and 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:where in which , , is the angular frequency, and is a function only of defined aswith the dominant frequency of the source, the thickness of PML, the distance from the point inside PML and the interface between the interior region and PML region, and a constant , and 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, . We call it the Helmholtz equation with PML. Moreover, (2) is also valid for variable (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 Section 2, 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 Figure 1(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.
(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 Figure 1(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 and , respectively, represent the step size for the fine grid and the coarse grid, and define and . For each and , we set and . 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 , , 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. Letfor , andfor . Upon introducing two notationswe approximate as Similar treatment with respect to leads to whereBy setting , we construct an approximation for the Laplacian term as We then approximate by a weighted average: whereWe now establish a 25-point difference approximation for the Helmholtz equation as
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 , , and , 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) satisfythen the 25-point finite difference approximation (14) is pointwise consistent with the Helmholtz equation (1) and is a second order scheme.
Proof. Assume that and . It follows from the Taylor theorem thatwhere Therefore, when , , and hold, we have that the left hand side of the 25-point finite difference approximation (14) is equivalent towhere 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 Figure 1(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 with () yield the 25-point finite difference equationwhere the coefficients are given by
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, , , , , , , , and . This indicates that the parameters in the scheme (19) may vary according to the wavenumber . 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 Figure 1(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 , where is the propagation angle from the -axis. The following analysis is based on the assumption that the wavenumber is a positive constant.
Substituting into (19), dividing both sides by the factor , and finally applying the Euler formula lead to the following dispersion equation:where 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,Combining (15) with (23) yieldsInserting (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 , , , , , , , , , , , , , and . 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 notationswe, respectively, let the real part and the imaginary part of the dispersion equation be zero, which lead to the following two equations:
Note that (26) are functions of 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 , where is the wavenumber for the second kind node and is the step size for the fine grid.
Step 2. Let and choose , , where is a given positive integer.
Step 3. Inserting above 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).
Let be the number of gridpoints per wavelength and there holds [24]. According to , 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 , , , and vary smoothly, while the parameter is very sensitive to the number of gridpoints per wavelength , which may be the effect of the ghost nodes.
(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 Figure 1(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 and there are four velocities: from the top, m/s, 2800 m/s, 3000 m/s, and 3500 m/s. A point source is located at the point , where is the Ricker wavelet in the frequency domain defined by whose dominant frequency is . In addition, all the receivers are placed on the top of the ground; that is, they are in the line of . In the following computation, the time sampling is 8 ms and the highest frequency we compute is 60 Hz. For the refined subgridding 25p, and denote the coarse and fine step size, respectively. In this experiment, we let and . Specifically, we only use fine grids for the top layer whose velocity is .
(a)
(b)
(c)
(d)
(e)
(f)
In Figure 3(b), the monofrequency wavefield (real part) for obtained by the refined subgridding 25p with 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 and the refined subgridding 25p with 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 decrease in memory and a 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 and that generated by the refined subgridding 25p with 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 point and source point with , we consider the solution of the Helmholtz equation in a two-layered medium in which satisfieswith continuity conditionswith the wavenumber in which and are given positive numbers satisfying . Without loss of generality, we assume that holds for this problem. We remark that the condition that indicates the velocity in the layer with is slower than that in the layer with , 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),(ii), wherein which is the Hankel function of the first kind with order zero, , andIn 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 . Furthermore, (32) represents the transmitted waves, which presents that there are only transmitted waves in the layer with . Visual explanations will be given in Figure 4.
(a)
(b)
(c)
(d)
We then consider a two-layered problem (28) with the source point fixed at , and . The bounded domain 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 in -norm, which is defined as follows: for any complex vector , where is the complex modulus of . The -norm measures the error related to the amplitude of the waves. In addition, let represent the number of all unknowns in the two-dimensional computational domain in one test; that is, is the sum of the number of the unknowns in the interior area and that of the unknowns in the PML. In addition, let represent the number of all unknowns which lie in the interior domain in the layer with . Each decay exponent is computed by , where denotes the numerical solution’s error in the layer with for one test. For the refined subgridding 25p, and denote the coarse and fine step size, respectively. In this experiment, the coarse grids are refined only once in the layer with bigger wavenumber . Therefore, holds.
In Tables 1 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 is much smaller than that by the refined uniform 25p with , while it is a little bigger than that by the refined uniform 25p with . However, there is a sharp decrease in both the number of unknowns and the computational time for the refined subgridding 25p with , when compared with the refined uniform 25p with . It is also found that, for the two methods, the error for the layer with 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 as the far-field contribution.
We next consider the decay exponent of the refined subgridding 25p in the far-field contribution. Table 3 presents the decay exponent for the refined subgridding 25p. We find that the decay exponent is one. As denotes the number of the unknowns which lie in the layer with , the decay exponent being one indicates that the refined subgridding 25p is second order in accuracy.
For further comparison, in Figures 4(a)–4(c) we present the real part of the solution for the two-layered problem in the domain by the analytical method, the refined uniform 25p with , and the refined subgridding 25p with , 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 , respectively, by the analytical method, the refined uniform 25p with , the refined subgridding 25p with , and the refined subgridding 25p with . From this figure, it is obvious that the solution by the refined uniform 25p with 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 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.
Appendix
We next present choice strategies for optimal parameters for the third kind node’s stencil.
To obtain a second order scheme, we firstly let the parameters in (19) satisfy (15). We then let all the parameters of ghost nodes in Figure 1(d) be zero:Combining (A.1) with (15) yields
We next present two equations related to the real and imaginary part of the dispersion equation (21) for the third kind node. Inserting (A.2) into the dispersion equation (21) to eliminate the parameters , , , , , , , , , , , and and letting the dispersion equation’s real part and imaginary part, respectively, be zero lead to the following two equations:Based on (A.2) and (A.3), we come to the refined choice strategy for the third kind node. For simplicity, we omit the details here.
Furthermore, based on the number of gridpoints per wavelength , some optimal parameters for the third kind of transitional nodes are presented in Figure 2(b). This figure only shows five parameters in (A.3). The parameters , vary smoothly, while the parameters , , and are sensitive to the number of gridpoints per wavelength , which may be the effect of the ghost nodes.
Competing Interests
The authors declare that there are no competing interests regarding the publication of this paper.
Acknowledgments
This research is partially supported by the Natural Science Foundation of China under Grants 11301310, 11501106, 11226303, 10771224, and 11071264, the Youth Foundation of School of Mathematical Sciences at Shandong Normal University, the Training Programme Foundation for Excellent Young Scholar of Guangdong Province of China under Grant YQ2015164, and Guangdong Provincial Government of China through the “Computational Science Innovative Research Team” program.