We present full volume cosmological simulations using the moving-mesh code arepo to study the coevolution of dust and galaxies . We extend the dust model in arepo to include thermal sputtering of grains and investigate the evolution of the dust mass function , the cosmic distribution of dust beyond the interstellar medium , and the dependence of dust-to-stellar mass ratio on galactic properties . The simulated dust mass function is well-described by a Schechter fit and lies closest to observations at z = 0 . The radial scaling of projected dust surface density out to distances of 10 \text { Mpc } around galaxies with magnitudes 17 < i < 21 is similar to that seen in Sloan Digital Sky Survey data , albeit with a lower normalisation . At z = 0 , the predicted dust density of \Omega _ { \text { dust } } \approx 1.3 \times 10 ^ { -6 } lies in the range of \Omega _ { \text { dust } } values seen in low-redshift observations . We find that dust-to-stellar mass ratio anti-correlates with stellar mass for galaxies living along the star formation main sequence . Moreover , we estimate the 850 \mu \text { m } number density functions for simulated galaxies and analyse the relation between dust-to-stellar flux and mass ratios at z = 0 . At high redshift , our model fails to produce enough dust-rich galaxies , and this tension is not alleviated by adopting a top-heavy initial mass function . We do not capture a decline in \Omega _ { \text { dust } } from z = 2 to z = 0 , which suggests that dust production mechanisms more strongly dependent on star formation may help to produce the observed number of dusty galaxies near the peak of cosmic star formation .