Remember when we implemented the Barnes-Hut method in Python, but it was too slow to be of any practical use?

Well, I pushed that project much farther than it was ever meant to go, and ended up making it run 100 times faster, with optional OpenMP parallelism to boot, and it’s still written in 100% Python syntax. The trick is JIT-compiling it all with Numba. I give you pykdgrav.

Benchmarking against GIZMO on an identical problem, I generally find that GIZMO is about 4 times faster in computing the forces. Mind you, GIZMO is state-of-the-art, written in pure C, and has benefited from 20 years of GADGET development, so I was actually expecting much worse.

So why did I code this when there are so many other highly-optimized N-body codes out there? Because it’s fun, and because I wanted to show that C-like performance in Python numerical computing is becoming a reality, even for non-trivial, non-vectorizable algorithms. This project would not have been possible even 6 months ago, since it relies on Numba’s jitclass.

Like any Python package, it is also 100% portable and trivial to integrate into your own code. It is actually efficient enough that I have started using it in my day-to-day simulation analysis of star clusters.