Skip to content

Commit

Permalink
Surface navigation now split in finding the frame and relocation
Browse files Browse the repository at this point in the history
  • Loading branch information
@agheata
agheata committed Sep 7, 2024
1 parent 2526590 commit 97c1ef4
Show file tree
Hide file tree
Showing 8 changed files with 182 additions and 128 deletions.
5 changes: 3 additions & 2 deletions include/AdePT/core/AdePTTransport.icc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -67,7 +67,7 @@ template <typename IntegrationLayer>
bool AdePTTransport<IntegrationLayer>::InitializeGeometry(const vecgeom::cxx::VPlacedVolume *world)
{
auto &cudaManager = vecgeom::cxx::CudaManager::Instance();
bool success = true;
bool success = true;
#ifdef ADEPT_USE_SURF
#ifdef ADEPT_USE_SURF_SINGLE
using BrepHelper = vgbrep::BrepHelper<float>;
Expand All @@ -83,9 +83,10 @@ bool AdePTTransport<IntegrationLayer>::InitializeGeometry(const vecgeom::cxx::VP
cudaManager.SynchronizeNavigationTable();
#else
// Upload solid geometry to GPU.
cudaDeviceSetLimit(cudaLimitStackSize, 1024 * 8);
cudaManager.LoadGeometry(world);
auto world_dev = cudaManager.Synchronize();
success = world_dev != nullptr;
success = world_dev != nullptr;
#endif
// Initialize BVH
InitBVH();
Expand Down
17 changes: 13 additions & 4 deletions include/AdePT/kernels/electrons.cuh
View file
Original file line number Diff line number Diff line change
Expand Up @@ -171,14 +171,20 @@ static __device__ __forceinline__ void TransportElectrons(adept::TrackManager<Tr

// Check if there's a volume boundary in between.
bool propagated = true;
long hitsurf_index = -1;
double geometryStepLength;
vecgeom::NavigationState nextState;
if (BzFieldValue != 0) {
geometryStepLength = fieldPropagatorBz.ComputeStepAndNextVolume<AdePTNavigator>(
eKin, restMass, Charge, geometricalStepLengthFromPhysics, pos, dir, navState, nextState, propagated, safety);
eKin, restMass, Charge, geometricalStepLengthFromPhysics, pos, dir, navState, nextState, hitsurf_index, propagated, safety);
} else {
geometryStepLength = AdePTNavigator::ComputeStepAndNextVolume(pos, dir, geometricalStepLengthFromPhysics,
navState, nextState, kPush);
#ifdef ADEPT_USE_SURF
double geometryStepLength = AdePTNavigator::ComputeStepAndNextVolume(pos, dir, geometricalStepLengthFromPhysics,
navState, nextState, hitsurf_index, kPush);
#else
double geometryStepLength = AdePTNavigator::ComputeStepAndNextVolume(pos, dir, geometricalStepLengthFromPhysics,
navState, nextState, kPush);
#endif
pos += geometryStepLength * dir;
}

Expand Down Expand Up @@ -306,8 +312,11 @@ static __device__ __forceinline__ void TransportElectrons(adept::TrackManager<Tr

// Kill the particle if it left the world.
if (!nextState.IsOutside()) {
#ifdef ADEPT_USE_SURF
AdePTNavigator::RelocateToNextVolume(pos, dir, hitsurf_index, nextState);
#else
AdePTNavigator::RelocateToNextVolume(pos, dir, nextState);

#endif
// Move to the next boundary.
navState = nextState;
// Check if the next volume belongs to the GPU region and push it to the appropriate queue
Expand Down
18 changes: 15 additions & 3 deletions include/AdePT/kernels/gammas.cuh
View file
Original file line number Diff line number Diff line change
Expand Up @@ -95,8 +95,16 @@ __global__ void TransportGammas(adept::TrackManager<Track> *gammas, Secondaries

// Check if there's a volume boundary in between.
vecgeom::NavigationState nextState;
#ifdef ADEPT_USE_SURF
long hitsurf_index = -1;
double geometryStepLength = AdePTNavigator::ComputeStepAndNextVolume(pos, dir, geometricalStepLengthFromPhysics,
navState, nextState, hitsurf_index, kPush);
#else
double geometryStepLength = AdePTNavigator::ComputeStepAndNextVolume(pos, dir, geometricalStepLengthFromPhysics,
navState, nextState, kPush);
#endif
// printf("pvol=%d step=%g onboundary=%d pos={%g, %g, %g} dir={%g, %g, %g}\n", navState.TopId(), geometryStepLength,
// nextState.IsOnBoundary(), pos[0], pos[1], pos[2], dir[0], dir[1], dir[2]);
pos += geometryStepLength * dir;

// Set boundary state in navState so the next step and secondaries get the
Expand All @@ -121,15 +129,19 @@ __global__ void TransportGammas(adept::TrackManager<Track> *gammas, Secondaries

// Kill the particle if it left the world.
if (!nextState.IsOutside()) {
#ifdef ADEPT_USE_SURF
AdePTNavigator::RelocateToNextVolume(pos, dir, hitsurf_index, nextState);
#else
AdePTNavigator::RelocateToNextVolume(pos, dir, nextState);

#endif
// Move to the next boundary.
navState = nextState;
//printf(" -> pvol=%d pos={%g, %g, %g} \n", navState.TopId(), pos[0], pos[1], pos[2]);
// Check if the next volume belongs to the GPU region and push it to the appropriate queue
#ifndef ADEPT_USE_SURF
const int nextlvolID = navState.Top()->GetLogicalVolume()->id();
const int nextlvolID = navState.Top()->GetLogicalVolume()->id();
#else
const int nextlvolID = navState.GetLogicalId();
const int nextlvolID = navState.GetLogicalId();
#endif
VolAuxData const &nextauxData = auxDataArray[nextlvolID];
if (nextauxData.fGPUregion > 0)
Expand Down
200 changes: 101 additions & 99 deletions include/AdePT/magneticfield/fieldPropagatorConstBz.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -32,8 +32,9 @@ class fieldPropagatorConstBz {
__host__ __device__ Precision ComputeStepAndNextVolume(double kinE, double mass, int charge, Precision physicsStep,
Vector3D &position, Vector3D &direction,
vecgeom::NavigationState const &current_state,
vecgeom::NavigationState &new_state, bool &propagated,
const Precision safety = 0.0, const int max_iteration = 100);
vecgeom::NavigationState &new_state, long &hitsurf_index,
bool &propagated, const Precision safety = 0.0,
const int max_iteration = 100);

private:
Precision BzValue;
Expand Down Expand Up @@ -87,7 +88,8 @@ template <class Navigator>
__host__ __device__ Precision fieldPropagatorConstBz::ComputeStepAndNextVolume(
double kinE, double mass, int charge, Precision physicsStep, vecgeom::Vector3D<vecgeom::Precision> &position,
vecgeom::Vector3D<vecgeom::Precision> &direction, vecgeom::NavigationState const &current_state,
vecgeom::NavigationState &next_state, bool &propagated, const vecgeom::Precision safetyIn, const int max_iterations)
vecgeom::NavigationState &next_state, long &hitsurf_index, bool &propagated, const vecgeom::Precision safetyIn,
const int max_iterations)
{
using Precision = vecgeom::Precision;
#ifdef VECGEOM_FLOAT_PRECISION
Expand All @@ -108,107 +110,107 @@ __host__ __device__ Precision fieldPropagatorConstBz::ComputeStepAndNextVolume(
const Precision epsilon_step = 1.0e-7 * physicsStep; // Ignore remainder if < e_s * PhysicsStep
int chordIters = 0;

if (charge == 0) {
stepDone = Navigator::ComputeStepAndNextVolume(position, direction, remains, current_state, next_state, kPush);
position += stepDone * direction;
} else {
bool continueIteration = false;

Precision safety = safetyIn;
Vector3D safetyOrigin = position;
// Prepare next_state in case we skip navigation inside the safety sphere.
current_state.CopyTo(&next_state);
next_state.SetBoundaryState(false);

Precision maxNextSafeMove = safeLength;

bool lastWasZero = false;
// Locate the intersection of the curved trajectory and the boundaries of the current
// volume (including daughters).
do {
Vector3D endPosition = position;
Vector3D endDirection = direction;
Precision safeMove = min(remains, maxNextSafeMove);

helixBz.DoStep<Vector3D, Precision, int>(position, direction, charge, momentumMag, safeMove, endPosition,
endDirection);

Vector3D chordVec = endPosition - position;
Precision chordLen = chordVec.Length();
Vector3D chordDir = (1 / chordLen) * chordVec;

Precision currentSafety = safety - (position - safetyOrigin).Length();
Precision move;
if (currentSafety > chordLen) {
move = chordLen;
} else {
Precision newSafety = 0;
if (stepDone > 0) {
newSafety = Navigator::ComputeSafety(position, current_state);
}
if (newSafety > chordLen) {
move = chordLen;
safetyOrigin = position;
safety = newSafety;
} else {
move = Navigator::ComputeStepAndNextVolume(position, chordDir, chordLen, current_state, next_state, kPush);
}
}
bool continueIteration = false;

static constexpr Precision ReduceFactor = 0.1;
static constexpr int ReduceIters = 6;
Precision safety = safetyIn;
Vector3D safetyOrigin = position;
// Prepare next_state in case we skip navigation inside the safety sphere.
current_state.CopyTo(&next_state);
next_state.SetBoundaryState(false);

if (lastWasZero && chordIters >= ReduceIters) {
lastWasZero = false;
}
Precision maxNextSafeMove = safeLength;

if (move == chordLen) {
position = endPosition;
direction = endDirection;
move = safeMove;
// We want to try the maximum step in the next iteration.
maxNextSafeMove = safeLength;
continueIteration = true;
} else if (move <= kPush + Navigator::kBoundaryPush && stepDone == 0) {
// FIXME: Even for zero steps, the Navigator will return kPush + possibly
// Navigator::kBoundaryPush instead of a real 0.
move = 0;
lastWasZero = true;

// Reduce the step attempted in the next iteration to navigate around
// boundaries where the chord step may end in a volume we just left.
maxNextSafeMove = ReduceFactor * safeMove;
continueIteration = chordIters < ReduceIters;

if (!continueIteration) {
// Let's move to the other side of this boundary -- this side we cannot progress !!
move = Navigator::kBoundaryPush;
// printf("fieldProp-ConstBz: pushing by %10.4g \n ", move );
}
bool lastWasZero = false;
// Locate the intersection of the curved trajectory and the boundaries of the current
// volume (including daughters).
do {
Vector3D endPosition = position;
Vector3D endDirection = direction;
Precision safeMove = min(remains, maxNextSafeMove);

helixBz.DoStep<Vector3D, Precision, int>(position, direction, charge, momentumMag, safeMove, endPosition,
endDirection);

Vector3D chordVec = endPosition - position;
Precision chordLen = chordVec.Length();
Vector3D chordDir = (1 / chordLen) * chordVec;

Precision currentSafety = safety - (position - safetyOrigin).Length();
Precision move;
if (currentSafety > chordLen) {
move = chordLen;
} else {
Precision newSafety = 0;
if (stepDone > 0) {
newSafety = Navigator::ComputeSafety(position, current_state);
}
if (newSafety > chordLen) {
move = chordLen;
safetyOrigin = position;
safety = newSafety;
} else {
// Accept the intersection point on the surface. This means that
// the point at the boundary will be on the 'straight'-line chord,
// not the curved trajectory.
// ( This involves a bias -- relevant for muons in trackers.
// Currently it's controlled/limited by the acceptable step size ie. 'safeLength' )
position = position + move * chordDir;

// Primitive approximation of end direction and move to the crossing point ...
Precision fraction = chordLen > 0 ? move / chordLen : 0;
direction = direction * (1.0 - fraction) + endDirection * fraction;
direction = direction.Unit();
// safeMove is how much the track would have been moved if not hitting the boundary
// We approximate the actual reduction along the curved trajectory to be the same
// as the reduction of the full chord due to the boundary crossing.
move = fraction * safeMove;
continueIteration = false;
#ifdef ADEPT_USE_SURF
move = Navigator::ComputeStepAndNextVolume(position, chordDir, chordLen, current_state, next_state,
hitsurf_index, kPush);
#else
move = Navigator::ComputeStepAndNextVolume(position, chordDir, chordLen, current_state, next_state, kPush);
#endif
}
stepDone += move;
remains -= move;
chordIters++;

} while (continueIteration && (remains > epsilon_step) && (chordIters < max_iterations));
}
}

static constexpr Precision ReduceFactor = 0.1;
static constexpr int ReduceIters = 6;

if (lastWasZero && chordIters >= ReduceIters) {
lastWasZero = false;
}

if (move == chordLen) {
position = endPosition;
direction = endDirection;
move = safeMove;
// We want to try the maximum step in the next iteration.
maxNextSafeMove = safeLength;
continueIteration = true;
} else if (move <= kPush + Navigator::kBoundaryPush && stepDone == 0) {
// FIXME: Even for zero steps, the Navigator will return kPush + possibly
// Navigator::kBoundaryPush instead of a real 0.
move = 0;
lastWasZero = true;

// Reduce the step attempted in the next iteration to navigate around
// boundaries where the chord step may end in a volume we just left.
maxNextSafeMove = ReduceFactor * safeMove;
continueIteration = chordIters < ReduceIters;

if (!continueIteration) {
// Let's move to the other side of this boundary -- this side we cannot progress !!
move = Navigator::kBoundaryPush;
// printf("fieldProp-ConstBz: pushing by %10.4g \n ", move );
}
} else {
// Accept the intersection point on the surface. This means that
// the point at the boundary will be on the 'straight'-line chord,
// not the curved trajectory.
// ( This involves a bias -- relevant for muons in trackers.
// Currently it's controlled/limited by the acceptable step size ie. 'safeLength' )
position = position + move * chordDir;

// Primitive approximation of end direction and move to the crossing point ...
Precision fraction = chordLen > 0 ? move / chordLen : 0;
direction = direction * (1.0 - fraction) + endDirection * fraction;
direction = direction.Unit();
// safeMove is how much the track would have been moved if not hitting the boundary
// We approximate the actual reduction along the curved trajectory to be the same
// as the reduction of the full chord due to the boundary crossing.
move = fraction * safeMove;
continueIteration = false;
}
stepDone += move;
remains -= move;
chordIters++;

} while (continueIteration && (remains > epsilon_step) && (chordIters < max_iterations));

propagated = (chordIters < max_iterations);
return stepDone;
Expand Down
20 changes: 15 additions & 5 deletions include/AdePT/magneticfield/fieldPropagatorRungeKutta.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -24,8 +24,8 @@ class fieldPropagatorRungeKutta {
static inline __host__ __device__ __host__ __device__ Real_t ComputeStepAndNextVolume(
Field_t const &magneticField, double kinE, double mass, int charge, double physicsStep,
vecgeom::Vector3D<Real_t> &position, vecgeom::Vector3D<Real_t> &direction,
vecgeom::NavigationState const &current_state, vecgeom::NavigationState &next_state, bool &propagated,
const Real_t & /*safety*/, const int max_iterations, int &iterDone, int threadId);
vecgeom::NavigationState const &current_state, vecgeom::NavigationState &next_state, long &hitsurf_index,
bool &propagated, const Real_t & /*safety*/, const int max_iterations, int &iterDone, int threadId);
// Move the track,
// updating 'position', 'direction', the next state and returning the length moved.
protected:
Expand Down Expand Up @@ -192,9 +192,9 @@ inline __host__ __device__ Real_t
fieldPropagatorRungeKutta<Field_t, RkDriver_t, Real_t, Navigator_t>::ComputeStepAndNextVolume(
Field_t const &magField, double kinE, double mass, int charge, double physicsStep,
vecgeom::Vector3D<Real_t> &position, vecgeom::Vector3D<Real_t> &direction,
vecgeom::NavigationState const &current_state, vecgeom::NavigationState &next_state, bool &propagated,
const Real_t & /*safety*/, // eventually In/Out ?
const int max_iterations, int &itersDone // useful for now - to monitor and report -- unclear if needed later
vecgeom::NavigationState const &current_state, vecgeom::NavigationState &next_state, long &hitsurf_index,
bool &propagated, const Real_t & /*safety*/, // eventually In/Out ?
const int max_iterations, int &itersDone // useful for now - to monitor and report -- unclear if needed later
,
int indx)
{
Expand Down Expand Up @@ -225,7 +225,12 @@ fieldPropagatorRungeKutta<Field_t, RkDriver_t, Real_t, Navigator_t>::ComputeStep
bool found_end = false;

if (inZeroFieldRegion) {
#ifdef ADEPT_USE_SURF
stepDone = Navigator_t::ComputeStepAndNextVolume(position, direction, remains, current_state, next_state,
hitsurf_index, kPush);
#else
stepDone = Navigator_t::ComputeStepAndNextVolume(position, direction, remains, current_state, next_state, kPush);
#endif
position += stepDone * direction;
} else {
bool continueIteration = false;
Expand Down Expand Up @@ -254,8 +259,13 @@ fieldPropagatorRungeKutta<Field_t, RkDriver_t, Real_t, Navigator_t>::ComputeStep

// Check Intersection
//-- vecgeom::Vector3D<Real_t> ChordDir= (1.0/chordDist) * ChordVec;
#ifdef ADEPT_USE_SURF
Real_t linearStep = Navigator_t::ComputeStepAndNextVolume(position, chordVec, chordDist, current_state,
next_state, hitsurf_index, kPush);
#else
Real_t linearStep =
Navigator_t::ComputeStepAndNextVolume(position, chordVec, chordDist, current_state, next_state, kPush);
#endif
Real_t curvedStep;

if (lastWasZero && chordIters >= ReduceIters) {
Expand Down
Loading

0 comments on commit 97c1ef4

Please sign in to comment.