Optimizing gravitation calculation for particles i

Question

I've have created a small visualisation of particles in python. I'm caclulation the movement of particels in a 2D space with zero gravity. As each particle attracts all other particles based on the particle mass and distance.

I made a visualsation in pygame and everything works as plan (with the caluclation), however i need to optimize the calculation extreamly. Today the system can calculate about 100-150 particles in a deacent framerate. I put all the calculation in a seperate thread that gave me some more but not nearly what i want.

I've looked at scipy and numpy but since I'm no scientist or mathguru i just get confused. It looks like I'm on the right track but i have no clue howto.

I need to calculate all the attraction on all particles i have to a loop in a loop. And since I need to find if any have collided, i have to do the same all over again.

It breaks my heart to write that kind of code....

Numpy has the ability to calculate array with array, however i haven't found any what to calculate all items in array with all the items from same/another array. Is there one? If so i could create and couple of arrays and calculate much faster and there must be a function to get index from to 2 arrays where their values match (Collitiondetect iow)

Here is todays attraction/collsion calculation:

class Particle:
    def __init__(self):
        self.x = random.randint(10,790)
        self.y = random.randint(10,590)
        self.speedx = 0.0
        self.speedy = 0.0
        self.mass = 4

#Attraction    
for p in Particles:
    for p2 in Particles:
        if p != p2:
            xdiff = P.x - P2.x
            ydiff = P.y - P2.y
            dist = math.sqrt((xdiff**2)+(ydiff**2))
            force = 0.125*(p.mass*p2.mass)/(dist**2)
            acceleration = force / p.mass
            xc = xdiff/dist
            yc = ydiff/dist
            P.speedx -= acceleration * xc
            P.speedy -= acceleration * yc
for p in Particles:
    p.x += p.speedx
    p.y += p.speedy

#Collision
for P in Particles:
   for P2 in Particles:
        if p != P2:
            Distance = math.sqrt(  ((p.x-P2.x)**2)  +  ((p.y-P2.y)**2)  )
            if Distance < (p.radius+P2.radius):
                p.speedx = ((p.mass*p.speedx)+(P2.mass*P2.speedx))/(p.mass+P2.mass)
                p.speedy = ((p.mass*p.speedy)+(P2.mass*P2.speedy))/(p.mass+P2.mass)
                p.x = ((p.mass*p.x)+(P2.mass*P2.x))/(p.mass+P2.mass)
                p.y = ((p.mass*p.y)+(P2.mass*P2.y))/(p.mass+P2.mass)
                p.mass += P2.mass
                p.radius = math.sqrt(p.mass)
                Particles.remove(P2)

Answer 1

I've worked on this previously, and one of the things I've seen in the past to accelerate collision calculations is to actually store a list of nearby particles.

Basically, the idea is inside of your gravity calculation you do something like:

for (int i = 0; i < n; i++)
{
    for (int j = i + 1; j < n; j++)
    {
        DoGravity(Particle[i], Particle[j]);
        if (IsClose(Particle[i], Particle[j]))
        {
            Particle[i].AddNeighbor(Particle[j]);
            Particle[j].AddNeighbor(Particle[i]);
        }
    }
}

Then, you simply pass over all particles and you do collision detection on each on in turn. This is usually something like O(n) in best case, but it can easily degrade to O(n^2) in the worst case.

Another alternative is to try placing your particles inside of a Octree. Building one up is something like O(n), then you can query it to see if anything is near each other. At that point you'd just do collision detection on the pairs. Doing this is O(n log n) I believe.

Not only that, but you can use the Octree to accelerate the gravity calculation as well. Instead of O(n^2) behavior, it drops down to O(n log n) as well. Most Octree implementations include an "opening parameter" that controls the speed vs accuracy trade off you'll be making. So Octrees tend to be less accurate than a direct pairwise calculation and complicated to code up, but they also make large scale simulations possible.

If you use the Octree in this manner, you'll do what's known as a Barnes-Hut Simulation.

Note: Since you're working in 2D, the 2D analogue to an Octree is known as a Quadtree. See the following Wikipedia article for more information: http://en.wikipedia.org/wiki/Quadtree

Answer 2

You can first try to work with complex numbers: the relevant gravitation and dynamics formulas are very simple in this formalism, and can also be quite fast (because NumPy can do the calculation internally, instead of you handling separately x and y coordinates). For instance, the force between two particules at z and z' is simply:

(z-z')/abs(z-z')**3

NumPy can calculate such a quantity very quickly, for all z/z' pairs. For instance, the matrix of all z-z' values is simply obtained from the 1D array Z of coordinates as Z-Z[:, numpy.newaxis] (the diagonal terms [z=z'] do require some special care, when calculating 1/abs(z-z')**3: they should be set to zero).

As for the time evolution, you can certainly use SciPy's fast differential equation routines: they are much more precise than the step by step Euler integration.

In any case, delving into NumPy would be very useful, especially if you plan to do scientific calculations, as NumPy is very fast.

Answer 3

