Opticks : GPU Optical Photon Simulation for Particle Physics with NVIDIA OptiX

Opticks Benefits

Outline

• Short video
• Problem and approach to solving
  • Optical Photon Simulation problem
  • Hybrid solution : Geant4 + Opticks
  • NVIDIA OptiX Ray Tracing Engine
• Opticks : optical photon simulation with NVIDIA OptiX
  • geometry workflow : geocache, GPU textures
  • large geometry techniques, geometry modelling
  • Geant4/Opticks event workflow
• Validation, performance comparisons
• Summary

Optical Photon Simulation Problem...

JPMT Before Contact 2

Ray Traced Image Synthesis ≈ Optical Photon Simulation

OptiX Pixel Calculation

http://on-demand.gputechconf.com/siggraph/2013/presentation/SG3106-Building-Ray-Tracing-Applications-OptiX.pdf

Geometry, light sources, optical physics ->

• pixel values at image plane
• photon parameters at detectors (eg PMTs)

Ray tracing has many applications :

• advertising, design, entertainment, games,...
• BUT : most ray tracers just render images

NVIDIA OptiX had foresight to be less specific:

• general geometry intersection API
• OptiX is to ray tracing what OpenGL is to rasterization

rasterization

project 3D primitives onto 2D image plane, combine fragments into pixel values

ray tracing

cast rays thru image pixels into scene, recursively reflect/refract with geometry intersected, combine returns into pixel values

NVIDIA OptiX 1

NVIDIA OptiX 2

https://research.nvidia.com/publication/optix-general-purpose-ray-tracing-engine

Opticks Geometry : Recreate Geant4 "Context" on GPU

Export Tesselated Geometry

• Geant4 -> G4DAE XML https://bitbucket.org/simoncblyth/g4dae

Load G4DAE XML

• label primitives with boundary indices
• find repeated geometry "instances"
• write NumPy serialization .npy geocache

Geocache -> few seconds startup, instead of few minutes

Load/upload geocache

• load geocache from file
• upload OptiX textures to GPU
• upload geometry buffers to GPU

GPU geometry buffers read by:

• OptiX bounding box and intersection programs
• OpenGL GLSL shaders for visualization

Large Geometry Techniques : Instancing Mandatory

Geometry analysed to find instances

JUNO: ~90M --> 0.1M triangles

• 18k 20" PMTs
• 36k 3" PMTs

OpenGL/OptiX instancing

Multiple Renderers

• global, full instances, bbox instances
• rapid switching controls

Geometry Modelling : Tesselated vs Analytic Photomultiplier Tubes

Analytic : more realistic, faster, less memory, much more effort

For Dayabay PMT:

• partition CSG solids into 12 single primitive parts (instead of 2928 triangles)
• splitting at geometrical intersections avoids implementing general CSG boolean handling
• geometry provided to OptiX in form of ray intersection and bounding box code

Aim : analytic description of geometry on critical optical path, remainder tesselated

Hybrid Geant4/Opticks Event Workflow

GPU Resident Photons

Seeded on GPU

associate photons -> gensteps (via seed buffer)

Generated on GPU, using gensteps:

~Stacks collected before photon generation:

• number of photons to generate
• start/end position of step
• other quantities needed for GPU generation

Propagated on GPU

Only photons hitting PMTs copied to CPU

Thrust: high level C++ access to CUDA

• https://developer.nvidia.com/Thrust

Geant4/Detector Simulation

• modified Scintillation, Cerenkov processes
  • collect genstep
  • skip optical photon generation

Opticks (OptiX/Thrust GPU interoperation)

• OptiX : upload gensteps
• Thrust : seeding, distribute genstep indices to photons
• OptiX : launch photon generation and propagation
• Thrust : pullback photons that hit PMTs
• Thrust : index photon step sequences (optional)

Geant4/Detector Simulation

• populate standard hit collections
• subsequent electronics simulation proceeds unaltered

Multi-event handling

• reuse/resize OptiX buffers for each event

Compare Opticks/Geant4 Simulations with Simple Lights/Geometries

1M Photons -> Water Sphere (S-Polarized)

