2022年11月 – Python in the box

11月 17, 2022

The process of brute force calculation by Coulomb force with multiprocessing and numpy. TEST3

Each result returned by the distributed processing is not made into a PN-sized matrix, but is placed as is. This lightens the processing.

import random
import math
import time
import gc
import numpy as np

from multiprocessing import Pool

np.set_printoptions(threshold=np.inf)
random.seed(1)

XX = 0
YY = 1
ZZ = 2

core = 6

#total particle num
PN = 10000
#---------------------------------------------------
#
# parameter : divide_N
#
# The number of particles per side of a small square 
# divided for the splitting process when creating a 
# large PN^2 region with all particles PN.
#
# The smaller the setting, the more it will be divided
# into many smaller jobs.
#---------------------------------------------------
divide_N = 1000
xyzV = [[0 for i in range(3)] for j in range(PN)]
xyzP = np.zeros((PN,3))
xyzF = np.zeros((PN,3))


def wrapper(args):
  return calc_XY_start_end(*args)

def find_pair():
  global core
  global PN
  global xyzF
  global xyzP

  worklist = []
  thread   = 0

  remainder = PN %  divide_N
  quotient  = PN // divide_N

  if remainder != 0:                 #ex PN=11, divide_N=4, remainder=3, quotient=2
      quotient = quotient + 1        #   2 -> 3
  extra = divide_N - remainder       #   1

  X_st = 0
  X_ed = 0
  Y_st = 0
  Y_ed = 0

  Thead_num = 0
  sum_f = np.zeros((PN,PN,3))

  for Xn in range(0,quotient):
    for Yn in range(0,Xn+1):

      worklist.append([Xn, Yn, quotient, remainder, divide_N, extra])
      Thead_num = Thead_num + 1
      #----------------------------------------
      #local_force = calc_XY_start_end(Xn, Yn, quotient, remainder, divide_N, extra)
      #----------------------------------------


  #start thread. results is callback in array.
  p = Pool(core)
  callback = p.map(wrapper, worklist)
  p.close()


  t1 = 0
  for Th in range(0,Thead_num):

    xn_l = callback[Th][0]
    yn_l = callback[Th][1]
    xxst = callback[Th][2]
    xxed = callback[Th][3]
    yyst = callback[Th][4]
    yyed = callback[Th][5]
    local_force = callback[Th][6]

    time_sta = time.time()

    sb = local_force.shape
    ten   = -1 * local_force.transpose(1,0,2)
    tensb = ten.shape

    sum_f[yyst:yyst+sb[0], xxst:xxst+sb[1]] += local_force

    if xn_l != yn_l:
      sum_f[xxst:xxst+tensb[0], yyst:yyst+tensb[1]] += ten

    time_end = time.time()
    t1 =t1 + (time_end  - time_sta)


  print ("=====t1=",t1)
  print(np.sum(sum_f))
  #mmm = sum_f.sum(axis=1) 
  #xyzF = xyzF + mmm 


def calc_XY_start_end(Xn, Yn, quotient, remainder, divide_N, extra):
  global PN
  global xyzP

  if Xn != quotient-1 and remainder != 0 or remainder == 0:
    X_st = Xn*divide_N
    X_ed = (Xn+1)*divide_N-1
  else:
    X_st = Xn*divide_N - (divide_N - remainder)
    X_ed = (Xn+1)*divide_N  - (divide_N - remainder) -1
  if Yn != quotient-1 and remainder != 0 or remainder == 0:
    Y_st = Yn*divide_N
    Y_ed = (Yn+1)*divide_N-1
  else:
    Y_st = Yn*divide_N - (divide_N - remainder)
    Y_ed = (Yn+1)*divide_N  -(divide_N - remainder) -1

  Xa = np.arange(X_st,X_ed+1)
  Ya = np.arange(Y_st,Y_ed+1)
  mx1,my1 = np.meshgrid(Xa,Ya)
  local_xyz  = xyzP[mx1]-xyzP[my1]
  distance_a = np.linalg.norm(local_xyz, axis=2)

  tmp_d = distance_a[:,:,np.newaxis]
  tmp_f = local_xyz/tmp_d
  tmp_f = np.nan_to_num(tmp_f,0)

  #extra clear
  if remainder != 0:
    if Xn + 1 == quotient:
      for dd in range(0,extra):
        tmp_f[:,dd] = 0
    if Yn + 1 == quotient:
      for dd in range(0,extra):
        tmp_f[[dd]] = 0

  return Xn,Yn,X_st,X_ed,Y_st,Y_ed,tmp_f



