Diffusion simulation. Works, but how to make it fast

[PhysChem]

Hi all!

I wrote a -primitive- diffusion simulator, although it works there is one problem: It is damn slow. (30000 diffusion-step took almost a day )


How it's look like:

There is a 1000*1000*1000 box with 50000 particles. The x and y koordinats are homogeneously distributed (and random of course), BUT the z coordinates are not. Let's cut the box "horizontally" in 1000 plane (so each plane has a single z-coordinate). The planes on the top are very populated, while the planes with low z-coordinate are barely populated. (The population of each horizontal plane can be described with a 2-degree polynom, but it is not important.)
The particles can move only along the x,y or z axis. I don't want to allow the particles to leave the box. Therefore, if a particle's z-coordinate equals 1000/or equals null, it's movement is restricted (it cannot go upwards/ or downwards).
X and Y koordinates are restricted to, I was lazy to make too much if condition, so i used the np.remainder-trick (but only for x and y coordinates). That means if one of my particles can "teleport back" to the other side of my box.

Now let's see the the important part of the code: (I didn't paste here, how the matrix of the original coordinates (k numpy-matrix) was made. It is long and unimportant...)

svv=[[0, 0, 1], [0, 0, -1], [0, 1, 0], [0, -1, 0], [1, 0, 0], [-1, 0, 0]]
felsvv=[[0, 0, 1], [0, 1, 0], [0, -1, 0], [1, 0, 0], [-1, 0, 0]]
lesvv=[[0, 0, -1], [0, 1, 0], [0, -1, 0], [1, 0, 0], [-1, 0, 0]]
lista=[]
o=np.matrix([1000, 1000, 1001])
for j in range(30000):
   for i in range(50000):
      if k[i,2]==0:
         l=np.matrix([random.choice(felsvv)])
      elif k[i,2]==1000:
         l=np.matrix([random.choice(lesvv)])
      else:   
         l=np.matrix([random.choice(svv)])    
      ujsor=k[i,:]+l
      lista.append(ujsor)
   k=np.remainder(np.matrix(np.array(lista)), o)
   lista=[]
print(k)

Obviously a for cycle inside an other for cycle is not good for speed. As you can see, the particles are moving one-by one.
It would be faster, if i could create a big L-matrix (stepping-matrix) with the shape of 50000*3, so all the particles could move at once.
Like this way: L=np.matrix([random.choice(svv) for _in range(50000)])
But in that case, I cannot restrict the z-koordinates, and the particles would escape on the top and the bottom of the box.

[PhysChem]

I solved it :-)

The old code worked the following way:
1.step) the python read the first row of the k-matrix (which contains the coordinates),
2.step) then -according to the last coordinate- it decided what kind of l-matrix (shaped 1*3) could be generated. ( for example if the last coordinate of the k-matrix=999, then l cannot be [0,0,1])
3.step) the new coordinates of the row was calculated (ujsor=k[i,:]+l)
4. step) the new coordinates was put into a list (called lista)
These four steps was repeated until every row in the k-matrix was modified. That was the i-cycle. Then -in order to kept the x and y coordinates between 1 and 1000- np.remainder function was used, Finally the list was cleared, and everything was repeated again and again 29999 times (the j-cycle)

The new code:

svv=[[0,0,1], [0,0,-1], [0,1,0], [0,-1,0], [1,0,0], [-1,0,0]]
lesvv=[[0,0,-1], [0,1,0], [0,-1,0], [1,0,0], [-1,0,0]]
felsvv=[[0,0,1], [0,1,0], [0,-1,0], [1,0,0], [-1,0,0]]
adatlista=[]
for j in range(30000):
 for i in range(50000):
  if k[i,2]==1:
     c=random.choice(felsvv)
     adatlista.append(c)
  elif k[i,2]==1000:
     c=random.choice(lesvv)
     adatlista.append(c)
  else:
     c=random.choice(svv)
     adatlista.append(c)
 L=np.matrix(adatlista)
 k=k+L
 adatlista=[]
print(k)

Instead of row-vector addition, the complete L-matrix (50000*3 shaped) was made, but the python interpreter paid regard of the corresponding row of the k-matrix. After that, the L-matrix added to the k-matrix.

While the old code was able to change the position of every particles about 1500 times per hour, the new code can do the same about 9500 times per hour.

I hope i was understandable, since my english language knowledge is quite sux.

[j.crater]

Hey! It's great to see you got your code optimized well :) And thanks for postng your findings and solution!