Context Navigation

G4MagIntegratorDriver.cc@ 1347

Last change on this file since 1347 was 1340, checked in by garnier, 15 years ago
update ti head
File size: 35.6 KB

Rev	Line
[831]	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	//
[1340]	27	// $Id: G4MagIntegratorDriver.cc,v 1.57 2010/07/14 10:00:36 gcosmo Exp $
	28	// GEANT4 tag $Name: field-V09-03-03 $
[831]	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
[1231]	55	// The (default) maximum number of steps is Base
	56	// divided by the order of Stepper
[831]	57	//
	58	const G4int G4MagInt_Driver::fMaxStepBase = 250; // Was 5000
	59
	60	#ifndef G4NO_FIELD_STATISTICS
	61	#define G4FLD_STATS 1
	62	#endif
	63
[1231]	64	// ---------------------------------------------------------
	65
[831]	66	// Constructor
	67	//
	68	G4MagInt_Driver::G4MagInt_Driver( G4double hminimum,
[1231]	69	G4MagIntegratorStepper *pItsStepper,
	70	G4int numComponents,
	71	G4int statisticsVerbose)
[831]	72	: fSmallestFraction( 1.0e-12 ),
	73	fNoIntegrationVariables(numComponents),
	74	fMinNoVars(12),
	75	fNoVars( std::max( fNoIntegrationVariables, fMinNoVars )),
	76	fVerboseLevel(0),
[1231]	77	fNoTotalSteps(0), fNoBadSteps(0), fNoSmallSteps(0),
	78	fNoInitialSmallSteps(0), fDyerr_max(0.0), fDyerr_mx2(0.0),
[831]	79	fDyerrPos_smTot(0.0), fDyerrPos_lgTot(0.0), fDyerrVel_lgTot(0.0),
[1340]	80	fSumH_sm(0.0), fSumH_lg(0.0),
[831]	81	fStatisticsVerboseLevel(statisticsVerbose)
	82	{
[1231]	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
[831]	85
[1231]	86	RenewStepperAndAdjust( pItsStepper );
	87	fMinimumStep= hminimum;
	88	fMaxNoSteps = fMaxStepBase / pIntStepper->IntegratorOrder();
[831]	89	#ifdef G4DEBUG_FIELD
[1231]	90	fVerboseLevel=2;
[831]	91	#endif
	92
[1231]	93	if( (fVerboseLevel > 0) \|\| (fStatisticsVerboseLevel > 1) )
	94	{
	95	G4cout << "MagIntDriver version: Accur-Adv: "
	96	<< "invE_nS, QuickAdv-2sqrt with Statistics "
[831]	97	#ifdef G4FLD_STATS
[1231]	98	<< " enabled "
[831]	99	#else
[1231]	100	<< " disabled "
[831]	101	#endif
[1231]	102	<< G4endl;
	103	}
[831]	104	}
	105
[1231]	106	// ---------------------------------------------------------
	107
[831]	108	// Destructor
	109	//
	110	G4MagInt_Driver::~G4MagInt_Driver()
	111	{
[1231]	112	if( fStatisticsVerboseLevel > 1 )
	113	{
	114	PrintStatisticsReport();
	115	}
[831]	116	}
	117
[1231]	118	// To add much printing for debugging purposes, uncomment the following
	119	// and set verbose level to 1 or higher value !
[831]	120	// #define G4DEBUG_FIELD 1
	121
[1231]	122	// ---------------------------------------------------------
	123
