A faulty optical process can make the position or direction of the

photon into NaN on the GPU.  Now we abort these photons rather than
let them lock up the intersect_mesh() method.  There is a new history
bit (#31) that indicates when a NAN_ABORT has occurred, and this bit
is checked for by GPU.propagate().  If set for any of the photons, a
warning message is printed.

While not as good as preventing the NaN problem in the first place,
this at least ensures we are aware of the problem.

3 files changed, 13 insertions, 1 deletions

diff --git a/gpu.py b/gpu.py
index c9e6269..5ad195d 100644
--- a/gpu.py
+++ b/gpu.py
@@ -12,6 +12,7 @@ from tools import timeit
 import src
 from geometry import standard_wavelengths
 from color import map_to_color
+import sys
 
 cuda.init()
 
@@ -308,6 +309,9 @@ class GPU(object):
                 output_queue_gpu = temp
                 output_queue_gpu[:1].set(np.uint32(1))
                 nphotons = input_queue_gpu[:1].get()[0]
+        
+        if gpuarray.max(self.histories_gpu).get() & (1 << 31):
+            print >>sys.stderr, "WARNING: ABORTED PHOTONS IN THIS EVENT"
 
         if 'profile' in __builtins__:
             self.context.synchronize()
diff --git a/src/kernel.cu b/src/kernel.cu
index cc567b8..307412e 100644
--- a/src/kernel.cu
+++ b/src/kernel.cu
@@ -345,6 +345,13 @@ __global__ void propagate(int first_photon, int nthreads,
 
 		int command;
 
+		// check for NaN and fail
+		if (isnan(p.direction.x*p.direction.y*p.direction.z*p.position.x*p.position.y*p.position.z)) 
+		{
+			p.history |= NO_HIT | NAN_ABORT;
+			break;
+		}
+
 		fill_state(s, p);
 
 		if (p.last_hit_triangle == -1)
diff --git a/src/photon.h b/src/photon.h
index f471866..4203a50 100644
--- a/src/photon.h
+++ b/src/photon.h
@@ -44,7 +44,8 @@ enum
 	SURFACE_ABSORB   = 0x1 << 3,
 	RAYLEIGH_SCATTER = 0x1 << 4,
 	REFLECT_DIFFUSE  = 0x1 << 5,
-	REFLECT_SPECULAR = 0x1 << 6
+	REFLECT_SPECULAR = 0x1 << 6,
+	NAN_ABORT        = 0x1 << 31
 }; // processes
 
 enum {BREAK, CONTINUE, PASS}; // return value from propagate_to_boundary