| Geant4 Cross Reference |
1 // ******************************************* 1
2 // * License and Disclaimer
3 // *
4 // * The Geant4 software is copyright of th
5 // * the Geant4 Collaboration. It is provided
6 // * conditions of the Geant4 Software License
7 // * LICENSE and available at http://cern.ch/
8 // * include a list of copyright holders.
9 // *
10 // * Neither the authors of this software syst
11 // * institutes,nor the agencies providing fin
12 // * work make any representation or warran
13 // * regarding this software system or assum
14 // * use. Please see the license in the file
15 // * for the full disclaimer and the limitatio
16 // *
17 // * This code implementation is the result
18 // * technical work of the GEANT4 collaboratio
19 // * By using, copying, modifying or distri
20 // * any work based on the software) you ag
21 // * use in resulting scientific publicati
22 // * acceptance of all terms of the Geant4 Sof
23 // *******************************************
24 //
25 // G4BorisDriver
26 //
27 // Class description:
28 //
29 // G4BorisDriver is a driver class using the
30 // method to integrate the equation of motion.
31 //
32 //
33 // Author: Divyansh Tiwari, Google Summer of C
34 // Supervision: John Apostolakis,Renee Fatemi,
35 // -------------------------------------------
36 #include <cassert>
37
38 #include "G4BorisDriver.hh"
39
40 #include "G4SystemOfUnits.hh"
41 #include "G4LineSection.hh"
42 #include "G4FieldUtils.hh"
43
44 const G4double G4BorisDriver::fErrorConstraint
45 fMaxSteppingDecrease / fSafetyFact
46
47 const G4double G4BorisDriver::fErrorConstraint
48 fMaxSteppingIncrease / fSafetyFact
49
50 // -------------------------------------------
51
52 G4BorisDriver::
53 G4BorisDriver( G4double hminimum, G4BorisSchem
54 G4int numberOfComponents, G4boo
55 : fMinimumStep(hminimum),
56 fVerbosity(verbosity),
57 boris(Boris)
58 // , interval_sequence{2,4}
59 {
60 assert(boris->GetNumberOfVariables() == nu
61
62 if(boris->GetNumberOfVariables() != number
63 {
64 std::ostringstream msg;
65 msg << "Disagreement in number of variab
66 << boris->GetNumberOfVariables()
67 << " vs no of components = " << numb
68 G4Exception("G4BorisDriver Constructor:"
69 "GeomField1001", FatalExcept
70 }
71
72 }
73
74 // -------------------------------------------
75
76 G4bool G4BorisDriver::AccurateAdvance( G4Field
77 G4doubl
78 G4doubl
79 G4doubl
80 {
81 // Specification: Driver with adaptive step
82 // Integrate starting values at y_current o
83 // On output 'track' is replaced by values
84
85 // Ensure that hstep > 0
86 if(hstep == 0)
87 {
88 std::ostringstream message;
89 message << "Proposed step is zero; hstep
90 G4Exception("G4BorisDriver::AccurateAdva
91 "GeomField1001", JustWarning
92 return true;
93 }
94 if(hstep < 0)
95 {
96 std::ostringstream message;
97 message << "Invalid run condition." << G
98 << "Proposed step is negative; h
99 << "Requested step cannot be neg
100 G4Exception("G4BorisDriver::AccurateAdva
101 "GeomField0003", EventMustBe
102 return false;
103 }
104
105 if( hinitial == 0.0 ) { hinitial = hstep; }
106 if( hinitial < 0.0 ) { hinitial = std::fabs
107 // G4double htrial = std::min( hstep, hinit
108 G4double htrial = hstep;
109 // Decide first step size
110
111 // G4int noOfSteps = h/hstep;
112
113 // integration variables
114 //
115 track.DumpToArray(yCurrent);
116
117 const G4double restMass = track.GetRestMass
118 const G4double charge = track.GetCharge()*e
119 const G4int nvar= GetNumberOfVariables()
120
121 // copy non-integration variables to out ar
122 //
123 std::memcpy(yOut + nvar,
124 yCurrent + nvar,
125 sizeof(G4double)*(G4FieldTrack:
126
127 G4double curveLength = track.GetCurveLength
128 const G4double endCurveLength = curveLength
129
130 // -- Initial version: Did it in one step -
131 // G4FieldTrack yFldTrk(track);
132 // yFldTrk.LoadFromArray(yCurrent, G4FieldT
133 // yFldTrk.SetCurveLength(curveLength);
134 // G4double dchord_step, dyerr_len;
135 // QuickAdvance(yFldTrk, dydxCurrent, htria
136
137 const G4double hThreshold =
138 std::max(epsilon * hstep, fSmallestFra
139
140 G4double htry= htrial;
141
142 for (G4int nstp = 0; nstp < fMaxNoSteps; ++
143 {
144 G4double hdid= 0.0, hnext=0.0;
145
146 OneGoodStep(yCurrent, curveLength, htry,
147 //*********
148
149 // Simple check: move (distance of displ
150 const G4ThreeVector StartPos = field_uti
151 const G4ThreeVector EndPos = field_utils
152 CheckStep(EndPos, StartPos, hdid);
153
154 // Check 1. for finish and 2. *av
155 if (curveLength >= endCurveLength || htr
156 {
157 break;
158 }
159
160 htry = std::max(hnext, fMinimumStep);
161 if (curveLength + htry > endCurveLength)
162 {
163 htry = endCurveLength - curveLength;
164 }
165 }
166
167 // upload new state
168 track.LoadFromArray(yCurrent, G4FieldTrack:
169 track.SetCurveLength(curveLength);
170
171 return true;
172 }
173
174 // -------------------------------------------
175
176 void G4BorisDriver::OneGoodStep(G4double y[],
177 G4double& curv
178 G4double htry
179 G4double epsi
180 G4double rest
181 G4double char
182 G4double& hdid
183 G4double& hnex
184 {
185 G4double error2 = DBL_MAX;
186 G4double yerr[G4FieldTrack::ncompSVEC], yt
187
188 G4double h = htry;
189
190 const G4int max_trials = 100;
191
192 for (G4int iter = 0; iter < max_trials; ++
193 {
194 boris->StepWithErrorEstimate(y, restMa
195
196 error2 = field_utils::relativeError2(y
197 e
198 if (error2 <= 1.0)
199 {
200 break;
201 }
202
203 h = ShrinkStepSize2(h, error2);
204
205 G4double xnew = curveLength + h;
206 if(xnew == curveLength)
207 {
208 std::ostringstream message;
209 message << "Stepsize underflow in
210 << "- Step's start x=" <<
211 << " and end x= " << xnew
212 << " are equal !! " << G4e
213 << " Due to step-size= "
214 << ". Note that input step
215 G4Exception("G4IntegrationDriver::
216 "GeomField1001", JustW
217 break;
218 }
219 }
220
221 hnext = GrowStepSize2(h, error2);
222 curveLength += (hdid = h);
223
224 field_utils::copy(y, ytemp, GetNumberOfVar
225 }
226
227 // ===========--------------------------------
228
229 G4bool G4BorisDriver::
230 QuickAdvance( G4FieldTrack& track, const G4dou
231 G4double hstep, G4double& missDi
232 {
233 const auto nvar = boris->GetNumberOfVariab
234
235 track.DumpToArray(yIn);
236 const G4double curveLength = track.GetCurv
237
238 // call the boris method for step length h
239 G4double restMass = track.GetRestMass();
240 G4double charge = track.GetCharge()*e_SI;
241
242 // boris->DoStep(restMass, charge, yIn, y
243 // boris->DoStep(restMass, charge, yMid, y
244 boris->StepWithMidAndErrorEstimate(yIn, re
245 yMid, y
246 // Same, and also return mid-point evalu
247
248 // How to calculate chord length??
249 const auto mid = field_utils::makeVector(y
250 field_utils::Value3D::Pos
251 const auto in = field_utils::makeVector(y
252 field_utils::Value3D::Pos
253 const auto out = field_utils::makeVector(y
254 field_utils::Value3D::Pos
255
256 missDist = G4LineSection::Distline(mid, in
257
258 dyerr = field_utils::absoluteError(yOut, y
259
260 // copy non-integrated variables to output
261 //
262 std::memcpy(yOut + nvar, yIn + nvar,
263 sizeof(G4double) * (G4FieldTra
264
265 // set new state
266 //
267 track.LoadFromArray(yOut, G4FieldTrack::nc
268 track.SetCurveLength(curveLength + hstep)
269
270 return true;
271 }
272
273 // -------------------------------------------
274
275 void G4BorisDriver::
276 GetDerivatives( const G4FieldTrack& yTrack, G4
277 {
278 G4double EBfieldValue[6];
279 GetDerivatives(yTrack, dydx, EBfieldValue);
280 }
281
282 // -------------------------------------------
283
284 void G4BorisDriver::
285 GetDerivatives( const G4FieldTrack& yTrack, G4
286 G4double EBfieldValue[]) const
287 {
288 // G4Exception("G4BorisDriver::GetDerivativ
289 // "GeomField0003", FatalExcept
290 G4double ytemp[G4FieldTrack::ncompSVEC];
291 yTrack.DumpToArray(ytemp);
292 GetEquationOfMotion()->EvaluateRhsReturnB(y
293 }
294
295 // -------------------------------------------
296
297 G4double G4BorisDriver::ShrinkStepSize2(G4doub
298 {
299 if (error2 > fErrorConstraintShrink * fErro
300 {
301 return fMaxSteppingDecrease * h;
302 }
303 return fSafetyFactor * h * std::pow(error2,
304 }
305
306 // -------------------------------------------
307
308 G4double G4BorisDriver::GrowStepSize2(G4double
309 // Given the square of the relative error,
310 {
311 if (error2 < fErrorConstraintGrow * fError
312 {
313 return fMaxSteppingIncrease * h;
314 }
315 return fSafetyFactor * h * std::pow(error2
316 }
317
318 // -------------------------------------------
319
320 void G4BorisDriver::SetEquationOfMotion(G4Equa
321 {
322 G4Exception("G4BorisDriver::SetEquationOfMo
323 "This method is not implemented
324 }
325
326 // -------------------------------------------
327
328 void
329 G4BorisDriver::StreamInfo( std::ostream& os )
330 {
331 os << "State of G4BorisDriver: " << std::end
332 os << " Method is implemented, but gives n
333 }
334