We present a stable and convergent numerical method to simulate unsteady incompressible viscoelastic flow in two dimensional complex geometry. The incompressible viscous momentum equation, coupled with the multiple mode Giesekus constitutive equation for viscoelastic stress, are used to model viscoelastic fluids. We cast the hyperbolic part of the equations as a conservative finite difference method and leave the elliptic part as a source term computed by implicit methods. A projection method is used to enforce the velocity incompressibility constraint and to update pressure. Irregular computational domains are discretized with a Cartesian grid embedded boundary method. We demonstrate our numerical method by computing flow passing a sphere, flow in a well-rounded contraction channel, and in an abrupt contraction channel, with component Weissenberg numbers exceeding 300, and with a range of relaxation times spanning four decades. Published by Elsevier B.V.