Dear Cosima,
Regarding solving for thousands of geodesics, a faster alternative would be to solve them as a PDE, i.e. as a dust fluid. This is analogous to solving the Burgers equation (the pressureless Euler equation) for nonrelativistic flows. The Einstein Toolkit is naturally designed to handle such scenarios. But the details depend on what you are trying to simulate.
For example, for timelike geodesics, this could be done by evolving eq. (10) in
https://arxiv.org/abs/1410.7777with p_a given by eq. (8) and H given by eq. (16), setting the specific enthalpy to h=1.
Some caveats:
1) For null geodesics, you can do something similar, but need a different foliation.
2) If you try to use the Valencia formulation to evolve dust, you must evolve the continuity equation for the density \rho, along with the Euler equation for the velocity. This system is ill-posed and unstable for dust (pressure = 0), because you are evolving an equation that is no-longer needed: One needs to drop the continuity equation, eliminate \rho completely from the equations, and evolve only 3 equations for the spatial components of the 4-velocity. Then the system is well-posed (strongly hyperbolic), just like the Burgers equation (whose relativistic generalization you are solving). Eq. (10) in the above paper does this.
3) As Roland Haas pointed out, because you are evolving dust as a pressureless fluid, even though the flow is geodesic, caustics or shocks may form. If that happens, the way to handle it depends on what you are trying to simulate with geodesics. If you need to simulate the actual geodesics, and not dust, you can use the method of characteristics and disregard shock solutions. (For null geodesics you won't have this problem.)
4) If you need to know where each geodesic particle is at a given time, along with eq. (10) above, you can solve the other Hamilton equation for the position x^a conjugate to p_a,
$ dx^a/dt=\partial H / \partial p_a $,
written in Eulerian PDE form. This evolves a 'label' of the initial position of each particle.
For any clarification, please feel free to contact me or Roland Haas.
Best regards,
Haris Markakis
________________________________________________
Charalampos Markakis
Postdoctoral Research Associate
National Center for Supercomputing Applications
University of Illinois at Urbana-Champaign
1205 West Clark St, Room 4022, MC-257, Urbana, IL 61801
markakis@illinois.edu
________________________________________________