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

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.