| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // G4BogackiShampine23 implementation << 26 // Bogacki-Shampine - 4 - 3(2) non-FSAL implementation by Somnath Banerjee
>> 27 // Supervision / code review: John Apostolakis
27 // 28 //
28 // Bogacki-Shampine - 4 - 3(2) non-FSAL imple << 29 // Sponsored by Google in Google Summer of Code 2015.
29 // <<
30 // Implementation of the method proposed in t <<
31 // "A 3(2) pair of Runge - Kutta formulas" <<
32 // by P. Bogacki and L. F. Shampine, <<
33 // Appl. Math. Lett., vol. 2, no. 4, pp. 32 <<
34 // 30 //
35 // The Bogacki shampine method has the followi << 31 // First version: 20 May 2015
36 // <<
37 // 0 | <<
38 // 1/2|1/2 <<
39 // 3/4|0 3/4 <<
40 // 1 |2/9 1/3 4/9 <<
41 // ------------------- <<
42 // |2/9 1/3 4/9 0 <<
43 // |7/24 1/4 1/3 1/8 <<
44 // 32 //
45 // Created: Somnath Banerjee, Google Summer of << 33 // History
46 // Supervision: John Apostolakis, CERN << 34 // -----------------------------
47 // ------------------------------------------- << 35 // Created by Somnath Banerjee on 20 May 2015
>> 36 ///////////////////////////////////////////////////////////////////////////////
>> 37
>> 38
>> 39 /*
>> 40
>> 41 This contains the stepper function of the G4BogackiShampine23 class
>> 42
>> 43 The Bogacki shampine method has the following Butcher's tableau
>> 44
>> 45 0 |
>> 46 1/2|1/2
>> 47 3/4|0 3/4
>> 48 1 |2/9 1/3 4/9
>> 49 -------------------
>> 50 |2/9 1/3 4/9 0
>> 51 |7/24 1/4 1/3 1/8
>> 52
>> 53 */
48 54
49 #include "G4BogackiShampine23.hh" 55 #include "G4BogackiShampine23.hh"
50 #include "G4LineSection.hh" 56 #include "G4LineSection.hh"
51 #include "G4FieldUtils.hh" <<
52 57
53 using namespace field_utils; << 58 using namespace std;
54 59
55 G4BogackiShampine23::G4BogackiShampine23(G4Equ << 60 //Constructor
56 G4int << 61 G4BogackiShampine23::G4BogackiShampine23(G4EquationOfMotion *EqRhs,
57 : G4MagIntegratorStepper(EqRhs, integrationV << 62 G4int noIntegrationVariables,
>> 63 G4bool primary)
>> 64 : G4MagIntegratorStepper(EqRhs, noIntegrationVariables),
>> 65 fLastStepLength(0.), fAuxStepper(0)
58 { 66 {
59 SetIntegrationOrder(3); << 67 const G4int numberOfVariables = noIntegrationVariables;
60 SetFSAL(true); << 68
>> 69 ak2 = new G4double[numberOfVariables] ;
>> 70 ak3 = new G4double[numberOfVariables] ;
>> 71 ak4 = new G4double[numberOfVariables] ;
>> 72
>> 73 pseudoDydx_for_DistChord = new G4double[numberOfVariables];
>> 74
>> 75 const G4int numStateVars = std::max(noIntegrationVariables,
>> 76 GetNumberOfStateVariables() );
>> 77
>> 78 yTemp = new G4double[numberOfVariables] ;
>> 79 yIn = new G4double[numberOfVariables] ;
>> 80
>> 81 fLastInitialVector = new G4double[numStateVars] ;
>> 82 fLastFinalVector = new G4double[numStateVars] ;
>> 83 fLastDyDx = new G4double[numStateVars];
>> 84
>> 85 fMidVector = new G4double[numStateVars];
>> 86 fMidError = new G4double[numStateVars];
>> 87 if( primary )
>> 88 {
>> 89 fAuxStepper = new G4BogackiShampine23(EqRhs, numberOfVariables, !primary);
>> 90 }
61 } 91 }
62 92
63 void G4BogackiShampine23::makeStep(const G4dou << 93
64 const G4dou << 94 //Destructor
65 const G4dou << 95 G4BogackiShampine23::~G4BogackiShampine23()
66 G4double yO <<
67 G4double* d <<
68 G4double* y <<
69 { 96 {
>> 97 delete[] ak2;
>> 98 delete[] ak3;
>> 99 delete[] ak4;
>> 100
>> 101 delete[] yTemp;
>> 102 delete[] yIn;
>> 103
>> 104 delete[] fLastInitialVector;
>> 105 delete[] fLastFinalVector;
>> 106 delete[] fLastDyDx;
>> 107 delete[] fMidVector;
>> 108 delete[] fMidError;
70 109
71 G4double yTemp[G4FieldTrack::ncompSVEC]; << 110 delete fAuxStepper;
72 for(G4int i = GetNumberOfVariables(); i < Ge << 111 }
73 { <<
74 yOutput[i] = yTemp[i] = yInput[i]; <<
75 } <<
76 112
77 G4double ak2[G4FieldTrack::ncompSVEC], << 113 //******************************************************************************
78 ak3[G4FieldTrack::ncompSVEC]; << 114 //
>> 115 // Given values for n = 4 variables yIn[0,...,n-1]
>> 116 // known at x, use the 3rd order Bogacki Shampine method
>> 117 // to advance the solution over an interval Step
>> 118 // and return the incremented variables as yOut[0,...,n-1]. Also
>> 119 // return an estimate of the local truncation error yErr[] using the
>> 120 // embedded 2nd order method. The user supplies routine
>> 121 // RightHandSide(y,dydx), which returns derivatives dydx for y .
>> 122
>> 123
>> 124 //******************************************************************************
>> 125
>> 126
>> 127 void
>> 128 G4BogackiShampine23::Stepper( const G4double yInput[],
>> 129 const G4double DyDx[],
>> 130 G4double Step,
>> 131 G4double yOut[],
>> 132 G4double yErr[])
>> 133 {
>> 134
>> 135 G4int i;
79 136
80 const G4double b21 = 0.5 , << 137 const G4double b21 = 0.5 ,
81 b31 = 0., b32 = 3.0 / 4.0, << 138 b31 = 0. , b32 = 3.0/4.0 ,
82 b41 = 2.0 / 9.0, b42 = 1.0 / << 139 b41 = 2.0/9.0, b42 = 1.0/3.0 , b43 = 4.0/9.0;
83 << 140
84 const G4double dc1 = b41 - 7.0 / 24.0, dc2 << 141
85 dc3 = b43 - 1.0 / 3.0, dc4 << 142 const G4double dc1 = b41 - 7.0/24.0 , dc2 = b42 - 1.0/4.0 ,
86 << 143 dc3 = b43 - 1.0/3.0 , dc4 = - 0.125 ;
87 // RightHandSide(yInput, dydx); << 144
88 for(G4int i = 0; i < GetNumberOfVariables(); <<
89 { <<
90 yTemp[i] = yInput[i] + b21 * hstep * dydx[ <<
91 } <<
92 145
93 RightHandSide(yTemp, ak2); <<
94 for(G4int i = 0; i < GetNumberOfVariables(); <<
95 { <<
96 yTemp[i] = yInput[i] + hstep * (b31 * dydx <<
97 } <<
98 146
99 RightHandSide(yTemp, ak3); << 147 // Initialise time to t0, needed when it is not updated by the integration.
100 for(G4int i = 0; i < GetNumberOfVariables(); << 148 // [ Note: Only for time dependent fields (usually electric)
101 { << 149 // is it neccessary to integrate the time.]
102 yOutput[i] = yInput[i] + hstep * (b41*dydx << 150 yOut[7] = yTemp[7] = yIn[7];
103 } << 151
>> 152 const G4int numberOfVariables= this->GetNumberOfVariables(); // The number of variables to be integrated over
>> 153
>> 154 // Saving yInput because yInput and yOut can be aliases for same array
>> 155
>> 156 for(i=0;i<numberOfVariables;i++)
>> 157 {
>> 158 yIn[i]=yInput[i];
>> 159 }
>> 160 // RightHandSide(yIn, dydx) ; // 1st Step --Not doing, getting passed
104 161
105 if ((dydxOutput != nullptr) && (yError != nu << 162 for(i=0;i<numberOfVariables;i++)
106 { <<
107 RightHandSide(yOutput, dydxOutput); <<
108 for(G4int i = 0; i < GetNumberOfVariables( <<
109 { 163 {
110 yError[i] = hstep * (dc1 * dydx[i] + dc2 << 164 yTemp[i] = yIn[i] + b21*Step*DyDx[i] ;
111 dc3 * ak3[i] + dc4 <<
112 } 165 }
113 } << 166 RightHandSide(yTemp, ak2) ; // 2nd Step
>> 167
>> 168 for(i=0;i<numberOfVariables;i++)
>> 169 {
>> 170 yTemp[i] = yIn[i] + Step*(b31*DyDx[i] + b32*ak2[i]) ;
>> 171 }
>> 172 RightHandSide(yTemp, ak3) ; // 3rd Step
>> 173
>> 174 for(i=0;i<numberOfVariables;i++)
>> 175 {
>> 176 yOut[i] = yIn[i] + Step*(b41*DyDx[i] + b42*ak2[i] + b43*ak3[i]) ;
>> 177 }
>> 178 RightHandSide(yOut, ak4) ; // 4th Step
>> 179
>> 180 for(i=0;i<numberOfVariables;i++)
>> 181 {
>> 182 // yOut[i] = yIn[i] + Step*(c1*DyDx[i]+ c2*ak2[i] + c3*ak3[i] + c4*ak4[i]);
>> 183
>> 184 yErr[i] = Step*(dc1*DyDx[i] + dc2*ak2[i] + dc3*ak3[i] +
>> 185 dc4*ak4[i] ) ;
>> 186
>> 187
>> 188 // Store Input and Final values, for possible use in calculating chord
>> 189 fLastInitialVector[i] = yIn[i] ;
>> 190 fLastFinalVector[i] = yOut[i];
>> 191 fLastDyDx[i] = DyDx[i];
>> 192 }
>> 193 // NormaliseTangentVector( yOut ); // Not wanted
>> 194
>> 195 fLastStepLength =Step;
>> 196
>> 197 return ;
114 } 198 }
115 199
116 void G4BogackiShampine23::Stepper(const G4doub << 200 G4double G4BogackiShampine23::DistChord() const
117 const G4doub <<
118 G4double hst <<
119 G4double yOu <<
120 G4double yEr <<
121 { 201 {
122 copy(fyIn, yInput); << 202 G4double distLine, distChord;
123 copy(fdydx, dydx); << 203 G4ThreeVector initialPoint, finalPoint, midPoint;
124 fhstep = hstep; <<
125 <<
126 makeStep(fyIn, fdydx, fhstep, fyOut, fdydxOu <<
127 204
128 copy(yOutput, fyOut); << 205 // Store last initial and final points (they will be overwritten in self-Stepper call!)
129 } << 206 initialPoint = G4ThreeVector( fLastInitialVector[0],
>> 207 fLastInitialVector[1], fLastInitialVector[2]);
>> 208 finalPoint = G4ThreeVector( fLastFinalVector[0],
>> 209 fLastFinalVector[1], fLastFinalVector[2]);
130 210
131 void G4BogackiShampine23::Stepper(const G4doub << 211 // Do half a step using StepNoErr
132 const G4doub <<
133 G4double hst <<
134 G4double yOu <<
135 G4double yEr <<
136 G4double dyd <<
137 { <<
138 copy(fyIn, yInput); <<
139 copy(fdydx, dydx); <<
140 fhstep = hstep; <<
141 212
142 makeStep(fyIn, fdydx, fhstep, fyOut, fdydxOu << 213 fAuxStepper->Stepper( fLastInitialVector, fLastDyDx, 0.5 * fLastStepLength,
>> 214 fMidVector, fMidError );
143 215
144 copy(yOutput, fyOut); << 216 midPoint = G4ThreeVector( fMidVector[0], fMidVector[1], fMidVector[2]);
145 copy(dydxOutput, fdydxOut); <<
146 } <<
147 217
148 G4double G4BogackiShampine23::DistChord() cons << 218 // Use stored values of Initial and Endpoint + new Midpoint to evaluate
149 { << 219 // distance of Chord
150 G4double yMid[G4FieldTrack::ncompSVEC]; <<
151 makeStep(fyIn, fdydx, fhstep / 2., yMid); <<
152 220
153 const G4ThreeVector begin = makeVector(fyIn, <<
154 const G4ThreeVector mid = makeVector(yMid, V <<
155 const G4ThreeVector end = makeVector(fyOut, <<
156 221
157 return G4LineSection::Distline(mid, begin, e << 222 if (initialPoint != finalPoint)
>> 223 {
>> 224 distLine = G4LineSection::Distline( midPoint, initialPoint, finalPoint );
>> 225 distChord = distLine;
>> 226 }
>> 227 else
>> 228 {
>> 229 distChord = (midPoint-initialPoint).mag();
>> 230 }
>> 231 return distChord;
158 } 232 }
>> 233
>> 234 //------Verified-------
159 235