# \[ANN\] AstroNbodySim.jl: Unitful and differentiable gravitational N-body simulation code in Julia

**URL:** <https://discourse.julialang.org/t/ann-astronbodysim-jl-unitful-and-differentiable-gravitational-n-body-simulation-code-in-julia/74782>\
**Category:** Package Announcements\
**Tags:** package\
**Created:** [January 18, 2022, 6:10am UTC](https://discourse.julialang.org/t/ann-astronbodysim-jl-unitful-and-differentiable-gravitational-n-body-simulation-code-in-julia/74782 "2022-01-18T06:10:50Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![islent](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/islent/32/11326_2.png) [@islent](https://discourse.julialang.org/u/islent)\
**Post date:** [January 18, 2022, 6:10am UTC](https://discourse.julialang.org/t/ann-astronbodysim-jl-unitful-and-differentiable-gravitational-n-body-simulation-code-in-julia/74782/1 "2022-01-18T06:10:50Z")

</div>

> **[GitHub - JuliaAstroSim/AstroNbodySim.jl: Unitful and differentiable...](https://github.com/JuliaAstroSim/AstroNbodySim.jl)**
>
> Unitful and differentiable gravitational N-body simulation code in Julia - GitHub - JuliaAstroSim/AstroNbodySim.jl: Unitful and differentiable gravitational N-body simulation code in Julia

It is a package for scientific gravitational n-body simulations, with following features:

- Compute with units
- User-friendly
  - Well documented
  - Readable programming
  - Vectorized array operations
  - Dispatch on types for various simulation settings
  - Float16, Float32, Float64, Int128, BigFloat, Measurement, etc.

- Cross-platform: Linux, Windows, MacOS. Easy to deploy
- Hybrid Parallelism: multi-threading, distributed parallelism, GPU acceleration
- Modularity and Versatility: 9 packages, designed for general purposes, highly extentable
- Realtime visualzation (interactive)
- Auto-test workflow

Documentation: [Home · AstroNbodySim.jl](https://juliaastrosim.github.io/AstroNbodySim.jl/dev)

Here’s some figures and animations to demonstrate what `AstroNbodySim.jl` can do:

1. Realtime visualization of simulations on GPU  
 ![Plummer](https://global.discourse-cdn.com/julialang/original/3X/c/3/c3673364572e64d47fc587b4f8e71590857200a6.gif)

2. Galactic collision  
 ![GalacticCollision](https://global.discourse-cdn.com/julialang/original/3X/8/7/878f7e9b5f3c1f96786a22a5f5689a0f0bc8081a.gif)

3. Uncertainty propagation  

4. Autodiff of background potential field  

5. User-difined pipeline: Tidal disruption event (TDE)

6. Lagrange radii and scale radius

7. Solar System  
 ![SolarSystem](https://global.discourse-cdn.com/julialang/original/3X/2/8/286c7480a6fc6a837f3adeaa96ebfdc3d8118189.gif)

Package ecosystem:

- Basic data structure: [PhysicalParticles.jl](https://github.com/JuliaAstroSim/PhysicalParticles.jl)
- File I/O: [AstroIO.jl](https://github.com/JuliaAstroSim/AstroIO.jl)
- Initial Condition: [AstroIC.jl](https://github.com/JuliaAstroSim/AstroIC.jl)
- Parallelism: [ParallelOperations.jl](https://github.com/JuliaAstroSim/ParallelOperations.jl)
- Trees: [PhysicalTrees.jl](https://github.com/JuliaAstroSim/PhysicalTrees.jl)
- Meshes: [PhysicalMeshes.jl](https://github.com/JuliaAstroSim/PhysicalMeshes.jl)
- Plotting: [AstroPlot.jl](https://github.com/JuliaAstroSim/AstroPlot.jl)
- Simulation: [AstroNbodySim](https://github.com/JuliaAstroSim/AstroNbodySim.jl)
- Benchmark: [BenchmarkPlots](https://github.com/JuliaAstroSim/BenchmarkPlots.jl)
- Parameter space exploration: [ParameterSpace](https://github.com/JuliaAstroSim/ParameterSpace.jl)

Issues and PRs are welcomed!

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [January 18, 2022, 6:21am UTC](https://discourse.julialang.org/t/ann-astronbodysim-jl-unitful-and-differentiable-gravitational-n-body-simulation-code-in-julia/74782/2 "2022-01-18T06:21:11Z")

</div>

Just wondering, is there a reason this isn’t using `DifferentialEquations` for the actual simulation?

---

<div class="post-metadata">

**Author:** ![islent](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/islent/32/11326_2.png) [@islent](https://discourse.julialang.org/u/islent)\
**Post date:** [January 18, 2022, 6:41am UTC](https://discourse.julialang.org/t/ann-astronbodysim-jl-unitful-and-differentiable-gravitational-n-body-simulation-code-in-julia/74782/3 "2022-01-18T06:41:50Z")

</div>

1. We aim for heavy-load Nbody problems (\> 10^6 particles) on HPC, so that the parallelization algorithm has to be customized
2. Same as above, we have many specific simulation steps for scientific goals
3. About force solvers:
  - `DifferentialEquations` does not have a distributed octree force solver
  - Direct summation method does not need `DifferentialEquations`
  - `DifferentialEquations` might be integrated as a sub-solver for particle-mesh solvers in the future

4. The code architecture is evolved from [Gadget-2](https://wwwmpa.mpa-garching.mpg.de/gadget/). (We have many functions designed for Gadget-2 users)

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [January 18, 2022, 11:05am UTC](https://discourse.julialang.org/t/ann-astronbodysim-jl-unitful-and-differentiable-gravitational-n-body-simulation-code-in-julia/74782/4 "2022-01-18T11:05:42Z")

</div>

That doesn’t make much sense though. There’s nothing specific to time steppers that would be customized for parallelization outside of the array interface. It seems like there’s an abstraction break here where time steppers, parallel/unitful arrays, and fast multipole method summations are put together where in practice the parallel efficiency of the three do not necessarily require interwoven implementation in order to achieve maximal efficiency. Pulling the three apart would ultimately be more efficient due ot how it can interact with other parts of the ecosystem.

---

<div class="post-metadata">

**Author:** ![xzackli](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/xzackli/32/38301_2.png) [@xzackli](https://discourse.julialang.org/u/xzackli)\
**Post date:** [January 18, 2022, 2:40pm UTC](https://discourse.julialang.org/t/ann-astronbodysim-jl-unitful-and-differentiable-gravitational-n-body-simulation-code-in-julia/74782/5 "2022-01-18T14:40:01Z")

</div>

This looks great! I love the diverse range of supporting packages, but it seems like your docs are light on the scientific details (presumably you’re putting those in a paper).

1. Do you have tests of i.e. cluster evolution with respect to standard codes like gadget?
2. It’s not clear to me what your integrator is. What’s the scheme for time stepping?
3. What’s the tree scheme? BH?

---

<div class="post-metadata">

**Author:** ![carstenbauer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/carstenbauer/32/4981_2.png) [@carstenbauer](https://discourse.julialang.org/u/carstenbauer)\
**Post date:** [January 18, 2022, 3:00pm UTC](https://discourse.julialang.org/t/ann-astronbodysim-jl-unitful-and-differentiable-gravitational-n-body-simulation-code-in-julia/74782/6 "2022-01-18T15:00:56Z")

</div>

> [@islent](#):
>
> We aim for heavy-load Nbody problems (\> 10^6 10610^6 particles) on HPC, so that the parallelization algorithm has to be customized

From taking a quick look, It seems like you are using Julia’s built-in distributed computing facilities. I am curious

1. why you chose it over MPI,
2. how well it works for you, and
3. whether you plan to compare to MPI at some point.

---

<div class="post-metadata">

**Author:** ![carstenbauer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/carstenbauer/32/4981_2.png) [@carstenbauer](https://discourse.julialang.org/u/carstenbauer)\
**Post date:** [January 18, 2022, 3:27pm UTC](https://discourse.julialang.org/t/ann-astronbodysim-jl-unitful-and-differentiable-gravitational-n-body-simulation-code-in-julia/74782/7 "2022-01-18T15:27:04Z")

</div>

> [@ChrisRackauckas](#):
>
> Pulling the three apart would ultimately be more efficient due ot how it can interact with other parts of the ecosystem.

While I agree, improving generic array abstractions seems like a very different and potentially much harder development effort to me.

Also, what are suggesting specifically? Should one use `DArray`s in combination with DiffEq for distributed ODE solving? If I just read the [DiffEq docs](https://diffeq.sciml.ai/dev/basics/faq/#GPUs,-multithreading-and-distributed-computation-support), I would think that I shouldn’t go down this route.

> GPUArrays.jl (CuArrays.jl), ArrayFire.jl, DistributedArrays.jl have been tested and work in various forms, where **the last one is still not recommended for common use yet**.

So, I guess I should develop my own custom array type then? Is this mentioned / explained somewhere in the docs? I think an example showcasing something like this would be great. (Maybe I’m just not aware of it?)

---

<div class="post-metadata">

**Author:** ![islent](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/islent/32/11326_2.png) [@islent](https://discourse.julialang.org/u/islent)\
**Post date:** [January 18, 2022, 4:03pm UTC](https://discourse.julialang.org/t/ann-astronbodysim-jl-unitful-and-differentiable-gravitational-n-body-simulation-code-in-julia/74782/8 "2022-01-18T16:03:33Z")

</div>

The code project started from Julia 0.7, in 2018. At that time, `DifferentialEquations` was not our preferred choice. However, I would add DiffEq as one of the solvers in future development.

---

<div class="post-metadata">

**Author:** ![islent](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/islent/32/11326_2.png) [@islent](https://discourse.julialang.org/u/islent)\
**Post date:** [January 18, 2022, 4:15pm UTC](https://discourse.julialang.org/t/ann-astronbodysim-jl-unitful-and-differentiable-gravitational-n-body-simulation-code-in-julia/74782/9 "2022-01-18T16:15:04Z")

</div>

Yes, I’m preparing a paper, which will be submitted in one or two weeks. All figures are generated by codes in `examples/`

1. In `test/`, we have auto-tests validating the algorithms. For example, the evolution of energy and momentum are compared with Gadget-2. (Gadget-2 need unix-like platform, thus the comparison is not operated by default)
2. We have an explicit Euler integrator and an KDK Leapfrog integrator, both supporting fixed time-step and adaptive time-step (enough for collisionless n-body problems).
3. Basically same with Gadget-2: Hilbert-Peano octree.

---

<div class="post-metadata">

**Author:** ![islent](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/islent/32/11326_2.png) [@islent](https://discourse.julialang.org/u/islent)\
**Post date:** [January 18, 2022, 4:33pm UTC](https://discourse.julialang.org/t/ann-astronbodysim-jl-unitful-and-differentiable-gravitational-n-body-simulation-code-in-julia/74782/10 "2022-01-18T16:33:26Z")

</div>

1. Easier to understand and debug
2. `ParallelOperations` is type-stable and fast. Data communication is realized in minimum cost. We have benchmarks in `Benchmark/` (tens of seconds for one step of 10^6 particles). Optimization works on distributed-memories are in progress. The multi-threading method is still a problem, because `LoopVecterization` do not support user-defined structs.
3. Generally speaking, direct summation method (CPU) and tree method (CPU only) in `AstroNbodySim` is slower than Gadget-2. One reason is that we rebuild the tree in each step. There are lots of optimization works to do. Fortunately, the simple GPU implementation of direct summation method is comparably fast (when using Float32).

---

<div class="post-metadata">

**Author:** ![islent](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/islent/32/11326_2.png) [@islent](https://discourse.julialang.org/u/islent)\
**Post date:** [January 18, 2022, 4:43pm UTC](https://discourse.julialang.org/t/ann-astronbodysim-jl-unitful-and-differentiable-gravitational-n-body-simulation-code-in-julia/74782/11 "2022-01-18T16:43:05Z")

</div>

Our implementation of distributed simulation data is quite similar to `DistributedArrays`. Particle data are handled by `StructArrays` for efficient memory accessing.

I will mention this in docs.

---

<div class="post-metadata">

**Author:** ![xzackli](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/xzackli/32/38301_2.png) [@xzackli](https://discourse.julialang.org/u/xzackli)\
**Post date:** [January 18, 2022, 4:47pm UTC](https://discourse.julialang.org/t/ann-astronbodysim-jl-unitful-and-differentiable-gravitational-n-body-simulation-code-in-julia/74782/12 "2022-01-18T16:47:28Z")

</div>

Maybe I should just wait for the paper, but I’m curious if uncertainty propagation actually works? Maybe it’s just for small systems. Pure speculation: at kpc scales and long timescales, I would expect a particle to be destined to enter a particular halo, but once it enters the halo the N-body calculation becomes Monte Carlo. So the uncertainty on position should just be the radius of the halo?

---

<div class="post-metadata">

**Author:** ![islent](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/islent/32/11326_2.png) [@islent](https://discourse.julialang.org/u/islent)\
**Post date:** [January 18, 2022, 5:02pm UTC](https://discourse.julialang.org/t/ann-astronbodysim-jl-unitful-and-differentiable-gravitational-n-body-simulation-code-in-julia/74782/13 "2022-01-18T17:02:41Z")

</div>

uncertainty propagation is supported by Measurements.jl, which tracks all uncertainties during computation. I don’t fully understand your question, but if you have all uncertainty measurements in your initial conditions, you can get a reliable uncertainty prrdiction (assuming time integration is accurate), and differentiate the result over initial conditions.

See examples/01-binary.jl

---

<div class="post-metadata">

**Author:** ![KZiemian](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kziemian/32/9020_2.png) [@KZiemian](https://discourse.julialang.org/u/KZiemian)\
**Post date:** [January 18, 2022, 5:06pm UTC](https://discourse.julialang.org/t/ann-astronbodysim-jl-unitful-and-differentiable-gravitational-n-body-simulation-code-in-julia/74782/14 "2022-01-18T17:06:30Z")

</div>

It looks amazing. Where I can find information about your planes for the future?

---

<div class="post-metadata">

**Author:** ![lmiq](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lmiq/32/18314_2.png) [@lmiq](https://discourse.julialang.org/u/lmiq)\
**Post date:** [January 18, 2022, 5:23pm UTC](https://discourse.julialang.org/t/ann-astronbodysim-jl-unitful-and-differentiable-gravitational-n-body-simulation-code-in-julia/74782/15 "2022-01-18T17:23:12Z")

</div>

> [@islent](#):
>
> but if you have all uncertainty measurements in your initial conditions, you can get a reliable uncertainty prrdiction

See this comment and some associated discussion: [Notebook on particle simulations - #3 by ChrisRackauckas](https://discourse.julialang.org/t/notebook-on-particle-simulations/68496/3)

A curiosity: in a typical simulation with periodic boundary conditions, how big (system size) and how many particles (10^6?) are we talking about?

Can these calculations be performed with a cutoff (great enough such that the gravitational potential is small yet much smaller than the system size?). How many time-steps/second is a good benchmark? (I’m curious how would the trees compare with [cell lists](https://github.com/m3g/CellListMap.jl) with a sufficiently large cutoff if that is possible).

---

<div class="post-metadata">

**Author:** ![giordano](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/giordano/32/2166_2.png) [@giordano](https://discourse.julialang.org/u/giordano)\
**Post date:** [January 18, 2022, 5:27pm UTC](https://discourse.julialang.org/t/ann-astronbodysim-jl-unitful-and-differentiable-gravitational-n-body-simulation-code-in-julia/74782/16 "2022-01-18T17:27:33Z")

</div>

Note that the `Measurements.jl` follows [linear error propagation theory](https://en.wikipedia.org/wiki/Propagation_of_uncertainty), which is based on some assumptions: errors are normally distributed, the first derivative is enough to locally describe the variation of the function (simple functions like `sin` and `cos` already fail this condition locally), etc… If some of these assumptions don’t hold, you get basically garbage numbers out, which I think is what @xzackli was suggesting above

---

<div class="post-metadata">

**Author:** ![yzhu](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/yzhu/32/32910_2.png) [@yzhu](https://discourse.julialang.org/u/yzhu)\
**Post date:** [January 18, 2022, 7:47pm UTC](https://discourse.julialang.org/t/ann-astronbodysim-jl-unitful-and-differentiable-gravitational-n-body-simulation-code-in-julia/74782/17 "2022-01-18T19:47:28Z")

</div>

Thanks! We will publish a paper describing the code and future works soon. The plans will be also updated on GitHub afterwards.

---

<div class="post-metadata">

**Author:** ![KZiemian](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kziemian/32/9020_2.png) [@KZiemian](https://discourse.julialang.org/u/KZiemian)\
**Post date:** [January 18, 2022, 9:04pm UTC](https://discourse.julialang.org/t/ann-astronbodysim-jl-unitful-and-differentiable-gravitational-n-body-simulation-code-in-julia/74782/18 "2022-01-18T21:04:58Z")

</div>

I wish you good luck. I will try to watch this project regular, to stay on board with developments.

---

<div class="post-metadata">

**Author:** ![islent](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/islent/32/11326_2.png) [@islent](https://discourse.julialang.org/u/islent)\
**Post date:** [January 19, 2022, 2:08am UTC](https://discourse.julialang.org/t/ann-astronbodysim-jl-unitful-and-differentiable-gravitational-n-body-simulation-code-in-julia/74782/19 "2022-01-19T02:08:33Z")

</div>

> Where I can find information about your planes for the future?

Thanks! Check them here:

- [dev · GitHub](https://github.com/orgs/JuliaAstroSim/projects/4)
- [Roadmap of AstroNbodySim.jl · Issue #5 · JuliaAstroSim/AstroNbodySim.jl · GitHub](https://github.com/JuliaAstroSim/AstroNbodySim.jl/issues/5)

---

<div class="post-metadata">

**Author:** ![islent](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/islent/32/11326_2.png) [@islent](https://discourse.julialang.org/u/islent)\
**Post date:** [January 19, 2022, 2:28am UTC](https://discourse.julialang.org/t/ann-astronbodysim-jl-unitful-and-differentiable-gravitational-n-body-simulation-code-in-julia/74782/20 "2022-01-19T02:28:54Z")

</div>

> A curiosity: in a typical simulation with periodic boundary conditions, how big (system size) and how many particles (10^6?) are we talking about?

Dynamic simulation under periodic boundary conditions is not supported by now. It will be supported in cosmological simulations. In the case of galactic dynamics, we only support vacuum boundary conditions.

> Can these calculations be performed with a cutoff (great enough such that the gravitational potential is small yet much smaller than the system size?)

Yes, the gravitational potential is softened, or in other words, we solve collisionless Boltzmann equations. Check sec 2.9 in [Galactic Dynamics](https://books.google.com.hk/books?hl=en&lr=&id=6mF4CKxlbLsC&oi=fnd&pg=PP1&dq=galactic+dynamics+binney&ots=YZGxGbGcij&sig=yTuuLidoanBDWXjY8xo3BIT74G0&redir_esc=y&hl=zh-CN&sourceid=cndr#v=onepage&q=galactic%20dynamics%20binney&f=false) for details.

> How many time-steps/second is a good benchmark?

See the sections about time integration in [Gadget-2 paper](https://academic.oup.com/mnras/article/364/4/1105/1042826)

- In the case of adaptive time-steps, \Delta t is determined by the accelerations in the last time-step
- For fixed time-steps, it is reasonable to choose the smallest time-step in adaptive time-steps.

[Next page](https://discourse.julialang.org/t/ann-astronbodysim-jl-unitful-and-differentiable-gravitational-n-body-simulation-code-in-julia/74782.md?page=2)