source: trunk/source/geometry/magneticfield/src/G4MagIntegratorDriver.cc@ 1199

Last change on this file since 1199 was 921, checked in by garnier, 17 years ago

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