The Arithmetic Mean Solver in Lagged Diffusivity Method for Nonlinear Diffusion Equations

Th is paper deals with the solution of nonlinear system arising from finite difference discretization of nonlinear diffusion convection equations by the lagged diffusivity functional iteration method combined with d ifferent inner iterative solvers. The analysis of the whole procedure with the splitt ing methods of the Arithmet ic Mean (AM) and of the Alternating Group Explicit (AGE) has been developed. A comparison in terms of number of iterations has been done with the BiCG-STAB algorithm. Some numerical experiments have been carried out and they seem to show the effectiveness of the lagged diffusivity procedure with the Arithmet ic Mean method as inner solver.


Introduction
We consider a nonlinear diffusion convection equation where the diffusion coefficient, denoted by σ, depends on the solution.
When we use a fin ite difference discretization, this elliptic equation supplemented by a suitable boundary condition, can be transcribed into a nonlinear system of algebraic equations.
We wish to compute a solution of this system of nonlinear equations with a common iterat ive procedure in which the nonlinear term, corresponding to the discretizat ion of the diffusivity σ, may be evaluated at the previous iteration (see [22]). In literature, this approach of nonlinearity lagging in the diffusivity term is denoted as Lagged Diffusivity Fixed Point Iteration or Lagged Diffusivity Functional Iteration.
In Section 2, a model problem described by a nonlinear diffusion convection equation subject to homogeneous Dirich let boundary conditions is presented and a finite difference discretization is described. Then, the lagged diffusivity procedure for the solution of this nonlinear difference system is stated.
Since, a pu rpose here is to re-examin e the lagged diffusivity procedure for solving the system of nonlinear difference equat ions of elliptic type in the context of Parallel Co mputing, the linear difference system that arises at each iterat ion of the lagged diffusiv ity procedure is solved with the iterative splitting methods of the Arith met ic Mean introduced in [20][21] and of the Alternating Group Exp licit (A GE) introduced by Evans (see [1][2][3]).
Thus, the outer iterates of the lagged diffusivity procedure are the approximate solutions of the linear systems computed with an inner iterative solver; a criterion for acceptability of these approximate solutions is given. A stopping rule for the lagged diffusivity procedure is also given.
Section 3 is devoted to remind the Arith met ic Mean method and AGE method for the solution of linear systems with block tridiagonal coefficient matrix.
In Sect ion 4, the convergence of the lagged diffusivity iteration method is analysed under mild and reasonable assumptions imposed on the diffusivity σ using well known standard techniques (see [16]).
In Section 5, numerical experiments show the behaviour of the inner-outer iterations of the procedure. In this section a comparison of the lagged diffusivity iterat ion method with different iterat ive solvers is also presented. The Arithmet ic Mean and the AGE methods are co mpared, in terms of number of iterat ions, with BiCG-STA B method (see [23]).

