Multigrid method based on transformation-free high-order scheme for solving 2D Helmholtz equation on nonuniform grids

High-order compact difference schemes can achieve higher-order accuracy on uniform grids. However, in some cases these may not achieve the desired accuracy. Therefore, we propose a multigrid method based on high-order compact difference scheme on nonuniform grids. We will use interpolation and restriction operators developed by Ge and Cao (J. Comput. Phys. 230:4051-4070, 2011). The suggested scheme has up to fourth-order accuracy. Lastly, some numerical experiments are given to show the accuracy and performance of the proposed scheme.


Introduction
Two-and three-dimensional elliptic partial differential equations (PDEs) play a pivotal role in different fields of science and technology. High-order compact schemes (HOC) are used for the solution of the Helmholtz equation and other elliptic PDEs [, ]. Consider the two-dimensional (D) Helmholtz equation where is a rectangular domain and k is a wave number. The forcing function f (x, y) and the solution u(x, y) have the required continuous differentiability up to a specific order. The equation has many real-world applications like elasticity, electromagnetic waves, acoustic wave scattering, weather and climate prediction, water wave propagation, noise reduction in silencers, and radar scattering. In this paper, we use a finite-difference approximation on nonuniform grids in discrete domain to obtain a scheme up to fourth-order accuracy. We also considered the Helmholtz equation with constant value of k. Equation () has been solved by different techniques such as finite-difference method (FDM) [], fast-Fourier-transform-based (FFT) methods [], finite-element method (FEM) [], the spectral-element method [], compact finite-difference method [], and multigrid methods []. The multigrid method based on HOC schemes is among the most efficient iterative techniques for solving PDEs [, ].
In FDM the number of mesh points will be enlarged to increase the accuracy; however, it will also increase the computational time. The Helmholtz equation is solved by FEM and spectral-element method, but the limitations of these methods are of high computational cost []. Many iterative techniques for the Helmholtz equation suffer due to their slow convergence. The investigation on fast iterative methods to efficiently solve the large algebraic systems arising from high-order difference schemes for PDEs is more attractive. Multigrid methods together with the HOC schemes on uniform mesh sizes are developed in [-]. In most cases where sudden changes occur in a flow, the step sizes have to be rectified over the entire domain. Under these situations, where points are concentrated in the regions of sharp variation, local mesh refinement procedures [, , , -] are necessary, thus dramatically reducing the computational time and computer storage. Ge and Cao [, ] developed a multigrid method with HOC scheme on nonuniform grids for solving D convection diffusion equation and D Poisson equation. This paper is based on approach that an interpolation operator and a projection operator that are suited for a HOC scheme using nonuniform mesh are represented by a transformation-free HOC scheme on nonuniform grids. The main focus in this paper is to develop a multigrid method based on a HOC scheme on nonuniform grids for solving of the D Helmholtz equation. To the best of our knowledge, the D Helmholtz equation is not solved by a multigrid method based on a HOC scheme on nonuniform grids.

HOC scheme on nonuniform grids
Discretization is performed on twodimensional nonuniform gird points. The interval [a  , a  ] is divided into subintervals In the x-direction, consider h x = a  -a  N x , and the forward and backward step sizes are given by , then the grids turn to be uniform. The approximate values of a function u(x, y) at interior grid points (x i , y j ) are represented by u  , and the estimated values of other eight neighboring points are determined by u i , i = , , , . . . , , as in Figure . The Taylor series expansion is performed for appropriate description of a sufficiently smooth function u(x, y) in the given domain at points  and , which are Multiplying equation () by θ bx and () by θ fx , then adding and solving for the second-order derivative, which gives where , and the second-order central difference operator along the x-direction is defined as if θ fx = θ bx = , then equation () reduces to uniform grids of the central difference operator. Hence, the second-order derivative for the x-direction is and the approximation of the second-order derivative for the variable y can be find accordingly. Therefore, the central difference (CD) scheme for the Helmholtz equation can be discretized as where τ  is the truncation error and is defined as The equations H  , H  , H  and L  , L  , L  are defined as If τ  is dropped off from equation (), then the CD scheme for nonuniform grids becomes According to the definition of δ  x , δ  y , the CD scheme can be written as In equation (), only five grid points are involved. From the definition of τ  we can see that when h fx = h bx and h fy = h by , then equation () is of second-order accuracy. In order to improve the accuracy, we consider and Applying the central difference scheme to equation (), we have Similarly, equation () will be Through central difference schemes, the first-and second-order derivatives in equations (), () can be approximated. Now combining equations () and () with equations () and (), the nine-point HOC scheme on nonuniform mesh points for two-dimensional Helmholtz equation () can be written as The coefficients of the LHS in equation () are given as It is easier to know that this scheme has third to fourth order of accuracy from expansion of τ  .

