|
ExaFMM 1
Fast-multipole Method for exascale systems
|
wrapper/coulomb.py
Go to the documentation of this file.
00001 from numpy import * 00002 import ctypes 00003 solver=ctypes.CDLL("libcoulomb.so") 00004 00005 n = 100 00006 x = random.rand(3*n) # coordinates 00007 q = ones(n) / n # charges 00008 p = zeros(n) # potential 00009 f = zeros(3*n) # force 00010 solver.FMMcalccoulomb(n, x.ctypes.data, q.ctypes.data, p.ctypes.data, f.ctypes.data, 0) # FMM coulomb solver 00011 00012 # compare with direct summation in python 00013 diffp = normp = difff = normf = 0 00014 for i in range(n): 00015 pd = fx = fy = fz = 0 00016 for j in range(n): 00017 dx = x[3*i+0] - x[3*j+0] 00018 dy = x[3*i+1] - x[3*j+1] 00019 dz = x[3*i+2] - x[3*j+2] 00020 R2 = dx * dx + dy * dy + dz * dz 00021 if R2 == 0: 00022 invR = 0 00023 else: 00024 invR = 1 / sqrt(R2) 00025 invR3 = q[j] * invR * invR * invR 00026 pd += q[j] * invR 00027 fx += dx * invR3 00028 fy += dy * invR3 00029 fz += dz * invR3 00030 diffp = (p[i] - pd) * (p[i] - pd) 00031 normp = pd * pd 00032 difff = (f[3*i+0] - fx) * (f[3*i+0] - fx) \ 00033 + (f[3*i+1] - fy) * (f[3*i+1] - fy) \ 00034 + (f[3*i+2] - fz) * (f[3*i+2] - fz) 00035 normf = fx * fx + fy * fy + fz * fz 00036 print ['potential error : ',sqrt(diffp/normp)] 00037 print ['force error : ',sqrt(difff/normf)]
