Numerical methods: I built a CUDA-accelerated black hole raytracer for my master’s thesis – would love feedback and collaborators
 |
After a long time of work, I'm making my master's thesis project public: a CUDA-accelerated numerical relativity raytracer for Schwarzschild black holes that compares seven different integration methods, including two novel integrators I derived specifically for this problem. My thesis was inspired by some of the results in the Nasa, Orbits, Flight Book 1963 What it does: You can render accretion disks with full gravitational redshift and relativistic Doppler beaming, run real-time webcam gravitational lensing (point your camera at yourself and see light bend around a black hole), and interactively control every physical parameter through a PyQt5 GUI. What might actually interesting in this sub: The geodesic equation d²u/dφ² + u = (3/2) r_s u², which was the Binet-equation for Schwarzschild Blackholes, was solved with seven methods side by side: RK4, Euler, Adams-Bashforth, Adams-Bashforth 4, Adams-Moulton 4, and two I derived myself. The Bowie method (4th-order explicit Taylor-series, no first-order system splitting) and the Obrechkoff method (4th-order implicit with an analytical 2×2 Jacobian and Newton iteration) are new to this problem. The interesting result: at the photon sphere, every standard method eventually diverges outward and that very fast. The Bowie method stays stable much longer and is much faster. The Obrechkoff method spirals inward, which was strange, the only one that does it. Same equation, same initial conditions, fundamentally different error propagation just from the integration scheme. I did not expect that. I'm still unsure why this is happening. The Bowie method held the photon sphere for 7.3+ full orbits (best stability) with angular momentum conservation at 7×10⁻¹³ % relative error. Obrechkoff held 7.4 orbits with energy conservation at 2×10⁻⁹ %. For rendering, both are 6-8x faster than RK4 at the same step size because they need far fewer φ-integration steps for weakly-deflected paths (or in other words they are calculated in r-distances and then re-transformed, and are still faster). What's in the repo: • Full source code (Python, Numba CUDA, PyQt5) • All seven integrators with fixed and adaptive step size variants • Complete accretion disk, Milky Way background, and webcam renderers • Two massive analysis scripts: trajectory_analysis.py (4300 lines) benchmarks every method against analytical Schwarzschild solutions across 38 impact parameters (1k+ plots), and integration_analysis.py produces multi-method comparison grids, difference maps, radar charts, and redshift profiles • My full compiled thesis as a PDF with all derivations, convergence proofs, and truncation error analysis You can run it in three commands:
Requirements: Python 3.10+, CUDA-capable GPU (CUDA needs to be installed already [11-13 versions]) helps but it falls back to CPU. Why I'm posting this: My thesis supervisor passed away at the end of this project, and I never got the chance to publish the novel integrators or get feedback from anybody, or colleagues or a community. I'm releasing this because I think the methods might be genuinely useful, because the Bowie and Obrechkoff integrators outperform RK4 for this class of ODE and I suspect they'd work well for other second-order equations with the same structure. I'd love feedback on the integrators, the rendering approach, the CUDA implementation, or really anything. If you want to take the code in a new direction, please do. I was thinking about doing a PhD in Kerr metrics about different integrators and better visualizations for this problem but I couldn't find a professor for this task in europe. That's why I'm putting it out there and I will just work in the industry from now on. Repo: github.com/al-sca/blackhole-raytracer (https://github.com/al-sca/blackhole-raytracer) See Mr. John Cooks Blog for a reference about the "newly" discovered numerical method by me: https://www.johndcook.com/blog/2025/12/23/bowie-integrator-and-the-nonlinear-pendulum/ Happy to answer questions in the comments. [EDIT] Here an overview of the app, when running the "uv run main.py": submitted by /u/a_new_rusty_crab |
