N body simulation in C++

Question

I am trying to implement an OpenMP version of the 2-dimensional n-body simulation.

But there is a problem: I assume each particle's initial velocity and acceleration are zero. When the particles first gather together, they would disperse out in a high speed, and do not gather again.

This doesn't seem to be consistent with the Newton Law, right?

Can someone explain why that happens and how I can fix the error?

Here is part of my code:

/* update one frame */
void update() {
    int i, j;

    omp_set_num_threads(8);
    # pragma omp parallel private(j)
    {

    # pragma omp for schedule(static)
        for ( i = 0; i < g_N; ++i ) {
            g_pv[i].a_x = 0.0;
            g_pv[i].a_y = 0.0;
            for ( j = 0; j < g_N; ++j ) {
                if (i == j)
                    continue;
                double r_2 = pow((g_pv[i].pos_x - g_pv[j].pos_x),2) + pow((g_pv[i].pos_y - g_pv[j].pos_y),2);
                g_pv[i].a_x += (-1) * G * g_pv[j].m * (g_pv[i].pos_x - g_pv[j].pos_x) / (pow(r_2 + e,1.5));
                g_pv[i].a_y += (-1) * G * g_pv[j].m * (g_pv[i].pos_y - g_pv[j].pos_y) / (pow(r_2 + e,1.5));             
            } 
            g_pv[i].v_x += period * g_pv[i].a_x;
            g_pv[i].v_y += period * g_pv[i].a_y;
        }
    # pragma omp for schedule(static)
        for ( int i = 0; i < g_N; ++i ) {
            g_pv[i].pos_x += g_pv[i].v_x * period;
            g_pv[i].pos_y += g_pv[i].v_y * period;
        }
    }
}

Don't worry about the OpenMP, just treat it as an sequential version. OpenMP would not affect the outcome much.

Edit: to clarify, here is the entire code (there might be some errors in this part, but the problem I described should occur in the above code section)

# include <iostream>
# include <fstream>
# include <iomanip>
# include <cmath>
# include <vector>
# include <cstdlib>
# include <omp.h>

# include <GL/glew.h>
# include <GL/freeglut.h>
# include <GL/gl.h>

using namespace std;

/* the size of the opengl window */
# define WIDTH 2000
# define HEIGHT 2000

/* define the global constants */
const double G = 6.67 * pow(10, -11);
// const double G = 6.67;
const double e = 0.00001;
const double period = 1;

/* define the structure of particle */
struct particle
{
    double m;
    double pos_x;
    double pos_y;
    double v_x;
    double v_y;
    double a_x;
    double a_y;

    particle(double m = 0, double pos_x = 0, double pos_y = 0, 
            double v_x = 0, double v_y = 0, double a_x = 0, double a_y = 0)
    {
        this->m         = m;
        this->pos_x     = pos_x;
        this->pos_y     = pos_y;
        this->v_x       = v_x;
        this->v_y       = v_y;
        this->a_x       = a_x;
        this->a_y       = a_y;
    }
};

/* define the global data */
int g_N;                        // number of particles
vector<particle> g_pv;          // particle vector

void setUp();
void update();
void display();

int main(int argc, char ** argv) {

    /* set up the window */
    glutInit(&argc, argv);
    glutInitDisplayMode (GLUT_SINGLE | GLUT_RGB);
    glutInitWindowSize (WIDTH, HEIGHT);
    glutInitWindowPosition (100, 100);
    glutCreateWindow ("openmp");

    /* initialize */
    setUp();

    glutDisplayFunc(display);
    glutMainLoop();

    return 0;
}

