Skip to content

Commit 61e33bc

Browse files
hahnjohahnjo
committed
example19: Use interaction queues instead of SOAData
This revives an old concept from the toy examples 1, 2, and 4 (deleted in commit 9230c4a) to have one queue per discrete interaction. This appears to be equally fast as the quite complex scheduling kernel based on SOAData.
1 parent 1bc7c3d commit 61e33bc

File tree

4 files changed

+300
-365
lines changed

4 files changed

+300
-365
lines changed

‎examples/Example19/electrons.cu‎

Lines changed: 117 additions & 121 deletions
Original file line numberDiff line numberDiff line change
@@ -31,8 +31,8 @@
3131
template <bool IsElectron>
3232
static __device__ __forceinline__ void TransportElectrons(Track *electrons, const adept::MParray *active,
3333
Secondaries &secondaries, adept::MParray *activeQueue,
34-
GlobalScoring *globalScoring,
35-
ScoringPerVolume *scoringPerVolume, SOAData soaData)
34+
Interactions interactions, GlobalScoring *globalScoring,
35+
ScoringPerVolume *scoringPerVolume)
3636
{
3737
#ifdef VECGEOM_FLOAT_PRECISION
3838
const Precision kPush = 10 * vecgeom::kTolerance;
@@ -65,9 +65,6 @@ static __device__ __forceinline__ void TransportElectrons(Track *electrons, cons
6565
if (push) activeQueue->push_back(globalSlot);
6666
};
6767

68-
// Signal that this globalSlot doesn't undergo an interaction (yet)
69-
soaData.nextInteraction[i] = -1;
70-
7168
// Init a track with the needed data to call into G4HepEm.
7269
G4HepEmElectronTrack elTrack;
7370
G4HepEmTrack *theTrack = elTrack.GetTrack();
@@ -291,158 +288,157 @@ static __device__ __forceinline__ void TransportElectrons(Track *electrons, cons
291288
continue;
292289
}
293290

294-
soaData.nextInteraction[i] = winnerProcessIndex;
291+
interactions.queues[winnerProcessIndex]->push_back(globalSlot);
295292

296293
survive(false);
297294
}
298295
}
299296

300297
// Instantiate kernels for electrons and positrons.
301298
__global__ void TransportElectrons(Track *electrons, const adept::MParray *active, Secondaries secondaries,
302-
adept::MParray *activeQueue, GlobalScoring *globalScoring,
303-
ScoringPerVolume *scoringPerVolume, SOAData soaData)
299+
adept::MParray *activeQueue, Interactions interactions, GlobalScoring *globalScoring,
300+
ScoringPerVolume *scoringPerVolume)
304301
{
305-
TransportElectrons</*IsElectron*/ true>(electrons, active, secondaries, activeQueue, globalScoring, scoringPerVolume,
306-
soaData);
302+
TransportElectrons</*IsElectron*/ true>(electrons, active, secondaries, activeQueue, interactions, globalScoring,
303+
scoringPerVolume);
307304
}
308305
__global__ void TransportPositrons(Track *positrons, const adept::MParray *active, Secondaries secondaries,
309-
adept::MParray *activeQueue, GlobalScoring *globalScoring,
310-
ScoringPerVolume *scoringPerVolume, SOAData soaData)
306+
adept::MParray *activeQueue, Interactions interactions, GlobalScoring *globalScoring,
307+
ScoringPerVolume *scoringPerVolume)
311308
{
312-
TransportElectrons</*IsElectron*/ false>(positrons, active, secondaries, activeQueue, globalScoring, scoringPerVolume,
313-
soaData);
309+
TransportElectrons</*IsElectron*/ false>(positrons, active, secondaries, activeQueue, interactions, globalScoring,
310+
scoringPerVolume);
314311
}
315312

