argument '-j' to threadtest.py now specifies a list of device ids to be used. GPUThread objects now shallow copy the geometry so that threads are not writing to the same memory when the geometry is loaded onto the gpu. the model number for the 12" Hamamatsu PMT is R11708, not r7081 (which is the model for the 10" PMT); all references to the 12" PMT were changed accordingly. only allocate space for 20 materials and 20 surfaces on the gpu instead of 100 to save some space. started to modify track.py to build its own photons and module since the GPUThread object only copies photon hit times back from the gpu (not track information), but I am waiting to find out if pycuda GPUArrays can be used with vector types.

12 files changed, 153 insertions, 133 deletions

diff --git a/detectors/lbne.py b/detectors/lbne.py
index 9f74975..2833ebe 100644
--- a/detectors/lbne.py
+++ b/detectors/lbne.py
@@ -19,9 +19,9 @@ class LBNE(Geometry):
         super(LBNE, self).__init__()
 
         if cut_pmt:
-            pmt_mesh = r7081_cut
+            pmt_mesh = r11708_cut
         else:
-            pmt_mesh = r7081
+            pmt_mesh = r11708
 
         # outer cylinder
         cylinder_mesh = \
diff --git a/geometry.py b/geometry.py
index 2405e07..f0f5c20 100644
--- a/geometry.py
+++ b/geometry.py
@@ -165,8 +165,8 @@ def morton_order(mesh, bits):
                             np.max(mesh[:,:,1]),
                             np.max(mesh[:,:,2])])
 
-    if bits <= 0 or bits > 12:
-        raise Exception('number of bits must be in the range (0,12].')
+    if bits <= 0 or bits > 10:
+        raise Exception('number of bits must be in the range (0,10].')
 
     max_value = 2**bits - 1
 
@@ -403,31 +403,43 @@ class Geometry(object):
         self.node_map_gpu = cuda.to_device(self.node_map)
         self.node_length_gpu = cuda.to_device(self.node_length)
 
-        print 'Device usage:'
-        print 'vertices:', vertices.nbytes
-        print 'triangles:', triangles.nbytes
-        print 'lower_bounds:', lower_bounds.nbytes
-        print 'upper_bounds:', upper_bounds.nbytes
-        print 'node_map:', self.node_map.nbytes
-        print 'node_length:', self.node_length.nbytes
+        def format_size(size):
+            if size < 1e3:
+                return '%.1f%s' % (size, ' ')
+            elif size < 1e6:
+                return '%.1f%s' % (size/1e3, 'K')
+            elif size < 1e9:
+                return '%.1f%s' % (size/1e6, 'M')
+            else:
+                return '%.1f%s' % (size/1e9, 'G')
+
+        def format_array(name, array):
+            return '%-15s %6s %6s' % (name, format_size(len(array)), format_size(array.nbytes))
+
+        print 'Device Usage:'
+        print format_array('vertices', vertices)
+        print format_array('triangles', triangles)
+        print format_array('lower_bounds', self.lower_bounds)
+        print format_array('upper_bounds', self.upper_bounds)
+        print format_array('node_map', self.node_map)
+        print format_array('node_length', self.node_length)
+        print '%-15s %6s %6s' % ('total', '', format_size(vertices.nbytes + triangles.nbytes + self.lower_bounds.nbytes + self.upper_bounds.nbytes + self.node_map.nbytes + self.node_length.nbytes))
 
         set_pointer = module.get_function('set_pointer')
         set_pointer(self.triangles_gpu, self.vertices_gpu, 
                     block=(1,1,1), grid=(1,1))
 
-        lower_bounds_tex = module.get_texref('lower_bounds')
-        upper_bounds_tex = module.get_texref('upper_bounds')
-        node_map_tex = module.get_texref('node_map')
-        node_length_tex = module.get_texref('node_length')
+        self.lower_bounds_tex = module.get_texref('lower_bounds')
+        self.upper_bounds_tex = module.get_texref('upper_bounds')
+        self.node_map_tex = module.get_texref('node_map')
+        self.node_length_tex = module.get_texref('node_length')
 
