Skip to content

Commit 822af79

Browse files
bernhardmgruberbernhardmgruber
committed
Keep some Track members in registers
The track variables are loaded at the beginning of the kernel and written back when the track survives. The signature of InitAsSecondary() had to be changed to not load from the track again.
1 parent f18b1b2 commit 822af79

File tree

3 files changed

+71
-53
lines changed

3 files changed

+71
-53
lines changed

‎examples/Example19/electrons.cu‎

Lines changed: 39 additions & 31 deletions
Original file line numberDiff line numberDiff line change
@@ -47,13 +47,21 @@ static __device__ __forceinline__ void TransportElectrons(Track *electrons, cons
4747
for (int i = blockIdx.x * blockDim.x + threadIdx.x; i < activeSize; i += blockDim.x * gridDim.x) {
4848
const int globalSlot = (*active)[i];
4949
Track &currentTrack = electrons[globalSlot];
50-
const auto volume = currentTrack.navState.Top();
50+
auto energy = currentTrack.energy;
51+
auto pos = currentTrack.pos;
52+
auto dir = currentTrack.dir;
53+
auto navState = currentTrack.navState;
54+
const auto volume = navState.Top();
5155
const int volumeID = volume->id();
5256
// the MCC vector is indexed by the logical volume id
5357
const int lvolID = volume->GetLogicalVolume()->id();
5458
const int theMCIndex = MCIndex[lvolID];
5559

5660
auto survive = [&](bool push = true) {
61+
currentTrack.energy = energy;
62+
currentTrack.pos = pos;
63+
currentTrack.dir = dir;
64+
currentTrack.navState = navState;
5765
if (push) activeQueue->push_back(globalSlot);
5866
};
5967

@@ -63,9 +71,9 @@ static __device__ __forceinline__ void TransportElectrons(Track *electrons, cons
6371
// Init a track with the needed data to call into G4HepEm.
6472
G4HepEmElectronTrack elTrack;
6573
G4HepEmTrack *theTrack = elTrack.GetTrack();
66-
theTrack->SetEKin(currentTrack.energy);
74+
theTrack->SetEKin(energy);
6775
theTrack->SetMCIndex(theMCIndex);
68-
theTrack->SetOnBoundary(currentTrack.navState.IsOnBoundary());
76+
theTrack->SetOnBoundary(navState.IsOnBoundary());
6977
theTrack->SetCharge(Charge);
7078
G4HepEmMSCTrackData *mscData = elTrack.GetMSCTrackData();
7179
mscData->fIsFirstStep = currentTrack.initialRange < 0;
@@ -80,8 +88,8 @@ static __device__ __forceinline__ void TransportElectrons(Track *electrons, cons
8088

8189
// Compute safety, needed for MSC step limit.
8290
double safety = 0;
83-
if (!currentTrack.navState.IsOnBoundary()) {
84-
safety = BVHNavigator::ComputeSafety(currentTrack.pos, currentTrack.navState);
91+
if (!navState.IsOnBoundary()) {
92+
safety = BVHNavigator::ComputeSafety(pos, navState);
8593
}
8694
theTrack->SetSafety(safety);
8795

@@ -102,9 +110,9 @@ static __device__ __forceinline__ void TransportElectrons(Track *electrons, cons
102110

103111
bool restrictedPhysicalStepLength = false;
104112
if (BzFieldValue != 0) {
105-
const double momentumMag = sqrt(currentTrack.energy * (currentTrack.energy + 2.0 * Mass));
113+
const double momentumMag = sqrt(energy * (energy + 2.0 * Mass));
106114
// Distance along the track direction to reach the maximum allowed error
107-
const double safeLength = fieldPropagatorBz.ComputeSafeLength(momentumMag, Charge, currentTrack.dir);
115+
const double safeLength = fieldPropagatorBz.ComputeSafeLength(momentumMag, Charge, dir);
108116

109117
constexpr int MaxSafeLength = 10;
110118
double limit = MaxSafeLength * safeLength;
@@ -145,22 +153,20 @@ static __device__ __forceinline__ void TransportElectrons(Track *electrons, cons
145153
vecgeom::NavStateIndex nextState;
146154
if (BzFieldValue != 0) {
147155
geometryStepLength = fieldPropagatorBz.ComputeStepAndNextVolume<BVHNavigator>(
148-
currentTrack.energy, Mass, Charge, geometricalStepLengthFromPhysics, currentTrack.pos, currentTrack.dir,
149-
currentTrack.navState, nextState, propagated, safety);
156+
energy, Mass, Charge, geometricalStepLengthFromPhysics, pos, dir, navState, nextState, propagated, safety);
150157
} else {
151-
geometryStepLength =
152-
BVHNavigator::ComputeStepAndNextVolume(currentTrack.pos, currentTrack.dir, geometricalStepLengthFromPhysics,
153-
currentTrack.navState, nextState, kPush);
154-
currentTrack.pos += geometryStepLength * currentTrack.dir;
158+
geometryStepLength = BVHNavigator::ComputeStepAndNextVolume(pos, dir, geometricalStepLengthFromPhysics, navState,
159+
nextState, kPush);
160+
pos += geometryStepLength * dir;
155161
}
156162

157163
// Set boundary state in navState so the next step and secondaries get the
158-
// correct information (currentTrack.navState = nextState only if relocated
164+
// correct information (navState = nextState only if relocated
159165
// in case of a boundary; see below)
160-
currentTrack.navState.SetBoundaryState(nextState.IsOnBoundary());
166+
navState.SetBoundaryState(nextState.IsOnBoundary());
161167

162168
// Propagate information from geometrical step to MSC.
163-
theTrack->SetDirection(currentTrack.dir.x(), currentTrack.dir.y(), currentTrack.dir.z());
169+
theTrack->SetDirection(dir.x(), dir.y(), dir.z());
164170
theTrack->SetGStepLength(geometryStepLength);
165171
theTrack->SetOnBoundary(nextState.IsOnBoundary());
166172

@@ -169,7 +175,7 @@ static __device__ __forceinline__ void TransportElectrons(Track *electrons, cons
169175

170176
// Collect the direction change and displacement by MSC.
171177
const double *direction = theTrack->GetDirection();
172-
currentTrack.dir.Set(direction[0], direction[1], direction[2]);
178+
dir.Set(direction[0], direction[1], direction[2]);
173179
if (!nextState.IsOnBoundary()) {
174180
const double *mscDisplacement = mscData->GetDisplacement();
175181
vecgeom::Vector3D<Precision> displacement(mscDisplacement[0], mscDisplacement[1], mscDisplacement[2]);
@@ -186,18 +192,18 @@ static __device__ __forceinline__ void TransportElectrons(Track *electrons, cons
186192
// Apply displacement, depending on how close we are to a boundary.
187193
// 1a. Far away from geometry boundary:
188194
if (reducedSafety > 0.0 && dispR <= reducedSafety) {
189-
currentTrack.pos += displacement;
195+
pos += displacement;
190196
} else {
191197
// Recompute safety.
192-
safety = BVHNavigator::ComputeSafety(currentTrack.pos, currentTrack.navState);
198+
safety = BVHNavigator::ComputeSafety(pos, navState);
193199
reducedSafety = sFact * safety;
194200

195201
// 1b. Far away from geometry boundary:
196202
if (reducedSafety > 0.0 && dispR <= reducedSafety) {
197-
currentTrack.pos += displacement;
203+
pos += displacement;
198204
// 2. Push to boundary:
199205
} else if (reducedSafety > kGeomMinLength) {
200-
currentTrack.pos += displacement * (reducedSafety / dispR);
206+
pos += displacement * (reducedSafety / dispR);
201207
}
202208
// 3. Very small safety: do nothing.
203209
}
@@ -209,7 +215,7 @@ static __device__ __forceinline__ void TransportElectrons(Track *electrons, cons
209215
atomicAdd(&scoringPerVolume->chargedTrackLength[volumeID], elTrack.GetPStepLength());
210216

211217
// Collect the changes in energy and deposit.
212-
currentTrack.energy = theTrack->GetEKin();
218+
energy = theTrack->GetEKin();
213219
double energyDeposit = theTrack->GetEnergyDeposit();
214220
atomicAdd(&globalScoring->energyDeposit, energyDeposit);
215221
atomicAdd(&scoringPerVolume->energyDeposit[volumeID], energyDeposit);
@@ -234,13 +240,13 @@ static __device__ __forceinline__ void TransportElectrons(Track *electrons, cons
234240
double sinPhi, cosPhi;
235241
sincos(phi, &sinPhi, &cosPhi);
236242

237-
gamma1.InitAsSecondary(/*parent=*/currentTrack);
243+
gamma1.InitAsSecondary(pos, navState);
238244
newRNG.Advance();
239245
gamma1.rngState = newRNG;
240246
gamma1.energy = copcore::units::kElectronMassC2;
241247
gamma1.dir.Set(sint * cosPhi, sint * sinPhi, cost);
242248

243-
gamma2.InitAsSecondary(/*parent=*/currentTrack);
249+
gamma2.InitAsSecondary(pos, navState);
244250
// Reuse the RNG state of the dying track.
245251
gamma2.rngState = currentTrack.rngState;
246252
gamma2.energy = copcore::units::kElectronMassC2;
@@ -256,10 +262,10 @@ static __device__ __forceinline__ void TransportElectrons(Track *electrons, cons
256262

257263
// Kill the particle if it left the world.
258264
if (nextState.Top() != nullptr) {
259-
BVHNavigator::RelocateToNextVolume(currentTrack.pos, currentTrack.dir, nextState);
265+
BVHNavigator::RelocateToNextVolume(pos, dir, nextState);
260266

261267
// Move to the next boundary.
262-
currentTrack.navState = nextState;
268+
navState = nextState;
263269
survive();
264270
}
265271
continue;
@@ -313,7 +319,9 @@ __device__ void ElectronInteraction(int const globalSlot, SOAData const & /*soaD
313319
GlobalScoring *globalScoring, ScoringPerVolume *scoringPerVolume)
314320
{
315321
Track &currentTrack = particles[globalSlot];
316-
const auto volume = currentTrack.navState.Top();
322+
const auto pos = currentTrack.pos;
323+
const auto navState = currentTrack.navState;
324+
const auto volume = navState.Top();
317325
// the MCC vector is indexed by the logical volume id
318326
const int lvolID = volume->GetLogicalVolume()->id();
319327
const int theMCIndex = MCIndex[lvolID];
@@ -338,7 +346,7 @@ __device__ void ElectronInteraction(int const globalSlot, SOAData const & /*soaD
338346
Track &secondary = secondaries.electrons.NextTrack();
339347
atomicAdd(&globalScoring->numElectrons, 1);
340348

341-
secondary.InitAsSecondary(/*parent=*/currentTrack);
349+
secondary.InitAsSecondary(pos, navState);
342350
secondary.rngState = newRNG;
343351
secondary.energy = deltaEkin;
344352
secondary.dir.Set(dirSecondary[0], dirSecondary[1], dirSecondary[2]);
@@ -362,7 +370,7 @@ __device__ void ElectronInteraction(int const globalSlot, SOAData const & /*soaD
362370
Track &gamma = secondaries.gammas.NextTrack();
363371
atomicAdd(&globalScoring->numGammas, 1);
364372

365-
gamma.InitAsSecondary(/*parent=*/currentTrack);
373+
gamma.InitAsSecondary(pos, navState);
366374
gamma.rngState = newRNG;
367375
gamma.energy = deltaEkin;
368376
gamma.dir.Set(dirSecondary[0], dirSecondary[1], dirSecondary[2]);
@@ -382,12 +390,12 @@ __device__ void ElectronInteraction(int const globalSlot, SOAData const & /*soaD
382390
Track &gamma2 = secondaries.gammas.NextTrack();
383391
atomicAdd(&globalScoring->numGammas, 2);
384392

385-
gamma1.InitAsSecondary(/*parent=*/currentTrack);
393+
gamma1.InitAsSecondary(pos, navState);
386394
gamma1.rngState = newRNG;
387395
gamma1.energy = theGamma1Ekin;
388396
gamma1.dir.Set(theGamma1Dir[0], theGamma1Dir[1], theGamma1Dir[2]);
389397

390-
gamma2.InitAsSecondary(/*parent=*/currentTrack);
398+
gamma2.InitAsSecondary(pos, navState);
391399
// Reuse the RNG state of the dying track.
392400
gamma2.rngState = currentTrack.rngState;
393401
gamma2.energy = theGamma2Ekin;

‎examples/Example19/example.cuh‎

Lines changed: 4 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -36,7 +36,8 @@ struct Track {
3636

3737
__host__ __device__ double Uniform() { return rngState.Rndm(); }
3838

39-
__host__ __device__ void InitAsSecondary(const Track &parent)
39+
__host__ __device__ void InitAsSecondary(const vecgeom::Vector3D<Precision> &parentPos,
40+
const vecgeom::NavStateIndex &parentNavState)
4041
{
4142
// The caller is responsible to branch a new RNG state and to set the energy.
4243
this->numIALeft[0] = -1.0;
@@ -49,8 +50,8 @@ struct Track {
4950

5051
// A secondary inherits the position of its parent; the caller is responsible
5152
// to update the directions.
52-
this->pos = parent.pos;
53-
this->navState = parent.navState;
53+
this->pos = parentPos;
54+
this->navState = parentNavState;
5455
}
5556
};
5657

‎examples/Example19/gammas.cu‎

Lines changed: 28 additions & 19 deletions
Original file line numberDiff line numberDiff line change
@@ -32,13 +32,19 @@ __global__ void TransportGammas(Track *gammas, const adept::MParray *active, Sec
3232
for (int i = blockIdx.x * blockDim.x + threadIdx.x; i < activeSize; i += blockDim.x * gridDim.x) {
3333
const int slot = (*active)[i];
3434
Track &currentTrack = gammas[slot];
35-
const auto volume = currentTrack.navState.Top();
35+
const auto energy = currentTrack.energy;
36+
auto pos = currentTrack.pos;
37+
const auto dir = currentTrack.dir;
38+
auto navState = currentTrack.navState;
39+
const auto volume = navState.Top();
3640
const int volumeID = volume->id();
3741
// the MCC vector is indexed by the logical volume id
3842
const int lvolID = volume->GetLogicalVolume()->id();
3943
const int theMCIndex = MCIndex[lvolID];
4044

4145
auto survive = [&](bool push = true) {
46+
currentTrack.pos = pos;
47+
currentTrack.navState = navState;
4248
if (push) activeQueue->push_back(slot);
4349
};
4450

@@ -48,7 +54,7 @@ __global__ void TransportGammas(Track *gammas, const adept::MParray *active, Sec
4854
// Init a track with the needed data to call into G4HepEm.
4955
G4HepEmGammaTrack gammaTrack;
5056
G4HepEmTrack *theTrack = gammaTrack.GetTrack();
51-
theTrack->SetEKin(currentTrack.energy);
57+
theTrack->SetEKin(energy);
5258
theTrack->SetMCIndex(theMCIndex);
5359

5460
// Sample the `number-of-interaction-left` and put it into the track.
@@ -71,15 +77,15 @@ __global__ void TransportGammas(Track *gammas, const adept::MParray *active, Sec
7177

7278
// Check if there's a volume boundary in between.
7379
vecgeom::NavStateIndex nextState;
74-
double geometryStepLength = BVHNavigator::ComputeStepAndNextVolume(
75-
currentTrack.pos, currentTrack.dir, geometricalStepLengthFromPhysics, currentTrack.navState, nextState, kPush);
76-
currentTrack.pos += geometryStepLength * currentTrack.dir;
80+
double geometryStepLength =
81+
BVHNavigator::ComputeStepAndNextVolume(pos, dir, geometricalStepLengthFromPhysics, navState, nextState, kPush);
82+
pos += geometryStepLength * dir;
7783
atomicAdd(&globalScoring->neutralSteps, 1);
7884

7985
// Set boundary state in navState so the next step and secondaries get the
80-
// correct information (currentTrack.navState = nextState only if relocated
86+
// correct information (navState = nextState only if relocated
8187
// in case of a boundary; see below)
82-
currentTrack.navState.SetBoundaryState(nextState.IsOnBoundary());
88+
navState.SetBoundaryState(nextState.IsOnBoundary());
8389

8490
// Propagate information from geometrical step to G4HepEm.
8591
theTrack->SetGStepLength(geometryStepLength);
@@ -99,10 +105,10 @@ __global__ void TransportGammas(Track *gammas, const adept::MParray *active, Sec
99105

100106
// Kill the particle if it left the world.
101107
if (nextState.Top() != nullptr) {
102-
BVHNavigator::RelocateToNextVolume(currentTrack.pos, currentTrack.dir, nextState);
108+
BVHNavigator::RelocateToNextVolume(pos, dir, nextState);
103109

104110
// Move to the next boundary.
105-
currentTrack.navState = nextState;
111+
navState = nextState;
106112
survive();
107113
}
108114
continue;
@@ -128,12 +134,15 @@ __device__ void GammaInteraction(int const globalSlot, SOAData const &soaData, i
128134
ScoringPerVolume *scoringPerVolume)
129135
{
130136
Track &currentTrack = particles[globalSlot];
131-
const auto volume = currentTrack.navState.Top();
137+
const auto energy = currentTrack.energy;
138+
const auto pos = currentTrack.pos;
139+
const auto dir = currentTrack.dir;
140+
const auto navState = currentTrack.navState;
141+
const auto volume = navState.Top();
132142
const int volumeID = volume->id();
133143
// the MCC vector is indexed by the logical volume id
134144
const int lvolID = volume->GetLogicalVolume()->id();
135145
const int theMCIndex = MCIndex[lvolID];
136-
const auto energy = currentTrack.energy;
137146

138147
auto survive = [&] { activeQueue->push_back(globalSlot); };
139148

@@ -152,7 +161,7 @@ __device__ void GammaInteraction(int const globalSlot, SOAData const &soaData, i
152161
G4HepEmGammaInteractionConversion::SampleKinEnergies(&g4HepEmData, energy, logEnergy, theMCIndex, elKinEnergy,
153162
posKinEnergy, &rnge);
154163

155-
double dirPrimary[] = {currentTrack.dir.x(), currentTrack.dir.y(), currentTrack.dir.z()};
164+
double dirPrimary[] = {dir.x(), dir.y(), dir.z()};
156165
double dirSecondaryEl[3], dirSecondaryPos[3];
157166
G4HepEmGammaInteractionConversion::SampleDirections(dirPrimary, dirSecondaryEl, dirSecondaryPos, elKinEnergy,
158167
posKinEnergy, &rnge);
@@ -162,12 +171,12 @@ __device__ void GammaInteraction(int const globalSlot, SOAData const &soaData, i
162171
atomicAdd(&globalScoring->numElectrons, 1);
163172
atomicAdd(&globalScoring->numPositrons, 1);
164173

165-
electron.InitAsSecondary(/*parent=*/currentTrack);
174+
electron.InitAsSecondary(pos, navState);
166175
electron.rngState = newRNG;
167176
electron.energy = elKinEnergy;
168177
electron.dir.Set(dirSecondaryEl[0], dirSecondaryEl[1], dirSecondaryEl[2]);
169178

170-
positron.InitAsSecondary(/*parent=*/currentTrack);
179+
positron.InitAsSecondary(pos, navState);
171180
// Reuse the RNG state of the dying track.
172181
positron.rngState = currentTrack.rngState;
173182
positron.energy = posKinEnergy;
@@ -181,7 +190,7 @@ __device__ void GammaInteraction(int const globalSlot, SOAData const &soaData, i
181190
survive();
182191
return;
183192
}
184-
const double origDirPrimary[] = {currentTrack.dir.x(), currentTrack.dir.y(), currentTrack.dir.z()};
193+
const double origDirPrimary[] = {dir.x(), dir.y(), dir.z()};
185194
double dirPrimary[3];
186195
const double newEnergyGamma =
187196
G4HepEmGammaInteractionCompton::SamplePhotonEnergyAndDirection(energy, dirPrimary, origDirPrimary, &rnge);
@@ -193,10 +202,10 @@ __device__ void GammaInteraction(int const globalSlot, SOAData const &soaData, i
193202
Track &electron = secondaries.electrons.NextTrack();
194203
atomicAdd(&globalScoring->numElectrons, 1);
195204

196-
electron.InitAsSecondary(/*parent=*/currentTrack);
205+
electron.InitAsSecondary(pos, navState);
197206
electron.rngState = newRNG;
198207
electron.energy = energyEl;
199-
electron.dir = energy * currentTrack.dir - newEnergyGamma * newDirGamma;
208+
electron.dir = energy * dir - newEnergyGamma * newDirGamma;
200209
electron.dir.Normalize();
201210
} else {
202211
atomicAdd(&globalScoring->energyDeposit, energyEl);
@@ -227,11 +236,11 @@ __device__ void GammaInteraction(int const globalSlot, SOAData const &soaData, i
227236
Track &electron = secondaries.electrons.NextTrack();
228237
atomicAdd(&globalScoring->numElectrons, 1);
229238

230-
double dirGamma[] = {currentTrack.dir.x(), currentTrack.dir.y(), currentTrack.dir.z()};
239+
double dirGamma[] = {dir.x(), dir.y(), dir.z()};
231240
double dirPhotoElec[3];
232241
G4HepEmGammaInteractionPhotoelectric::SamplePhotoElectronDirection(photoElecE, dirGamma, dirPhotoElec, &rnge);
233242

234-
electron.InitAsSecondary(/*parent=*/currentTrack);
243+
electron.InitAsSecondary(pos, navState);
235244
electron.rngState = newRNG;
236245
electron.energy = photoElecE;
237246
electron.dir.Set(dirPhotoElec[0], dirPhotoElec[1], dirPhotoElec[2]);

0 commit comments

Comments
 (0)