def init_lattice():
  global xyzF
  global xyzP

  pnum = 0
  while pnum < PN:
    xyzP[pnum][XX] = random.uniform(-1,1)
    xyzP[pnum][YY] = random.uniform(-1,1)
    xyzP[pnum][ZZ] = random.uniform(-1,1)
    xyzF[pnum][XX] = random.uniform(-1,1)
    xyzF[pnum][YY] = random.uniform(-1,1)
    xyzF[pnum][ZZ] = random.uniform(-1,1)
    pnum += 1


def results_sum():
  global xyzF
  pnum = 0
  total_F = 0
#  while pnum < PN:
#    total_F = total_F + xyzF[pnum][XX]
#    total_F = total_F + xyzF[pnum][YY]
#    total_F = total_F + xyzF[pnum][ZZ]
#    pnum += 1
#
#  print (total_F)


if __name__ == "__main__":

  init_lattice()
  find_pair()
  results_sum()

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

144

145

146

147

148

149

150

151

152

153

154

155

156

157

158

159

160

161

162

163

164

165

166

167

168

169

170

171

172

173

174

175

176

177

178

179

180

181

182

183

184

import random

import math

import time

import gc

import numpy as np

from multiprocessing import Pool

np.set_printoptions(threshold=np.inf)

random.seed(1)

XX = 0

YY = 1

ZZ = 2

core = 6

#total particle num

PN = 10000

#---------------------------------------------------

# parameter : divide_N

# The number of particles per side of a small square

# divided for the splitting process when creating a

# large PN^2 region with all particles PN.

# The smaller the setting, the more it will be divided

# into many smaller jobs.

#---------------------------------------------------

divide_N = 1000

xyzV = [[0 for i in range(3)] for j in range(PN)]

xyzP = np.zeros((PN,3))

xyzF = np.zeros((PN,3))

def wrapper(args):

return calc_XY_start_end(*args)

def find_pair():

global core

global PN

global xyzF

global xyzP

worklist = []

thread = 0

remainder = PN % divide_N

quotient = PN // divide_N

if remainder != 0: #ex PN=11, divide_N=4, remainder=3, quotient=2

quotient = quotient + 1 # 2 -> 3

extra = divide_N - remainder # 1

X_st = 0

X_ed = 0

Y_st = 0

Y_ed = 0

Thead_num = 0

sum_f = np.zeros((PN,PN,3))

for Xn in range(0,quotient):

for Yn in range(0,Xn+1):

worklist.append([Xn, Yn, quotient, remainder, divide_N, extra])

Thead_num = Thead_num + 1

#----------------------------------------

#local_force = calc_XY_start_end(Xn, Yn, quotient, remainder, divide_N, extra)

#----------------------------------------

#start thread. results is callback in array.

p = Pool(core)

callback = p.map(wrapper, worklist)

p.close()

t1 = 0

for Th in range(0,Thead_num):

xn_l = callback[Th][0]

yn_l = callback[Th][1]

xxst = callback[Th][2]

xxed = callback[Th][3]

yyst = callback[Th][4]

yyed = callback[Th][5]

local_force = callback[Th][6]

time_sta = time.time()

sb = local_force.shape

ten = -1 * local_force.transpose(1,0,2)

tensb = ten.shape

sum_f[yyst:yyst+sb[0], xxst:xxst+sb[1]] += local_force

if xn_l != yn_l:

sum_f[xxst:xxst+tensb[0], yyst:yyst+tensb[1]] += ten

time_end = time.time()

t1 =t1 + (time_end - time_sta)

print ("=====t1=",t1)

print(np.sum(sum_f))

#mmm = sum_f.sum(axis=1)

#xyzF = xyzF + mmm

def calc_XY_start_end(Xn, Yn, quotient, remainder, divide_N, extra):

global PN

global xyzP

if Xn != quotient-1 and remainder != 0 or remainder == 0:

X_st = Xn*divide_N

X_ed = (Xn+1)*divide_N-1

else:

X_st = Xn*divide_N - (divide_N - remainder)

