We develop a technique for calculating three-dimensional classical partition functions using projected entangled-pair states (PEPS). Our method is based on variational PEPS optimization algorithms for two-dimensional quantum spin systems, and allows us to compute free energies directly in the thermodynamic limit. The main focus of this work is classical frustration in three-dimensional many-body systems leading to an extensive ground-state degeneracy. We provide high-accuracy results for the residual entropy of the dimer model on the cubic lattice, water-ice $I_h$ and water-ice $I_c$. In addition, we show that these systems are in a Coulomb phase by computing the dipolar form of the correlation functions. As a further benchmark of our methods, we calculate the critical temperature and exponents of the Ising model on the cubic lattice.