The paper presents a finite element scheme combined with a stream line integration scheme to solve viscoelastic fluid flows with multi-mode differential models. The numerical scheme is applied to the flow through axisymmetrical contractions. An eight-mode Giesekus model, fitted for an LDPE melt, has been used. Results are given for the size of the vortices in front of the contraction. These compare well with experimental values up to very high Deborah numbers. Results for entrance corrections and vortex intensities are also presented. The vortex intensity appears to have a maximum as a function of the Deborah number. The velocity profiles in the contraction are over-developed and the velocity on the axis relative to the mean velocity reaches a maximum at the start of the vort! ex. It is argued that the size of the vortex depends on the ratio of elongational properties versus shear properties. This is further supported by a parameter study on a one-mode Giesekus model.