G4MagInt_Driver Class Reference

#include <G4MagIntegratorDriver.hh>


Public Member Functions

G4bool AccurateAdvance (G4FieldTrack &y_current, G4double hstep, G4double eps, G4double hinitial=0.0)
G4bool QuickAdvance (G4FieldTrack &y_val, const G4double dydx[], G4double hstep, G4double &dchord_step, G4double &dyerr)
G4bool QuickAdvance (G4FieldTrack &y_posvel, const G4double dydx[], G4double hstep, G4double &dchord_step, G4double &dyerr_pos_sq, G4double &dyerr_mom_rel_sq)
 G4MagInt_Driver (G4double hminimum, G4MagIntegratorStepper *pItsStepper, G4int numberOfComponents=6, G4int statisticsVerbosity=1)
 ~G4MagInt_Driver ()
G4double GetHmin () const
G4double Hmin () const
G4double GetSafety () const
G4double GetPshrnk () const
G4double GetPgrow () const
G4double GetErrcon () const
void GetDerivatives (const G4FieldTrack &y_curr, G4double dydx[])
void RenewStepperAndAdjust (G4MagIntegratorStepper *pItsStepper)
void ReSetParameters (G4double new_safety=0.9)
void SetSafety (G4double valS)
void SetPshrnk (G4double valPs)
void SetPgrow (G4double valPg)
void SetErrcon (G4double valEc)
G4double ComputeAndSetErrcon ()
void SetChargeMomentumMass (G4double particleCharge, G4double MomentumXc, G4double Mass)
const G4MagIntegratorStepperGetStepper () const
void OneGoodStep (G4double ystart[], const G4double dydx[], G4double &x, G4double htry, G4double eps, G4double &hdid, G4double &hnext)
G4double ComputeNewStepSize (G4double errMaxNorm, G4double hstepCurrent)
G4double ComputeNewStepSize_WithinLimits (G4double errMaxNorm, G4double hstepCurrent)
G4int GetMaxNoSteps () const
void SetMaxNoSteps (G4int val)
void SetHmin (G4double newval)
void SetVerboseLevel (G4int newLevel)
G4double GetVerboseLevel () const
G4double GetSmallestFraction () const
void SetSmallestFraction (G4double val)

Protected Member Functions

void WarnSmallStepSize (G4double hnext, G4double hstep, G4double h, G4double xDone, G4int noSteps)
void WarnTooManyStep (G4double x1start, G4double x2end, G4double xCurrent)
void WarnEndPointTooFar (G4double endPointDist, G4double hStepSize, G4double epsilonRelative, G4int debugFlag)
void PrintStatus (const G4double *StartArr, G4double xstart, const G4double *CurrentArr, G4double xcurrent, G4double requestStep, G4int subStepNo)
void PrintStatus (const G4FieldTrack &StartFT, const G4FieldTrack &CurrentFT, G4double requestStep, G4int subStepNo)
void PrintStat_Aux (const G4FieldTrack &aFieldTrack, G4double requestStep, G4double actualStep, G4int subStepNo, G4double subStepSize, G4double dotVelocities)
void PrintStatisticsReport ()


Detailed Description

Definition at line 48 of file G4MagIntegratorDriver.hh.


Constructor & Destructor Documentation

G4MagInt_Driver::G4MagInt_Driver ( G4double  hminimum,
G4MagIntegratorStepper pItsStepper,
G4int  numberOfComponents = 6,
G4int  statisticsVerbosity = 1 
)

Definition at line 69 of file G4MagIntegratorDriver.cc.

References G4cout, G4endl, G4MagIntegratorStepper::IntegratorOrder(), and RenewStepperAndAdjust().

00073   : fSmallestFraction( 1.0e-12 ), 
00074     fNoIntegrationVariables(numComponents), 
00075     fMinNoVars(12), 
00076     fNoVars( std::max( fNoIntegrationVariables, fMinNoVars )),
00077     fStatisticsVerboseLevel(statisticsVerbose),
00078     fNoTotalSteps(0),  fNoBadSteps(0), fNoSmallSteps(0),
00079     fNoInitialSmallSteps(0), 
00080     fDyerr_max(0.0), fDyerr_mx2(0.0), 
00081     fDyerrPos_smTot(0.0), fDyerrPos_lgTot(0.0), fDyerrVel_lgTot(0.0), 
00082     fSumH_sm(0.0), fSumH_lg(0.0),
00083     fVerboseLevel(0)
00084 {  
00085   // In order to accomodate "Laboratory Time", which is [7], fMinNoVars=8
00086   // is required. For proper time of flight and spin,  fMinNoVars must be 12
00087 
00088   RenewStepperAndAdjust( pStepper );
00089   fMinimumStep= hminimum;
00090   fMaxNoSteps = fMaxStepBase / pIntStepper->IntegratorOrder();
00091 #ifdef G4DEBUG_FIELD
00092   fVerboseLevel=2;
00093 #endif
00094 
00095   if( (fVerboseLevel > 0) || (fStatisticsVerboseLevel > 1) )
00096   {
00097     G4cout << "MagIntDriver version: Accur-Adv: "
00098            << "invE_nS, QuickAdv-2sqrt with Statistics "
00099 #ifdef G4FLD_STATS
00100            << " enabled "
00101 #else
00102            << " disabled "
00103 #endif
00104            << G4endl;
00105   }
00106 }

G4MagInt_Driver::~G4MagInt_Driver (  ) 

Definition at line 112 of file G4MagIntegratorDriver.cc.

References PrintStatisticsReport().

00113 { 
00114   if( fStatisticsVerboseLevel > 1 )
00115   {
00116     PrintStatisticsReport();
00117   }
00118 }


Member Function Documentation

G4bool G4MagInt_Driver::AccurateAdvance ( G4FieldTrack y_current,
G4double  hstep,
G4double  eps,
G4double  hinitial = 0.0 
)

Definition at line 127 of file G4MagIntegratorDriver.cc.

References ComputeNewStepSize(), G4MagIntegratorStepper::ComputeRightHandSide(), G4FieldTrack::DumpToArray(), EventMustBeAborted, FatalException, G4cerr, G4cout, G4endl, G4Exception(), G4FieldTrack::GetCurveLength(), Hmin(), JustWarning, G4FieldTrack::LoadFromArray(), G4FieldTrack::ncompSVEC, OneGoodStep(), PrintStatus(), QuickAdvance(), G4FieldTrack::SetCurveLength(), WarnEndPointTooFar(), WarnSmallStepSize(), and WarnTooManyStep().