316313
template <bool IsElectron, int ProcessIndex>
317-
__device__ void ElectronInteraction(int const globalSlot, SOAData const & /*soaData*/, int const /*soaSlot*/,
318-
Track *particles, Secondaries secondaries, adept::MParray *activeQueue,
319-
GlobalScoring *globalScoring, ScoringPerVolume *scoringPerVolume)
314+
__device__ void ElectronInteraction(Track *particles, const adept::MParray *active, Secondaries secondaries,
315+
adept::MParray *activeQueue, GlobalScoring *globalScoring,
316+
ScoringPerVolume *scoringPerVolume)
320317
{
321-
Track &currentTrack = particles[globalSlot];
322-
auto energy = currentTrack.energy;
323-
const auto pos = currentTrack.pos;
324-
auto dir = currentTrack.dir;
325-
const auto navState = currentTrack.navState;
326-
const auto volume = navState.Top();
327-
// the MCC vector is indexed by the logical volume id
328-
const int lvolID = volume->GetLogicalVolume()->id();
329-
const int theMCIndex = MCIndex[lvolID];
330-
331-
auto survive = [&] {
332-
currentTrack.dir = dir;
333-
currentTrack.energy = energy;
334-
activeQueue->push_back(globalSlot);
335-
};
336-
337-
const double theElCut = g4HepEmData.fTheMatCutData->fMatCutData[theMCIndex].fSecElProdCutE;
338-
339-
RanluxppDouble newRNG{currentTrack.rngState.Branch()};
340-
G4HepEmRandomEngine rnge{&currentTrack.rngState};
341-
342-
if constexpr (ProcessIndex == 0) {
343-
// Invoke ionization (for e-/e+):
344-
double deltaEkin = (IsElectron) ? G4HepEmElectronInteractionIoni::SampleETransferMoller(theElCut, energy, &rnge)
345-
: G4HepEmElectronInteractionIoni::SampleETransferBhabha(theElCut, energy, &rnge);
346-
347-
double dirPrimary[] = {dir.x(), dir.y(), dir.z()};
348-
double dirSecondary[3];
349-
G4HepEmElectronInteractionIoni::SampleDirections(energy, deltaEkin, dirSecondary, dirPrimary, &rnge);
350-
351-
Track &secondary = secondaries.electrons.NextTrack();
352-
atomicAdd(&globalScoring->numElectrons, 1);
353-
354-
secondary.InitAsSecondary(pos, navState);
355-
secondary.rngState = newRNG;
356-
secondary.energy = deltaEkin;
357-
secondary.dir.Set(dirSecondary[0], dirSecondary[1], dirSecondary[2]);
358-
359-
energy -= deltaEkin;
360-
dir.Set(dirPrimary[0], dirPrimary[1], dirPrimary[2]);
361-
survive();
362-
} else if constexpr (ProcessIndex == 1) {
363-
// Invoke model for Bremsstrahlung: either SB- or Rel-Brem.
364-
double logEnergy = std::log(energy);
365-
double deltaEkin = energy < g4HepEmPars.fElectronBremModelLim
366-
? G4HepEmElectronInteractionBrem::SampleETransferSB(&g4HepEmData, energy, logEnergy,
367-
theMCIndex, &rnge, IsElectron)
368-
: G4HepEmElectronInteractionBrem::SampleETransferRB(&g4HepEmData, energy, logEnergy,
369-
theMCIndex, &rnge, IsElectron);
370-
371-
double dirPrimary[] = {dir.x(), dir.y(), dir.z()};
372-
double dirSecondary[3];
373-
G4HepEmElectronInteractionBrem::SampleDirections(energy, deltaEkin, dirSecondary, dirPrimary, &rnge);
374-
375-
Track &gamma = secondaries.gammas.NextTrack();
376-
atomicAdd(&globalScoring->numGammas, 1);
377-
378-
gamma.InitAsSecondary(pos, navState);
379-
gamma.rngState = newRNG;
380-
gamma.energy = deltaEkin;
381-
gamma.dir.Set(dirSecondary[0], dirSecondary[1], dirSecondary[2]);
382-
383-
energy -= deltaEkin;
384-
dir.Set(dirPrimary[0], dirPrimary[1], dirPrimary[2]);
385-
survive();
386-
} else if constexpr (ProcessIndex == 2) {
387-
// Invoke annihilation (in-flight) for e+
388-
double dirPrimary[] = {dir.x(), dir.y(), dir.z()};
389-
double theGamma1Ekin, theGamma2Ekin;
390-
double theGamma1Dir[3], theGamma2Dir[3];
391-
G4HepEmPositronInteractionAnnihilation::SampleEnergyAndDirectionsInFlight(
392-
energy, dirPrimary, &theGamma1Ekin, theGamma1Dir, &theGamma2Ekin, theGamma2Dir, &rnge);
393-
394-
Track &gamma1 = secondaries.gammas.NextTrack();
395-
Track &gamma2 = secondaries.gammas.NextTrack();
396-
atomicAdd(&globalScoring->numGammas, 2);
397-
398-
gamma1.InitAsSecondary(pos, navState);
399-
gamma1.rngState = newRNG;
400-
gamma1.energy = theGamma1Ekin;
401-
gamma1.dir.Set(theGamma1Dir[0], theGamma1Dir[1], theGamma1Dir[2]);
402-
403-
gamma2.InitAsSecondary(pos, navState);
404-
// Reuse the RNG state of the dying track.
405-
gamma2.rngState = currentTrack.rngState;
406-
gamma2.energy = theGamma2Ekin;
407-
gamma2.dir.Set(theGamma2Dir[0], theGamma2Dir[1], theGamma2Dir[2]);
408-
409-
// The current track is killed by not enqueuing into the next activeQueue.
318+
int activeSize = active->size();
319+
for (int i = blockIdx.x * blockDim.x + threadIdx.x; i < activeSize; i += blockDim.x * gridDim.x) {
320+
const int globalSlot = (*active)[i];
321+
Track &currentTrack = particles[globalSlot];
322+
auto energy = currentTrack.energy;
323+
const auto pos = currentTrack.pos;
324+
auto dir = currentTrack.dir;
325+
const auto navState = currentTrack.navState;
326+
const auto volume = navState.Top();
327+
// the MCC vector is indexed by the logical volume id
328+
const int lvolID = volume->GetLogicalVolume()->id();
329+
const int theMCIndex = MCIndex[lvolID];
330+
331+
auto survive = [&] {
332+
currentTrack.dir = dir;
333+
currentTrack.energy = energy;
334+
activeQueue->push_back(globalSlot);
335+
};
336+
337+
const double theElCut = g4HepEmData.fTheMatCutData->fMatCutData[theMCIndex].fSecElProdCutE;
338+
339+
RanluxppDouble newRNG{currentTrack.rngState.Branch()};
340+
G4HepEmRandomEngine rnge{&currentTrack.rngState};
341+
342+
if constexpr (ProcessIndex == 0) {
343+
// Invoke ionization (for e-/e+):
344+
double deltaEkin = (IsElectron) ? G4HepEmElectronInteractionIoni::SampleETransferMoller(theElCut, energy, &rnge)
345+
: G4HepEmElectronInteractionIoni::SampleETransferBhabha(theElCut, energy, &rnge);
346+
347+
double dirPrimary[] = {dir.x(), dir.y(), dir.z()};
348+
double dirSecondary[3];
349+
G4HepEmElectronInteractionIoni::SampleDirections(energy, deltaEkin, dirSecondary, dirPrimary, &rnge);
350+
351+
Track &secondary = secondaries.electrons.NextTrack();
352+
atomicAdd(&globalScoring->numElectrons, 1);
353+
354+
secondary.InitAsSecondary(pos, navState);
355+
secondary.rngState = newRNG;
356+
secondary.energy = deltaEkin;
357+
secondary.dir.Set(dirSecondary[0], dirSecondary[1], dirSecondary[2]);
358+
359+
energy -= deltaEkin;
360+
dir.Set(dirPrimary[0], dirPrimary[1], dirPrimary[2]);
361+
survive();
362+
} else if constexpr (ProcessIndex == 1) {
363+
// Invoke model for Bremsstrahlung: either SB- or Rel-Brem.
364+
double logEnergy = std::log(energy);
365+
double deltaEkin = energy < g4HepEmPars.fElectronBremModelLim
366+
? G4HepEmElectronInteractionBrem::SampleETransferSB(&g4HepEmData, energy, logEnergy,
367+
theMCIndex, &rnge, IsElectron)
368+
: G4HepEmElectronInteractionBrem::SampleETransferRB(&g4HepEmData, energy, logEnergy,
369+
theMCIndex, &rnge, IsElectron);
370+
371+
double dirPrimary[] = {dir.x(), dir.y(), dir.z()};
372+
double dirSecondary[3];
373+
G4HepEmElectronInteractionBrem::SampleDirections(energy, deltaEkin, dirSecondary, dirPrimary, &rnge);
374+
375+
Track &gamma = secondaries.gammas.NextTrack();
376+
atomicAdd(&globalScoring->numGammas, 1);
377+
378+
gamma.InitAsSecondary(pos, navState);
379+
gamma.rngState = newRNG;
380+
gamma.energy = deltaEkin;
381+
gamma.dir.Set(dirSecondary[0], dirSecondary[1], dirSecondary[2]);
382+
383+
energy -= deltaEkin;
384+
dir.Set(dirPrimary[0], dirPrimary[1], dirPrimary[2]);
385+
survive();
386+
} else if constexpr (ProcessIndex == 2) {
387+
// Invoke annihilation (in-flight) for e+
388+
double dirPrimary[] = {dir.x(), dir.y(), dir.z()};
389+
double theGamma1Ekin, theGamma2Ekin;
390+
double theGamma1Dir[3], theGamma2Dir[3];
391+
G4HepEmPositronInteractionAnnihilation::SampleEnergyAndDirectionsInFlight(
392+
energy, dirPrimary, &theGamma1Ekin, theGamma1Dir, &theGamma2Ekin, theGamma2Dir, &rnge);
393+
394+
Track &gamma1 = secondaries.gammas.NextTrack();
395+
Track &gamma2 = secondaries.gammas.NextTrack();
396+
atomicAdd(&globalScoring->numGammas, 2);
397+
398+
gamma1.InitAsSecondary(pos, navState);
399+
gamma1.rngState = newRNG;
400+
gamma1.energy = theGamma1Ekin;
401+
gamma1.dir.Set(theGamma1Dir[0], theGamma1Dir[1], theGamma1Dir[2]);
402+
403+
gamma2.InitAsSecondary(pos, navState);
404+
// Reuse the RNG state of the dying track.
405+
gamma2.rngState = currentTrack.rngState;
406+
gamma2.energy = theGamma2Ekin;
407+
gamma2.dir.Set(theGamma2Dir[0], theGamma2Dir[1], theGamma2Dir[2]);
408+
409+
// The current track is killed by not enqueuing into the next activeQueue.
410+
}
410411
}
411412
}
412413

413414
__global__ void IonizationEl(Track *particles, const adept::MParray *active, Secondaries secondaries,
414415
adept::MParray *activeQueue, GlobalScoring *globalScoring,
415-
ScoringPerVolume *scoringPerVolume, SOAData const soaData)
416+
ScoringPerVolume *scoringPerVolume)
416417
{
417-
InteractionLoop<0>(&ElectronInteraction<true, 0>, active, soaData, particles, secondaries, activeQueue, globalScoring,
418-
scoringPerVolume);
418+
ElectronInteraction<true, 0>(particles, active, secondaries, activeQueue, globalScoring, scoringPerVolume);
419419
}
420420
__global__ void BremsstrahlungEl(Track *particles, const adept::MParray *active, Secondaries secondaries,
421421
adept::MParray *activeQueue, GlobalScoring *globalScoring,
422-
ScoringPerVolume *scoringPerVolume, SOAData const soaData)
422+
ScoringPerVolume *scoringPerVolume)
423423
{
424-
InteractionLoop<1>(&ElectronInteraction<true, 1>, active, soaData, particles, secondaries, activeQueue, globalScoring,
425-
scoringPerVolume);
424+
ElectronInteraction<true, 1>(particles, active, secondaries, activeQueue, globalScoring, scoringPerVolume);
426425
}
427426

