diff options
| -rw-r--r-- | detectors/__init__.py | 10 | ||||
| -rw-r--r-- | detectors/lbne.py | 10 | ||||
| -rw-r--r-- | geometry.py | 21 | ||||
| -rw-r--r-- | solids/r7081.py | 19 | ||||
| -rw-r--r-- | src/kernel.cu | 29 | ||||
| -rw-r--r-- | threadtest.py | 26 | ||||
| -rw-r--r-- | track.py | 14 | ||||
| -rwxr-xr-x | view.py | 22 |
8 files changed, 92 insertions, 59 deletions
diff --git a/detectors/__init__.py b/detectors/__init__.py index 8bcf69b..7a36101 100644 --- a/detectors/__init__.py +++ b/detectors/__init__.py @@ -1 +1,11 @@ from lbne import LBNE + +# from lbne document 94 +radius = 25.0 +height = 50.0 +nstrings = 324 +pmts_per_string = 102 +endcap_spacing = .485 + +lbne = LBNE(radius, height, nstrings, pmts_per_string, endcap_spacing) +minilbne = LBNE(radius/10, height/10, nstrings//10, pmts_per_string//10, endcap_spacing) diff --git a/detectors/lbne.py b/detectors/lbne.py index 7f1bb41..04751e6 100644 --- a/detectors/lbne.py +++ b/detectors/lbne.py @@ -13,16 +13,8 @@ from transform import rotate, make_rotation_matrix from itertools import product import make -endcap_spacing = .485 - -radius = 25.0/10.0 -height = 50.0/10.0 - -nstrings = 324//10 -pmts_per_string = 102//10 - class LBNE(Geometry): - def __init__(self): + def __init__(self, radius, height, nstrings, pmts_per_string, endcap_spacing): super(LBNE, self).__init__() # outer cylinder diff --git a/geometry.py b/geometry.py index be31e4b..04be460 100644 --- a/geometry.py +++ b/geometry.py @@ -336,6 +336,9 @@ class Geometry(object): set_material = module.get_function('set_material') for i, material in enumerate(self.materials): if material is None: + if color is False: + print 'warning: some triangles have no associated material.' + continue refractive_index = np.interp(standard_wavelengths, material.refractive_index[:,0], material.refractive_index[:,1]).astype(np.float32) @@ -363,22 +366,6 @@ class Geometry(object): set_surface(np.int32(i), surface.absorption_gpu, surface.reflection_diffuse_gpu, surface.reflection_specular_gpu, block=(1,1,1), grid=(1,1)) - self.material1_index_tex = module.get_texref('material1_lookup') - self.material2_index_tex = module.get_texref('material2_lookup') - self.surface_index_tex = module.get_texref('surface_lookup') - - self.material1_index_gpu = cuda.to_device(self.material1_index) - self.material2_index_gpu = cuda.to_device(self.material2_index) - self.surface_index_gpu = cuda.to_device(self.surface_index) - - self.material1_index_tex.set_address(self.material1_index_gpu, self.material1_index.nbytes) - self.material2_index_tex.set_address(self.material2_index_gpu, self.material2_index.nbytes) - self.surface_index_tex.set_address(self.surface_index_gpu, self.surface_index.nbytes) - - self.material1_index_tex.set_format(cuda.array_format.SIGNED_INT32, 1) - self.material2_index_tex.set_format(cuda.array_format.SIGNED_INT32, 1) - self.surface_index_tex.set_format(cuda.array_format.SIGNED_INT32, 1) - vertices = np.empty(len(self.mesh.vertices), dtype=gpuarray.vec.float4) vertices['x'] = self.mesh.vertices[:,0] vertices['y'] = self.mesh.vertices[:,1] @@ -393,7 +380,7 @@ class Geometry(object): if color: triangles['w'] = self.colors else: - triangles['w'] = (self.material1_index << 24) | (self.material2_index << 16) | (self.surface_index << 8) + triangles['w'] = ((self.material1_index & 0xff) << 24) | ((self.material2_index & 0xff) << 16) | ((self.surface_index & 0xff) << 8) lower_bounds = np.empty(self.lower_bounds.shape[0], dtype=gpuarray.vec.float4) lower_bounds['x'] = self.lower_bounds[:,0] diff --git a/solids/r7081.py b/solids/r7081.py index af9c2f0..ed677d8 100644 --- a/solids/r7081.py +++ b/solids/r7081.py @@ -27,3 +27,22 @@ r7081_inner_solid = Solid(r7081_inner_mesh, vacuum, glass, inner_surface, 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/src/kernel.cu b/src/kernel.cu index 6c3ef1b..405f06c 100644 --- a/src/kernel.cu +++ b/src/kernel.cu @@ -25,11 +25,6 @@ enum texture<float4, 1, cudaReadModeElementType> vertices; __device__ uint4 *triangles; -/* material/surface index lookup for each triangle */ -texture<int, 1, cudaReadModeElementType> material1_lookup; -texture<int, 1, cudaReadModeElementType> material2_lookup; -texture<int, 1, cudaReadModeElementType> surface_lookup; - /* lower/upper bounds for the bounding box associated with each node/leaf */ texture<float4, 1, cudaReadModeElementType> upper_bounds; texture<float4, 1, cudaReadModeElementType> lower_bounds; @@ -43,6 +38,14 @@ __device__ float3 make_float3(const float4 &a) return make_float3(a.x, a.y, a.z); } +__device__ int convert(int c) +{ + if (c & 0x80) + return (0xFFFFFF00 | c); + else + return c; +} + /* Test the intersection between a ray starting at `origin` traveling in the direction `direction` and the bounding box around node `i`. If the ray intersects the bounding box return true, else return false. */ @@ -204,7 +207,7 @@ __global__ void ray_trace(int nthreads, float3 *positions, float3 *directions, i if (triangle_index == -1) { - pixels[id] = 0; + pixels[id] = 0x000000; } else { @@ -263,9 +266,9 @@ __global__ void propagate(int nthreads, float3 *positions, float3 *directions, f float3 v1 = make_float3(tex1Dfetch(vertices, triangle_data.y)); float3 v2 = make_float3(tex1Dfetch(vertices, triangle_data.z)); - int material_in_index = tex1Dfetch(material1_lookup, last_hit_triangle); - int material_out_index = tex1Dfetch(material2_lookup, last_hit_triangle); - int surface_index = tex1Dfetch(surface_lookup, last_hit_triangle); + int material_in_index = convert(0xFF & (triangle_data.w >> 24)); + int material_out_index = convert(0xFF & (triangle_data.w >> 16)); + int surface_index = convert(0xFF & (triangle_data.w >> 8)); float3 surface_normal = cross(v1-v0,v2-v1); surface_normal /= norm(surface_normal); @@ -408,9 +411,13 @@ __global__ void propagate(int nthreads, float3 *positions, float3 *directions, f float reflection_probability = reflection_coefficient*reflection_coefficient; if ((curand_uniform(&rng) < reflection_probability) || isnan(refracted_angle)) + { direction = rotate(surface_normal, incident_angle, p); + } else + { direction = rotate(surface_normal, PI-refracted_angle, p); + } polarization = p; @@ -422,9 +429,13 @@ __global__ void propagate(int nthreads, float3 *positions, float3 *directions, f float reflection_probability = reflection_coefficient*reflection_coefficient; if ((curand_uniform(&rng) < reflection_probability) || isnan(refracted_angle)) + { direction = rotate(surface_normal, incident_angle, p); + } else + { direction = rotate(surface_normal, PI-refracted_angle, p); + } polarization = cross(p, direction); diff --git a/threadtest.py b/threadtest.py index 91a1857..d7b9d7c 100644 --- a/threadtest.py +++ b/threadtest.py @@ -1,6 +1,5 @@ from gputhread import * from Queue import Queue -from detectors import LBNE from sample import uniform_sphere import numpy as np from pycuda import gpuarray @@ -11,7 +10,7 @@ from histogram import Histogram jobs = Queue() output = Queue() -def create_job(position=(0,0,0), nphotons=1000, max_steps=10): +def create_job(position=(0,0,0), nphotons=5000, max_steps=10): positions = np.tile(position, nphotons).reshape((nphotons, 3)) directions = uniform_sphere(nphotons) polarizations = uniform_sphere(nphotons) @@ -23,14 +22,14 @@ def create_job(position=(0,0,0), nphotons=1000, max_steps=10): return Job(positions, directions, wavelengths, polarizations, times, states, last_hit_triangles, max_steps) -def generate_event(position=(0,0,0), nphotons=1000): +def generate_event(detector, position=(0,0,0), nphotons=5000): jobs.put(create_job(position, nphotons)) jobs.join() job = output.get() last_hit_triangles = job.last_hit_triangles - solids = lbne.solid_id[last_hit_triangles] + solids = detector.solid_id[last_hit_triangles] solids[last_hit_triangles == -1] = -1 surface_absorbed = job.states == 2 @@ -38,7 +37,7 @@ def generate_event(position=(0,0,0), nphotons=1000): print 'state %2i %i' % (i, len(job.states[job.states == i])) event_times = [] - for i in lbne.pmtids: + for i in detector.pmtids: photons = np.logical_and(solids == i, surface_absorbed) hit_times = job.times[photons] @@ -50,7 +49,7 @@ def generate_event(position=(0,0,0), nphotons=1000): return event_times -def likelihood(event_times, position=(0,0,0), nphotons=1000, neval=100): +def likelihood(detector, event_times, position=(0,0,0), nphotons=5000, neval=100): for i in range(neval): jobs.put(create_job(position, nphotons)) jobs.join() @@ -60,11 +59,11 @@ def likelihood(event_times, position=(0,0,0), nphotons=1000, neval=100): job = output.get() last_hit_triangles = job.last_hit_triangles - solids = lbne.solid_id[job.last_hit_triangles] + solids = detector.solid_id[job.last_hit_triangles] solids[last_hit_triangles == -1] = -1 surface_absorbed = job.states == 2 - for j in lbne.pmtids: + for j in detector.pmtids: pmt_photons = solids == j photons = np.logical_and(pmt_photons, surface_absorbed) @@ -108,22 +107,23 @@ if __name__ == '__main__': parser.add_option('-n', type='int', dest='nblocks', default=64) options, args = parser.parse_args() - lbne = LBNE() - lbne.build(bits=options.nbits) + from detectors import minilbne + + minilbne.build(bits=options.nbits) cuda.init() gputhreads = [] for i in range(options.ndevices): - gputhreads.append(GPUThread(i, lbne, jobs, output, options.nblocks)) + gputhreads.append(GPUThread(i, minilbne, jobs, output, options.nblocks)) gputhreads[-1].start() try: - event_times = generate_event() + event_times = generate_event(minilbne) for z in np.linspace(-1.0, 1.0, 100): t0 = time.time() - log_likelihood = likelihood(event_times, (z,0,0)) + log_likelihood = likelihood(minilbne, event_times, (z,0,0)) elapsed = time.time() - t0 print 'z = %5.2f, likelihood = %s, elapsed %.2f sec' % (z, log_likelihood, elapsed) finally: @@ -1,7 +1,7 @@ import numpy as np import pycuda.driver as cuda from gputhread import GPUThread -from detectors import LBNE +from detectors import minilbne from Queue import Queue from threadtest import create_job import matplotlib.pyplot as plt @@ -10,20 +10,18 @@ from color import map_wavelength from solids import r7081 from geometry import Geometry -nphotons = 100000 +nphotons = 1000 jobs = Queue() output = Queue() -#geometry = LBNE() -geometry = Geometry() -geometry.add_solid(r7081, displacement=(0,-1,0)) +geometry = minilbne geometry.build(bits=8) cuda.init() try: - gputhread = GPUThread(0, geometry, jobs, output, 64) + gputhread = GPUThread(5, geometry, jobs, output, 64) gputhread.start() job = create_job((0,0,0), nphotons) @@ -46,7 +44,7 @@ finally: gputhread.stop() gputhread.join() -mask = job.states != 0 +mask = job.states == 2 rgb = (map_wavelength(job.wavelengths[mask])*255).astype(np.uint32) @@ -56,5 +54,5 @@ def format_hex_string(s): colors = map(format_hex_string, map(hex, rgb[:,0] << 16 | rgb[:,1] << 8 | rgb[:,2])) plt.figure() -plt.plot(*roundrobin(x[mask,:,0], x[mask,:,1], colors)) +plt.plot(*roundrobin(x[mask,:,0], x[mask,:,2], colors)) plt.show() @@ -1,4 +1,5 @@ #!/usr/bin/env python +import os import numpy as np import pygame @@ -171,6 +172,22 @@ def view(geometry, name=''): done = True break + if event.key == K_F12: + if name == '': + root, ext = 'screenshot', 'png' + else: + root, ext = name, 'png' + + filename = '.'.join([root, ext]) + + i = 1 + while os.path.exists(filename): + filename = '.'.join([root + str(i), ext]) + i += 1 + + pygame.image.save(screen, filename) + print 'image saved to %s' % filename + if event.type == KEYUP: if event.key == K_LSHIFT or event.key == K_RSHIFT: shift = False @@ -207,9 +224,8 @@ if __name__ == '__main__': members = dict(inspect.getmembers(detectors) + inspect.getmembers(solids)) if args[0] in members: - if type(members[args[0]]) is type and \ - issubclass(members[args[0]], Geometry): - geometry = members[args[0]]() + if isinstance(members[args[0]], Geometry): + geometry = members[args[0]] geometry.build(options.bits) view(geometry, args[0]) elif isinstance(members[args[0]], Solid): |