-        lower_bounds_tex.set_address(self.lower_bounds_gpu, lower_bounds.nbytes)
-        upper_bounds_tex.set_address(self.upper_bounds_gpu, upper_bounds.nbytes)
-        node_map_tex.set_address(self.node_map_gpu, self.node_map.nbytes)
-        node_length_tex.set_address(self.node_length_gpu, self.node_length.nbytes)
+        self.lower_bounds_tex.set_address(self.lower_bounds_gpu, lower_bounds.nbytes)
+        self.upper_bounds_tex.set_address(self.upper_bounds_gpu, upper_bounds.nbytes)
+        self.node_map_tex.set_address(self.node_map_gpu, self.node_map.nbytes)
+        self.node_length_tex.set_address(self.node_length_gpu, self.node_length.nbytes)
 
-        lower_bounds_tex.set_format(cuda.array_format.FLOAT, 4)
-        upper_bounds_tex.set_format(cuda.array_format.FLOAT, 4)
-        node_map_tex.set_format(cuda.array_format.UNSIGNED_INT32, 1)
-        node_length_tex.set_format(cuda.array_format.UNSIGNED_INT32, 1)
-
-        return [lower_bounds_tex, upper_bounds_tex, node_map_tex, node_length_tex]
+        self.lower_bounds_tex.set_format(cuda.array_format.FLOAT, 4)
+        self.upper_bounds_tex.set_format(cuda.array_format.FLOAT, 4)
+        self.node_map_tex.set_format(cuda.array_format.UNSIGNED_INT32, 1)
+        self.node_length_tex.set_format(cuda.array_format.UNSIGNED_INT32, 1)
diff --git a/gputhread.py b/gputhread.py
index ba0bb4b..06921d5 100644
--- a/gputhread.py
+++ b/gputhread.py
@@ -1,4 +1,5 @@
 import numpy as np
+from copy import copy
 import time
 import pycuda.driver as cuda
 from pycuda.characterize import sizeof
@@ -13,7 +14,7 @@ class GPUThread(threading.Thread):
         threading.Thread.__init__(self)
 
         self.device_id = device_id
-        self.geometry = geometry
+        self.geometry = copy(geometry)
         self.jobs = jobs
         self.output = output
         self.nblocks = nblocks