X_ed = (Xn+1)*divide_N - (divide_N - remainder) -1

if Yn != quotient-1 and remainder != 0 or remainder == 0:

Y_st = Yn*divide_N

Y_ed = (Yn+1)*divide_N-1

else:

Y_st = Yn*divide_N - (divide_N - remainder)

Y_ed = (Yn+1)*divide_N -(divide_N - remainder) -1

Xa = np.arange(X_st,X_ed+1)

Ya = np.arange(Y_st,Y_ed+1)

mx1,my1 = np.meshgrid(Xa,Ya)

local_xyz = xyzP[mx1]-xyzP[my1]

distance_a = np.linalg.norm(local_xyz, axis=2)

tmp_d = distance_a[:,:,np.newaxis]

tmp_f = local_xyz/tmp_d

tmp_f = np.nan_to_num(tmp_f,0)

#extra clear

if remainder != 0:

if Xn + 1 == quotient:

for dd in range(0,extra):

tmp_f[:,dd] = 0

if Yn + 1 == quotient:

for dd in range(0,extra):

tmp_f[[dd]] = 0

return Xn,Yn,X_st,X_ed,Y_st,Y_ed,tmp_f

def init_lattice():

global xyzF

global xyzP

pnum = 0

while pnum < PN:

xyzP[pnum][XX] = random.uniform(-1,1)

xyzP[pnum][YY] = random.uniform(-1,1)

xyzP[pnum][ZZ] = random.uniform(-1,1)

xyzF[pnum][XX] = random.uniform(-1,1)

xyzF[pnum][YY] = random.uniform(-1,1)

xyzF[pnum][ZZ] = random.uniform(-1,1)

pnum += 1

def results_sum():

global xyzF

pnum = 0

total_F = 0

# while pnum < PN:

# total_F = total_F + xyzF[pnum][XX]

# total_F = total_F + xyzF[pnum][YY]

# total_F = total_F + xyzF[pnum][ZZ]

# pnum += 1

# print (total_F)

if __name__ == "__main__":

init_lattice()

find_pair()

results_sum()

11月 16, 202211月 16, 2022

The process of brute force calculation by Coulomb force with multiprocessing and numpy. TEST2

The split process reduced the matrix size as much as possible, but it still eats up too much memory.
Multiprocessing consumes a lot of memory.

import random
import math
import time
import gc
import numpy as np

from multiprocessing import Pool

np.set_printoptions(threshold=np.inf)
random.seed(1)

XX = 0
YY = 1
ZZ = 2

#total particle num
PN = 1200
#---------------------------------------------------
#
# parameter : divide_N
#
# The number of particles per side of a small square 
# divided for the splitting process when creating a 
# large PN^2 region with all particles PN.
#
# The smaller the setting, the more it will be divided
# into many smaller jobs.
#---------------------------------------------------
divide_N = 50
core = 2
xyzV = [[0 for i in range(3)] for j in range(PN)]
xyzP = np.zeros((PN,3))
xyzF = np.zeros((PN,3))

def wrapper(args):
  return calc_XY_start_end(*args)

def find_pair():
  global PN
  global xyzF
  global xyzP

  worklist = []
  thread   = 0

  remainder = PN %  divide_N
  quotient  = PN // divide_N

  if remainder != 0:                 #ex PN=11, divide_N=4, remainder=3, quotient=2
      quotient = quotient + 1        #   2 -> 3
  extra = divide_N - remainder       #   1

  X_st = 0
  X_ed = 0
  Y_st = 0
  Y_ed = 0

  Thead_num = 0
  sum_f = np.zeros((PN,PN,3))

  for Xn in range(0,quotient):
    for Yn in range(0,Xn+1):

      worklist.append([Xn, Yn, quotient, remainder, divide_N, extra])
      Thead_num = Thead_num + 1
      #----------------------------------------
      #local_force = calc_XY_start_end(Xn, Yn, quotient, remainder, divide_N, extra)
      #----------------------------------------

  p = Pool(core)

  #start thread. results is callback in array.
  callback = p.map(wrapper, worklist)
  p.close()


  for Th in range(0,Thead_num):

    xn_l = callback[Th][0]
    yn_l = callback[Th][1]
    xxst = callback[Th][2]
    xxed = callback[Th][3]
    yyst = callback[Th][4]
    yyed = callback[Th][5]
    local_force = callback[Th][6]

    if xxst != 0:
      r2 = np.zeros((xxst,divide_N,3))
      local_force        = np.insert(local_force     ,0,r2,axis=1)
    if xxed+1 != PN:
      r2 = np.zeros((PN - xxed - 1,divide_N,3))
      local_force      = np.insert(local_force     ,xxed + 1,r2,axis=1)

    if yyst != 0:
      r2 = np.zeros((yyst,PN,3))
      local_force      = np.insert(local_force     ,0,r2,axis=0)
    if yyed+1 != PN:
      r2 = np.zeros((PN - yyed - 1,PN,3))
      local_force      = np.insert(local_force     ,yyed + 1,r2,axis=0)

    del r2

    if xn_l != yn_l:
      sum_f = sum_f + local_force
      sum_f = sum_f - local_force.transpose(1,0,2)
    else:
      sum_f = sum_f + local_force

    del local_force
    gc.collect()

  #print (sum_f)

  print("=====================================")
  mmm = sum_f.sum(axis=1)
  xyzF = xyzF + mmm


