| Geant4 Cross Reference |
1 // ********************************************************************
2 // * License and Disclaimer *
3 // * *
4 // * The Geant4 software is copyright of the Copyright Holders of *
5 // * the Geant4 Collaboration. It is provided under the terms and *
6 // * conditions of the Geant4 Software License, included in the file *
7 // * LICENSE and available at http://cern.ch/geant4/license . These *
8 // * include a list of copyright holders. *
9 // * *
10 // * Neither the authors of this software system, nor their employing *
11 // * institutes,nor the agencies providing financial support for this *
12 // * work make any representation or warranty, express or implied, *
13 // * regarding this software system or assume any liability for its *
14 // * use. Please see the license in the file LICENSE and URL above *
15 // * for the full disclaimer and the limitation of liability. *
16 // * *
17 // * This code implementation is the result of the scientific and *
18 // * technical work of the GEANT4 collaboration. *
19 // * By using, copying, modifying or distributing the software (or *
20 // * any work based on the software) you agree to acknowledge its *
21 // * use in resulting scientific publications, and indicate your *
22 // * acceptance of all terms of the Geant4 Software license. *
23 // ********************************************************************
24 //
25 // G4BorisDriver
26 //
27 // Class description:
28 //
29 // G4BorisDriver is a driver class using the second order Boris
30 // method to integrate the equation of motion.
31 //
32 //
33 // Author: Divyansh Tiwari, Google Summer of Code 2022
34 // Supervision: John Apostolakis,Renee Fatemi, Soon Yung Jun
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::fErrorConstraintShrink = std::pow(
45 fMaxSteppingDecrease / fSafetyFactor, 1. / fPowerShrink);
46
47 const G4double G4BorisDriver::fErrorConstraintGrow = std::pow(
48 fMaxSteppingIncrease / fSafetyFactor, 1. / fPowerGrow);
49
50 // --------------------------------------------------------------------------
51
52 G4BorisDriver::
53 G4BorisDriver( G4double hminimum, G4BorisScheme* Boris,
54 G4int numberOfComponents, G4bool verbosity )
55 : fMinimumStep(hminimum),
56 fVerbosity(verbosity),
57 boris(Boris)
58 // , interval_sequence{2,4}
59 {
60 assert(boris->GetNumberOfVariables() == numberOfComponents);
61
62 if(boris->GetNumberOfVariables() != numberOfComponents)
63 {
64 std::ostringstream msg;
65 msg << "Disagreement in number of variables = "
66 << boris->GetNumberOfVariables()
67 << " vs no of components = " << numberOfComponents;
68 G4Exception("G4BorisDriver Constructor:",
69 "GeomField1001", FatalException, msg);
70 }
71
72 }
73
74 // --------------------------------------------------------------------------
75
76 G4bool G4BorisDriver::AccurateAdvance( G4FieldTrack& track,
77 G4double hstep,
78 G4double epsilon,
79 G4double hinitial )
80 {
81 // Specification: Driver with adaptive stepsize control.
82 // Integrate starting values at y_current over hstep x2 with (relative) accuracy 'eps'.
83 // On output 'track' is replaced by values at the end of the integration interval.
84
85 // Ensure that hstep > 0
86 if(hstep == 0)
87 {
88 std::ostringstream message;
89 message << "Proposed step is zero; hstep = " << hstep << " !";
90 G4Exception("G4BorisDriver::AccurateAdvance()",
91 "GeomField1001", JustWarning, message);
92 return true;
93 }
94 if(hstep < 0)
95 {
96 std::ostringstream message;
97 message << "Invalid run condition." << G4endl
98 << "Proposed step is negative; hstep = " << hstep << G4endl
99 << "Requested step cannot be negative! Aborting event.";
100 G4Exception("G4BorisDriver::AccurateAdvance()",
101 "GeomField0003", EventMustBeAborted, message);
102 return false;
103 }
104
105 if( hinitial == 0.0 ) { hinitial = hstep; }
106 if( hinitial < 0.0 ) { hinitial = std::fabs( hinitial ); }
107 // G4double htrial = std::min( hstep, hinitial );
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_SI;
119 const G4int nvar= GetNumberOfVariables();
120
121 // copy non-integration variables to out array
122 //
123 std::memcpy(yOut + nvar,
124 yCurrent + nvar,
125 sizeof(G4double)*(G4FieldTrack::ncompSVEC-nvar));
126
127 G4double curveLength = track.GetCurveLength(); // starting value
128 const G4double endCurveLength = curveLength + hstep;
129
130 // -- Initial version: Did it in one step -- did not account for errors !!!
131 // G4FieldTrack yFldTrk(track);
132 // yFldTrk.LoadFromArray(yCurrent, G4FieldTrack::ncompSVEC);
133 // yFldTrk.SetCurveLength(curveLength);
134 // G4double dchord_step, dyerr_len;
135 // QuickAdvance(yFldTrk, dydxCurrent, htrial, dchord_step, dyerr_len);
136
137 const G4double hThreshold =
138 std::max(epsilon * hstep, fSmallestFraction * curveLength);
139
140 G4double htry= htrial;
141
142 for (G4int nstp = 0; nstp < fMaxNoSteps; ++nstp)
143 {
144 G4double hdid= 0.0, hnext=0.0;
145
146 OneGoodStep(yCurrent, curveLength, htry, epsilon, restMass, charge, hdid, hnext);
147 //*********
148
149 // Simple check: move (distance of displacement) is smaller than length along curve!
150 const G4ThreeVector StartPos = field_utils::makeVector(yCurrent, field_utils::Value3D::Position);
151 const G4ThreeVector EndPos = field_utils::makeVector(yCurrent, field_utils::Value3D::Position);
152 CheckStep(EndPos, StartPos, hdid);
153
154 // Check 1. for finish and 2. *avoid* numerous small last steps
155 if (curveLength >= endCurveLength || htry < hThreshold)
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::ncompSVEC);
169 track.SetCurveLength(curveLength);
170
171 return true;
172 }
173
174 // --------------------------------------------------------------------------
175
176 void G4BorisDriver::OneGoodStep(G4double y[], // InOut
177 G4double& curveLength, // InOut
178 G4double htry,
179 G4double epsilon_rel,
180 G4double restMass,
181 G4double charge,
182 G4double& hdid, // Out
183 G4double& hnext) // Out
184 {
185 G4double error2 = DBL_MAX;
186 G4double yerr[G4FieldTrack::ncompSVEC], ytemp[G4FieldTrack::ncompSVEC];
187
188 G4double h = htry;
189
190 const G4int max_trials = 100;
191
192 for (G4int iter = 0; iter < max_trials; ++iter)
193 {
194 boris->StepWithErrorEstimate(y, restMass, charge, h, ytemp, yerr);
195
196 error2 = field_utils::relativeError2(y, yerr, std::max(h, fMinimumStep),
197 epsilon_rel);
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 Stepper !" << G4endl
210 << "- Step's start x=" << curveLength
211 << " and end x= " << xnew
212 << " are equal !! " << G4endl
213 << " Due to step-size= " << h
214 << ". Note that input step was " << htry;
215 G4Exception("G4IntegrationDriver::OneGoodStep()",
216 "GeomField1001", JustWarning, message);
217 break;
218 }
219 }
220
221 hnext = GrowStepSize2(h, error2);
222 curveLength += (hdid = h);
223
224 field_utils::copy(y, ytemp, GetNumberOfVariables());
225 }
226
227 // ===========------------------------------------------------------===========
228
229 G4bool G4BorisDriver::
230 QuickAdvance( G4FieldTrack& track, const G4double /*dydx*/[],
231 G4double hstep, G4double& missDist, G4double& dyerr)
232 {
233 const auto nvar = boris->GetNumberOfVariables();
234
235 track.DumpToArray(yIn);
236 const G4double curveLength = track.GetCurveLength();
237
238 // call the boris method for step length hstep
239 G4double restMass = track.GetRestMass();
240 G4double charge = track.GetCharge()*e_SI;
241
242 // boris->DoStep(restMass, charge, yIn, yMid, hstep*0.5);
243 // boris->DoStep(restMass, charge, yMid, yOut, hstep*0.5); // Use mid-point !!
244 boris->StepWithMidAndErrorEstimate(yIn, restMass, charge, hstep,
245 yMid, yOut, yError);
246 // Same, and also return mid-point evaluation
247
248 // How to calculate chord length??
249 const auto mid = field_utils::makeVector(yMid,
250 field_utils::Value3D::Position);
251 const auto in = field_utils::makeVector(yIn,
252 field_utils::Value3D::Position);
253 const auto out = field_utils::makeVector(yOut,
254 field_utils::Value3D::Position);
255
256 missDist = G4LineSection::Distline(mid, in, out);
257
258 dyerr = field_utils::absoluteError(yOut, yError, hstep);
259
260 // copy non-integrated variables to output array
261 //
262 std::memcpy(yOut + nvar, yIn + nvar,
263 sizeof(G4double) * (G4FieldTrack::ncompSVEC - nvar));
264
265 // set new state
266 //
267 track.LoadFromArray(yOut, G4FieldTrack::ncompSVEC);
268 track.SetCurveLength(curveLength + hstep);
269
270 return true;
271 }
272
273 // --------------------------------------------------------------------------------
274
275 void G4BorisDriver::
276 GetDerivatives( const G4FieldTrack& yTrack, G4double dydx[]) const
277 {
278 G4double EBfieldValue[6];
279 GetDerivatives(yTrack, dydx, EBfieldValue);
280 }
281
282 // --------------------------------------------------------------------------------
283
284 void G4BorisDriver::
285 GetDerivatives( const G4FieldTrack& yTrack, G4double dydx[],
286 G4double EBfieldValue[]) const
287 {
288 // G4Exception("G4BorisDriver::GetDerivatives()",
289 // "GeomField0003", FatalException, "This method is not implemented.");
290 G4double ytemp[G4FieldTrack::ncompSVEC];
291 yTrack.DumpToArray(ytemp);
292 GetEquationOfMotion()->EvaluateRhsReturnB(ytemp, dydx, EBfieldValue);
293 }
294
295 // --------------------------------------------------------------------------------
296
297 G4double G4BorisDriver::ShrinkStepSize2(G4double h, G4double error2) const
298 {
299 if (error2 > fErrorConstraintShrink * fErrorConstraintShrink)
300 {
301 return fMaxSteppingDecrease * h;
302 }
303 return fSafetyFactor * h * std::pow(error2, 0.5 * fPowerShrink);
304 }
305
306 // --------------------------------------------------------------------------------
307
308 G4double G4BorisDriver::GrowStepSize2(G4double h, G4double error2) const
309 // Given the square of the relative error,
310 {
311 if (error2 < fErrorConstraintGrow * fErrorConstraintGrow)
312 {
313 return fMaxSteppingIncrease * h;
314 }
315 return fSafetyFactor * h * std::pow(error2, 0.5 * fPowerGrow);
316 }
317
318 // --------------------------------------------------------------------------------
319
320 void G4BorisDriver::SetEquationOfMotion(G4EquationOfMotion* /*equation*/ )
321 {
322 G4Exception("G4BorisDriver::SetEquationOfMotion()", "GeomField0003", FatalException,
323 "This method is not implemented. BorisDriver/Stepper should keep its equation");
324 }
325
326 // --------------------------------------------------------------------------------
327
328 void
329 G4BorisDriver::StreamInfo( std::ostream& os ) const
330 {
331 os << "State of G4BorisDriver: " << std::endl;
332 os << " Method is implemented, but gives no information. " << std::endl;
333 }
334