Merge pull request #4 from PyCOMPLETE/release/GPUupdate

aoeftiger · web-flow · commit 39067257a816 · 2017-06-14T17:56:23.000+02:00
Release/gpuupdate
diff --git a/Bassetti_Erskine.py b/Bassetti_Erskine.py
@@ -8,7 +8,7 @@
 #     This file is part of the code:
 #                                                                      		    
 # 
-#		           PyPIC Version 2.2.0                     
+#		           PyPIC Version 2.2.1                     
 #                  
 #                                                                       
 #     Author and contact:   Giovanni IADAROLA 
diff --git a/FFT_OpenBoundary.py b/FFT_OpenBoundary.py
@@ -8,7 +8,7 @@
 #     This file is part of the code:
 #                                                                      		    
 # 
-#		           PyPIC Version 2.2.0                     
+#		           PyPIC Version 2.2.1                     
 #                  
 #                                                                       
 #     Author and contact:   Giovanni IADAROLA 
diff --git a/FFT_OpenBoundary_SquareGrid.py b/FFT_OpenBoundary_SquareGrid.py
@@ -8,7 +8,7 @@
 #     This file is part of the code:
 #                                                                      		    
 # 
-#		           PyPIC Version 2.2.0                     
+#		           PyPIC Version 2.2.1                     
 #                  
 #                                                                       
 #     Author and contact:   Giovanni IADAROLA 
diff --git a/FFT_PEC_Boundary_SquareGrid.py b/FFT_PEC_Boundary_SquareGrid.py
@@ -8,7 +8,7 @@
 #     This file is part of the code:
 #                                                                      		    
 # 
-#		           PyPIC Version 2.2.0                     
+#		           PyPIC Version 2.2.1                     
 #                  
 #                                                                       
 #     Author and contact:   Giovanni IADAROLA 
diff --git a/FiniteDifferences_ShortleyWeller_SquareGrid.py b/FiniteDifferences_ShortleyWeller_SquareGrid.py
@@ -8,7 +8,7 @@
 #     This file is part of the code:
 #                                                                      		    
 # 
-#		           PyPIC Version 2.2.0                     
+#		           PyPIC Version 2.2.1                     
 #                  
 #                                                                       
 #     Author and contact:   Giovanni IADAROLA 
diff --git a/FiniteDifferences_ShortleyWeller_SquareGrid_extrapolation.py b/FiniteDifferences_ShortleyWeller_SquareGrid_extrapolation.py
@@ -8,7 +8,7 @@
 #     This file is part of the code:
 #                                                                      		    
 # 
-#		           PyPIC Version 2.2.0                     
+#		           PyPIC Version 2.2.1                     
 #                  
 #                                                                       
 #     Author and contact:   Giovanni IADAROLA 
diff --git a/FiniteDifferences_Staircase_SquareGrid.py b/FiniteDifferences_Staircase_SquareGrid.py
@@ -8,7 +8,7 @@
 #     This file is part of the code:
 #                                                                      		    
 # 
-#		           PyPIC Version 2.2.0                     
+#		           PyPIC Version 2.2.1                     
 #                  
 #                                                                       
 #     Author and contact:   Giovanni IADAROLA 
diff --git a/GPU/__init__.py b/GPU/__init__.py
@@ -0,0 +1 @@
+from _version import __version__
diff --git a/GPU/_version.py b/GPU/_version.py
@@ -0,0 +1,16 @@
+import os, subprocess
+worktree = os.path.dirname(os.path.abspath(__file__))
+gitdir = worktree + '/.git/'
+try:
+    __version__ = subprocess.check_output(
+        'git --git-dir=' + gitdir + ' --work-tree=' +
+        worktree + ' describe --long --dirty --abbrev=10 --tags', shell=True)
+    __version__ = __version__.rstrip() # remove trailing \n
+    __version__ = __version__[1:] # remove leading v
+    # remove commit hash to conform to PEP440:
+    split_ = __version__.split('-')
+    __version__ = split_[0]
+    if split_[1] != '0':
+        __version__ += '.' + split_[1]
+except:
+    __version__ = '(no git available to determine version)'
diff --git a/GPU/itest/PyPIC GPU Example improved_p2m.ipynb b/GPU/itest/PyPIC GPU Example improved_p2m.ipynb
diff --git a/GPU/itest/PyPIC GPU Example.ipynb b/GPU/itest/PyPIC GPU Example.ipynb
diff --git a/GPU/m2p/m2p.py b/GPU/m2p/m2p.py
@@ -10,31 +10,55 @@ def mesh_to_particles_CPU_3d(mesh, mesh_quantity, indices, weights):
     ip, jp, kp = indices
     stridex = mesh.nx
     stridey = mesh.ny
