Simulation of viscoelastic fluid flows in complex geometry at high Weissenburg (Wi) number is still a challenging problem in computational rheology. In this thesis, parallel computing toolbox available in OpenFOAM has been analysed in details. The scalability of parallel viscoelastic flow solver has been critically evaluated under benchmark flow problems, including 2D and 3D 4:1 contraction flow, 2D flow past a cylinder with 50% blockage, using up to 4096 CPU cores and 55 million computational grids. Areas for further improvements in parallel computational rheology are discussed. A new monitoring and preserving molecular conformation method is proposed to overcome the unphysical artefact problem in simulation of the FENE-CD-JS model under benchmark flow problems. It greatly enhances the stability of parallel viscoelastic flow solver in simulation of high Wi number flows and, for the first time, successfully captures elastic turbulence in flow past a cylinder problem. A new constitutive model, named as FENE-CD2-JS model, is proposed to overcome the existing shortcomings of the original FENE dumbbell model. The model accounts for high order interactions between non-equilibrium polymer chain molecules and could reproduce the asymmetric oscillatory vortex dynamics in the 16:1 contraction flow geometry at high Wi number (up to 48) flow observed in the experiments. For the first time, the analysis of the spatial distribution of non-equilibrium polymer conformation, through the conformational tensor, in strong complex flow and the results of their power-law scaling are also presented.