We present the first 2D hydrodynamical finite-volume simulations in which dust is fully coupled with the gas, including its back-reaction onto it, and at the same time the dust size is evolving according to coagulation and fragmentation based on a subgrid model. The aim of this analysis is to present the differences occurring when dust evolution is included relative to simulations with fixed dust size, with and without an embedded Jupiter-mass planet that triggers gap formation. We use the two-fluid polar Godunov-type code RoSSBi developed by Surville et al. combined with a new local subgrid method for dust evolution based on the model by Birnstiel et al. We find striking differences between simulations with variable and fixed dust sizes. The timescales for dust depletion differ significantly and yield a completely different evolution of the dust surface density. In general, sharp features such as pileups of dust in the inner disk and near gap edges, when a massive planet is present, become much weaker. This has important implications for the interpretation of observed substructure in disks, suggesting that the presence of a massive planet does not necessarily cause sharp gaps and rings in the dust component. Also, particles with different dust sizes show a different distribution, pointing to the importance of multiwavelength synthetic observations in order to compare with observations by ALMA and other instruments. We also find that simulations adopting fixed intermediate particle sizes, in the range of 10−2 to 10−1 cm, best approximate the surface density evolution seen in simulations with dust evolution.