Quasi-Newton preconditioners for the solution of large nonlinear systems in porous media

Newton's method for the solution of systems of nonlinear equations requires the solution of a number of linear systems with the Jacobian J as the coefficient matrix. When J is large and sparse, for example for problems arising from the discretization of a nonlinear PDE, which is our case as we describe below, the preconditioned Krylov based iterative schemes can be employed for the solution of the linear system, so that two nested iterative procedures need to be implemented [2,3]. A crucial issue, for the reduction of the total linear iterations, is to use efficient preconditioning techniques. In general, ILU-type preconditioners can be employed and calculated at each nonlinear iteration. In this case, a large cost of the calculation of the preconditioners has been paid. In this work preconditioners for solving the linear systems of the Newton method in each nonlinear iteration are studied. The preconditioner is defined by means of a Broyden-type rank-one update of a given initial preconditioner, at each nonlinear iteration. Our sequence of preconditioners built like this, are bounded in the sense that the norm of the matrix obtained from the identity matrix minus the preconditioner times the Jacobian, in each outer iteration [1]. The approach proposed in this paper is to solve the inner systems of the Newton method with an iterative Krylov subspace method, starting with the ILU(0) preconditioner, computed from the initial Jacobian, and to update this preconditioner using a rank one sum. A sequence of preconditioners Pk can thus be defined by imposing the secant condition, as used in the implementation of quasi-Newton methods. With this approach an algorithm has been constructed where the preconditioners are updated until a kmax number of nonlinear iterations. Afterwards a new ILU(0) preconditioner is computed as a restarted procedure. We have applied the obtained algorithm to the solution of the nonlinear system of algebraic equations arising from the discretization of the highly nonlinear equations governing the two phase model in porous media [4,5]. The discretization has been done using linear finite elements (triangles in two and tetrahedra in three dimensions) yielding a system of first order differential equations integrated in time via backward Euler forward difference method [6]. The numerical results show an improvement both in terms of iteration number and CPU time with respect to the ILU(0) preconditioner computed in each nonlinear iteration. In addition these improvements are mainly when the values of the kmax parameter are small.