GRHD: Make `cons2prim` work with 2D velocities

I am glad I have the reftests in place, because this was very delicate to implement correctly. Part of this struggle was also to not introduce extra allocations from boxing of the closures, due to Julia sharing variables between scopes ...

This PR also applies a general clean up to the cons2prim function.