Referenced by G4ChordFinder::AdvanceChordLimited(), G4ChordFinder::ApproxCurvePointS(), G4ChordFinder::ApproxCurvePointV(), G4MultiLevelLocator::EstimateIntersectionPoint(), G4BrentLocator::EstimateIntersectionPoint(), and G4VIntersectionLocator::ReEstimateEndpoint().

00131 {
00132   // Runge-Kutta driver with adaptive stepsize control. Integrate starting
00133   // values at y_current over hstep x2 with accuracy eps. 
00134   // On output ystart is replaced by values at the end of the integration 
00135   // interval. RightHandSide is the right-hand side of ODE system. 
00136   // The source is similar to odeint routine from NRC p.721-722 .
00137 
00138   G4int nstp, i, no_warnings=0;
00139   G4double x, hnext, hdid, h;
00140 
00141 #ifdef G4DEBUG_FIELD
00142   static G4int dbg=1;
00143   static G4int nStpPr=50;   // For debug printing of long integrations
00144   G4double ySubStepStart[G4FieldTrack::ncompSVEC];
00145   G4FieldTrack  yFldTrkStart(y_current);
00146 #endif
00147 
00148   G4double y[G4FieldTrack::ncompSVEC], dydx[G4FieldTrack::ncompSVEC];
00149   G4double ystart[G4FieldTrack::ncompSVEC], yEnd[G4FieldTrack::ncompSVEC]; 
00150   G4double  x1, x2;
00151   G4bool succeeded = true, lastStepSucceeded;
00152 
00153   G4double startCurveLength;
00154 
00155   G4int  noFullIntegr=0, noSmallIntegr = 0 ;
00156   static G4int  noGoodSteps =0 ;  // Bad = chord > curve-len 
00157   const  G4int  nvar= fNoVars;
00158 
00159   G4FieldTrack yStartFT(y_current);
00160 
00161   //  Ensure that hstep > 0
00162   //
00163   if( hstep <= 0.0 )
00164   { 
00165     if(hstep==0.0)
00166     {
00167       std::ostringstream message;
00168       message << "Proposed step is zero; hstep = " << hstep << " !";
00169       G4Exception("G4MagInt_Driver::AccurateAdvance()", 
00170                   "GeomField1001", JustWarning, message);
00171       return succeeded; 
00172     }
00173     else
00174     { 
00175       std::ostringstream message;
00176       message << "Invalid run condition." << G4endl
00177               << "Proposed step is negative; hstep = " << hstep << "." << G4endl
00178               << "Requested step cannot be negative! Aborting event.";
00179       G4Exception("G4MagInt_Driver::AccurateAdvance()", 
00180                   "GeomField0003", EventMustBeAborted, message);
00181       return false;
00182     }
00183   }
00184 
00185   y_current.DumpToArray( ystart );
00186 
00187   startCurveLength= y_current.GetCurveLength();
00188   x1= startCurveLength; 
00189   x2= x1 + hstep;
00190 
00191   if ( (hinitial > 0.0) && (hinitial < hstep)
00192     && (hinitial > perMillion * hstep) )
00193   {
00194      h = hinitial;
00195   }
00196   else  //  Initial Step size "h" defaults to the full interval
00197   {
00198      h = hstep;
00199   }
00200 
00201   x = x1;
00202 
00203   for (i=0;i<nvar;i++)  { y[i] = ystart[i]; }
00204 
00205   G4bool lastStep= false;
00206   nstp=1;
00207 
00208   do
00209   {
00210     G4ThreeVector StartPos( y[0], y[1], y[2] );
00211 
00212 #ifdef G4DEBUG_FIELD
00213     G4double xSubStepStart= x; 
00214     for (i=0;i<nvar;i++)  { ySubStepStart[i] = y[i]; }
00215     yFldTrkStart.LoadFromArray(y, fNoIntegrationVariables);
00216     yFldTrkStart.SetCurveLength(x);
00217 #endif
00218 
00219     // Old method - inline call to Equation of Motion
00220     //   pIntStepper->RightHandSide( y, dydx );
00221     // New method allows to cache field, or state (eg momentum magnitude)
00222     pIntStepper->ComputeRightHandSide( y, dydx );
00223     fNoTotalSteps++;
00224 
00225     // Perform the Integration
00226     //      
00227     if( h > fMinimumStep )
00228     { 
00229       OneGoodStep(y,dydx,x,h,eps,hdid,hnext) ;
00230       //--------------------------------------
00231       lastStepSucceeded= (hdid == h);   
00232 #ifdef G4DEBUG_FIELD
00233       if (dbg>2) {
00234         PrintStatus( ySubStepStart, xSubStepStart, y, x, h,  nstp); // Only
00235       }
00236 #endif
00237     }
00238     else
00239     {
00240       G4FieldTrack yFldTrk( G4ThreeVector(0,0,0), 
00241                             G4ThreeVector(0,0,0), 0., 0., 0., 0. );
00242       G4double dchord_step, dyerr, dyerr_len;   // What to do with these ?
00243       yFldTrk.LoadFromArray(y, fNoIntegrationVariables); 
00244       yFldTrk.SetCurveLength( x );
00245 
00246       QuickAdvance( yFldTrk, dydx, h, dchord_step, dyerr_len ); 
00247       //-----------------------------------------------------
00248 
00249       yFldTrk.DumpToArray(y);    
00250 
00251 #ifdef G4FLD_STATS
00252       fNoSmallSteps++; 
00253       if ( dyerr_len > fDyerr_max)  { fDyerr_max= dyerr_len; }
00254       fDyerrPos_smTot += dyerr_len;
00255       fSumH_sm += h;  // Length total for 'small' steps
00256       if (nstp<=1)  { fNoInitialSmallSteps++; }
00257 #endif
00258 #ifdef G4DEBUG_FIELD
00259       if (dbg>1)
00260       {
00261         if(fNoSmallSteps<2) { PrintStatus(ySubStepStart, x1, y, x, h, -nstp); }
00262         G4cout << "Another sub-min step, no " << fNoSmallSteps 
00263                << " of " << fNoTotalSteps << " this time " << nstp << G4endl; 
00264         PrintStatus( ySubStepStart, x1, y, x, h,  nstp);   // Only this
00265         G4cout << " dyerr= " << dyerr_len << " relative = " << dyerr_len / h 
00266                << " epsilon= " << eps << " hstep= " << hstep 
00267                << " h= " << h << " hmin= " << fMinimumStep << G4endl;
00268       }
00269 #endif        
00270       if( h == 0.0 )
00271       { 
00272         G4Exception("G4MagInt_Driver::AccurateAdvance()",
00273                     "GeomField0003", FatalException,
00274                     "Integration Step became Zero!"); 
00275       }
00276       dyerr = dyerr_len / h;
00277       hdid= h;
00278       x += hdid;
00279 
00280       // Compute suggested new step
00281       hnext= ComputeNewStepSize( dyerr/eps, h);
00282 
00283       // .. hnext= ComputeNewStepSize_WithinLimits( dyerr/eps, h);
00284       lastStepSucceeded= (dyerr<= eps);
00285     }
00286 
00287     if (lastStepSucceeded)  { noFullIntegr++; }
00288     else                    { noSmallIntegr++; }
00289 
00290     G4ThreeVector EndPos( y[0], y[1], y[2] );
00291 
00292 #ifdef  G4DEBUG_FIELD
00293     if( (dbg>0) && (dbg<=2) && (nstp>nStpPr))
00294     {
00295       if( nstp==nStpPr )  { G4cout << "***** Many steps ****" << G4endl; }
00296       G4cout << "MagIntDrv: " ; 
00297       G4cout << "hdid="  << std::setw(12) << hdid  << " "
00298              << "hnext=" << std::setw(12) << hnext << " " 
00299              << "hstep=" << std::setw(12) << hstep << " (requested) " 
00300              << G4endl;
00301       PrintStatus( ystart, x1, y, x, h, (nstp==nStpPr) ? -nstp: nstp); 
00302     }
00303 #endif
00304 
00305     // Check the endpoint
00306     G4double endPointDist= (EndPos-StartPos).mag(); 
00307     if ( endPointDist >= hdid*(1.+perMillion) )
00308     {
00309       fNoBadSteps++;
00310 
00311       // Issue a warning only for gross differences -
00312       // we understand how small difference occur.
00313       if ( endPointDist >= hdid*(1.+perThousand) )
00314       { 
00315 #ifdef G4DEBUG_FIELD
00316         if (dbg)
00317         {
00318           WarnEndPointTooFar ( endPointDist, hdid, eps, dbg ); 
00319           G4cerr << "  Total steps:  bad " << fNoBadSteps
00320                  << " good " << noGoodSteps << " current h= " << hdid
00321                  << G4endl;
00322           PrintStatus( ystart, x1, y, x, hstep, no_warnings?nstp:-nstp);  
00323         }
00324 #endif
00325         no_warnings++;
00326       }
00327     }
00328     else
00329     {
00330       noGoodSteps ++;
00331     } 
00332 // #endif
00333 
00334     //  Avoid numerous small last steps
00335     if( (h < eps * hstep) || (h < fSmallestFraction * startCurveLength) )
00336     {
00337       // No more integration -- the next step will not happen
00338       lastStep = true;  
00339     }
00340     else
00341     {
00342       // Check the proposed next stepsize
00343       if(std::fabs(hnext) <= Hmin())
00344       {
00345 #ifdef  G4DEBUG_FIELD
00346         // If simply a very small interval is being integrated, do not warn
00347         if( (x < x2 * (1-eps) ) &&        // The last step can be small: OK
00348             (std::fabs(hstep) > Hmin()) ) // and if we are asked, it's OK
00349         {
00350           if(dbg>0)
00351           {
00352             WarnSmallStepSize( hnext, hstep, h, x-x1, nstp );  
00353             PrintStatus( ystart, x1, y, x, hstep, no_warnings?nstp:-nstp);
00354           }
00355           no_warnings++;
00356         }
00357 #endif
00358         // Make sure that the next step is at least Hmin.
00359         h = Hmin();
00360       }
00361       else
00362       {
00363         h = hnext;
00364       }
00365 
00366       //  Ensure that the next step does not overshoot
00367       if ( x+h > x2 )
00368       {                // When stepsize overshoots, decrease it!
00369         h = x2 - x ;   // Must cope with difficult rounding-error
00370       }                // issues if hstep << x2
00371 
00372       if ( h == 0.0 )
00373       {
00374         // Cannot progress - accept this as last step - by default
00375         lastStep = true;
00376 #ifdef G4DEBUG_FIELD
00377         if (dbg>2)
00378         {
00379           int prec= G4cout.precision(12); 
00380           G4cout << "Warning: G4MagIntegratorDriver::AccurateAdvance"
00381                  << G4endl
00382                  << "  Integration step 'h' became "
00383                  << h << " due to roundoff. " << G4endl
00384                  << " Calculated as difference of x2= "<< x2 << " and x=" << x
00385                  << "  Forcing termination of advance." << G4endl;
00386           G4cout.precision(prec);
00387         }          
00388 #endif
00389       }
00390     }
00391   } while ( ((nstp++)<=fMaxNoSteps) && (x < x2) && (!lastStep) );
00392      // Have we reached the end ?
00393      // --> a better test might be x-x2 > an_epsilon
00394 
00395   succeeded=  (x>=x2);  // If it was a "forced" last step
00396 
00397   for (i=0;i<nvar;i++)  { yEnd[i] = y[i]; }
00398 
00399   // Put back the values.
00400   y_current.LoadFromArray( yEnd, fNoIntegrationVariables );
00401   y_current.SetCurveLength( x );
00402 
00403   if(nstp > fMaxNoSteps)
00404   {
00405     no_warnings++;
00406     succeeded = false;
00407 #ifdef G4DEBUG_FIELD
00408     if (dbg)
00409     {
00410       WarnTooManyStep( x1, x2, x );  //  Issue WARNING
00411       PrintStatus( yEnd, x1, y, x, hstep, -nstp);
00412     }
00413 #endif
00414   }
00415 
00416 #ifdef G4DEBUG_FIELD
00417   if( dbg && no_warnings )
00418   {
00419     G4cerr << "G4MagIntegratorDriver exit status: no-steps " << nstp <<G4endl;
00420     PrintStatus( yEnd, x1, y, x, hstep, nstp);
00421   }
00422 #endif
00423 
00424   return succeeded;
00425 }  // end of AccurateAdvance ...........................