[831]	124	G4bool
	125	G4MagInt_Driver::AccurateAdvance(G4FieldTrack& y_current,
[1231]	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 .
[831]	135
[1231]	136	G4int nstp, i, no_warnings=0;
	137	G4double x, hnext, hdid, h;
[831]	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
[1231]	155	const G4int nvar= fNoVars;
[831]	156
	157	G4FieldTrack yStartFT(y_current);
	158
[1231]	159	// Ensure that hstep > 0
	160	//
[831]	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
[1231]	188	if ( (hinitial > 0.0) && (hinitial < hstep)
	189	&& (hinitial > perMillion * hstep) )
	190	{
[831]	191	h = hinitial;
[1231]	192	}
	193	else // Initial Step size "h" defaults to the full interval
	194	{
[831]	195	h = hstep;
	196	}
	197
	198	x = x1;
	199
[1231]	200	for (i=0;i<nvar;i++) { y[i] = ystart[i]; }
[831]	201
[1231]	202	G4bool lastStep= false;
[831]	203	nstp=1;
	204
[1231]	205	do
	206	{
	207	G4ThreeVector StartPos( y[0], y[1], y[2] );
[831]	208
[1231]	209	#ifdef G4DEBUG_FIELD
[1340]	210	G4double xSubStepStart= x;
[1231]	211	for (i=0;i<nvar;i++) { ySubStepStart[i] = y[i]; }
	212	yFldTrkStart.LoadFromArray(y, fNoIntegrationVariables);
	213	yFldTrkStart.SetCurveLength(x);
	214	#endif
[831]	215
[1231]	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
[1340]	230	if (dbg>2) {
	231	PrintStatus( ySubStepStart, xSubStepStart, y, x, h, nstp); // Only
[1231]	232	}
[831]	233	#endif
[1231]	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 );
[831]	242
[1231]	243	QuickAdvance( yFldTrk, dydx, h, dchord_step, dyerr_len );
	244	//-----------------------------------------------------
[831]	245
[1231]	246	yFldTrk.DumpToArray(y);
[831]	247
	248	#ifdef G4FLD_STATS
[1231]	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++; }
[831]	254	#endif
	255	#ifdef G4DEBUG_FIELD
[1231]	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;
[831]	276
[1231]	277	// Compute suggested new step
	278	hnext= ComputeNewStepSize( dyerr/eps, h);
[831]	279
[1231]	280	// .. hnext= ComputeNewStepSize_WithinLimits( dyerr/eps, h);
	281	lastStepSucceeded= (dyerr<= eps);
	282	}
[831]	283
[1231]	284	if (lastStepSucceeded) { noFullIntegr++; }
	285	else { noSmallIntegr++; }
	286
	287	G4ThreeVector EndPos( y[0], y[1], y[2] );
	288
[831]	289	#ifdef G4DEBUG_FIELD
[1231]	290	if( (dbg>0) && (dbg<=2) && (nstp>nStpPr))
	291	{
	292	if( nstp==nStpPr ) { G4cout << "*** Many steps **" << G4endl; }
[1340]	293	G4cout << "MagIntDrv: " ;
[1231]	294	G4cout << "hdid=" << std::setw(12) << hdid << " "
[1340]	295	<< "hnext=" << std::setw(12) << hnext << " "
	296	<< "hstep=" << std::setw(12) << hstep << " (requested) "
	297	<< G4endl;
[1231]	298	PrintStatus( ystart, x1, y, x, h, (nstp==nStpPr) ? -nstp: nstp);
	299	}
[831]	300	#endif
	301
[1231]	302	// Check the endpoint
	303	G4double endPointDist= (EndPos-StartPos).mag();
	304	if ( endPointDist >= hdid*(1.+perMillion) )
	305	{
	306	fNoBadSteps++;
[831]	307
[1231]	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)
[831]	314	{
[1231]	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	}
[831]	321	#endif
[1231]	322	no_warnings++;
	323	}
	324	}
	325	else
	326	{
	327	noGoodSteps ++;
	328	}
	329	// #endif
[831]	330
[1231]	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	{
[831]	342	#ifdef G4DEBUG_FIELD
[1231]	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	}
[831]	354	#endif
[1231]	355	// Make sure that the next step is at least Hmin.
	356	h = Hmin();
	357	}
	358	else
	359	{
	360	h = hnext;
	361	}