/* read the input data */
void setUp() {
    glMatrixMode (GL_PROJECTION);
    glLoadIdentity ();
    /* Sets a 2d projection matrix
     * (0,0) is the lower left corner (WIDTH, HEIGHT) is the upper right */
    glOrtho (0, WIDTH, 0, HEIGHT, 0, 1);
    glDisable(GL_DEPTH_TEST);

    ifstream inFile;
    inFile.open("input_25.txt");

    inFile >> g_N;
    g_pv.resize(g_N);
    for ( int i = 0; i < g_N; ++i )
    {
        inFile >> g_pv[i].m >> g_pv[i].pos_x >> g_pv[i].pos_y
               >> g_pv[i].v_x >> g_pv[i].v_y >> g_pv[i].a_x >> g_pv[i].a_y; 
    }
    inFile.close();
}

/* display in openGL */
void display(void) {
    glClear(GL_COLOR_BUFFER_BIT);
    for(int i = 0; i < g_pv.size(); ++i) {
        /* Get the ith particle */
        particle p = g_pv[i];

        /* Draw the particle as a little square. */
        glBegin(GL_QUADS);
        glColor3f (1.0, 1.0, 1.0);
        glVertex2f(p.pos_x + 2, p.pos_y + 2);
        glVertex2f(p.pos_x - 2, p.pos_y + 2);
        glVertex2f(p.pos_x - 2, p.pos_y - 2);
        glVertex2f(p.pos_x + 2, p.pos_y - 2);
        glEnd();
    }   
    update();
    glutPostRedisplay();
    glFlush ();
}

/* update one frame */
void update() {
    int i, j;

    omp_set_num_threads(8);
    # pragma omp parallel private(j)
    {
        /* compute the force */
        # pragma omp for schedule(static)
            for ( i = 0; i < g_N; ++i ) {
                g_pv[i].a_x = 0.0;
                g_pv[i].a_y = 0.0;
                for ( j = 0; j < g_N; ++j ) {
                    if (i == j)
                        continue;
                    double r_2 = pow((g_pv[i].pos_x - g_pv[j].pos_x),2) + pow((g_pv[i].pos_y - g_pv[j].pos_y),2);
                    g_pv[i].a_x += (-1) * G * g_pv[j].m * (g_pv[i].pos_x - g_pv[j].pos_x) / (pow(r_2 + e,1.5));
                    g_pv[i].a_y += (-1) * G * g_pv[j].m * (g_pv[i].pos_y - g_pv[j].pos_y) / (pow(r_2 + e,1.5));             
                } 
                g_pv[i].v_x += period * g_pv[i].a_x;
                g_pv[i].v_y += period * g_pv[i].a_y;
            }

        /* compute the velocity */
        # pragma omp for schedule(static)
            for ( int i = 0; i < g_N; ++i ) {
                g_pv[i].pos_x += g_pv[i].v_x * period;
                g_pv[i].pos_y += g_pv[i].v_y * period;
            }
    }
}

Answer 1

The first evident problem is that you are using a simple "Taylor integration scheme". The point is that since you are approximating the infinitely small dt, to a finite time difference Δt that is your period, you are expanding the well known motion equations in a Taylor expansion, with infinite terms.... that you are truncating.

In general:

x_t+Δt = x_t + x′_t Δt + 1/2 x″_t Δt² + 1/6 x′″_t Δt³ + O(Δt⁴)

being x′_t the first derivative at time t, the velocity v; x″_t the second derivative, the acceleration a,...; O(Δt⁴) is the order of the error - in this example, we are truncating the expansion at the 3^rd order, and we have a local error of the 4^th.

In your case you are using (Euler method):

x_t+Δt = x_t + v_t Δt + O(Δt²)

and

v_t+Δt = v_t + a_t Δt + O(Δt²)

In this case, since you are stopping the expansion to the first term, this is a first order approximation. you are going to get a local error of order O(Δt²), that corresponds to a global error of order O(Δt). That's how the truncation error behaves - see this reference for more details.

I know very well of two different integration scheme that improve the Euler method:

• Verlet Integration (often used in Molecular Dynamics codes)
• Leapfrog Integration.

and I know of Runge-Kutta methods (find the references from the two cited articles) that are even of higher order, but I have never used them personally.

