99{
100
101
102
103
104
105
108
109#ifdef G4DEBUG_FIELD
110 G4int no_warnings = 0;
111 static G4int dbg = 1;
114#endif
115
120
122
123 const G4int nvar = fNoVars;
124
126
127
128
129 if( hstep <= 0.0 )
130 {
131 if( hstep == 0.0 )
132 {
133 std::ostringstream message;
134 message << "Proposed step is zero; hstep = " << hstep << " !";
135 G4Exception(
"G4OldMagIntDriver::AccurateAdvance()",
137 return succeeded;
138 }
139 else
140 {
141 std::ostringstream message;
142 message <<
"Invalid run condition." <<
G4endl
143 <<
"Proposed step is negative; hstep = " << hstep <<
"." <<
G4endl
144 << "Requested step cannot be negative! Aborting event.";
145 G4Exception(
"G4OldMagIntDriver::AccurateAdvance()",
147 return false;
148 }
149 }
150
152
154 x1= startCurveLength;
155 x2= x1 + hstep;
156
157 if ( (hinitial > 0.0) && (hinitial < hstep)
158 && (hinitial > perMillion * hstep) )
159 {
160 h = hinitial;
161 }
162 else
163 {
164 h = hstep;
165 }
166
167 x = x1;
168
169 for ( i=0; i<nvar; ++i) { y[i] = ystart[i]; }
170
171#ifdef G4DEBUG_FIELD
172
173 G4cout <<
"IDriver::AccurAdv called. "
174 <<
" Input: hstep = " << hstep <<
" hinitial= " << hinitial <<
" , current: h = " << h <<
G4endl;
175#endif
176
178 nstp = 1;
179
180 do
181 {
183
184#ifdef G4DEBUG_FIELD
186 for (i=0; i<nvar; ++i) { ySubStepStart[i] = y[i]; }
187 yFldTrkStart.LoadFromArray(y, fNoIntegrationVariables);
188 yFldTrkStart.SetCurveLength(x);
189 if(dbg)
191#endif
192
194 ++fNoTotalSteps;
195
196
197
198 if( h > fMinimumStep )
199 {
201
202
203#ifdef G4DEBUG_FIELD
204 if (dbg)
205 {
206 G4cout <<
"IntegrationDriver -- after OneGoodStep / requesting step = " << h <<
G4endl;
207
209 }
210#endif
211 }
212 else
213 {
216 G4double dchord_step, dyerr, dyerr_len;
217 yFldTrk.LoadFromArray(y, fNoIntegrationVariables);
218 yFldTrk.SetCurveLength( x );
219
220 QuickAdvance( yFldTrk, dydx, h, dchord_step, dyerr_len );
221
222
223 yFldTrk.DumpToArray(y);
224
225#ifdef G4FLD_STATS
226 ++fNoSmallSteps;
227 if ( dyerr_len > fDyerr_max ) { fDyerr_max = dyerr_len; }
228 fDyerrPos_smTot += dyerr_len;
229 fSumH_sm += h;
230 if (nstp<=1) { ++fNoInitialSmallSteps; }
231#endif
232#ifdef G4DEBUG_FIELD
233 if (dbg>1)
234 {
235 if(fNoSmallSteps<2) {
PrintStatus(ySubStepStart, x1, y, x, h, -nstp); }
236 G4cout <<
"Another sub-min step, no " << fNoSmallSteps
237 <<
" of " << fNoTotalSteps <<
" this time " << nstp <<
G4endl;
239 G4cout <<
" dyerr= " << dyerr_len <<
" relative = " << dyerr_len / h
240 << " epsilon= " << eps << " hstep= " << hstep
241 <<
" h= " << h <<
" hmin= " << fMinimumStep <<
G4endl;
242 }
243#endif
244 if( h == 0.0 )
245 {
246 G4Exception(
"G4OldMagIntDriver::AccurateAdvance()",
248 "Integration Step became Zero!");
249 }
250 dyerr = dyerr_len / h;
251 hdid = h;
252 x += hdid;
253
254
256
257
258 }
259
261
262#if (G4DEBUG_FIELD>1)
263 static G4int nStpPr = 50;
264 if( (dbg>0) && (dbg<=2) && (nstp>nStpPr))
265 {
266 if( nstp==nStpPr ) {
G4cout <<
"***** Many steps ****" <<
G4endl; }
268 G4cout <<
"hdid=" << std::setw(12) << hdid <<
" "
269 << "hnext=" << std::setw(12) << hnext << " "
270 << "hstep=" << std::setw(12) << hstep << " (requested) "
272 PrintStatus( ystart, x1, y, x, h, (nstp==nStpPr) ? -nstp: nstp);
273 }
274#endif
275
276
277 G4double endPointDist= (EndPos-StartPos).mag();
278 if ( endPointDist >= hdid*(1.+perMillion) )
279 {
280 ++fNoBadSteps;
281
282
283
284 if ( endPointDist >= hdid*(1.+perThousand) )
285 {
286#ifdef G4DEBUG_FIELD
287 if (dbg)
288 {
290 G4cerr <<
" Total steps: bad " << fNoBadSteps
291 <<
" current h= " << hdid <<
G4endl;
292 PrintStatus( ystart, x1, y, x, hstep, no_warnings?nstp:-nstp);
293 }
294 ++no_warnings;
295#endif
296 }
297 }
298
299
300 if( (h < eps * hstep) || (h < fSmallestFraction * startCurveLength) )
301 {
302
303 lastStep = true;
304 }
305 else
306 {
307
308 if(std::fabs(hnext) <=
Hmin())
309 {
310#ifdef G4DEBUG_FIELD
311
312 if( (x < x2 * (1-eps) ) &&
313 (std::fabs(hstep) >
Hmin()) )
314 {
315 if(dbg>0)
316 {
318 PrintStatus( ystart, x1, y, x, hstep, no_warnings?nstp:-nstp);
319 }
320 ++no_warnings;
321 }
322#endif
323
325 }
326 else
327 {
328 h = hnext;
329 }
330
331
332 if ( x+h > x2 )
333 {
334 h = x2 - x ;
335 }
336
337 if ( h == 0.0 )
338 {
339
340 lastStep = true;
341#ifdef G4DEBUG_FIELD
342 if (dbg>2)
343 {
344 int prec=
G4cout.precision(12);
345 G4cout <<
"Warning: G4MagIntegratorDriver::AccurateAdvance"
347 << " Integration step 'h' became "
348 << h <<
" due to roundoff. " <<
G4endl
349 << " Calculated as difference of x2= "<< x2 << " and x=" << x
350 <<
" Forcing termination of advance." <<
G4endl;
352 }
353#endif
354 }
355 }
356 } while ( ((nstp++)<=fMaxNoSteps) && (x < x2) && (!lastStep) );
357
358
359
360
361
362 succeeded = (x>=x2);
363
364 for (i=0; i<nvar; ++i) { yEnd[i] = y[i]; }
365
366
369
370 if(nstp > fMaxNoSteps)
371 {
372 succeeded = false;
373#ifdef G4DEBUG_FIELD
374 ++no_warnings;
375 if (dbg)
376 {
379 }
380#endif
381 }
382
383#ifdef G4DEBUG_FIELD
384 if( dbg && no_warnings )
385 {
386 G4cerr <<
"G4MagIntegratorDriver exit status: no-steps " << nstp <<
G4endl;
388 }
389#endif
390
391 return succeeded;
392}
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)