Using a C Heartbeat function¶

The heartbeat function of a REBOUND simulation gets called after every timestep. There are many different things you can do in a heartbeat function, for example creating outputs, adding particles, adjusting parameters that depend on time, etc. REBOUND supports heartbeat functions in both its C and python interface.

A python heartbeat function can sometimes become the bottleneck of a simulation because it gets called every single timestep. This tutorial shows you how to implement the heartbeat function in C, then link it to REBOUND using python. Note that alternatively you can of course always just use the C version of REBOUND directly and never bother with python at all.

We start by creating a REBOUND simulation which contains the planets of our Solar System as a test case.

In [1]:

import rebound
sim = rebound.Simulation()
rebound.data.add_solar_system(sim)
sim.integrator = "whfast"
sim.dt = sim.particles[1].P/30.13 # About 30 steps for each Mercury Orbit

Let us first create a simple heartbeat function in python. It simply calculates the eccentricity of Mercury (you could do something with it, here we just calculate it and then ignore it).

In [3]:

# We use a global variable to store the value of the eccentricity
e = 0 
def heartbeat(sim_pointer):
    global e
    # The function argument is a pointer to the simulation:
    # Here we get its contents:
    sim = sim_pointer.contents 
    e = sim.particles[1].e
sim.heartbeat = heartbeat
sim.integrate(sim.t+1)
print("Eccentricity: %f " %e)

Eccentricity: 0.205636

Let's measure how long it takes to integrate 1000 orbits:

In [4]:

import time
start = time.time()
sim.integrate(sim.t + sim.particles[1].P*1000)
stop = time.time()
print("Runtime: %f s"%(stop-start))

Runtime: 0.474900 s

We now implement this in C. For this to work, it is best to download and work with a full REBOUND repository (download the package from github, rather than just installing the python package with pip install).

We first write our heartbeat function in C. The following cell writes to a new file in the current directory, heartbeat.c (you can also use an external editor and terminal window to do the same without the jupyter magic commands):

In [6]:

%%writefile heartbeat.c
#include "rebound.h"
double e =0; // global variable
void heartbeat(struct reb_simulation* sim_pointer){
    struct reb_orbit orbit = reb_orbit_from_particle(sim_pointer->G, sim_pointer->particles[0], sim_pointer->particles[1]);
    e = orbit.e;
}

Overwriting heartbeat.c

Before we compile and link our heartbeat function as a shared library, we need the REBOUND header file and the shared library file. Different operating systems and compilers handles the paths to shared libraries differently. This can quickly get rather frustrating. If you're familiar with C, by all means go ahead and do it the proper way. A hack to get around most of these difficulties is to simply copy the REBOUND header and library to the current folder.

In [7]:

!cp ../src/librebound.so .
!cp ../src/rebound.h .

If you installed REBOUND with pip, you can look up the paths to the two files in python and e.g. create symlinks to them into your working directory.

In [8]:

from pathlib import Path
print(rebound.__libpath__)
print(Path(rebound.__file__).parent / "rebound.h")

/path/to/your/venv/lib/python3.9/site-packages/rebound/../librebound.cpython-39-x86_64-linux-gnu.so
/path/to/your/venv/lib/python3.9/site-packages/rebound/rebound.h

Next we can compile and link the code. -fPIC instructs gcc to create Position Independent Code, which might not be needed on your operating system.

In [9]:

!gcc -c -O3 -fPIC heartbeat.c -o heartbeat.o

In [10]:

!gcc -L. -shared heartbeat.o -o heartbeat.so -lrebound

Now we load the library using ctypes.

In [11]:

from ctypes import cdll
clibheartbeat = cdll.LoadLibrary("heartbeat.so")

We can now finally set the function pointer in our simulation to the new heartbeat function and then run the simulation.

In [12]:

sim.heartbeat = clibheartbeat.heartbeat

In [13]:

start = time.time()
sim.integrate(sim.t + sim.particles[1].P*1000)
stop = time.time()
print("Runtime: %f s"%(stop-start))

Runtime: 0.153589 s

Note that the simulation runs significantly faster using the C heartbeat function as we avoid all the python overhead.

We can print out the value of the global variable e in the heartbeat library (using global variables in a shared library is not the best way to store data - all simulations will see the same variable and you could end up with unexpected behaviour if you are running multiple simulations in one python program).

In [14]:

from ctypes import c_double
print(c_double.in_dll(clibheartbeat,"e").value)

0.2057293010348337

In [ ]: