Viscoelasticity is a characteristic of many complex fluids like polymer melts, petroleum, blood, etc. The investigation of viscoelastic flow mechanism has practical significance in both scientific and engineering field. Owing to strongly nonlinear, numerical method becomes a practical way to solve viscoelastic flow problem. In the study, the mathematical model of three-dimensional flow of viscoelastic fluids is established. The planar contraction flow as a benchmark problem for the numerical investigation of viscoelastic flow is solved by using the penalty finite element method with a decoupled algorithm. The multi-mode finitely extensible nonlinear elastic dumbbell with a Peterlin closure approximation (FENE-P) constitutive model is used to describe the viscoelastic rheological properties. The discrete elastic viscous split stress formulation in cooperating with the inconsistent streamline upwind scheme is employed to improve the computation stability. The numerical methods proposed in the study can be well used to predict complex flow patterns of viscoelastic fluids.