99{
100
101
102
103
104
105
106 G4int nstp, i, no_warnings = 0;
108
109#ifdef G4DEBUG_FIELD
110 static G4int dbg = 1;
113#endif
114
118 G4bool succeeded =
true, lastStepSucceeded;
119
121
122 G4int noFullIntegr = 0, noSmallIntegr = 0;
124 const G4int nvar = fNoVars;
125
127
128
129
130 if( hstep <= 0.0 )
131 {
132 if( hstep == 0.0 )
133 {
134 std::ostringstream message;
135 message << "Proposed step is zero; hstep = " << hstep << " !";
136 G4Exception(
"G4OldMagIntDriver::AccurateAdvance()",
138 return succeeded;
139 }
140 else
141 {
142 std::ostringstream message;
143 message <<
"Invalid run condition." <<
G4endl
144 <<
"Proposed step is negative; hstep = " << hstep <<
"." <<
G4endl
145 << "Requested step cannot be negative! Aborting event.";
146 G4Exception(
"G4OldMagIntDriver::AccurateAdvance()",
148 return false;
149 }
150 }
151
153
155 x1= startCurveLength;
156 x2= x1 + hstep;
157
158 if ( (hinitial > 0.0) && (hinitial < hstep)
159 && (hinitial > perMillion * hstep) )
160 {
161 h = hinitial;
162 }
163 else
164 {
165 h = hstep;
166 }
167
168 x = x1;
169
170 for ( i=0; i<nvar; ++i) { y[i] = ystart[i]; }
171
172#ifdef G4DEBUG_FIELD
173
174 G4cout <<
"IDriver::AccurAdv called. "
175 <<
" Input: hstep = " << hstep <<
" hinitial= " << hinitial <<
" , current: h = " << h <<
G4endl;
176#endif
177
179 nstp = 1;
180
181 do
182 {
184
185#ifdef G4DEBUG_FIELD
187 for (i=0; i<nvar; ++i) { ySubStepStart[i] = y[i]; }
188 yFldTrkStart.LoadFromArray(y, fNoIntegrationVariables);
189 yFldTrkStart.SetCurveLength(x);
190 if(dbg)
192#endif
193
195 ++fNoTotalSteps;
196
197
198
199 if( h > fMinimumStep )
200 {
202
203 lastStepSucceeded = (hdid == h);
204#ifdef G4DEBUG_FIELD
205 if (dbg)
206 {
207 G4cout <<
"IntegrationDriver -- after OneGoodStep / requesting step = " << h <<
G4endl;
208
210 }
211#endif
212 }
213 else
214 {
217 G4double dchord_step, dyerr, dyerr_len;
218 yFldTrk.LoadFromArray(y, fNoIntegrationVariables);
219 yFldTrk.SetCurveLength( x );
220
221 QuickAdvance( yFldTrk, dydx, h, dchord_step, dyerr_len );
222
223
224 yFldTrk.DumpToArray(y);
225
226#ifdef G4FLD_STATS
227 ++fNoSmallSteps;
228 if ( dyerr_len > fDyerr_max ) { fDyerr_max = dyerr_len; }
229 fDyerrPos_smTot += dyerr_len;
230 fSumH_sm += h;
231 if (nstp<=1) { ++fNoInitialSmallSteps; }
232#endif
233#ifdef G4DEBUG_FIELD
234 if (dbg>1)
235 {
236 if(fNoSmallSteps<2) {
PrintStatus(ySubStepStart, x1, y, x, h, -nstp); }
237 G4cout <<
"Another sub-min step, no " << fNoSmallSteps
238 <<
" of " << fNoTotalSteps <<
" this time " << nstp <<
G4endl;
240 G4cout <<
" dyerr= " << dyerr_len <<
" relative = " << dyerr_len / h
241 << " epsilon= " << eps << " hstep= " << hstep
242 <<
" h= " << h <<
" hmin= " << fMinimumStep <<
G4endl;
243 }
244#endif
245 if( h == 0.0 )
246 {
247 G4Exception(
"G4OldMagIntDriver::AccurateAdvance()",
249 "Integration Step became Zero!");
250 }
251 dyerr = dyerr_len / h;
252 hdid = h;
253 x += hdid;
254
255
257
258
259 lastStepSucceeded = (dyerr<= eps);
260 }
261
262 if (lastStepSucceeded) { ++noFullIntegr; }
263 else { ++noSmallIntegr; }
264
266
267#if (G4DEBUG_FIELD>1)
268 static G4int nStpPr = 50;
269 if( (dbg>0) && (dbg<=2) && (nstp>nStpPr))
270 {
271 if( nstp==nStpPr ) {
G4cout <<
"***** Many steps ****" <<
G4endl; }
273 G4cout <<
"hdid=" << std::setw(12) << hdid <<
" "
274 << "hnext=" << std::setw(12) << hnext << " "
275 << "hstep=" << std::setw(12) << hstep << " (requested) "
277 PrintStatus( ystart, x1, y, x, h, (nstp==nStpPr) ? -nstp: nstp);
278 }
279#endif
280
281
282 G4double endPointDist= (EndPos-StartPos).mag();
283 if ( endPointDist >= hdid*(1.+perMillion) )
284 {
285 ++fNoBadSteps;
286
287
288
289 if ( endPointDist >= hdid*(1.+perThousand) )
290 {
291#ifdef G4DEBUG_FIELD
292 if (dbg)
293 {
295 G4cerr <<
" Total steps: bad " << fNoBadSteps
296 << " good " << noGoodSteps << " current h= " << hdid
298 PrintStatus( ystart, x1, y, x, hstep, no_warnings?nstp:-nstp);
299 }
300#endif
301 ++no_warnings;
302 }
303 }
304 else
305 {
306 ++noGoodSteps;
307 }
308
309
310
311 if( (h < eps * hstep) || (h < fSmallestFraction * startCurveLength) )
312 {
313
314 lastStep = true;
315 }
316 else
317 {
318
319 if(std::fabs(hnext) <=
Hmin())
320 {
321#ifdef G4DEBUG_FIELD
322
323 if( (x < x2 * (1-eps) ) &&
324 (std::fabs(hstep) >
Hmin()) )
325 {
326 if(dbg>0)
327 {
329 PrintStatus( ystart, x1, y, x, hstep, no_warnings?nstp:-nstp);
330 }
331 ++no_warnings;
332 }
333#endif
334
336 }
337 else
338 {
339 h = hnext;
340 }
341
342
343 if ( x+h > x2 )
344 {
345 h = x2 - x ;
346 }
347
348 if ( h == 0.0 )
349 {
350
351 lastStep = true;
352#ifdef G4DEBUG_FIELD
353 if (dbg>2)
354 {
355 int prec=
G4cout.precision(12);
356 G4cout <<
"Warning: G4MagIntegratorDriver::AccurateAdvance"
358 << " Integration step 'h' became "
359 << h <<
" due to roundoff. " <<
G4endl
360 << " Calculated as difference of x2= "<< x2 << " and x=" << x
361 <<
" Forcing termination of advance." <<
G4endl;
363 }
364#endif
365 }
366 }
367 } while ( ((nstp++)<=fMaxNoSteps) && (x < x2) && (!lastStep) );
368
369
370
371
372
373 succeeded = (x>=x2);
374
375 for (i=0; i<nvar; ++i) { yEnd[i] = y[i]; }
376
377
380
381 if(nstp > fMaxNoSteps)
382 {
383 ++no_warnings;
384 succeeded = false;
385#ifdef G4DEBUG_FIELD
386 if (dbg)
387 {
390 }
391#endif
392 }
393
394#ifdef G4DEBUG_FIELD
395 if( dbg && no_warnings )
396 {
397 G4cerr <<
"G4MagIntegratorDriver exit status: no-steps " << nstp <<
G4endl;
399 }
400#endif
401
402 return succeeded;
403}
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
CLHEP::Hep3Vector G4ThreeVector
G4GLOB_DLL std::ostream G4cerr
static void PrintStatus(const G4double *StartArr, G4double xstart, const G4double *CurrentArr, G4double xcurrent, G4double requestStep, unsigned int subStepNo, unsigned int noIntegrationVariables)
G4double GetCurveLength() const
void SetCurveLength(G4double nCurve_s)
void DumpToArray(G4double valArr[ncompSVEC]) const
void LoadFromArray(const G4double valArr[ncompSVEC], G4int noVarsIntegrated)
void RightHandSide(const G4double y[], G4double dydx[]) const
virtual G4bool QuickAdvance(G4FieldTrack &y_val, const G4double dydx[], G4double hstep, G4double &dchord_step, G4double &dyerr) override
void PrintStatus(const G4double *StartArr, G4double xstart, const G4double *CurrentArr, G4double xcurrent, G4double requestStep, G4int subStepNo)
void OneGoodStep(G4double ystart[], const G4double dydx[], G4double &x, G4double htry, G4double eps, G4double &hdid, G4double &hnext)
void WarnSmallStepSize(G4double hnext, G4double hstep, G4double h, G4double xDone, G4int noSteps)
virtual G4double ComputeNewStepSize(G4double errMaxNorm, G4double hstepCurrent) override
void WarnEndPointTooFar(G4double endPointDist, G4double hStepSize, G4double epsilonRelative, G4int debugFlag)
void WarnTooManyStep(G4double x1start, G4double x2end, G4double xCurrent)