The flow of a Newtonian fluid is known to become unstable when the viscosity does not dominate its dynamics. This behaviour has traditionally been characterised by the non-dimensional Reynolds number, which measures the ratio between inertial and viscous forces. However, in some complex fluids, instabilities may be driven by an elastic mechanism that is determined by the evolution of the fluid microstructure. Molecular dynamics simulations offer a methodology for studying the dynamics of molecular fluids at the microscale. Macroscopic-type flow instabilities are examined with novel molecular dynamics simulations of shear flow between two concentric rotating cylinders. The basic flow of a Newtonian fluid bifurcates at a critical Reynolds number within 3% of the theoretical prediction, where beyond this value counter-rotating vortices form in the Taylor-Couette flow configuration. A spontaneous development of waviness in the vortices is observed at higher Reynolds numbers, and further simulations with polymers in solution as the sheared fluid are performed. Molecular dynamics simulations, however, become prohibitively expensive for large macroscopic flows. The present work addresses this problem for the context of planar shear flow of a Newtonian solvent over polymers grafted to a solid substrate, using a new software library developed for performing massively-parallel continuum-molecular hybrid simulations.