-    mesh_quantity = np.ravel(mesh_quantity)
-    particles_quantity = (mesh_quantity[jp   + stridex*ip     + stridex*stridey*kp    ] * weights[0]
-                        + mesh_quantity[jp   + stridex*(ip+1) + stridex*stridey*kp    ] * weights[1]
-                        + mesh_quantity[jp+1 + stridex*ip     + stridex*stridey*kp    ] * weights[2]
-                        + mesh_quantity[jp+1 + stridex*(ip+1) + stridex*stridey*kp    ] * weights[3]
-                        + mesh_quantity[jp   + stridex*ip     + stridex*stridey*(kp+1)] * weights[4]
-                        + mesh_quantity[jp   + stridex*(ip+1) + stridex*stridey*(kp+1)] * weights[5]
-                        + mesh_quantity[jp+1 + stridex*ip     + stridex*stridey*(kp+1)] * weights[6]
-                        + mesh_quantity[jp+1 + stridex*(ip+1) + stridex*stridey*(kp+1)] * weights[7])
+    mq = np.ravel(mesh_quantity)
+
+    @np.vectorize
+    def check_outside(ip, jp, kp):
+        outside_idx = (jp < 0 or jp >= mesh.nx - 1 or
+                       ip < 0 or ip >= mesh.ny - 1 or
+                       kp < 0 or kp >= mesh.nz - 1)
+        return outside_idx
+    outside_idx = check_outside(ip, jp, kp)
+    inside_idx = ~outside_idx
+    ip, jp, kp = ip[inside_idx], jp[inside_idx], kp[inside_idx]
+    weights = map(lambda w: w[inside_idx], weights)
+
+    particles_quantity = np.empty(len(indices[0]), dtype=mesh_quantity.dtype)
+    particles_quantity[inside_idx] = (
+        mq[jp   + stridex*ip     + stridex*stridey*kp    ] * weights[0]
+      + mq[jp   + stridex*(ip+1) + stridex*stridey*kp    ] * weights[1]
+      + mq[jp+1 + stridex*ip     + stridex*stridey*kp    ] * weights[2]
+      + mq[jp+1 + stridex*(ip+1) + stridex*stridey*kp    ] * weights[3]
+      + mq[jp   + stridex*ip     + stridex*stridey*(kp+1)] * weights[4]
+      + mq[jp   + stridex*(ip+1) + stridex*stridey*(kp+1)] * weights[5]
+      + mq[jp+1 + stridex*ip     + stridex*stridey*(kp+1)] * weights[6]
+      + mq[jp+1 + stridex*(ip+1) + stridex*stridey*(kp+1)] * weights[7])
+
+    particles_quantity[outside_idx] = 0
     return particles_quantity
 
 def mesh_to_particles_CPU_2d(mesh, mesh_quantity, indices, weights):
     """CPU kernel for 3d mesh to particles quantity interpolation"""
     ip, jp = indices
+    stridex = mesh.nx
+    mesh_quantity = np.ravel(mesh_quantity)
+
+    @np.vectorize
     def check_outside(ip, jp):
-        outside_idx = ip < 0 and jp < 0 and ip >= mesh.ny and jp >= mesh.nx
+        outside_idx = (jp < 0 or jp >= mesh.nx - 1 or
+                       ip < 0 or ip >= mesh.ny - 1)
         return outside_idx
-    check_out = np.vectorize(check_outside)
-    outside_idx = check_out(ip,jp)
+    outside_idx = check_outside(ip, jp)
+    inside_idx = ~outside_idx
+    ip, jp = ip[inside_idx], jp[inside_idx]
+    weights = map(lambda w: w[inside_idx], weights)
+
+    particles_quantity = np.empty(len(indices[0]), dtype=mesh_quantity.dtype)
+    particles_quantity[inside_idx] = (
+        mesh_quantity[jp   + stridex*ip    ] * weights[0]
+      + mesh_quantity[jp   + stridex*(ip+1)] * weights[1]
+      + mesh_quantity[jp+1 + stridex*ip    ] * weights[2]
+      + mesh_quantity[jp+1 + stridex*(ip+1)] * weights[3])
 
-    mesh_quantity = np.ravel(mesh_quantity)
-    stridex = mesh.nx
-    particles_quantity = (mesh_quantity[jp   + stridex*ip    ] * weights[0]
-                        + mesh_quantity[jp   + stridex*(ip+1)] * weights[1]
-                        + mesh_quantity[jp+1 + stridex*ip    ] * weights[2]
-                        + mesh_quantity[jp+1 + stridex*(ip+1)] * weights[3])
     particles_quantity[outside_idx] = 0
     return particles_quantity
diff --git a/GPU/p2m/p2m.py b/GPU/p2m/p2m.py
@@ -16,6 +16,10 @@ def particles_to_mesh_CPU_3d(mesh, n_macroparticles, mesh_indices, weights):
         ip = mesh_indices[0][p]
         jp = mesh_indices[1][p]
         kp = mesh_indices[2][p]
+        if     (jp < 0 or jp >= mesh.nx - 1 or
+                ip < 0 or ip >= mesh.ny - 1 or
+                kp < 0 or kp >= mesh.nz - 1):
+            continue
         mesh_density [jp   + stridex*ip     + stridex*stridey*kp    ] += weights[0][p]
         mesh_density [jp   + stridex*(ip+1) + stridex*stridey*kp    ] += weights[1][p]
         mesh_density [jp+1 + stridex*ip     + stridex*stridey*kp    ] += weights[2][p]
