Hydrodynamical shocks are a manifestation of the non-linearity of the Euler equations and play a fundamental role in cosmological gas dynamics. In this work, we identify and analyse shocks in the Illustris simulation, and contrast the results with those of non-radiative runs.We show that simulations with more comprehensive physical models of galaxy formation pose new challenges for shock finding algorithms due to radiative cooling and star-forming processes,prompting us to develop a number of methodology improvements. We find in Illustris a total shock surface area which is about 1.4 times larger at the present epoch compared to nonradiative runs, and an energy dissipation rate at shocks which is higher by a factor of around 7.Remarkably, shocks with Mach numbers above and below M ≈ 10 contribute about equally to the total dissipation across cosmic time. This is in sharp contrast to non-radiative simulations,and we demonstrate that a large part of the difference arises due to strong black hole radiomode feedback in Illustris. We also provide an overview of the large diversity of shock morphologies, which includes complex networks of halo-internal shocks, shocks on to cosmic sheets, feedback shocks due to black holes and galactic winds, as well as ubiquitous accretion shocks. In high-redshift systems more massive than 1012M⊙, we discover the existence of a double accretion shock pattern in haloes. They are created when gas streams along filaments without being shocked at the outer accretion shock, but then forms a second, roughly spherical accretion shock further inside.