G4double G4MagInt_Driver::ComputeAndSetErrcon (  )  [inline]

Definition at line 74 of file G4MagIntegratorDriver.icc.

References GetPgrow(), and GetSafety().

Referenced by ReSetParameters(), SetPgrow(), and SetSafety().

00075 {
00076       errcon = std::pow(max_stepping_increase/GetSafety(),1.0/GetPgrow());
00077       return errcon;
00078 } 

G4double G4MagInt_Driver::ComputeNewStepSize ( G4double  errMaxNorm,
G4double  hstepCurrent 
)

Definition at line 754 of file G4MagIntegratorDriver.cc.

References GetPgrow(), GetPshrnk(), and GetSafety().

Referenced by AccurateAdvance(), G4ChordFinderSaf::FindNextChord(), G4ChordFinder::FindNextChord(), and QuickAdvance().

00757 {
00758   G4double hnew;
00759 
00760   // Compute size of next Step for a failed step
00761   if(errMaxNorm > 1.0 )
00762   {
00763     // Step failed; compute the size of retrial Step.
00764     hnew = GetSafety()*hstepCurrent*std::pow(errMaxNorm,GetPshrnk()) ;
00765   } else if(errMaxNorm > 0.0 ) {
00766     // Compute size of next Step for a successful step
00767     hnew = GetSafety()*hstepCurrent*std::pow(errMaxNorm,GetPgrow()) ;
00768   } else {
00769     // if error estimate is zero (possible) or negative (dubious)
00770     hnew = max_stepping_increase * hstepCurrent; 
00771   }
00772 
00773   return hnew;
00774 }

