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

Last change on this file since 1292 was 1231, checked in by garnier, 16 years ago

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