Multigrid method
The multigrid method is one of the most efficient and fastest methods for solving PDEs.
In the multigrid method, the rate of convergence is independent of the mesh size. This method is more effective for solving large-scale sparse linear systems obtained from the discretization of elliptic PDEs [, , -]. The main principle of the multigrid method is to smoothen the error on coarse grid level using basic iterative methods such as Jacobi or Gauss-Seidel method, etc. The multigrid method consists of three important components that are relaxation, restriction, and interpolation operators. These are applied as 'a single iteration of a multigrid cycle comprised of manipulating the error by the application of relaxation method, fixing the residuals on the coarse grid level, solving the error equation on the coarse grid and adjusting the correction of coarse grid up to the fine grid level' . Some specific methods have been applied for the solution of the D and D Helmholtz equations with HOC schemes on uniform grids [-, , , , ]. A full weighting restriction operator and the standard bilinear interpolation operator are used as the inter-grid transfer operators. But in the case of nonuniform grids, these restriction and interpolation operators cannot be used; so new restriction and interpolation operators for nonuniform grids are proposed by Ge and Cao [] by using the area law developed by Liu []. In the following section, we give out the derivation of the two operators for the completeness.

Restriction operator
The principle of developing restriction operator is based on the evaluation of the residuals on the coarse grid level with the use of residuals on the fine grid level. In the multigrid method, Liu developed a law for the restriction of the residual [], known as the area law.
For every point on the coarse grid level, there are corresponding eight fine grid points surrounding it. On the coarse grid, there is a contribution of different degree between the reference grid points and the corresponding surrounding grid points on the fine grids, and a full weighting restriction operator for nonuniform grids is constructed on the base of area law. These points are shown for convenience in Figure . The basic idea for getting the full weighting restriction operator of each grid point is to analyze the weighting coefficients of the residuals. On the coarse grids, the reference point (i, j) of the fine grids have the major contribution to it, so the corresponding weighting coefficient is evaluated by a  /a. At that instant, we noticed that the grid points near the reference point (i, j) have much more contributions than those far away from it. For instance, the weighting coefficient of the point (i + , j) is given by a  /a, that of the point (i -, j) by a  /a, and so on. Now suppose that r ij is the residual at the fine grid point (i, j) and thatr¯i ,j is the corresponding residual at the coarse grid point (ī,j). It is very simple to see that i = ī and j = j; thus, the full weighting restriction operator on nonuniform grids can be written as in []: in which If the step size reduces to equal size, then the total area is divided into sixteen equal small parts by the grid lines and half-grid lines. Denoting the area of each part by a, we obtain that a  = ā, a  = a  = a  = a  = ā, and a  = a  = a  = a  =ā. Due to this situation, the restriction operator will reduced to the full weighting operator on equal mesh sizes []:r¯i ,j =   r i,j + (r i-,j + r i+,j + r i,j+ + r i,j- ) + (r i+,j+ + r i-,j+ + r i+,j- + r i-,j- ) .

Interpolation operator
For the construction of an interpolation operator, we use a similar strategy. We observed that when grid points are shifted from coarse level to the fine level, at that instant, the grids points on the coarse level are the grid points on fine level. These grid points are shifted directly from the coarse grid level to the fine grid level. The interpolation operator is expressed as r i,j =r¯i ,j . Thus, the points on the fine grid are interpolated with their own neighboring points on the coarse level. The formula for error correction along the x-and y-directions are interpolated as [] In case of central grid points, we use four grid points around them on the coarse grid level to interpolate as follows []: When the grid sizes become equal, then the interpolation operator reduces to the bilinear interpolation on equal step sizes []: (r¯i -,j- +r¯i ,j- +r¯i ,j +r¯i -,j ).

Relaxation operator (smoother)
In the multigrid method, the relaxation operator is an important operator. Its work is not to remove the errors, but to damp the high-frequency components of the errors on the present grid level. A simple smoother (Gauss-Seidel relaxation) method can efficiently remove the errors in all directions for simple isotropic problems [, ], but in case of anisotropic and boundary layer problems, the line Gauss-Seidel [, ] and alternating line Gauss-Seidel methods [, -] are shown to be more robust smoothers. In this paper, we use three relaxations to smooth the residuals on each coarse grid such as the line Gauss-Seidel relaxation, natural Gauss-Seidel relaxation, and Red-black Gauss-Seidel relaxation.