G4double G4MagInt_Driver::ComputeNewStepSize_WithinLimits ( G4double  errMaxNorm,
G4double  hstepCurrent 
)

Definition at line 784 of file G4MagIntegratorDriver.cc.

References GetPgrow(), GetPshrnk(), and GetSafety().

00787 {
00788   G4double hnew;
00789 
00790   // Compute size of next Step for a failed step
00791   if (errMaxNorm > 1.0 )
00792   {
00793     // Step failed; compute the size of retrial Step.
00794     hnew = GetSafety()*hstepCurrent*std::pow(errMaxNorm,GetPshrnk()) ;
00795   
00796     if (hnew < max_stepping_decrease*hstepCurrent)
00797     {
00798       hnew = max_stepping_decrease*hstepCurrent ;
00799                          // reduce stepsize, but no more
00800                          // than this factor (value= 1/10)
00801     }
00802   }
00803   else
00804   {
00805     // Compute size of next Step for a successful step
00806     if (errMaxNorm > errcon)
00807      { hnew = GetSafety()*hstepCurrent*std::pow(errMaxNorm,GetPgrow()); }
00808     else  // No more than a factor of 5 increase
00809      { hnew = max_stepping_increase * hstepCurrent; }
00810   }
00811   return hnew;
00812 }

void G4MagInt_Driver::GetDerivatives ( const G4FieldTrack y_curr,
G4double  dydx[] 
) [inline]

Definition at line 144 of file G4MagIntegratorDriver.icc.

References G4FieldTrack::DumpToArray(), and G4FieldTrack::ncompSVEC.

00146 { 
00147   G4double  tmpValArr[G4FieldTrack::ncompSVEC];
00148   y_curr.DumpToArray( tmpValArr  );
00149   pIntStepper -> RightHandSide( tmpValArr , dydx );
00150 }

G4double G4MagInt_Driver::GetErrcon (  )  const [inline]

Definition at line 62 of file G4MagIntegratorDriver.icc.

00063 {
00064       return errcon;
00065 }

G4double G4MagInt_Driver::GetHmin (  )  const [inline]

Definition at line 32 of file G4MagIntegratorDriver.icc.

00033 {
00034       return fMinimumStep;
00035 } 

G4int G4MagInt_Driver::GetMaxNoSteps (  )  const [inline]

Definition at line 132 of file G4MagIntegratorDriver.icc.

00133 {
00134   return fMaxNoSteps;
00135 }

G4double G4MagInt_Driver::GetPgrow (  )  const [inline]

Definition at line 56 of file G4MagIntegratorDriver.icc.

Referenced by ComputeAndSetErrcon(), ComputeNewStepSize(), ComputeNewStepSize_WithinLimits(), and OneGoodStep().

00057 {
00058       return pgrow;
00059 }

G4double G4MagInt_Driver::GetPshrnk (  )  const [inline]

Definition at line 50 of file G4MagIntegratorDriver.icc.

Referenced by ComputeNewStepSize(), ComputeNewStepSize_WithinLimits(), and OneGoodStep().

00051 {
00052       return pshrnk;
00053 } 

G4double G4MagInt_Driver::GetSafety (  )  const [inline]

Definition at line 44 of file G4MagIntegratorDriver.icc.

Referenced by ComputeAndSetErrcon(), ComputeNewStepSize(), ComputeNewStepSize_WithinLimits(), and OneGoodStep().

00045 {
00046       return safety;
00047 }

G4double G4MagInt_Driver::GetSmallestFraction (  )  const [inline]

Definition at line 165 of file G4MagIntegratorDriver.icc.

00166 {
00167       return fSmallestFraction; 
00168 } 

const G4MagIntegratorStepper * G4MagInt_Driver::GetStepper (  )  const [inline]

Definition at line 126 of file G4MagIntegratorDriver.icc.

Referenced by G4ErrorPropagatorManager::InitFieldForBackwards().

00127 {
00128   return pIntStepper;
00129 }

G4double G4MagInt_Driver::GetVerboseLevel (  )  const [inline]

Definition at line 153 of file G4MagIntegratorDriver.icc.

00154 {
00155       return fVerboseLevel;
00156 } 

G4double G4MagInt_Driver::Hmin (  )  const [inline]

Definition at line 38 of file G4MagIntegratorDriver.icc.

Referenced by AccurateAdvance(), and WarnSmallStepSize().

00039 {
00040       return fMinimumStep;
00041 }

void G4MagInt_Driver::OneGoodStep ( G4double  ystart[],
const G4double  dydx[],
G4double x,
G4double  htry,
G4double  eps,
G4double hdid,
G4double hnext 
)

Definition at line 515 of file G4MagIntegratorDriver.cc.

References G4cerr, G4endl, GetPgrow(), GetPshrnk(), GetSafety(), G4FieldTrack::ncompSVEC, and sqr().

Referenced by AccurateAdvance().

