The paper presents a parallel algorithm for numerical simulation of the two-phase flow in the digital rock samples. The algorithm is based on a finite-difference approach, including the solution of the Poisson equation with spectral preconditioning. The phase-field method is used as the interface-capturing method. Since the algorithm is oriented on multi-GPU systems, the domain decomposition approach is implemented. To simplify parallelization, we investigate the CUDA unified memory model, which allows the transfer of halo data from different domains between GPUs automatically. The implementation performance with unified memory is examined in the case of the Cahn-Hilliard equation, which is the basis of the phase-field method, and compared with another implementation based on MPI data transfers. Using the developed algorithm, we perform drainage simulation, which is the underlying process in many lab measurements, such as an unsteady-state method for determining relative permeability curves and the centrifuge method for capillary pressure estimation.