The mathematical models in reservoir simulation are usually discretized into large linear equations, and solving them needs lots of time. Taking into account the mathematical characteristics of the black oil model, a multilevel preconditioning solution method is designed to deal with the algebraic equations in reservoir numerical simulation. Takes into account some of the properties of pressure, saturation, and implicit well variables in flow model, the multilevel preconditioner is comprised of several different iterative methods, such as algebraic multigrid method, Incomplete LU factorization, Gauss-Seidel iteration with downwind ordering and crosswind blocks and et al. The efficiency and robustness of multilevel preconditioner is proved by a million-cell benchmark problem and a real-world matured reservoir with high heterogeneity, high water-cut, geological faults, and complex well scheduling. The numerical results indicate that the proposed method is not only robust with respect to the heterogeneity, anisotropy, and number of wells but also efficient method that can solve large Jacobian system in reservoir simulation quickly and precisely.