[831]	362
[1231]	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
[831]	368
[1231]	369	if ( h == 0.0 )
	370	{
	371	// Cannot progress - accept this as last step - by default
	372	lastStep = true;
	373	#ifdef G4DEBUG_FIELD
[1340]	374	if (dbg>2)
[1231]	375	{
[1340]	376	int prec= G4cout.precision(12);
[1231]	377	G4cout << "Warning: G4MagIntegratorDriver::AccurateAdvance"
	378	<< G4endl
	379	<< " Integration step 'h' became "
[1340]	380	<< h << " due to roundoff. " << G4endl
	381	<< " Calculated as difference of x2= "<< x2 << " and x=" << x
[1231]	382	<< " Forcing termination of advance." << G4endl;
[1340]	383	G4cout.precision(prec);
[1231]	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
[831]	392	succeeded= (x>=x2); // If it was a "forced" last step
	393
[1231]	394	for (i=0;i<nvar;i++) { yEnd[i] = y[i]; }
[831]	395
	396	// Put back the values.
	397	y_current.LoadFromArray( yEnd, fNoIntegrationVariables );
	398	y_current.SetCurveLength( x );
	399
[1231]	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	}
[831]	410	#endif
	411	}
	412
	413	#ifdef G4DEBUG_FIELD
[1231]	414	if( dbg && no_warnings )
	415	{
	416	G4cerr << "G4MagIntegratorDriver exit status: no-steps " << nstp <<G4endl;
	417	PrintStatus( yEnd, x1, y, x, hstep, nstp);
[831]	418	}
	419	#endif
	420
	421	return succeeded;
	422	} // end of AccurateAdvance ...........................
	423
[1231]	424	// ---------------------------------------------------------
	425
[831]	426	void
	427	G4MagInt_Driver::WarnSmallStepSize( G4double hnext, G4double hstep,
[1231]	428	G4double h, G4double xDone,
	429	G4int nstp)
[831]	430	{
	431	static G4int noWarningsIssued =0;
	432	const G4int maxNoWarnings = 10; // Number of verbose warnings
[1231]	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;
[831]	447	G4cerr << " this sub-step " << h
[1231]	448	<< " req_tot_len " << hstep
	449	<< " done " << xDone
	450	<< " min " << Hmin()
	451	<< G4endl;
[831]	452	}
	453	noWarningsIssued++;
	454	}
	455
[1231]	456	// ---------------------------------------------------------
	457
[831]	458	void
	459	G4MagInt_Driver::WarnTooManyStep( G4double x1start,
[1231]	460	G4double x2end,
	461	G4double xCurrent)
[831]	462	{
[1231]	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;
[831]	469	}
	470
[1231]	471	// ---------------------------------------------------------
	472
[831]	473	void
	474	G4MagInt_Driver::WarnEndPointTooFar (G4double endPointDist,
[1231]	475	G4double h ,
	476	G4double eps,
	477	G4int dbg)
[831]	478	{
[1231]	479	static G4double maxRelError=0.0, maxRelError_last_printed=0.0;
	480	G4bool isNewMax, prNewMax;
[831]	481
[1231]	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; }
[831]	486
[1231]	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;
[831]	497
[1231]	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	}
[831]	516	}
[1231]	517
[831]	518	// ---------------------------------------------------------
	519
	520	void
	521	G4MagInt_Driver::OneGoodStep( G4double y[], // InOut
[1231]	522	const G4double dydx[],
	523	G4double& x, // InOut
	524	G4double htry,
	525	G4double eps_rel_max,
	526	G4double& hdid, // Out
	527	G4double& hnext ) // Out
[831]	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	{
[1231]	543	G4double errmax_sq;
	544	G4double h, htemp, xnew ;
[831]	545
[1231]	546	G4double yerr[G4FieldTrack::ncompSVEC], ytemp[G4FieldTrack::ncompSVEC];
[831]	547
[1231]	548	h = htry ; // Set stepsize to the initial trial value
[831]	549
[1231]	550	G4double inv_eps_vel_sq = 1.0 / (eps_rel_max*eps_rel_max);
[831]	551
[1231]	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
[831]	555
[1231]	556	G4int iter;
[831]	557
[1231]	558	static G4int tot_no_trials=0;
	559	const G4int max_trials=100;
[831]	560
[1231]	561	G4ThreeVector Spin(y[9],y[10],y[11]);
	562	G4bool hasSpin= (Spin.mag2() > 0.0);
[831]	563
[1231]	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);
[831]	571
[1231]	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
[831]	576
[1231]	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
[831]	582
[1231]	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;
[1340]	589	errmax_sq = std::max( errmax_sq, errspin_sq );
	590	}