def calc_XY_start_end(Xn, Yn, quotient, remainder, divide_N, extra):
  global PN
  global xyzP

  if Xn != quotient-1 and remainder != 0 or remainder == 0:
    X_st = Xn*divide_N
    X_ed = (Xn+1)*divide_N-1
  else:
    X_st = Xn*divide_N - (divide_N - remainder)
    X_ed = (Xn+1)*divide_N  - (divide_N - remainder) -1
  if Yn != quotient-1 and remainder != 0 or remainder == 0:
    Y_st = Yn*divide_N
    Y_ed = (Yn+1)*divide_N-1
  else:
    Y_st = Yn*divide_N - (divide_N - remainder)
    Y_ed = (Yn+1)*divide_N  -(divide_N - remainder) -1

  Xa = np.arange(X_st,X_ed+1)
  Ya = np.arange(Y_st,Y_ed+1)
  mx1,my1 = np.meshgrid(Xa,Ya)
  local_xyz  = xyzP[mx1]-xyzP[my1]
  distance_a = np.linalg.norm(local_xyz, axis=2)

  tmp_d = distance_a[:,:,np.newaxis]
  tmp_f = local_xyz/tmp_d
  tmp_f = np.nan_to_num(tmp_f,0)

  #extra clear
  if remainder != 0:
    if Xn + 1 == quotient:
      for dd in range(0,extra):
        tmp_f[:,dd] = 0
    if Yn + 1 == quotient:
      for dd in range(0,extra):
        tmp_f[[dd]] = 0

  return Xn,Yn,X_st,X_ed,Y_st,Y_ed,tmp_f



def init_lattice():
  global xyzF
  global xyzP

  pnum = 0
  while pnum < PN:
    xyzP[pnum][XX] = random.uniform(-1,1)
    xyzP[pnum][YY] = random.uniform(-1,1)
    xyzP[pnum][ZZ] = random.uniform(-1,1)
    xyzF[pnum][XX] = random.uniform(-1,1)
    xyzF[pnum][YY] = random.uniform(-1,1)
    xyzF[pnum][ZZ] = random.uniform(-1,1)
    pnum += 1


def results_sum():
  global xyzF
  pnum = 0
  total_F = 0
  while pnum < PN:
    total_F = total_F + xyzF[pnum][XX]
    total_F = total_F + xyzF[pnum][YY]
    total_F = total_F + xyzF[pnum][ZZ]
    pnum += 1

  print (total_F)


if __name__ == "__main__":

  init_lattice()
  find_pair()
  results_sum()

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

144

145

146

147

148

149

150

151

152

153

154

155

156

157

158

159

160

161

162

163

164

165

166

167

168

169

170

171

172

173

174

175

176

177

178

179

180

181

182

183

184

185

186

187

188

189

190

191

import random

import math

import time

import gc

import numpy as np

from multiprocessing import Pool

np.set_printoptions(threshold=np.inf)

random.seed(1)

XX = 0

YY = 1

ZZ = 2

#total particle num

PN = 1200

#---------------------------------------------------

# parameter : divide_N

# The number of particles per side of a small square

# divided for the splitting process when creating a

# large PN^2 region with all particles PN.

# The smaller the setting, the more it will be divided

# into many smaller jobs.

