| Geant4 Cross Reference |
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 // G4OldMagIntDriver
27 //
28 // Class description:
29 //
30 // Provides a driver that talks to the Integrator Stepper, and insures that
31 // the error is within acceptable bounds.
32
33 // V.Grichine, 07.10.1996 - Created
34 // W.Wander, 28.01.1998 - Added ability for low order integrators
35 // J.Apostolakis, 08.11.2001 - Respect minimum step in AccurateAdvance
36 // --------------------------------------------------------------------
37 #ifndef G4OLD_MAGINT_DRIVER_HH
38 #define G4OLD_MAGINT_DRIVER_HH
39
40 #include "G4VIntegrationDriver.hh"
41 #include "G4MagIntegratorStepper.hh"
42 #include "G4ChordFinderDelegate.hh"
43
44 class G4OldMagIntDriver : public G4VIntegrationDriver,
45 public G4ChordFinderDelegate<G4OldMagIntDriver>
46 {
47 public:
48
49 G4OldMagIntDriver(G4double hminimum,
50 G4MagIntegratorStepper* pItsStepper,
51 G4int numberOfComponents = 6,
52 G4int statisticsVerbosity = 0);
53 ~G4OldMagIntDriver() override;
54 // Constructor, destructor.
55
56 G4OldMagIntDriver(const G4OldMagIntDriver&) = delete;
57 G4OldMagIntDriver& operator=(const G4OldMagIntDriver&) = delete;
58
59 inline G4double AdvanceChordLimited(G4FieldTrack& track,
60 G4double stepMax,
61 G4double epsStep,
62 G4double chordDistance) override;
63
64 inline void OnStartTracking() override;
65 inline void OnComputeStep(const G4FieldTrack* = nullptr) override {}
66 inline G4bool DoesReIntegrate() const override { return true; }
67
68 G4bool AccurateAdvance(G4FieldTrack& y_current,
69 G4double hstep,
70 G4double eps, // Requested y_err/hstep
71 G4double hinitial = 0.0) override;
72 // Above drivers for integrator (Runge-Kutta) with stepsize control.
73 // Integrates ODE starting values y_current
74 // from current s (s=s0) to s=s0+h with accuracy eps.
75 // On output ystart is replaced by value at end of interval.
76 // The concept is similar to the odeint routine from NRC p.721-722.
77
78 G4bool QuickAdvance(G4FieldTrack& y_val, // INOUT
79 const G4double dydx[],
80 G4double hstep,
81 G4double& dchord_step,
82 G4double& dyerr) override;
83 // QuickAdvance just tries one Step - it does not ensure accuracy.
84
85 G4bool QuickAdvance( G4FieldTrack& y_posvel, // INOUT
86 const G4double dydx[],
87 G4double hstep, // IN
88 G4double& dchord_step,
89 G4double& dyerr_pos_sq,
90 G4double& dyerr_mom_rel_sq);
91 // QuickAdvance that also just tries one Step (so also does not
92 // ensure accuracy), but does return the errors in position and
93 // momentum (normalised: Delta_Integration(p^2)/(p^2) ).
94
95 inline G4double GetHmin() const;
96 inline G4double Hmin() const; // Obsolete
97 inline G4double GetSafety() const;
98 inline G4double GetPshrnk() const;
99 inline G4double GetPgrow() const;
100 inline G4double GetErrcon() const;
101 void GetDerivatives(const G4FieldTrack& y_curr, // INput
102 G4double dydx[]) const override; // OUTput
103
104 void GetDerivatives(const G4FieldTrack& track,
105 G4double dydx[],
106 G4double field[]) const override;
107 // Accessors
108
109 G4EquationOfMotion* GetEquationOfMotion() override;
110 void SetEquationOfMotion(G4EquationOfMotion* equation) override;
111
112 void RenewStepperAndAdjust(G4MagIntegratorStepper* pItsStepper) override;
113 // Sets a new stepper pItsStepper for this driver. Then it calls
114 // ReSetParameters to reset its parameters accordingly.
115
116 inline void ReSetParameters(G4double new_safety = 0.9);
117 // i) sets the exponents (pgrow & pshrnk),
118 // using the current Stepper's order,
119 // ii) sets the safety
120 // ii) calculates "errcon" according to the above values.
121
122 inline void SetSafety(G4double valS);
123 inline void SetPshrnk(G4double valPs);
124 inline void SetPgrow (G4double valPg);
125 inline void SetErrcon(G4double valEc);
126 // When setting safety or pgrow, errcon will be set to a compatible value.
127
128 inline G4double ComputeAndSetErrcon();
129
130 const G4MagIntegratorStepper* GetStepper() const override;
131 G4MagIntegratorStepper* GetStepper() override;
132
133 void OneGoodStep( G4double ystart[], // Like old RKF45step()
134 const G4double dydx[],
135 G4double& x,
136 G4double htry,
137 G4double eps, // memb variables ?
138 G4double& hdid,
139 G4double& hnext) ;
140 // This takes one Step that is as large as possible while
141 // satisfying the accuracy criterion of:
142 // yerr < eps * |y_end-y_start|
143
144 G4double ComputeNewStepSize(G4double errMaxNorm, // normalised
145 G4double hstepCurrent) override;
146 // Taking the last step's normalised error, calculate
147 // a step size for the next step.
148 // Do not limit the next step's size within a factor of the
149 // current one.
150
151 void StreamInfo( std::ostream& os ) const override;
152
153 G4double ComputeNewStepSize_WithinLimits(G4double errMaxNorm, // normalised
154 G4double hstepCurrent);
155 // Taking the last step's normalised error, calculate
156 // a step size for the next step.
157 // Limit the next step's size within a range around the current one.
158
159 inline G4int GetMaxNoSteps() const;
160 inline void SetMaxNoSteps(G4int val);
161 // Modify and Get the Maximum number of Steps that can be
162 // taken for the integration of a single segment -
163 // (i.e. a single call to AccurateAdvance).
164
165 inline void SetHmin(G4double newval);
166 void SetVerboseLevel(G4int newLevel) override;
167 G4int GetVerboseLevel() const override;
168
169 inline G4double GetSmallestFraction() const;
170 void SetSmallestFraction( G4double val );
171
172 protected:
173
174 void WarnSmallStepSize(G4double hnext, G4double hstep,
175 G4double h, G4double xDone,
176 G4int noSteps);
177
178 void WarnTooManyStep(G4double x1start, G4double x2end, G4double xCurrent);
179 void WarnEndPointTooFar(G4double endPointDist,
180 G4double hStepSize ,
181 G4double epsilonRelative,
182 G4int debugFlag);
183 // Issue warnings for undesirable situations
184
185 void PrintStatus(const G4double* StartArr,
186 G4double xstart,
187 const G4double* CurrentArr,
188 G4double xcurrent,
189 G4double requestStep,
190 G4int subStepNo);
191 void PrintStatus(const G4FieldTrack& StartFT,
192 const G4FieldTrack& CurrentFT,
193 G4double requestStep,
194 G4int subStepNo);
195 void PrintStat_Aux(const G4FieldTrack& aFieldTrack,
196 G4double requestStep,
197 G4double actualStep,
198 G4int subStepNo,
199 G4double subStepSize,
200 G4double dotVelocities);
201 // Verbose output for debugging
202
203 void PrintStatisticsReport();
204 // Report on the number of steps, maximum errors etc.
205
206 #ifdef QUICK_ADV_TWO
207 G4bool QuickAdvance( G4double yarrin[], // In
208 const G4double dydx[],
209 G4double hstep,
210 G4double yarrout[], // Out
211 G4double& dchord_step, // Out
212 G4double& dyerr ); // in length
213 #endif
214
215 private:
216
217 // ---------------------------------------------------------------
218 // INVARIANTS
219
220 G4double fMinimumStep = 0.0;
221 // Minimum Step allowed in a Step (in absolute units)
222 G4double fSmallestFraction = 1.0e-12; // Expected range 1e-12 to 5e-15
223 // Smallest fraction of (existing) curve length - in relative units
224 // below this fraction the current step will be the last
225
226 const G4int fNoIntegrationVariables = 0; // Variables in integration
227 const G4int fMinNoVars = 12; // Minimum number for FieldTrack
228 const G4int fNoVars = 0; // Full number of variable
229
230 G4int fMaxNoSteps;
231 G4int fMaxStepBase = 250; // was 5000
232 // Default maximum number of steps is Base divided by the order of Stepper
233
234 G4double safety;
235 G4double pshrnk; // exponent for shrinking
236 G4double pgrow; // exponent for growth
237 G4double errcon;
238 // Parameters used to grow and shrink trial stepsize.
239
240 G4int fStatisticsVerboseLevel = 0;
241
242 // ---------------------------------------------------------------
243 // DEPENDENT Objects
244
245 G4MagIntegratorStepper* pIntStepper = nullptr;
246
247 // ---------------------------------------------------------------
248 // STATE
249
250 unsigned long fNoTotalSteps=0, fNoBadSteps=0;
251 unsigned long fNoSmallSteps=0, fNoInitialSmallSteps=0, fNoCalls=0;
252 G4double fDyerr_max=0.0, fDyerr_mx2=0.0;
253 G4double fDyerrPos_smTot=0.0, fDyerrPos_lgTot=0.0, fDyerrVel_lgTot=0.0;
254 G4double fSumH_sm=0.0, fSumH_lg=0.0;
255 // Step Statistics
256
257 G4int fVerboseLevel = 0; // Verbosity level for printing (debug, ..)
258 // Could be varied during tracking - to help identify issues
259
260 using ChordFinderDelegate = G4ChordFinderDelegate<G4OldMagIntDriver>;
261 };
262
263 #include "G4OldMagIntDriver.icc"
264
265 #endif
266