diff --git a/solids/__init__.py b/solids/__init__.py
index 8fc66a4..f4c5812 100644
--- a/solids/__init__.py
+++ b/solids/__init__.py
@@ -1,2 +1,2 @@
-from r7081 import r7081
-from r7081_cut import r7081_cut
+from r11708 import r11708
+from r11708_cut import r11708_cut
diff --git a/solids/r11708.py b/solids/r11708.py
new file mode 100644
index 0000000..a0476d0
--- /dev/null
+++ b/solids/r11708.py
@@ -0,0 +1,48 @@
+import os
+import sys
+import numpy as np
+
+dir = os.path.split(os.path.realpath(__file__))[0]
+sys.path.append(dir + '/..')
+
+import models
+from stl import mesh_from_stl
+from geometry import *
+from materials import *
+
+r11708_outer_mesh = mesh_from_stl(models.dir + '/hamamatsu_12inch_outer.stl')
+r11708_inner_mesh = mesh_from_stl(models.dir + '/hamamatsu_12inch_inner.stl')
+
+photocathode_triangles = np.mean(r11708_inner_mesh[:], axis=1)[:,1] > 0
+
+inner_color = np.empty(len(r11708_inner_mesh.triangles), np.uint32)
+inner_color[photocathode_triangles] = 0xff0000
+inner_color[~photocathode_triangles] = 0x00ff00
+
+inner_surface = np.empty(len(r11708_inner_mesh.triangles), np.object)
+inner_surface[photocathode_triangles] = black_surface
+inner_surface[~photocathode_triangles] = shiny_surface
+
+r11708_inner_solid = Solid(r11708_inner_mesh, vacuum, glass, inner_surface, color=inner_color)
+r11708_outer_solid = Solid(r11708_outer_mesh, glass, lightwater_sno)
+
+r11708 = r11708_inner_solid + r11708_outer_solid
+
+if __name__ == '__main__':
+    from view import view
+    from copy import deepcopy
+
+    r11708_outer_mesh_cutaway = deepcopy(r11708_outer_mesh)
+    r11708_outer_mesh_cutaway.triangles = \
+        r11708_outer_mesh_cutaway.triangles[\
+        np.mean(r11708_outer_mesh_cutaway[:], axis=1)[:,0] > 0]
+
+    r11708_outer_solid_cutaway = Solid(r11708_outer_mesh_cutaway, glass, lightwater_sno)
+
+    r11708_cutaway = r11708_inner_solid + r11708_outer_solid_cutaway
+
+    geometry = Geometry()
+    geometry.add_solid(r11708_cutaway)
+    geometry.build(bits=8)
+
+    view(geometry, 'r11708_cutaway')
diff --git a/solids/r11708_cut.py b/solids/r11708_cut.py
new file mode 100644
index 0000000..a9ffeba
--- /dev/null
+++ b/solids/r11708_cut.py
@@ -0,0 +1,48 @@
+import os
+import sys
+import numpy as np
+
+dir = os.path.split(os.path.realpath(__file__))[0]
+sys.path.append(dir + '/..')
+
+import models
+from stl import mesh_from_stl
+from geometry import *
+from materials import *
+
+r11708_outer_mesh = mesh_from_stl(models.dir + '/hamamatsu_12inch_outer_cut.stl')
+r11708_inner_mesh = mesh_from_stl(models.dir + '/hamamatsu_12inch_inner_cut.stl')
+
+photocathode_triangles = np.mean(r11708_inner_mesh[:], axis=1)[:,1] > 0
+
+inner_color = np.empty(len(r11708_inner_mesh.triangles), np.uint32)
+inner_color[photocathode_triangles] = 0xff0000
+inner_color[~photocathode_triangles] = 0x00ff00
+
+inner_surface = np.empty(len(r11708_inner_mesh.triangles), np.object)
+inner_surface[photocathode_triangles] = black_surface
+inner_surface[~photocathode_triangles] = shiny_surface
+
+r11708_inner_solid = Solid(r11708_inner_mesh, vacuum, glass, inner_surface, color=inner_color)
+r11708_outer_solid = Solid(r11708_outer_mesh, glass, lightwater_sno)
+
+r11708_cut = r11708_inner_solid + r11708_outer_solid
+
+if __name__ == '__main__':
+    from view import view
+    from copy import deepcopy
+
+    r11708_outer_mesh_cutaway = deepcopy(r11708_outer_mesh)
+    r11708_outer_mesh_cutaway.triangles = \
+        r11708_outer_mesh_cutaway.triangles[\
+        np.mean(r11708_outer_mesh_cutaway[:], axis=1)[:,0] > 0]
+
+    r11708_outer_solid_cutaway = Solid(r11708_outer_mesh_cutaway, glass, lightwater_sno)
+
+    r11708_cutaway = r11708_inner_solid + r11708_outer_solid_cutaway
+
+    geometry = Geometry()
+    geometry.add_solid(r11708_cutaway)
+    geometry.build(bits=8)
+
+    view(geometry, 'r11708_cutaway')
diff --git a/solids/r7081.py b/solids/r7081.py
deleted file mode 100644
index ed677d8..0000000
--- a/solids/r7081.py
+++ /dev/null
@@ -1,48 +0,0 @@
-import os
-import sys
-import numpy as np
-
-dir = os.path.split(os.path.realpath(__file__))[0]
-sys.path.append(dir + '/..')
-
-import models
-from stl import mesh_from_stl
-from geometry import *
-from materials import *
-
-r7081_outer_mesh = mesh_from_stl(models.dir + '/hamamatsu_12inch_outer.stl')
-r7081_inner_mesh = mesh_from_stl(models.dir + '/hamamatsu_12inch_inner.stl')
-
-photocathode_triangles = np.mean(r7081_inner_mesh[:], axis=1)[:,1] > 0
-
-inner_color = np.empty(len(r7081_inner_mesh.triangles), np.uint32)
-inner_color[photocathode_triangles] = 0xff0000
-inner_color[~photocathode_triangles] = 0x00ff00
-
-inner_surface = np.empty(len(r7081_inner_mesh.triangles), np.object)
-inner_surface[photocathode_triangles] = black_surface
-inner_surface[~photocathode_triangles] = shiny_surface
-
-r7081_inner_solid = Solid(r7081_inner_mesh, vacuum, glass, inner_surface, color=inner_color)
-r7081_outer_solid = Solid(r7081_outer_mesh, glass, lightwater_sno)
-
-r7081 = r7081_inner_solid + r7081_outer_solid
-
-if __name__ == '__main__':
-    from view import view
-    from copy import deepcopy
-
-    r7081_outer_mesh_cutaway = deepcopy(r7081_outer_mesh)
-    r7081_outer_mesh_cutaway.triangles = \
-        r7081_outer_mesh_cutaway.triangles[\
-        np.mean(r7081_outer_mesh_cutaway[:], axis=1)[:,0] > 0]
-
-    r7081_outer_solid_cutaway = Solid(r7081_outer_mesh_cutaway, glass, lightwater_sno)
-
-    r7081_cutaway = r7081_inner_solid + r7081_outer_solid_cutaway
-
-    geometry = Geometry()
-    geometry.add_solid(r7081_cutaway)
-    geometry.build(bits=8)
-
-    view(geometry, 'r7081_cutaway')
diff --git a/solids/r7081_cut.py b/solids/r7081_cut.py
deleted file mode 100644
index d38d54b..0000000
--- a/solids/r7081_cut.py
+++ /dev/null
@@ -1,48 +0,0 @@
-import os
-import sys
-import numpy as np
-
-dir = os.path.split(os.path.realpath(__file__))[0]
-sys.path.append(dir + '/..')
-
-import models
-from stl import mesh_from_stl
-from geometry import *
-from materials import *
-
-r7081_outer_mesh = mesh_from_stl(models.dir + '/hamamatsu_12inch_outer_cut.stl')
-r7081_inner_mesh = mesh_from_stl(models.dir + '/hamamatsu_12inch_inner_cut.stl')
-
-photocathode_triangles = np.mean(r7081_inner_mesh[:], axis=1)[:,1] > 0
-
-inner_color = np.empty(len(r7081_inner_mesh.triangles), np.uint32)
-inner_color[photocathode_triangles] = 0xff0000
-inner_color[~photocathode_triangles] = 0x00ff00
-
-inner_surface = np.empty(len(r7081_inner_mesh.triangles), np.object)
-inner_surface[photocathode_triangles] = black_surface
-inner_surface[~photocathode_triangles] = shiny_surface
-
-r7081_inner_solid = Solid(r7081_inner_mesh, vacuum, glass, inner_surface, color=inner_color)
-r7081_outer_solid = Solid(r7081_outer_mesh, glass, lightwater_sno)
-
-r7081_cut = r7081_inner_solid + r7081_outer_solid
-
-if __name__ == '__main__':
-    from view import view
-    from copy import deepcopy
-
-    r7081_outer_mesh_cutaway = deepcopy(r7081_outer_mesh)
-    r7081_outer_mesh_cutaway.triangles = \
-        r7081_outer_mesh_cutaway.triangles[\
-        np.mean(r7081_outer_mesh_cutaway[:], axis=1)[:,0] > 0]
-
-    r7081_outer_solid_cutaway = Solid(r7081_outer_mesh_cutaway, glass, lightwater_sno)
-
-    r7081_cutaway = r7081_inner_solid + r7081_outer_solid_cutaway
-
-    geometry = Geometry()
-    geometry.add_solid(r7081_cutaway)
-    geometry.build(bits=8)
-
-    view(geometry, 'r7081_cutaway')
diff --git a/src/materials.h b/src/materials.h
index 77c9f43..6115396 100644
--- a/src/materials.h
+++ b/src/materials.h
@@ -20,8 +20,8 @@ struct Surface
 	float *reflection_specular;
 };
 