#---------------------------------------------------

divide_N = 50

core = 2

xyzV = [[0 for i in range(3)] for j in range(PN)]

xyzP = np.zeros((PN,3))

xyzF = np.zeros((PN,3))

def wrapper(args):

return calc_XY_start_end(*args)

def find_pair():

global PN

global xyzF

global xyzP

worklist = []

thread = 0

remainder = PN % divide_N

quotient = PN // divide_N

if remainder != 0: #ex PN=11, divide_N=4, remainder=3, quotient=2

quotient = quotient + 1 # 2 -> 3

extra = divide_N - remainder # 1

X_st = 0

X_ed = 0

Y_st = 0

Y_ed = 0

Thead_num = 0

sum_f = np.zeros((PN,PN,3))

for Xn in range(0,quotient):

for Yn in range(0,Xn+1):

worklist.append([Xn, Yn, quotient, remainder, divide_N, extra])

Thead_num = Thead_num + 1

#----------------------------------------

#local_force = calc_XY_start_end(Xn, Yn, quotient, remainder, divide_N, extra)

#----------------------------------------

p = Pool(core)

#start thread. results is callback in array.

callback = p.map(wrapper, worklist)

p.close()

for Th in range(0,Thead_num):

xn_l = callback[Th][0]

yn_l = callback[Th][1]

xxst = callback[Th][2]

xxed = callback[Th][3]

yyst = callback[Th][4]

yyed = callback[Th][5]

local_force = callback[Th][6]

if xxst != 0:

r2 = np.zeros((xxst,divide_N,3))

local_force = np.insert(local_force ,0,r2,axis=1)

if xxed+1 != PN:

r2 = np.zeros((PN - xxed - 1,divide_N,3))

local_force = np.insert(local_force ,xxed + 1,r2,axis=1)

if yyst != 0:

r2 = np.zeros((yyst,PN,3))

local_force = np.insert(local_force ,0,r2,axis=0)

if yyed+1 != PN:

r2 = np.zeros((PN - yyed - 1,PN,3))

local_force = np.insert(local_force ,yyed + 1,r2,axis=0)

del r2

if xn_l != yn_l:

sum_f = sum_f + local_force

sum_f = sum_f - local_force.transpose(1,0,2)

else:

sum_f = sum_f + local_force

del local_force

gc.collect()

#print (sum_f)

print("=====================================")

mmm = sum_f.sum(axis=1)

xyzF = xyzF + mmm

def calc_XY_start_end(Xn, Yn, quotient, remainder, divide_N, extra):

global PN

global xyzP

if Xn != quotient-1 and remainder != 0 or remainder == 0:

X_st = Xn*divide_N

X_ed = (Xn+1)*divide_N-1

else:

X_st = Xn*divide_N - (divide_N - remainder)

X_ed = (Xn+1)*divide_N - (divide_N - remainder) -1

if Yn != quotient-1 and remainder != 0 or remainder == 0:

Y_st = Yn*divide_N

Y_ed = (Yn+1)*divide_N-1

else:

Y_st = Yn*divide_N - (divide_N - remainder)

Y_ed = (Yn+1)*divide_N -(divide_N - remainder) -1

Xa = np.arange(X_st,X_ed+1)

Ya = np.arange(Y_st,Y_ed+1)

mx1,my1 = np.meshgrid(Xa,Ya)

local_xyz = xyzP[mx1]-xyzP[my1]

distance_a = np.linalg.norm(local_xyz, axis=2)

tmp_d = distance_a[:,:,np.newaxis]

tmp_f = local_xyz/tmp_d

tmp_f = np.nan_to_num(tmp_f,0)

#extra clear

if remainder != 0:

if Xn + 1 == quotient:

for dd in range(0,extra):

tmp_f[:,dd] = 0

if Yn + 1 == quotient:

for dd in range(0,extra):

tmp_f[[dd]] = 0

return Xn,Yn,X_st,X_ed,Y_st,Y_ed,tmp_f

def init_lattice():

global xyzF

global xyzP

pnum = 0

while pnum < PN:

xyzP[pnum][XX] = random.uniform(-1,1)

xyzP[pnum][YY] = random.uniform(-1,1)

xyzP[pnum][ZZ] = random.uniform(-1,1)

xyzF[pnum][XX] = random.uniform(-1,1)

