Research Article | Open Access

Tingting Wu, Zhongying Chen, Jian Chen, "Optimal 25-Point Finite-Difference Subgridding Techniques for the 2D Helmholtz Equation", *Mathematical Problems in Engineering*, vol. 2016, Article ID 1719846, 16 pages, 2016. https://doi.org/10.1155/2016/1719846

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

**Academic Editor:**Yan-Wu Wang

#### 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