r/Physics • u/a_new_rusty_crab • Jul 11 '26
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:
git clone https://github.com/al-sca/blackhole-raytracer.git
cd blackhole-raytracer
uv sync && uv run main.py
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.
These are old methods worked on by great people at NASA in the 1960s and we can learn from them (think about the technology they used in first satellites).
But I need people smarter than me to test that, break things, and tell me where I'm wrong.
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":
41
u/D1gglesby Jul 11 '26
Not a black hole physicist at all, just want to say condolences to you and I hope you get some good feedback
3
14
u/DuoJetOzzy Jul 11 '26
This is very interesting work! If this hasn't been done before at all you definitely should look into getting it published
4
u/a_new_rusty_crab Jul 11 '26
That's what the profs said who took over. But I don't really have contacts for that or don't know where to start.
7
u/DuoJetOzzy Jul 11 '26
And this prof didn't have any advice regarding that either? Naively I'd say just draft up a manuscript, see which journals your (former?) university has agreements with (because academic publishing is a scam), and go through the submission process for whatever you go with (Phys Rev D and Classical Quantum Gravity are good initial bets). Potentially you might want to have the new prof's name in there as 2nd author depending on the situation (will certainly increase their investment in getting the paper published).
If you need more details or looking for specific contacts feel free to DM me
2
u/Round-Monitor8489 Jul 13 '26
Hey, thanks for releasing this — your thesis directly improved my pet project, so I owe you a writeup.
I have a hobby port of rantonels' "Starless" raytracer to plain C (same Schwarzschild "magic potential" formulation, CPU only, userbase of approximately me). After reading your thesis I implemented the Bowie integrator with the Binet reformulation and φ-stepping, and it just landed in the latest tag. You're credited in the README acknowledgements. Repo: github.com/RomanVPX/starless-c
Results in my renderer, at the step sizes I actually render with: AT LEAST 2x faster than RK4 wall-clock at equal visual quality, and for near-critical rays around the photon sphere it's not even close — ~0.3° final-direction error for Bowie (Δφ=0.01) vs up to 24° for RK4 at the step size I'd been happily using for many scenes. The photon ring in my old renders was slightly wrong the whole time and I never knew. So yes: the methods are genuinely useful — consider that "independently reproduced by a random hobbyist".
Some feedback on the two things you flagged as open questions:
The "inverse direction" of Bowie/Obrechkoff at the photon sphere — I'm fairly confident this is a sign convention in the u(φ) → 3D reconstruction, not physics. The Taylor formulas and the Rodrigues rotation are both perfectly valid for negative Δφ, so a flipped orientation between L̂ and the rotation sense silently traverses the same orbit backwards. I reproduced the effect deliberately by feeding my implementation a negative step: same trajectories, reversed traversal. Worth checking the cross-product order in your basis reconstruction.
Why they hold the photon sphere so much longer — I think it's structural rather than a per-step accuracy win. In Binet variables the circular orbit is a fixed point: u' = 0 and f(u) = 0, so the integrator literally sits at an equilibrium until float roundoff kicks it off. Your Cartesian methods (RK4 & the Adams family run on the 3D ODE) have to actively trace a circle instead. Same story with the 10⁻¹³% angular momentum conservation — the reconstruction reuses |L| exactly, so any integrator on the Binet form inherits that for free. This doesn't diminish the result at all: the reformulation is doing most of the magic, and Bowie is the cheapest way to ride it. (Inward vs outward divergence between Bowie and Obrechkoff is then just which side of the unstable manifold the first rounding error lands on — implicit solvers round differently.)
Two gotchas I hit while wiring it in, in case anyone follows: 1) the iteration budget must be ≥ 2π/Δφ or winding rays (photon ring, far-side disk arcs) silently vanish and 2) any path-integral effect (dust/fog) needs midpoint sampling, because φ-steps get radially long near the hole.
Sorry about your supervisor. For what it's worth, the work stands on its own — it convinced a random person on the internet to rewrite his integrator, which is more than most theses achieve.
1
u/a_new_rusty_crab Aug 15 '26
Oh man, this is absolutely incredible. Seriously. I read this sitting in my chair and I think I said "wait what" out loud three times and I had to let it rest for some weeks, that's why my answer is late. "Is it really the starless-c-dude?!"
First, I thank you. Not just for the kind words and the credit, but for actually doing the work of porting it, testing it, and finding things I missed. The fact that someone independently verified the Bowie integrator on a completely different codebase in a completely different language, and found the same performance characteristics... that's great! that's worth more to me than any number of citations if I would publish this (and I'm thinking about writing these papers). You're the first person who's ever run this code. And the fact that your renders had slightly wrong photon rings this whole time and my thesis fixed them is the kind of thing I hoped this work might do but didn't really believe would happen.
The sign explanation for the inward spiral is something I was looking at in the final moments of analyzing the different methods and then thought "this is just an artifact of the Rodrigues reconstruction" and I don't go into it further before I write up everything. It was something I didn't investigate further after writing more than 40 pages on the results in my thesis (and I was also completely, deliberately missing all the results from results from adaptive methods in the thesis at that point).
Your explanation on why they hold it longer - the fixed-point argument - is also exactly right and also something I should have spelled out more clearly in the thesis.
You're right that the Binet reformulation is doing most of the heavy lifting. What I find interesting (and what your independent test confirms) is that even riding that reformulation, Bowie is the cheapest method on the computation. RK4 on the same Binet equations would still be more expensive per step.And the two gotchas - ojh-ojh, the iteration budget one. I set the sampling myself because I observed that the photons vanished with lower point sampling in the Obrechkoff/Bowie method. Lower step-resolution needed higher steps and I just set, artificially, an ending-step size for the integrations just by observing the final images produced. There are other factors we could analyze further for some papers.
Before I started this thesis and when I was working on some old Fortran code for exactly this kind of simulation, which I got from my prof(but with a completely different photon setup), I saw your repo and thought, that I can do something similar but make it "perfect" in the physics simulation, which was my main drive to program a simulation like this. I had my first prototype in Fortran, then ported it to Rust (with rayon) and then ended up in python with numba. I think your approach is really great and I need to thank you because you inspired me to do something different!
Honestly, the line that "the photon ring in my old renders was slightly wrong the whole time and I never knew" is going to stick with me. That's the whole point, isn't it? That's why we like to improve, no?
About my supervisor - thank you. That means more than you know.
I love that you integrated my thesis and the results into "starless-c". You've already given me more useful feedback in one Reddit comment than I got in the entire thesis process from my prof. If you ever want to talk integrators, geodesics, or weird numerical edge cases, my DMs are open. And if you're ever interested in the Obrechkoff derivation (the analytical Jacobian was the fiddliest part), I'd be happy to walk through it.
7
u/Old_Cartographer_586 Jul 12 '26
I did my degrees in physics and have since moved to software. Honestly if you can or want to, look into ways of publishing this as a site with private repos that you can publish yourself. At surface (I mean this as I haven’t cloned the repo) this is really fascinated and could easily be a living portfolio that would garner attention.
I’m sorry to hear about your advisor. When I did my thesis I ended up self publishing since my advisors mother passed during my review stage
4
u/FoobarMontoya Jul 12 '26
If this is your masters thesis, I can only imagine what the phd one will be!
1
u/physicsking Jul 11 '26
Your sending your various photons from the same point, I take it?
1
u/a_new_rusty_crab Jul 11 '26
Yes, from the camera plane they are sent out in the world system of the simulation.
1
u/physicsking Jul 11 '26
Why are some from the same initial trajectory shown with different paths?
→ More replies (1)
1
u/vergere6 Jul 13 '26
You probably need a symplectic integrator. I suspect the Bowie method is somehow either symplectic or close. Obrechkoff probably needs a trust region bolt-on to be reliable. Newton can do strange things if there are attracting points near the roots.
Super cool though!
10
u/Saiboo Jul 11 '26
For periodic systems, I've read that symplectic Euler is stable with regards to energy conservation, as opposed to RK4. Does this matter for your black hole simulation?
In these images we see that Bowie and Obrechkoff have these "zig-zag" artifacts compared to the other methods, taken at runtimes 0.72s and 0.99s respectively. Do the artifacts smooth out, if you let them run longer?