Multiscale methods have been developed as an alternative approach to upscaling and to accelerate reservoir simulation. The key idea of all these methods is to construct a set of prolongation operators that map between unknowns associated with cells in a fine grid holding petrophysical properties and unknowns on a coarser grid used for dynamic simulation. Herein, we extend one such method-the multiscale restricted-smoothed basis (MsRSB) method-to polymer flooding including shear-thinning (and thickening) effects, which gives highly nonlinear fluid models that are challenging to simulate. To this end, we first formulate a sequentially implicit solution procedure for polymer models with non-Newtonian rheology. By treating the implicit velocity dependence of the viscosities in an inner iteration loop, we obtain a formulation that appears to be more robust and stable than the standard fully-implicit approach. We then use a general algebraic multiscale framework to formulate an efficient and versatile multiscale solver. The unique feature of the MsRSB method is how the prolongation operators are constructed. By using restricted smoothing much in the same way as in smoothed aggregation multigrid methods, one gets a robust and flexible method that enables coarse partitions and prolongation operators to be constructed in an semi-automated manner even for highly complex geo-cellular models with high media contrasts and unstructured cell connections. By setting iterative tolerances appropriately, the resulting iterative multiscale solver can be set to compute mass-conservative approximations to the sequential or fully-implicit solution to arbitrary accuracy and hence be used to trade accuracy for efficiency. We first verify the sequential solution procedure and multiscale solver against a well-established commercial simulator on a test case with simple geometry, highly heterogeneous media properties, and strongly nonlinear fluid behavior. Next, the sequential fine-scale and multiscale solvers are validated on a synthetic simulation model of a shallow-marine reservoir. Here, the computational time is dominated by the pressure solves, and 5-8 times speedup is observed when replacing the fine-scale pressure solver by the iterative multiscale method. We also demonstrate the flexibility of the method by applying it to model with unstructured polyhedral cells that adapt to well positions and faults.