[831]	591
[1231]	592	if ( errmax_sq <= 1.0 ) { break; } // Step succeeded.
[831]	593
[1231]	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
[831]	613	#ifdef G4FLD_STATS
[1231]	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;
[831]	618	#endif
	619
[1231]	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);
[831]	630
[1231]	631	for(G4int k=0;k<fNoIntegrationVariables;k++) { y[k] = ytemp[k]; }
[831]	632
[1231]	633	return;
[831]	634	} // end of OneGoodStep .............................
	635
	636	//----------------------------------------------------------------------
[1231]	637
[831]	638	// QuickAdvance just tries one Step - it does not ensure accuracy
	639	//
	640	G4bool G4MagInt_Driver::QuickAdvance(
[1231]	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 )
[831]	647	{
[1231]	648	G4Exception("G4MagInt_Driver::QuickAdvance()", "NotImplemented",
	649	FatalException, "Not yet implemented.");
[831]	650
[1231]	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;
[831]	655	}
	656
[1231]	657	//----------------------------------------------------------------------
	658
[831]	659	G4bool G4MagInt_Driver::QuickAdvance(
[1231]	660	G4FieldTrack& y_posvel, // INOUT
	661	const G4double dydx[],
	662	G4double hstep, // In
	663	G4double& dchord_step,
	664	G4double& dyerr )
[831]	665	{
[1231]	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;
[831]	671
[1231]	672	static G4int no_call=0;
	673	no_call ++;
[831]	674
[1231]	675	// Move data into array
	676	y_posvel.DumpToArray( yarrin ); // yarrin <== y_posvel
	677	s_start = y_posvel.GetCurveLength();
[831]	678
[1231]	679	// Do an Integration Step
	680	pIntStepper-> Stepper(yarrin, dydx, hstep, yarrout, yerr_vec) ;
	681	// *******
[831]	682
[1231]	683	// Estimate curve-chord distance
	684	dchord_step= pIntStepper-> DistChord();
	685	// *********
[831]	686
[1231]	687	// Put back the values. yarrout ==> y_posvel
	688	y_posvel.LoadFromArray( yarrout, fNoIntegrationVariables );
	689	y_posvel.SetCurveLength( s_start + hstep );
[831]	690
	691	#ifdef G4DEBUG_FIELD
[1231]	692	if(fVerboseLevel>2)
	693	{
	694	G4cout << "G4MagIntDrv: Quick Advance" << G4endl;
	695	PrintStatus( yarrin, s_start, yarrout, s_start+hstep, hstep, 1);
	696	}
[831]	697	#endif
	698
[1231]	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;
[831]	706
[1231]	707	// Calculate also the change in the momentum squared also ???
	708	// G4double veloc_square = y_posvel.GetVelocity().mag2();
	709	// ...
[831]	710
	711	#ifdef RETURN_A_NEW_STEP_LENGTH
[1231]	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 ;
[831]	715
[1231]	716	// Look at the velocity deviation ?
	717	// sqr(yerr_vec[3])+sqr(yerr_vec[4])+sqr(yerr_vec[5]));
[831]	718
[1231]	719	// Set suggested new step
	720	hstep= ComputeNewStepSize( dyerr_len, hstep);
[831]	721	#endif
	722
[1231]	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	}
[831]	732
[1231]	733	return true;
[831]	734	}
	735
[1231]	736	// --------------------------------------------------------------------------
	737
[831]	738	#ifdef QUICK_ADV_ARRAY_IN_AND_OUT
	739	G4bool G4MagInt_Driver::QuickAdvance(
[1231]	740	G4double yarrin[], // In
	741	const G4double dydx[],
	742	G4double hstep, // In
	743	G4double yarrout[],
	744	G4double& dchord_step,
	745	G4double& dyerr ) // In length
[831]	746	{
[1231]	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];
[831]	751	}
	752	#endif
	753
	754	// --------------------------------------------------------------------------