0.5M Photons -> Dayabay PMT

Photon step records

128 bit per step : highly compressed position, time, wavelength, polarization vector, material/history codes

Photon flag sequence

16x 4-bit step flags recorded in uint64 sequence, indexed using Thrust GPU sort (1M indexed ~0.040s)

Sequence index -> interactive OpenGL selection of photons by flag sequence

Opticks/Geant4 Rainbow Step Sequence Comparison

Flags:

• BT/BR: boundary transmit/reflect
• TO/SC/SA: torch/scatter/surface absorb

Statistically consistent photon histories in the two simulations : Multiple orders of rainbow apparent

 64-bit uint  Opticks    Geant4    chi2                                      (tag:5,-5)

        8ccd   819160    819654    0.15  [4 ] TO BT BT SA                    (cross droplet)
         8bd   102087    101615    1.09  [3 ] TO BR SA                       (external reflect)
       8cbcd    61869     61890    0.00  [5 ] TO BT BR BT SA                 (bow 1)
      8cbbcd     9618      9577    0.09  [6 ] TO BT BR BR BT SA              (bow 2)
     8cbbbcd     2604      2687    1.30  [7 ] TO BT BR BR BR BT SA           (bow 3)
    8cbbbbcd     1056      1030    0.32  [8 ] TO BT BR BR BR BR BT SA        (bow 4)
       86ccd     1014      1000    0.10  [5 ] TO BT BT SC SA
   8cbbbbbcd      472       516    1.96  [9 ] TO BT BR BR BR BR BR BT SA     (bow 5)
         86d      498       473    0.64  [3 ] TO SC SA
  bbbbbbbbcd      304       294    0.17  [10] TO BT BR BR BR BR BR BR BR BR  (bow 8+ truncated)
  8cbbbbbbcd      272       247    1.20  [10] TO BT BR BR BR BR BR BR BT SA  (bow 6)
  cbbbbbbbcd      183       161    1.41  [10] TO BT BR BR BR BR BR BR BR BT  (bow 7 truncated)

1M Rainbow S-Polarized, Comparison Opticks/Geant4

Deviation angle(degrees) of 1M parallel monochromatic photons in disc shaped beam incident on water sphere. Numbered bands are visible range expectations of first 11 rainbows. S-Polarized intersection (E field perpendicular to plane of incidence) arranged by directing polarization radially.

PMT Opticks/Geant4 step distribution comparison TO BT [SD]

Good agreement reached, after several fixes: geometry, total internal reflection, group velocity

position(xyz), time(t)

polarization(abc), radius(r)

PMT Opticks/Geant4 step distribution comparison : chi2/ndf

Photon Propagation Times Geant4 cf Opticks

• Opticks > 200X Geant4 with only 384 core mobile GPU[1] (multi-GPU workstation up to 20x more cores)

• Interop uses OpenGL buffers allowing visualization, Compute uses OptiX buffers
• Interop/Compute : perfectly identical results, monitored by digest

[1] MacBook Pro (2013), NVIDIA GeForce GT 750M, 2048 MB, 384 cores

Opticks Event Buffers

Required GPU Buffers

genstep

(n_genstep,6,4) float32: parameters of Cerenkov, Scintillation or Torch genstep. Steps collected from Geant4 are copied to GPU.

seed

(n_photon, 1) uint32, genstep_id for each photon : associates gensteps to photons

photon

(n_photon,4,4) float32: position, time, direction, polarization, wavelength, flags

Debug GPU Buffers

sequence

(n_photon, 1, 2) uint64 : up to 16 flag and material codes packed into 2*64bits

record

(n_photon, 16, 2, 4) int16 : up to 16 propagation points : highly compressed position, time, polarization, wavelength, flags

CPU Buffers

hit

(n_hit, 4,4) float32 : only photons that hit PMTs copied back to CPU

Time/Memory profile of multi-event production mode running

• First event has long prelaunch of 0.37s when GPU buffers are created, subsequent events reuse the buffer, resizing as needed.
• GPU launches every ~0.11s taking ~0.07s (for 0.5M photons) overhead of ~0.04s for downloadHits unusually large from high hit fraction as directly shooting PMT