00530                                                      :
00531 // Numerical Recipes in C: The Art of Scientific Computing (NRC), Second
00532 // Edition, by William H. Press, Saul A. Teukolsky, William T.
00533 // Vetterling, and Brian P. Flannery (Cambridge University Press 1992),
00534 // 16.2 Adaptive StepSize Control for Runge-Kutta, p. 719
00535 
00536 {
00537   G4double errmax_sq;
00538   G4double h, htemp, xnew ;
00539 
00540   G4double yerr[G4FieldTrack::ncompSVEC], ytemp[G4FieldTrack::ncompSVEC];
00541 
00542   h = htry ; // Set stepsize to the initial trial value
00543 
00544   G4double inv_eps_vel_sq = 1.0 / (eps_rel_max*eps_rel_max);
00545 
00546   G4double errpos_sq=0.0;    // square of displacement error
00547   G4double errvel_sq=0.0;    // square of momentum vector difference
00548   G4double errspin_sq=0.0;   // square of spin vector difference
00549 
00550   G4int iter;
00551 
00552   static G4int tot_no_trials=0; 
00553   const G4int max_trials=100; 
00554 
00555   G4ThreeVector Spin(y[9],y[10],y[11]);
00556   G4bool     hasSpin= (Spin.mag2() > 0.0); 
00557 
00558   for (iter=0; iter<max_trials ;iter++)
00559   {
00560     tot_no_trials++;
00561     pIntStepper-> Stepper(y,dydx,h,ytemp,yerr); 
00562     //            *******
00563     G4double eps_pos = eps_rel_max * std::max(h, fMinimumStep); 
00564     G4double inv_eps_pos_sq = 1.0 / (eps_pos*eps_pos); 
00565 
00566     // Evaluate accuracy
00567     //
00568     errpos_sq =  sqr(yerr[0]) + sqr(yerr[1]) + sqr(yerr[2]) ;
00569     errpos_sq *= inv_eps_pos_sq; // Scale relative to required tolerance
00570 
00571     // Accuracy for momentum
00572     errvel_sq =  (sqr(yerr[3]) + sqr(yerr[4]) + sqr(yerr[5]) )
00573                / (sqr(y[3]) + sqr(y[4]) + sqr(y[5]) );
00574     errvel_sq *= inv_eps_vel_sq;
00575     errmax_sq = std::max( errpos_sq, errvel_sq ); // Square of maximum error
00576 
00577     if( hasSpin )
00578     { 
00579       // Accuracy for spin
00580       errspin_sq =  ( sqr(yerr[9]) + sqr(yerr[10]) + sqr(yerr[11]) )
00581                  /  ( sqr(y[9]) + sqr(y[10]) + sqr(y[11]) );
00582       errspin_sq *= inv_eps_vel_sq;
00583       errmax_sq = std::max( errmax_sq, errspin_sq ); 
00584    }
00585 
00586     if ( errmax_sq <= 1.0 )  { break; } // Step succeeded. 
00587 
00588     // Step failed; compute the size of retrial Step.
00589     htemp = GetSafety()*h* std::pow( errmax_sq, 0.5*GetPshrnk() );
00590 
00591     if (htemp >= 0.1*h)  { h = htemp; }  // Truncation error too large,
00592     else  { h = 0.1*h; }                 // reduce stepsize, but no more
00593                                          // than a factor of 10
00594     xnew = x + h;
00595     if(xnew == x)
00596     {
00597       G4cerr << "G4MagIntegratorDriver::OneGoodStep:" << G4endl
00598              << "  Stepsize underflow in Stepper " << G4endl ;
00599       G4cerr << "  Step's start x=" << x << " and end x= " << xnew 
00600              << " are equal !! " << G4endl
00601              <<"  Due to step-size= " << h 
00602              << " . Note that input step was " << htry << G4endl;
00603       break;
00604     }
00605   }
00606 
00607 #ifdef G4FLD_STATS
00608   // Sum of squares of position error // and momentum dir (underestimated)
00609   fSumH_lg += h; 
00610   fDyerrPos_lgTot += errpos_sq;
00611   fDyerrVel_lgTot += errvel_sq * h * h; 
00612 #endif
00613 
00614   // Compute size of next Step
00615   if (errmax_sq > errcon*errcon)
00616   { 
00617     hnext = GetSafety()*h*std::pow(errmax_sq, 0.5*GetPgrow());
00618   }
00619   else
00620   {
00621     hnext = max_stepping_increase*h ; // No more than a factor of 5 increase
00622   }
00623   x += (hdid = h);
00624 
00625   for(G4int k=0;k<fNoIntegrationVariables;k++) { y[k] = ytemp[k]; }
00626 
00627   return;
00628 }   // end of  OneGoodStep .............................

void G4MagInt_Driver::PrintStat_Aux ( const G4FieldTrack aFieldTrack,
G4double  requestStep,
G4double  actualStep,
G4int  subStepNo,
G4double  subStepSize,
G4double  dotVelocities 
) [protected]

Definition at line 919 of file G4MagIntegratorDriver.cc.

References G4cout, G4endl, G4FieldTrack::GetCurveLength(), G4FieldTrack::GetKineticEnergy(), G4FieldTrack::GetMomentumDir(), and G4FieldTrack::GetPosition().

Referenced by PrintStatus().

00926 {
00927     const G4ThreeVector Position=      aFieldTrack.GetPosition();
00928     const G4ThreeVector UnitVelocity=  aFieldTrack.GetMomentumDir();
00929  
00930     if( subStepNo >= 0)
00931     {
00932        G4cout << std::setw( 5) << subStepNo << " ";
00933     }
00934     else
00935     {
00936        G4cout << std::setw( 5) << "Start" << " ";
00937     }
00938     G4double curveLen= aFieldTrack.GetCurveLength();
00939     G4cout << std::setw( 7) << curveLen;
00940     G4cout << std::setw( 9) << Position.x() << " "
00941            << std::setw( 9) << Position.y() << " "
00942            << std::setw( 9) << Position.z() << " "
00943            << std::setw( 8) << UnitVelocity.x() << " "
00944            << std::setw( 8) << UnitVelocity.y() << " "
00945            << std::setw( 8) << UnitVelocity.z() << " ";
00946     G4int oldprec= G4cout.precision(3);
00947     G4cout << std::setw( 8) << UnitVelocity.mag2()-1.0 << " ";
00948     G4cout.precision(6);
00949     G4cout << std::setw(10) << dotVeloc_StartCurr << " ";
00950     G4cout.precision(oldprec);
00951     G4cout << std::setw( 7) << aFieldTrack.GetKineticEnergy();
00952     G4cout << std::setw(12) << step_len << " ";
00953 
00954     static G4double oldCurveLength= 0.0;
00955     static G4double oldSubStepLength= 0.0;
00956     static G4int oldSubStepNo= -1;
00957 
00958     G4double subStep_len=0.0;
00959     if( curveLen > oldCurveLength )
00960     {
00961       subStep_len= curveLen - oldCurveLength;
00962     }
00963     else if (subStepNo == oldSubStepNo)
00964     {
00965       subStep_len= oldSubStepLength;
00966     }
00967     oldCurveLength= curveLen;
00968     oldSubStepLength= subStep_len;
00969 
00970     G4cout << std::setw(12) << subStep_len << " "; 
00971     G4cout << std::setw(12) << subStepSize << " "; 
00972     if( requestStep != -1.0 )
00973     {
00974       G4cout << std::setw( 9) << requestStep << " ";
00975     }
00976     else
00977     {
00978        G4cout << std::setw( 9) << " InitialStep " << " ";
00979     }
00980     G4cout << G4endl;
00981 }

void G4MagInt_Driver::PrintStatisticsReport (  )  [protected]

Definition at line 985 of file G4MagIntegratorDriver.cc.

References G4cout, and G4endl.

Referenced by ~G4MagInt_Driver().

