Earlier quoted context omitted.
That's mentioned on the "Conclusions" page of TFA: > Large-scale simulation: So far, we have only focused on systems with a few objects. What about large-scale systems with thousands or millions of objects? Turns out it is not so easy because the computation of gravity scales as . Have a look at Barnes-Hut algorithm to see how to speed up the simulation. In fact, we have documentations about it on this website as wel…
C or C++? Ha! I've implemented FMM in Fortran 77. (Not my choice; it was a summer internship and the "boss" wanted it that way.) It was a little painful.
Writing N-body gravity simulations code in Python
31–40 of 51 posts
Re: Writing N-body gravity simulations code in Python
#32Earlier quoted context omitted.
Not true for regular GPUs like RTX5090. They have atrocious float64 performance compared to CPU. You need a special GPU designed for scientific computations (many float64 cores)
This is confusing: GPUs seem great for scientific computations. You often want f64 for scientific computation. GPUs aren't good with f64. I'm trying to evaluate if I can get away with f32 for GPU use for my molecular docking software. Might be OK, but I've hit cases broadly where f64 is fine, but f32 is not. I suppose this is because the dominant uses games and AI/ML use f32, or for AI even less‽
Re: Writing N-body gravity simulations code in Python
#33https://benchmarksgame-team.pages.debian.net/benchmarksgame/...
Re: Writing N-body gravity simulations code in Python
#34Earlier quoted context omitted.
One way of doing that is by each time step delta perform two integrations, one which give you x_1, and another x_2, with different accuracies. The error is estimated by the difference |x_1 - x_2| and you make the error match your tolerance by adjusting the time step. Naturally, this difference becomes big when two objects are close together, since the acceleration will induce a large change of velocity, and a lower t…
I have run into the problem where a constant time step can suddenly result in bodies getting flung out of the simulation because they go very close. Your solution sounds interesting, but isn't it only practical when you have a small number of bodies?
If you see bodies flung out after close passes, three solutions are available: reduce the time step, use a higher order time integrator, and (the most common method) add regularization. Regularization (often called "softening") removes the singularity by adding a constant to the squared distance. So 1 over zero becomes one over a small-ish and finite number.
Re: Writing N-body gravity simulations code in Python
#35My favorite thing about this kind of code is that people are constantly inventing new techniques to do time integration. It's the sort of thing you'd think was a solved problem but then when you read about time integration strategies you realize how rich the space of solutions are. And they are more like engineering problems than pure math, with tradeoffs and better fits based on the kind of system.
They're more like engineering than pure math because pure math hasn't solved the n-body problem.
Re: Writing N-body gravity simulations code in Python
#36Next logical optimization: Barnes Hut? Groups source bodies using a recursive tree of cubes. Gives huge speedups with high body counts. FMM is a step after, which also groups target bodies. Much more complicated to implement.
That's mentioned on the "Conclusions" page of TFA: > Large-scale simulation: So far, we have only focused on systems with a few objects. What about large-scale systems with thousands or millions of objects? Turns out it is not so easy because the computation of gravity scales as . Have a look at Barnes-Hut algorithm to see how to speed up the simulation. In fact, we have documentations about it on this website as wel…
Re: Writing N-body gravity simulations code in Python
#37My favorite thing about this kind of code is that people are constantly inventing new techniques to do time integration. It's the sort of thing you'd think was a solved problem but then when you read about time integration strategies you realize how rich the space of solutions are. And they are more like engineering problems than pure math, with tradeoffs and better fits based on the kind of system.
Decades ago ... I think it was Computer Recreations column in "Scientific American" ... the strategy, since computers were less-abled then, was to advance all the bodies by some amount of time — call it a "tick" — when bodies got close, the "tick" got smaller: therefore the calculations more nuanced, precise. Further apart and you could run the solar system on generalities.
Re: Writing N-body gravity simulations code in Python
#38For comparison, the famous Debian language benchmark competition: https://benchmarksgame-team.pages.debian.net/benchmarksgame/...
2022 "N-Body Performance With a kD-Tree: Comparing Rust to Other Languages"
https://ieeexplore.ieee.org/document/10216574
2025 "Parallel N-Body Performance Comparison: Julia, Rust, and More"
https://link.springer.com/chapter/10.1007/978-3-031-85638-9_...
Re: Writing N-body gravity simulations code in Python
#39Earlier quoted context omitted.
I have run into the problem where a constant time step can suddenly result in bodies getting flung out of the simulation because they go very close. Your solution sounds interesting, but isn't it only practical when you have a small number of bodies?
Yes, the author uses a globally-adaptive time stepper, which is only efficient for very small N. There are adaptive time step methods that are local, and those are used for large systems. If you see bodies flung out after close passes, three solutions are available: reduce the time step, use a higher order time integrator, and (the most common method) add regularization. Regularization (often called "softening") remo…
IIRC that is what I did in the end. It is fudge, but it works.
Re: Writing N-body gravity simulations code in Python
#40Once you have the matrix implementation in Step 2 (Implementation 3) it's rather straightforward to extend your N-body simulator to run on a GPU with Jax --- you can just add `import jax.numpy as jnp` and replace all the `np.`s with `jnp`s. For a few-body system (e.g., the Solar System) this probably won't provide any speedup. But once you get to ~100 bodies you should start to see substantial speedups by running the…
> rather straightforward For the programmer, yes it is easy enough. But there is a lot of complexity hidden behind that change. If you care about how your tools work that might be a problem (I'm not judging).
Obviously if you want maximum performance you'll probably still have to roll up your sleeves and write CUDA yourself, but there are a lot of situations where you can still get most of the benefit of using a GPU and never have to worry yourself over how to write CUDA.