Perturbation Methods and Formal Modeling for Dynamic SystemsView this Special Issue
Research Article | Open Access
Multigrid Method for Solution of 3D Helmholtz Equation Based on HOC Schemes
A higher order compact difference (HOC) scheme with uniform mesh sizes in different coordinate directions is employed to discretize a two- and three-dimensional Helmholtz equation. In case of two dimension, the stencil is of 9 points while in three-dimensional case, the scheme has 27 points and has fourth- to fifth-order accuracy. Multigrid method using Gauss-Seidel relaxation is designed to solve the resulting sparse linear systems. Numerical experiments were conducted to test the accuracy of the sixth-order compact difference scheme with Multigrid method and to compare it with the standard second-order finite-difference scheme and fourth-order 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 .
The struggle for computing accurate solution using different grid sizes has increased researchers’ curiosity for developing high order difference schemes. Compact finite-difference scheme is widely used in vast area of computational problems, such as the Helmholtz equations and other elliptic equations [1, 2]. We seek high-accuracy numerical solution of the three-dimensional 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 finite-difference 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 . The second-order central-difference 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 finite-difference method , the finite-element method , the spectral-element method , and the compact finite-difference method .
In finite-difference 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 finite-element methods, but the disadvantage of finite-element 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 spectral-element method, it shows that it requires fever grid nodes per wavelength as compared to the finite-element methods for the Helmholtz equation . But due to the less sparse resultant matrix, compared to the resulting finite-element matrix, computational time of both methods is the same . Turkel et al. solved Helmholtz equation in 2D and 3D domain for variable wave number . 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 high-frequency 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 . They used Dirichlet and/or Neumann boundary conditions for the development of a fourth-order compact finite-difference method. Later on a sixth-order finite-difference method was developed by Nabavi et al. for solving Helmholtz equation in one-dimensional and two-dimensional domain with Neumann boundary conditions . It is shown that sixth-order accuracy is the best that can be achieved for Poisson and Helmholtz equation in 2D case . In this work the basic issue discussed is to develop a sixth-order compact finite-difference method for solving three-dimensional Helmholtz equation with Multigrid method. The present study is the first that uses sixth-order compact finite-difference 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 .
2. Higher Order Compact Scheme
2.1. Two-Dimensional Case
Consider the two-dimensional 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 second-order 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 fourth-order 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 sixth-order 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 sixth-order approximation of the two-dimensional 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 right-hand side of (18) can be expressed as follows: where therefore (18) can be written as follows: This is the compact sixth-order approximation of two-dimensional Helmholtz equation, which can be written in the form , where is a sparse, symmetric, and block tridiagonal matrix.
2.2. Three-Dimensional Case
Consider the three-dimensional 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 sixth-order scheme of 3D problem consider the following discretization: now (32) becomes
Consequently, the compact sixth-order approximation of the three-dimensional 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 left-hand side of (42) is where Also the right-hand side of (42) can be expressed as follows: where Thus the compact sixth-order approximation of three-dimensional 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 right-hand 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 fourth-order accurate approximation of and a second-order accurate approximation of and . Fourth-order 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 sixth-order accurate method for two- and three-dimensional 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. Two-Dimensional Case
Considering (4), using (5), we have From (4), we have In above equation (50), we need fourth-order approximation of ; that is, Making use of (51) in (50), we have Now the second-order approximation of is and also the second-order 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. Three-Dimensional 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 two-dimensional case while in three-dimensional 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 finite-difference scheme, fourth-order compact finite-difference scheme, and LU decomposition. We are solving this equation by Multigrid method with sixth-order compact finite-difference scheme using Gauss-Seidel method as a smoother.
4. Multigrid Method
The results of second-order finite-difference scheme and sixth-order compact finite-difference 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 fourth-order compact difference schemes. It is shown that the sixth-order compact difference scheme is the most cost effective compared to the second-order finite-difference schemes and the corresponding conventional second-order central-difference scheme with Multigrid methods. The standard Multigrid method with a point Gauss-Seidel type relaxation and standard mesh coarsening strategy does not work very well with unequal mesh sizes discretized Poisson equation .
The strategy used is the line relaxation to replace point relaxation. That is, line Gauss-Seidel 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 sixth-order 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 right-hand 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 LU-decompositon or Gauss-Seidel relaxations, we use bilinear interpolation to transfer correction from a coarse grid to a fine grid, and we also use a full-weighting 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 V-cycling and for W-cycling), here we use V-cycle with .
FAS Multigrid Cycle(1)if is the coarsest grid, then solve (48) using a time marching technique and then stop. Else do the pre-smoothing 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 two-dimensional Helmholtz equation (1) on the unit square domain . The right hand-side 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 sixth-order compact finite-difference scheme is compared against the standard second-order central difference scheme and fourth-order compact difference scheme, in terms of solution of accuracy, Multigrid convergence rate, and CPU timing. One pre-smoothing and post-smoothing are applied at each level. The iteration stops when the Euclidean norm (2-norm) 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).