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, in one hand, to use different types of elements representing the physics as for instance the capillary forces and mass conservation and, on the other hand, to show the flexibility of the geometric multigrid methods to handle efficiently the resulting linear systems. We numerically solve the V-VP, Vorticity-Velocity-Pressure, and S-V-P, Stress-Velocity-Pressure, formulations of the incompressible Navier-Stokes equations based on the least squares principles using different types of finite elements, conforming, nonconforming and discontinuous of low as well as high 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 efficient solver, geometric multigrid, for several prototypical benchmark configurations (driven cavity, flow around obstacles), and we investigate the effects of pressure jumps for the capillary force in multiphase flow simulations (static bubble configuration).