Live data from Hacker News

Writing N-body gravity simulations code in Python

alvinng4.github.io

31–40 of 51 posts

Re: Writing N-body gravity simulations code in Python

#31
post #23
post #12

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.

Oh wow! I implemented BH in rust not long ago. Was straightforward. Set it up with grpahics so I could see the cubes etc. Then looked into FMM... I couldn't figure out where to start! Looked very formidable. multipole seems to be coming up in everything scientific I look at these days...

Re: Writing N-body gravity simulations code in Python

#32
post #13

Earlier 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‽

You need cards like the A100 / H100 / H200 / B200. GPUs aren’t fundamentally worse at f64. Nvidia just makes it worse in many cards for market segmentation.

Re: Writing N-body gravity simulations code in Python

#34

Earlier 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?

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") 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

#35
post #10

My 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.

There are countless applications other than the n-body problem for which my point holds

Re: Writing N-body gravity simulations code in Python

#36
post #12

Next 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…

Supercomputers will simulate trillions of masses. The HACC code, commonly used to verify the performance of these machines, uses a uniform grid (interpolation and a 3D FFT) and local corrections to compute the motion of ~8 trillion bodies.

Re: Writing N-body gravity simulations code in Python

#37
post #10

My 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.

Cut the time interval shorter and shorter, and then BAM you've independently discovered calculus.

Re: Writing N-body gravity simulations code in Python

#38

For comparison, the famous Debian language benchmark competition: https://benchmarksgame-team.pages.debian.net/benchmarksgame/...

And for those who have access

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

#39

Earlier 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…

>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.

IIRC that is what I did in the end. It is fudge, but it works.

Re: Writing N-body gravity simulations code in Python

#40
post #7

Once 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).

Oh for sure there is a ton that is going on behind the scenes when you substitute `np` --> `jnp`. But as someone who worked on gravitational dynamics a decade ago, it's really incredible that so much work has been done to make running numerical calculations on a GPU so straightforward for the programmer.

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.

Post reply on HN