Localized Linear Systems in Implicit Time Stepping for Advection-Diffusion-Reaction Equations

Implicit reservoir simulation models offer improved robustness compared to semi-implicit or explicit alternatives. The implicit treatment gives rise to a large nonlinear algebraic system of equations that must be solved at each time-step. Newton-like iterative methods are often employed in order to solve these nonlinear systems. At each nonlinear iteration, large, sparse linear systems must be solved to obtain the Newton update vector. It is observed that these computed Newton updates are often sparse, even though the sum of the Newton updates over a converged time-step may not be. Sparsity in the Newton update suggests the presence of a spatially localized propagation of corrections along the nonlinear iteration sequence. Substantial computational savings may be realized by restricting the linear solution process to obtain only the nonzero update elements. This requires an a priori identification of the set of nonzero update elements. To preserve the convergence behavior of the original Newton-like process, it is necessary to avoid missing any nonzero element in the identification procedure. This ensures that the localized and full linear computations result in the same solution. As a first step towards the development of such a localization method for general fully-implicit simulation, the focus is on sequential-implicit methods for general two-phase flow. We develop a mathematically sound framework to predict this sparsity pattern before the system is solved. The development first mathematically relates the Newton update in functional space to that of the discrete system. Next, the Newton update formula in functional space is homogenized and solved in such a way that it results in conservative estimates of the numerical Newton update. The cost of evaluating the estimates is linear in the number of nonzero components. Upon projection onto the discrete mesh, the analytical estimates produce a conservative indication of the update's sparsity pattern. The estimates are used to label the components of the solution vector that will be nonzero, and the corresponding submatrix is solved. The computed result is guaranteed to be identical to the one obtained by solving the entire system.;When applied to various simulations of two- and three-phase flow recovery processes in the full SPE10 geological model, the observed reduction in computational effort ranges between four to tenfold depending on the level of total compressibility in the system and on the time step size. We propose, apply, and test a novel algorithm to resolve a system of hyperbolic equations obtained from an Equation of State (EOS) based compositional simulator. When applied to various fully-implicit flow and multicomponent transport simulations, involving six thermodynamic species, on the full SPE10 geological model, the observed reduction in computational effort ranges between four to twelvefold depending on the level of locality present in the system. We apply this algorithm to several injection and depletion scenarios with and without gravity to investigate the adaptivity and robustness of the proposed method to the underlying heterogeneity and complexity. We demonstrate that the algorithm enables efficient and robust full-resolution fully-implicit simulation without the errors introduced by adaptive discretization methods or the stability concerns of semi-implicit approaches.