Localized Computation of Newton Updates for General Fully-implicit Reservoir Simulation

Reported observations suggest that Newton updates that are computed during the course of a fully-implicit time step are often sparse. The level of sparsity can vary dramatically from nonlinear iteration to the next, and across time steps. Reported observations suggest that the level of sparsity can be as large as 95%. This work develops an algorithm that accurately predetermines the nonzero elements of the Newton update, and subsequently, can compute it by only solving a truncated linear system. Several alternative ad hoc sparsity prediction strategies have been proposed. Due to their inability to consistently and accurately predetermine the sparsity set, the resulting Newton updates that are computed are inaccurate, leading to a severe degradation of the nonlinear convergence rate. An exact strategy based on an analysis of the sparsity graph of the Jacobian matrix was also proposed for two phase incompressible flow without gravity. Although exact, the proposed strategy cannot be generalized to more complex physics or numerical approximations. Recently, a theoretically sharp and conservative estimate for the sparsity set was derived specifically for the pressure and saturation variables in two-phase sequential-implicit simulation. In this strategy, the discrete Newton update was related to analytical solutions of linear Partial Differential Equations for flow and transport independently. The analytical solutions were evaluated and projected onto the computational domain, thereby providing an estimate of the sparsity set. The theoretically reliable algorithm was demonstrated to reduce the sequential-implicit simulation time for general two phase flow in the full SPE 10 comparative geological model by 5 fold. In this work, the approach is extended to general fully-implicit simulation of coupled flow and multicomponent transport. This is accomplished by considering a canonical functional form of the equations for flow and a system of transported quantities. The analytical estimate is derived by solving the system of linear differential equations using the Schur complement decomposition in functional space. When applied to various simulations of three-phase flow recovery processes in the full SPE 10 model, the observed reduction in computational effort ranged between four and tenfold depending on the level of total compressibility in the system and on the time step size. To investigate the scalability of the algorithm, we applied it to refined models of the SPE 10 case and to multicomponent problems. The improvement in computational speed scales strongly with the number of transport components, and to a lesser degree with problem size.