Context Navigation

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

Last change on this file since 1315 was 1315, checked in by garnier, 14 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	//
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	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
420	void
421	G4MagInt_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
452	void
453	G4MagInt_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
467	void
468	G4MagInt_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
514	void
515	G4MagInt_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()hstd::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	//
633	G4bool 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
652	G4bool 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
732	G4bool 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	//
752	G4double
753	G4MagInt_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()hstepCurrentstd::pow(errMaxNorm,GetPshrnk()) ;
764	} else if(errMaxNorm > 0.0 ) {
765	// Compute size of next Step for a successful step
766	hnew = GetSafety()hstepCurrentstd::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	//
782	G4double
783	G4MagInt_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()hstepCurrentstd::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()hstepCurrentstd::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
815	void 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
840	void 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
918	void 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
984	void 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
1032	void 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.

Download in other formats: