Harmonic function with gradient of constant norm in hyperbolic 3-space
The main paper you want to consult is
É. Cartan, Familles de surfaces isoparamétriques dans les espaces à courbure constante, Annali di Mat. 17 (1938), 177–191.
In this paper, Cartan considers the problem of studying the functions $f$ defined on an open set in a space $M$ of constant curvature that satisfy two equations of the form $|\nabla f|^2 = A(f)$ and $\Delta f = B(f)$ for some functions, $A>0$ and $B$, of one variable. He shows that the (hypersurface) level sets of $f$ are what he called isoparametric, i.e., they have all of their principal curvatures constant.
This had been considered somewhat earlier earlier by Levi-Civita for the case $n=3$ in flat space and by Segre for all $n$ (again, in flat space). They showed that, up to rigid motion and homothety, a connected isoparametric hypersurface must be an open subset of $S^k\times \mathbb{R}^{n-k-1}$ for some $0\le k\le n{-}1$.
Cartan showed that, when the ambient sectional curvature is negative, i.e., in hyperbolic $n$-space, then, again, a connected isoparametric hypersurface must either be totally umbilic (i.e., an open subset of either a totally geodesic hyperplane, an equidistant hypersurface, a horosphere, or a hypersphere) or be an open subset of a hypersurface of the form $S^k\times H^{n-k-1}$, i.e., it must be a totally geodesic hypersurface or horosphere or else the set of points of constant distance from a totally geodesic submanifold of lower dimension. Thus, if you did have a solution of $|\nabla f|^2=1$ and $\Delta f = 0$ (even a local one in a neighborhood of $p$) in hyperbolic $n$-space, then $f-f(p)$ would be the oriented distance from the level set $f=f(p)$, and it is easy to calculate that no such function is harmonic, so you have a contradiction.
By the way, there is a much greater variety of isoparametric hypersurfaces in spaces of constant positive curvature. In fact, even though Cartan found many nontrivial examples, and there has been a lot of work since then constructing more examples and classifying them, the classification of these is still not complete.
Such a harmonic function does not exist. Here's a sketch of a proof. Consider a level set $\Sigma_c$ of $f=c$. The gradient $\nabla f$ is a constant length vector perpendicular to $\Sigma_c$ at each point of $\Sigma_c$ (let's assume $|\nabla f|=1$). Also, the vector field $\nabla f$ is divergence-free, since $f$ is harmonic, and thus the flow by $\nabla f$ is volume preserving. Now, flow by the vector field $\nabla f$ for time $t$ takes $\Sigma_c$ into $\Sigma_{c+t}$. Since the vectors are constant length, this gives a local orthogonal coordinate system about $\Sigma_c$, which therefore is Fermi coordinates. So the flow lines of $\nabla f$ are geodesics by the Gauss lemma. Since the flow is volume preserving and $|\nabla f|=1$, it also preserves the area of $\Sigma_{c+t}$. This implies that $\Sigma_c$ is a minimal surface, since the derivative of the variation of area under orthogonal deformation is the trace of the second fundamental form. One then computes that the principal curvatures of $\Sigma_c$ at each point are $\pm 1$. This is because under the orthogonal flow for time $t$, the second fundamental form at time $t$ is determined uniquely by the second fundamental form at time $0$. The only way that it can remain trace $0$ is if the principal curvatures are $\pm 1$. Thus, by Gauss' equation, $\Sigma_c$ is isometric to $\mathbb{H}^2_{-2}$. However, by a result of Doug Moore (generalizing a classic result of Hilbert), there is no isometric immersion of $\mathbb{H}^2_{-2}$ into $\mathbb{H}^3_{-1}$.
Addendum: Given GB's comment, I should add some explanation of why the principal curvatures are $\pm 1$. If we take a region $\Omega \subset \Sigma_c$, and let $A_t$ denote the area of the image of $\Omega$ at time $t$ under the flow, then $\partial A_t/\partial t = \int_\Omega H dA$, where $H$ is the mean curvature. Since $\nabla f$ preserves the area of $\Omega$, we see that $H=0$, so $\Sigma_c$ is a minimal surface.
The second variation of area is given by $\partial^2 A_t/\partial t^2 = \int_\Omega (2 Det(B)-Ric(\nabla f)) dA =0$, where $B$ is the second fundamental form, and $Ric$ is the Ricci curvature. Since $\Sigma_c$ is minimal, $Det(B)=-p^2$, where $p$ is a principal curvature. Also, $Ric(\nabla f)=-2$, since $|\nabla f|=1$. So we get that $p=1$, and the principal curvatures are $\pm 1$.
One can also see more directly that the principal curvatures of $\Sigma_t$ evolve like the curvature of plane hyperbolic curves (this can be seen by considering osculating spheres). The only way that the principal curvatures can remain opposite is if they are $\pm 1$.