WHFast512 tutorial¶

This tutorial shows how to use WHFast512, the highly optimized version of the symplectic Wisdom-Holman integrator. WHFast512 is written in x86 assembly and requires a CPU with AVX512 support. Currently only Linux is supported.

WHFast512 is very fast, up to 8 times faster than WHFast. However, it comes with some restrictions. Most importantly, it only works if the number of particles is less than or equal to 9 (= 8 planets + 1 star). It does not support test particles or variational particles. And you must work in units where $G=1$. But since gravity is scale invariant, you can always rescale your system such that $G=1$ without changing physics (see Units.ipynb.

First WHFast512 integration

Setting up a WHFast512 simulation works the same way as setting up any other REBOUND simulation:

In [1]:

import rebound
import numpy as np
sim = rebound.Simulation()
sim.add("solar system")  # built-in initial conditions for testing
sim.integrator = "whfast512"

import rebound import numpy as np sim = rebound.Simulation() sim.add("solar system") # built-in initial conditions for testing sim.integrator = "whfast512"

Since we work in units in which $G=1$ and we assume a length unit in the code to be one astronomical unit, the time unit in the code is $yr/(2\pi)$. We can thus set a timestep of 5 days the following way:

In [2]:

sim.dt = 5.0/365.25*np.pi*2

sim.dt = 5.0/365.25*np.pi*2

We also need to set the exact_finish_time flag because WHFast512 does not adjust the timestep we set:

In [3]:

sim.exact_finish_time = False

sim.exact_finish_time = False

Let's integrate this forward in time one step:

In [4]:

sim.steps(1)

sim.steps(1)

If we print out the time, you will notice that it is not 5 days but rather 13000 years:

In [5]:

sim.t/(np.pi*2.0)

sim.t/(np.pi*2.0)

Out[5]:

13689.253935660507

This is because by default WHFast512 concatenates $10^6$ timesteps to achieve high performance:

In [6]:

sim.integrator.concatenate_steps

sim.integrator.concatenate_steps

Out[6]:

1000000

You can change that number to a smaller or larger value depending on how fine-grained you required your output to be.

Speed

We next compare the speed of WHFast512 compared to WHFast. To do that let's copy the WHFast512 simulation, then integrate both foward in time for 200k years.

In [7]:

sim_whfast = sim.copy()
sim.integrate(2e5*np.pi*2, exact_finish_time=False)
sim_whfast.integrator = "whfast"
sim_whfast.integrate(sim.t, exact_finish_time=False) # integrating to same time as WHFast512

sim_whfast = sim.copy() sim.integrate(2e5*np.pi*2, exact_finish_time=False) sim_whfast.integrator = "whfast" sim_whfast.integrate(sim.t, exact_finish_time=False) # integrating to same time as WHFast512

For this test case, we get a speed up of:

In [8]:

print("%.1fx"%(sim_whfast.walltime/sim.walltime))

print("%.1fx"%(sim_whfast.walltime/sim.walltime))

8.2x

Note that a longer run would provide a more accurate runtime comparison.

Advanced options

WHFast512 is using Jacobi coordinates for the Hamiltonian Splitting and support symplectic correctors of order 17. In addition, WHFast512 supports general relativistic corrections. You can turn these features on with the following flags:

In [9]:

sim.integrator.gr_potential = True
sim.integrator.corrector = 17

sim.integrator.gr_potential = True sim.integrator.corrector = 17

Integrating multiple systems simultaneously

If your planetary systems contains 2, 3, or 4 planets, then you can combine the integration of multiple systems as WHFast512 is able to integrate the orbits of 8 planets in parallel. Here is an example that integrates four slightly different 2-planet systems in parallel with WHFast512 by setting the N_systems variable:

In [10]:

sim_combined = rebound.Simulation()
for i in range(4):
    sim = rebound.Simulation()
    sim.add(m=1) # star
    sim.add(m=1e-3, a=1)
    sim.add(m=1e-3, a=2+i/2.0, e=0.1)
    # adding all particles to the combined simulation
    for p in sim.particles:
        sim_combined.add(p)
# integrating the combined simulation
sim_combined.exact_finish_time = False
sim_combined.dt = 5.0/365.25*np.pi*2
sim_combined.integrator = "whfast512"
sim_combined.integrator.N_systems = 4
sim_combined.steps(1)

sim_combined = rebound.Simulation() for i in range(4): sim = rebound.Simulation() sim.add(m=1) # star sim.add(m=1e-3, a=1) sim.add(m=1e-3, a=2+i/2.0, e=0.1) # adding all particles to the combined simulation for p in sim.particles: sim_combined.add(p) # integrating the combined simulation sim_combined.exact_finish_time = False sim_combined.dt = 5.0/365.25*np.pi*2 sim_combined.integrator = "whfast512" sim_combined.integrator.N_systems = 4 sim_combined.steps(1)

We can then untangle the final result:

In [11]:

for i in range(4):
    sim = rebound.Simulation()
    for j in range(i*3, (i+1)*3):
        sim.add(sim_combined.particles[j])
    print("Semi-major axis of planet 2 in simulation %d: %f", i, sim.particles[2].a)

for i in range(4): sim = rebound.Simulation() for j in range(i*3, (i+1)*3): sim.add(sim_combined.particles[j]) print("Semi-major axis of planet 2 in simulation %d: %f", i, sim.particles[2].a)

Semi-major axis of planet 2 in simulation %d: %f 0 2.0017178895698557
Semi-major axis of planet 2 in simulation %d: %f 1 2.5002496749523435
Semi-major axis of planet 2 in simulation %d: %f 2 3.000244732224506
Semi-major axis of planet 2 in simulation %d: %f 3 3.5001849012088213