tpmt-;tpmt--  --compute --multievent 10  --production    # 10x 0.5M photons incident on PMT

2016-11-30 23:00:51.446 INFO  [224583] [TimesTable::dump@103] Opticks::postpropagate filter: NONE

2.656           2.656          2.656      30501.000      30501.000 : OpticksRun::createEvent.BEG_0
0.008           2.664          0.008      30504.000          3.000 : OpSeeder::seedPhotonsFromGenstepsViaOptiX_0
0.367           3.031          0.367      30675.000        171.000 : OPropagator::prelaunch_0
0.070           3.102          0.070      30675.000          0.000 : OPropagator::launch_0
0.031           3.133          0.031      30682.000          7.000 : OEvent::downloadHits_0
0.008           3.141          0.008      30682.000          0.000 : OpticksRun::resetEvent.END_0

0.070           3.211          0.070      30682.000          0.000 : OPropagator::launch_1
0.031           3.242          0.031      30689.000          7.000 : OEvent::downloadHits_1
0.008           3.250          0.008      30689.000          0.000 : OpticksRun::resetEvent.END_1

0.070           3.320          0.070      30689.000          0.000 : OPropagator::launch_2
0.023           3.344          0.023      30696.000          7.000 : OEvent::downloadHits_2
0.008           3.352          0.008      30696.000          0.000 : OpticksRun::saveEvent.END_2
0.008           3.359          0.008      30696.000          0.000 : OpticksRun::resetEvent.END_2

...

Time/Memory profile of multi-event production mode running (2)

• timings of 10 GPU launches, fairly stable
• time until first launch: 3.1s (loading geocache from file, uploading to GPU, creating event buffers)
• memory profiling suggests leaking ~7MB per event, (Needs Investigation)

2016-11-30 23:00:51.447 INFO  [224583] [TimesTable::dump@103] Opticks::postpropagate filter: OPropagator::launch

3.102           3.102          0.070      30675.000          0.000 : OPropagator::launch_0
0.109           3.211          0.070      30682.000          0.000 : OPropagator::launch_1
0.109           3.320          0.070      30689.000          0.000 : OPropagator::launch_2
0.109           3.430          0.070      30696.000          0.000 : OPropagator::launch_3
0.109           3.539          0.070      30703.000          0.000 : OPropagator::launch_4
0.102           3.641          0.070      30710.000          0.000 : OPropagator::launch_5
0.109           3.750          0.070      30716.000          0.000 : OPropagator::launch_6
0.102           3.852          0.070      30723.000          0.000 : OPropagator::launch_7
0.102           3.953          0.062      30730.000          0.000 : OPropagator::launch_8
0.102           4.055          0.062      30737.000          0.000 : OPropagator::launch_9

tconcentric : test geometry configured via boundaries

• Idealized geometry for testing optical physics in situation with no geometry issues.
• Boundaries specify : outer-material/outer-surface/inner-surface/inner-material

tconcentric-testconfig(){
    local test_config=(
        mode=BoxInBox
        analytic=1
        shape=sphere boundary=StainlessSteel///Acrylic         parameters=0,0,0,$(( 5000 + 5 ))
        shape=sphere boundary=Acrylic//RSOilSurface/MineralOil parameters=0,0,0,$(( 5000 - 5 ))
        shape=sphere boundary=MineralOil///Acrylic             parameters=0,0,0,$(( 4000 + 5 ))
        shape=sphere boundary=Acrylic///LiquidScintillator     parameters=0,0,0,$(( 4000 - 5 ))
        shape=sphere boundary=LiquidScintillator///Acrylic     parameters=0,0,0,$(( 3000 + 5 ))
        shape=sphere boundary=Acrylic///GdDopedLS              parameters=0,0,0,$(( 3000 - 5 ))
    )
    echo "$(join _ ${test_config[@]})"
}

https://bitbucket.org/simoncblyth/opticks/src/tip/tests/tconcentric.bash

tconcentric : fixed polarization "Torch" light source

Torch source: 1M 430nm photons (+Y polarized) shot down X axis

