Links

Content Skeleton

This Page

Previous topic

Chroma Geometry Source Overview

Next topic

Photon Propagation

Chroma Materials

Material handling in Chroma

(chroma_env)delta:chroma blyth$ grep -l material *.py
detector.py   # passing mention, hands parameter to Geometry subclass of Detector
sim.py        # passing mention, hands detector.detector_material to G4ParallelGenerator
pmt.py        # pmt Solid construction using material arguments to build funcs
geometry.py   # implementation of Solid and Material
(chroma_env)delta:chroma blyth$ find . -name '*.py' -exec grep -l material {} \;
./demo/optics.py
./demo/pmt.py
./detector.py
./generator/g4gen.py    # using Material to create G4Material instances, for photon generation within the material
./generator/photon.py
./geometry.py
./gpu/geometry.py
./pmt.py
./sim.py

chroma/geometry.py

Solid

  • constant or iterator inner/outer material parameters

Material

  • constant or over standard wavelengths properties
015 # all material/surface properties are interpolated at these
016 # wavelengths when they are sent to the gpu
017 standard_wavelengths = np.arange(60, 810, 20).astype(np.float32)
...
200 class Material(object):
201     """Material optical properties."""
202     def __init__(self, name='none'):
203         self.name = name
204
205         self.refractive_index = None
206         self.absorption_length = None
207         self.scattering_length = None
208         self.set('reemission_prob', 0)
209         self.set('reemission_cdf', 0)
210         self.density = 0.0 # g/cm^3
211         self.composition = {} # by mass
212
213     def set(self, name, value, wavelengths=standard_wavelengths):
214         if np.iterable(value):
215             if len(value) != len(wavelengths):
216                 raise ValueError('shape mismatch')
217         else:
218             value = np.tile(value, len(wavelengths))
219
220         self.__dict__[name] = np.array(zip(wavelengths, value), dtype=np.float32)
221
222 # Empty material
223 vacuum = Material('vacuum')
224 vacuum.set('refractive_index', 1.0)
225 vacuum.set('absorption_length', 1e6)
226 vacuum.set('scattering_length', 1e6)

Surface

Wavelength dependent surface properties handled similarly:

229 class Surface(object):
230     """Surface optical properties."""
231     def __init__(self, name='none', model=0):
232         self.name = name
233         self.model = model
234
235         self.set('detect', 0)
236         self.set('absorb', 0)
237         self.set('reemit', 0)
238         self.set('reflect_diffuse', 0)
239         self.set('reflect_specular', 0)
240         self.set('eta', 0)
241         self.set('k', 0)
242         self.set('reemission_cdf', 0)
243
244         self.thickness = 0.0
245         self.transmissive = 0

chroma/gpu/geometry.py

Material properties passed over to GPU side

  1. refractive_index
  2. absorption_length
  3. scattering_length
  4. reemission_prob
  5. reemission_cdf

Looks like (from chroma/sim.py) have just one GPUGeometry.

13 class GPUGeometry(object):
14     def __init__(self, geometry, wavelengths=None, print_usage=False, min_free_gpu_mem=300e6):
15         if wavelengths is None:
16             wavelengths = standard_wavelengths
17
18         try:
19             wavelength_step = np.unique(np.diff(wavelengths)).item()
20         except ValueError:
21             raise ValueError('wavelengths must be equally spaced apart.')
22
23         geometry_source = get_cu_source('geometry_types.h')
24         material_struct_size = characterize.sizeof('Material', geometry_source)
25         surface_struct_size = characterize.sizeof('Surface', geometry_source)
26         geometry_struct_size = characterize.sizeof('Geometry', geometry_source)
27
28         self.material_data = []
29         self.material_ptrs = []
30
31         def interp_material_property(wavelengths, property):
32             # note that it is essential that the material properties be
33             # interpolated linearly. this fact is used in the propagation
34             # code to guarantee that probabilities still sum to one.
35             return np.interp(wavelengths, property[:,0], property[:,1]).astype(np.float32)
36
37         for i in range(len(geometry.unique_materials)):
38             material = geometry.unique_materials[i]
39
40             if material is None:
41                 raise Exception('one or more triangles is missing a material.')
42
43             refractive_index = interp_material_property(wavelengths, material.refractive_index)
44             refractive_index_gpu = ga.to_gpu(refractive_index)
45             absorption_length = interp_material_property(wavelengths, material.absorption_length)
46             absorption_length_gpu = ga.to_gpu(absorption_length)
47             scattering_length = interp_material_property(wavelengths, material.scattering_length)
48             scattering_length_gpu = ga.to_gpu(scattering_length)
49             reemission_prob = interp_material_property(wavelengths, material.reemission_prob)
50             reemission_prob_gpu = ga.to_gpu(reemission_prob)
51             reemission_cdf = interp_material_property(wavelengths, material.reemission_cdf)
52             reemission_cdf_gpu = ga.to_gpu(reemission_cdf)
..