@@ -36,10 +40,12 @@ def particles_to_mesh_CPU_2d(mesh, n_macroparticles, mesh_indices, weights):
     for p in xrange(n_macroparticles):
         ip = mesh_indices[0][p]
         jp = mesh_indices[1][p]
-        if ip >= 0 and jp >= 0 and ip < mesh.ny and jp < mesh.nx:
-            mesh_density [jp   + stridex*ip    ] += weights[0][p]
-            mesh_density [jp   + stridex*(ip+1)] += weights[1][p]
-            mesh_density [jp+1 + stridex*ip    ] += weights[2][p]
-            mesh_density [jp+1 + stridex*(ip+1)] += weights[3][p]
+        if     (jp < 0 or jp >= mesh.nx - 1 or
+                ip < 0 or ip >= mesh.ny - 1):
+            continue
+        mesh_density [jp   + stridex*ip    ] += weights[0][p]
+        mesh_density [jp   + stridex*(ip+1)] += weights[1][p]
+        mesh_density [jp+1 + stridex*ip    ] += weights[2][p]
+        mesh_density [jp+1 + stridex*(ip+1)] += weights[3][p]
     mesh_density = mesh_density.reshape(mesh.shape)
     return mesh_density
diff --git a/GPU/p2m/p2m_kernels.cu b/GPU/p2m/p2m_kernels.cu
@@ -3,14 +3,24 @@
    @author: Stefan Hegglin, Adrian Oeftiger
 */
 
+#include <cuda.h>
+
 // implementation from: http://docs.nvidia.com/cuda/cuda-c-programming-guide/#atomicadd
-// very slow, for testing purposes
-__device__ double atomicAdd(double* address, double val)
+// very slow, for <NVIDIA P100 purposes where double atomicAdd does not exist yet
+// 2017-05-22 edit based on:
+// http://stackoverflow.com/questions/39274472/error-function-atomicadddouble-double-has-already-been-defined
+
+#if !defined(__CUDA_ARCH__) || __CUDA_ARCH__ >= 600
+
+#else
+static __inline__ __device__ double atomicAdd(double* address, double val)
 {
     unsigned long long int* address_as_ull =
                               (unsigned long long int*)address;
     unsigned long long int old = *address_as_ull, assumed;
 
+    if (val==0.0)
+        return __longlong_as_double(old);
     do {
         assumed = old;
         old = atomicCAS(address_as_ull, assumed,
@@ -22,6 +32,7 @@ __device__ double atomicAdd(double* address, double val)
 
     return __longlong_as_double(old);
 }
+#endif
 
 
 extern "C" {
diff --git a/PyPIC_Scatter_Gather.py b/PyPIC_Scatter_Gather.py
@@ -8,7 +8,7 @@
 #     This file is part of the code:
 #                                                                      		    
 # 
-#		           PyPIC Version 2.2.0                     
+#		           PyPIC Version 2.2.1                     
 #                  
 #                                                                       
 #     Author and contact:   Giovanni IADAROLA 
@@ -66,7 +66,7 @@ class PyPIC_Scatter_Gather(object):
     def __init__(self, x_aper=None, y_aper=None, dx=None, dy=None, xg=None, yg=None, 
                 x_min=None, x_max=None, y_min=None, y_max=None, *args, **kwargs):
 
-        print 'PyPIC Version 2.2.0'
+        print 'PyPIC Version 2.2.1'
         
         if xg!=None and yg!=None:
             assert(x_aper==None and y_aper==None and dx==None and dy==None)
diff --git a/change_version_number.py b/change_version_number.py
@@ -9,7 +9,7 @@
         with open(filename) as fid:
             content=fid.read()
         if 'giovanni.iadarola@cern.ch' in content:
-            content=content.replace('PyPIC Version 2.2.0', 'PyPIC Version 2.2.0')
+            content=content.replace('PyPIC Version 2.2.1', 'PyPIC Version 2.2.1')
             with open(filename,'w') as fid:
                 fid.write(content)
 
diff --git a/geom_impact_poly.py b/geom_impact_poly.py
@@ -7,7 +7,7 @@
 #     
 #     This file is part of the code:
 #                                                                                          
-#                   PyPIC Version 2.2.0                   
+#                   PyPIC Version 2.2.1                   
 #                  
 #                                                                       
 #     Author and contact:   Giovanni IADAROLA 

Original file line number	Diff line number	Diff line change
`@@ -8,7 +8,7 @@`
`8`	`8`	`# This file is part of the code:`
`9`	`9`	`#`
`10`	`10`	`#`
`11`		`-# PyPIC Version 2.2.0`
	`11`	`+# PyPIC Version 2.2.1`
`12`	`12`	`#`
`13`	`13`	`#`
`14`	`14`	`# Author and contact: Giovanni IADAROLA`