00986 {
00987   G4int noPrecBig= 6;
00988   G4int oldPrec= G4cout.precision(noPrecBig);
00989 
00990   G4cout << "G4MagInt_Driver Statistics of steps undertaken. " << G4endl;
00991   G4cout << "G4MagInt_Driver: Number of Steps: "
00992          << " Total= " <<  fNoTotalSteps
00993          << " Bad= "   <<  fNoBadSteps 
00994          << " Small= " <<  fNoSmallSteps 
00995          << " Non-initial small= " << (fNoSmallSteps-fNoInitialSmallSteps)
00996          << G4endl;
00997 
00998 #ifdef G4FLD_STATS
00999   G4cout << "MID dyerr: " 
01000          << " maximum= " << fDyerr_max 
01001          << " Sum small= " << fDyerrPos_smTot 
01002          << " std::sqrt(Sum large^2): pos= " << std::sqrt(fDyerrPos_lgTot)
01003          << " vel= " << std::sqrt( fDyerrVel_lgTot )
01004          << " Total h-distance: small= " << fSumH_sm 
01005          << " large= " << fSumH_lg
01006          << G4endl;
01007 
01008 #if 0
01009   G4int noPrecSmall=4; 
01010   // Single line precis of statistics ... optional
01011   G4cout.precision(noPrecSmall);
01012   G4cout << "MIDnums: " << fMinimumStep
01013          << "   " << fNoTotalSteps 
01014          << "  "  <<  fNoSmallSteps
01015          << "  "  << fNoSmallSteps-fNoInitialSmallSteps
01016          << "  "  << fNoBadSteps         
01017          << "   " << fDyerr_max
01018          << "   " << fDyerr_mx2 
01019          << "   " << fDyerrPos_smTot 
01020          << "   " << fSumH_sm
01021          << "   " << fDyerrPos_lgTot
01022          << "   " << fDyerrVel_lgTot
01023          << "   " << fSumH_lg
01024          << G4endl;
01025 #endif 
01026 #endif 
01027 
01028  G4cout.precision(oldPrec);
01029 }

void G4MagInt_Driver::PrintStatus ( const G4FieldTrack StartFT,
const G4FieldTrack CurrentFT,
G4double  requestStep,
G4int  subStepNo 
) [protected]

Definition at line 841 of file G4MagIntegratorDriver.cc.

References G4cout, G4endl, G4FieldTrack::GetCurveLength(), G4FieldTrack::GetMomentumDir(), G4FieldTrack::GetPosition(), and PrintStat_Aux().

00846 {
00847     G4int verboseLevel= fVerboseLevel;
00848     static G4int noPrecision= 5;
00849     G4int oldPrec= G4cout.precision(noPrecision);
00850     // G4cout.setf(ios_base::fixed,ios_base::floatfield);
00851 
00852     const G4ThreeVector StartPosition=       StartFT.GetPosition();
00853     const G4ThreeVector StartUnitVelocity=   StartFT.GetMomentumDir();
00854     const G4ThreeVector CurrentPosition=     CurrentFT.GetPosition();
00855     const G4ThreeVector CurrentUnitVelocity= CurrentFT.GetMomentumDir();
00856 
00857     G4double  DotStartCurrentVeloc= StartUnitVelocity.dot(CurrentUnitVelocity);
00858 
00859     G4double step_len= CurrentFT.GetCurveLength() - StartFT.GetCurveLength();
00860     G4double subStepSize = step_len;
00861      
00862     if( (subStepNo <= 1) || (verboseLevel > 3) )
00863     {
00864        subStepNo = - subStepNo;        // To allow printing banner
00865 
00866        G4cout << std::setw( 6)  << " " << std::setw( 25)
00867               << " G4MagInt_Driver: Current Position  and  Direction" << " "
00868               << G4endl; 
00869        G4cout << std::setw( 5) << "Step#" << " "
00870               << std::setw( 7) << "s-curve" << " "
00871               << std::setw( 9) << "X(mm)" << " "
00872               << std::setw( 9) << "Y(mm)" << " "  
00873               << std::setw( 9) << "Z(mm)" << " "
00874               << std::setw( 8) << " N_x " << " "
00875               << std::setw( 8) << " N_y " << " "
00876               << std::setw( 8) << " N_z " << " "
00877               << std::setw( 8) << " N^2-1 " << " "
00878               << std::setw(10) << " N(0).N " << " "
00879               << std::setw( 7) << "KinEner " << " "
00880               << std::setw(12) << "Track-l" << " "   // Add the Sub-step ??
00881               << std::setw(12) << "Step-len" << " " 
00882               << std::setw(12) << "Step-len" << " " 
00883               << std::setw( 9) << "ReqStep" << " "  
00884               << G4endl;
00885     }
00886 
00887     if( (subStepNo <= 0) )
00888     {
00889       PrintStat_Aux( StartFT,  requestStep, 0., 
00890                        0,        0.0,         1.0);
00891       //*************
00892     }
00893 
00894     if( verboseLevel <= 3 )
00895     {
00896       G4cout.precision(noPrecision);
00897       PrintStat_Aux( CurrentFT, requestStep, step_len, 
00898                      subStepNo, subStepSize, DotStartCurrentVeloc );
00899       //*************
00900     }
00901 
00902     else // if( verboseLevel > 3 )
00903     {
00904        //  Multi-line output
00905        
00906        // G4cout << "Current  Position is " << CurrentPosition << G4endl 
00907        //    << " and UnitVelocity is " << CurrentUnitVelocity << G4endl;
00908        // G4cout << "Step taken was " << step_len  
00909        //    << " out of PhysicalStep= " <<  requestStep << G4endl;
00910        // G4cout << "Final safety is: " << safety << G4endl;
00911        // G4cout << "Chord length = " << (CurrentPosition-StartPosition).mag()
00912        //        << G4endl << G4endl; 
00913     }
00914     G4cout.precision(oldPrec);
00915 }

void G4MagInt_Driver::PrintStatus ( const G4double StartArr,
G4double  xstart,
const G4double CurrentArr,
G4double  xcurrent,
G4double  requestStep,
G4int  subStepNo 
) [protected]

Definition at line 816 of file G4MagIntegratorDriver.cc.

References G4FieldTrack::LoadFromArray(), and G4FieldTrack::SetCurveLength().

Referenced by AccurateAdvance(), and QuickAdvance().

