In both bubbly and porous media flow, the jumps in coefficients may yield an ill-conditioned linear system. The solution of this system using an iterative technique like the conjugate gradient (CG) is delayed because of the presence of small eigenvalues in the spectrum of the coefficient matrix. To accelerate the convergence, we use two levels of preconditioning. For the first level, we choose between out-of-the-box incomplete LU decomposition, sparse approximate inverse, and truncated Neumann series-based preconditioner. For the second level, we use deflation. Through our experiments, we show that it is possible to achieve a computationally fast solver on a graphics processing unit. The preconditioners discussed in this work exhibit fine-grained parallelism. We show that the graphics processing unit version of the two-level preconditioned CG can be up to two times faster than a dual quad core CPU implementation@en