We present results for the equation of state of the two-dimensional Hubbard model on an isotropic square lattice as obtained from a controlled and numerically exact large-cluster dynamical mean field simulation. Our results are obtained for large but finite systems and are extrapolated to infinite system size using a known finite-size scaling relation, and are supplemented by reliable error bars accounting for all sources of errors. We establish the importance of examining the decay of spatial spin correlations to determine whether a sufficiently large cluster has been used and with this in mind we present the energy, entropy, double occupancy, and nearest-neighbor spin correlations extrapolated to the thermodynamic limit. We discuss the implications of these calculations on pseudogap physics of the 2D Hubbard model away from half filling, where we find a strong behavioral shift in energy below a temperature ${T}^{*}$ which becomes more pronounced for larger clusters. Finally, we provide reference calculations and tables for the equation of state for values of doping away from half filling which are of interest to cold-atom experiments.