The order I'm speaking of here is the order of the truncation of the approximation, strictly related to the order of the error performed with the truncation itself.

The Leapfrog method is of second order.

The Verlet method is of third order, with local error O(Δt⁴). This is a third order method even if you don't see any third order derivative because they cancel out in the derivation, as shown in the wikipedia reference.

An higher order integration scheme let you obtain more precise trajectories without shrinking the time-step, Δt.

However the most important properties these integration methods have, missing on the simple Taylor forward integration, whatever is the truncation order, are:

• reversibility
• they are Symplectic, i.e., they do conserve energy

From the references you will find a lot of material to learn it, but a good N-body book would let you learn it in a more structured way.

Just a final note about the very important difference between the two schemes: Verlet method does not give automatically the velocities, that need to evaluated (straightforwardly) as a subsequent step. On the other hand Leapfrog method does evaluate both positions and velocities, but at different times - with this scheme you cannot evaluate both quantities at the same time.

Once you have a good integration scheme you need to worry about what's the maximum error you are able to tolerate, some error analysis is needed at this point, and then you will need to implement multiple-time-scale methods to have accurate integration even in case two objects are too close even for O(Δt⁴) or larger (thinking to gravitational slings).

They can be global (that is, reduce period everywhere) or local (do it just for some partition of the system, where the particles are too close)... but at this point you should be able to go and find more references by yourself.

Answer 2

I'm expanding my comment to an answer (as suggested by Z boson) with several suggestions as to how you might address the problem.

This question really belongs to the Computational Science.SE, as I don't think there is anything wrong with the code per se, but rather the algorithm seems faulty: As the particles get close you can get a force of G / pow(e,1.5) ~ G * 1e7. This is large. Very, very large (compared to your time step). Why? Suppose you have two planets, one massive sitting at (0, 0) and a small one at (0, 1), where the latter gets a very large acceleration towards the former. On the next step the small planet will be at (0, -100), or whatever, and the force on it zero, which means that it will never return and now has a considerable velocity, too. Your simulation has blown up.

As you said, this does not jive well with Newton's laws, so this indicates that your numerical scheme has failed. Worry not, there are several ways to fix this. You already anticipated as much, as you added e. Just make it larger, say 10, and you ought to be fine. Alternatively, set a very small timestep, period. You could also just not compute the force if the particles get too close or have some heuristic as to what should happen if planets collide (maybe the explode and vanish, or throw an exception). Or have a repulsive force say with r ^{- 2} potential or something: g_pv[i].a_y += (-1) * G * g_pv[j].m * (g_pv[i].pos_y - g_pv[j].pos_y) * (1 / pow(r_2 + e,1.5) - 1 / pow(r_2 + e,2.5)); This is similar to how the phenomenological Lennard-Jones interaction incorporates the repulsion arising from the Pauli exclusion principle. Note that you can increase the sharpness of the repulsion (change 2.5 to 12.5, if you like) which means that the repulsion has less of an effect far away (good), but that it needs to be resolved more accurately, leading to a smaller period (bad). Ideally you would start with an initial configuration that does not lead to collisions, but this is nigh-impossible to predict. In the end, you'll probably want to use a combination of the methods that were listed above.

Using OpenMP might lead to a slight speed-up, but you should really be using algorithms for long-range forces, such as Barnes-Hut simulation. For more, see for example the 2013 Summer School on Fast Methods for Long Range Interactions in Complex Particle Systems and the booklet (available for free) they published on the most recent developments. You also might not want to display every time step: Scientific simulations usually save every 5000th step or so. If you want nice looking movies, you can interpolate with some moving averages to smooth out the noise (your simulation does not have temperature or anything of the sort, so you will probably be fine without averaging). Also I'm not sure if your data structures are optimized for the job or if you might be running into cache misses. I'm not really an expert, so maybe someone else might want to weigh in on this. Finally, consider not using pow, but rather the fast inverse square root or similar methods.