Not sure if this will help you out much but it's part of a solution I've been working on for the same problem. I didn't notice a huge gain in performance doing it this way, still starting to grind to a halt around 200 particles, but maybe it'll give you some ideas.

C++ module for calculating the x and y components of gravitational attraction on a 2d plane:

#include <Python.h>
#include <math.h>

double _acceleration(double &Vxa, double &Vya, double &Vxb, double &Vyb, double xa, double ya, double xb, double yb, double massa, double massb)
{
   double xdiff = xa - xb;
   double ydiff = ya - yb;
   double distance = sqrt(xdiff*xdiff + ydiff*ydiff) * pow(10, 5);

   if (distance <= 0)
      distance = 1;

   double force = (6.674 * pow(10, -11))*(massa*massb)/(distance*distance);

   double acca = force / massa;
   double accb = force / massb;
   double xcomponent = xdiff/distance;
   double ycomponent = ydiff/distance;

   Vxa -= acca * xcomponent;
   Vya -= acca * ycomponent;
   Vxb += accb * xcomponent;
   Vyb += accb * ycomponent;

   return distance;
}

static PyObject* gforces(PyObject* self, PyObject* args)
{
   double Vxa, Vya, Vxb, Vyb, xa, ya, xb, yb, massa, massb, distance;

   if (!PyArg_ParseTuple(args, "ffffdffffdffffdd", &Vxa, &Vya, &Vxb, &Vyb, &xa, &ya, &xb, &yb, &massa, &massb))
      return NULL;

   distance = _acceleration(Vxa, Vya, Vxb, Vyb, xa, ya, xb, yb, massa, massb);

   return Py_BuildValue("ffffddd", Vxa, Vya, Vxb, Vyb, distance);
}

static PyMethodDef GForcesMethods[] = {
{"gforces", gforces, METH_VARARGS, "Calculate the x and y acceleration of two masses and the distance between the two."},
{NULL, NULL, 0, NULL}
};

PyMODINIT_FUNC
initgforces(void)
{
(void) Py_InitModule("gforces", GForcesMethods);
}

If you compile this as a pyd file you should get a python object that you can import. You do have to get all your compiler and linker options correct though. I'm using dev-C++ and have my compiler options set to -shared -o gforces.pyd and linker set to -lpython27 (make sure you use the same version you have installed) and add your python directory path to the includes and libraries tabs.

The object takes the arguments ( p1.speedx, p1.speedy, p2.speedx, p2.speedy, p1.x, p1.y, p2.x, p2.y, p1.mass, p2.mass ) and returns the new p1.speedx, p1.speedy, p2.speedx, p2.speedy, and distance between p1 p2.

Using the above module, I've also tried to cut out a few steps for the collision detection by comparing the returned distance to the sum of the particles' radii as such:

def updateForces(self):         #part of a handler class for the particles
    prevDone = []
    for i in self.planetGroup:  #planetGroup is a sprite group from pygame
        prevDone.append(i.ID)
        for j in self.planetGroup:
            if not (j.ID in prevDone):               #my particles have unique IDs
                distance = i.calcGForce( j )         #calcGForce just calls the above  
                if distance <= i.radius + j.radius:  #object and assigns the returned 
                    #collision handling goes here    #values for the x and y speed
                                                     #components to the particles

Hope this helps a bit. Any further advice or pointing out of gross errors on my part is welcome, I'm new to this as well.

Answer 4

(This may should go in a comment but I don't have the needed reputation to do that)

I don't see how you do the time stepping. You have

P.speedx -= acceleration * xc
P.speedy -= acceleration * yc

but to get the new speed at time t+delta_t you would do

P.speedx -= acceleration * xc * delta_t
P.speedy -= acceleration * yc * delta_t

and then update the position like so:

P.x = P.x + P.speedx * delta_t
P.y = P.y + P.speedy * delta_t

Then to your speed concern. Maybe it would be better to store the particle information not in a class but in numpy arrays? But I don't think you can avoid loops.

Also, have you looked at wikipedia, there it describes some methods to speed up the calculation.

(edited due to Mike's comment)

Answer 5

to do fast calculation, you need to store x, y, speedx, speedy, m in numpy arrays. For example:

import numpy as np

p = np.array([
    (0,0),
    (1,0),
    (0,1),
    (1,1),
    (2,2),
], dtype = np.float)

p is a 5x2 array which store x, y position of particles. To calculate the distance between each pair, you can use:

print np.sqrt(np.sum((p[:, np.newaxis] - p[np.newaxis, :])**2, axis=-1))

the output is:

[[ 0.          1.          1.          1.41421356  2.82842712]
 [ 1.          0.          1.41421356  1.          2.23606798]
 [ 1.          1.41421356  0.          1.          2.23606798]
 [ 1.41421356  1.          1.          0.          1.41421356]
 [ 2.82842712  2.23606798  2.23606798  1.41421356  0.        ]]

or you can use cdist from scipy:

from scipy.spatial.distance import cdist
print cdist(p, p)