xyzF[pnum][YY] = random.uniform(-1,1)

xyzF[pnum][ZZ] = random.uniform(-1,1)

pnum += 1

def results_sum():

global xyzF

pnum = 0

total_F = 0

while pnum < PN:

total_F = total_F + xyzF[pnum][XX]

total_F = total_F + xyzF[pnum][YY]

total_F = total_F + xyzF[pnum][ZZ]

pnum += 1

print (total_F)

if __name__ == "__main__":

init_lattice()

find_pair()

results_sum()

11月 16, 202211月 16, 2022

The process of brute force calculation by Coulomb force with multiprocessing and numpy. TEST1

The small distributed matrices also eat too much memory because they are only PNxPN in size.
Too slow.

###########################
# multiprocessing test ####
###########################
import random
import math
import time
import scipy.special as scm
from multiprocessing import Pool

random.seed(1)

PX = 0;PY = 1;PZ = 2;
VX = 3;VY = 4;VZ = 5;
FX = 6;FY = 7;FZ = 8;

#number of particles in a line
line_num = 15

#total particle num
PN = line_num * line_num * line_num

#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)]

#Number of combinations of coulomb force calculation
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
  #each thread work list
  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]

  #make thread core num
  p = Pool(core)

  #start thread. results is callback in array.
  callback = p.map(wrapper, worklist)
  p.close()

  #summation each thread results
  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 init_lattice():
  global xyz

  pnum = 0
  while pnum < PN:
    xyz[pnum][PX] = random.uniform(-1,1)
    xyz[pnum][PY] = random.uniform(-1,1)
    xyz[pnum][PZ] = random.uniform(-1,1)
    xyz[pnum][FX] = random.uniform(-1,1)
    xyz[pnum][FY] = random.uniform(-1,1)
    xyz[pnum][FZ] = random.uniform(-1,1)
    pnum += 1

if __name__ == "__main__":
  init_lattice()
  find_pair()

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

###########################

# multiprocessing test ####

###########################

import random

import math

import time

import scipy.special as scm

from multiprocessing import Pool

random.seed(1)

PX = 0;PY = 1;PZ = 2;

VX = 3;VY = 4;VZ = 5;

FX = 6;FY = 7;FZ = 8;

#number of particles in a line

line_num = 15

#total particle num

PN = line_num * line_num * line_num

#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)]

#Number of combinations of coulomb force calculation

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

#each thread work list

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]

#make thread core num

p = Pool(core)

#start thread. results is callback in array.

callback = p.map(wrapper, worklist)

p.close()

#summation each thread results

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 init_lattice():

global xyz

pnum = 0

while pnum < PN:

xyz[pnum][PX] = random.uniform(-1,1)

xyz[pnum][PY] = random.uniform(-1,1)

xyz[pnum][PZ] = random.uniform(-1,1)

xyz[pnum][FX] = random.uniform(-1,1)

xyz[pnum][FY] = random.uniform(-1,1)

xyz[pnum][FZ] = random.uniform(-1,1)

pnum += 1

if __name__ == "__main__":

init_lattice()

find_pair()

11月 15, 2022

The process of brute force calculation by Coulomb force

The process of brute force calculation by Coulomb force is calculated as a matrix using numpy. In the calculation, each element is calculated in a small matrix and summed up at the end to enable parallel processing later.

import random
import math
import time
import numpy as np

np.set_printoptions(threshold=np.inf)
random.seed(1)

XX = 0
YY = 1
ZZ = 2

#total particle num
PN = 100
#---------------------------------------------------
#
# parameter : divide_N
#
# The number of particles per side of a small square 
# divided for the splitting process when creating a 
# large PN^2 region with all particles PN.
#
# The smaller the setting, the more it will be divided
# into many smaller jobs.
#---------------------------------------------------
divide_N = 3

xyzV = [[0 for i in range(3)] for j in range(PN)]
xyzP = np.zeros((PN,3))
xyzF = np.zeros((PN,3))