• propagation times can be calculated from known material group velocities

tconcentric-torchconfig(){
    local torch_config=(
                 type=point
                 photons=1000000
                 mode=fixpol         polarization=0,1,0
                 source=0,0,0        target=1,0,0
                 material=GdDopedLS  wavelength=430
                 weight=1.0
                 time=0.1
                 zenithazimuth=0,1,0,1
                 radius=0
               )
    echo "$(join _ ${torch_config[@]})"
}

https://bitbucket.org/simoncblyth/opticks/src/tip/tests/tconcentric.bash

tconcentric : spherical GdLS/LS/MineralOil

tconcentric : Opticks/Geant4 chi2 comparison

.      seqhis_ana  1:concentric   -1:concentric     c2
.                       1000000   1000000       373.13/356 =  1.05  (pval:0.256 prob:0.744)
0000           8ccccd    669843    670001             0.02  [6 ] TO BT BT BT BT SA
0001               4d     83950     84149             0.24  [2 ] TO AB
0002          8cccc6d     45490     44770             5.74  [7 ] TO SC BT BT BT BT SA
0003           4ccccd     28955     28718             0.97  [6 ] TO BT BT BT BT AB
0004             4ccd     23187     23170             0.01  [4 ] TO BT BT AB
0005          8cccc5d     20238     20140             0.24  [7 ] TO RE BT BT BT BT SA
0006          8cc6ccd     10214     10357             0.99  [7 ] TO BT BT SC BT BT SA
0007          86ccccd     10176     10318             0.98  [7 ] TO BT BT BT BT SC SA
0008          89ccccd      7540      7710             1.90  [7 ] TO BT BT BT BT DR SA
0009         8cccc55d      5976      5934             0.15  [8 ] TO RE RE BT BT BT BT SA
0010              45d      5779      5766             0.01  [3 ] TO RE AB
0011  8cccccccc9ccccd      5339      5269             0.46  [15] TO BT BT BT BT DR BT BT BT BT BT BT BT BT SA
0012          8cc5ccd      5111      4940             2.91  [7 ] TO BT BT RE BT BT SA
0013              46d      4797      4886             0.82  [3 ] TO SC AB
0014      8cccc9ccccd      4494      4469             0.07  [11] TO BT BT BT BT DR BT BT BT BT SA
0015      8cccccc6ccd      3317      3302             0.03  [11] TO BT BT SC BT BT BT BT BT BT SA
0016         8cccc66d      2670      2675             0.00  [8 ] TO SC SC BT BT BT BT SA
0017          49ccccd      2432      2383             0.50  [7 ] TO BT BT BT BT DR AB
0018          4cccc6d      2043      1991             0.67  [7 ] TO SC BT BT BT BT AB
0019            4cc6d      1755      1826             1.41  [5 ] TO SC BT BT AB

Top 20 chart above, (category 100 down to ~100 photons for propagation of 1M photons)

tconcentric : Opticks/Geant4 history counts chi2/df ~ 1.0

Match achieved between two optical simulations, run from single executable

CPU

Daya Bay detsim optical physics updated to G4 10.2

GPU

above ported to CUDA/OptiX

Match achieved only after many bug fixes, including:

• scattering, was comparing different implementations
• diffuse reflection, ported G4 approach
• reemission "rejoining"
• missing optical surfaces
• truncation recording discrepancy
• group velocity of wrong material after refraction (G4 issue 1275)
• interpolation mismatches

Work guided by the next largest chi2 contributor

https://bugzilla-geant4.kek.jp/show_bug.cgi?id=1275

Group Velocity problems -> Time shifts

Group velocity approximated using "E dn/dE" dispersion term.

• G4 "bin shifting" for dispersion incompatible with GPU texture fixed domains,
• stepLength/groupvel -> propagation time, small groupvel mismatches -> time shifts
• following refraction G4 uses groupvel of wrong material -> time shifts : https://bugzilla-geant4.kek.jp/show_bug.cgi?id=1275

            c                         dn  dE          dn          n1 - n0           delta_n
vg =  --------------- ,    E dn/dE =  --  --     ~  --------  ~  -----------   ~   -----------
        n + E dn/dE                   dE  dlogE      d(logE)     logE1 - logE0      log(E1/E0)

///         0  +  1  +  2  +  3  +  4  +  5        <--- 6 original values
///         |    /     /     /     /      |
///         |   /     /     /     /       |
///         0  1     2     3     4        5        <--- still 6, dispersion term domain

Find photon tracks require UseGivenVelocity=true to get proposed velocity to be used.

// eg from the G4TrackAction
const_cast(track)->UseGivenVelocity(true)

tconcentric : "Unbiased" Scattering : TO SC ..

In [1]: ab.sel = "TO SC .." ; print ab.his
Out[1]:
.        seqhis_ana  1:concentric   -1:concentric      c2
.                         66489     66126        73.48/67 =  1.10  (pval:0.274 prob:0.726)
0000          8cccc6d     45490     44770             5.74  [7 ] TO SC BT BT BT BT SA
0001              46d      4797      4886             0.82  [3 ] TO SC AB
0002         8cccc66d      2670      2675             0.00  [8 ] TO SC SC BT BT BT BT SA
0003          4cccc6d      2043      1991             0.67  [7 ] TO SC BT BT BT BT AB
0004            4cc6d      1755      1826             1.41  [5 ] TO SC BT BT AB
0005         8cccc56d      1158      1215             1.37  [8 ] TO SC RE BT BT BT BT SA
0006         8cc6cc6d       817       739             3.91  [8 ] TO SC BT BT SC BT BT SA
0007         86cccc6d       709       756             1.51  [8 ] TO SC BT BT BT BT SC SA
0008         89cccc6d       538       486             2.64  [8 ] TO SC BT BT BT BT DR SA
0009         8cc5cc6d       381       372             0.11  [8 ] TO SC BT BT RE BT BT SA
0010        8cccc556d       321       369             3.34  [9 ] TO SC RE RE BT BT BT BT SA
0011 8cccccccc9cccc6d       348       367             0.50  [16] TO SC BT BT BT BT DR BT BT BT BT BT BT BT BT SA
0012             456d       292       330             2.32  [4 ] TO SC RE AB
0013     8cccc9cccc6d       316       319             0.01  [12] TO SC BT BT BT BT DR BT BT BT BT SA
0014             466d       277       302             1.08  [4 ] TO SC SC AB
0015     8cccccc6cc6d       255       243             0.29  [12] TO SC BT BT SC BT BT BT BT BT BT SA
0016         49cccc6d       162       158             0.05  [8 ] TO SC BT BT BT BT DR AB
0017        8cccc666d       143       156             0.57  [9 ] TO SC SC SC BT BT BT BT SA
0018     8cccc6cccc6d       121       135             0.77  [12] TO SC BT BT BT BT SC BT BT BT BT SA
0019 8cccccccc6cccc6d       126       108             1.38  [16] TO SC BT BT BT BT SC BT BT BT BT BT BT BT BT SA

tconcentric : "Unbiased" Scattering : TO SC .. polarization

Outgoing polarization

• perpendicular to outgoing direction (photon stays transverse)
• within plane of initial polarization and outgoing direction

tconcentric : Opticks/Geant4 distrib chi2/df ~ 1.0

• Top 100 history categories correspond to ~900 propagation points
• 8 quantities at each point : ~7200 histograms pairs to chi2 compare
• selecting discrepant points : distchi2 > 1.1 (yields 41 out of 900 points)

XYZT:position/time ABCW:polarization/wavelength

OpticksDocs

import numpy as np
a = np.load("photons.npy")

Summary

Opticks enables Geant4 based simulations to benefit from optical photon simulation taking effectively zero time and zero CPU memory, thanks to massive parallelism made accessible by NVIDIA OptiX.

• The more photons the bigger the overall speedup (99% -> 100x)
• Drastic speedup -> better detector understanding -> greater precision
• Large PMT based neutrino experiments, such as JUNO, can benefit the most

https://bitbucket.org/simoncblyth/opticks

YouTube Video

https://www.youtube.com/watch?v=CBpOha4RzIs