NB all five material properties are linearly interpolated onto the same standard (or provided) wavelengths.

Same class also passes surface properties

  1. detect
  2. absorb
  3. reemit
  4. reflect_diffuse
  5. reflect_specular
  6. eta
  7. k
  8. reemission_cdf
78         for i in range(len(geometry.unique_surfaces)):
79             surface = geometry.unique_surfaces[i]
80
81             if surface is None:
82                 # need something to copy to the surface array struct
83                 # that is the same size as a 64-bit pointer.
84                 # this pointer will never be used by the simulation.
85                 self.surface_ptrs.append(np.uint64(0))
86                 continue
87
88             detect = interp_material_property(wavelengths, surface.detect)
89             detect_gpu = ga.to_gpu(detect)
90             absorb = interp_material_property(wavelengths, surface.absorb)
91             absorb_gpu = ga.to_gpu(absorb)
92             reemit = interp_material_property(wavelengths, surface.reemit)
93             reemit_gpu = ga.to_gpu(reemit)
94             reflect_diffuse = interp_material_property(wavelengths, surface.reflect_diffuse)
95             reflect_diffuse_gpu = ga.to_gpu(reflect_diffuse)
96             reflect_specular = interp_material_property(wavelengths, surface.reflect_specular)
97             reflect_specular_gpu = ga.to_gpu(reflect_specular)
98             eta = interp_material_property(wavelengths, surface.eta)
99             eta_gpu = ga.to_gpu(eta)
100             k = interp_material_property(wavelengths, surface.k)
101             k_gpu = ga.to_gpu(k)
102             reemission_cdf = interp_material_property(wavelengths, surface.reemission_cdf)
103             reemission_cdf_gpu = ga.to_gpu(reemission_cdf)
104
105             self.surface_data.append(detect_gpu)
106             self.surface_data.append(absorb_gpu)
107             self.surface_data.append(reemit_gpu)
108             self.surface_data.append(reflect_diffuse_gpu)
109             self.surface_data.append(reflect_specular_gpu)
110             self.surface_data.append(eta_gpu)
111             self.surface_data.append(k_gpu)
112             self.surface_data.append(reemission_cdf_gpu)
113
114             surface_gpu = \
115                 make_gpu_struct(surface_struct_size,
116                                 [detect_gpu, absorb_gpu, reemit_gpu,
117                                  reflect_diffuse_gpu,reflect_specular_gpu,
118                                  eta_gpu, k_gpu, reemission_cdf_gpu,
119                                  np.uint32(surface.model),
120                                  np.uint32(len(wavelengths)),
121                                  np.uint32(surface.transmissive),
122                                  np.float32(wavelength_step),
123                                  np.float32(wavelengths[0]),
124                                  np.float32(surface.thickness)])
125
126             self.surface_ptrs.append(surface_gpu)

chroma/cuda/propagate.cu

propagate

Within the propagate stepping the fill_state(s, p, g) state, photon, geometry sets material props with the state.

