GR scaling with large number of test particles

[daelee6]

Hello!

I was running simulations with a small number of active particles and a large number of test particles and noticed that using "gr" scaled as N^2 rather than N_active * N.

Here is a small example:

import rebound
import reboundx
from reboundx import constants
import numpy as np
import matplotlib.pyplot as plt
import time

rng = np.random.default_rng()

# source particle, not added to sim
source = rebound.Particle(x=1., y=-8., z=-2., 
                          vx=-0.8, vy=-2.4, vz=-0.2)

def add_particles(sim):
    p = rebound.Particle()
    p.x = rng.uniform(-0.1,0.1)
    p.y = rng.uniform(-0.1,0.1)
    p.z = rng.uniform(-0.1,0.1)
    p.vx = rng.uniform(-0.01,0.01)
    p.vy = rng.uniform(-0.01,0.01)
    p.vz = rng.uniform(-0.01,0.01)
    
    p += source
    sim.add(p)

# choice of GR routine
forces = [None, "gr_potential", "gr"]

# number of test particles
nparticles = np.array([10,40,160,640,2560])

t = np.zeros((len(forces), len(nparticles)))

for i, force in enumerate(forces):
    for j, n in enumerate(nparticles):
        # create simulation with Sun and Earth
        sim = rebound.Simulation()
        sim.units = ('yr', 'AU', 'Msun')
        sim.add(m=1., hash="Sun")
        sim.add(m=0.000001, a=1., hash="Earth")

        # two active particles
        sim.move_to_com()
        sim.N_active = 2
        
        # add in GR
        rebx = reboundx.Extras(sim)
        if force is not None:
            gr = rebx.load_force(force)
            rebx.add_force(gr)
            gr.params["c"] = constants.C * np.sqrt(sim.G) # set speed of light
            sim.particles["Sun"].params["gr_source"] = 1
        
        # add in test particles
        for __ in range(n):
            add_particles(sim)
        
        # integrate and check runtime
        start_time = time.time()
        sim.integrate(1.)
        end_time = time.time()
        t[i,j] = end_time - start_time
        print(n, t[i,j])
        sim = None
        rebx = None
    
    # plot runtime
    label = force
    if label is None:
        label = "none"
    plt.scatter(nparticles, t[i], label=label)

plt.plot(nparticles, t[0,-1]*(nparticles/nparticles[-1]),'k--',label=r"$O(n)$")
plt.plot(nparticles, t[2,-1]*(nparticles/nparticles[-1])**2,'k-.',label=r"$O(n^2)$")
plt.legend()

plt.xscale('log')
plt.yscale('log')
plt.xlabel("Number of particles")
plt.ylabel("Runtime (s)")

plt.show()

The result I get is

It looks like this can be adjusted by changing lines 77, 99, and 153 in "gr.c" to use N_active:

65a66,67
>     const int N_active = sim->N_active;
>
77c79
<     for(int i=0; i<N; i++){
---
>     for(int i=0; i<N_active; i++){
99c101
<     reb_particles_transform_inertial_to_jacobi_posvelacc(ps, ps_j, ps, N, N);
---
>     reb_particles_transform_inertial_to_jacobi_posvelacc(ps, ps_j, ps, N, N_active);
153c155
<     reb_particles_transform_jacobi_to_inertial_acc(ps, ps_j, ps, N, N);
---
>     reb_particles_transform_jacobi_to_inertial_acc(ps, ps_j, ps, N, N_active);

I've implemented it as a separate force and it seems to give the correct scaling while still giving identical simulations to "gr".

Thank you!

[dtamayo]

This is wonderful, thank you so much! I hadn't thought about the fact that the coordinate transformations could cause problems for large numbers of particles.

Would you like to create a pull request so it shows your contribution to the codebase, or would you prefer for me to add this in. Happy to do either!

[daelee6]

If you could add it, that would be very helpful! I'm not exactly sure how to treat N_active. I currently have it set to
const int N_active = sim->N_active > 0 ? sim->N_active : N;
to cover cases when N_active = -1; ... could there be any cases when N_active = 0?