def find_pair():
  global PN
  global xyzF
  global xyzP

  remainder = PN %  divide_N
  quotient  = PN // divide_N

  if remainder != 0:                 #ex PN=11, divide_N=4, remainder=3, quotient=2
      quotient = quotient + 1        #   2 -> 3
  extra = divide_N - remainder       #   1

  X_st = 0
  X_ed = 0
  Y_st = 0
  Y_ed = 0

  sum_f = np.zeros((PN,PN,3))

  for Xn in range(0,quotient):
    for Yn in range(0,Xn+1):

      #----------------------------------------
      local_force = calc_XY_start_end(Xn, Yn, quotient, remainder, divide_N, extra)
      #----------------------------------------

      if Xn != Yn:
        sum_f = sum_f + local_force
        sum_f = sum_f - local_force.transpose(1,0,2)
      else:
        sum_f = sum_f + local_force

  print (sum_f)

  print("=====================================")
  mmm = sum_f.sum(axis=1)
  xyzF = xyzF + mmm


def calc_XY_start_end(Xn, Yn, quotient, remainder, divide_N, extra):
  global PN
  global xyzP

  if Xn != quotient-1 and remainder != 0 or remainder == 0:
    X_st = Xn*divide_N
    X_ed = (Xn+1)*divide_N-1
  else:
    X_st = Xn*divide_N - (divide_N - remainder)
    X_ed = (Xn+1)*divide_N  - (divide_N - remainder) -1
  if Yn != quotient-1 and remainder != 0 or remainder == 0:
    Y_st = Yn*divide_N
    Y_ed = (Yn+1)*divide_N-1
  else:
    Y_st = Yn*divide_N - (divide_N - remainder)
    Y_ed = (Yn+1)*divide_N  -(divide_N - remainder) -1

  Xa = np.arange(X_st,X_ed+1)
  Ya = np.arange(Y_st,Y_ed+1)
  mx1,my1 = np.meshgrid(Xa,Ya)
  local_xyz  = xyzP[mx1]-xyzP[my1]
  distance_a = np.linalg.norm(local_xyz, axis=2)

  tmp_d = distance_a[:,:,np.newaxis]
  tmp_f = local_xyz/tmp_d
  tmp_f = np.nan_to_num(tmp_f,0)

  #extra clear
  if remainder != 0:
    if Xn + 1 == quotient:
      for dd in range(0,extra):
        tmp_f[:,dd] = 0
    if Yn + 1 == quotient:
      for dd in range(0,extra):
        tmp_f[[dd]] = 0

  if X_st != 0:
    r2 = np.zeros((X_st,divide_N,3))
    tmp_f        = np.insert(tmp_f     ,0,r2,axis=1)
  if X_ed+1 != PN:
    r2 = np.zeros((PN - X_ed - 1,divide_N,3))
    tmp_f      = np.insert(tmp_f     ,X_ed + 1,r2,axis=1)

  if Y_st != 0:
    r2 = np.zeros((Y_st,PN,3))
    tmp_f      = np.insert(tmp_f     ,0,r2,axis=0)
  if Y_ed+1 != PN:
    r2 = np.zeros((PN - Y_ed - 1,PN,3))
    tmp_f      = np.insert(tmp_f     ,Y_ed + 1,r2,axis=0)

  #print (tmp_f)
  #print("===================================================")
  return tmp_f



def init_lattice():
  global xyzF
  global xyzP

  pnum = 0
  while pnum < PN:
    xyzP[pnum][XX] = random.uniform(-1,1)
    xyzP[pnum][YY] = random.uniform(-1,1)
    xyzP[pnum][ZZ] = random.uniform(-1,1)
    xyzF[pnum][XX] = random.uniform(-1,1)
    xyzF[pnum][YY] = random.uniform(-1,1)
    xyzF[pnum][ZZ] = random.uniform(-1,1)
    pnum += 1


def results_sum():
  global xyzF
  pnum = 0
  total_F = 0
  while pnum < PN:
    total_F = total_F + xyzF[pnum][XX]
    total_F = total_F + xyzF[pnum][YY]
    total_F = total_F + xyzF[pnum][ZZ]
    pnum += 1

  print (total_F)


if __name__ == "__main__":

  init_lattice()
  find_pair()
  results_sum()

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

144

145

146

147

148

149

150

151

152

153

154

155

156

157

158

159

import random

import math

import time

import numpy as np

np.set_printoptions(threshold=np.inf)

random.seed(1)

XX = 0

YY = 1

ZZ = 2

#total particle num

PN = 100

#---------------------------------------------------

# parameter : divide_N

# The number of particles per side of a small square

# divided for the splitting process when creating a

# large PN^2 region with all particles PN.

# The smaller the setting, the more it will be divided

# into many smaller jobs.

#---------------------------------------------------

divide_N = 3

xyzV = [[0 for i in range(3)] for j in range(PN)]

xyzP = np.zeros((PN,3))

xyzF = np.zeros((PN,3))

def find_pair():

global PN

global xyzF

global xyzP

remainder = PN % divide_N

quotient = PN // divide_N

if remainder != 0: #ex PN=11, divide_N=4, remainder=3, quotient=2

quotient = quotient + 1 # 2 -> 3

extra = divide_N - remainder # 1

X_st = 0

X_ed = 0

Y_st = 0

Y_ed = 0

sum_f = np.zeros((PN,PN,3))

for Xn in range(0,quotient):

for Yn in range(0,Xn+1):

#----------------------------------------

local_force = calc_XY_start_end(Xn, Yn, quotient, remainder, divide_N, extra)

#----------------------------------------

if Xn != Yn:

sum_f = sum_f + local_force

sum_f = sum_f - local_force.transpose(1,0,2)

else:

sum_f = sum_f + local_force

print (sum_f)

print("=====================================")

mmm = sum_f.sum(axis=1)

xyzF = xyzF + mmm

def calc_XY_start_end(Xn, Yn, quotient, remainder, divide_N, extra):

global PN

global xyzP

if Xn != quotient-1 and remainder != 0 or remainder == 0:

X_st = Xn*divide_N

X_ed = (Xn+1)*divide_N-1

else:

X_st = Xn*divide_N - (divide_N - remainder)

X_ed = (Xn+1)*divide_N - (divide_N - remainder) -1

if Yn != quotient-1 and remainder != 0 or remainder == 0:

Y_st = Yn*divide_N

Y_ed = (Yn+1)*divide_N-1

else:

Y_st = Yn*divide_N - (divide_N - remainder)

Y_ed = (Yn+1)*divide_N -(divide_N - remainder) -1

Xa = np.arange(X_st,X_ed+1)

Ya = np.arange(Y_st,Y_ed+1)

mx1,my1 = np.meshgrid(Xa,Ya)

local_xyz = xyzP[mx1]-xyzP[my1]

distance_a = np.linalg.norm(local_xyz, axis=2)

tmp_d = distance_a[:,:,np.newaxis]

tmp_f = local_xyz/tmp_d

tmp_f = np.nan_to_num(tmp_f,0)

#extra clear

if remainder != 0:

if Xn + 1 == quotient:

for dd in range(0,extra):

tmp_f[:,dd] = 0

if Yn + 1 == quotient:

for dd in range(0,extra):

tmp_f[[dd]] = 0

if X_st != 0:

r2 = np.zeros((X_st,divide_N,3))

tmp_f = np.insert(tmp_f ,0,r2,axis=1)

if X_ed+1 != PN:

r2 = np.zeros((PN - X_ed - 1,divide_N,3))

tmp_f = np.insert(tmp_f ,X_ed + 1,r2,axis=1)

if Y_st != 0:

r2 = np.zeros((Y_st,PN,3))

tmp_f = np.insert(tmp_f ,0,r2,axis=0)

if Y_ed+1 != PN:

r2 = np.zeros((PN - Y_ed - 1,PN,3))

tmp_f = np.insert(tmp_f ,Y_ed + 1,r2,axis=0)

#print (tmp_f)

#print("===================================================")

return tmp_f

def init_lattice():

global xyzF

global xyzP

pnum = 0

while pnum < PN:

xyzP[pnum][XX] = random.uniform(-1,1)

xyzP[pnum][YY] = random.uniform(-1,1)

xyzP[pnum][ZZ] = random.uniform(-1,1)

xyzF[pnum][XX] = random.uniform(-1,1)

xyzF[pnum][YY] = random.uniform(-1,1)

xyzF[pnum][ZZ] = random.uniform(-1,1)

pnum += 1

def results_sum():

global xyzF

pnum = 0

total_F = 0

while pnum < PN:

total_F = total_F + xyzF[pnum][XX]

total_F = total_F + xyzF[pnum][YY]

total_F = total_F + xyzF[pnum][ZZ]

pnum += 1

print (total_F)

if __name__ == "__main__":

init_lattice()

find_pair()

results_sum()