Memory allocation and computations for Laplace ’ s equation of 3-D arbitrary boundary problems

Computation iteration schemes and memory allocation technique for finite difference method were presented in this paper. The transformed form of a groundwater flow problem in the generalized curvilinear coordinates was taken to be the illustrating example and a 3-dimensional second order accurate 19-point scheme was presented. Traditional element-byelement methods (e.g. SOR) are preferred since it is simple and memory efficient but time consuming in computation. For efficient memory allocation, an index method was presented to store the sparse non-symmetric matrix of the problem. For computations, conjugate-gradientlike methods were reported to be computationally efficient. Among them, using incomplete Choleski decomposition as preconditioner was reported to be good method for iteration convergence. In general, the developed index method in this paper has the following advantages: (1) adaptable to various governing and boundary conditions, (2) flexible for higher order approximation, (3) independence of problem dimension, (4) efficient for complex problems when global matrix is not symmetric, (5) convenience for general sparse matrices, (6) computationally efficient in the most time consuming procedure of matrix multiplication, and (7) applicable to any developed matrix solver.


Introduction
Computation iteration schemes of the transformed form of a groundwater flow problem in generalized curvilinear coordinates are studied in this paper.The groundwater flow governing equation [1] der ived from Darcy's Law and equation of continuity is where φ = piezometric head and K = hydraulic conductivity tensor.As the free surface of the unconfined aquifer is especially interested in this research, the governing is transformed in the computational domain for a more accurate solution.Hence, the governing equation becomes [2] ( ) where and are axes of the generalized curvilinear coordinates, k ij Γ is the Christoffel tensor, and g ij is contravariant metric tensor.In addition to the complex of the governing equation in the generalized curvilinear coordinate, the free surface is boundary condition is highly non-linear which can not be simplified except small free gradient [1].
The numerical solution of this governing partial differential equation, subject to the initial and boundary conditions, can be solved by dividing the domain into a collection of points or elemental cells, discretizing the equation to a set of algebraic equations on this collection, and solving the resulting linear system of equations.Traditional element-by-element methods (e.g.SOR) are preferred since it is simple and memory efficient but time consuming in computation.In recent years, conjugate-gradient-like methods were reported to be computationally efficient [3,4].Among them, using incomplete Choleski decomposition as preconditioner was reported to be a much better method for iteration convergence [5].Conjugate gradient like methods were applicable to many fields [6][7][8][9].In this paper, an index method about memory allocation was developed and applied to a 3-dimensional arbitrary boundary problem.Some numerical experiments had been accomplished for ground water mounding problem using conjugate-gradient-like methods in the generalized curvilinear coordinates.

The governing equation and boundary condition
The experimental law of Darcy can be expressed as where q is the specific flux vector.For impressible fluid The governing equation (Eq.( 1)) results from substituting Darcy's law into the equation of continuity.The free surface boundary condition of groundwater includes (1) zero pressure condition and (2) prescribed flux condition.It is assumed that the free surface separates a 100% saturation from a region of zero saturation.Applying the zero pressure condition on the free surface to the definition of the piezometric head gives z = φ (5) which shows that the piezometric head on the free surface is equal to the elevation head.The unsteady free surface boundary condition is where N is the vertical recharge rate, e n is the effective porosity and F is the free surface function.