00822                                  :  
00823   //                                 <dydx>           - as Initial Force
00824   //                                 stepTaken(hdid)  - last step taken
00825   //                                 nextStep (hnext) - proposal for size
00826 {
00827    G4FieldTrack  StartFT(G4ThreeVector(0,0,0),
00828                  G4ThreeVector(0,0,0), 0., 0., 0., 0. );
00829    G4FieldTrack  CurrentFT (StartFT);
00830 
00831    StartFT.LoadFromArray( StartArr, fNoIntegrationVariables); 
00832    StartFT.SetCurveLength( xstart);
00833    CurrentFT.LoadFromArray( CurrentArr, fNoIntegrationVariables); 
00834    CurrentFT.SetCurveLength( xcurrent );
00835 
00836    PrintStatus(StartFT, CurrentFT, requestStep, subStepNo ); 
00837 }

G4bool G4MagInt_Driver::QuickAdvance ( G4FieldTrack y_posvel,
const G4double  dydx[],
G4double  hstep,
G4double dchord_step,
G4double dyerr_pos_sq,
G4double dyerr_mom_rel_sq 
)

Definition at line 634 of file G4MagIntegratorDriver.cc.

References FatalException, G4Exception(), and G4FieldTrack::GetPosition().

00641 {
00642   G4Exception("G4MagInt_Driver::QuickAdvance()", "GeomField0001",
00643               FatalException, "Not yet implemented."); 
00644 
00645   // Use the parameters of this method, to please compiler
00646   dchord_step = dyerr_pos_sq = hstep * hstep * dydx[0]; 
00647   dyerr_mom_rel_sq = y_posvel.GetPosition().mag2();
00648   return true;
00649 }

G4bool G4MagInt_Driver::QuickAdvance ( G4FieldTrack y_val,
const G4double  dydx[],
G4double  hstep,
G4double dchord_step,
G4double dyerr 
)

Definition at line 653 of file G4MagIntegratorDriver.cc.

References ComputeNewStepSize(), G4FieldTrack::DumpToArray(), G4cout, G4endl, G4FieldTrack::GetCurveLength(), G4FieldTrack::LoadFromArray(), G4FieldTrack::ncompSVEC, PrintStatus(), G4FieldTrack::SetCurveLength(), and sqr().

Referenced by AccurateAdvance(), G4ChordFinderSaf::FindNextChord(), and G4ChordFinder::FindNextChord().

00659 {
00660   G4double dyerr_pos_sq, dyerr_mom_rel_sq;  
00661   G4double yerr_vec[G4FieldTrack::ncompSVEC],
00662            yarrin[G4FieldTrack::ncompSVEC], yarrout[G4FieldTrack::ncompSVEC]; 
00663   G4double s_start;
00664   G4double dyerr_mom_sq, vel_mag_sq, inv_vel_mag_sq;
00665 
00666   static G4int no_call=0; 
00667   no_call ++; 
00668 
00669   // Move data into array
00670   y_posvel.DumpToArray( yarrin );      //  yarrin  <== y_posvel 
00671   s_start = y_posvel.GetCurveLength();
00672 
00673   // Do an Integration Step
00674   pIntStepper-> Stepper(yarrin, dydx, hstep, yarrout, yerr_vec) ; 
00675   //            *******
00676 
00677   // Estimate curve-chord distance
00678   dchord_step= pIntStepper-> DistChord();
00679   //                         *********
00680 
00681   // Put back the values.  yarrout ==> y_posvel
00682   y_posvel.LoadFromArray( yarrout, fNoIntegrationVariables );
00683   y_posvel.SetCurveLength( s_start + hstep );
00684 
00685 #ifdef  G4DEBUG_FIELD
00686   if(fVerboseLevel>2)
00687   {
00688     G4cout << "G4MagIntDrv: Quick Advance" << G4endl;
00689     PrintStatus( yarrin, s_start, yarrout, s_start+hstep, hstep,  1); 
00690   }
00691 #endif
00692 
00693   // A single measure of the error   
00694   //      TO-DO :  account for  energy,  spin, ... ? 
00695   vel_mag_sq   = ( sqr(yarrout[3])+sqr(yarrout[4])+sqr(yarrout[5]) );
00696   inv_vel_mag_sq = 1.0 / vel_mag_sq; 
00697   dyerr_pos_sq = ( sqr(yerr_vec[0])+sqr(yerr_vec[1])+sqr(yerr_vec[2]));
00698   dyerr_mom_sq = ( sqr(yerr_vec[3])+sqr(yerr_vec[4])+sqr(yerr_vec[5]));
00699   dyerr_mom_rel_sq = dyerr_mom_sq * inv_vel_mag_sq;
00700 
00701   // Calculate also the change in the momentum squared also ???
00702   // G4double veloc_square = y_posvel.GetVelocity().mag2();
00703   // ...
00704 
00705 #ifdef RETURN_A_NEW_STEP_LENGTH
00706   // The following step cannot be done here because "eps" is not known.
00707   dyerr_len = std::sqrt( dyerr_len_sq ); 
00708   dyerr_len_sq /= eps ;
00709 
00710   // Look at the velocity deviation ?
00711   //  sqr(yerr_vec[3])+sqr(yerr_vec[4])+sqr(yerr_vec[5]));
00712 
00713   // Set suggested new step
00714   hstep= ComputeNewStepSize( dyerr_len, hstep);
00715 #endif
00716 
00717   if( dyerr_pos_sq > ( dyerr_mom_rel_sq * sqr(hstep) ) )
00718   {
00719     dyerr = std::sqrt(dyerr_pos_sq);
00720   }
00721   else
00722   {
00723     // Scale it to the current step size - for now
00724     dyerr = std::sqrt(dyerr_mom_rel_sq) * hstep;
00725   }
00726 
00727   return true;
00728 }

void G4MagInt_Driver::RenewStepperAndAdjust ( G4MagIntegratorStepper pItsStepper  )  [inline]

Definition at line 110 of file G4MagIntegratorDriver.icc.

References ReSetParameters().

Referenced by G4MagInt_Driver().

00111 {  
00112       pIntStepper = pItsStepper; 
00113       ReSetParameters();
00114 }

void G4MagInt_Driver::ReSetParameters ( G4double  new_safety = 0.9  )  [inline]

Definition at line 81 of file G4MagIntegratorDriver.icc.

References ComputeAndSetErrcon(), and G4MagIntegratorStepper::IntegratorOrder().

Referenced by RenewStepperAndAdjust().

00082 {
00083       safety = new_safety;
00084       pshrnk = -1.0 / pIntStepper->IntegratorOrder();
00085       pgrow  = -1.0 / (1.0 + pIntStepper->IntegratorOrder());
00086       ComputeAndSetErrcon();
00087 }

void G4MagInt_Driver::SetChargeMomentumMass ( G4double  particleCharge,
G4double  MomentumXc,
G4double  Mass 
) [inline]

Definition at line 117 of file G4MagIntegratorDriver.icc.