The Lagged Diffusivity Functional Iteration Method
Consider the model problem described by the nonlinear diffusion convection equation where u= u(x,y) is the density function at the point (x,y) of a diffusion mediu m R, σ =σ(u)>0 is the diffusion coefficient or diffusivity and is dependent on the solution u, q=q(x,y) ≥ 0 is the absorption term, the velocity vector p=(p 1 , p 2 ) T is assumed to be constant and f(x,y) is a real valued sufficiently smooth function.
In the boundary ∂R o f R, equation (1) can be supplemented by a homogeneous Dirichlet boundary condition of the form ( , ) = 0.
(2) In the following, we suppose R to be a rectangular domain with boundary ∂R and we assume that the functions σ, q and f satisfy the "smoothness" conditions: (i) the function σ = σ(u) is continuous in u; the functions q(x,y) and f(x,y) are continuous in x,y respectively; (ii) there exist two positive constants σmin and σmax such that 0 < ≤ ( ) ≤ , uniformly in u; in addition, q(x,y) ≥ q min ≥ 0; (iii) for fixed (x,y)∈R, the function σ(u) satisfies Lipschitz condition in u with constant Γ (uniformly in x, y), The nonlinearity introduced by the u-dependence of the coefficient σ(u) requires that, in general, the solution of equation (1) is approximated by numerical methods. We superimpose on R∪∂R a grid of points R h ∪∂R h ; the set of the internal points R h of the grid are the mesh points ( , ), for i = 1, ..., N and j = 1, ..., M, with uniform mesh size h along x and y directions respectively, i.e. +1 = + ℎ and +1 = + ℎ for i = 0, ..., N, j = 0, ..., M.
In the case of diffusion equation ( = ), the matrix A(u) is also symmetric; then it is a sy mmetric and positive defin ite matrix ( [24, p. 91]).
We remark that while the grid function is defined on the whole mesh region R h ∪∂R h , the vector ∈ ℝ represents the grid function � � defined only on the interior mesh points R h .
Here, we suppose that a solution * of the system (4) exists in ℬ (see Section 4).
For solving the nonlinear system (4) the easiest and maybe the most co mmon method is to lag the nonlinear term in (4) generating an iterative procedure denoted as Lagged Diffusivity Functional Iteration.
With this iterative procedure the nonlinear system (4) can be solved via a sequence of systems of linear equations.
Specifically, g iven a sequence of positive numbers { } such that → 0 as → ∞ and an initial estimate (0 ) of the solution * of the system (4), we generate a sequence of iterates � ( ) �, = 0, 1, ..., with the following rule for the transition from a current iterat ion ( ) to the new iterate ( +1) : • Find an appro ximate solution ( +1) of the linear system with the criterion for acceptability of the solution on the norm � � ( +1) �� ≤ +1 .
Then, the lagged diffusivity procedure is co mposed by an outer iteration that generates the sequence � ( ) � and by an inner iterative solver of the linear system (8). This solver must be particularly well suited for implementation on parallel co mputers.

