# Simulation of interacting particles [closed]

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?

• Just in passing, I'll note that l0 is one of the worst names you could choose - for many readers, it's easily misread as 10. – Toby Speight Oct 19 '17 at 15:44
• This is NOT a working code. I doesn't even compile! It lacks the declaration of many essential variables (coord arrays x[] and y[], to mention the two most prominent), it pretends to use variables which were never assigned (theta[]), it has references to variables outside their scope (i and j used past the and of for(i=...) and for(j=...)), etc, etc, etc... – CiaPan Oct 19 '17 at 20:41
• @OliverRange: c, x, y, e_x, e_y, dx, dy and theta still have no type (and no sample set of values) and there's a closing brace too much. Please provide an example we can compile. – hoffmale Oct 20 '17 at 12:21
• @hoffmale As this question now has answers, we prefer if a new question is asked instead of this question being edited. – Simon Forsberg Oct 20 '17 at 13:03
• Follow-up question – 200_success Oct 21 '17 at 8:52

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.
• I have provided a runnable code. Could you please give me suggestions to optimize it? – Oliver Range Oct 20 '17 at 12:31
• They have removed my updates. Sorry – Oliver Range Oct 20 '17 at 13:08

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...

• Could you please check my updates in the question? I tried to describe it better and complete the code. Could you please give more suggestions for the updated question? – Oliver Range Oct 20 '17 at 12:13
• They have removed my updates. Sorry – Oliver Range Oct 20 '17 at 13:08
• Could you please check my new question here ?codereview.stackexchange.com/questions/178379/… – Oliver Range Oct 20 '17 at 13:16