A have a problem with the speed of the program written in Python. The program is "Simulation of ferromagnetic particles in a magnetic field", and more specifically in the magnetically inert liquid. The program works, but very slowly compared to the same program written in C++, but I have it written in Python, on a project to study.
Overall, the program code is based on loops, with lots of floating point calculations. There is a random number of particles (generated randomly position) which interact with each other under the influence of a magnetic field.
This is the initial position:
https://i.stack.imgur.com/T15Bb.jpg
And the final:
https://i.stack.imgur.com/0nU5D.jpg
The main loop (in SymMain.py, with k variable) iterations are time steps are calculated coordinates of the particles in them at the time and the forces acting on it (the attractive force and a small repulsive force). In order to speed up, I wanted to use parallel processing to calculate the number of iterations at the same time, but this is not possible, because the calculation in one iteration depend on the calculations in the previous one.
I did not know that Python is so slower than C + +. As an example, the simulation of 529 particles in one time step is calculated (on my computer):
C + + ~ 0.5s
Python ~ 50s
So, how to speed up the program? Please help.
This code is only calculate, not drawing particles as in the pictures. The number of simulate partcile is on SymConst.py, this is nrH*nrL.
SymMain.py
#coding:windows-1250
from os import system
from SymCalc import *
from SymParticle import *
if __name__ == '__main__':
App = SymCalc()
App.MainLoop()
SymParticle.py
#coding:windows-1250
from random import randint
from math import *
from SymConst import *
from SymParticle import *
class SymCalc(object):
def __init__(self):
# declaration lists containing the properties of the particles
ParticleList = []
ParticleListTemp = []
t = 0.0
# the initial values of particle
for x in range(0, nParticle):
ParticleList.append(Particle(x+1, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 8e-6))
ParticleListTemp.append(Particle(x+1, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 8e-6))
# generating random coordinates X, Y
for x in range(0, nParticle):
self.Rand(ParticleList[x], x)
# time steps
for k in range(0, k0):
print "Time step = {0}".format(k)
# calculation of forces
for i in range(0, nParticle):
for j in range(0, i-1):
self.ParticleCalculate(ParticleList[i], ParticleList[j], 1)
for j in range(i+1, nParticle):
self.ParticleCalculate(ParticleList[i], ParticleList[j], 1)
# display data
if(k%Constant == 0):
for i in range(0, nParticle):
self.Print(k, i, ParticleList[i], t)
# changing the position of the particle
for i in range(0, nParticle):
self.ChangeOfPosition (ParticleList[i], dt)
# reset forces
for j in range(0, nParticle):
self.ParticleCalculate(ParticleList[i], ParticleList[j], 0)
self.ParticleCalculate(ParticleList[i], ParticleListTemp[j], 0)
# next time step
t += dt
# random coordinates of the particles
def Rand(self, part, lp):
l = lp%nrL # współrzędna X pola
part.x = (l-nrL/2)*(L0/nrL) + ((randint(0,32767)%100)-50)*(L0/nrL-2*part.r)/99 # X
h = (lp+1-(lp)%nrL)/nrL # Y
part.y = (h-nrH/2)*(H0/nrH) + ((randint(0,32767)%100)-50)*(H0/nrH-2*part.r)/99 # współrzędna Y cząstki
# calculating function of the force acting on the particles
def ParticleCalculate(self, part, part_tmp, p):
# auxiliary variables
# r_4, dx, dy, rr, sum, sx, sy, chi_eff, tmp, frepx, frepy, f0
md = []
if(p == 0):
part.frx = 0
part.fry = 0
part.fx = 0
part.fy = 0
if(p == 1):
# versor coordinates connecting the geometrical center of the particle
dx = part.x - part_tmp.x
dy = part.y - part_tmp.y
# the distance between the geometric means of the particle
rr = sqrt(dx*dx + dy*dy)
if(rr < 0.85*(part.r + part_tmp.r)):
print "ERROR: Invalid distance between the particles! Simulation aborted..."
if(rr >= 10*part.r):
# magnetic dipoles
chi_eff = (3.*(MI_P - 1.))/((MI_P - 1) + 3.)
md.append((4.*MI_0*pi*part.r*part.r*part.r*chi_eff*M_H0)/3.0)
md.append((4.*MI_0*pi*part_tmp.r*part_tmp.r*part_tmp.r*chi_eff*M_H0)/3.0)
tmp = pow(rr,7)
# first member
sum = (5.*(md[0]*part.nx*dx + md[0]*part.ny*dy) * (md[1]*part_tmp.nx*dx + md[1]*part_tmp.ny*dy)) / tmp
sx = sum*dx
sy = sum*dy
tmp = tmp / (rr*rr)
# second member
sum = (md[0]*part.nx*md[1]*part_tmp.nx + md[0]*part.ny*md[1]*part_tmp.ny) / tmp
sx -= sum*dx
sy -= sum*dy
# third member
sx -= (md[0]*(md[1]*part_tmp.nx*dx + md[1]*part_tmp.ny*dy)*part.nx + md[1]*(md[0]*part.nx*dx + md[0]*part.ny*dy)*part_tmp.nx)/tmp
sy -= (md[0]*(md[1]*part_tmp.nx*dx + md[1]*part_tmp.ny*dy)*part.ny + md[1]*(md[0]*part.nx*dx + md[0]*part.ny*dy)*part_tmp.ny)/tmp
# finally
tmp = (-3./(4*pi*MI_0))
sx *= tmp
sy *= tmp
part.fx += sx
part.fy += sy
# short-range repulsive force
tmp = pow(rr,15)
r_4 = pow((part.r+part.r),4)
f0 = 3.*fabs(md[0])*fabs(md[1])/(4.*pi*MI_0*r_4)
frepx = pow((part.r+part.r),15)*dx/(tmp*rr)*f0
frepy = pow((part.r+part.r),15)*dy/(tmp*rr)*f0
part.frx += frepx;
part.fry += frepy;
# change the position of the particle
def ChangeOfPosition(self, part, dt):
part.ddx = 0
part.ddy = 0
part.ddx = 1/(6*pi*part.r*eta)*(part.fx+part.frx)*dt
part.ddy = 1/(6*pi*part.r*eta)*(part.fy+part.fry)*dt
# particles new position value
part.x += part.ddx
part.y += part.ddy
if(part.x < -L0/2):
part.x += L0
elif(part.x > L0/2):
part.x -= L0
elif(part.y < -H0/2):
part.y += H0
elif(part.y > H0/2):
part.y -= H0
# display data
def Print(self, k, i, part, t):
print "-"*50
print "\nParticle {0}".format(i+1)
print "The resultant magnetostatic force fx = {0}".format(part.fx - part.frx)
print "The resultant magnetostatic force fy = {0}\n".format(part.fy - part.fry)
if(i == nParticle-1):
print "\n\t\t---t={0}[s]---".format(t)
print "-"*50
SymParticle.py
#coding:windows-1250
class Particle(object):
# generating a particle properties
def __init__(self, num, x, y, fx, fy, frx, fry, nx, ny, ddx, ddy, r):
self.num = num
self.x = x
self.y = y
self.fx = fx
self.fy = fy
self.frx = frx
self.fry = fry
self.nx = nx
self.ny = ny
self.ddx = ddx
self.ddy = ddy
self.r = r
SymConst.py
#coding:windows-1250
### Constant
M_H0 = 3e4
MI_0 = 12.56e-7
MI_P = 2000
eta = 0.1
k0 = 10001
dt = 1e-6 # time step
H0 = 5.95e-4
L0 = 5.95e-4
nrH = 4
nrL = 4
Constant = 100
nParticle = nrH*nrL # number of particle