152     if (p.history & (NO_HIT | BULK_ABSORB | SURFACE_DETECT | SURFACE_ABSORB | NAN_ABORT))
153     return;
///
///     FLAGGED AS A DEAD PHOTON ALREADY, NOTHING TO DO
///
154
155     State s;
156
157     int steps = 0;
158     while (steps < max_steps) {
159     steps++;
160
161     int command;
162
163     // check for NaN and fail
164     if (isnan(p.direction.x*p.direction.y*p.direction.z*p.position.x*p.position.y*p.position.z)) {
165         p.history |= NO_HIT | NAN_ABORT;
166         break;
167     }
168
169     fill_state(s, p, g);
170
171     if (p.last_hit_triangle == -1)
172         break;
173
174     command = propagate_to_boundary(p, s, rng, use_weights, scatter_first);
175     scatter_first = 0; // Only use the scatter_first value once
176
177     if (command == BREAK)
178         break;
179
180     if (command == CONTINUE)
181         continue;
182
183     if (s.surface_index != -1) {
184       command = propagate_at_surface(p, s, rng, g, use_weights);
185
186         if (command == BREAK)
187         break;
188
189         if (command == CONTINUE)
190         continue;
191     }
192
193     propagate_at_boundary(p, s, rng);
194
195     } // while (steps < max_steps)

chroma/cuda/photon.h

State

30 struct State
31 {
32     bool inside_to_outside;
33
34     float3 surface_normal;
35
36     float refractive_index1, refractive_index2;
37     float absorption_length;
38     float scattering_length;
39     float reemission_prob;
40     Material *material1;
41
42     int surface_index;
43
44     float distance_to_boundary;
45 };

fill_state

79 __device__ void
80 fill_state(State &s, Photon &p, Geometry *g)
81 {
82     p.last_hit_triangle = intersect_mesh(p.position, p.direction, g,
83                                          s.distance_to_boundary,
84                                          p.last_hit_triangle);
85
86     if (p.last_hit_triangle == -1) {
87         p.history |= NO_HIT;
88         return;
89     }
90
91     Triangle t = get_triangle(g, p.last_hit_triangle);
92
93     unsigned int material_code = g->material_codes[p.last_hit_triangle];
94
95     int inner_material_index = convert(0xFF & (material_code >> 24));
96     int outer_material_index = convert(0xFF & (material_code >> 16));
97     s.surface_index = convert(0xFF & (material_code >> 8));
98
99     float3 v01 = t.v1 - t.v0;
100     float3 v12 = t.v2 - t.v1;
101
102     s.surface_normal = normalize(cross(v01, v12));
103
104     Material *material1, *material2;
105     if (dot(s.surface_normal,-p.direction) > 0.0f) {
106         // outside to inside
107         material1 = g->materials[outer_material_index];
108         material2 = g->materials[inner_material_index];
109
110         s.inside_to_outside = false;
111     }
112     else {
113         // inside to outside
114         material1 = g->materials[inner_material_index];
115         material2 = g->materials[outer_material_index];
116         s.surface_normal = -s.surface_normal;
117
118         s.inside_to_outside = true;
119     }
120
121     s.refractive_index1 = interp_property(material1, p.wavelength,
122                                           material1->refractive_index);
123     s.refractive_index2 = interp_property(material2, p.wavelength,
124                                           material2->refractive_index);
125     s.absorption_length = interp_property(material1, p.wavelength,
126                                           material1->absorption_length);
127     s.scattering_length = interp_property(material1, p.wavelength,
128                                           material1->scattering_length);
129     s.reemission_prob = interp_property(material1, p.wavelength,
130                                         material1->reemission_prob);
131
132     s.material1 = material1;
133 } // fill_state
  1. for COLLADA integration need to implement the GDML G4 material (and surface) wavelength array properties into COLLADA extra tags

material_codes

packed in chroma/geometry.py 4 bytes, high to low: material1/material2/surface/blank

material_codes = (((geometry.material1_index & 0xff) << 24) |
                  ((geometry.material2_index & 0xff) << 16) |
                  ((geometry.surface_index & 0xff) << 8)).astype(np.uint32)
self.material_codes = ga.to_gpu(material_codes)

unpacked:

92
93     unsigned int material_code = g->material_codes[p.last_hit_triangle];
94
95     int inner_material_index = convert(0xFF & (material_code >> 24));
96     int outer_material_index = convert(0xFF & (material_code >> 16));
97     s.surface_index = convert(0xFF & (material_code >> 8));

suggests correspondence:

  • material1 is inner
  • material2 is outer