Context Navigation

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

Last change on this file since 956 was 921, checked in by garnier, 15 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	//
52	const G4double G4MagInt_Driver::max_stepping_increase = 5.0;
53	const 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	//
57	const G4int G4MagInt_Driver::fMaxStepBase = 250; // Was 5000
58
59	#ifndef G4NO_FIELD_STATISTICS
60	#define G4FLD_STATS 1
61	#endif
62
63	// Constructor
64	//
65	G4MagInt_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	//
104	G4MagInt_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
117	G4bool
118	G4MagInt_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
387	void
388	G4MagInt_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
414	void
415	G4MagInt_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
427	void
428	G4MagInt_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
472	void
473	G4MagInt_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()hstd::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	//
583	G4bool 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
602	G4bool 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
677	G4bool 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
698	G4double
699	G4MagInt_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()hstepCurrentstd::pow(errMaxNorm,GetPshrnk()) ;
710	}else if(errMaxNorm > 0.0 ){
711	// Compute size of next Step for a successful step
712	hnew = GetSafety()hstepCurrentstd::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
727	G4double
728	G4MagInt_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()hstepCurrentstd::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()hstepCurrentstd::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
756	void 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
780	void 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
860	void 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
914	void 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
961	void 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.

Download in other formats: