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