Mathematical modelling and numerical simulations of phase separation becomes much more involved if one component is a macromolecular compound. In this case, the large molecular relaxation time gives rise to a dynamic coupling between intra-molecular processes and the unmixing on experimentally relevant time scales, with interesting new phenomena, for which the name â viscoelastic phase separationâ has been coined. Our model of viscoelastic phase separation describes time evolution of the volume fraction of a polymer and the bulk stress leading to a strongly coupled (possibly degenerate) cross-diffusion system. The evolution of volume fraction is governed by the Cahn-Hilliard type equation, while the bulk stress is a parabolic relaxation equation. The system is further combined with the Navier-Stokes-Peterlin system, describing time evolution of the velocity and (elastic) conformation tensor. Under some physically relevant assumptions on boundedness of model parameters we have proved that global in time weak solutions exist. Further, we have derived a suitable notion of the relative energy taking into account the non-convex nature of the energy law for the viscoelastic phase separation. This allows us to prove the weak-strong uniqueness principle and consequently the uniqueness of a weak solution in special cases. Our extensive numerical simulations confirm robustness of the analysed model and the convergence of a suitable numerical scheme with respect to the relative energy.