Mathematical Problems in Engineering

Volume 2016, Article ID 1719846, 16 pages

http://dx.doi.org/10.1155/2016/1719846

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

^{1}School of Mathematical Sciences, Shandong Normal University, Jinan 250014, China^{2}Guangdong Province Key Laboratory of Computational Science, Sun Yat-sen University, Guangzhou 510275, China^{3}Department of Mathematics, Foshan University, Foshan 528000, China

Received 3 December 2015; Revised 20 February 2016; Accepted 3 March 2016

Academic Editor: Yan-Wu Wang

Copyright © 2016 Tingting Wu et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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