Field Propagation using Runge-Kutta integration and example15 as demo (

#166)

* First version of classes for Runge-Kutta integration of track in field

- MagneticFieldEquation
- DormandPrince45 - Stepper
- RkIntegrationDriver
- fieldConstants.h with kB2C, delta_intersection parameters
- fieldPropagatorRungeKutta : new class for use by electrons.cu, used in new Example14

More:

Example15: field propagation using Runge-Kutta.  Based on Example13 otherwise.

Tests - unit test test_magfieldRK.cpp
  Simple checks for equation, stepper and driver classes
  Driver: check version of Advance and V1 (old)
  Checks driver vs helix results (on cpu).

Details:
--------
RkIntegrationDriver
   Introduced new Advance method.  Using AdvanceV1 as backward compatible.
   Enabled used of RK classes with double integrands (field remains float.)

 PrintFieldVectors: auxiliary methods, initially host-only
 Changed VECCORE_ATT_HOST_DEVICE to __host__ __device__

example15: Added ability to print Track info
           check result of RK integration in electrons.cu
           report differences
           Optional argument for Bz field value.
           Use TrackML as default geometry.

* Electrons.cu + RK: return iterations; print good & bad steps.

* Runge-Kutta integration: changes to identify cause of errors. (For debugging)

* Fix for compilation - due to change to use fieldConstants.h

* electrons.cu: Adapted to changes in G4HepEM (using Example13)

Took revisions from Example13/electrons.cu

* Example15: Fixed multiple parts using changes in Example13

  - electrons.cu
  - example15.cpp
  - example15.cu
  - example15.cuh
  - example15.h

Changed fieldPropagatorRungeKutta.h to suppress 'id' in argument

* Fix to pass world volume to G4TransporationManager

* 15/electrons: Improved printing of differences RK vs helix

* fieldPropagatorRungeKutta: added checks vs Helix

* Example15 / magfield : revised prints, new header file.

* fieldPropagatorRungeKutta: fix update of momentum in substep of ComputeStepAndNextVolume

* fieldPropagatorRungeKutta: Refine intersection for comparisons

* Example15 / magneticfield: checks first refined, then disabled.

* fieldPropagatorRungeKutta: IntegrateToEnd reports if loopCt > 1

RkIntegrationDriver: small cleanup

* fieldPropagatorRungeKutta: fix condition for ending - adding condition (curvedStep > 0)

* Example15: improvements to printing changes (investigating).

* fieldPropagatorRungeKutta: added step reduction for zero steps (copying from fpConstBz)

* example15.cu: fix call to ensure that ClearQueue occurs (not just when printing.)

* example15/electrons.cu: first cleanup, removed max_step=0.25mm, always RK.

* fieldPropagatorRungeKutta.h: sharper step reduction for stuck tracks, flip side if failing

After a maximum number of failed (zero) steps, reducing the step size each time,
accept that the track will not cross the boundary, and flip it to the other volume,
which it is apparently entering (or never left.)

* example15.cu: Added optional printing and info printing ;  Revised magfield/inc/CompareResponses

* example15: reset cms2018.gdml as default geometryf

* example15.cu: removed most (optional) verbosity

* fieldPropagatorConstBz::ComputeNextStepAndVolume: moved slot after safety in args

* Example13/electrons.cu : added arguments to ComputeStepAndNextVolume

* fieldPropagatorConstBz: After multiple zero steps move to other side of
boundary;

fieldPropagatorConstBz::ComputeStepAndNextVolume
- Optional defaul values in interface (at compile time)
- Protected some verbosity using if(verbose)
- Added some optional verbosity

Both of these are to be trimmed / refined.

* fieldPropagatorConstBz: cleaned up several printouts; defaults in ComputeNextStepAndVolume

Enable defaults for the last 3 arguments of ComputeNextStepAndVolume

* fieldPropagatorRungeKutta: removed optional prints & cleanup

* fieldPropagatorConstBz::ComputeStepAndNextVolume enforce that safety arg is const

* fieldPropagatorConstBz.h : clean-up; added const in arg

Safety was const in ComputeStepAndNextVolume implementation only - added it to declaration.

* fieldPropagatorConstBz: fix vecgeom::Precision; cosmetics

* Example13: set B=3.8T ; electrons.cu: took out slot in args of fieldProp...

* CompareResponces.h: Added Copyright, reuse

* Example15/electrons.cu : cleanup. Fix to remove slot.

* example15.cu: cleanup

* example15.cpp: Set stack limit to 8K

examples/CMakeLists.txt

@@ -17,6 +17,7 @@ add_subdirectory(Example11)
 add_subdirectory(Example12)
 add_subdirectory(Example13)
 add_subdirectory(Example14)
+add_subdirectory(Example15)
 add_subdirectory(Example16)
 add_subdirectory(Example17)
 add_subdirectory(Example18)

examples/Example11/README.md

@@ -19,7 +19,7 @@ Results are reproducible using one RANLUX++ state per track.
 ### Kernels
 
 This example uses one stream per particle type to launch kernels asynchronously.
-They are synchronized via a forth stream using CUDA events.
+They are synchronized via a fourth stream using CUDA events.
 
 #### `TransportElectrons<bool IsElectron>`
 

examples/Example12/README.md

@@ -22,7 +22,7 @@ Additionally, the kernels score energy deposit and the charged track length per
 ### Kernels
 
 This example uses one stream per particle type to launch kernels asynchronously.
-They are synchronized via a forth stream using CUDA events.
+They are synchronized via a fourth stream using CUDA events.
 
 #### `TransportElectrons<bool IsElectron>`
 

examples/Example13/README.md

@@ -24,7 +24,7 @@ Additionally, the kernels score energy deposit and the charged track length per
 ### Kernels
 
 This example uses one stream per particle type to launch kernels asynchronously.
-They are synchronized via a forth stream using CUDA events.
+They are synchronized via a fourth stream using CUDA events.
 
 #### `TransportElectrons<bool IsElectron>`
 

examples/Example13/electrons.cu

@@ -123,9 +123,12 @@ static __device__ __forceinline__ void TransportElectrons(Track *electrons, cons
     bool propagated = true;
     double geometryStepLength;
     vecgeom::NavStateIndex nextState;
+    constexpr int max_iters = 10;
+    // printf("max_iters= %3d\n", max_iters);
     if (BzFieldValue != 0) {
       geometryStepLength = fieldPropagatorBz.ComputeStepAndNextVolume<BVHNavigator>(
-          energy, Mass, Charge, geometricalStepLengthFromPhysics, pos, dir, navState, nextState, propagated, safety);
+          energy, Mass, Charge, geometricalStepLengthFromPhysics, pos, dir, navState, nextState,
+          propagated, safety, max_iters );
     } else {
       geometryStepLength = BVHNavigator::ComputeStepAndNextVolume(pos, dir, geometricalStepLengthFromPhysics, navState,
                                                                   nextState, kPush);

examples/Example13/example13.cuh

@@ -135,7 +135,7 @@ extern __constant__ __device__ struct G4HepEmData g4HepEmData;
 
 extern __constant__ __device__ int *MCIndex;
 
-// constexpr vecgeom::Precision BzFieldValue = 0.1 * copcore::units::tesla;
-constexpr vecgeom::Precision BzFieldValue = 0;
+constexpr vecgeom::Precision BzFieldValue = 3.8 * copcore::units::tesla;
+// constexpr vecgeom::Precision BzFieldValue = 0;
 
 #endif

examples/Example15/CMakeLists.txt

@@ -0,0 +1,42 @@
+# SPDX-FileCopyrightText: 2021 CERN
+# SPDX-License-Identifier: Apache-2.0
+
+if(NOT TARGET G4HepEm::g4HepEm)
+  message(STATUS "Disabling example15 (needs G4HepEm)")
+  return()
+endif()
+
+if(Geant4_FOUND)
+  if(NOT Geant4_gdml_FOUND)
+    message(STATUS "Disabling example15 (needs Geant4 with GDML support)")
+    return()
+  endif()
+else()
+  message(STATUS "Disabling example15 (needs Geant4)")
+  return()
+endif()
+
+# example15 is the AdePT demo example using material cuts as defined in the input gdml file
+add_executable(example15
+  example15.cpp
+  example15.cu
+  electrons.cu
+  gammas.cu)
+target_link_libraries(example15
+  PRIVATE
+    AdePT
+    CopCore::CopCore
+    VecGeom::vecgeom
+    VecGeom::vecgeomcuda_static
+    VecGeom::vgdml
+    ${Geant4_LIBRARIES}
+    G4HepEm::g4HepEmData
+    G4HepEm::g4HepEmInit
+    G4HepEm::g4HepEmRun
+    CUDA::cudart)
+
+set_target_properties(example15 PROPERTIES CUDA_SEPARABLE_COMPILATION ON CUDA_RESOLVE_DEVICE_SYMBOLS ON)
+
+# Tests
+add_test(NAME example15
+  COMMAND $<TARGET_FILE:example15> -gdml_file ${CMAKE_BINARY_DIR}/cms2018.gdml)

examples/Example15/README.md

@@ -0,0 +1,64 @@
+<!--
+SPDX-FileCopyrightText: 2021 CERN
+SPDX-License-Identifier: CC-BY-4.0
+-->
+
+## Example 14
+
+Example demonstrating particle transportation on GPUs with a non-uniform magnetic field - in arbitrary geometry read from a GDML file.
+
+New feature is
+ - definition of a non-uniform magnetic field
+ - integration of the equation of motion of charged particles using embedded Runge-Kutta method, Dormand Prince 4(5).
+ Caveats (at present) : 
+   - fixed accuracy of integration (constant in source)
+
+The features carried over from Example13 are:
+ * arbitrary geometry via gdml file (tested with cms2018.gdml from VecGeom persistency/gdml/gdmls folder) and optionally a magnetic field with constant Bz,
+ * geometry read as Geant4 geometry, reading in regions and cuts, to initialize G4HepEm data
+ * geometry read then into VecGeom, and synchronized to GPU
+ * G4HepEm material-cuts couple indices mapped to VecGeom logical volume id's
+ * physics processes for e-/e+ (including MSC) and gammas using G4HepEm.
+ * scoring per placed volume, no sensitive detector feature
+
+Electrons, positrons, and gammas are stored in separate containers in device memory.
+Free positions in the storage are handed out with monotonic slot numbers, slots are not reused.
+Active tracks are passed via three queues per particle type (see `struct ParticleQueues`).
+Results are reproducible using one RANLUX++ state per track.
+
+Additionally, the kernels score energy deposit and the charged track length per volume.
+
+### Kernels
+
+This example uses one stream per particle type to launch kernels asynchronously.
+They are synchronized via a fourth stream using CUDA events.
+
+#### `TransportElectrons<bool IsElectron>`
+
+1. Obtain safety unless the track is currently on a boundary.
+2. Determine physics step limit, including conversion to geometric step length according to MSC.
+3. Query geometry (or optionally magnetic field) to get geometry step length.
+4. Convert geometry to true step length according to MSC, apply net direction change and discplacement.
+5. Apply continuous effects; kill track if stopped.
+6. If the particle reaches a boundary, perform relocation.
+7. If not, and if there is a discrete process:
+ 1. Sample the final state.
+ 2. Update the primary and produce secondaries.
+
+#### `TransportGammas`
+
+1. Determine the physics step limit.
+2. Query VecGeom to get geometry step length (no magnetic field for neutral particles!).
+3. If the particle reaches a boundary, perform relocation.
+4. If not, and if there is a discrete process:
+ 1. Sample the final state.
+ 2. Update the primary and produce secondaries.
+
+#### `FinishIteration`
+
+Clear the queues and return the tracks in flight.
+This kernel runs after all secondary particles were produced.
+
+#### `InitPrimaries` and `InitParticleQueues`
+
+Used to initialize multiple primary particles with separate seeds.