REBOUND Simulation and PC Benchmarking Tests

The journey of a thousand lines of code, begins w/ a single CMake.txt file... that doesn't want to work.

—Me, attempting to compile anything from source

Abstract:

Project Github: https://github.com/0xVRGL/GalaxySim

REBOUND Website: https://rebound.hanno-rein.de/

REBOUND is an N-Body integrator used to mathematically and computationally model the movement of particles under the influence of gravity.

As all objects with mass will exert a gravitational pull, an absolute model of a solar system must account for all objects of mass regardless if their gravitational pull relatively small. For example, it is usually sufficient to only calculate the gravitational influence of the Sun when constructing a model of the solar system, but as each planet itself has mass, the planets too exert a gravitational force on the sun, thus pulling the sun off its own axis ever so slightly.

As a larger number of particles are added to the simulation, the number of gravitational interactions -and requisite equations needed- grow O(N²), resulting in ever inevitably computational times.

Henceforth the idea to use such as a PC benchmarking test, albeit a rudimentary one.

Rebound Simulation Code:

REBOUND's base code is written in C, whereas this project is written in C++ and compiled through Microsoft's vscode compiler.

As a result restrict header tags must be redefined to __restrict -as for rebounds library be understood by the compiler-, and _USE_MATH_DEFINES, must be defined before the C's main library is allowed to run.

#define _USE_MATH_DEFINES

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>

#define restrict __restrict
extern "C" {
#include "rebound.h"
}
#undef restrict

class TimeLogger {
    public:
        time_t simstart;
        time_t simend;
        time_t simdiff;

        int N;

        void SetSimStartTime() {
            time(&simstart);
        }

        void SetSimEndTime() {
            time(&simend);

            simdiff = simend - simstart;
        }

        void WriteTimeToFile() {
            FILE *fptimelog;
            char fname[64];

            snprintf(fname, sizeof(fname),"%d_Particle_RSim_Log.txt", N);
            FILE *fp_log = fopen(fname, "w");

            fprintf(fp_log, "UNIX Time Start: %ld\n", (long)simstart);
            fprintf(fp_log, "UNIX Time End: %ld\n", (long)simend);
            fprintf(fp_log, "UNIX Time Diff: %ld\n", (long)simdiff);

            fclose(fp_log);
        }
};

void heartbeat(struct reb_simulation* rsim) {
    if (reb_simulation_output_check(rsim, 1.0)){
        //Built-in progress line: N_tot, t, dt, cpu time, and % of tmax done
        reb_simulation_output_timing(rsim, 100.0);

        //printf("%f ", rsim->particles[2].ax);
        //printf("%f\n", rsim->particles[2].ay);


    }
}

int main() {
    int N;

    TimeLogger timelogger{};

    struct reb_simulation* rsim = reb_simulation_create();
    reb_simulation_set_integrator(rsim, "whfast");
    rsim->heartbeat = &heartbeat;

    //Central Mass Particle.
    struct reb_particle central = {0};
    central.m = 1.0;    //I'm assuming this is Solar Mass Units.
    reb_simulation_add(rsim, central);

    //Orbiting Massless Paricles
    printf("Awaiting Input:");
    scanf("%d", &N);
    timelogger.N = N;

    for (int i = 0; i < N; i++) {
        double a = reb_random_uniform(rsim, 0.4,20.);
        double e = reb_random_uniform(rsim, 0.01,0.2);
        double omega = reb_random_uniform(rsim, 0.,2.*M_PI);
        double f = reb_random_uniform(rsim, 0.,2.*M_PI);

        struct reb_particle p = reb_particle_from_orbit(1.,rsim->particles[0],0.,a,e,0.,0.,omega,f);
        reb_simulation_add(rsim, p);

        printf("Adding particle %d\n", i);
    }

    printf("\n");

    timelogger.SetSimStartTime();

    //Sim Integration for 100s: dt=0.01s thus 10000 steps
    reb_simulation_integrate(rsim, 100);
    reb_simulation_move_to_com(rsim);

    timelogger.SetSimEndTime();

    printf("\n\n");
    printf("Writing Time Info to Log");
    timelogger.WriteTimeToFile();

    printf("\n");
    printf("UNIX Time Start: %ld\n", (long)timelogger.simstart);
    printf("UNIX Time End: %ld\n", (long)timelogger.simend);
    printf("Total Time: %lld\n\n", timelogger.simdiff);

    reb_simulation_output_orbits(rsim, "output");

    reb_simulation_free(rsim);

    system("pause");
}

CMakeLists:

This project uses the rebound library as compiled from source in IntellaJ CLion rather than as provided by rebound's github, as CLion uses CMakeLists instead of Make, new compilation files had to be written.

GalaxySim/src/CMakeLists.txt:

cmake_minimum_required(VERSION 3.31)
project(rebound C)

set(SOURCES rebound.c
        tree.c
        particle.c
        gravity.c
        integrator_whfast.c
        integrator_whfast512.c
        integrator_saba.c
        integrator_ias15.c
        integrator_sei.c
        integrator_bs.c
        integrator_leapfrog.c
        integrator_mercurius.c
        integrator_trace.c
        integrator_eos.c
        boundary.c
        binarydata.c
        output.c
        collision.c
        communication_mpi.c
        display.c
        tools.c
        rotations.c
        simulation.c
        derivatives.c
        simulationarchive.c
        glad.c
        integrator_janus.c
        transformations.c
        fmemopen.c
        server.c
        frequency_analysis.c)

#object file build target.
add_library(rebound OBJECT ${SOURCES})
target_compile_definitions(rebound PRIVATE BUILDINGLIBREBOUND _GNU_SOURCE)

#dll file build target.
add_library(librebound SHARED $<TARGET_OBJECTS:rebound>)
target_include_directories(librebound PUBLIC ${CMAKE_CURRENT_SOURCE_DIR})
target_compile_definitions(librebound PRIVATE _USRDLL _WINDLL)

GalaxySim/CMakeLists.txt:

cmake_minimum_required(VERSION 3.31)
project(GalaxySim)

add_subdirectory(${PROJECT_SOURCE_DIR}/src)

add_executable(GalaxySim GalaxySim.cpp)
target_link_libraries(GalaxySim PRIVATE librebound)
#again, have to manually copy over .dll so the exe works.

Benchmark Testing and Results:

Initial benchmark tests returned results consistent with an exponential rise in computing time, quadratic regression was then used to a curve fit, as provided by MATLAB's polyfit and polyval functions.

The program was then run on a friend's considerably more powerful PC and results were compared.

Amount UNIX Time Start UNIX Time End Total Time (s) Amount Est Time (s) 50p 1784779993 1784779999 6 10,000p 120429 (1.39 days) 100p 1784780020 1784780037 17 100,000p 12004281 (4.56 months) 150p 1784840470 1784840504 34 1,000,000p 1200042801 (38 yrs) 200p 1784840585 1784840643 58 250p 1784840800 1784840888 88 300p 1784840917 1784841040 123 350p 1784853451 1784853612 161 400p 1784853102 1784853310 208 450p 1784852796 1784853057 261 500p 1784780056 1784780380 324 550p 1784852388 1784852773 385 600p 1784851882 1784852335 453 650p 1784851333 1784851862 529 700p 1784850690 1784851298 608 750p 1784849967 1784850660 693 800p 1784849060 1784849858 798 850p 1784847886 1784848777 891 900p 1784842807 1784843812 1005 950p 1784841104 1784842246 1142 1000p 1784782477 1784783685 1208 2000p 1784834623 1784839457 4834

Blue Line: My PC

Orange Line: Friends PC

Data within table collected from my PC, see matlab code below for friends data.

Matlab Quadratic Regression Code:

x_dat_b = [0,50,100,150,200,250,300,350,400,450,500,550,600,650,700,750,800,850,900,950,1000,2000];
y_dat_b = [0,6,17,34,58,88,123,161,208,261,324,385,453,529,608,693,798,891,1005,1142,1208,4834];

x_dat_j = [0, 100, 200, 300, 500, 700];
y_dat_j = [0, 11, 27, 60, 144, 288];

x_dat_jo = [0,100,200,300,2000];
y_dat_jo = [0,8,26,54,2105];

coeffs = polyfit(x_dat_jo,y_dat_jo,2);

yvect = polyval(coeffs, x_dat_jo);

sse = sum((y_dat_jo-yvect).^2);
sst = sum((y_dat_jo-mean(y_dat_jo)).^2);

r2 = 1 - sse/sst;

x_fun_b = 0:0.01:2500;
y_fun_b = 0.0012*x_fun_b.^2+0.0428*x_fun_b+1.3282;

x_fun_j = 0:0.01:2500;
y_fun_j = 0.00055945*x_fun_j.^2+0.0142*x_fun_j+2.0172;

plot(x_fun_b,y_fun_b);
xlabel('# of Particles');
ylabel('Time (s)');
hold on;

plot(x_fun_j,y_fun_j);

plot(x_dat_b,y_dat_b, 'blax');

plot(x_dat_j,y_dat_j, 'x', 'Color', '#47c4ed');

hold off;