We report on simulations of an unsteady three dimensional viscoelastic fluid flow through a model porous medium, employing a finite volume methodology (FVM) with a staggered grid. Boundary conditions at the walls of the porous structures are imposed using a second order immersed boundary method (IBM), allowing for accurate simulations using a relatively coarse grid. We compare the viscoelastic stresses obtained using this new IBM technique with those published in literature and find good correspondence. Next, we applied this methodology to model viscoelastic fluids with a FENE-P constitutive model flowing through closely spaced cylinders. Using periodic boundary conditions, we modeled the flow behavior for Newtonian and viscoelastic fluids for successive contractions and expansions. We observe the presence of counter-rotating vortices in between the closely spaced cylinders. The viscoelastic flow structure is symmetric for lower Deborah (De) number, but onset of an asymmetry occurs after a critical De for an infinite array of cylinders. In the presence of side walls, we observe that the onset of flow asymmetry happens at a much lower De, which can be related to higher viscoelastic stresses normal to the flow direction and larger extensional viscosities which affect the curved streamlines. The three-dimensional flow characteristics for viscoelastic flow at higher De number are quite different in comparison with Newtonian flow behavior.