Numerical experiments
In order to check the effectiveness of the present method, some problems are chosen. The V-cycle multigrid method is used with zero initial guess, and the process is stopped when the Euclidean norm of the residual vector is reduced by  - on the finest grid level. The effectiveness of the multigrid method with HOC scheme and CD scheme () is presented. The reported errors are the l  -norms of the errors between the computed solution and the exact solution on finest grid. The order of accuracy for a difference scheme is defined as where Error(N  ) and Error(N  ) are the maximum absolute errors approximated for two different grids with N  +  and N  +  points in both direction, whereas N  is half of N  . We use the l  -norm for comparison of the numerical solution and the exact solution, which is defined as where e i,j is the error vector defined as, e i,j = u i,jv i,j , and v i,j is the discrete approximation of u i,j which implies that u i,j = v i,j + O(h)  . First, we use different grid sizes from   to   to compute the accuracy order. The l  -norms of the error and accuracy order for the same value of λ and different values of N, k are presented in Tables  and . We consider the case where N =  and N =  for the accuracy order of the scheme. We also examined the behavior of the scheme for different values of k. The scheme is sensitive for  ≤ k ≤ . If we increase the value of k = ,, then the error does not decrease further. The scheme behaves robustly with respect to the wave number k. However, for any value of N and k, overall, the error does Table 1 The number of multigrid V-cycles with two schemes and different values of 10, 50, 100, 500, 1,000 and λ = -0.9 the boundary conditions are given by the analytic solution, that is, This problem has a steep boundary layer along x = ; therefore, we are using nonuniform grids along the x-axis, which are accumulating near x = , and uniform grids along the y-axis with the following stretching function []: where λ is a stretching parameter and controls the tightness of the grid points in the x-direction. When λ < , more grid points are accumulated to the boundary x =  and to 1.0661e -4 6.3918e -5 4.1235e -5 0.101 -0.4 9.1100e -5 6.1818e -5 3.2561e -5 1.9760e -5 0.559 -0.6 7.1398e -5 4.5502e -5 1.6021e -5 6.2030e -6 0.650 -0.8 4.7341e -5 8.8061e -6 7.6660e -6 3.7210e -7 0.242 -0.9 7.2104e - 5 1.1268e -5 9.9271e -6 9.3180e -7 2.677  the boundary x =  for λ > . If λ = , then the grids reduced to be uniform. When λ = -. and the grid numbers are   , the grid distribution in the xy-plane is shown in Figure . The estimated accuracy and maximum absolute error with different stretching parameter λ are presented in Table . We see that when λ = , the results are very poor. A more accurate solution and order of convergence are obtained from HOC and CD schemes with decreasing stretching parameter λ on nonuniform grids. We observe that when λ = -., the solution obtained with HOC scheme is more accurate, but when λ further decreases to -., the accuracy decreases. This situation is not wondering because putting more grids in the boundary layer area will necessarily cause lack of mesh points in the other regions of the domain. Figure  indicates the configuration of solution in the xy-plane. Table  shows the l  -norm of the error, CPU timing, and the order of accuracy for different stretching parameters λ in problem . It is also obvious from the results that the line Gauss-Seidel relaxation is the most efficient smoother with the least multigrid V-cycle numbers for such type of problems. (a) shows the exact solution, (b) the solution obtained from HOC scheme on uniform grids, (c) the computed solution obtained from a HOC scheme on nonuniform grids, and (d) the computed solution of CD scheme on nonuniform grids.
Example  Consider the PDE with a source term f (x, y), Its analytic solution is The source function is determined by the analytic solution with the boundary layers along x =  and y = . Hence, nonuniform grids along the coordinate directions with accumulation near x = , y =  is used by the following stretching formula: When λ gets closer to , more grids are accumulated near x = , y = . When λ = . and the grids size is   , the grids distribution is given in Figure . Table  indicates the error norms and order of accuracy for different stretching parameters λ for problem . The value of λ changes from . to .. We observe that in nonuniform grids with increasing the stretching parameter λ, more and more grids accumulate into the boundary layers; consequently, more accurate results are obtained from HOC and CD schemes. The rate of convergence continuously increases with the increase of λ. We observe that when λ = ., a considerably most accurate solution is obtained with the HOC scheme, but when λ increases to ., it leads to decrease in accuracy. Figure  shows the contours of the exact solution in the xy-plane. Table  and Table       We observe that the exact solution does not show high variations; therefore, nonuniform grids are not necessary. Uniform grids are used for this problem to check the effectiveness of the multigrid method. The results obtained from the HOC and CD schemes are presented. The reported error is the error norm over the discretized grid points on the finest grid level. Table  lists the number of multigrid V-cycles and the corresponding CPU time in seconds for solving problem  on the   ,   ,   , and   grids. We can see that, for this problem, the multigrid method is very efficient and all the smoothers work well.

Conclusion
In this paper, we have proposed a transformation-free high-order compact finite-difference scheme on nonuniform grids for solution of the D Helmholtz equation to get up to fourth-order accuracy. Furthermore, we have applied the multigrid method based on the HOC scheme on nonuniform grids, which solved the resulting system efficiently. In the case of boundary layer problems with suitable grid stretching ratios, the accuracy is up to fourth order for the HOC scheme and second order for the CD scheme. Numerical results show that the multigrid method with HOC scheme has the required accuracy and is faster than the CD scheme.