From d7c27304aa45aa10a0c2675a4562f5cfd31715b6 Mon Sep 17 00:00:00 2001 From: John Apostolakis Date: Sat, 26 Feb 2022 08:06:19 +0100 Subject: [PATCH] Runge-Kutta integration: changes to identify cause of errors. (For debugging) --- examples/Example15/electrons.cu | 14 +-- examples/Example15/example15.cu | 2 +- magneticfield/inc/fieldPropagatorRungeKutta.h | 100 +++++++++++++++--- 3 files changed, 94 insertions(+), 22 deletions(-) diff --git a/examples/Example15/electrons.cu b/examples/Example15/electrons.cu index 313411e83..37d60a1c5 100644 --- a/examples/Example15/electrons.cu +++ b/examples/Example15/electrons.cu @@ -228,7 +228,7 @@ static __device__ __forceinline__ void TransportElectrons(Track *electrons, cons magneticFieldB, currentTrack.energy, Mass, Charge, stepLengthFromPhysics, currentTrack.pos, currentTrack.dir, currentTrack.navState, nextState, - propagated, /*lengthDone,*/ safety, max_iterations, iterDone + propagated, /*lengthDone,*/ safety, max_iterations, iterDone, i ); #ifdef CHECK_RESULTS constexpr Precision thesholdDiff=3.0e-5; @@ -241,18 +241,18 @@ static __device__ __forceinline__ void TransportElectrons(Track *electrons, cons CompareResponseVector3D( i, startDirection, directionHx, currentTrack.dir, "Direction", thesholdDiff ); const char* Outcome[2]={ "Good", " Bad" }; - // if( badPosition || badDirection) { - printf("%4s track (id= %3d) e_kin= %8.4g stepLen= %7.3g iters= %5d\n ", Outcome[badPosition||badDirection], - i, currentTrack.energy, stepLengthFromPhysics, iterDone); - currentTrack.print(i, /* verbose= */ true ); - // } + if( badPosition || badDirection) { + printf("%4s track (id= %3d) e_kin= %8.4g stepLen= %7.3g iters= %5d\n ", Outcome[badPosition||badDirection], + i, currentTrack.energy, stepLengthFromPhysics, iterDone); + currentTrack.print(i, /* verbose= */ true ); + } #endif #else fieldPropagatorConstBz fieldPropagatorBz(BzFieldValue); geometryStepLength = fieldPropagatorBz.ComputeStepAndNextVolume( currentTrack.energy, Mass, Charge, stepLengthFromPhysics, currentTrack.pos, currentTrack.dir, - currentTrack.navState, nextState, propagated); + currentTrack.navState, nextState, propagated, i); #endif } else { geometryStepLength = diff --git a/examples/Example15/example15.cu b/examples/Example15/example15.cu index 7d8ebf8e1..4633d99ef 100644 --- a/examples/Example15/example15.cu +++ b/examples/Example15/example15.cu @@ -492,7 +492,7 @@ void example15(int numParticles, double energy, int batch, const int *MCIndex_ho // std::cout << " Tracks AFTER Swap \n"; // std::cout << " Remaining: e- " << numElectrons << " g " << numGammas << " e+ " << numPositrons << "\n"; - printActiveTracks( electrons, positrons, gammas, verbTrk ); + // printActiveTracks( electrons, positrons, gammas, verbTrk ); if (numElectrons == previousElectrons && numPositrons == previousPositrons && numGammas == 0) { loopingNo++; diff --git a/magneticfield/inc/fieldPropagatorRungeKutta.h b/magneticfield/inc/fieldPropagatorRungeKutta.h index 47d501128..dbffd36c8 100644 --- a/magneticfield/inc/fieldPropagatorRungeKutta.h +++ b/magneticfield/inc/fieldPropagatorRungeKutta.h @@ -24,7 +24,9 @@ class fieldPropagatorRungeKutta { int charge, Real_t step, vecgeom::Vector3D &position, - vecgeom::Vector3D &direction); + vecgeom::Vector3D &direction + , int id // For debugging + ); static inline __host__ __device__ @@ -38,7 +40,8 @@ class fieldPropagatorRungeKutta { bool & propagated, const Real_t & /*safety*/, const int max_iterations - , int & iterDone + , int & iterDone + , int threadId ); // Move the track, // updating 'position', 'direction', the next state and returning the length moved. @@ -46,7 +49,10 @@ class fieldPropagatorRungeKutta { static inline __host__ __device__ void IntegrateTrackToEnd( Field_t const & magField, // RkDriver_t &integrationDriverRK, vecgeom::Vector3D & position, vecgeom::Vector3D & momentumVec, - int charge, Real_t stepLength ); + int charge, Real_t stepLength + , int id // Temporary - for debugging + , bool verbose = true // >>2 + ); static inline __host__ __device__ bool IntegrateTrackForProgress( Field_t const & magField, // RkDriver_t &integrationDriverRK, @@ -79,7 +85,9 @@ void fieldPropagatorRungeKutta int charge, Real_t step, vecgeom::Vector3D &position, - vecgeom::Vector3D &direction ) + vecgeom::Vector3D &direction + , int id // For debugging + ) { Real_t momentumMag = sqrt(kinE * (kinE + 2.0 * mass)); Real_t invMomentumMag = 1.0 / momentumMag; @@ -87,11 +95,13 @@ void fieldPropagatorRungeKutta // Only charged particles ( e-, e+, any .. ) can be propagated vecgeom::Vector3D positionVec = position; vecgeom::Vector3D momentumVec = momentumMag * direction; - IntegrateTrackToEnd /**/ ( magField, - positionVec, - momentumVec, - charge, - step ); + IntegrateTrackToEnd( magField, + positionVec, + momentumVec, + charge, + step + , id, true // For debugging + ); position = positionVec; direction = invMomentumMag * momentumVec; // Deviation of magnitude of direction from unit indicates integration error @@ -108,6 +118,8 @@ void fieldPropagatorRungeKutta::Integr vecgeom::Vector3D & momentumVec, int charge, Real_t stepLength + , int id + , bool verbose ) { // Version 1. Finish the integration of lanes ... @@ -115,28 +127,89 @@ void fieldPropagatorRungeKutta::Integr const unsigned int trialsPerCall = vecCore::Min( 30U, fMaxTrials / 2) ; // Parameter that can be tuned unsigned int totalTrials=0; + static int callNum = 0; + callNum ++; Real_t lenRemains = stepLength; Real_t hTry = stepLength; // suggested 'good' length per integration step bool unfinished = true; + + Real_t totLen = 0.0; + // unsigned int loopCt=0; do { Real_t hAdvanced = 0; // length integrated this iteration (of do-while) Real_t dydx_end[Nvar]; + + const vecgeom::Vector3D posBegin= position; // For print + const vecgeom::Vector3D momBegin= momentumVec; // >> bool done= RkDriver_t::Advance( position, momentumVec, charge, lenRemains, magField, hTry, dydx_end, - hAdvanced, totalTrials, trialsPerCall); + hAdvanced, totalTrials, + id, // Temporary ? + trialsPerCall); // Runge-Kutta single call ( number of steps <= trialsPerCall ) lenRemains -= hAdvanced; - unfinished = !done; // (lenRemains > 0.0); + unfinished = lenRemains > 0.0; /// Was = !done ... for debugging ???? + totLen+= hAdvanced; + + // if( !done || (loopCt++>0) ){ + if( id == 1 ) { + const vecgeom::Vector3D deltaPos= position - posBegin; + const vecgeom::Vector3D deltaMomentum= momentumVec - momBegin; + + printf(" id %3d call %4d iters %3d hdid= %9.5g " // totLen= %9.5g lenRemains= %9.5g " + " ret= %1d #= pos = %9.6g %9.6g %9.6g momemtumV= %14.9g %14.9g %14.9g hTry= %7.4g remains= %7.4g " + " Delta-pos= %7.3g %7.3g %7.3g (mag= %7.3g) Delta-mom= %7.3g %7.3g %7.3g (mag= %7.3g) " + " \n", + id, callNum, totalTrials, hAdvanced, // totLen, lenRemains, + done, + position[0], position[1], position[2], + momentumVec[0], momentumVec[1], momentumVec[2], + hTry, lenRemains + , deltaPos[0], deltaPos[1], deltaPos[2], deltaPos.Mag() + , deltaMomentum[0], deltaMomentum[1], deltaMomentum[2], deltaMomentum.Mag() + ); + } // sumAdvanced += hAdvanced; // Gravy .. } while ( unfinished && (totalTrials < fMaxTrials) ); } +#if 0 +template +inline __host__ __device__ +bool fieldPropagatorRungeKutta::IntegrateOneStep( + Field_t const & magField, + vecgeom::Vector3D & position, + vecgeom::Vector3D & momentumVec, + int charge, + Real_t & lenRemains, + Real_t & hTry, // Trial step size + Real_t & hAdvanced, // Length advanced + int & totalTrials + ) +// Return whether full length was completed +{ + Real_t hTry = stepLength; // suggested 'good' length per integration step + bool unfinished = true; + + hAdvanced = 0; // length integrated this iteration + Real_t dydx_end[Nvar]; + + bool done= + RkDriver_t::Advance( position, momentumVec, charge, lenRemains, magField, hTry, dydx_end, + hAdvanced, totalTrials, trialsPerCall); + // Runge-Kutta single call ( number of steps <= trialsPerCall ) + + lenRemains -= hAdvanced; + return done; +} +#endif + // ---------------------------------------------------------------------------- template @@ -259,6 +332,7 @@ fieldPropagatorRungeKutta ::ComputeSte 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 threadId ) { using copcore::units::MeV; @@ -285,8 +359,6 @@ fieldPropagatorRungeKutta ::ComputeSte Real_t ratioOverFld = (bmag>0) ? momentumVec.Dot(B0fieldVec) / (bmag * bmag) : 0.0; vecgeom::Vector3D PtransB = momentumVec - ratioOverFld * B0fieldVec; - // Real_t curv = fabs(Track::kB2C * charge * bmag / ( PtransB.Mag() + tiny)); - // Calculate inverse curvature instead - save a division Real_t inv_curv = fabs( PtransB.Mag() / ( Track::kB2C * charge * bmag + tiny) @@ -324,7 +396,7 @@ fieldPropagatorRungeKutta ::ComputeSte vecgeom::Vector3D endMomentumVec = momentumVec; // momentumMag * direction; Real_t safeArc = min(remains, safeLength); - IntegrateTrackToEnd( magField, endPosition, endMomentumVec, charge, safeArc); + IntegrateTrackToEnd( magField, endPosition, endMomentumVec, charge, safeArc, threadId); //----------------- vecgeom::Vector3D chordVec = endPosition - position; Real_t chordLen = chordVec.Length();