Anisotropic flow coefficients v(1)-v(5) in heavy ion collisions are computed by combining a classical Yang-Mills description of the early time Glasma flow with the subsequent relativistic viscous hydrodynamic evolution of matter through the quark-gluon plasma and hadron gas phases. The Glasma dynamics, as realized in the impact parameter dependent Glasma (IP-Glasma) model, takes into account event-by-event geometric fluctuations in nucleon positions and intrinsic subnucleon scale color charge fluctuations; the preequilibrium flow of matter is then matched to the music algorithm describing viscous hydrodynamic flow and particle production at freeze-out. The IP-Glasma+MUSIC model describes well both transverse momentum dependent and integrated v(n) data measured at the Large Hadron Collider and the Relativistic Heavy Ion Collider. The model also reproduces the event-by-event distributions of v(2), v(3) and v(4) measured by the ATLAS Collaboration. The implications of our results for better understanding of the dynamics of the Glasma and for the extraction of transport properties of the quark-gluon plasma are outlined.