Research Article  Open Access
Multigrid Method for Solution of 3D Helmholtz Equation Based on HOC Schemes
Abstract
A higher order compact difference (HOC) scheme with uniform mesh sizes in different coordinate directions is employed to discretize a two and threedimensional Helmholtz equation. In case of two dimension, the stencil is of 9 points while in threedimensional case, the scheme has 27 points and has fourth to fifthorder accuracy. Multigrid method using GaussSeidel relaxation is designed to solve the resulting sparse linear systems. Numerical experiments were conducted to test the accuracy of the sixthorder compact difference scheme with Multigrid method and to compare it with the standard secondorder finitedifference scheme and fourthorder compact difference scheme. Performance of the scheme is tested through numerical examples. Accuracy and efficiency of the new scheme are established by using the errors norms .
1. Introduction
The struggle for computing accurate solution using different grid sizes has increased researchers’ curiosity for developing high order difference schemes. Compact finitedifference scheme is widely used in vast area of computational problems, such as the Helmholtz equations and other elliptic equations [1, 2]. We seek highaccuracy numerical solution of the threedimensional Helmholtz equation as follows: where is a cubic solid domain and is a wave number. The above equation is an elliptic partial differential equation. This equation has broad application in physical phenomena, such as elasticity, electromagnetic waves, potential in time harmonic acoustic and electromagnetic fields, acoustic wave scattering, water wave propagation, noise reduction in silencers, radar scattering, and membrane vibration which are governed in frequency domain. The solution and the forcing function are assumed to have the required continuous partial derivatives up to specific orders and are sufficiently smooth. In this study we analyze finitedifference approximation on uniform grids in , , and directions. In this study we also use a constant value of to obtain a scheme with sixth order of accuracy [3]. The secondorder centraldifference operator at , denoted by , is defined as follows: Difference operators and are defined similarly. Equation (1) can be discretized at given grid points , , and as follows: where denotes the truncated terms on the order . Helmholtz equation has been numerically solved by different techniques and different approaches are developed such as the finitedifference method [4], the finiteelement method [5], the spectralelement method [6], and the compact finitedifference method [7].
In finitedifference methods, the stencil of grid points needs to be enlarged, in order to increase the order of accuracy of approximation, but this is not desirable. Helmholtz equation is extensively solved by finiteelement methods, but the disadvantage of finiteelement method is the high computational cost; another disadvantage is the pollution effect; that is, its results are less accurate solution at higher wave number for the given nodes per wavelength. In spectralelement method, it shows that it requires fever grid nodes per wavelength as compared to the finiteelement methods for the Helmholtz equation [6]. But due to the less sparse resultant matrix, compared to the resulting finiteelement matrix, computational time of both methods is the same [6]. Turkel et al. solved Helmholtz equation in 2D and 3D domain for variable wave number [8]. Many iterative techniques for the Helmholtz equation suffer because of their slow convergence, when high frequencies are required. The struggle for fast iterative methods for highfrequency Helmholtz equations becomes the focus of research.
In general to obtain more accurate numerical solution add more nodes and use smaller mesh sizes, which needs more computational time and storage space. Singer and Turkel have conducted an important work in this regard [4]. They used Dirichlet and/or Neumann boundary conditions for the development of a fourthorder compact finitedifference method. Later on a sixthorder finitedifference method was developed by Nabavi et al. for solving Helmholtz equation in onedimensional and twodimensional domain with Neumann boundary conditions [7]. It is shown that sixthorder accuracy is the best that can be achieved for Poisson and Helmholtz equation in 2D case [4]. In this work the basic issue discussed is to develop a sixthorder compact finitedifference method for solving threedimensional Helmholtz equation with Multigrid method. The present study is the first that uses sixthorder compact finitedifference scheme for three dimensions and the designing of specialized Multigrid method. This scheme is different from [1, 7] but the difference is of high order, because the term in 3 dimensions has no counterpart in 2 dimensions. It is also different when is kept constant in [8].
2. Higher Order Compact Scheme
2.1. TwoDimensional Case
Consider the twodimensional Helmholtz equation: In order to achieve an appropriate description, from Taylor series expansion, we have where the grid space is ; adding above expressions (5) and solving for the secondorder derivative which gives using (2) and (6), we have similarly we can find the approximation for the variable . Therefore the central difference scheme for Helmholtz equation can be written as follows: where The higher order derivatives are expressed through mixed derivatives by successively differentiating (4); this process also includes derivative of the forcing function . Applying the appropriate derivatives of (4), we get and putting the above equation in (9), we get in order to find the fourthorder approximation of in above equation, which can be obtained from Taylor series expansion such that putting (12) in (11), we get where is laplace operator; for compact sixthorder approximation, differentiating (4), we have where is biharmonic operator; putting (10) and (14) in (15), we get and using (14), (16), and (12), we have Consequently, the compact sixthorder approximation of the twodimensional Helmholtz equation can be written in modified form which leads to where is the approximate solution of (4), such that . Since using (19), (18) can be written as follows: where Also the righthand side of (18) can be expressed as follows: where therefore (18) can be written as follows: This is the compact sixthorder approximation of twodimensional Helmholtz equation, which can be written in the form , where is a sparse, symmetric, and block tridiagonal matrix.
2.2. ThreeDimensional Case
Consider the threedimensional Helmholtz equation: In order to achieve an appropriate description, from Taylor series expansion, we have and adding the above expressions and solving for the second derivative which gives using (2) and (27), we have similarly we can find the approximation for the variables and . Therefore the central difference scheme for Helmholtz equation can be written as follows: where For simplicity we use and . The higher order derivatives are expressed through mixed derivatives by successively differentiating (25); this process also includes derivative of the forcing function . Differentiating (25) twice with respect to and then solving it for , , and , respectively, we have therefore (30) becomes where is laplace operator; for compact sixthorder scheme of 3D problem consider the following discretization: now (32) becomes
Hence
For compact sixthorder approximation, from (25), we have and also using (31) and (36) in (37), we get and using (19) and (39), we have Finally
Consequently, the compact sixthorder approximation of the threedimensional Helmholtz equation can be written in modified form which leads to the following: where is the approximate solution to (25) and ; also is the biharmonic operator. Since using (43) the lefthand side of (42) is where Also the righthand side of (42) can be expressed as follows: where Thus the compact sixthorder approximation of threedimensional Helmholtz equation in simplest form is as follows: Equation (48) can be written in the form . In this equation we assume that the derivative of the righthand side function can be determined numerically. In case where is not known analytically, then one can approximate the derivatives with high order finite differences. We need only a fourthorder accurate approximation of and a secondorder accurate approximation of and . Fourthorder accurate approximation to is obtained by using .
3. Boundary Condition
When a Dirichlet boundary condition is applied, then the formula in (18) can be used for all interior points. In case of Neumann boundary condition, we develop a sixthorder accurate method for two and threedimensional cases; that is, and . Similar formulas hold in the other directions. To be specific, we consider the coordinate line , in (18) and introduce a ghost point . At the boundary , we specify both the Helmholtz equation and the Neumann boundary condition.
3.1. TwoDimensional Case
Considering (4), using (5), we have From (4), we have In above equation (50), we need fourthorder approximation of ; that is, Making use of (51) in (50), we have Now the secondorder approximation of is and also the secondorder approximation of is Using the derivatives of (4), we have Using (52), (54), and (55) in (49), we have for ; in (56) the approximation used is where is arbitrary constant; using (57) in (56), we have where , − , and . Putting in (24), we have Multiplying (58) by and adding (59), we get formula for boundary nodes where in order to eliminate all the elements with .
3.2. ThreeDimensional Case
Consider (25), by using (26), we have To be specific, we consider the coordinate line and introduce a ghost point ; from (25), we have Using (62) in (63), we have Now using center difference discretization at , (25) will become Equations (60) and (65) yield that for Neumann boundary points the stencil is of five points for twodimensional case while in threedimensional case the stencil has seven points. The matrix is inverted including the extra line . Applying Neumann boundary conditions accuracy and CPU timing remain the same in all examples.
Different methods are used to solve (48) such as finitedifference scheme, fourthorder compact finitedifference scheme, and LU decomposition. We are solving this equation by Multigrid method with sixthorder compact finitedifference scheme using GaussSeidel method as a smoother.
4. Multigrid Method
The results of secondorder finitedifference scheme and sixthorder compact finitedifference scheme are in sparse linear systems; that can be solved efficiently by Multigrid methods. In order to remove high frequency error, Multigrid method utilizes some relaxation methods. Multigrid method makes use of coarse grid correction to smooth the errors. Efficient Multigrid method is implemented in [9, 10] for solving 2D Poisson equation discretized by the standard fourthorder compact difference schemes. It is shown that the sixthorder compact difference scheme is the most cost effective compared to the secondorder finitedifference schemes and the corresponding conventional secondorder centraldifference scheme with Multigrid methods. The standard Multigrid method with a point GaussSeidel type relaxation and standard mesh coarsening strategy does not work very well with unequal mesh sizes discretized Poisson equation [11].
The strategy used is the line relaxation to replace point relaxation. That is, line GaussSeidel relaxation can be shown to be very effective in removing high frequency errors. Consider in this paper for the particular problem the dominant direction is always either the or direction, but not both. Let be the dominant direction; we only perform line relaxation along direction on each successive grid. The coefficient matrix of the sixthorder compact finite difference with this ordering can be written in block tridiagonal matrix of the block order . The order of the coefficient matrix is of , where , where and are both symmetric tridiagonal submatrices of order . They represent the submatrix of of each grid line along one direction, where is the part of the solution vector representing the grid points on each jth line, and is the corresponding part on the righthand side vector. On each level needs only one factorization. The factorization cost of is negligible because matrix B has constant blocks which does not change from one grid line to another grid line.
In Multigrid method with LUdecompositon or GaussSeidel relaxations, we use bilinear interpolation to transfer correction from a coarse grid to a fine grid, and we also use a fullweighting scheme to update the residual on a coarse grid.
Algorithm 1 (Multigrid Algorithm). Assuming that we set up these Multigrid parameters:pre smoothing steps on each level. post smoothing steps on each level. is the number of Multigrid cycles on each level ( for Vcycling and for Wcycling), here we use Vcycle with .
FAS Multigrid Cycle(1)if is the coarsest grid, then solve (48) using a time marching technique and then stop. Else do the presmoothing step:
(2)Restriction: consider
(3)Interpolation: consider
(4)One has
Here the restriction operator is by full weighting and the interpolation is by bilinear operator.
5. Numerical Calculations
In order to show the efficiency and applicability of the Multigrid method, numerical experiments are conducted to solve a twodimensional Helmholtz equation (1) on the unit square domain . The right handside function and the pure Dirichlet boundary conditions are applied on all side of a unit square and unit cubic domain.
5.1. 2D Case
Example 1. Consider the following:
The exact solution is . The sixthorder compact finitedifference scheme is compared against the standard secondorder central difference scheme and fourthorder compact difference scheme, in terms of solution of accuracy, Multigrid convergence rate, and CPU timing. One presmoothing and postsmoothing are applied at each level. The iteration stops when the Euclidean norm (2norm) of the residual vector is reduced by . The maximum absolute error reported is the maximum absolute error between the computed solution and the exact solution over the entire fine grid points. In order to compare the numerical solution and the exact solution, we use norm. The matrix norm of the error vector is defined as follows: where the error vector and the residual , is the number of nodes, is the wave number, and furthermore is the MG with 2nd order, is MG with 4th order, and is MG with 6th order. The norm of the error for and the data in Table 1, Table 2 indicates the behavior of for different values (see Figure 1).