-__device__ Material materials[100];
-__device__ Surface surfaces[100];
+__device__ Material materials[20];
+__device__ Surface surfaces[20];
 
 __device__ float interp_property(const float &x, const float *fp)
 {
diff --git a/threadtest.py b/threadtest.py
index 5dffb98..92b0c81 100644
--- a/threadtest.py
+++ b/threadtest.py
@@ -61,7 +61,7 @@ if __name__ == '__main__':
 
     parser = optparse.OptionParser('%prog')
     parser.add_option('-b', type='int', dest='nbits', default=8)
-    parser.add_option('-j', type='int', dest='ndevices', default=1)
+    parser.add_option('-j', type='string', dest='devices', default=1)
     parser.add_option('-n', type='int', dest='nblocks', default=64)
     options, args = parser.parse_args()
 
@@ -74,7 +74,7 @@ if __name__ == '__main__':
     cuda.init()
 
     gputhreads = []
-    for i in range(options.ndevices):
+    for i in [int(s) for s in options.devices.split(',')]:
         gputhreads.append(GPUThread(i, detector, jobs, output, options.nblocks))
         gputhreads[-1].start()
 
diff --git a/track.py b/track.py
index 6f53c09..101e08e 100644
--- a/track.py
+++ b/track.py
@@ -7,8 +7,9 @@ from threadtest import create_job
 import matplotlib.pyplot as plt
 from itertoolset import roundrobin
 from color import map_wavelength
-from solids import r7081
+from solids import r11708
 from geometry import Geometry
+import src
 
 nphotons = 1000
 
@@ -21,8 +22,14 @@ geometry.build(bits=8)
 cuda.init()
 
 try:
-    gputhread = GPUThread(5, geometry, jobs, output, 64)
-    gputhread.start()
+    #gputhread = GPUThread(5, geometry, jobs, output, 64)
+    #gputhread.start()
+
+    device = cuda.Device(5)
+    context = device.make_context()
+    module = SourceModule(src.kernel, options=['-I' + src.dir], no_extern_c=True, cache_dir=False)
+
+    propagate = module.get_function('propagate')
 
     job = create_job((0,0,0), nphotons)
 
diff --git a/view.py b/view.py
index 952964f..2abb76d 100755
--- a/view.py
+++ b/view.py
@@ -55,7 +55,7 @@ def view(viewable, name='', bits=8):
 
     module = SourceModule(src.kernel, options=['-I' + src.dir],
                           no_extern_c=True, cache_dir=False)
-    texrefs = geometry.load(module, color=True)
+    geometry.load(module, color=True)
     cuda_raytrace = module.get_function('ray_trace')
     cuda_rotate = module.get_function('rotate')
     cuda_translate = module.get_function('translate')
@@ -100,7 +100,7 @@ def view(viewable, name='', bits=8):
     def render():
         """Render the mesh and display to screen."""
         t0 = time.time()
-        cuda_raytrace(np.int32(pixels.size), origins_gpu, directions_gpu, np.int32(geometry.node_map.size-1), np.int32(geometry.first_node), pixels_gpu, texrefs=texrefs, **gpu_kwargs)
+        cuda_raytrace(np.int32(pixels.size), origins_gpu, directions_gpu, np.int32(geometry.node_map.size-1), np.int32(geometry.first_node), pixels_gpu, **gpu_kwargs)
         cuda.Context.synchronize()
         elapsed = time.time() - t0