Impossible slow simulation switching from 10k to 100k particles using hardsphere. Is this normal?

[Francyrad]

Dear @hannorein

I did some experiments for my moon formation scenario and i found this pattern:

Basically when i run my scenario using 10k particles, it takes usually 0.2 seconds with my cpu power of 20000. If i run 100k particles for the same problem, it's gonna cost up to 1 minute to do 10 s of simulation.

My question is: is this normal with hardsphere? I would expect just a 10x of execution time (2s), my particles are not even touching in the beginning of my simulation, or my code as somekind of a bottleneck using mpi?

Still thank you for your help and best regards

#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>
#include <math.h>
#include <time.h>
#include <sys/time.h>
#include <sys/stat.h>
#include <string.h>
#include "rebound.h"
#include "tools.h"
#include "output.h"



void heartbeat(struct reb_simulation* const r);
void heartbeat(struct reb_simulation* const r){
    if (reb_output_check(r,10.0*r->dt)){
        reb_output_timing(r,10);
    }
}



const double G = 6.67430e-11;
const double M_earth = 5.972e24;
const double R_earth = 6371000.0;
const double rho_particle = 3300.0; // kg/m^3
const double r_inner = 6372000.0;
const double r_outer = 25484000.0;
const double boxsize =  3 * r_outer;
const int N_particles = 10000;

const double M_moon = 7.342e22; // kg (massa lunare)
const double total_mass = M_moon; // massa totale delle particelle = 1 massa lunare
double random_uniform(double min, double max) {
    return min + (max - min) * ((double) rand() / (double) RAND_MAX);
}
double coefficient_of_restitution_constant(const struct reb_simulation* const r, double v) {
    // v is the normal impact velocity.
    // Here, we just use a constant coefficient of restitution
    return 0.6;
}

int main(int argc, char* argv[]){
    struct reb_simulation* r = reb_create_simulation();
    // Setup constants
    r->integrator    = REB_INTEGRATOR_LEAPFROG;
    r->gravity    = REB_GRAVITY_TREE;
    r->boundary    = REB_BOUNDARY_OPEN;
    r->collision = REB_COLLISION_TREE;
    r->coefficient_of_restitution = coefficient_of_restitution_constant;
    r->collision_resolve = reb_collision_resolve_hardsphere;
    r->opening_angle2    = 1.5;        // This constant determines the accuracy of the tree code gravity estimate.
    r->G         = G;
    r->softening     = 0.02;        // Gravitational softening length
    r->dt         = 1;        // Timestep
    
    // Setup root boxes for gravity tree.
    // Here, we use 2x2=4 root boxes (each with length 'boxsize')
    // This allows you to use up to 8 MPI nodes.
    reb_configure_box(r, 3* boxsize, 4, 2, 1);
    
    // Initialize MPI
    // This can only be done after reb_configure_box.
    reb_mpi_init(r);
    
    // Setup particles only on master node
    
    
    struct reb_particle Earth = {0};
    Earth.m = M_earth;
    Earth.r = R_earth;
    
    
    double M_particle = total_mass / N_particles;
    double R_particle = pow((3.0 * M_particle) / (4.0 * M_PI * rho_particle), 1.0 / 3.0) / 1000.0; // Raggio in km
        
    if (r->mpi_id==0){
        reb_add(r, Earth);
        for (int i=0;i<N_particles;i++){
            double radius = random_uniform(r_inner, r_outer);
            double theta = random_uniform(0, 2 * M_PI);

            struct reb_particle p = {0};
            p.m = M_particle;
            p.r = R_particle * 1000.0; // Raggio in metri
            p.x = radius * cos(theta);
            p.y = radius * sin(theta);
            p.z = 0;
            double v_kep = sqrt(G * M_earth / radius);
            p.vx = -v_kep * sin(theta);
            p.vy = v_kep * cos(theta);
            p.vz = 0;

            reb_add(r, p);
                        }
                    }

    r->heartbeat = heartbeat;
    reb_simulationarchive_automate_interval(r, "archive.bin", 100.);
    reb_integrate(r, INFINITY);

    reb_mpi_finalize(r);
    reb_free_simulation(r);
        
  
        
    }

[hannorein]

The tree code scaling is not $O(N)$ but $O(N log(N))$. Regardless, the issue with your setup is that you initialize your particles in 2D but scale their radii as if the particles were 3D. As a result you will have significantly more collisions between particles as you increase the number of particles.