The aim of this work is to investigate a Galerkin least-squares (GLS) multifield formulation for inelastic non-Newtonian fluid flows. We present the mechanical modeling of isochoric flows combining mass and momentum balance laws in continuum mechanics with an inelastic constitutive equation for the stress tensor. For the latter, we use the generalized Newtonian liquid model, which may predict either shear-thinning or shear-thickening. We employ a finite element formulation stabilized via a GLS scheme in three primal variables: extra stress, velocity, and pressure. This formulation keeps the inertial terms and has the capability of predicting viscosity dependency on the strain rate. The GLS method circumvents the compatibility conditions that arise in mixed formulations between the approximation functions of pressure and velocity and, in the multifield case, of extra stress and velocity. The GLS terms are added elementwise, as functions of the grid Reynolds number, so as to add artificial diffusivity selectively to diffusion and advection dominant flow regions—an important feature in the case of variable viscosity fluids. We present numerical results for the lid-driven cavity flow of shear-thinning and shear-thickening fluids, using the power-law viscosity function for Reynolds numbers between 50 and 500 and power-law exponents from 0.25 to 1.5. We also present results concerning flows of shear-thinning Carreau fluids through abrupt planar and axisymmetric contractions. We study ranges of Carreau numbers from 1 to 100, Reynolds numbers from 1 to 100, and power-law exponents equal to 0.1 and 0.5. Besides accounting for inertia effects in the flow, the GLS method captures some interesting features of shear-thinning flows, such as the reduction of the fluid stresses, the flattening of the velocity profile in the contraction plane, and the separation of the boundary layer downstream the contraction.