[1231]	755
[831]	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)
[1231]	762	G4double hstepCurrent) // current step size
[831]	763	{
	764	G4double hnew;
	765
	766	// Compute size of next Step for a failed step
[1231]	767	if(errMaxNorm > 1.0 )
	768	{
[831]	769	// Step failed; compute the size of retrial Step.
	770	hnew = GetSafety()hstepCurrentstd::pow(errMaxNorm,GetPshrnk()) ;
[1231]	771	} else if(errMaxNorm > 0.0 ) {
[831]	772	// Compute size of next Step for a successful step
	773	hnew = GetSafety()hstepCurrentstd::pow(errMaxNorm,GetPgrow()) ;
[1231]	774	} else {
[831]	775	// if error estimate is zero (possible) or negative (dubious)
	776	hnew = max_stepping_increase * hstepCurrent;
	777	}
	778
	779	return hnew;
	780	}
	781
[1231]	782	// ---------------------------------------------------------------------------
	783
	784	// This method computes new step sizes limiting changes within certain factors
[831]	785	//
[1231]	786	// It shares its logic with AccurateAdvance.
	787	// They are kept separate currently for optimisation.
	788	//
[831]	789	G4double
	790	G4MagInt_Driver::ComputeNewStepSize_WithinLimits(
	791	G4double errMaxNorm, // max error (normalised)
[1231]	792	G4double hstepCurrent) // current step size
[831]	793	{
	794	G4double hnew;
	795
	796	// Compute size of next Step for a failed step
[1231]	797	if (errMaxNorm > 1.0 )
	798	{
[831]	799	// Step failed; compute the size of retrial Step.
	800	hnew = GetSafety()hstepCurrentstd::pow(errMaxNorm,GetPshrnk()) ;
	801
[1231]	802	if (hnew < max_stepping_decrease*hstepCurrent)
	803	{
	804	hnew = max_stepping_decrease*hstepCurrent ;
[831]	805	// reduce stepsize, but no more
	806	// than this factor (value= 1/10)
[1231]	807	}
	808	}
	809	else
	810	{
[831]	811	// Compute size of next Step for a successful step
[1231]	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; }
[831]	816	}
	817	return hnew;
	818	}
	819
[1231]	820	// ---------------------------------------------------------------------------
[831]	821
	822	void G4MagInt_Driver::PrintStatus( const G4double* StartArr,
[1231]	823	G4double xstart,
	824	const G4double* CurrentArr,
	825	G4double xcurrent,
	826	G4double requestStep,
	827	G4int subStepNo)