The generalized curvilinear coordinates
Transforming the problem and solving it in a regular domain is a technique that has been widely used to solve problems with complex shapes of boundaries such as ducts, engine intakes, and complete aircraft or automobiles.The procedure involves generating a grid in the physical domain, then transforming the physical domain (Cartesian coordinate) into the computational domain (generalized curvilinear coordinate), and finally solving the problem in the computational domain.Grid generation for the physical domain is a key step to obtaining accurate, reasonable solutions.In the development of aerospace engineering, conformal mapping has played an important role.The basic idea of conformal mapping is to transform an arbitrary contour into a circle using analytical functions so that problems of potential flow can be solved; this method utilizes the invariance of solutions to Laplace's equation mapped by analytic functions.This technique is suitable and efficient; yet applications to fluid dynamics were rare before 1970's because of the computational difficulties.Conformal mapping was also applied in the field of ground water mechanics.Many closed-form solutions for simplified ground water flow problems were developed before computers were widely used [10].The field of computational fluid dynamics has rapidly developed with the use of computers.Since the computational difficulties can be overcome, the equations of motion can be solved numerically considering pressure, viscous shear and compressibility.To undertake accurate calculations, a wellorganized computational grid is necessary.In the early 1970's, the theory of numerically-generated, boundary-conforming coordinate systems was developed in order to meet the need for a wellorganized grids of arbitrarily-shaped geometry.
For the finite difference method, the imposition of boundary conditions on a domain of complex geometry requires a complicated interpolation of the data on local grid lines [11] and typically, a local loss of accuracy in the computational solution.Such difficulties motivated the introduction of a mapping or transformation from the physical domain (x, y, z) (rectangular coordinate domain) to a computational domain (ξ, η, ζ) (generalized coordinate domain).A distorted region in the physical domain is mapped into a rectangular region in generalized coordinate space (Figure 1).It is assumed that there is a unique, single-valued relation between the physical domain and the generalized coordinate domain; for a three-dimensional problem, these relations are 3 with the inverse relations The relation between the physical domain and the computational domain can be expressed by the well-known Jacobian matrix, J , and its inverse, 1 J − , defined as and The governing equation transform is performed by the relations between physical domain and the generalized curvilinear domain [12][13][14].In tensor notation, a point in Cartesian coordinates is given by ( ) and in generalized curvilinear coordinates is given by ( ) The directed line segment dr, joining two points an incremental distance apart (Figure 2), is written which defines the base vector , i g , in curvilinear coordinates.IMETI 2016 where i g and j δ are covariant base vectors in curvilinear and Cartesian coordinates, respectively.
i g can be written as a linear combination of j δ , are tangent to the coordinate curves but are not of unit length.As the covariant base vectors are not necessarily orthogonal in generalized curvilinear coordinates, contravariant base vectors, i g , are introduced by The contravariant vectors, i g , are perpendicular to the surface i ξ = constant.From Eqs. ( 12) and ( 13), where j i δ is the Kronecker delta, which is 0 if i j ≠ and 1 if i = j.The covariant metric tensor is defined by and the contravariant metric tensor is given by The determinant of ij g is called g and . The cross (vector) product of any two base vectors is where Applying the relations between Cartesian coordinates and the generalized curvilinear coordinates, the governing equation can be derived as Eq. ( 2).
The steady state free surface boundary condition can be transformed from Eq. ( 6) to be This numerical model was verified by Hele-Shaw model and showed good agreements [2].

Numerical implementation and the sparse matrix
The governing equation for ground water mounding domain is a Laplacian type equation.Assuming that every block of the aquifer is homogeneous, the governing equation in the generalized curvilinear domain is Eq. ( 2).In two-dimensional problem, the discretization is carried out by a second-order accurate, nine-point scheme (Figure 3a).Information at point (I, J) is approximated by its surrounding eight points (i.e. points (I-1, J-1), (I-1, J), ... , and (I+1, J+1)).For a three-dimensional problem, a second-order accurate scheme is used as well.Approximate information at (I, J) is from its surrounding 18 points (Figure 3b).
According to the governing equation and boundary condition of the numerical model stated above, a set of linear system of equations can be found Figure 3a.Two-dimensional: second order accurate nine-point scheme.Figure 3b.Three-dimensional: second order accurate 19-point scheme.
For a 5 by 4 grid problem after transformed in the generalized curvilinear coordinates (Figure 4a), matrix A (Figure 4b) is apparently to be sparse and non-symmetric if the problem domain contains two types of soil property.The soil properties are different at either side of node 3, 8, 13, and 18.In the following matrix A, c denotes value from constant head boundary condition, g denotes value from governing equation, b denotes value from no flow or free surface boundary condition, and h denotes value from flow through heterogeneity condition.The numerical implementation flow chart is described in Figure 5. Systems of linear equations of "guessed free surface" can be solved by conjugate-gradient-like methods and its extensions, which were proposed in this research.

