Skip to content

overlap safety call in MSC with relocation #688

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
SeverinDiederichs wants to merge 1 commit into
base: master
Choose a base branch
from
Open

overlap safety call in MSC with relocation #688

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
55 changes: 37 additions & 18 deletions include/AdePT/transport/AdePTTransport.cuh
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -125,6 +125,9 @@ __global__ void InitParticleQueues(ParticleQueues queues, size_t CapacityTranspo
adept::MParray::MakeInstanceAt(CapacityTransport, queues.nextActive);
#ifdef ADEPT_USE_SPLIT_KERNELS
adept::MParray::MakeInstanceAt(CapacityTransport, queues.propagation);
if (queues.setupInteractions) {
adept::MParray::MakeInstanceAt(CapacityTransport, queues.setupInteractions);
}
for (int i = 0; i < ParticleQueues::numSplitQueues; i++) {
adept::MParray::MakeInstanceAt(CapacityTransport, queues.splitQueues[i]);
}
Expand Down Expand Up @@ -254,6 +257,9 @@ __global__ void FinishIteration(AllParticleQueues all, Stats *stats, TracksAndSl
all.queues[i].initiallyActive->clear();
#ifdef ADEPT_USE_SPLIT_KERNELS
all.queues[i].propagation->clear();
if (all.queues[i].setupInteractions) {
all.queues[i].setupInteractions->clear();
}
for (int j = 0; j < ParticleQueues::numSplitQueues; j++) {
all.queues[i].splitQueues[j]->clear();
}
Expand Down Expand Up @@ -394,6 +400,9 @@ __global__ void ClearAllQueues(AllParticleQueues all)
if (i == GPUQueueIndex::GammaWDT) return;
#ifdef ADEPT_USE_SPLIT_KERNELS
all.queues[i].propagation->clear();
if (all.queues[i].setupInteractions) {
all.queues[i].setupInteractions->clear();
}
for (int j = 0; j < ParticleQueues::numSplitQueues; j++) {
all.queues[i].splitQueues[j]->clear();
}
Expand Down Expand Up @@ -664,7 +673,12 @@ std::unique_ptr<GPUstate, GPUstateDeleter> InitializeGPU(int trackCapacity, int
particleType.queues.nextActive = static_cast<adept::MParray *>(gpuPtr);
#ifdef ADEPT_USE_SPLIT_KERNELS
gpuMalloc(gpuPtr, sizeOfQueueStorage);
particleType.queues.propagation = static_cast<adept::MParray *>(gpuPtr);
particleType.queues.propagation = static_cast<adept::MParray *>(gpuPtr);
particleType.queues.setupInteractions = nullptr;
if (particleIndex != GPUQueueIndex::Gamma) {
gpuMalloc(gpuPtr, sizeOfQueueStorage);
particleType.queues.setupInteractions = static_cast<adept::MParray *>(gpuPtr);
}
for (int j = 0; j < ParticleQueues::numSplitQueues; j++) {
gpuMalloc(gpuPtr, sizeOfQueueStorage);
particleType.queues.splitQueues[j] = static_cast<adept::MParray *>(gpuPtr);
Expand Down Expand Up @@ -898,21 +912,24 @@ void TransportLoop(int trackCapacity, int stepCapacity, int numThreads, TrackBuf
};
const AllParticleQueues allParticleQueues = {{electrons.queues, positrons.queues, gammas.queues, woodcockQueues}};
#ifdef ADEPT_USE_SPLIT_KERNELS
const SplitQueues gammaSplitQueues = {{gammas.queues.splitQueues[ParticleQueues::gammaConversion],
gammas.queues.splitQueues[ParticleQueues::gammaCompton],
gammas.queues.splitQueues[ParticleQueues::gammaPhotoelectric],
gammas.queues.splitQueues[ParticleQueues::gammaWoodcock],
gammas.queues.splitQueues[ParticleQueues::relocation]}};
const SplitQueues gammaSplitQueues = {{gammas.queues.splitQueues[ParticleQueues::gammaConversion],
gammas.queues.splitQueues[ParticleQueues::gammaCompton],
gammas.queues.splitQueues[ParticleQueues::gammaPhotoelectric],
gammas.queues.splitQueues[ParticleQueues::gammaWoodcock],
gammas.queues.splitQueues[ParticleQueues::relocation]}};
// Slots left as nullptr are never indexed: ElectronSetupInteractions only dereferences
// queues[winnerProcessIndex] for winnerProcessIndex in [0, 2] (3 is filtered out as it
// continues on the host) plus positronStoppedAnnihilation. The relocation slot is unused
// for charged particles, whose relocation queue is filled by ElectronMSC and passed to
// ElectronRelocation directly.
const SplitQueues electronSplitQueues = {{electrons.queues.splitQueues[ParticleQueues::chargedIonization],
electrons.queues.splitQueues[ParticleQueues::chargedBremsstrahlung],
nullptr, nullptr,
electrons.queues.splitQueues[ParticleQueues::relocation]}};
nullptr, nullptr, nullptr}};
const SplitQueues positronSplitQueues = {
{positrons.queues.splitQueues[ParticleQueues::chargedIonization],
positrons.queues.splitQueues[ParticleQueues::chargedBremsstrahlung],
positrons.queues.splitQueues[ParticleQueues::positronAnnihilation],
positrons.queues.splitQueues[ParticleQueues::positronStoppedAnnihilation],
positrons.queues.splitQueues[ParticleQueues::relocation]}};
positrons.queues.splitQueues[ParticleQueues::positronStoppedAnnihilation], nullptr}};
#endif
const TracksAndSlots tracksAndSlots = {electrons.tracks,
positrons.tracks,
Expand Down Expand Up @@ -1033,15 +1050,16 @@ void TransportLoop(int trackCapacity, int stepCapacity, int numThreads, TrackBuf
ElectronPropagation<true><<<blocks, threads, 0, electrons.stream>>>(
electrons.tracks, gpuState.hepEmBuffers_d.electronsHepEm, electrons.queues.propagation);
ElectronMSC<true><<<blocks, threads, 0, electrons.stream>>>(
electrons.tracks, gpuState.hepEmBuffers_d.electronsHepEm, electrons.queues.propagation);
ElectronSetupInteractions<true><<<blocks, threads, 0, electrons.stream>>>(
gpuState.hepEmBuffers_d.electronsHepEm, electrons.queues.propagation, particleManager, electronSplitQueues,
kernelOptions);
electrons.tracks, gpuState.hepEmBuffers_d.electronsHepEm, electrons.queues.propagation,
electrons.queues.setupInteractions, electrons.queues.splitQueues[ParticleQueues::relocation]);
ADEPT_DEVICE_API_CALL(EventRecord(electrons.setupEvent, electrons.stream));
ADEPT_DEVICE_API_CALL(StreamWaitEvent(electrons.auxiliaryStream, electrons.setupEvent, 0));
ElectronRelocation<true><<<blocks, threads, 0, electrons.stream>>>(
gpuState.hepEmBuffers_d.electronsHepEm, particleManager,
electrons.queues.splitQueues[ParticleQueues::relocation], kernelOptions);
ElectronSetupInteractions<true><<<blocks, threads, 0, electrons.auxiliaryStream>>>(
gpuState.hepEmBuffers_d.electronsHepEm, electrons.queues.setupInteractions, particleManager,
electronSplitQueues, kernelOptions);
ElectronIonization<true><<<blocks, threads, 0, electrons.auxiliaryStream>>>(
gpuState.hepEmBuffers_d.electronsHepEm, particleManager,
electrons.queues.splitQueues[ParticleQueues::chargedIonization], kernelOptions);
Expand Down Expand Up @@ -1075,15 +1093,16 @@ void TransportLoop(int trackCapacity, int stepCapacity, int numThreads, TrackBuf
ElectronPropagation<false><<<blocks, threads, 0, positrons.stream>>>(
positrons.tracks, gpuState.hepEmBuffers_d.positronsHepEm, positrons.queues.propagation);
ElectronMSC<false><<<blocks, threads, 0, positrons.stream>>>(
positrons.tracks, gpuState.hepEmBuffers_d.positronsHepEm, positrons.queues.propagation);
ElectronSetupInteractions<false><<<blocks, threads, 0, positrons.stream>>>(
gpuState.hepEmBuffers_d.positronsHepEm, positrons.queues.propagation, particleManager, positronSplitQueues,
kernelOptions);
positrons.tracks, gpuState.hepEmBuffers_d.positronsHepEm, positrons.queues.propagation,
positrons.queues.setupInteractions, positrons.queues.splitQueues[ParticleQueues::relocation]);
ADEPT_DEVICE_API_CALL(EventRecord(positrons.setupEvent, positrons.stream));
ADEPT_DEVICE_API_CALL(StreamWaitEvent(positrons.auxiliaryStream, positrons.setupEvent, 0));
ElectronRelocation<false><<<blocks, threads, 0, positrons.stream>>>(
gpuState.hepEmBuffers_d.positronsHepEm, particleManager,
positrons.queues.splitQueues[ParticleQueues::relocation], kernelOptions);
ElectronSetupInteractions<false><<<blocks, threads, 0, positrons.auxiliaryStream>>>(
gpuState.hepEmBuffers_d.positronsHepEm, positrons.queues.setupInteractions, particleManager,
positronSplitQueues, kernelOptions);
ElectronIonization<false><<<blocks, threads, 0, positrons.auxiliaryStream>>>(
gpuState.hepEmBuffers_d.positronsHepEm, particleManager,
positrons.queues.splitQueues[ParticleQueues::chargedIonization], kernelOptions);
Expand Down
112 changes: 57 additions & 55 deletions include/AdePT/transport/kernels/electrons_split.cuh
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -360,7 +360,8 @@ __global__ void ElectronPropagation(ChargedTrack *electronsOrPositrons, G4HepEmE
}

template <bool IsElectron>
__global__ void ElectronMSC(ChargedTrack *electrons, G4HepEmElectronTrack *hepEMTracks, const adept::MParray *active)
__global__ void ElectronMSC(ChargedTrack *electrons, G4HepEmElectronTrack *hepEMTracks, const adept::MParray *active,
adept::MParray *setupQueue, adept::MParray *relocatingQueue)
{
constexpr double restMass = copcore::units::kElectronMassC2;

Expand All @@ -369,57 +370,19 @@ __global__ void ElectronMSC(ChargedTrack *electrons, G4HepEmElectronTrack *hepEM
const int slot = (*active)[i];

ChargedTrack &currentTrack = electrons[slot];
// the MCC vector is indexed by the logical volume id
const int lvolID = currentTrack.navState.GetLogicalId();

// Retrieve HepEM track
G4HepEmElectronTrack &elTrack = hepEMTracks[slot];
G4HepEmTrack *theTrack = elTrack.GetTrack();

G4HepEmMSCTrackData *mscData = elTrack.GetMSCTrackData();
G4HepEmRandomEngine rnge(&currentTrack.rngState);

// Apply continuous effects.
currentTrack.stopped = G4HepEmElectronManager::PerformContinuous(&g4HepEmData, &g4HepEmPars, &elTrack, &rnge);

// Collect the direction change and displacement by MSC.
// Collect the direction change by MSC.
const double *direction = theTrack->GetDirection();
currentTrack.dir.Set(direction[0], direction[1], direction[2]);
if (!currentTrack.nextState.IsOnBoundary()) {
const double *mscDisplacement = mscData->GetDisplacement();
vecgeom::Vector3D<double> displacement(mscDisplacement[0], mscDisplacement[1], mscDisplacement[2]);
const double dLength2 = displacement.Length2();
constexpr double kGeomMinLength = 0.05 * copcore::units::nm; // 0.05 [nm]
constexpr double kGeomMinLength2 = kGeomMinLength * kGeomMinLength; // (0.05 [nm])^2
if (dLength2 > kGeomMinLength2) {
const double dispR = std::sqrt(dLength2);
// Estimate safety by subtracting the geometrical step length.
double safety = currentTrack.GetSafety(currentTrack.pos);
constexpr double sFact = 0.99;
double reducedSafety = sFact * safety;

// Apply displacement, depending on how close we are to a boundary.
// 1a. Far away from geometry boundary:
if (reducedSafety > 0.0 && dispR <= reducedSafety) {
currentTrack.pos += displacement;
} else {
// Recompute safety.
// Use maximum accuracy only if safety is smaller than physicalStepLength
safety = AdePTNavigator::ComputeSafety(currentTrack.pos, currentTrack.navState, dispR);
currentTrack.SetSafety(currentTrack.pos, safety);
reducedSafety = sFact * safety;

// 1b. Far away from geometry boundary:
if (reducedSafety > 0.0 && dispR <= reducedSafety) {
currentTrack.pos += displacement;
// 2. Push to boundary:
} else if (reducedSafety > kGeomMinLength) {
currentTrack.pos += displacement * (reducedSafety / dispR);
}
// 3. Very small safety: do nothing.
}
}
}

// Collect the charged step length (might be changed by MSC). Collect the changes in energy and deposit.
currentTrack.eKin = theTrack->GetEKin();
Expand All @@ -432,14 +395,24 @@ __global__ void ElectronMSC(ChargedTrack *electrons, G4HepEmElectronTrack *hepEM
currentTrack.globalTime += deltaTime;
currentTrack.localTime += deltaTime;
currentTrack.properTime += deltaTime * (restMass / (restMass + currentTrack.eKin));

// Push particles to the appropriate queue for the next kernel, either
// ElectronRelocation or ElectronSetupInteractions.
// Particles on boundary need to be relocated, all others undergo the setup for interactions,
// which includes the displacement of the particle by MSC.
if (currentTrack.nextState.IsOnBoundary() && !currentTrack.stopped) {
relocatingQueue->push_back(slot);
} else {
setupQueue->push_back(slot);
}
}
}

/***
* @brief Adds tracks to interaction and relocation queues depending on their state
* @brief Applies MSC displacement and adds tracks to interaction queues depending on their state
*/
template <bool IsElectron>
__global__ void ElectronSetupInteractions(G4HepEmElectronTrack *hepEMTracks, const adept::MParray *propagationQueue,
__global__ void ElectronSetupInteractions(G4HepEmElectronTrack *hepEMTracks, const adept::MParray *setupQueue,
ParticleManager particleManager, SplitQueues splitQueues,
const TransportKernelOptions options)
{
Expand All @@ -449,34 +422,63 @@ __global__ void ElectronSetupInteractions(G4HepEmElectronTrack *hepEMTracks, con
const bool returnAllSteps = options.returnAllSteps;
const bool returnLastStep = options.returnLastStep;

int activeSize = propagationQueue->size();
int activeSize = setupQueue->size();
for (int i = blockIdx.x * blockDim.x + threadIdx.x; i < activeSize; i += blockDim.x * gridDim.x) {
const int slot = (*propagationQueue)[i];
const int slot = (*setupQueue)[i];
ChargedTrack &currentTrack = electronsOrPositrons.TrackAt(slot);

// Retrieve HepEM track
G4HepEmElectronTrack &elTrack = hepEMTracks[slot];
G4HepEmTrack *theTrack = elTrack.GetTrack();

// the MCC vector is indexed by the logical volume id
const int lvolID = currentTrack.navState.GetLogicalId();

VolAuxData const &auxData = adept::transport::gVolAuxData[lvolID];

bool trackSurvives = true;

// Retrieve HepEM track
G4HepEmElectronTrack &elTrack = hepEMTracks[slot];
G4HepEmTrack *theTrack = elTrack.GetTrack();
if (!currentTrack.nextState.IsOnBoundary()) {
G4HepEmMSCTrackData *mscData = elTrack.GetMSCTrackData();
const double *mscDisplacement = mscData->GetDisplacement();
vecgeom::Vector3D<double> displacement(mscDisplacement[0], mscDisplacement[1], mscDisplacement[2]);
const double dLength2 = displacement.Length2();
constexpr double kGeomMinLength = 0.05 * copcore::units::nm; // 0.05 [nm]
constexpr double kGeomMinLength2 = kGeomMinLength * kGeomMinLength; // (0.05 [nm])^2
if (dLength2 > kGeomMinLength2) {
const double dispR = std::sqrt(dLength2);
// Estimate safety by subtracting the geometrical step length.
double safety = currentTrack.GetSafety(currentTrack.pos);
constexpr double sFact = 0.99;
double reducedSafety = sFact * safety;

G4HepEmMSCTrackData *mscData = elTrack.GetMSCTrackData();
// Apply displacement, depending on how close we are to a boundary.
// 1a. Far away from geometry boundary:
if (reducedSafety > 0.0 && dispR <= reducedSafety) {
currentTrack.pos += displacement;
} else {
// Recompute safety.
// Use maximum accuracy only if safety is smaller than the displacement.
safety = AdePTNavigator::ComputeSafety(currentTrack.pos, currentTrack.navState, dispR);
currentTrack.SetSafety(currentTrack.pos, safety);
reducedSafety = sFact * safety;

// 1b. Far away from geometry boundary:
if (reducedSafety > 0.0 && dispR <= reducedSafety) {
currentTrack.pos += displacement;
// 2. Push to boundary:
} else if (reducedSafety > kGeomMinLength) {
currentTrack.pos += displacement * (reducedSafety / dispR);
}
// 3. Very small safety: do nothing.
}
}
}

double energyDeposit = theTrack->GetEnergyDeposit();

bool reached_interaction = true;

// Set Non-stopped, on-boundary tracks for relocation
if (currentTrack.nextState.IsOnBoundary() && !currentTrack.stopped) {
// Add particle to relocation queue
splitQueues.queues[ParticleQueues::relocation]->push_back(slot);
continue;
}

auto winnerProcessIndex = theTrack->GetWinnerProcessIndex();

// Now check whether the non-relocating tracks reached an interaction
Expand Down
3 changes: 3 additions & 0 deletions include/AdePT/transport/queues/ParticleQueues.cuh
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -43,6 +43,9 @@ struct ParticleQueues {
adept::MParray *initiallyActive;
#ifdef ADEPT_USE_SPLIT_KERNELS
adept::MParray *propagation;
// Filled by charged-particle MSC for tracks that do not require relocation.
// This queue is allocated only for electrons and positrons.
adept::MParray *setupInteractions;
adept::MParray *splitQueues[numSplitQueues];
#endif

Expand Down
Loading