Geant4-11
G4FSALIntegrationDriver.icc
Go to the documentation of this file.
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// G4FSALIntegrationDriver inline implementation
27//
28// Created: D.Sorokin, 2017
29// --------------------------------------------------------------------
30
31#include "G4FieldUtils.hh"
32
33template <class T>
34G4FSALIntegrationDriver<T>::
35G4FSALIntegrationDriver ( G4double hminimum, T* pStepper,
36 G4int numComponents, G4int statisticsVerbose )
37 : Base(pStepper),
38 fMinimumStep(hminimum),
39 fSmallestFraction(1e-12),
40 fVerboseLevel(statisticsVerbose),
41 fNoQuickAvanceCalls(0),
42 fNoAccurateAdvanceCalls(0),
43 fNoAccurateAdvanceBadSteps(0),
44 fNoAccurateAdvanceGoodSteps(0)
45{
46 if (numComponents != Base::GetStepper()->GetNumberOfVariables())
47 {
48 std::ostringstream message;
49 message << "Driver's number of integrated components "
50 << numComponents
51 << " != Stepper's number of components "
52 << pStepper->GetNumberOfVariables();
53 G4Exception("G4FSALIntegrationDriver","GeomField0002",
54 FatalException, message);
55 }
56}
57
58template <class T>
59G4FSALIntegrationDriver<T>::~G4FSALIntegrationDriver()
60{
61#ifdef G4VERBOSE
62 if( fVerboseLevel > 0 )
63 G4cout << "G4FSALIntegration Driver Stats: "
64 << "#QuickAdvance " << fNoQuickAvanceCalls
65 << " - #AccurateAdvance " << fNoAccurateAdvanceCalls << G4endl
66 << "#good steps " << fNoAccurateAdvanceGoodSteps << " "
67 << "#bad steps " << fNoAccurateAdvanceBadSteps << G4endl;
68#endif
69}
70
71// Runge-Kutta driver with adaptive stepsize control. Integrate starting
72// values at y_current over hstep x2 with accuracy eps.
73// On output ystart is replaced by values at the end of the integration
74// interval. RightHandSide is the right-hand side of ODE system.
75// The source is similar to odeint routine from NRC p.721-722 .
76//
77template <class T>
78G4bool G4FSALIntegrationDriver<T>::
79AccurateAdvance( G4FieldTrack& track, G4double hstep,
80 G4double eps, G4double hinitial )
81{
82 ++fNoAccurateAdvanceCalls;
83
84 if (hstep < GetMinimumStep())
85 {
86 G4double dchord_step = 0.0, dyerr = 0.0;
87 G4double dydx[G4FieldTrack::ncompSVEC];
88 Base::GetDerivatives(track, dydx);
89 return QuickAdvance(track, dydx, hstep, dchord_step, dyerr);
90 }
91
92 G4bool succeeded = false;
93
94 G4double hnext, hdid;
95
96 G4double y[G4FieldTrack::ncompSVEC], dydx[G4FieldTrack::ncompSVEC];
97
98 track.DumpToArray(y);
99
100 // hstep somtimes is too small. No need to add large curveLength.
101 G4double curveLength = 0.0;
102 G4double endCurveLength = hstep;
103
104
105 G4double h = hstep;
106 if (hinitial > CLHEP::perMillion * hstep && hinitial < hstep)
107 {
108 h = hinitial;
109 }
110
111 Base::GetStepper()->RightHandSide(y, dydx);
112
113 for (G4int iter = 0; iter < Base::GetMaxNoSteps(); ++iter)
114 {
115 const G4ThreeVector StartPos =
116 field_utils::makeVector(y, field_utils::Value3D::Position);
117
118 OneGoodStep(y, dydx, curveLength, h, eps, hdid, hnext);
119
120 const G4ThreeVector EndPos =
121 field_utils::makeVector(y, field_utils::Value3D::Position);
122
123 CheckStep(EndPos, StartPos, hdid);
124
125 G4double restCurveLength = endCurveLength - curveLength;
126 if (restCurveLength < GetSmallestFraction() * hstep)
127 {
128 succeeded = true;
129 break;
130 }
131 h = std::min(hnext, restCurveLength);
132 }
133
134
135 if (succeeded)
136 {
137 track.LoadFromArray(y, Base::GetStepper()->GetNumberOfVariables());
138 track.SetCurveLength(track.GetCurveLength() + curveLength);
139 }
140
141 return succeeded;
142}
143
144// Driver for one Runge-Kutta Step with monitoring of local truncation error
145// to ensure accuracy and adjust stepsize. Input are dependent variable
146// array y[0,...,5] and its derivative dydx[0,...,5] at the
147// starting value of the independent variable x . Also input are stepsize
148// to be attempted htry, and the required accuracy eps. On output y and x
149// are replaced by their new values, hdid is the stepsize that was actually
150// accomplished, and hnext is the estimated next stepsize.
151// This is similar to the function rkqs from the book:
152// Numerical Recipes in C: The Art of Scientific Computing (NRC), Second
153// Edition, by William H. Press, Saul A. Teukolsky, William T.
154// Vetterling, and Brian P. Flannery (Cambridge University Press 1992),
155// 16.2 Adaptive StepSize Control for Runge-Kutta, p. 719
156//
157template <class T>
158void G4FSALIntegrationDriver<T>::
159OneGoodStep(G4double y[],
160 G4double dydx[],
161 G4double& curveLength, // InOut
162 G4double htry,
163 G4double eps_rel_max,
164 G4double& hdid, // Out
165 G4double& hnext) // Out
166{
167 G4double error2 = DBL_MAX;
168
169 G4double yError[G4FieldTrack::ncompSVEC],
170 yOut[G4FieldTrack::ncompSVEC],
171 dydxOut[G4FieldTrack::ncompSVEC];
172
173 // Set stepsize to the initial trial value
174 G4double hstep = htry;
175
176 static G4ThreadLocal G4int tot_no_trials = 0;
177 const G4int max_trials = 100;
178
179 for (G4int iter = 0; iter < max_trials; ++iter)
180 {
181 ++tot_no_trials;
182
183 Base::GetStepper()->Stepper(y, dydx, hstep, yOut, yError, dydxOut);
184 error2 = field_utils::relativeError2(y, yError, hstep, eps_rel_max);
185
186 // Step succeeded.
187 if (error2 <= 1) break;
188
189 hstep = Base::ShrinkStepSize2(hstep, error2);
190 }
191
192 hnext = Base::GrowStepSize2(hstep, error2);
193 curveLength += (hdid = hstep);
194
195 for(G4int k = 0; k < Base::GetStepper()->GetNumberOfVariables(); ++k)
196 {
197 y[k] = yOut[k];
198 dydx[k] = dydxOut[k];
199 }
200}
201
202template <class T>
203G4bool G4FSALIntegrationDriver<T>::
204QuickAdvance( G4FieldTrack& fieldTrack, const G4double dydxIn[],
205 G4double hstep,
206 G4double& dchord_step, G4double& dyerr )
207{
208 ++fNoQuickAvanceCalls;
209
210 if (hstep == 0)
211 {
212 std::ostringstream message;
213 message << "Proposed step is zero; hstep = " << hstep << " !";
214 G4Exception("G4FSALIntegrationDriver ::QuickAdvance()",
215 "GeomField1001", JustWarning, message);
216 return true;
217 }
218 if (hstep < 0)
219 {
220 std::ostringstream message;
221 message << "Invalid run condition." << G4endl
222 << "Proposed step is negative; hstep = "
223 << hstep << "." << G4endl
224 << "Requested step cannot be negative! Aborting event.";
225 G4Exception("G4FSALIntegrationDriver ::QuickAdvance()",
226 "GeomField0003", EventMustBeAborted, message);
227 return false;
228 }
229
230 G4double yError[G4FieldTrack::ncompSVEC],
231 yIn[G4FieldTrack::ncompSVEC],
232 yOut[G4FieldTrack::ncompSVEC],
233 dydxOut[G4FieldTrack::ncompSVEC];
234
235 fieldTrack.DumpToArray(yIn);
236
237 Base::GetStepper()->Stepper(yIn, dydxIn, hstep, yOut, yError, dydxOut);
238 dchord_step = Base::GetStepper()->DistChord();
239
240 fieldTrack.LoadFromArray(yOut, Base::GetStepper()->GetNumberOfVariables());
241 fieldTrack.SetCurveLength(fieldTrack.GetCurveLength() + hstep);
242
243 dyerr = field_utils::absoluteError(yOut, yError, hstep);
244
245 return true;
246}
247
248template <class T>
249void G4FSALIntegrationDriver<T>::SetSmallestFraction(G4double newFraction)
250{
251 if( newFraction > 1.e-16 && newFraction < 1e-8 )
252 {
253 fSmallestFraction = newFraction;
254 }
255 else
256 {
257 std::ostringstream message;
258 message << "Smallest Fraction not changed. " << G4endl
259 << " Proposed value was " << newFraction << G4endl
260 << " Value must be between 1.e-8 and 1.e-16";
261 G4Exception("G4FSALIntegrationDriver::SetSmallestFraction()",
262 "GeomField1001", JustWarning, message);
263 }
264}
265
266template <class T>
267void G4FSALIntegrationDriver<T>::CheckStep(
268 const G4ThreeVector& posIn, const G4ThreeVector& posOut, G4double hdid)
269{
270 const G4double endPointDist = (posOut - posIn).mag();
271 if (endPointDist >= hdid * (1. + CLHEP::perMillion))
272 {
273 ++fNoAccurateAdvanceBadSteps;
274#ifdef G4DEBUG_FIELD
275 // Issue a warning only for gross differences -
276 // we understand how small difference occur.
277 if (endPointDist >= hdid * (1. + perThousand))
278 {
279 G4Exception("G4FSALIntegrationDriver::CheckStep()",
280 "GeomField1002", JustWarning,
281 "endPointDist >= hdid!");
282 }
283#endif
284 }
285 else
286 {
287 ++fNoAccurateAdvanceGoodSteps;
288 }
289}
290
291template <class T>
292inline G4double G4FSALIntegrationDriver<T>::GetMinimumStep() const
293{
294 return fMinimumStep;
295}
296
297template <class T>
298void G4FSALIntegrationDriver<T>::SetMinimumStep(G4double minimumStepLength)
299{
300 fMinimumStep = minimumStepLength;
301}
302
303template <class T>
304G4int G4FSALIntegrationDriver<T>::GetVerboseLevel() const
305{
306 return fVerboseLevel;
307}
308
309template <class T>
310void G4FSALIntegrationDriver<T>::SetVerboseLevel(G4int newLevel)
311{
312 fVerboseLevel = newLevel;
313}
314
315template <class T>
316G4double G4FSALIntegrationDriver<T>::GetSmallestFraction() const
317{
318 return fSmallestFraction;
319}
320
321template <class T>
322void G4FSALIntegrationDriver<T>::IncrementQuickAdvanceCalls()
323{
324 ++fNoQuickAvanceCalls;
325}
326
327template <class T>
328G4double
329G4FSALIntegrationDriver<T>::AdvanceChordLimited(G4FieldTrack& track,
330 G4double hstep,
331 G4double eps,
332 G4double chordDistance)
333{
334 return ChordFinderDelegate::AdvanceChordLimitedImpl(track, hstep,
335 eps, chordDistance);
336}
337
338template <class T>
339void G4FSALIntegrationDriver<T>::StreamInfo( std::ostream& os ) const
340{
341// Write out the parameters / state of the driver
342 os << "State of G4IntegrationDriver: " << std::endl;
343 os << "--Base state (G4RKIntegrationDriver): " << std::endl;
344 Base::StreamInfo( os );
345 os << "--Own state (G4IntegrationDriver<>): " << std::endl;
346 os << " fMinimumStep = " << fMinimumStep << std::endl;
347 os << " Smallest Fraction = " << fSmallestFraction << std::endl;
348
349 os << " verbose level = " << fVerboseLevel << std::endl;
350 os << " Reintegrates = " << DoesReIntegrate() << std::endl;
351 os << "--Chord Finder Delegate state: " << std::endl;
352 ChordFinderDelegate::StreamDelegateInfo( os );
353}
Generated by doxygen 1.9.3