The frequency response functions (FRFs) are of critical interest to dynamic problems, however, they suffer from computational challenges for viscoelastically damped systems. In this paper, a physics-space-based reduction method is proposed for predicting the FRFs of large-scale viscoelastically damped systems involving the standard linear solid model. A physics-dimension subspace is constructed based on original system matrices and viscoelastic parameters, which can be easily generated by using a recursive manner. A projection basis generation algorithm is then developed to generate a standard orthonormal basis within the physics-dimension subspace. With the help of the standard orthonormal basis and the moment-matching-based reduction method, a physics-space-based reduction method is proposed for efficiently predicting the FRFs of large-scale viscoelastically damped systems. Unlike the widely used state-space reduction method, the reduced system of the proposed method can preserve system’s physical structure so that the physical meaning can be captured. Using both theoretical and numerical analyses, the proposed physics-space-based method is more accurate and efficient than the state-space-based reduction method.