428427
__global__ void IonizationPos(Track *particles, const adept::MParray *active, Secondaries secondaries,
429428
adept::MParray *activeQueue, GlobalScoring *globalScoring,
430-
ScoringPerVolume *scoringPerVolume, SOAData const soaData)
429+
ScoringPerVolume *scoringPerVolume)
431430
{
432-
InteractionLoop<0>(&ElectronInteraction<false, 0>, active, soaData, particles, secondaries, activeQueue,
433-
globalScoring, scoringPerVolume);
431+
ElectronInteraction<false, 0>(particles, active, secondaries, activeQueue, globalScoring, scoringPerVolume);
434432
}
435433
__global__ void BremsstrahlungPos(Track *particles, const adept::MParray *active, Secondaries secondaries,
436434
adept::MParray *activeQueue, GlobalScoring *globalScoring,
437-
ScoringPerVolume *scoringPerVolume, SOAData const soaData)
435+
ScoringPerVolume *scoringPerVolume)
438436
{
439-
InteractionLoop<1>(&ElectronInteraction<false, 1>, active, soaData, particles, secondaries, activeQueue,
440-
globalScoring, scoringPerVolume);
437+
ElectronInteraction<false, 1>(particles, active, secondaries, activeQueue, globalScoring, scoringPerVolume);
441438
}
442439
__global__ void AnnihilationPos(Track *particles, const adept::MParray *active, Secondaries secondaries,
443440
adept::MParray *activeQueue, GlobalScoring *globalScoring,
444-
ScoringPerVolume *scoringPerVolume, SOAData const soaData)
441+
ScoringPerVolume *scoringPerVolume)
445442
{
446-
InteractionLoop<2>(&ElectronInteraction<false, 2>, active, soaData, particles, secondaries, activeQueue,
447-
globalScoring, scoringPerVolume);
443+
ElectronInteraction<false, 2>(particles, active, secondaries, activeQueue, globalScoring, scoringPerVolume);
448444
}

0 commit comments

Comments
 (0)