Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/geometry/magneticfield/src/G4DoLoMcPriRK34.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

  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 // G4DoLoMcPriRK34 implementation
 27 //
 28 // Created: Somnath Banerjee, Google Summer of Code 2015, 7 July 2015
 29 // Supervision: John Apostolakis, CERN
 30 // --------------------------------------------------------------------
 31 
 32 #include "G4DoLoMcPriRK34.hh"
 33 #include "G4LineSection.hh"
 34 
 35 // Constructor
 36 //
 37 G4DoLoMcPriRK34::G4DoLoMcPriRK34(G4EquationOfMotion* EqRhs,
 38                                   G4int noIntegrationVariables,
 39                                   G4bool primary)
 40    : G4MagIntegratorStepper(EqRhs, noIntegrationVariables)
 41 {
 42    const G4int numberOfVariables = noIntegrationVariables;
 43    
 44    // New Chunk of memory being created for use by the Stepper
 45    
 46    // aki - for storing intermediate RHS
 47    //
 48    ak2 = new G4double[numberOfVariables];
 49    ak3 = new G4double[numberOfVariables];
 50    ak4 = new G4double[numberOfVariables];
 51    ak5 = new G4double[numberOfVariables];
 52    ak6 = new G4double[numberOfVariables];
 53    
 54    yTemp = new G4double[numberOfVariables] ;
 55    yIn = new G4double[numberOfVariables] ;
 56       
 57    fLastInitialVector = new G4double[numberOfVariables] ;
 58    fLastFinalVector = new G4double[numberOfVariables] ;
 59    fLastDyDx = new G4double[numberOfVariables];
 60    
 61    fMidVector = new G4double[numberOfVariables];
 62    fMidError = new G4double[numberOfVariables];
 63    if( primary )
 64    {
 65      fAuxStepper = new G4DoLoMcPriRK34(EqRhs, numberOfVariables, !primary);
 66    }
 67 }
 68 
 69 // Destructor
 70 //
 71 G4DoLoMcPriRK34::~G4DoLoMcPriRK34()
 72 {
 73    // clear all previously allocated memory for Stepper and DistChord
 74 
 75    delete [] ak2;
 76    delete [] ak3;
 77    delete [] ak4;
 78    delete [] ak5;
 79    delete [] ak6;
 80    
 81    delete [] yTemp;
 82    delete [] yIn;
 83    
 84    delete [] fLastInitialVector;
 85    delete [] fLastFinalVector;
 86    delete [] fLastDyDx;
 87    delete [] fMidVector;
 88    delete [] fMidError;
 89    
 90    delete fAuxStepper;
 91 }
 92 
 93 // Stepper
 94 //
 95 // Passing in the value of yInput[],the first time dydx[] and Step length
 96 // Giving back yOut and yErr arrays for output and error respectively
 97 //
 98 void G4DoLoMcPriRK34::Stepper(const G4double yInput[],
 99                               const G4double DyDx[],
100                                     G4double Step,
101                                     G4double yOut[],
102                                     G4double yErr[] )
103 {
104    G4int i;
105    
106    // The various constants defined on the basis of butcher tableu
107    //
108    const G4double b21 = 7.0/27.0 ,
109                   b31 = 7.0/72.0 , 
110                   b32 = 7.0/24.0 ,
111    
112                   b41 = 3043.0/3528.0 , 
113                   b42 = -3757.0/1176.0 ,
114                   b43 = 1445.0/441.0,
115    
116                   b51 = 17617.0/11662.0 ,
117                   b52 = -4023.0/686.0 , 
118                   b53 = 9372.0/1715.0 ,
119                   b54 = -66.0/595.0 ,
120    
121                   b61 = 29.0/238.0 ,
122                   b62 = 0.0 , 
123                   b63 = 216.0/385.0 ,
124                   b64 = 54.0/85.0 ,
125                   b65 = -7.0/22.0 ,
126 
127                   dc1 = 363.0/2975.0 - b61 ,
128                   dc2 = 0.0 - b62 ,
129                   dc3 = 981.0/1750.0 - b63,
130                   dc4 = 2709.0/4250.0 - b64 ,
131                   dc5 = -3.0/10.0 - b65 ,
132                   dc6 = -1.0/50.0 ;        // end of declaration
133 
134    const G4int numberOfVariables = GetNumberOfVariables();
135    
136    // The number of variables to be integrated over
137    //
138    yOut[7] = yTemp[7]  = yIn[7];
139 
140    // Saving yInput because yInput and yOut can be aliases for same array
141    //
142    for(i=0; i<numberOfVariables; ++i)
143    {
144        yIn[i]=yInput[i];
145    }
146 
147    // RightHandSide(yIn, DyDx) ;   // 1st stage - Not doing, getting passed
148    
149    for(i=0; i<numberOfVariables; ++i)
150    {
151        yTemp[i] = yIn[i] + b21*Step*DyDx[i] ;
152    }
153    RightHandSide(yTemp, ak2) ;              // 2nd stage
154    
155    for(i=0; i<numberOfVariables; ++i)
156    {
157        yTemp[i] = yIn[i] + Step*(b31*DyDx[i] + b32*ak2[i]) ;
158    }
159    RightHandSide(yTemp, ak3) ;              // 3rd stage
160    
161    for(i=0; i<numberOfVariables; ++i)
162    {
163        yTemp[i] = yIn[i] + Step*(b41*DyDx[i] + b42*ak2[i] + b43*ak3[i]) ;
164    }
165    RightHandSide(yTemp, ak4) ;              // 4th stage
166    
167    for(i=0; i<numberOfVariables; ++i)
168    {
169        yTemp[i] = yIn[i] + Step*(b51*DyDx[i] + b52*ak2[i]
170                                + b53*ak3[i] + b54*ak4[i]) ;
171    }
172    RightHandSide(yTemp, ak5) ;              // 5th stage
173    
174    for(i=0; i<numberOfVariables; ++i)
175    {
176        yOut[i] = yIn[i] + Step*(b61*DyDx[i] + b62*ak2[i] + b63*ak3[i]
177                               + b64*ak4[i] + b65*ak5[i]) ;
178    }
179    RightHandSide(yOut, ak6) ;              // 6th and Final stage
180 
181    for(i=0; i<numberOfVariables; ++i)
182    {
183        
184        yErr[i] = Step*(dc1*DyDx[i] + dc2*ak2[i] + dc3*ak3[i] + dc4*ak4[i]
185                      + dc5*ak5[i] + dc6*ak6[i] ) ;
186 
187        // Store Input and Final values, for possible use in calculating chord
188        //
189        fLastInitialVector[i] = yIn[i] ;
190        fLastFinalVector[i]   = yOut[i];
191        fLastDyDx[i]          = DyDx[i];
192    }
193    
194    fLastStepLength = Step;
195    
196    return ;
197 }
198 
199 // DistChord
200 //
201 G4double  G4DoLoMcPriRK34::DistChord() const
202 {
203    G4double distLine, distChord;
204    G4ThreeVector initialPoint, finalPoint, midPoint;
205    
206    // Store last initial and final points
207    // (they will be overwritten in self-Stepper call!)
208    //
209    initialPoint = G4ThreeVector( fLastInitialVector[0],
210                                  fLastInitialVector[1], fLastInitialVector[2] );
211    finalPoint   = G4ThreeVector( fLastFinalVector[0],
212                                  fLastFinalVector[1],  fLastFinalVector[2] );
213    
214    // Do half a Step using StepNoErr
215    
216    fAuxStepper->Stepper( fLastInitialVector, fLastDyDx, 0.5 * fLastStepLength,
217                          fMidVector, fMidError );
218    
219    midPoint = G4ThreeVector( fMidVector[0], fMidVector[1], fMidVector[2]);
220    
221    // Use stored values of Initial and Endpoint + new Midpoint to evaluate
222    // distance of Chord
223    //
224    if (initialPoint != finalPoint)
225    {
226      distLine  = G4LineSection::Distline( midPoint, initialPoint, finalPoint );
227      distChord = distLine;
228    }
229    else
230    {
231      distChord = (midPoint-initialPoint).mag();
232    }
233    return distChord;
234 }
235 
236 void G4DoLoMcPriRK34::SetupInterpolation()
237 {
238 }
239 
240 void G4DoLoMcPriRK34::SetupInterpolate( const G4double /* yInput */ [] ,
241                                         const G4double  /* dydx */ [] ,
242                                         const G4double  /* Step */ )
243 {
244   // Do Nothing
245 }
246 
247 void G4DoLoMcPriRK34::Interpolate( G4double tau, 
248                                    G4double yOut[] )
249 {
250   Interpolate( fLastInitialVector, fLastDyDx, fLastStepLength, yOut, tau );
251 }
252 
253 // Function to evaluate the interpolation at tau fraction of the step
254 //
255 void G4DoLoMcPriRK34::Interpolate( const G4double yInput[],
256                                    const G4double dydx[],
257                                    const G4double Step,
258                                          G4double yOut[],
259                                          G4double tau )
260 {
261    G4double bf1, bf2, bf3, bf4, bf5, bf6;    
262 
263    const G4int numberOfVariables = GetNumberOfVariables();
264    
265    for(G4int i=0; i<numberOfVariables; ++i)
266    {
267      yIn[i]=yInput[i];
268    }
269    
270    G4double tau_2 = tau*tau, tau_3 = tau*tau_2;
271    
272    // Calculating the polynomials (coefficients for the respective stages)
273    //
274    bf1 = -(162.0*tau_3 - 504.0*tau_2 + 551.0*tau - 238.0)/238.0 ,
275    bf2 =  0.0 ,
276    bf3 =  27.0*tau*(27.0*tau_2 - 70.0*tau + 51.0 )/385.0 ,
277    bf4 = -27*tau*(27.0*tau_2 - 50.0*tau + 21.0)/85.0 ,
278    bf5 =  7.0*tau*(2232.0*tau_2 - 4166.0*tau + 1785.0 )/3278.0 ,
279    bf6 = tau*(tau - 1.0)*(387.0*tau - 238.0)/149.0 ;
280    
281    for( G4int i=0; i<numberOfVariables; ++i)
282    {
283      yOut[i] = yIn[i] + Step*tau*(bf1*dydx[i] + bf2*ak2[i] + bf3*ak3[i]
284              + bf4*ak4[i] + bf5*ak5[i] + bf6*ak6[i] ) ;
285    }
286 }
287