From c02b059a805ca301b986edbe55463ffc00f540e7 Mon Sep 17 00:00:00 2001 From: SeverinDiederichs Date: Wed, 12 Aug 2026 12:30:17 +0200 Subject: [PATCH] overlap safety call in MSC with relocation --- include/AdePT/transport/AdePTTransport.cuh | 55 ++++++--- .../transport/kernels/electrons_split.cuh | 112 +++++++++--------- .../AdePT/transport/queues/ParticleQueues.cuh | 3 + 3 files changed, 97 insertions(+), 73 deletions(-) diff --git a/include/AdePT/transport/AdePTTransport.cuh b/include/AdePT/transport/AdePTTransport.cuh index 04e4d6d4..80765b95 100644 --- a/include/AdePT/transport/AdePTTransport.cuh +++ b/include/AdePT/transport/AdePTTransport.cuh @@ -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]); } @@ -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(); } @@ -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(); } @@ -664,7 +673,12 @@ std::unique_ptr InitializeGPU(int trackCapacity, int particleType.queues.nextActive = static_cast(gpuPtr); #ifdef ADEPT_USE_SPLIT_KERNELS gpuMalloc(gpuPtr, sizeOfQueueStorage); - particleType.queues.propagation = static_cast(gpuPtr); + particleType.queues.propagation = static_cast(gpuPtr); + particleType.queues.setupInteractions = nullptr; + if (particleIndex != GPUQueueIndex::Gamma) { + gpuMalloc(gpuPtr, sizeOfQueueStorage); + particleType.queues.setupInteractions = static_cast(gpuPtr); + } for (int j = 0; j < ParticleQueues::numSplitQueues; j++) { gpuMalloc(gpuPtr, sizeOfQueueStorage); particleType.queues.splitQueues[j] = static_cast(gpuPtr); @@ -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, @@ -1033,15 +1050,16 @@ void TransportLoop(int trackCapacity, int stepCapacity, int numThreads, TrackBuf ElectronPropagation<<>>( electrons.tracks, gpuState.hepEmBuffers_d.electronsHepEm, electrons.queues.propagation); ElectronMSC<<>>( - electrons.tracks, gpuState.hepEmBuffers_d.electronsHepEm, electrons.queues.propagation); - ElectronSetupInteractions<<>>( - 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<<>>( gpuState.hepEmBuffers_d.electronsHepEm, particleManager, electrons.queues.splitQueues[ParticleQueues::relocation], kernelOptions); + ElectronSetupInteractions<<>>( + gpuState.hepEmBuffers_d.electronsHepEm, electrons.queues.setupInteractions, particleManager, + electronSplitQueues, kernelOptions); ElectronIonization<<>>( gpuState.hepEmBuffers_d.electronsHepEm, particleManager, electrons.queues.splitQueues[ParticleQueues::chargedIonization], kernelOptions); @@ -1075,15 +1093,16 @@ void TransportLoop(int trackCapacity, int stepCapacity, int numThreads, TrackBuf ElectronPropagation<<>>( positrons.tracks, gpuState.hepEmBuffers_d.positronsHepEm, positrons.queues.propagation); ElectronMSC<<>>( - positrons.tracks, gpuState.hepEmBuffers_d.positronsHepEm, positrons.queues.propagation); - ElectronSetupInteractions<<>>( - 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<<>>( gpuState.hepEmBuffers_d.positronsHepEm, particleManager, positrons.queues.splitQueues[ParticleQueues::relocation], kernelOptions); + ElectronSetupInteractions<<>>( + gpuState.hepEmBuffers_d.positronsHepEm, positrons.queues.setupInteractions, particleManager, + positronSplitQueues, kernelOptions); ElectronIonization<<>>( gpuState.hepEmBuffers_d.positronsHepEm, particleManager, positrons.queues.splitQueues[ParticleQueues::chargedIonization], kernelOptions); diff --git a/include/AdePT/transport/kernels/electrons_split.cuh b/include/AdePT/transport/kernels/electrons_split.cuh index d2088736..aadf3ce7 100644 --- a/include/AdePT/transport/kernels/electrons_split.cuh +++ b/include/AdePT/transport/kernels/electrons_split.cuh @@ -360,7 +360,8 @@ __global__ void ElectronPropagation(ChargedTrack *electronsOrPositrons, G4HepEmE } template -__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; @@ -369,57 +370,19 @@ __global__ void ElectronMSC(ChargedTrack *electrons, G4HepEmElectronTrack *hepEM const int slot = (*active)[i]; ChargedTrack ¤tTrack = 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(¤tTrack.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 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(); @@ -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 -__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) { @@ -449,10 +422,15 @@ __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 ¤tTrack = 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(); @@ -460,23 +438,47 @@ __global__ void ElectronSetupInteractions(G4HepEmElectronTrack *hepEMTracks, con 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 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 diff --git a/include/AdePT/transport/queues/ParticleQueues.cuh b/include/AdePT/transport/queues/ParticleQueues.cuh index d7ea00fa..4cdc82f1 100644 --- a/include/AdePT/transport/queues/ParticleQueues.cuh +++ b/include/AdePT/transport/queues/ParticleQueues.cuh @@ -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