The conjugate-gradient-like methods
The conjugate-gradient-like methods, which are based upon the method of descent, did not converge satisfactorily even though the computations were accelerated.Therefore, preconditioning was used.Preconditioning solves the system of linear equations, transformed from Eq. ( 19), in which the eigenvalues of A are more clustered than those of A [15].Two preconditioners were used: diagonal scaling and incomplete Choleski decomposition.Diagonal scaling is very simple and easy to program since the matrices in the preconditioning algorithm are divided by the diagonal matrix of A. In the incomplete Choleski decomposition, a modification by lumping the neglected values of the off-diagonal elements to the diagonal elements is also used [16].The idea of preconditioning is illustrated in the following description.From Eqs. ( 19) and ( 20 Let M = BB T , then M -1 is the preconditioner.B is usually evaluated from A. For example, A=BB T for a symmetric matrix.This matrix is useful because M -1 A is closer to an identity matrix than A (Strikwerda, 1989).Algorithms for conjugate-gradient-like methods can be found in many places [5,17,18].

Results
The conjugate-gradient-like methods were used as a pilot study of solving Eqs. ( 2) and (18).Figure 6 is a sketch of the two-dimensional flow field and boundary conditions for steady recharge and flow through a two-layered aquifer with a horizontal, impervious bottom and finite, constant head lateral boundaries.2) and (18) Computations of the steady state ground water mound location by different methods were compared by the number of iterations for various schemes to solve the system of linear equations to achieve an error less than 10 -6 .The procedure requires an initial guess for the water table elevation, computation of φ in the flow domain, and an iteration scheme.This example had a recharge width of 10 m with its center located 100 m from the left boundary and steadily supplying water at a rate (N) of 0.001 cm/sec to an aquifer with 40 m initial saturated thickness (H 0 ) and 200 m long (L).The hydraulic conductivity of the upper aquifer (K 1 ) was 0.002 cm/sec and lower aquifer (K 2 ) is 0.003 cm/sec, respectively.The computation grid for the problem is 57 by 57.All the computation is set to be double precision and a convergence criterion of 10 -6 is used.Four iterative schemes and the sky-line method were compared for the above example: that IGCR (10) and IGCR(20) were able to solve this non-symmetric matrix while other methods were not suitable (i.e., SOR, GCR(30), DGCR (10), DGCR(20) and DRCR(30)) or were divergent (i.e., ICG and CG) as CG and ICG only works for solving systems of linear equations with symmetric coefficient matrices.It is evident that the application of preconditioner reduces the iterations (comparing GCR(n), DGCR(n) and IGCR(n)).In this case, using incomplete Choleski decomposition showed a great improvement.The conjugate-gradient-like and Sky-Line (SL) methods were used to solve a three-dimensional ground water mounding model.This example used a rate of 0.05 cm/sec recharge on an area 1m 2 for a homogeneous aquifer (K=0.1 cm/sec) with initial thickness 15 m.The aquifer has horizontal dimensions of 28 m by 28 m (i.e., 784 m 2 area).Three methods are tested by comparing their iterations vs. error.ICG is not valid as shown in Figure 8.For error convergence, IGCR(n) is the best (most efficient) method.Since the three-dimensional global matrix has a very wide band width, there are a large number of zero terms stored when SL is used.LU decomposition of SL is no longer efficient because there are many more operations than for IGCR(n).IGCR(n) displays advantages in computational efficiency and memory efficiency over SL for the three-dimensional finite difference model.

Figure 1 .
Figure 1.Correspondence of the physical and generalized coordinate domain[12].

Figure 4b .Figure 5 .
Figure 4b.Matrix A of a 5 by 4 grid domain in generalized curvilinear coordinates.

( 1 )Figure 7 .
Figure 7. Error vs. iterations for various iterative schemes applied on a 57 by 57 grid Figure 7 shows number of iterations vs. their error by various iterative schemes for a 57 by 57 grid problem.It is indicated that IGCR(30) is the best method from error convergence.Figure 7 also shown

Figure 8 .
Figure 8. Error vs. iterations for various iterative schemes applied on a 22 by 22 by 22 grid