Iterative Parallel Solution of the Linear Systems
In this section, we remind the block form of the Alternating Group Exp licit (A GE) and of the Arith met ic Mean methods for the solution of the linear system = , (11) when the n×n matrix A has the block tridiagonal form Here, each square block , i = 1, ..., M, is a nonsingular N×N matrix and the blocks and (i = 1, ..., M; . The AGE method consists in considering the following splitting of the matrix A = 1 + 2 , where 1 and 2 are the following matrices Thus, starting from a vector (0) , the AGE method generates a sequence of iterates � ( ) � as follows; for ≥ 0 and = 0, 1, ..., until convergence: is the identity matrix of order n. The A GE method is convergent when the matrices 1 and 2 are symmetric and positive definite with ≥ 0. In this case it is proved that the optimal choice for r is = √ where 0 < ≤ ( 1 ) ≤ , 0 < ≤ ( 2 ) ≤ , and ( ) , = 1, 2, are the eigenvalues of the matrix (see, e.g. [1]).
Furthermore, if the matrix A is irreducib ly (or strictly) diagonally do minant with positive d iagonal elements, > 0 , and nonpositive off diagonal elements, ≤ 0, ≠ , and then the AGE method is convergent. Indeed, from the hypotheses we have that the matrix A is an ℳ -matrix. Set = + and = − , = 1, 2. The choice o f r in (14) yields that the matrices are strictly (or irreducib ly) d iagonally do minant and have positive diagonal elements and nonpositive off diagonal elements, then are ℳ-matrices with −1 ≥ 0. Fro m the hypotheses on the matrix A, the matrices are nonnegative � ≥ 0, = 1, 2�.

Thus, the iteration matrix T of the A GE method is
, can be seen as a splitting method ( , ) , i.e., = − , by , = 1, 2, and = 1 + 2 , then we have the expression (15) fo r −1 . Now, the proof runs as that of the Regular Splitting Theorem in [18, p. 119].
The Arithmet ic Mean method uses the following two splittings of the matrix = 1 + 1 , = 2 + 2 , where 1 and 2 are the following matrices We suppose M even. If M is odd, we can proceed in a similar way.
In the paper [21], it is proved that the block form of the Arith metic Mean method above described, is convergent when the matrix A is: • irreducibly (strictly) d iagonally do minant with positive diagonal entries and nonpositive off diagonal elements with > 0 ( ≥ 0); Here ( ) denotes an eigenvalue of a matrix V and = − , s=1,2, and ̃ is the symmetric positive definite matrix ̃= + ; • symmetric positive defin ite with ≥ 0.
At each iteration k of the A GE method, we have to solve M /2 linear systems of order 2N, i = 1, 3, 5, ..., M -1, to obtain the vector ( +1 /2 ) . The solution of (17) can be seen as a block partit ioned vector , where each block has N components.
These M / 2 systems can be solved simultaneously (in parallel).
Then, in order to obtain the new iterate ( +1 ) of the A GE method, we have to solve, in parallel, M / 2-1 systems as (17) with i = 2,4, 6, ..., M -2 and two linear systems of order N � 1 ′ + � 1 = 1 , � ′ + � = , that can be solved in parallel, as well. The A GE method has an intrinsic parallelism.
Furthermore, we have to solve two linear systems of order N � 1 + where 1 and indicate the first and the last block of the vector [2] and 1 , the corresponding r.h.s. of the second equation of (16).
These systems can be solved in parallel. The Arith metic Mean method introduces also an explicit parallelism in order to increase the degree of multiprogramming, that is the number o f processes that can be executed simultaneously ( [13, p. 87]).
When the system (11) arises from the fin ite difference discretizat ion of the problem (1)-(2), the diagonal blocks of A in (12) are tridiagonal while the sub and superdiagonal blocks are diagonal.
Thus, the systems (17) and (18) can be solved directly as in [9,8] or iteratively, generating a t wo-stage iterative method as in [12]. A direct solution of these systems can be performed by cyclic reduction solvers ( [17, p. 125

Analysis of the Convergence of the Lagged Diffusivity Functional Iteration Method
In this section, we prove the convergence of the sequence � ( ) � generated by the lagged diffusivity iterat ion method for the solution of the system (4) under the smoothness assumptions (i)-(iii).
We define and , i = 0, ..., N +1 and j = 0, ..., M +1, the grid functions defined on R h ∪∂R h and satisfying the Dirich let boundary condition on ∂R h .
For grid functions � � and � � of this type, the discrete 2 ( ℎ ) inner product and norm are defined by the formulas respectively. We say that the grid functions � � defined on R h ∪∂R h and vanishing on ∂R h satisfy Property A if they are uniformly bounded and have uniformly bounded backward difference quotients ∇ and ∇ at each mesh point ( , ) of R h ∪∂R h . The set of all grid functions which satisfy Property A is denoted by ℬ. Thus, ℬ is the set of all grid functions � � for wh ich there exist t wo positive constants and such that ∥ ∥ ℎ ≤ , The constant is independent of h; also is independent of h but it depends on ‖ ‖ (see [16]).
The iterate ( +1) of the lagged diffusivity functional iteration method satisfies the system (8) with the acceptability criterion (9), that is ( +1) is the solution of the linear diffusion convection equation whose diffusivity depends on the previous iterate ( ) with inhomogeneous term − − ( ( +1 ) ).
We assume that all the iterates ( ) satisfy Property A. Thus in particular, by inequality (20), the backward difference quotients of each grid function are bounded and they depend on the inhomogeneous term; then, we have that there exist two constants > 0 and 0 > 0 such that Assume that the mapping ( ) is uniformly monotone in ℬ.
Suppose that { } is a sequence of positive numbers such that → 0 as → ∞. Let ( ) ∈ ℬ be arbitrary and let ( +1) be the solution of ( ) = satisfying the condition (9) with ( ) as in (8).
If all the vectors ( ) belong to ℬ and satisfy Property A with (21) instead of (20), then the sequence � ( ) � converges to * .
Proof. The proof runs as that in [6,Theor. 1]  A more general result on the convergence of the lagged diffusivity functional iterat ion method can be obtained by defining the mapping u=G(u) where ( ) = ( ) −1 , for all ∈ ℝ . A solution of the system (4) is a fixed point of the mapping ( ) .
Since the s moothness condition (iii) on ( ) , we can have that the matrix A(u) satisfies the Lipschitz-continuity condition for every bounded subset Ω of ℝ with a Lipschitz constant . Then we can write, Here, ∥• ∥ indicates an arbitrary vector and matrix norm. The last inequality assures that the mapping ( ) satisfies a Lipschitz condition on a bounded subset Ω of ℝ .

Numerical Experiments
In this section, we consider a numerical experimentation of the lagged diffusivity functional iterat ion method for the solution of the nonlinear system (4) generated by the finite difference discretization above described, of the ellipt ic problem (1)- (2). Indeed, we have to solve the system ( ) = . In these experiments, the vector solution * is prefixed and it is composed by the values of the prescribed function 3 ∶ ( ) = 2 ( 1 + 2 − 2 2 ) ; where, in the case of 2, we have 2_1 for a=3, b=2, c=1; 2_2 for a=1.5, b=0.1, c=0.9 and 2_3 for a=1, b =0.01, c=0.99.
The vector is co mputed as ≡ * = ( * ) * , where the matrix ( * ) of order n, has elements as in (5) and (6) with N=M and = 2 . In all the experiments we have N=256 and 1 = 2 ≡ .
At each iteration , = 0,1, …, of the lagged diffusivity procedure, we have to solve the linear system of order n � ( ) � = , with the splitting method of the Arithmet ic Mean or of the Alternating Group Exp licit described in a previous section (see e.g. [9], [10] for an evaluation of the block form o f the Arith metic Mean method on different parallel architectures and [7] for a description of the Fortran code imp lementing the method). These methods are compared with BiCG-STAB method imp lemented as in [14, p. 50].
In the tables, we report the number of iterat ions * and, in brackets, the total number of iterations of the inner solver , i.e., = � . *
The writing max it. indicates that at a certain iteration , the ma ximu m nu mber of iterat ions of the inner solver has been reached. The maximu m nu mber o f inner iterat ions is set equal to 20000.
The writ ing n.c. indicates that at a certain iteration , the condition (7) is not satisfied.
The symbol "  " close to the number of inner and outer iterations denotes that at a certain iteration , the condition (7) is not satisfied and, in these cases, the lagged diffusivity iteration method generates the iteration by performing a prefixed nu mber (equal to 20) of iterations of the inner iterative solver.
In the tables, we indicate with lag-AM, lag-AGE and lag-B, the lagged diffusivity functional iteration method with the Arith metic Mean, the Alternating Group Exp licit and the BiCG-STA B method, respectively, as inner solver.
Furthermore, we observe that, since +1 decreases, for increasing, as (22) and the lagged diffusivity functional iteration method stops at the iteration * when the criterion for * +1 in (10) is satisfied, we have * +1 = Indeed, in the experiments we obtain * = ⌈ log 2 � 1 � ⌉ .

Conclusions
Fro m the nu merical experiments the following conclusions can be drawn: • the outer residual has the same order of and the error in the discrete 2 ( ℎ ) norm has, in wo rst cases, order ℎ ; • the AM method gives better results when the ratio between the maximu m value of and the smallest component of is small, that is the coefficient matrix of the linear system is strongly asymmetric (see [20]) or the deviation from asy mmetry is decreasing. We define as the deviation fro m asymmetry of a matrix the difference between the Frobenius norms of the symmetric and nonsymmetric parts of the matrix. Furthermo re, we can observe the same behaviour of the AGE method with the one of the AM method, respect to the deviation fro m asymmetry of the coefficient matrix that occurs at each step of the lagged diffusivity procedure; • the lagged diffusivity functional iteration method combined with the AM or the A GE method breaks down when the coefficient matrix, at a certain iteration , is not an ℳ -matrix (n.c.). In some cases, the AGE inner solver, requires a large nu mber of inner iterat ions especially for nearly symmetric matrices; • the behaviour of the residual when we use AM or AGE method as inner solver, is always monotone, except in one case where the lagged diffusivity procedure breaks down (the coefficient matrix is not an ℳ-matrix) but it has been possible to force the convergence by running a few number of iterations of the inner iterat ive solver. The nonmonotonicity of the residual happens at these "forced" iterations. This technique of forcing convergence is successful when condition (7) is "almost satisfied".
• the lagged diffusivity functional iteration method combined with the BiCG-STA B method, when it does not break down, requires a nu mber of inner iterations that is not large and seems to be independent from the deviation fro m asymmetry of the coefficient matrix. We observe that the failure of the lagged diffusivity p rocedure with the BiCG-STA B method, in most cases, happens when the decreasing of the residual is not monotone; • the behaviour of the lagged diffusivity procedure with AM or A GE method as iterat ive solver depends on the choice of the initial vector. The choice (0 ) = in the cases ( ) = 2_2 or ( ) = 2_3 yields to negative values for ( ) for a certain (tipically = 1 or = 2 ) so that condition (7) is not satisfied; that is the initial iterate is too close to a region where a sufficient condition for determining the iterates of the outer procedure is not satisfied. Then, the control on condition (7) is a detection to start from another initial vector.