Padi truncation of the thermodynamic perturbation theory is used to calculate the solvation chemical potential of a dipolar solute in a model fluid of dipolar hard spheres. Monte Carlo simulations of the solvation thermodynamics are carried out over a wide range of solute and solvent dipoles in order to address the following major issues: (i) testing the performance of the Padi perturbation theory against simulations, (ii) understanding the mechanism of nonlinear solvation, and (iii) elucidating the fundamental limitations of the dielectric continuum picture of dipole solvation. The Pade form of the solvation chemical potential constructed in the paper agrees with the whole body of simulation results within an accuracy of 3%. Internal energy and entropy of solvation are also accurately described by the perturbation treatment. Simulations show a complex nonlinear solvation mechanism in dipolar liquids: At low solvent polarities the solvation nonlinearity is due to orientational saturation that switches to the electrostriction mechanism at higher dipolar strengths of the solvent. We find that the optimum cavity radius of the Onsager reaction-field theory of solvation depends substantially on solvent polarity. A general method of testing the performance of linear solvation theories is proposed. It shows that the fundamental failure of continuum theories consists in their inaccurate description of the internal energy and entropy of solvation.