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

Last change on this file since 1326 was 1315, checked in by garnier, 15 years ago

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