Alternating Asymmetric Iterative Algorithm Based on Domain Decomposition for 3D Poisson Problem

Poisson equation is a widely used partial differential equation. It is very important to study its numerical solution. Based on the strategy of domain decomposition, the alternating asymmetric iterative algorithm for 3D Poisson equation is provided. The solution domain is divided into several sub-domains, and eight asymmetric iterative schemes with the relaxation factor for 3D Poisson equation are constructed. When the numbers of iteration are odd or even, the computational process of the presented iterative algorithm are proposed respectively. In the calculation of the inner interfaces, the group explicit method is used, which makes the algorithm to be performed fast and in parallel, and avoids the difficulty of solving large-scale linear equations. Furthermore, the convergence of the algorithm is analyzed theoretically. Finally, by comparing with the numerical experimental results of Jacobi and Gauss Seidel iterative algorithms, it is shown that the alternating asymmetric iterative algorithm based on domain decomposition has shorter computation time, fewer iteration numbers and good parallelism.


Introduction
Poisson equation is an elliptic partial differential equation, which frequently appears in many fields such as fluid dynamics, heat transfer, electromagnetics, acoustics, electrostatics mechanical engineering and so on. Many researches on studding the numerical techniques to approximate the solution of Poisson equation have been made in the past few decades. The application of finite difference methods for solving Poisson equation will normally lead to a large, block, and sparse system of equations. Direct methods and iterative methods [1] are normally considered as common approaches for solving such system of equations. Several high precision multigrid and compact difference methods are given in [2][3][4][5][6]. Romao et al. [7,8] provides the Galerkin and least-squares finite element methods in the solution of 3D Poisson equation. In [9,10], the Haar wavelet methods are given. Speyer et al. [11] provide a preconditioned bi-conjugate gradient stabilized method which is efficient, albeit nonmonotonic and convergent.
With the continuous improvement of computer hardware, people are more and more focused on solving large-scale scientific and engineering problems quickly and efficiently on parallel computers. Therefore, people wish to find some direct methods and iterative methods, which have the characteristic of much better solving elliptic equations and easier parallel implementation. In recent years, parallel algorithms are also constantly emerging. Several new parallel methods of direct solution are proposed. P. Valero-Lara and A. Pinelli et al. [12] provide the implementation of a fast solver based on a block cyclic reduction algorithm for the linear systems of a three dimensional separable elliptic problem. And they also study on the parallel characteristics of an algorithm for the direct solution of linear systems with a block-tridiagonal coefficient matrix (BLKTRI problem) [13]. C. P. Stone et al. [14] analyze the performance of a block tridiagonal benchmark. Many authors have given the implementation of scalar tridiagonal solver on GPUs [15][16][17]. Y. Zhang [16] also illustrates several methods to solve tridiagonal systems on GPUs. Because of the direct method to solve large-scale sparse and block diagonal equation systems, when the coefficient matrix is close to singularity, the calculation will often stop or make mistakes. So people also find iterative methods which can be solved by constructing some efficient iterative schemes to approximate the problem itself, so that the iteration can reach a certain accuracy. In [18][19][20][21], a class of efficient parallel finite difference iterative algorithms for Poisson equation were also proposed.
In addition, the domain decomposition method [22] is also a powerful tool for parallel implementation, which studies parallelization from the model level of physical problems. This kind of method can decompose scale problem into small-scale problem and solve serial problem into parallel problem. The explicit-implicit domain decomposition method is proposed by Kuznetsov [23]. Because the numerical boundary conditions on the internal boundary are often not the same as those of the original mathematical model or the corresponding physical problems, different methods to obtain the internal boundary information form different explicit-implicit domain decomposition (EIDD) methods. This leads to the idea of parallel implementation for iterative method based on domain decomposition. In [24], the authors have proposed a kind of finite difference parallel iterative algorithm for two-dimensional Poisson problem, and verified its efficiency and accuracy.
This paper extends to the study of the domain decomposition method for three-dimensional Poisson problem. Several finite difference asymmetric iterative schemes are constructed, and each asymmetric iterative schemes are used to solve the sub-domains alternatively and in parallel; in the processing of inner interfaces, group explicit (GE) method [25,26] is used. The calculation on the whole solution domain is explicit but using the implicit iterative schemes, which greatly avoids the difficulty of solving linear equations and improves the calculation speed and accuracy. When the number of iteration is odd or even, the iterative process of the presented algorithm is given respectively, and a kind of efficient iterative algorithm is established based on domain decomposition for solving three-dimensional Poisson equation.
This paper is outlined as follows. In Section 2, we present several asymmetric iterative schemes. Section 3 gives the alternating asymmetric iterative algorithm. And the convergence and the optimal relaxation factor are obtained in Section 4. In Section 5, we perform the numerical experiments to examine the presented algorithm. Finally we give the conclusion of this paper in Section 6.

Asymmetric Iterative Schemes
Consider the three-dimensional Poisson problem, with the boundary condition, u(x, y, z) = g(x, y, z), (x, y, z) ∈ ∂Ω. where , and ∂Ω is the boundary of the domain Ω. We divide the solution domain Ω into uniform grid, the space step h x = L/l in x direction, h y = M/m in y direction and h z = K/s in z direction. For implicity, the space steps are assumed equal that h x = h y = h z = h. Denote x i = ih, i = 0, 1, ..., l; y j = jh, j = 0, 1, ..., m; z k = kh, k = 0, 1, ..., s; u (n) i,j,k as numerical solution on the nth iteration level at the grid node (x i , y j , z k ). We can give the classical difference discretization in Equation (3) for the 3D Poisson Equation (1), namely, Then we construct eight asymmetric iterative schemes by the difference operator L with the relaxation factor ω as follows, Figure 1 represents the distribution of unknown solution at the (n + 1)th iteration level for the eight asymmetric iterative schemes (5)- (12).

The Domain Decomposition
We can divide the 3D solution domain Ω into multi-subdomains. For simplicity, we use six grid planes x = p, x = p + 1, y = q, y = q + 1, z = r, z = r + 1 to discrete the solution domain Ω into eight sub-domains, and note Ω i , i = 1, 2, ..., 8 as subsets of grid points, while p, q, r are positive integers with Denote π i is the interfaces of the sub-domain Ω i . The sorting order of sub-domains is as follows: the subspaces above z = r + 1 are sorted anticlockwise starting from the upper right sub-domain of the inner layer, and the sub-domains under z = r are sorted anticlockwise starting from the lower right subspace (as shown in Figure 2). The specific description is as follows:

Algorithm Implementation
In this subsection, we provide a new alternating asymmetric iterative (AAI) algorithm based on domain decomposition for 3D Poisson problem (1) and (2). We give different computational processes in each sub-domains at the odd iteration layers and even iteration layers respectively, and use the asymmetric iterative schemes alternatively. The Group Explicit (GE) method is used to solve the inner interfaces, which makes the algorithm to be computed fast and in parallel, and avoids the difficulty of solving large-scale linear equations.

Implementation of Odd Level Iteration
When the iteration number are odd, namely, n = 2a + 1, a = 0, 1, ..., we solve the grid nodes from the boundaries to the inner interfaces step by step, that is, using the asymmetric iterative schemes (5)- (12) to solve the grid points in Ω i , i = 1, 2, ..., 8 respectively: Obviously, the numerical solution can be obtained independently in parallel when the iteration numbers are odd, which saves a lot of computational time compared with the full-implicit iteration case. In addition, although the asymmetric iterative schemes are implicit, the computational process can be transformed into explicit, which can obviously improve the calculation speed and avoid solving large and complex linear equations.
Therefore, it can be seen that the interfaces π = Interfaces I ∪ Interfaces II ∪ Interfaces III. When the interfaces π are solved, the inner grid nodes in the domain Ω can be solved in order like the odd case of iteration numbers in Section 3.2.1. We give the computational procedures in detail as follows.
(1) The solution to Interfaces I We use the asymmetric iterative schemes (5)- (12) to solve the Interfaces I, then the following linear equations can be obtained: where matrices M 1 and N 1 are represented as bellow, Then we just solve the above eight-order sparse linear Equation (14) to obtain the numerical solution of the inter f ace I.

