We compute rigorously the scaling limit of multi-point energy correlations in the critical Ising model on a torus. For the one-point function, averaged between horizontal and vertical edges of the square lattice, this result has been known since the 1969 work of Ferdinand and Fischer. We propose an alternative proof, in a slightly greater generality, via a new exact formula in terms of determinants of discrete Laplacians. We also compute the main term of the asymptotics of the difference $\mathbb{E}(ε_{V}-ε_{H})$ of the energy density on a vertical and a horizontal edge, which is of order of $δ^{2}$, where $δ$ is the mesh size. The observable $ε_{V}-ε_{H}$ has been identified by Kadanoff and Ceva as (a component of) the stress-energy tensor. We then apply the discrete complex analysis methods of Smirnov and Hongler to compute the multi-point correlations. The fermionic observables are only periodic with doubled periods; by anti-symmetrization, this leads to contributions from four "sectors". The main new challenge arises in the doubly periodic sector, due to the existence of non-zero constant (discrete) analytic functions. We show that some additional input, namely the scaling limit of the one-point function and of relative contribution of sectors to the partition function, is sufficient to overcome this difficulty and successfully compute all correlations.