We use a diffusion Monte Carlo method to solve the many-body Schrodinger equation describing fully-heavy tetraquark systems. This approach allows to reduce the uncertainty of the numerical calculation at the percent level, accounts for multi-particle correlations in the physical observables, and avoids the usual quark-clustering assumed in other theoretical techniques applied to the same problem. The interaction between particles was modeled by the most general and accepted potential, i.e. a pairwise interaction including Coulomb, linear-confining and hyperfine spin-spin terms. This means that, in principle, our analysis should provide some rigorous statements about the mass location of the all-heavy tetraquark ground states, which is particularly timely due to the very recent observation made by the LHCb collaboration of some enhancements in the invariant mass spectra of $J/psi$-pairs. Our main results are: (i) the $ccbar cbar c$, $ccbar bbar b$ ($bbbar cbar c$) and $bbbar b bar b$ lowest-lying states are located well above their corresponding meson-meson thresholds; (ii) the $J^{PC}=0^{++}$ $ccbar cbar c$ ground state with preferred quark-antiquark pair configurations is compatible with the enhancement(s) observed by the LHCb collaboration; (iii) our results for the $ccbar cbar b$ and $bbbar cbar b$ sectors seem to indicate that the $0^+$ and $1^+$ ground states are almost degenerate with the $2^+$ located around $100,text{MeV}$ above them; (iv) smaller mass splittings for the $cbbar cbar b$ system are predicted, with absolute mass values in reasonable agreement with other theoretical works; (v) the $1^{++}$ $cbbar cbar b$ tetraquark ground state lies at its lowest $S$-wave meson-meson threshold and it is compatible with a molecular configuration.