ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4TsitourasRK45.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4TsitourasRK45.cc
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 // G4TsitourasRK45 implementation
27 //
28 // Author: Somnath Banerjee, Google Summer of Code 2015, 11.06.2015
29 // Supervision: John Apostolakis, CERN
30 // -------------------------------------------------------------------
31 
32 #include "G4TsitourasRK45.hh"
33 #include "G4LineSection.hh"
34 
36 //
37 // Constructor
38 //
39 G4TsitourasRK45::G4TsitourasRK45(G4EquationOfMotion *EqRhs,
40  G4int noIntegrationVariables,
41  G4bool primary)
42  : G4MagIntegratorStepper(EqRhs, noIntegrationVariables)
43 {
44  const G4int numberOfVariables = noIntegrationVariables;
45 
46  ak2 = new G4double[numberOfVariables] ;
47  ak3 = new G4double[numberOfVariables] ;
48  ak4 = new G4double[numberOfVariables] ;
49  ak5 = new G4double[numberOfVariables] ;
50  ak6 = new G4double[numberOfVariables] ;
51  ak7 = new G4double[numberOfVariables] ;
52  ak8 = new G4double[numberOfVariables] ;
53 
54  // Must ensure space extra 'state' variables exists - i.e. yIn[7]
55  //
56  const G4int numStateMax = std::max(GetNumberOfStateVariables(), 8);
57  const G4int numStateVars = std::max(noIntegrationVariables,
58  numStateMax );
59 
60  yTemp = new G4double[numStateVars] ;
61  yIn = new G4double[numStateVars] ;
62 
63  fLastInitialVector = new G4double[numberOfVariables] ;
64  fLastFinalVector = new G4double[numberOfVariables] ;
65 
66  fLastDyDx = new G4double[numberOfVariables];
67 
68  fMidVector = new G4double[numberOfVariables];
69  fMidError = new G4double[numberOfVariables];
70 
71  if( primary )
72  {
73  fAuxStepper = new G4TsitourasRK45(EqRhs, numberOfVariables, !primary);
74  }
75 }
76 
78 //
79 // Destructor
80 //
81 G4TsitourasRK45::~G4TsitourasRK45()
82 {
83  delete [] ak2;
84  delete [] ak3;
85  delete [] ak4;
86  delete [] ak5;
87  delete [] ak6;
88  delete [] ak7;
89  delete [] ak8;
90 
91  delete [] yTemp;
92  delete [] yIn;
93 
94  delete [] fLastInitialVector;
95  delete [] fLastFinalVector;
96  delete [] fLastDyDx;
97  delete [] fMidVector;
98  delete [] fMidError;
99 
100  delete fAuxStepper;
101 }
102 
103 // The following coefficients have been obtained from
104 // Table 1: The Coefficients of the new pair
105 //
106 // C. Tsitouras, "Runge–Kutta pairs of order 5(4) satisfying only
107 // the first column simplifying assumption"
108 // Computers & Mathematics with Applications, vol.62, no.2, pp.770-775, 2011.
109 //
110 // A corresponding matlab code was also found at:
111 // http://users.ntua.gr/tsitoura/new54.m
112 //
113 // Doing a step
114 //
115 void
116 G4TsitourasRK45::Stepper( const G4double yInput[],
117  const G4double dydx[],
118  G4double Step,
119  G4double yOut[],
120  G4double yErr[] )
121 {
122  const G4double b21 = 0.161 ,
123  b31 = -0.00848065549235698854 ,
124  b32 = 0.335480655492356989 ,
125 
126  b41 = 2.89715305710549343 ,
127  b42 = -6.35944848997507484 ,
128  b43 = 4.36229543286958141 ,
129 
130  b51 = 5.325864828439257,
131  b52 = -11.748883564062828,
132  b53 = 7.49553934288983621 ,
133  b54 = -0.09249506636175525,
134 
135  b61 = 5.8614554429464200,
136  b62 = -12.9209693178471093 ,
137  b63 = 8.1593678985761586 ,
138  b64 = -0.071584973281400997,
139  b65 = -0.0282690503940683829,
140 
141  b71 = 0.0964607668180652295 ,
142  b72 = 0.01,
143  b73 = 0.479889650414499575,
144  b74 = 1.37900857410374189,
145  b75 = -3.2900695154360807,
146  b76 = 2.32471052409977398,
147 
148  // c1 = 0.001780011052226 ,
149  // c2 = 0.000816434459657 ,
150  // c3 = -0.007880878010262 ,
151  // c4 = 0.144711007173263 ,
152  // c5 = -0.582357165452555 ,
153  // c6 = 0.458082105929187 ,
154  // c7 = 1.0/66.0 ;
155 
156  dc1 = 0.0935237485818927066 - b71 , // - 0.001780011052226,
157  dc2 = 0.00865288314156636761 - b72, // - 0.000816434459657,
158  dc3 = 0.492893099131431868 - b73 , // + 0.007880878010262,
159  dc4 = 1.14023541226785810 - b74 , // 0.144711007173263,
160  dc5 = - 2.3291801924393646 - b75, // + 0.582357165452555,
161  dc6 = 1.56887504931661552 - b76 , // - 0.458082105929187,
162  dc7 = 0.025; //- 1.0/66.0 ;
163 
164  // dc1 = -3.0/1280.0,
165  // dc2 = 0.0,
166  // dc3 = 6561.0/632320.0,
167  // dc4 = -343.0/20800.0,
168  // dc5 = 243.0/12800.0,
169  // dc6 = -1.0/95.0,
170  // dc7 = 0.0 ;
171 
172  const G4int numberOfVariables = GetNumberOfVariables();
173 
174  // The number of variables to be integrated over
175  //
176  yOut[7] = yTemp[7] = yIn[7] = yInput[7];
177 
178  // Saving yInput because yInput and yOut can be aliases for same array
179  //
180  for(G4int i=0; i<numberOfVariables; ++i)
181  {
182  yIn[i]=yInput[i];
183  }
184  // RightHandSide(yIn, dydx) ; // 1st Step - Not doing, getting passed
185 
186  for(G4int i=0; i<numberOfVariables; ++i)
187  {
188  yTemp[i] = yIn[i] + b21*Step*dydx[i] ;
189  }
190  RightHandSide(yTemp, ak2) ; // 2nd Stage
191 
192  for(G4int i=0; i<numberOfVariables; ++i)
193  {
194  yTemp[i] = yIn[i] + Step*(b31*dydx[i] + b32*ak2[i]) ;
195  }
196  RightHandSide(yTemp, ak3) ; // 3rd Stage
197 
198  for(G4int i=0; i<numberOfVariables; ++i)
199  {
200  yTemp[i] = yIn[i] + Step*(b41*dydx[i] + b42*ak2[i] + b43*ak3[i]) ;
201  }
202  RightHandSide(yTemp, ak4) ; // 4th Stage
203 
204  for(G4int i=0; i<numberOfVariables; ++i)
205  {
206  yTemp[i] = yIn[i] + Step*(b51*dydx[i] + b52*ak2[i] + b53*ak3[i] +
207  b54*ak4[i]) ;
208  }
209  RightHandSide(yTemp, ak5) ; // 5th Stage
210 
211  for(G4int i=0; i<numberOfVariables; ++i)
212  {
213  yTemp[i] = yIn[i] + Step*(b61*dydx[i] + b62*ak2[i] + b63*ak3[i] +
214  b64*ak4[i] + b65*ak5[i]) ;
215  }
216  RightHandSide(yTemp, ak6) ; // 6th Stage
217 
218  for(G4int i=0; i<numberOfVariables; ++i)
219  {
220  yOut[i] = yIn[i] + Step*(b71*dydx[i] + b72*ak2[i] + b73*ak3[i] +
221  b74*ak4[i] + b75*ak5[i] + b76*ak6[i]);
222  }
223  RightHandSide(yOut, ak7); // 7th Stage
224 
225  //Calculate the error in the step:
226  for(G4int i=0; i<numberOfVariables; ++i)
227  {
228  yErr[i] = Step*(dc1*dydx[i] + dc2*ak2[i] + dc3*ak3[i] + dc4*ak4[i] +
229  dc5*ak5[i] + dc6*ak6[i] + dc7*ak7[i] ) ;
230 
231  // Store Input and Final values, for possible use in calculating chord
232  //
233  fLastInitialVector[i] = yIn[i] ;
234  fLastFinalVector[i] = yOut[i];
235  fLastDyDx[i] = dydx[i];
236  }
237 
238  fLastStepLength = Step;
239 
240  return ;
241 }
242 
243 void G4TsitourasRK45::SetupInterpolation()
244  // (const G4double *yInput, const G4double *dydx, const G4double Step)
245 {
246  // Nothing to be done
247 }
248 
249 void G4TsitourasRK45::Interpolate(const G4double* yInput,
250  const G4double* dydx,
251  const G4double Step,
252  G4double* yOut,
253  G4double tau)
254 {
255  G4double bf1, bf2, bf3, bf4, bf5, bf6, bf7;
256  // Coefficients for all the seven stages.
257 
258  const G4int numberOfVariables = GetNumberOfVariables();
259 
260  G4double tau0 = tau;
261 
262  for(G4int i=0; i<numberOfVariables; ++i)
263  {
264  yIn[i] = yInput[i];
265  }
266 
267  G4double tau_2 = tau0*tau0 ;
268  // tau_3 = tau0*tau_2,
269  // tau_4 = tau_2*tau_2;
270 
271  bf1 = -1.0530884977290216*tau*(tau - 1.3299890189751412)*(tau_2 -
272  1.4364028541716351*tau + 0.7139816917074209);
273  bf2 = 0.1017*tau_2*(tau_2 - 2.1966568338249754*tau +
274  1.2949852507374631);
275  bf3 = 2.490627285651252793*tau_2*(tau_2 - 2.38535645472061657*tau
276  + 1.57803468208092486);
277  bf4 = -16.54810288924490272*(tau - 1.21712927295533244)*
278  (tau - 0.61620406037800089)*tau_2;
279  bf5 = 47.37952196281928122*(tau - 1.203071208372362603)*
280  (tau - 0.658047292653547382)*tau_2;
281  bf6 = -34.87065786149660974*(tau - 1.2)*(tau -
282  0.666666666666666667)*tau_2;
283  bf7 = 2.5*(tau - 1.0)*(tau - 0.6)*tau_2;
284 
285  // Putting together the coefficients calculated as the respective
286  // stage coefficients
287  //
288  for(G4int i=0; i<numberOfVariables; ++i)
289  {
290  yOut[i] = yIn[i] + Step*( bf1*dydx[i] + bf2*ak2[i] + bf3*ak3[i]
291  + bf4*ak4[i] + bf5*ak5[i] + bf6*ak6[i]
292  + bf7*ak7[i] ) ;
293  }
294 }
295 
296 G4double G4TsitourasRK45::DistChord() const
297 {
298  G4double distLine, distChord;
299  G4ThreeVector initialPoint, finalPoint, midPoint;
300 
301  // Store last initial and final points (they will be
302  // overwritten in self-Stepper call!)
303  //
304  initialPoint = G4ThreeVector( fLastInitialVector[0],
305  fLastInitialVector[1], fLastInitialVector[2]);
306  finalPoint = G4ThreeVector( fLastFinalVector[0],
307  fLastFinalVector[1], fLastFinalVector[2]);
308 
309  // Do half a step using StepNoErr
310  //
311  fAuxStepper->Stepper( fLastInitialVector, fLastDyDx, 0.5 * fLastStepLength,
312  fMidVector, fMidError );
313 
314  midPoint = G4ThreeVector( fMidVector[0], fMidVector[1], fMidVector[2]);
315 
316  // Use stored values of Initial and Endpoint + new Midpoint to evaluate
317  // distance of Chord
318  //
319  if (initialPoint != finalPoint)
320  {
321  distLine = G4LineSection::Distline( midPoint, initialPoint, finalPoint );
322  distChord = distLine;
323  }
324  else
325  {
326  distChord = (midPoint-initialPoint).mag();
327  }
328  return distChord;
329 }