Context Navigation

G4MagIntegratorDriver.cc @ 1342

Last change on this file since 1342 was 1340, checked in by garnier, 14 years ago
update ti head
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.57 2010/07/14 10:00:36 gcosmo Exp $
28	// GEANT4 tag $Name: field-V09-03-03 $
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
56	// divided by the order of Stepper
57	//
58	const 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	//
68	G4MagInt_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	//
110	G4MagInt_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
124	G4bool
125	G4MagInt_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
426	void
427	G4MagInt_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
458	void
459	G4MagInt_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
473	void
474	G4MagInt_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
520	void
521	G4MagInt_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()hstd::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	//
640	G4bool 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
659	G4bool 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
739	G4bool 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	//
759	G4double
760	G4MagInt_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()hstepCurrentstd::pow(errMaxNorm,GetPshrnk()) ;
771	} else if(errMaxNorm > 0.0 ) {
772	// Compute size of next Step for a successful step
773	hnew = GetSafety()hstepCurrentstd::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	//
789	G4double
790	G4MagInt_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()hstepCurrentstd::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()hstepCurrentstd::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
822	void 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
847	void 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
925	void 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
991	void 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
1039	void 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.

Context Navigation

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

Download in other formats: