Simulation of interacting particles [closed]

Question

I'm simulating a system of particles which interact with each other and move and rotate due to the interaction. There are 2000 particles in total but each particle interacts with other particles and their image (I'm using a periodic boundary condition) if particles and their images have distances larger thanl0. I just consider image of particles in 8 near boxes to the main box.

y[N],x[N] are the position of points in the plane; dx[n] and dy[N] their changes in each time step. c shows different time steps. e_x and e_y shows x and y component of direction of each particle respect to x and y-axis.

int main()
{
  unsigned short int N, l0 , R0;
  double rho, v0, eta, dt, a, b, D;
  long c_equilib, c_prod, c_display;
     while (c < 1000000)
            {
         c++;
        for(unsigned short int i=0; i<2000; i++){
          double S_theta=0.0, S_x=0.0, S_y=0.0, S;
                    for(unsigned short int j=0; j<2000; j++)
                    {
                        if (j==i) continue;
                        for(int ky=-1; ky<=1; ky++)
                        {
                            for(int kx=-1; kx<=1; kx++)
                            {
                                double r_x_ij=x[i]-(x[j]+kx*L),
                                    r_y_ij=y[i]-(y[j]+ky*L),
                                r_ij=sqrt(pow(r_x_ij,2.0)+pow(r_y_ij,2.0));
                                if (r_ij >= l0)
                                { 
                                    r_x_ij/=r_ij;
                                    r_y_ij/=r_ij;
                                    S_theta += (e_x[i]*e_x[j]+e_y[i]*e_y[j]) / pow(r_ij, 3) * ( 3.0* (r_x_ij*e_x[j]+r_y_ij*e_y[j]) *
                                                            pow(e_x[i]*r_y_ij - e_y[i]*r_x_ij, 2.0) -
                                                            e_x[i]*e_y[j] + e_y[i]*e_x[j] );

                                    S = ( 3.0* pow(r_x_ij*e_x[j]+r_y_ij*e_y[j], 2.0) -1.0 ) / pow(r_ij, 2.0);
                                    S_x += r_x_ij*S;
                                    S_y += r_y_ij*S;
                                }
                            }
                        }
                    }
        }
dx[i] += dt*a*S_x;
            dy[i] += dt*a*S_y;
        }
e_x[i]=cos(theta[i]);
e_y[i]=sin(theta[i]);
x[i]+=dx[i];
y[i]+=dy[i];
}
}

I'm looking to optimize the large loop of this C++ program. This is one part of my simulation (I have ignored some calculations inside the loop to make the question simpler). This part is extremely time-consuming. How can I make it faster?

Answer 1

I don't know how fast you expect this to run, but you need to understand that the inner part of your nested loop runs 36 trillion (1000000 * 2000 * 2000 * 3 * 3) times. This means that small changes in the run time of the inner loop will have a large impact on the total run time of the code, turning 1 hour of runtime to years.

Without knowing anything about the domain of your simulation, I can only propose this:

• Reduce the number of times the loops run
• Parallelize the code, so that you can have multiple threads working on different parts concurrently.
• Optimize the inner loop as much as possible. Perhaps there is some approximation of the calculation you're running that's significantly faster? Is there any way that you can make better use of processor caches or avoid branch mispredictions? Have a look at this latency chart.

Answer 2

You can note that i and j are constant in the loops controlled by kx and ky, so the first (micro)optimization could be calculating x-y differences once instead of nine times.

    for(unsigned short i= ...) {
      ...
      for(unsigned short j= ...) {
          if (j==i) continue;

          double delta_x = x[i]-x[j],
                 delta_y = y[i]-y[j];

          for(int ky=-1; ky<=1; ky++) {
              for(int kx=-1; kx<=1; kx++) {
                  double r_x_ij = delta_x - kx*L,
                         r_y_ij = delta_y - ky*L,
                         ....

Next, you calculate (possibly a quite expensive) square root just to discover (sometimes) it should be discarded. If so, compare squares first, then calculate a square root if necessary. Also, you can multiply instead of using pow(..., 2).

                  double r_x_ij = delta_x - kx*L,
                         r_y_ij = delta_y - ky*L,
                         r_ij_square = r_x_ij*r_x_ij + r_x_ij*r_x_ij;
                  if (r_ij_square >= l0_square) {
                      double r_ij = sqrt(r_ij_square);
                      ...

Then you'll no longer need pow(r_ij, 2.0) as the square is ready to use in r_ij_square. You can also replace pow(r_ij, 3) with a single multiplication of r_ij_square * r_ij.

However, all those fixes can not help for an enormous number of iterations which results from your simulated time span of \$1,000,000\$ steps and the square complexity of considering interaction between each two of \$2000\$ particles, which makes \$ 2000\times 2000 = 4,000,000\$ pairs...