In this talk we discuss numerical simulation techniques for incompressible fluids with complex rheology which means that local flow characteristics may differ significantly by several orders of magnitude, for instance due to non-isothermal behavior and pressure, resp., shear dependent viscosity. Such fluids usually include viscoplastic as well as viscoelastic effects which is typical for yield-stress fluids, granular material as well as polymer melts and kautschuk. Corresponding applications are relevant for polymer processing, but include also viscoplastic lubrication, fracking and macro encapsulation. In this talk, we present special discretization and solver techniques in which case the coupling between the velocity, pressure and additional variables for the stresses, which leads to restrictions for the choice of the FEM approximation spaces, and the (often) hyperbolic nature of the problem are handled with special Finite Element techniques including stabilization methods. The resulting linearized systems inside of outer Newton-like solvers are (special) nonsymmetric saddle point problems which are solved via geometrical multigrid approaches. We illustrate and analyze numerically the presented methodology for well-known benchmark configurations as well as protoypical industrial applications for several nonlinear flow models.