#include <G4MagIntegratorDriver.hh>
Definition at line 48 of file G4MagIntegratorDriver.hh.
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 }
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 }
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] |
G4double G4MagInt_Driver::GetHmin | ( | ) | const [inline] |
G4int G4MagInt_Driver::GetMaxNoSteps | ( | ) | const [inline] |
G4double G4MagInt_Driver::GetPgrow | ( | ) | const [inline] |
Definition at line 56 of file G4MagIntegratorDriver.icc.
Referenced by ComputeAndSetErrcon(), ComputeNewStepSize(), ComputeNewStepSize_WithinLimits(), and OneGoodStep().
G4double G4MagInt_Driver::GetPshrnk | ( | ) | const [inline] |
Definition at line 50 of file G4MagIntegratorDriver.icc.
Referenced by ComputeNewStepSize(), ComputeNewStepSize_WithinLimits(), and OneGoodStep().
G4double G4MagInt_Driver::GetSafety | ( | ) | const [inline] |
Definition at line 44 of file G4MagIntegratorDriver.icc.
Referenced by ComputeAndSetErrcon(), ComputeNewStepSize(), ComputeNewStepSize_WithinLimits(), and OneGoodStep().
G4double G4MagInt_Driver::GetSmallestFraction | ( | ) | const [inline] |
const G4MagIntegratorStepper * G4MagInt_Driver::GetStepper | ( | ) | const [inline] |
Definition at line 126 of file G4MagIntegratorDriver.icc.
Referenced by G4ErrorPropagatorManager::InitFieldForBackwards().
G4double G4MagInt_Driver::GetVerboseLevel | ( | ) | const [inline] |
G4double G4MagInt_Driver::Hmin | ( | ) | const [inline] |
Definition at line 38 of file G4MagIntegratorDriver.icc.
Referenced by AccurateAdvance(), and WarnSmallStepSize().
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] |
void G4MagInt_Driver::SetHmin | ( | G4double | newval | ) | [inline] |
void G4MagInt_Driver::SetMaxNoSteps | ( | G4int | val | ) | [inline] |
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().
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 }