Why not Chroma ?
Chroma Features
- Python/PyCUDA/NumPy based infrastructure for geometry/photon loading, kernel launch
- accelerated geometry intersection using BVH structure
- optical photon simulation CUDA kernels
My additions to Chroma
- G4DAE Geometry import
- G4Step transport and Cerenkov/Scintillation photon generation on GPU
- OpenGL/CUDA interop visualisations
Missing Features
- GPU workload scheduling
- Multi-GPU support
https://bitbucket.org/chroma/chroma
https://bitbucket.org/simoncblyth/chroma (my fork)
Parallels between Realistic Image Synthesis and Optical Simulation
Realistic image creation uses physically
based techniques and material definitions. Obvious parallels:
- ray traced rendering : image pixel calculation
- optical photon (OP) simulation : PMT hit calculation
Same rate determining step: geometry intersection
Applying techniques/hardware developed for fast ray tracing
can be hugely beneficial to optical photon simulation.
- expect OP simulation performance >100x Geant4
- OP processing time becomes effectively zero
Chroma Raycast with entire geometry in view
Render Split into 3x3 CUDA kernel launches, 1 thread per pixel, ~1.8s for 1.23M pixels, 2.4M tris (with )
OptiX Programming Model
OptiX provides a ray tracing pipeline analogous to OpenGL
rasterization pipeline.
- Blue: OptiX internals
- Yellow: User supplied CUDA Programs
Higher level API than pure CUDA, eg:
- shared host/device context system
Optical Photon Simulation port currently using:
- Ray Generation
- entry/exit, Cerenkov/Scintillation generation
- Intersection
- Triangle mesh intersection, boundary index lookup
- Closest Hit
- determine ray to boundary orientation
https://research.nvidia.com/sites/default/files/publications/Parker10Optix_1.pdf
OptiX Adoption Costs
Adoption of OptiX is compelling
- extremely fast intersection performance
- scales with CUDA cores across multiple GPUs
- improves with each release
- releases tune for new GPU architectures
Costs of adoption
- learn new tools: OptiX, Thrust
- develop C++ replacement for Chroma/G4DAE python
- port Chroma/G4 optical physics into new framework
- validate against Geant4
- develop memory efficient geometry representation
New C++ Packages Replacing Chroma
Basis packages
- NPY : host array handling, persistency
- NumpyServer : network IO of NPY arrays
- BCfg BRegex : configuration
Geometry packages
- GGeo : preparing geometry for GPU
- AssimpWrap : G4DAE geometry loading using forked Assimp
GPU library interface packages
- CUDAWrap : pseudo random numbers with cuRAND
- OptiXRap : geometry loading, OptiX launch control
- OGLWrap : OpenGL visualization
- ThrustRap : photon indexing, using CUDA Thrust
Main package
- GGeoView : optical physics OptiX programs, OpenGL shaders
https://bitbucket.org/simoncblyth/env/src/tip/graphics/ggeoview/
Selection of GPU development details
Some details of GPU developments described over the next pages
Optical Physics
- Porting Optical Physics from Geant4/Chroma into OptiX
- Optical Physics Implementation
Supplying the OptiX Programs
- Random Number Generation in OptiX programs
- Fast material/surface property lookup from boundary texture
- Reemission wavelength lookup from Inverted CDF texture
Handling Outputs
- Recording the steps of millions photons
- Indexing photon flag/material sequences
- Introducing CUDA Thrust
Porting Optical Physics from Geant4/Chroma into OptiX
Photon p ; State s ; PerRayData prd ;
while(bounce < bounce_max) // PSEUDO-CODE
{
bounce++
ray = optix::make_Ray(p.pos, p.dir,...)
rtTrace(geom, ray, prd)
if(!prd.boundary) break // MISS
cmd = propagate_to_boundary(p, s)
if(cmd == BREAK) break // ABSORB
if(cmd == CONTINUE) continue // REEMIT, SCATTER
// survivors pass to boundary
if(s.surface_index)
{
cmd = propagate_at_surface(p, s, g)
if(cmd == BREAK) break // SURFACE_ABSORB, SURFACE_DETECT
if(cmd == CONTINUE) continue // REFLECT_DIFFUSE, REFLECT_SPECULAR
}
else
{
propagate_at_boundary(p, s) // BOUNDARY_REFLECT BOUNDARY_TRANSMIT
}
}
https://bitbucket.org/simoncblyth/env/src/tip/graphics/ggeoview/cu/generate.cu
Random Number Generation in OptiX programs
cuRAND library from CUDA toolkit features:
- concurrent generation of reproducible pseudorandom number sequences
- sub-sequences are assigned to each CUDA thread, which maintains position in sub-sequence
- per-thread state is initialized within CUDA kernel
cuRAND Initialization demands large stack size
Stack sizes 10x typical for OptiX programs were needed,
resulting in slow OptiX running.
Workaround:
- use separate pure CUDA launches to initialize cuRAND
- copy curandState back to host and persist to file
- prior to OptiX launch, copy persisted curandState to GPU
- OptiX can then use cuRAND without having to initialize it
Packaged solution into CUDAWrap
Fast material/surface property lookup from boundary texture
AssimpWrap creates GGeo boundary instances and labels
triangles with boundary indices, boundaries contain:
- inner material : self
- outer material : parent
- inner surface : outwards going photons (self to parent)
- outer surface : inwards going photons (parent to self)
Properties are interpolated onto a common wavelength domain
- material : refractive_index, absorption_length, scattering_length, reemission_prob (float4)
- surface : detect, absorb, reflect_specular, reflect_diffuse (float4)
Interleaved properties used to create single boundary texture 2d (wavelength, qty line)
containing ~50 boundaries, 4 float4 each. CUDA tex2d property lookup:
float nmi = (nm - wavelength_domain.x)/wavelength_domain.z + 0.5f ;
float4 material1 = tex2D(wavelength_texture, nmi, line + 0.5f );
float refractive_index = material1.x ;
float absorption_length = material1.y ;
float scattering_length = material1.z ;
https://github.com/simoncblyth/assimp (my fork of Assimp)
Reemission wavelength lookup from Inverted CDF texture
Inverting Reemission CDF allows using texture lookup to obtain
reemission wavelength from uniform random throws.
Using 4096 probability bins.
Recording the steps of ~3 million photons
Up to 10 steps of the photon propagation are recorded.
Photon buffer : 4 * float4 = 512 bits/photon
- float4: position, time [32 * 4 = 128 bits]
- float4: direction, weight
- float4: polarization, wavelength
- float4: flags: material, boundary, history
Record buffer : 2 * short4 = 2*16*4 = 128 bits/record
- short4: position, time (snorm compressed) [4*16 = 64 bits]
- uchar4: polarization, wavelength (uchar compressed) [4*8 = 32 bits]
- uchar4: material, history flags [4*8 = 32 bits]
Compression uses known domains of position (geometry center, extent),
time (0:200ns), wavelength, polarization.
Union trickery allows recording ints into floats
GGeoView M1 Points
Indexing photon flag/material sequences
Selecting photons by flag/material sequences, requires indexing integer sequences.
- CK : Cerenkov
- BT : Boundary Transmit
- BR : Boundary Reflect
- RE : Reemission
- AB : Absorb
- BS : Bulk Scatter
Indexing history/material sequences for 3M photons:
- CPU STL map, sstream ~40s
- CUDA Thrust sorted sparse histogram <0.4s
Packaged indexing into ThrustRap ThrustIdx
GGeoView Flag Selection
Introducing CUDA Thrust
Distributed with CUDA
GPU performance without developing CUDA kernels
Mobile GPU Timings for Cerenkov and Scintillation photons
max_record:10
max_bounce:9
Cerenkov Ck*4.59 Scintillation
photons 0.61M 2.8M 2.8M
--(bytes)--------------------------------------
genstep size 736K 1.3M
photons size 37M 172M
records size 97M 430M
--(seconds)------------------------------------
createOpenGLCtx 0.692 - 0.599
loadGeometry 1.570 - 1.302
interpGeometry 0.211 - 0.190
initOptiX 4.216 - 6.521
loadGenstep 0.011 - 0.014
hostEvtAllocation ** 3.540 16.275 16.179
uploadEvt 0.232 1.066 0.552
generatePropagate ++ 1.404 6.453 7.907
evtDownload ** 0.348 1.602 1.780
evtSave ** 0.437 2.008 2.006
sequenceIndex 0.134 0.614 0.359
-----------------------------------------------
** scales by photon count
Operation with JUNO Geometry ?
The large number of PMTs may require a more memory
efficient geometry representation using OptiX features:
- parameterized geometry avoids tesselation, like Geant4
- geometry instancing avoids duplication
Memory access (not calculation) typically limits GPU performance,
improving memory efficiency expected to improve performance.
Next Steps
Test New Framework with IHEP 4-GPU workstation (together with Tao Lin)
- check performance scaling across 4-GPU cores
- attempt loading JUNO geometry
- investigate more memory efficient geometry techniques
Optical Photon Simulation
- Complete porting Optical Physics
- Instrument Geant4 optical photon propagation,
by recording photon steps into NPY array to enable step-by-step comparison
- Debugging to achieve match between Geant4 and GPU optical photon simulation,
G4DAE Geometry Exporter
- investigate issue inherited from GDML of a skipped
edge case (when a volume is shared between multiple volume pairs)
resulting in missing G4LogicalBorderSurface
- incorporate into Geant4 codebase
"Backup" : Details for Reference
On the following pages:
- GGeo/OptiX Generated Scintillation Photons cf Geant4
- GGeo/OptiX Generated Cerenkov Photons cf Geant4
- Cerenkov Photon Steps
- C++ Infrastructure : foundation packages
- C++ Infrastructure : domain packages
- propagate_to_boundary : ABSORB(REEMIT) / SCATTER / survive
- Comparison of GGeo/OptiX Generated Scintillation Photon Distributions
- Comparison of GGeo/OptiX Generated Cerenkov Photon Distributions
GGeo/OptiX Generated Scintillation Photons cf Geant4
GGeo/OptiX using inverted CDF reemission wavelength lookups (4096 bins)
GGeo/OptiX Generated Cerenkov Photons cf Geant4
Geant4/DetSim wavelength distribution has a blip at 200nm, corresponding to edge of water
refractive index properties.
GGeoView Cerenkov Geom M1
C++ Infrastructure : foundation packages
- NPY
Array persistency/manipulations inspired by NumPy,
using NPY serialization format
- 11 classes: G4StepNPY, PhotonsNPY, NPY, ...
- NumpyServer
Asynchronous IO of Geant4 Steps, Photons, Hits.
Communicates with remote G4DAEOpticks process, receiving
steps and replying with hits.
- 7 classes : numpydelegate, udp_server, ...
- CUDAWrap
cuRAND init and persist curandState (pure CUDA)
- avoids large stack size requirement of cuRAND init within OptiX
- 5 classes : cuRANDWrapper, LaunchSequence, LaunchCommon, ..
https://bitbucket.org/simoncblyth/env/src/tip/numerics/npy/
https://bitbucket.org/simoncblyth/env/src/tip/boost/basio/numpyserver/
https://bitbucket.org/simoncblyth/env/src/tip/cuda/cudawrap/
C++ Infrastructure : domain packages
- GGeo
GPU Geometry representation, NPY persistency
- 22 classes: GNode, GMaterial, GProperty, ...
- AssimpWrap
G4DAE -> GGeo geometry
- 7 classes : AssimpGGeo, AssimpTree, ...
- OptiXRap
GGeo -> OptiX geometry, OptiX launch control
- 7 classes : OptiXEngine, OptixGeometry, ...
- OGLRap
OpenGL shader based 3D visualization
- 29 classes : Scene, View, Camera, Rdr, Shdr, ...
https://bitbucket.org/simoncblyth/env/src/tip/optix/ggeo/
https://bitbucket.org/simoncblyth/env/src/tip/graphics/assimpwrap/
https://bitbucket.org/simoncblyth/env/src/tip/graphics/oglrap/
https://bitbucket.org/simoncblyth/env/src/tip/graphics/optixrap/
propagate_to_boundary : ABSORB(REEMIT) / SCATTER / survive
__device__ int propagate_to_boundary( Photon& p, State& s, curandState &rng)
{
float absorption_distance = -s.material1.y*logf(curand_uniform(&rng)); // .y:absorption_length
float scattering_distance = -s.material1.z*logf(curand_uniform(&rng)); // .z:scattering_length
if (absorption_distance <= scattering_distance)
{
if (absorption_distance <= s.distance_to_boundary)
{
p.time += absorption_distance/(SPEED_OF_LIGHT/s.material1.x); // .x:refractive_index
p.position += absorption_distance*p.direction;
if (curand_uniform(&rng) < s.material1.w) // .w:reemission_prob
{
// non-scintillators have zero reemission_prob
p.wavelength = reemission_lookup(curand_uniform(&rng));
p.direction = uniform_sphere(&rng);
p.polarization = normalize(cross(uniform_sphere(&rng), p.direction));
s.flag = BULK_REEMIT ;
return CONTINUE;
}
else
{
s.flag = BULK_ABSORB ;
return BREAK;
}
}
// otherwise sail to boundary
}
else
// scattering ..
Comparison of GGeo/OptiX Generated Scintillation Photon Distributions
Position, direction, polarization XYZ + time, wavelength, weight
Comparison of GGeo/OptiX Generated Cerenkov Photon Distributions
Position, direction, polarization XYZ + time, wavelength, weight
GGeoView