(2) The solution to Interfaces II
The computational procedure of the Interfaces I is depending on the use of GE method based on eight points per group. Similarly, we use the GE method based on four points per group to solve the Interfaces II between the inner boundaries of eight subspaces Ω i , i = 1, 2, ..., 8. Take one group of the Interfaces II for example to illustrate the order of the solution process. Figure 3 gives the direction of the iteration computation. Using the asymmetric iterative schemes (6), (7), (10), (11) to solve the grid nodes (i, q, r), (i, q, r + 1), (i, q + 1, r + 1), (i, q + 1, r)i = p + 2, p + 3, ..., l − 1 (shown in Figure 3), then we can provide the following fourth-order linear equations: where Then the numerical solution of such a set of inner boundary points can be calculated quickly only by solving the fourth order sparse linear Equation (15). In the same way, the points on the other five groups of inner boundary lines are also calculated by Group Explicit method, and we will not represent them one by one here.
(3) The solution to Interfaces III It can be seen from the calculation process of Interfaces I and II we solve the Interfaces III just depending on the group explicit method based on two points a group. The specific calculation process is the same as above (1) and (2), and we do not repeat it.
Finally, taking the above results as the interface conditions, we use the asymmetric iterative schemes different from the schemes at the odd levels to solve the inner points on Ω i , i = 1, 2, ..., 8.
Through the above specific implementation process of the domain decomposition iteration algorithm, we can see that the calculation of numerical solution can transform implicit iteration to explicit calculation no matter on the odd or even iteration layer. Combining with the domain decomposition method, the alternating asymmetric iterative (AAI) algorithm is well performed in parallel.

The Algorithm Convergence
In the last section, we propose a new AAI algorithm based domain decomposition for solving the three-dimensional Poisson problem (1)-(2), which can be written in the following matrix form, while T ω is the iterative matrix of AAI algorithm and b is the right-term. Then we can give the following theorem: Theorem 1. The sufficient and necessary conditions for the convergence of AAI algorithm are as follows: where ρ(T ω ) is the corresponding spectral radius of the iterative matrix T ω .
Consider the eigenvalue problem of Equation (17): Due to the asymmetry of the schemes (5) and (12) in the calculation direction, we take one of the iteration scheme (5) as an example, Firstly, we give the relationship between the eigenvalues µ of Jacobi iterative matrix B and the eigenvalues λ of the AAI iterative matrix T ω . Let V i,j,k be the eigenvectors of the Jacobi iterative matrix, then Taking Equation (21) into Equation (20), we can obtain, where µ = ± λ (λω + 1 − ω) 1 2 .
If λ is the eigenvalue of the matrix T ω , then Equation (24) determines that µ is eigenvalue of the matrix B, which is Jacobi iteration matrix of Poisson equation. On the contrary, if µ is eigenvalue of the matrix B, it can be determined only if there is a relationship between the eigenvalues λ of Jacobi iteration matrix and the eigenvalues µ of the given iteration matrix in the Equation (24).
In fact, if ω µ < ω < 2, Equation (33) is ture clearly; Otherwise 0 < ω < ω µ , and Therefore, ρ(T ω ) < 1, Equation (18) is proved and the presented AAI algorithm is convergent. Obviously, the spectrum radius ρ(T ω ) of the presented iterative matrix depends on the relaxation factor ω, so choosing approximate ω is important to the number of iterations and the convergence rate.
Since the optimal relaxation factor ω opt is obtained for 2D Poisson problem in [25], we can also provide the same computation for 3D case. When ρ(T ω opt ) obtains the minimum where ρ(B) is the spectrum radius of Jacobi iterative matrix, and ε is a positive, sufficiently small number. The optimal relaxation factor ω opt can be theoretically evaluated by Equation (36).

