分子動力学を使用したイオントラップシミュレーションとOpenGLによる可視化 で行っていたシミュレーション処理のクーロン力計算部分を、前回検討したmultiprocessingを使った並列処理版に置き換えて処理を行ってみた。シングルコアで動作させていた時より、よりスムーズにイオン粒子が動いているのが目視にて確認できる。
並列数4、粒子数は128個。
並列なしのシングルコア動作版の動画が下記(粒子数は同様に128個)
並列化した計算のサンプルソースコードが下記。
#######################################
# ion trap simulation ---------------
# with OpenGL/multiprocessing -------
#######################################
import random
import math
import itertools
import time
import scipy.special as scm
#import scipy.misc as scm
from OpenGL.GL import *
from OpenGL.GLU import *
from OpenGL.GLUT import *
from multiprocessing import Pool
#delta time
Dt = 0.01
step = 0
#Define xyz index
PX = 0;PY = 1;PZ = 2;
VX = 3;VY = 4;VZ = 5;
FX = 6;FY = 7;FZ = 8;
FF = [FX, FY, FZ]
PP = [PX, PY, PZ]
#Coulomb's Force parameter
CF = 1000.0
#mouse move
mx = 0.0
my = 0.0
#lattice size
lsize = 1.0
#number of particles in a line
line_num = 10
#total particle num
PN = line_num * line_num * line_num
th = math.pi * 0.5
ph = math.pi * 0.5
#ready to 9 parameters for particle
#(PX, PY, PZ, VX, VY, VZ, FX, FY, FZ)
xyz = [[0 for i in range(9)] for j in range(PN)]
combinum = int(scm.comb(PN, 2))
#thread number(local thread num)
core = 4
def find_pair_sub(prep,pend,thread):
global xyz
#local results array
xyzF = [[0 for i in range(3)] for j in range(PN)]
fx = 0; fy = 1; fz = 2
for i in range(prep,pend):
for j in range(i + 1, PN):
dx = xyz[i][PX] - xyz[j][PX]
dy = xyz[i][PY] - xyz[j][PY]
dz = xyz[i][PZ] - xyz[j][PZ]
r = math.sqrt(dx*dx + dy*dy + dz*dz)
xyzF[i][fx] = xyzF[i][fx] + dx/(r*r*r)
xyzF[i][fy] = xyzF[i][fy] + dy/(r*r*r)
xyzF[i][fz] = xyzF[i][fz] + dz/(r*r*r)
xyzF[j][fx] = xyzF[j][fx] - dx/(r*r*r)
xyzF[j][fy] = xyzF[j][fy] - dy/(r*r*r)
xyzF[j][fz] = xyzF[j][fz] - dz/(r*r*r)
return xyzF
def wrapper(args):
return find_pair_sub(*args)
def find_pair():
global PN
global combinum
pw = combinum // core
pl = combinum % core
localt = 0
thread = 0
pre = 0
worklist = []
ppp = pw
for i in range(PN) :
if core == 1:
worklist.append([pre,PN,thread])
break
localt = localt + (PN - i - 1)
if localt >= ppp:
worklist.append([pre,i,thread])
ppp += pw
thread += 1
pre = i
if i != pre:
prep = worklist[thread-1][0]
worklist[thread-1] = [prep,PN,thread-1]
p = Pool(core)
callback = p.map(wrapper, worklist)
p.close()
for j in range(core):
for i in range(PN):
xyz[i][FX] += callback[j][i][0]
xyz[i][FY] += callback[j][i][1]
xyz[i][FZ] += callback[j][i][2]
def update():
global step
step += 1
find_pair()
pnum = 0
while pnum < PN:
# set force that collects particles
xyz[pnum][FX] = xyz[pnum][FX] + CF*(lsize/2.0 - xyz[pnum][PX])
xyz[pnum][FY] = xyz[pnum][FY] + CF*(lsize/2.0 - xyz[pnum][PY])
xyz[pnum][FZ] = xyz[pnum][FZ] + CF*(lsize/2.0 - xyz[pnum][PZ])
# if too large force is making , reset.
CUT = 10000
for f in FF:
if xyz[pnum][f] > abs(CUT):
xyz[pnum][f] = 0
xyz[pnum][VX] = xyz[pnum][VX] + 0.5*Dt*xyz[pnum][FX]
xyz[pnum][VY] = xyz[pnum][VY] + 0.5*Dt*xyz[pnum][FY]
xyz[pnum][VZ] = xyz[pnum][VZ] + 0.5*Dt*xyz[pnum][FZ]
xyz[pnum][FX] = 0.0
xyz[pnum][FY] = 0.0
xyz[pnum][FZ] = 0.0
pnum += 1
#print "- step =", step," -"
pnum = 0
while pnum < PN:
#update position
xyz[pnum][PX] = xyz[pnum][PX] + Dt*xyz[pnum][VX]
xyz[pnum][PY] = xyz[pnum][PY] + Dt*xyz[pnum][VY]
xyz[pnum][PZ] = xyz[pnum][PZ] + Dt*xyz[pnum][VZ]
pnum += 1
def init_lattice():
#PN = line_num^3
delt = lsize/line_num;
xv_sum = 0
yv_sum = 0
zv_sum = 0
pnum = 0
k = 0
while k < line_num:
j = 0
while j < line_num:
i = 0
while i < line_num:
xyz[pnum][PX] = delt * i
xyz[pnum][PY] = delt * j
xyz[pnum][PZ] = delt * k
pnum += 1
i += 1
j += 1
k += 1
pnum = 0
while pnum < PN:
xyz[pnum][VX] = random.uniform(-1,1) # 1 x-velocity
xyz[pnum][VY] = random.uniform(-1,1) # 2 y-velocity
xyz[pnum][VZ] = random.uniform(-1,1) # 3 z-velocity
xv_sum += xyz[pnum][VX]
yv_sum += xyz[pnum][VY]
zv_sum += xyz[pnum][VZ]
pnum += 1
xv_sum = xv_sum / PN
yv_sum = yv_sum / PN
zv_sum = zv_sum / PN
pnum = 0
while pnum < PN:
xyz[pnum][VX] = xyz[pnum][VX] - xv_sum
xyz[pnum][VY] = xyz[pnum][VY] - yv_sum
xyz[pnum][VZ] = xyz[pnum][VZ] - zv_sum
pnum += 1
def draw():
cx = cy = cz = lsize * 0.5
glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT)
glClearColor(0.0, 0.0, 0.0, 1.0)
glLoadIdentity()
glColor3f(1.0,1.0,0.0)
gluLookAt(5*math.sin(th)*math.cos(ph)+cx, 5*math.cos(th)+cy, 5*math.sin(th)*math.sin(ph)+cz,\
cx, cy, cz, -math.cos(th)*math.cos(ph),math.sin(th),-math.cos(th)*math.sin(ph) )
glBegin(GL_LINE_LOOP)
glVertex3d(0,0,0)
glVertex3d(0,lsize,0)
glVertex3d(lsize,lsize,0)
glVertex3d(lsize,0,0)
glEnd()
glBegin(GL_LINE_LOOP)
glVertex3d(0,0,lsize)
glVertex3d(0,lsize,lsize)
glVertex3d(lsize,lsize,lsize)
glVertex3d(lsize,0,lsize)
glEnd()
glBegin(GL_LINES)
glVertex3d(0,0,0)
glVertex3d(0,0,lsize)
glEnd()
glBegin(GL_LINES)
glVertex3d(0,lsize,0)
glVertex3d(0,lsize,lsize)
glEnd()
glBegin(GL_LINES)
glVertex3d(lsize,0,0)
glVertex3d(lsize,0,lsize)
glEnd()
glBegin(GL_LINES)
glVertex3d(lsize,lsize,0)
glVertex3d(lsize,lsize,lsize)
glEnd()
pnum = 0
while pnum < PN:
glPointSize(3.0)
glColor3f(0.3, 0.3, 1.0)
glBegin(GL_POINTS)
glVertex3d(xyz[pnum][PX], xyz[pnum][PY], xyz[pnum][PZ])
glEnd()
#Reduce particle velocity
xyz[pnum][VX] *= 0.99
xyz[pnum][VY] *= 0.99
xyz[pnum][VZ] *= 0.99
pnum = pnum + 1
update()
glutSwapBuffers()
def init():
glClearColor(0.7, 0.7, 0.7, 0.7)
def idle():
glutPostRedisplay()
def reshape(w, h):
glViewport(0, 0,w,h)
glMatrixMode(GL_PROJECTION)
glLoadIdentity()
gluPerspective(30.0, w/h, 1.0, 100.0)
glMatrixMode (GL_MODELVIEW)
def motion(x, y):
global ph, th, mx, my
dltx = mx -x
dlty = my -y
#Invalid large motion
if 10 < abs(dltx):
dltx = 0
if 10 < abs(dlty):
dlty = 0
ph = ph - 0.01*dltx
th = th + 0.01*dlty
glutPostRedisplay()
glutSwapBuffers()
mx = x
my = y
if __name__ == "__main__":
init_lattice()
glutInitWindowPosition(50, 50);
glutInitWindowSize(500, 500);
glutInit(sys.argv)
glutInitDisplayMode(GLUT_RGBA | GLUT_DOUBLE )
glutCreateWindow("pyOpenGL TEST")
glutDisplayFunc(draw)
glutReshapeFunc(reshape)
glutMotionFunc(motion);
init()
glutIdleFunc(idle)
glutMainLoop()
I tried processing by replacing the Coulomb force calculation part of simulation processing which was done by ion trap simulation using molecular dynamics and visualization by OpenGL with a parallel processing version using multiprocessing discussed previously . It can be visually confirmed that ion particles are moving more smoothly than when operating with a single core.