We apply the multiscale FENE model to a 3D square-square contraction flow problem and to two 2D benchmark experiments. For this purpose, we couple the stochastic Brownian configuration field method (BCF) with our fully parallelized three-dimensional Navier-Stokes solver NaSt3DGPF. The robustness of the BCF method enables the numerical simulation of higher Deborah number flows for which most macroscopic methods suffer from stability issues. We validate our implementation by investigating the numerical error for a 2D viscoelastic Poiseuille flow that has an analytical solution. Furthermore, we compare the FENE model with the FENE-P closure for a two-dimensional 4: 1 contraction flow. We then compare the results of our 3D simulations with that of experimental measurements from literature and obtain a very good agreement. In particular, we are able to reproduce effects such as strong vortex enhancement, streamline divergence and flow inversion for highly elastic flows. Due to their computational complexity, our simulations require massively parallel computations. To this end, we use a domain decomposition approach with MPI. (C) 2014 Elsevier B.V. All rights reserved.