Hi Julia-people 
I am working with heat transfer, more specifically simulating thermal radiation heat transfer in participating media through Monte Carlo ray tracing. I am currently working with a 2D algorithm. The reason I am looking for inspiration and inputs is that I would like to port my code to be GPU code, and before I do that I wanted to ensure that my approach is optimal. I designed the following flow chart algorithm which should be general enough to work for both 1D, 2D and 3D.
By now I have designed two code implementations of it. The code is quite involved, so I will post only pseudo code here.
My first and original take on the code was a massive use of the @goto macro. I got this implementation to be very fast, ray tracing about 1 billion rays in five minutes on 16 threads in parallel. The code has been validated against analytical results. The pseudo code is as follows:
@label emitRay
emit ray() # emit ray
enclosure = whichEnclosure() # find out which enclosure we're in
walls = localWalls() # find the local walls
dist = distToSurface() # find the distance to the wall along direction
@label doesRayHitWall # label for use with @goto
if S < dist # if the distance to the wall is larger than the sampled distance
@scatter_or_absorb # then the ray is scattered or absorbed
point = updatePoint() # update the current point
if R_omega > omega # if a random number is larger than the scattering albedo
# then the ray is absorbed
enclosure = whichEnclosure()
absorb(enclosure) += 1
@goto emitRay
else # ray is scattered
scattering = sampledirection() # sample scattering direction
end
enclosure = whichEnclosure() # find out which enclosure we're in
walls = localWalls() # find the local walls
dist, wall_index = distToSurface() # find the distance to the wall along direction
if dist > S
@goto scatter_or_absorb # ray is scattered or absorbed
else
@goto rayHitsWall # go to label
end
else
@label rayHitsWall # ray hits a wall
solidwalls = solidwalls()
if wall_index == solidwalls # ray hits a solid wall
if R_alpha > alphaWall # then the ray is reflected by the wall
newDirection = sampleDirection()
else # the ray is absorbed by the wall
absorb(wall) += 1
end
@goto emitRay
else # ray hits an imaginary boundary
point = updatePoint()
enclosure = whichEnclosure()
walls = localWalls()
dist, index = distToSurface()
@goto doesRayHitWall
end
end
Even though I got the above implementation to execute very efficiently I recently learned that using @goto statements can be confusing and is an ‘oldschool’ way of coding. Therefore I also implemented the same flow chart in a recursion type code as described below:
function rayHitWall()
enclosure = whichEnclosure()
walls = localWalls()
dist, index = distToSurface()
if dist < S # if the remaining distance is shorted than the distance to the wall
solidwalls = solidWalls()
while index != solidwall && dist < S
point = updatePoint()
enclosure = whichEnclosure()
walls = localWalls()
dist, index = distToSurface()
solidwalls = solidWalls()
end
if dist > S # ray is absorbed or scattered
if R_omega > omega # ray is absorbed
enclosure = enclosure()
absorb(enclosure) += 1
return
else # ray is scattered
newdirection = sampleDirection()
rayHitWall() # call the function itself recursively
end
elseif R_alpha < alpha # ray is absorbed by the wall
enclosure = enclosure()
walls = localWalls()
dist, index = distToSurface()
solidwalls = solidWalls()
absorb(enclosure) += 1
return
else # ray is reflected
newdirection = sampleDirection()
rayHitWall() # call the function itself recursively
end
else
if R_omega > omega # ray is absorbed
enclosure = enclosure()
absorb(enclosure) += 1
return
else # ray is scattered
newdirection = sampleDirection()
rayHitWall() # call the function itself recursively
end
end
end
for i = 1:N_rays
emission = sampleEmission()
rayHitWall()
end
The above recursive code also works, however, it is significantly slower than the implementation based on @goto statements. Furthermore, there is a risk of getting ‘stackoverflow’ with this code.
But I was curious if someone else could come up with a more modern approach to implement the flowchart in (pseudo)code which would also be efficient, and possibly also applicable to GPU programming?
Any help would be appreciated.