In this talk we shall discuss the monolithic Newton multigrid FEM approach for the simulation of general nonlinear incompressible flow problems. The governing equations arise from model problems with non-isothermal behavior and pressure and shear dependent viscosity, in viscoelastic fluids. The well-known HWNP is overcome with LCR formulation wich guarantees the positive definiteness of the discrete conformation tensor. The coupling between the velocity gradient and the elastic stress, which leads to the restriction for the choice of FEM approximation spaces, and the hyperbolic nature of the problem are handled with edge-oriented stabilization. The resulting linearized system inside of an outer Newton solver is a typical nonsymmetric saddle point problem solved using the geometrical multigrid with a Vanka-like smoother. We validate quantitatively this approach for the well-known benchmark of flow around cylinder for two particular models namely Oldroyd-B and Giesekus.