The thermal structure of a protoplanetary disc is regulated by the opacity that dust grains provide. However, previous works have often considered simplified prescriptions for the dust opacity in hydrodynamical disc simulations, for example, by considering only a single particle size. In the present work, we perform 2D hydrodynamical simulations of protoplanetary discs where the opacity is self-consistently calculated for the dust population, taking into account the particle size, composition, and abundance. We first compared simulations utilizing single grain sizes to two different multi-grain size distributions at different levels of turbulence strengths, parameterized through the α-viscosity, and different gas surface densities. Assuming a single dust size leads to inaccurate calculations of the thermal structure of discs, because the grain size dominating the opacity increases with orbital radius. Overall the two grain size distributions, one limited by fragmentation only and the other determined from a more complete fragmentation-coagulation equilibrium, give comparable results for the thermal structure. We find that both grain size distributions give less steep opacity gradients that result in less steep aspect ratio gradients, in comparison to discs with only micrometer-sized dust. Moreover, in the discs with a grain size distribution, the innermost (<5 AU) outward migration region is removed and planets embedded in such discs experience lower migration rates. We also investigated the dependency of the water iceline position on the alpha-viscosity (α), the initial gas surface density (ςg,0) at 1 AU and the dust-to-gas ratio (fDG) and find rice α0.61ςg,00.8fDG0.37 independently of the distribution used in the disc. The inclusion of the feedback loop between grain growth, opacities, and disc thermodynamics allows for more self-consistent simulations of accretion discs and planet formation.