Hamiltonian Monte Carlo (HMC) samples from an unnormalized density by numerically integrating Hamiltonian dynamics.Girolami \& Calderhead (2011) extend HMC to Riemannian manifolds, but the resulting method faces Hamiltonian integration instability issues for practical usage. While previous works have been tackling this challenge by using more robust metric tensors than Fisher's information metric, our work focuses on designing numerically stable Hamiltonian dynamics.To do so, we start with the idea from Lu et al. (2017), which design momentum distributions to upper-bound the particle speed, and generalize this method to Riemannian manifolds.In our method, the upper bounds of velocity norm become \emph{position-dependent}, which intrinsically limits step sizes used in high curvature regions and, therefore, significantly reduces numeric errors.We also derive a more tractable algorithm to sample from relativistic momentum distributions without relying on the mean-field assumption.