Least-squares finite element methods are motivated, beside others, by the fact that in contrast to standard mixed finite element methods, the choice of the finite element spaces is not subject to the LBB stability condition and the corresponding discrete linear system is symmetric and positive definite. We intend to benefit from these two positive attractive features, on one hand, to use different types of elements representing the physics as for instance the jump in the pressure for multiphase flow and mass conservation and, on the other hand, to show the flexibility of the geometric multigrid methods to handle efficiently the resulting linear systems. With the aim to develop a solver for non-Newtonian problems, we introduce the stress as a new variable to recast the Navier-Stokes equations into first order systems of equations. We numerically solve S-V-P, Stress-Velocity-Pressure, formulation of the incompressible Navier-Stokes equations based on the least-squares principles using different types of finite elements of low as well as higher order. For the discrete systems, we use a conjugate gradient (CG) solver accelerated with a geometric multigrid preconditioner. In addition, we employ a Krylov space smoother which allows a parameter-free smoothing. Combining this linear solver with the Newton linearization results in a robust and efficient solver. We analyze the application of this general approach, of using different types of finite elements, and the efficiency of the solver, geometric multigrid, throughout the solution of the prototypical benchmark configuration ‘flow around cylinder’.