n-bodies: a parallel TBB solution: realizing addAcc(i,j)

Having settled the question of whether I should accumulate forces or accelerations last time, now it’s time to build the accumulation function.

    addAcc(int i, int j) {


i and j are indices selecting elements of the body array. First task is to compute the distance between them.

        double dx = body[i].pos[0]-body[j].pos[0];
        double dy = body[i].pos[1]-body[j].pos[1];
        double dz = body[i].pos[2]-body[j].pos[2];
        double distsq = dx*dx + dy*dy + dz*dz;


Pythagorean Theorem in three dimensions gets me the square of the hypotenuse, but before doing the square root, I’ll avoid the singularity:

        if (distsq < MINDIST) distsq = MINDIST;
        double dist = sqrt(distsq);


That is, if the point masses get too close together, act like they’re not. But wait! Why do I even need the square root, if I’m working with gravitation, an inverse-squared law? Well, because acceleration is a vector so I need the next step.

        double ud[3];
        ud[0] = dx/dist;
        ud[1] = dy/dist;
        ud[2] = dz/dist;


Array ud represents the unit vector (length 1 direction vector) pointing from body j to body i. I need just one more thing, the magnitude of those accelerations.

        double Gdivd = GFORCE/distsq;
        double ai = Gdivd*body[j].mass;
        double aj = Gdivd*body[i].mass;


All that’s left is to compute the acceleration vector components and apply them to the bodies.

        for (int k = 0; k < 3; ++k) {
            body[j].acc[k] += aj*ud[k];
            body[i].acc[k] -= ai*ud[k];


Next time: serial bodies test run

For more complete information about compiler optimizations, see our Optimization Notice.