I am currently studying the book "Riemannian Geometry" by Petersen.
Defintion: Let $(M, g)$ be Riemannian manifold and let $U \subset M$ be an open set. A function $r : U \to \mathbb{R}$ is said to be a distance function if $$|\nabla r| = 1. $$
In his book Petersen introduces the notation:
$$\partial_r = \nabla r.$$
Is it possible to introduce local coordinates $(x, \theta^1, \dots, \theta^{n-1}) $ on $U$ such that locally the points with fixed coordinate $x = r_0$ belong to the level set $r = r_0$ and such that $\partial_x = \nabla r$?
Then we could call the coordinate $x$ as $r$ and the notation of Petersen would be just an identity and everything would much clearer to me.
Yes sure. The level sets of $r$ are locally smooth hypersurfaces, because of the assumption on $\nabla r$. So, locally, you just have to look at the coordinate system defined by first choosing a local coordinate system on one of the level surfaces ($M_0$, say) and then use the normal distance as the $n-$th coordinate function. That this is in fact a coordinate system in a sufficiently close neighbourhood of the level surface follows since you basically restricting the exponentional map to the (normal part of the tangent spaces) of said level surface and since the derivative of the distance function in normal direction does not vanish. (You are looking at $(p, t) \mapsto \exp_p(t n)$ for $t$ in a neighbourhood of $0$ and $p$ taken from a level set and have to show that the differential of this map has full rank).
This is like geodesic normal coordinates not based on a point but based on a hypersurface.
That you get in fact coordinates for which the level sets coincide with the level sets of your distance function follows from the uniquenes of the integral curves of $\nabla r$. (Note that $r$ is constant along the level set from which the construction starts, so $\nabla r$ coincides with the normal up to sign, and since $\nabla r$ has length one it coincides with $n$ even up to length. It's not hard to see this remains true along the normal geodesics to the initial level surface, so the sets $r= const$ coincide locally with the sets of constant geodesic distance to $M_0$).
This construction breaks down once you reach a focal point of the level set or the cut locus.
(I guess you are after a kind of splitting theorem?).