So to make it simple. I currently use a method I found online to compute the gradient of a scalar field of a mesh.
To test how accurate this is, I made a sphere and followed the gradient direction of the geodesic field to construct a line. That is to say, pick $p$ on the sphere and define $\phi(X)$ the function that return the geodesic distance between $X$ and $p$.
By a property of the exponential map described by Manfredo Do Carmo, the geodesic gradient (the tangent) at $X$ uniquely describes the geodesic containing both $X$ and $P$.
In short, if you know the correct tangent at $X$ you can just follow that direction in the exponential map and eventually hit the source $p$.
So I went and coded exactly that:
In that image, the source $p$ is at the intersection of the 3 lines in the bottom part of the image. So my line passes very close to it, but it doesn't hit.
My assumption is, this happens because my gradient approximation is quite linearly based and as such even small errors quickly have big effects, especially as the points grow farther apart.
I am wondering if people know of methods or algorithms that can improve the accuracy of the gradient.
I got the formula for the discrete gradient at page 28 here.