[831]	828	// Potentially add as arguments:
	829	// <dydx> - as Initial Force
	830	// stepTaken(hdid) - last step taken
	831	// nextStep (hnext) - proposal for size
	832	{
[1231]	833	G4FieldTrack StartFT(G4ThreeVector(0,0,0),
	834	G4ThreeVector(0,0,0), 0., 0., 0., 0. );
[831]	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
[1231]	845	// ---------------------------------------------------------------------------
[831]	846
	847	void G4MagInt_Driver::PrintStatus(
	848	const G4FieldTrack& StartFT,
[1231]	849	const G4FieldTrack& CurrentFT,
[831]	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
[1231]	858	const G4ThreeVector StartPosition= StartFT.GetPosition();
	859	const G4ThreeVector StartUnitVelocity= StartFT.GetMomentumDir();
	860	const G4ThreeVector CurrentPosition= CurrentFT.GetPosition();
	861	const G4ThreeVector CurrentUnitVelocity= CurrentFT.GetMomentumDir();
[831]	862
	863	G4double DotStartCurrentVeloc= StartUnitVelocity.dot(CurrentUnitVelocity);
	864
[1231]	865	G4double step_len= CurrentFT.GetCurveLength() - StartFT.GetCurveLength();
[831]	866	G4double subStepSize = step_len;
	867
	868	if( (subStepNo <= 1) \|\| (verboseLevel > 3) )
	869	{
	870	subStepNo = - subStepNo; // To allow printing banner
	871
[1231]	872	G4cout << std::setw( 6) << " " << std::setw( 25)
	873	<< " G4MagInt_Driver: Current Position and Direction" << " "
	874	<< G4endl;
[831]	875	G4cout << std::setw( 5) << "Step#" << " "
[1231]	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;
[831]	891	}
	892
[1231]	893	if( (subStepNo <= 0) )
	894	{
	895	PrintStat_Aux( StartFT, requestStep, 0.,
	896	0, 0.0, 1.0);
	897	//*************
[831]	898	}
	899
	900	if( verboseLevel <= 3 )
	901	{
[1231]	902	G4cout.precision(noPrecision);
	903	PrintStat_Aux( CurrentFT, requestStep, step_len,
	904	subStepNo, subStepSize, DotStartCurrentVeloc );
	905	//*************
[831]	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;
[1231]	917	// G4cout << "Chord length = " << (CurrentPosition-StartPosition).mag()
	918	// << G4endl << G4endl;
[831]	919	}
	920	G4cout.precision(oldPrec);
	921	}
	922
[1231]	923	// ---------------------------------------------------------------------------
	924
[831]	925	void G4MagInt_Driver::PrintStat_Aux(
	926	const G4FieldTrack& aFieldTrack,
	927	G4double requestStep,
[1231]	928	G4double step_len,
[831]	929	G4int subStepNo,
[1231]	930	G4double subStepSize,
	931	G4double dotVeloc_StartCurr)
[831]	932	{
	933	const G4ThreeVector Position= aFieldTrack.GetPosition();
	934	const G4ThreeVector UnitVelocity= aFieldTrack.GetMomentumDir();
	935
	936	if( subStepNo >= 0)
[1231]	937	{
[831]	938	G4cout << std::setw( 5) << subStepNo << " ";
[1231]	939	}
[831]	940	else
[1231]	941	{
[831]	942	G4cout << std::setw( 5) << "Start" << " ";
[1231]	943	}
[831]	944	G4double curveLen= aFieldTrack.GetCurveLength();
	945	G4cout << std::setw( 7) << curveLen;
	946	G4cout << std::setw( 9) << Position.x() << " "
[1231]	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() << " ";
[831]	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;
[1231]	962	static G4int oldSubStepNo= -1;
[831]	963
	964	G4double subStep_len=0.0;
	965	if( curveLen > oldCurveLength )
[1231]	966	{
	967	subStep_len= curveLen - oldCurveLength;
	968	}
[831]	969	else if (subStepNo == oldSubStepNo)
[1231]	970	{
	971	subStep_len= oldSubStepLength;
	972	}
[831]	973	oldCurveLength= curveLen;
	974	oldSubStepLength= subStep_len;
	975
	976	G4cout << std::setw(12) << subStep_len << " ";
	977	G4cout << std::setw(12) << subStepSize << " ";
[1231]	978	if( requestStep != -1.0 )
	979	{
	980	G4cout << std::setw( 9) << requestStep << " ";
	981	}
[831]	982	else
[1231]	983	{
	984	G4cout << std::setw( 9) << " InitialStep " << " ";
	985	}
[831]	986	G4cout << G4endl;
	987	}
	988
[1231]	989	// ---------------------------------------------------------------------------
	990
[831]	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: "
[1231]	998	<< " Total= " << fNoTotalSteps
[831]	999	<< " Bad= " << fNoBadSteps
	1000	<< " Small= " << fNoSmallSteps
[1231]	1001	<< " Non-initial small= " << (fNoSmallSteps-fNoInitialSmallSteps)
	1002	<< G4endl;
[831]	1003
[1231]	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;
[831]	1013
	1014	#if 0
	1015	G4int noPrecSmall=4;
	1016	// Single line precis of statistics ... optional
	1017	G4cout.precision(noPrecSmall);
	1018	G4cout << "MIDnums: " << fMinimumStep
[1231]	1019	<< " " << fNoTotalSteps
	1020	<< " " << fNoSmallSteps
[831]	1021	<< " " << fNoSmallSteps-fNoInitialSmallSteps
[1231]	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
[831]	1033
	1034	G4cout.precision(oldPrec);
	1035	}
	1036
[1231]	1037	// ---------------------------------------------------------------------------
	1038
[831]	1039	void G4MagInt_Driver::SetSmallestFraction(G4double newFraction)
	1040	{
[1231]	1041	if( (newFraction > 1.e-16) && (newFraction < 1e-8) )
	1042	{
	1043	fSmallestFraction= newFraction;
[831]	1044	}
[1231]	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	}
[831]	1051	}

Note: See TracBrowser for help on using the repository browser.

Context Navigation

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

Download in other formats: