Proceeding with a multiplicative constant in light of the other 2 answers. I'll assume $u\in\mathbb R^d$, one can make an identification $\mathbb R^2 = \mathbb C$. These things like $f(a)-f(b)=g(a,b)|b-a|$ look like MVT which usually works out.
Set $f_j(u)=|u|^{p-1}u_j$, then $\partial_i f_j=(p-1)|u|^{p-2}\frac{u_j}{|u|}u_i + |u|^{p-1}\delta_{ij}$. Lazy calculations yield
$$|\nabla f|=\|\partial_i f_j\|_{\ell^2_{ij}}\le C\|\nabla f\|_{\ell^\infty_{ij}} \le C[(p-1) |u|^{p-3}\|u_i\|^2_{\ell^\infty_{i}} + |u|^{p-1}] \le Cp|u|^{p-1},$$
since $\|u_i\|_{\ell^\infty_i} \le \|u_i\|_{\ell^2_i} = |u|$. I'd guess a smarter calculation would give $C=1$, instead we have $C=C(d)$.
Now MVT gives
$$ | |u|^{p-1}u - |v|^{p-1}v | \le \sup_{\theta \in [0,1]} |\nabla f(\theta u + (1-\theta) v)| |u-v|$$
So the goal is just to bound the derivative on the line $[u,v]\subset \mathbb R^n$. $p=1$ is trivial. If $1<p<2$, then $s\mapsto s^{p-1}$ is concave, and $0^{p-1}=0$. This implies subadditivity, giving the result:
\begin{align}
|\nabla f(\theta u + (1-\theta) v)|
&= C(n,d,p)|\theta u + (1-\theta) v|^{p-1}
\\
&\le C(n,d,p)(|\theta u| + |(1-\theta)v|)^{p-1} \\
&\le C(n,d,p)(|\theta u|^{p-1} + |(1-\theta)v|^{p-1}) \\
&\le C(n,d,p)(|u|^{p-1} + |v|^{p-1}).
\end{align}
For $p\ge 2$, we can just directly use convexity, like the question body here
and with some more consideration of constants here.