Geant4  10.03.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4FSALIntegrationDriver.cc
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // $Id: FSALMagIntegratorDriver.cc
27 //
28 //
29 //
30 // Implementation for class FSALMagInt_Driver
31 // Tracking in space dependent magnetic field
32 //
33 // History of major changes: /To be filled/
34 
35 
36 #include <iomanip>
37 
38 #include "globals.hh"
39 #include "G4SystemOfUnits.hh"
40 #include "G4GeometryTolerance.hh"
42 #include "G4FieldTrack.hh"
43 
44 // Stepsize can increase by no more than 5.0
45 // and decrease by no more than 1/10. = 0.1
46 //
47 const G4double G4FSALIntegrationDriver::max_stepping_increase = 5.0; //Changed from 5.0 by [hackabot]
48 const G4double G4FSALIntegrationDriver::max_stepping_decrease = 0.1; //Changed from 0.1 by [hackabot]
49 
50 // The (default) maximum number of steps is Base
51 // divided by the order of Stepper
52 //
53 const G4int G4FSALIntegrationDriver::fMaxStepBase = 100000; // Was 5000, was 250
54 
55 #ifndef G4NO_FIELD_STATISTICS
56 #define G4FLD_STATS 1
57 #endif
58 
59 //#ifndef G4DEBUG_FIELD
60 //#define G4DEBUG_FIELD 1
61 //#endif
62 
63 // ---------------------------------------------------------
64 
65 // Constructor
66 //
68  G4VFSALIntegrationStepper *pStepper,
69  G4int numComponents,
70  G4int statisticsVerbose)
71 : fSmallestFraction( 1.0e-12 ),
72 fNoIntegrationVariables(numComponents),
73 fMinNoVars(12),
74 fNoVars( std::max( fNoIntegrationVariables, fMinNoVars )),
75 fStatisticsVerboseLevel(statisticsVerbose),
76 fNoTotalSteps(0), fNoBadSteps(0), fNoSmallSteps(0),
77 fNoInitialSmallSteps(0),
78 fDyerr_max(0.0), fDyerr_mx2(0.0),
79 fDyerrPos_smTot(0.0), fDyerrPos_lgTot(0.0), fDyerrVel_lgTot(0.0),
80 fSumH_sm(0.0), fSumH_lg(0.0),
81 fVerboseLevel(0),
82 TotalNoStepperCalls(0)
83 {
84  // In order to accomodate "Laboratory Time", which is [7], fMinNoVars=8
85  // is required. For proper time of flight and spin, fMinNoVars must be 12
86 
87  RenewStepperAndAdjust( pStepper );
88  fMinimumStep= hminimum;
89  fMaxNoSteps = fMaxStepBase / pIntStepper->IntegratorOrder();
90 #ifdef G4DEBUG_FIELD
91  fVerboseLevel=2;
92 #endif
93 
94  if( (fVerboseLevel > 0) || (fStatisticsVerboseLevel > 1) )
95  {
96  G4cout << "MagIntDriver version: Accur-Adv: "
97  << "invE_nS, QuickAdv-2sqrt with Statistics "
98 #ifdef G4FLD_STATS
99  << " enabled "
100 #else
101  << " disabled "
102 #endif
103  << G4endl;
104  }
105 }
106 
107 // ---------------------------------------------------------
108 
109 // Destructor
110 //
112 {
113  if( fStatisticsVerboseLevel > 1 )
114  {
116  }
117 }
118 
119 // To add much printing for debugging purposes, uncomment the following
120 // and set verbose level to 1 or higher value !
121 // #define G4DEBUG_FIELD 1
122 
123 // ---------------------------------------------------------
124 
125 G4bool
127  G4double hstep,
128  G4double eps,
129  G4double hinitial )
130 {
131  // Runge-Kutta driver with adaptive stepsize control. Integrate starting
132  // values at y_current over hstep x2 with accuracy eps.
133  // On output ystart is replaced by values at the end of the integration
134  // interval. RightHandSide is the right-hand side of ODE system.
135  // The source is similar to odeint routine from NRC p.721-722 .
136 
137  G4int nstp, i, no_warnings=0;
138  G4double x, hnext, hdid, h;
139 
140 #ifdef G4DEBUG_FIELD
141  static G4int dbg=10;
142  static G4int nStpPr=50; // For debug printing of long integrations
143  G4double ySubStepStart[G4FieldTrack::ncompSVEC];
144  G4FieldTrack yFldTrkStart(y_current);
145 #endif
146 
149  G4double x1, x2;
150  G4bool succeeded = true, lastStepSucceeded;
151 
152  G4double startCurveLength;
153 
154  G4int noFullIntegr=0, noSmallIntegr = 0 ;
155  static G4ThreadLocal G4int noGoodSteps =0 ; // Bad = chord > curve-len
156  const G4int nvar= fNoVars;
157 
158  G4FieldTrack yStartFT(y_current);
159 
160  // Ensure that hstep > 0
161  //
162  if( hstep <= 0.0 )
163  {
164  if(hstep==0.0)
165  {
166  std::ostringstream message;
167  message << "Proposed step is zero; hstep = " << hstep << " !";
168  G4Exception("G4FSALIntegrationDriver::AccurateAdvance()",
169  "GeomField1001", JustWarning, message);
170  return succeeded;
171  }
172  else
173  {
174  std::ostringstream message;
175  message << "Invalid run condition." << G4endl
176  << "Proposed step is negative; hstep = " << hstep << "." << G4endl
177  << "Requested step cannot be negative! Aborting event.";
178  G4Exception("G4FSALIntegrationDriver::AccurateAdvance()",
179  "GeomField0003", EventMustBeAborted, message);
180  return false;
181  }
182  }
183 
184  y_current.DumpToArray( ystart );
185 
186  startCurveLength= y_current.GetCurveLength();
187  x1= startCurveLength;
188  x2= x1 + hstep;
189 
190  if ( (hinitial > 0.0) && (hinitial < hstep)
191  && (hinitial > perMillion * hstep) )
192  {
193  h = hinitial;
194  }
195  else // Initial Step size "h" defaults to the full interval
196  {
197  h = hstep;
198  }
199 
200  x = x1;
201 
202  for (i=0;i<nvar;i++) { y[i] = ystart[i]; }
203 
204  G4bool lastStep= false;
205  nstp=1;
206 
207  pIntStepper->ComputeRightHandSide( y, dydx );
208 
209 
210  do
211  {
212  G4ThreeVector StartPos( y[0], y[1], y[2] );
213 
214 #ifdef G4DEBUG_FIELD
215  G4double xSubStepStart= x;
216  for (i=0;i<nvar;i++) { ySubStepStart[i] = y[i]; }
217  yFldTrkStart.LoadFromArray(y, fNoIntegrationVariables);
218  yFldTrkStart.SetCurveLength(x);
219 #endif
220 
221 
222 
223 
224  fNoTotalSteps++;
225 
226  // Perform the Integration
227  //
228  if( h > fMinimumStep )
229  {
230  OneGoodStep(y,dydx,x,h,eps,hdid,hnext) ;
231  //--------------------------------------
232  lastStepSucceeded= (hdid == h);
233 #ifdef G4DEBUG_FIELD
234  if (dbg>2) {
235  PrintStatus( ySubStepStart, xSubStepStart, y, x, h, nstp); // Only
236  }
237 #endif
238  }
239  else
240  {
241  G4FieldTrack yFldTrk( G4ThreeVector(0,0,0),
242  G4ThreeVector(0,0,0), 0., 0., 0., 0. );
243  G4double dchord_step, dyerr, dyerr_len; // What to do with these ?
244  yFldTrk.LoadFromArray(y, fNoIntegrationVariables);
245  yFldTrk.SetCurveLength( x );
246 
247  QuickAdvance( yFldTrk, dydx, h, dchord_step, dyerr_len );
248  //-----------------------------------------------------
249 
250  yFldTrk.DumpToArray(y);
251 
252 #ifdef G4FLD_STATS
253  fNoSmallSteps++;
254  if ( dyerr_len > fDyerr_max) { fDyerr_max= dyerr_len; }
255  fDyerrPos_smTot += dyerr_len;
256  fSumH_sm += h; // Length total for 'small' steps
257  if (nstp<=1) { fNoInitialSmallSteps++; }
258 #endif
259 #ifdef G4DEBUG_FIELD
260  if (dbg>1)
261  {
262  if(fNoSmallSteps<2) { PrintStatus(ySubStepStart, x1, y, x, h, -nstp); }
263  G4cout << "Another sub-min step, no " << fNoSmallSteps
264  << " of " << fNoTotalSteps << " this time " << nstp << G4endl;
265  PrintStatus( ySubStepStart, x1, y, x, h, nstp); // Only this
266  G4cout << " dyerr= " << dyerr_len << " relative = " << dyerr_len / h
267  << " epsilon= " << eps << " hstep= " << hstep
268  << " h= " << h << " hmin= " << fMinimumStep << G4endl;
269  }
270 #endif
271  if( h == 0.0 )
272  {
273  G4Exception("G4FSALIntegrationDriver::AccurateAdvance()",
274  "GeomField0003", FatalException,
275  "Integration Step became Zero!");
276  }
277  dyerr = dyerr_len / h;
278  hdid= h;
279  x += hdid;
280 
281  // Compute suggested new step
282  hnext= ComputeNewStepSize( dyerr/eps, h);
283 
284  // .. hnext= ComputeNewStepSize_WithinLimits( dyerr/eps, h);
285  lastStepSucceeded= (dyerr<= eps);
286  }
287 
288  if (lastStepSucceeded) { noFullIntegr++; }
289  else { noSmallIntegr++; }
290 
291  G4ThreeVector EndPos( y[0], y[1], y[2] );
292 
293 #ifdef G4DEBUG_FIELD
294  if( (dbg>0) && (dbg<=2) && (nstp>nStpPr))
295  {
296  if( nstp==nStpPr ) { G4cout << "***** Many steps ****" << G4endl; }
297  G4cout << "MagIntDrv: " ;
298  G4cout << "hdid=" << std::setw(12) << hdid << " "
299  << "hnext=" << std::setw(12) << hnext << " "
300  << "hstep=" << std::setw(12) << hstep << " (requested) "
301  << G4endl;
302  PrintStatus( ystart, x1, y, x, h, (nstp==nStpPr) ? -nstp: nstp);
303  }
304 #endif
305 
306  // Check the endpoint
307  G4double endPointDist= (EndPos-StartPos).mag();
308  if ( endPointDist >= hdid*(1.+perMillion) )
309  {
310  fNoBadSteps++;
311 
312  // Issue a warning only for gross differences -
313  // we understand how small difference occur.
314  if ( endPointDist >= hdid*(1.+perThousand) )
315  {
316 #ifdef G4DEBUG_FIELD
317  if (dbg)
318  {
319  WarnEndPointTooFar ( endPointDist, hdid, eps, dbg );
320  G4cerr << " Total steps: bad " << fNoBadSteps
321  << " good " << noGoodSteps << " current h= " << hdid
322  << G4endl;
323  PrintStatus( ystart, x1, y, x, hstep, no_warnings?nstp:-nstp);
324  }
325 #endif
326  no_warnings++;
327  }
328  }
329  else
330  {
331  noGoodSteps ++;
332  }
333  // #endif
334 
335  // Avoid numerous small last steps
336  if( (h < eps * hstep) || (h < fSmallestFraction * startCurveLength) )
337  {
338  // No more integration -- the next step will not happen
339  lastStep = true;
340  }
341  else
342  {
343  // Check the proposed next stepsize
344  if(std::fabs(hnext) <= Hmin())
345  {
346 #ifdef G4DEBUG_FIELD
347  // If simply a very small interval is being integrated, do not warn
348  if( (x < x2 * (1-eps) ) && // The last step can be small: OK
349  (std::fabs(hstep) > Hmin()) ) // and if we are asked, it's OK
350  {
351  if(dbg>0)
352  {
353  WarnSmallStepSize( hnext, hstep, h, x-x1, nstp );
354  PrintStatus( ystart, x1, y, x, hstep, no_warnings?nstp:-nstp);
355  }
356  no_warnings++;
357  }
358 #endif
359  // Make sure that the next step is at least Hmin.
360  h = Hmin();
361  }
362  else
363  {
364  h = hnext;
365  }
366 
367  // Ensure that the next step does not overshoot
368  if ( x+h > x2 )
369  { // When stepsize overshoots, decrease it!
370  h = x2 - x ; // Must cope with difficult rounding-error
371  } // issues if hstep << x2
372 
373  if ( h == 0.0 )
374  {
375  // Cannot progress - accept this as last step - by default
376  lastStep = true;
377 #ifdef G4DEBUG_FIELD
378  if (dbg>2)
379  {
380  int prec= G4cout.precision(12);
381  G4cout << "Warning: FSALMagIntegratorDriver::AccurateAdvance"
382  << G4endl
383  << " Integration step 'h' became "
384  << h << " due to roundoff. " << G4endl
385  << " Calculated as difference of x2= "<< x2 << " and x=" << x
386  << " Forcing termination of advance." << G4endl;
387  G4cout.precision(prec);
388  }
389 #endif
390  }
391  }
392  } while ( ((nstp++)<=fMaxNoSteps) && (x < x2) && (!lastStep) );
393  // Loop checking, 07.10.2016, J. Apostolakis
394 
395  // Have we reached the end ?
396  // --> a better test might be x-x2 > an_epsilon
397 
398  succeeded= (x>=x2); // If it was a "forced" last step
399 
400  for (i=0;i<nvar;i++) { yEnd[i] = y[i]; }
401 
402  // Put back the values.
403  y_current.LoadFromArray( yEnd, fNoIntegrationVariables );
404  y_current.SetCurveLength( x );
405 
406  if(nstp > fMaxNoSteps)
407  {
408  no_warnings++;
409  succeeded = false;
410 #ifdef G4DEBUG_FIELD
411  if (dbg)
412  {
413  WarnTooManyStep( x1, x2, x ); // Issue WARNING
414  PrintStatus( yEnd, x1, y, x, hstep, -nstp);
415  }
416 #endif
417  }
418 
419 #ifdef G4DEBUG_FIELD
420  if( dbg && no_warnings )
421  {
422  G4cerr << "FSALMagIntegratorDriver exit status: no-steps " << nstp <<G4endl;
423  PrintStatus( yEnd, x1, y, x, hstep, nstp);
424  }
425 #endif
426 
427  return succeeded;
428 } // end of AccurateAdvance ...........................
429 
430 // ---------------------------------------------------------
431 
432 void
434  G4double h, G4double xDone,
435  G4int nstp)
436 {
437  static G4ThreadLocal G4int noWarningsIssued =0;
438  const G4int maxNoWarnings = 10; // Number of verbose warnings
439  std::ostringstream message;
440  if( (noWarningsIssued < maxNoWarnings) || fVerboseLevel > 10 )
441  {
442  message << "The stepsize for the next iteration, " << hnext
443  << ", is too small - in Step number " << nstp << "." << G4endl
444  << "The minimum for the driver is " << Hmin() << G4endl
445  << "Requested integr. length was " << hstep << " ." << G4endl
446  << "The size of this sub-step was " << h << " ." << G4endl
447  << "The integrations has already gone " << xDone;
448  }
449  else
450  {
451  message << "Too small 'next' step " << hnext
452  << ", step-no: " << nstp << G4endl
453  << ", this sub-step: " << h
454  << ", req_tot_len: " << hstep
455  << ", done: " << xDone << ", min: " << Hmin();
456  }
457  G4Exception("G4FSALIntegrationDriver::WarnSmallStepSize()", "GeomField1001",
458  JustWarning, message);
459  noWarningsIssued++;
460 }
461 
462 // ---------------------------------------------------------
463 
464 void
466  G4double x2end,
467  G4double xCurrent)
468 {
469  std::ostringstream message;
470  message << "The number of steps used in the Integration driver"
471  << " (Runge-Kutta) is too many." << G4endl
472  << "Integration of the interval was not completed !" << G4endl
473  << "Only a " << (xCurrent-x1start)*100/(x2end-x1start)
474  << " % fraction of it was done.";
475  G4Exception("G4FSALIntegrationDriver::WarnTooManyStep()", "GeomField1001",
476  JustWarning, message);
477 }
478 
479 // ---------------------------------------------------------
480 
481 void
483  G4double h ,
484  G4double eps,
485  G4int dbg)
486 {
487  static G4ThreadLocal G4double maxRelError=0.0;
488  G4bool isNewMax, prNewMax;
489 
490  isNewMax = endPointDist > (1.0 + maxRelError) * h;
491  prNewMax = endPointDist > (1.0 + 1.05 * maxRelError) * h;
492  if( isNewMax ) { maxRelError= endPointDist / h - 1.0; }
493 
495  && ( (dbg>1) || prNewMax || (endPointDist >= h*(1.+eps) ) ) )
496  {
497  static G4ThreadLocal G4int noWarnings = 0;
498  std::ostringstream message;
499  if( (noWarnings ++ < 10) || (dbg>2) )
500  {
501  message << "The integration produced an end-point which " << G4endl
502  << "is further from the start-point than the curve length."
503  << G4endl;
504  }
505  message << " Distance of endpoints = " << endPointDist
506  << ", curve length = " << h << G4endl
507  << " Difference (curveLen-endpDist)= " << (h - endPointDist)
508  << ", relative = " << (h-endPointDist) / h
509  << ", epsilon = " << eps;
510  G4Exception("G4FSALIntegrationDriver::WarnEndPointTooFar()", "GeomField1001",
511  JustWarning, message);
512  }
513 }
514 
515 // ---------------------------------------------------------
516 
517 void
519  G4double dydx[],
520  G4double& x, // InOut
521  G4double htry,
522  G4double eps_rel_max,
523  G4double& hdid, // Out
524  G4double& hnext ) // Out
525 
526 // Driver for one Runge-Kutta Step with monitoring of local truncation error
527 // to ensure accuracy and adjust stepsize. Input are dependent variable
528 // array y[0,...,5] and its derivative dydx[0,...,5] at the
529 // starting value of the independent variable x . Also input are stepsize
530 // to be attempted htry, and the required accuracy eps. On output y and x
531 // are replaced by their new values, hdid is the stepsize that was actually
532 // accomplished, and hnext is the estimated next stepsize.
533 // This is similar to the function rkqs from the book:
534 // Numerical Recipes in C: The Art of Scientific Computing (NRC), Second
535 // Edition, by William H. Press, Saul A. Teukolsky, William T.
536 // Vetterling, and Brian P. Flannery (Cambridge University Press 1992),
537 // 16.2 Adaptive StepSize Control for Runge-Kutta, p. 719
538 
539 {
540  G4double errmax_sq;
541  G4double h, htemp, xnew ;
542 
544 
545  h = htry ; // Set stepsize to the initial trial value
546 
547  G4double inv_eps_vel_sq = 1.0 / (eps_rel_max*eps_rel_max);
548 
549  G4double errpos_sq=0.0; // square of displacement error
550  G4double errvel_sq=0.0; // square of momentum vector difference
551  G4double errspin_sq=0.0; // square of spin vector difference
552 
553  G4int iter;
554 
555  static G4ThreadLocal G4int tot_no_trials=0;
556  //-----------------------------------------------------------
557  //Made a change :
558  //max_trials = 100
559  const G4int max_trials=100000;
560  //[hackabot]
561  //------------------------------------------------------------
562  G4ThreeVector Spin(y[9],y[10],y[11]);
563  G4double spin_mag2 =Spin.mag2() ;
564  G4bool hasSpin= (spin_mag2 > 0.0);
565 
566  for (iter=0; iter<max_trials ;iter++)
567  {
568  tot_no_trials++;
569  pIntStepper-> Stepper(y,dydx,h,ytemp,yerr, dydx);
570  // *******
571  TotalNoStepperCalls++;
572  G4double eps_pos = eps_rel_max * std::max(h, fMinimumStep);
573  G4double inv_eps_pos_sq = 1.0 / (eps_pos*eps_pos);
574 
575  // Evaluate accuracy
576  //
577  errpos_sq = sqr(yerr[0]) + sqr(yerr[1]) + sqr(yerr[2]) ;
578  errpos_sq *= inv_eps_pos_sq; // Scale relative to required tolerance
579 
580  // Accuracy for momentum
581  G4double magvel_sq= sqr(y[3]) + sqr(y[4]) + sqr(y[5]) ;
582  G4double sumerr_sq = sqr(yerr[3]) + sqr(yerr[4]) + sqr(yerr[5]) ;
583  if( magvel_sq > 0.0 ) {
584  errvel_sq = sumerr_sq / magvel_sq;
585  }else{
586  G4cerr << "** G4MagIntegrationDriver: found case of zero momentum."
587  << " iteration= " << iter << " h= " << h << G4endl;
588  errvel_sq = sumerr_sq;
589  }
590  errvel_sq *= inv_eps_vel_sq;
591  errmax_sq = std::max( errpos_sq, errvel_sq ); // Square of maximum error
592 
593  if( hasSpin )
594  {
595  // Accuracy for spin
596  errspin_sq = ( sqr(yerr[9]) + sqr(yerr[10]) + sqr(yerr[11]) )
597  / spin_mag2; // ( sqr(y[9]) + sqr(y[10]) + sqr(y[11]) );
598  errspin_sq *= inv_eps_vel_sq;
599  errmax_sq = std::max( errmax_sq, errspin_sq );
600  }
601 
602  if ( errmax_sq <= 1.0 ) { break; } // Step succeeded.
603 
604  // Step failed; compute the size of retrial Step.
605  htemp = GetSafety()*h* std::pow( errmax_sq, 0.5*GetPshrnk() );
606 
607  if (htemp >= 0.1*h) { h = htemp; } // Truncation error too large,
608  else { h = 0.1*h; } // reduce stepsize, but no more
609  // than a factor of 10
610  xnew = x + h;
611  if(xnew == x)
612  {
613  G4cerr << "FSALMagIntegratorDriver::OneGoodStep:" << G4endl
614  << " Stepsize underflow in Stepper " << G4endl ;
615  G4cerr << " Step's start x=" << x << " and end x= " << xnew
616  << " are equal !! " << G4endl
617  <<" Due to step-size= " << h
618  << " . Note that input step was " << htry << G4endl;
619  break;
620  }
621  }
622 
623 #ifdef G4FLD_STATS
624  // Sum of squares of position error // and momentum dir (underestimated)
625  fSumH_lg += h;
626  fDyerrPos_lgTot += errpos_sq;
627  fDyerrVel_lgTot += errvel_sq * h * h;
628 #endif
629 
630  // Compute size of next Step
631  if (errmax_sq > errcon*errcon)
632  {
633  hnext = GetSafety()*h*std::pow(errmax_sq, 0.5*GetPgrow());
634  }
635  else
636  {
637  hnext = max_stepping_increase*h ; // No more than a factor of 5 increase
638  }
639  x += (hdid = h);
640 
641  for(G4int k=0;k<fNoIntegrationVariables;k++) { y[k] = ytemp[k]; }
642 
643  return;
644 } // end of OneGoodStep .............................
645 
646 //----------------------------------------------------------------------
647 
648 // QuickAdvance just tries one Step - it does not ensure accuracy
649 //
651  G4FieldTrack& y_posvel, // INOUT
652  G4double dydx[],
653  G4double hstep, // In
654  G4double& dchord_step,
655  G4double& dyerr_pos_sq,
656  G4double& dyerr_mom_rel_sq )
657 {
658  G4Exception("G4FSALIntegrationDriver::QuickAdvance()", "GeomField0001",
659  FatalException, "Not yet implemented.");
660 
661  // Use the parameters of this method, to please compiler
662  dchord_step = dyerr_pos_sq = hstep * hstep * dydx[0];
663  dyerr_mom_rel_sq = y_posvel.GetPosition().mag2();
664  return true;
665 }
666 
667 //----------------------------------------------------------------------
668 
670  G4FieldTrack& y_posvel, // INOUT
671  G4double dydx[],
672  G4double hstep, // In
673  G4double& dchord_step,
674  G4double& dyerr )
675 {
676  G4double dyerr_pos_sq, dyerr_mom_rel_sq;
679  G4double s_start;
680  G4double dyerr_mom_sq, vel_mag_sq, inv_vel_mag_sq;
681 
682  static G4ThreadLocal G4int no_call=0;
683  no_call ++;
684 
685  // Move data into array
686  y_posvel.DumpToArray( yarrin ); // yarrin <== y_posvel
687  s_start = y_posvel.GetCurveLength();
688 
689  // Do an Integration Step
690  pIntStepper-> Stepper(yarrin, dydx, hstep, yarrout, yerr_vec, dydx) ;
691  // *******
692 
693  // Estimate curve-chord distance
694  dchord_step= pIntStepper-> DistChord();
695  // *********
696 
697  // Put back the values. yarrout ==> y_posvel
698  y_posvel.LoadFromArray( yarrout, fNoIntegrationVariables );
699  y_posvel.SetCurveLength( s_start + hstep );
700 
701 #ifdef G4DEBUG_FIELD
702  if(fVerboseLevel>2)
703  {
704  G4cout << "G4MagIntDrv: Quick Advance" << G4endl;
705  PrintStatus( yarrin, s_start, yarrout, s_start+hstep, hstep, 1);
706  }
707 #endif
708 
709  // A single measure of the error
710  // TO-DO : account for energy, spin, ... ?
711  vel_mag_sq = ( sqr(yarrout[3])+sqr(yarrout[4])+sqr(yarrout[5]) );
712  inv_vel_mag_sq = 1.0 / vel_mag_sq;
713  dyerr_pos_sq = ( sqr(yerr_vec[0])+sqr(yerr_vec[1])+sqr(yerr_vec[2]));
714  dyerr_mom_sq = ( sqr(yerr_vec[3])+sqr(yerr_vec[4])+sqr(yerr_vec[5]));
715  dyerr_mom_rel_sq = dyerr_mom_sq * inv_vel_mag_sq;
716 
717  // Calculate also the change in the momentum squared also ???
718  // G4double veloc_square = y_posvel.GetVelocity().mag2();
719  // ...
720 
721 #ifdef RETURN_A_NEW_STEP_LENGTH
722  // The following step cannot be done here because "eps" is not known.
723  dyerr_len = std::sqrt( dyerr_len_sq );
724  dyerr_len_sq /= eps ;
725 
726  // Look at the velocity deviation ?
727  // sqr(yerr_vec[3])+sqr(yerr_vec[4])+sqr(yerr_vec[5]));
728 
729  // Set suggested new step
730  hstep= ComputeNewStepSize( dyerr_len, hstep);
731 #endif
732 
733  if( dyerr_pos_sq > ( dyerr_mom_rel_sq * sqr(hstep) ) )
734  {
735  dyerr = std::sqrt(dyerr_pos_sq);
736  }
737  else
738  {
739  // Scale it to the current step size - for now
740  dyerr = std::sqrt(dyerr_mom_rel_sq) * hstep;
741  }
742 
743  return true;
744 }
745 
746 // --------------------------------------------------------------------------
747 
748 #ifdef QUICK_ADV_ARRAY_IN_AND_OUT
750  G4double yarrin[], // In
751  const G4double dydx[],
752  G4double hstep, // In
753  G4double yarrout[],
754  G4double& dchord_step,
755  G4double& dyerr ) // In length
756 {
757  G4Exception("G4FSALIntegrationDriver::QuickAdvance()", "GeomField0001",
758  FatalException, "Not yet implemented.");
759  dyerr = dchord_step = hstep * yarrin[0] * dydx[0];
760  yarrout[0]= yarrin[0];
761 }
762 #endif
763 
764 // --------------------------------------------------------------------------
765 
766 // This method computes new step sizes - but does not limit changes to
767 // within certain factors
768 //
769 G4double
771  G4double errMaxNorm, // max error (normalised)
772  G4double hstepCurrent) // current step size
773 {
774  G4double hnew;
775 
776  // Compute size of next Step for a failed step
777  if(errMaxNorm > 1.0 )
778  {
779  // Step failed; compute the size of retrial Step.
780  hnew = GetSafety()*hstepCurrent*std::pow(errMaxNorm,GetPshrnk()) ;
781  } else if(errMaxNorm > 0.0 ) {
782  // Compute size of next Step for a successful step
783  hnew = GetSafety()*hstepCurrent*std::pow(errMaxNorm,GetPgrow()) ;
784  } else {
785  // if error estimate is zero (possible) or negative (dubious)
786  hnew = max_stepping_increase * hstepCurrent;
787  }
788 
789  return hnew;
790 }
791 
792 // ---------------------------------------------------------------------------
793 
794 // This method computes new step sizes limiting changes within certain factors
795 //
796 // It shares its logic with AccurateAdvance.
797 // They are kept separate currently for optimisation.
798 //
799 G4double
801  G4double errMaxNorm, // max error (normalised)
802  G4double hstepCurrent) // current step size
803 {
804  G4double hnew;
805 
806  // Compute size of next Step for a failed step
807  if (errMaxNorm > 1.0 )
808  {
809  // Step failed; compute the size of retrial Step.
810  hnew = GetSafety()*hstepCurrent*std::pow(errMaxNorm,GetPshrnk()) ;
811 
812  if (hnew < max_stepping_decrease*hstepCurrent)
813  {
814  hnew = max_stepping_decrease*hstepCurrent ;
815  // reduce stepsize, but no more
816  // than this factor (value= 1/10)
817  }
818  }
819  else
820  {
821  // Compute size of next Step for a successful step
822  if (errMaxNorm > errcon)
823  { hnew = GetSafety()*hstepCurrent*std::pow(errMaxNorm,GetPgrow()); }
824  else // No more than a factor of 5 increase
825  { hnew = max_stepping_increase * hstepCurrent; }
826  }
827  return hnew;
828 }
829 
830 // ---------------------------------------------------------------------------
831 
833  G4double xstart,
834  const G4double* CurrentArr,
835  G4double xcurrent,
836  G4double requestStep,
837  G4int subStepNo)
838 // Potentially add as arguments:
839 // <dydx> - as Initial Force
840 // stepTaken(hdid) - last step taken
841 // nextStep (hnext) - proposal for size
842 {
843  G4FieldTrack StartFT(G4ThreeVector(0,0,0),
844  G4ThreeVector(0,0,0), 0., 0., 0., 0. );
845  G4FieldTrack CurrentFT (StartFT);
846 
847  StartFT.LoadFromArray( StartArr, fNoIntegrationVariables);
848  StartFT.SetCurveLength( xstart);
849  CurrentFT.LoadFromArray( CurrentArr, fNoIntegrationVariables);
850  CurrentFT.SetCurveLength( xcurrent );
851 
852  PrintStatus(StartFT, CurrentFT, requestStep, subStepNo );
853 }
854 
855 // ---------------------------------------------------------------------------
856 
858  const G4FieldTrack& StartFT,
859  const G4FieldTrack& CurrentFT,
860  G4double requestStep,
861  G4int subStepNo)
862 {
863  G4int verboseLevel= fVerboseLevel;
864  static G4ThreadLocal G4int noPrecision= 5;
865  G4int oldPrec= G4cout.precision(noPrecision);
866  // G4cout.setf(ios_base::fixed,ios_base::floatfield);
867 
868  const G4ThreeVector StartPosition= StartFT.GetPosition();
869  const G4ThreeVector StartUnitVelocity= StartFT.GetMomentumDir();
870  const G4ThreeVector CurrentPosition= CurrentFT.GetPosition();
871  const G4ThreeVector CurrentUnitVelocity= CurrentFT.GetMomentumDir();
872 
873  G4double DotStartCurrentVeloc= StartUnitVelocity.dot(CurrentUnitVelocity);
874 
875  G4double step_len= CurrentFT.GetCurveLength() - StartFT.GetCurveLength();
876  G4double subStepSize = step_len;
877 
878  if( (subStepNo <= 1) || (verboseLevel > 3) )
879  {
880  subStepNo = - subStepNo; // To allow printing banner
881 
882  G4cout << std::setw( 6) << " " << std::setw( 25)
883  << " G4FSALIntegrationDriver: Current Position and Direction" << " "
884  << G4endl;
885  G4cout << std::setw( 5) << "Step#" << " "
886  << std::setw( 7) << "s-curve" << " "
887  << std::setw( 9) << "X(mm)" << " "
888  << std::setw( 9) << "Y(mm)" << " "
889  << std::setw( 9) << "Z(mm)" << " "
890  << std::setw( 8) << " N_x " << " "
891  << std::setw( 8) << " N_y " << " "
892  << std::setw( 8) << " N_z " << " "
893  << std::setw( 8) << " N^2-1 " << " "
894  << std::setw(10) << " N(0).N " << " "
895  << std::setw( 7) << "KinEner " << " "
896  << std::setw(12) << "Track-l" << " " // Add the Sub-step ??
897  << std::setw(12) << "Step-len" << " "
898  << std::setw(12) << "Step-len" << " "
899  << std::setw( 9) << "ReqStep" << " "
900  << G4endl;
901  }
902 
903  if( (subStepNo <= 0) )
904  {
905  PrintStat_Aux( StartFT, requestStep, 0.,
906  0, 0.0, 1.0);
907  //*************
908  }
909 
910  if( verboseLevel <= 3 )
911  {
912  G4cout.precision(noPrecision);
913  PrintStat_Aux( CurrentFT, requestStep, step_len,
914  subStepNo, subStepSize, DotStartCurrentVeloc );
915  //*************
916  }
917 
918  else // if( verboseLevel > 3 )
919  {
920  // Multi-line output
921 
922  // G4cout << "Current Position is " << CurrentPosition << G4endl
923  // << " and UnitVelocity is " << CurrentUnitVelocity << G4endl;
924  // G4cout << "Step taken was " << step_len
925  // << " out of PhysicalStep= " << requestStep << G4endl;
926  // G4cout << "Final safety is: " << safety << G4endl;
927  // G4cout << "Chord length = " << (CurrentPosition-StartPosition).mag()
928  // << G4endl << G4endl;
929  }
930  G4cout.precision(oldPrec);
931 }
932 
933 // ---------------------------------------------------------------------------
934 
936  const G4FieldTrack& aFieldTrack,
937  G4double requestStep,
938  G4double step_len,
939  G4int subStepNo,
940  G4double subStepSize,
941  G4double dotVeloc_StartCurr)
942 {
943  const G4ThreeVector Position= aFieldTrack.GetPosition();
944  const G4ThreeVector UnitVelocity= aFieldTrack.GetMomentumDir();
945 
946  if( subStepNo >= 0)
947  {
948  G4cout << std::setw( 5) << subStepNo << " ";
949  }
950  else
951  {
952  G4cout << std::setw( 5) << "Start" << " ";
953  }
954  G4double curveLen= aFieldTrack.GetCurveLength();
955  G4cout << std::setw( 7) << curveLen;
956  G4cout << std::setw( 9) << Position.x() << " "
957  << std::setw( 9) << Position.y() << " "
958  << std::setw( 9) << Position.z() << " "
959  << std::setw( 8) << UnitVelocity.x() << " "
960  << std::setw( 8) << UnitVelocity.y() << " "
961  << std::setw( 8) << UnitVelocity.z() << " ";
962  G4int oldprec= G4cout.precision(3);
963  G4cout << std::setw( 8) << UnitVelocity.mag2()-1.0 << " ";
964  G4cout.precision(6);
965  G4cout << std::setw(10) << dotVeloc_StartCurr << " ";
966  G4cout.precision(oldprec);
967  G4cout << std::setw( 7) << aFieldTrack.GetKineticEnergy();
968  G4cout << std::setw(12) << step_len << " ";
969 
970  static G4ThreadLocal G4double oldCurveLength= 0.0;
971  static G4ThreadLocal G4double oldSubStepLength= 0.0;
972  static G4ThreadLocal G4int oldSubStepNo= -1;
973 
974  G4double subStep_len=0.0;
975  if( curveLen > oldCurveLength )
976  {
977  subStep_len= curveLen - oldCurveLength;
978  }
979  else if (subStepNo == oldSubStepNo)
980  {
981  subStep_len= oldSubStepLength;
982  }
983  oldCurveLength= curveLen;
984  oldSubStepLength= subStep_len;
985 
986  G4cout << std::setw(12) << subStep_len << " ";
987  G4cout << std::setw(12) << subStepSize << " ";
988  if( requestStep != -1.0 )
989  {
990  G4cout << std::setw( 9) << requestStep << " ";
991  }
992  else
993  {
994  G4cout << std::setw( 9) << " InitialStep " << " ";
995  }
996  G4cout << G4endl;
997 }
998 
999 // ---------------------------------------------------------------------------
1000 
1002 {
1003  G4int noPrecBig= 6;
1004  G4int oldPrec= G4cout.precision(noPrecBig);
1005 
1006  G4cout << "G4FSALIntegrationDriver Statistics of steps undertaken. " << G4endl;
1007  G4cout << "G4FSALIntegrationDriver: Number of Steps: "
1008  << " Total= " << fNoTotalSteps
1009  << " Bad= " << fNoBadSteps
1010  << " Small= " << fNoSmallSteps
1011  << " Non-initial small= " << (fNoSmallSteps-fNoInitialSmallSteps)
1012  << G4endl;
1013 
1014 #ifdef G4FLD_STATS
1015  G4cout << "MID dyerr: "
1016  << " maximum= " << fDyerr_max
1017  << " Sum small= " << fDyerrPos_smTot
1018  << " std::sqrt(Sum large^2): pos= " << std::sqrt(fDyerrPos_lgTot)
1019  << " vel= " << std::sqrt( fDyerrVel_lgTot )
1020  << " Total h-distance: small= " << fSumH_sm
1021  << " large= " << fSumH_lg
1022  << G4endl;
1023 
1024 #if 0
1025  G4int noPrecSmall=4;
1026  // Single line precis of statistics ... optional
1027  G4cout.precision(noPrecSmall);
1028  G4cout << "MIDnums: " << fMinimumStep
1029  << " " << fNoTotalSteps
1030  << " " << fNoSmallSteps
1031  << " " << fNoSmallSteps-fNoInitialSmallSteps
1032  << " " << fNoBadSteps
1033  << " " << fDyerr_max
1034  << " " << fDyerr_mx2
1035  << " " << fDyerrPos_smTot
1036  << " " << fSumH_sm
1037  << " " << fDyerrPos_lgTot
1038  << " " << fDyerrVel_lgTot
1039  << " " << fSumH_lg
1040  << G4endl;
1041 #endif
1042 #endif
1043 
1044  G4cout.precision(oldPrec);
1045 }
1046 
1047 // ---------------------------------------------------------------------------
1048 
1050 {
1051  if( (newFraction > 1.e-16) && (newFraction < 1e-8) )
1052  {
1053  fSmallestFraction= newFraction;
1054  }
1055  else
1056  {
1057  G4cerr << "Warning: SmallestFraction not changed. " << G4endl
1058  << " Proposed value was " << newFraction << G4endl
1059  << " Value must be between 1.e-8 and 1.e-16" << G4endl;
1060  }
1061 }
1062 
1063 
1064 
1065 
1066 
1067 
1068 
1069 //--------------------------------------------------------------------------------
1070 // New functions introduced for more debugging purposes :
1071 
1072 
1073 
static constexpr double perMillion
Definition: G4SIunits.hh:334
void WarnEndPointTooFar(G4double endPointDist, G4double hStepSize, G4double epsilonRelative, G4int debugFlag)
G4double GetCurveLength() const
CLHEP::Hep3Vector G4ThreeVector
void WarnSmallStepSize(G4double hnext, G4double hstep, G4double h, G4double xDone, G4int noSteps)
void PrintStat_Aux(const G4FieldTrack &aFieldTrack, G4double requestStep, G4double actualStep, G4int subStepNo, G4double subStepSize, G4double dotVelocities)
double x() const
G4double GetKineticEnergy() const
double dot(const Hep3Vector &) const
void SetCurveLength(G4double nCurve_s)
G4double ComputeNewStepSize_WithinLimits(G4double errMaxNorm, G4double hstepCurrent)
G4double GetSurfaceTolerance() const
const G4ThreeVector & GetMomentumDir() const
static const G4double eps
G4bool AccurateAdvance(G4FieldTrack &y_current, G4double hstep, G4double eps, G4double hinitial=0.0)
tuple x
Definition: test.py:50
#define G4ThreadLocal
Definition: tls.hh:89
int G4int
Definition: G4Types.hh:78
double z() const
void WarnTooManyStep(G4double x1start, G4double x2end, G4double xCurrent)
G4double GetPshrnk() const
static const double prec
Definition: RanecuEngine.cc:58
G4ThreeVector GetPosition() const
G4GLOB_DLL std::ostream G4cout
G4FSALIntegrationDriver(G4double hminimum, G4VFSALIntegrationStepper *pItsStepper, G4int numberOfComponents=6, G4int statisticsVerbosity=1)
bool G4bool
Definition: G4Types.hh:79
void RenewStepperAndAdjust(G4VFSALIntegrationStepper *pItsStepper)
void SetSmallestFraction(G4double val)
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
void OneGoodStep(G4double ystart[], G4double dydx[], G4double &x, G4double htry, G4double eps, G4double &hdid, G4double &hnext)
void DumpToArray(G4double valArr[ncompSVEC]) const
G4double ComputeNewStepSize(G4double errMaxNorm, G4double hstepCurrent)
T max(const T t1, const T t2)
brief Return the largest of the two arguments
virtual G4int IntegratorOrder() const =0
void LoadFromArray(const G4double valArr[ncompSVEC], G4int noVarsIntegrated)
double y() const
double mag2() const
G4double Hmin() const
#define G4endl
Definition: G4ios.hh:61
G4bool QuickAdvance(G4FieldTrack &y_val, G4double dydx[], G4double hstep, G4double &dchord_step, G4double &dyerr)
static constexpr double perThousand
Definition: G4SIunits.hh:333
T sqr(const T &x)
Definition: templates.hh:145
virtual void ComputeRightHandSide(const G4double y[], G4double dydx[])
double G4double
Definition: G4Types.hh:76
G4double GetPgrow() const
void PrintStatus(const G4double *StartArr, G4double xstart, const G4double *CurrentArr, G4double xcurrent, G4double requestStep, G4int subStepNo)
static G4GeometryTolerance * GetInstance()
G4double GetSafety() const
G4GLOB_DLL std::ostream G4cerr