Numerical Experiments
In order to confirm the effectiveness of the AAI algorithm, the following experiments are carried out.
The initial iterative values u (0) i,j,k = 0 (i = 1, 2, · · · , l − 1; j = 1, 2, · · · , m − 1; k = 1, 2, ..., s − 1) is given. (1) Consider the 3D Laplace equation with the boundary condition, u(0, y, z) = sin(y + z); u(1, y, z) = exp( The exact solution of the 3D Poisson problem (37)-(38) is u(x, y, z) = exp( √ 2x)sin(y + z). Let u(x i , y j , z k ) be the exact solution and u (n) i,j,k the nth iterative solution, the errors are calculated in L ∞ -norm as: Moreover, the rate of convergence in space is calculated by where h 1 , h 2 are the space steps. Table 1 gives the errors E ∞ of the presented alternating asymmetric iteration algorithm based on domain decomposition for the 3D Laplace problem (37)-(38) with different values of ω when l = m = s = 31, h = 1/30, n = 150, we can obviously see that the errors is relatively smaller when the relaxation factor ω is about 1.9. we further see that the errors get the minimum when ω is about 1.82 shown in Figure 4a, which is match with the result of Equation (35). Figure 4b performs the errors with z = 0.5 when ω = 1.82, which illustrate the effectiveness of the AAI algorithm.
From Tables 2 and 3, we can see the iteration numbers of the AAI algorithm is the least during the Jacobi, Gauss-Seidel iterative methods under some error controls when h = 1/30, 1/50. In addition, the AAI algorithm obtains shorter times than the Jacobi and the Gauss-Seidel methods when the number of the grid nodes is in increasing. Table 4 gives the convergence rates and errors of the AAI algorithm. In computation of the rates of convergence in space, the spatial steps are taken as h = 1/(16 + 8d), d = 0, 1, · · · , 4. We can see the rates is of order 2 in space and the errors can up to 10 −5 . Figure 5 provides the errors, relative errors and numerical solutions at z = 1/3, 2/3 when l = m = s = 31, h = 1/30, n = 150, ω = 1.82, which shows the AAI algorithm is effect and accurate.          (2) Consider the 3D Poisson equation with the boundary condition, u(0, y, z) = 0; u(1, y, z) = sin(1)sin(y)sin(z); u(x, 0, z) = 0; u(x, 1, z) = sin(x)sin(1)sin(z); u(x, y, 0) = 0; u(x, y, 1) = sin(x)sin(y)sin (1).
The exact solution of the 3D Poisson problem (40)-(41) is u(x, y, z) = sin(x)sin(y)sin(z). Table 5 gives the errors E ∞ of AAI algorithm for the problem 2 with the different values of ω when l = m = s = 31, h = 1/30, n = 150, and Figure 6 shows the errors get nearly the minimum 2.2319 × 10 −4 while ω = 1.76. which show the effect of ω to the AAI algorithm. Tables 6-8 give the iteration numbers and times under some error controls, which also show obviously the presented algorithm has smaller iteration numbers than the Jacobi and Gauss-Seidel methods. The computational times are shorter with the grid points increasing. Figure 7 shows the errors, relative errors and numerical solutions of the AAI iterative algorithm based on domain decomposition for the problem 2 at z = 1/3, 2/3 when l = m = s = 51, h = 1/50, n = 150, ω = 1.96. All of the numerical experiments examine the effectiveness and accuracy of the presented AAI algorithm.    Since the presented AAI algorithm based on domain decomposition is constructed by the asymmetrical iterative schemes and GE method, which has interior parallelism and is easy to be implemented. During the process of the implementation of the AAI algorithm, there's no need to solve the large-scale sparse block tridiagonal matrices. When the iteration numbers are odd or even, we just to solve the 3D problems by the constructed iterative schemes, which are computed independently. Once the interfaces are solved by GE method, the other grid points can also be solved by the constructed iterative schemes directly and in parallel. Therefore, the key of parallel implementation lies in information transfer and time cost of the inner boundary. In fact, the whole computation is explicit but convergent, which save the most of the consuming time.
So In this paper, we use the Matlab software to implement the presented AAI algorithm. We extend this idea to solve the other time-dependent high-dimensional problems, and compare the times, speedup, caches and so on. The detailed parallel implementation and performance analysis are provided in [27].

Conclusions
In this paper, we provide a new alternating asymmetric iterative (AAI) algorithm for 3D Poisson problem based on domain decomposition. We use several asymmetrical iterative schemes to solve the sub-domains respectively. Meanwhile the asymmetrical iterative schemes are alternatively used on odd and even iteration levels to improve the accuracy. Moreover, we give the convergence of the algorithm and the optimal relaxation factor. Finally, several numerical experiments are taken to examine the effectiveness and accuracy of the presented algorithm.
The study will be extended to other high-dimensional diffusion problems and wave problems and so on, and also can be used to solve on more multi-subdomains, and the corresponding new algorithms will be designed. we will report these soon.