References G4MagIntegratorStepper::GetEquationOfMotion(), and G4EquationOfMotion::SetChargeMomentumMass().

00120 { 
00121       pIntStepper->GetEquationOfMotion()
00122                  ->SetChargeMomentumMass(particleCharge, MomentumXc, Mass); 
00123 }

void G4MagInt_Driver::SetErrcon ( G4double  valEc  )  [inline]

Definition at line 104 of file G4MagIntegratorDriver.icc.

00105 { 
00106       errcon=val;
00107 }

void G4MagInt_Driver::SetHmin ( G4double  newval  )  [inline]

Definition at line 68 of file G4MagIntegratorDriver.icc.

00069 {
00070       fMinimumStep = newval;
00071 } 

void G4MagInt_Driver::SetMaxNoSteps ( G4int  val  )  [inline]

Definition at line 138 of file G4MagIntegratorDriver.icc.

00139 {
00140   fMaxNoSteps= val;
00141 }

void G4MagInt_Driver::SetPgrow ( G4double  valPg  )  [inline]

Definition at line 97 of file G4MagIntegratorDriver.icc.

References ComputeAndSetErrcon().

00098 { 
00099       pgrow=val;
00100       ComputeAndSetErrcon(); 
00101 }

void G4MagInt_Driver::SetPshrnk ( G4double  valPs  )  [inline]

void G4MagInt_Driver::SetSafety ( G4double  valS  )  [inline]

Definition at line 90 of file G4MagIntegratorDriver.icc.

References ComputeAndSetErrcon().

00091 { 
00092       safety=val;
00093       ComputeAndSetErrcon();
00094 }

void G4MagInt_Driver::SetSmallestFraction ( G4double  val  ) 

Definition at line 1033 of file G4MagIntegratorDriver.cc.

References G4cerr, and G4endl.

01034 {
01035   if( (newFraction > 1.e-16) && (newFraction < 1e-8) )
01036   {
01037     fSmallestFraction= newFraction;
01038   }
01039   else
01040   { 
01041     G4cerr << "Warning: SmallestFraction not changed. " << G4endl
01042            << "  Proposed value was " << newFraction << G4endl
01043            << "  Value must be between 1.e-8 and 1.e-16" << G4endl;
01044   }
01045 }

void G4MagInt_Driver::SetVerboseLevel ( G4int  newLevel  )  [inline]

Definition at line 159 of file G4MagIntegratorDriver.icc.

Referenced by G4PropagatorInField::SetVerboseLevel().

00160 {
00161       fVerboseLevel= newLevel;
00162 }

void G4MagInt_Driver::WarnEndPointTooFar ( G4double  endPointDist,
G4double  hStepSize,
G4double  epsilonRelative,
G4int  debugFlag 
) [protected]

Definition at line 479 of file G4MagIntegratorDriver.cc.

References G4endl, G4Exception(), G4GeometryTolerance::GetInstance(), and JustWarning.

Referenced by AccurateAdvance().

00483 {
00484   static G4double maxRelError=0.0;
00485   G4bool isNewMax, prNewMax;
00486 
00487   isNewMax = endPointDist > (1.0 + maxRelError) * h;
00488   prNewMax = endPointDist > (1.0 + 1.05 * maxRelError) * h;
00489   if( isNewMax ) { maxRelError= endPointDist / h - 1.0; }
00490 
00491   if( dbg && (h > G4GeometryTolerance::GetInstance()->GetSurfaceTolerance()) 
00492           && ( (dbg>1) || prNewMax || (endPointDist >= h*(1.+eps) ) ) )
00493   { 
00494     static G4int noWarnings = 0;
00495     std::ostringstream message;
00496     if( (noWarnings ++ < 10) || (dbg>2) )
00497     {
00498       message << "The integration produced an end-point which " << G4endl
00499               << "is further from the start-point than the curve length."
00500               << G4endl;
00501     }
00502     message << "  Distance of endpoints = " << endPointDist
00503             << ", curve length = " << h << G4endl
00504             << "  Difference (curveLen-endpDist)= " << (h - endPointDist)
00505             << ", relative = " << (h-endPointDist) / h 
00506             << ", epsilon =  " << eps;
00507     G4Exception("G4MagInt_Driver::WarnEndPointTooFar()", "GeomField1001",
00508                 JustWarning, message);
00509   }
00510 }

void G4MagInt_Driver::WarnSmallStepSize ( G4double  hnext,
G4double  hstep,
G4double  h,
G4double  xDone,
G4int  noSteps 
) [protected]

Definition at line 430 of file G4MagIntegratorDriver.cc.

References G4endl, G4Exception(), Hmin(), and JustWarning.

Referenced by AccurateAdvance().

00433 {
00434   static G4int noWarningsIssued =0;
00435   const  G4int maxNoWarnings =  10;   // Number of verbose warnings
00436   std::ostringstream message;
00437   if( (noWarningsIssued < maxNoWarnings) || fVerboseLevel > 10 )
00438   {
00439     message << "The stepsize for the next iteration, " << hnext
00440             << ", is too small - in Step number " << nstp << "." << G4endl
00441             << "The minimum for the driver is " << Hmin()  << G4endl
00442             << "Requested integr. length was " << hstep << " ." << G4endl
00443             << "The size of this sub-step was " << h     << " ." << G4endl
00444             << "The integrations has already gone " << xDone;
00445   }
00446   else
00447   {
00448     message << "Too small 'next' step " << hnext
00449             << ", step-no: " << nstp << G4endl
00450             << ", this sub-step: " << h     
00451             << ",  req_tot_len: " << hstep 
00452             << ", done: " << xDone << ", min: " << Hmin();
00453   }
00454   G4Exception("G4MagInt_Driver::WarnSmallStepSize()", "GeomField1001",
00455               JustWarning, message);
00456   noWarningsIssued++;
00457 }

void G4MagInt_Driver::WarnTooManyStep ( G4double  x1start,
G4double  x2end,
G4double  xCurrent 
) [protected]

Definition at line 462 of file G4MagIntegratorDriver.cc.

References G4endl, G4Exception(), and JustWarning.

Referenced by AccurateAdvance().

00465 {
00466     std::ostringstream message;
00467     message << "The number of steps used in the Integration driver"
00468             << " (Runge-Kutta) is too many." << G4endl
00469             << "Integration of the interval was not completed !" << G4endl
00470             << "Only a " << (xCurrent-x1start)*100/(x2end-x1start)
00471             << " % fraction of it was done.";
00472     G4Exception("G4MagInt_Driver::WarnTooManyStep()", "GeomField1001",
00473                 JustWarning, message);
00474 }


The documentation for this class was generated from the following files:
Generated on Mon May 27 17:52:28 2013 for Geant4 by  doxygen 1.4.7