1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
|
//-*-c-*-
#include "linalg.h"
#include "mesh.h"
#include "geometry.h"
#include "photon.h"
extern "C"
{
__global__ void
propagate(int first_photon, int nthreads, unsigned int *input_queue,
unsigned int *output_queue, curandState *rng_states,
float3 *positions, float3 *directions,
float *wavelengths, float3 *polarizations,
float *times, unsigned int *histories,
int *last_hit_triangles, int max_steps,
Geometry *geometry)
{
int id = blockIdx.x*blockDim.x + threadIdx.x;
if (id >= nthreads)
return;
curandState rng = rng_states[id];
int photon_id = input_queue[first_photon + id];
Photon p;
p.position = positions[photon_id];
p.direction = directions[photon_id];
p.direction /= norm(p.direction);
p.polarization = polarizations[photon_id];
p.polarization /= norm(p.polarization);
p.wavelength = wavelengths[photon_id];
p.time = times[photon_id];
p.last_hit_triangle = last_hit_triangles[photon_id];
p.history = histories[photon_id];
if (p.history & (NO_HIT | BULK_ABSORB | SURFACE_DETECT | SURFACE_ABSORB))
return;
State s;
int steps = 0;
while (steps < max_steps) {
steps++;
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, geometry);
if (p.last_hit_triangle == -1)
break;
command = propagate_to_boundary(p, s, rng);
if (command == BREAK)
break;
if (command == CONTINUE)
continue;
if (s.surface_index != -1) {
command = propagate_at_surface(p, s, rng, geometry);
if (command == BREAK)
break;
if (command == CONTINUE)
continue;
}
propagate_at_boundary(p, s, rng);
} // while (steps < max_steps)
rng_states[id] = rng;
positions[photon_id] = p.position;
directions[photon_id] = p.direction;
polarizations[photon_id] = p.polarization;
wavelengths[photon_id] = p.wavelength;
times[photon_id] = p.time;
histories[photon_id] = p.history;
last_hit_triangles[photon_id] = p.last_hit_triangle;
// Not done, put photon in output queue
if ((p.history & (NO_HIT | BULK_ABSORB | SURFACE_DETECT | SURFACE_ABSORB)) == 0) {
int out_idx = atomicAdd(output_queue, 1);
output_queue[out_idx] = photon_id;
}
} // propagate
} // extern "C"
|
