Numerical simulations of thermal convection in a rotating spherical shell play an important role in dynamo models. In this paper, we present a parallel numerical model with high performance for the Earth's outer core convection based on a homegrown supercomputer. An approximate factorization fractional method combined with Crank-Nicolson scheme and second-order Adams-Bashford formula is employed for temporal integration. Spatial terms are discretized by a second-order finite volume scheme based on cubed-sphere grid. Two resultant large sparse linear algebraic equations are solved by preconditioned Krylov subspace iterative method. To accelerate convergence rate and improve parallel performance, linear solver is preconditioned with multilevel restricted additive Schwarz preconditioner based on domain decomposition and multigrid. The preconditioner reduces compute time and improve parallel performance of the solvers, which scales well to over ten thousand processor cores. Numerical results are in good agreement with reference solutions of the benchmark Case 0.