include/Geant4/G4FSALIntegrationDriver.icc

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 // G4FSALIntegrationDriver inline implementation
0027 //
0028 // Author: Dmitry Sorokin (CERN, Google Summer of Code 2017), 20.10.2017
0029 // --------------------------------------------------------------------
0030
0031 #include "G4FieldUtils.hh"
0032
0033 template <class T>
0034 G4FSALIntegrationDriver<T>::
0035 G4FSALIntegrationDriver ( G4double hminimum, T* pStepper,
0036                           G4int numComponents, G4int statisticsVerbose )
0037  : Base(pStepper),
0038    fMinimumStep(hminimum),
0039    fVerboseLevel(statisticsVerbose)
0040 {
0041   if (numComponents != Base::GetStepper()->GetNumberOfVariables())
0042   {
0043     std::ostringstream message;
0044     message << "Driver's number of integrated components "
0045             << numComponents
0046             << " != Stepper's number of components "
0047             << pStepper->GetNumberOfVariables();
0048     G4Exception("G4FSALIntegrationDriver","GeomField0002",
0049                 FatalException, message);
0050   }
0051 }
0052
0053 template <class T>
0054 G4FSALIntegrationDriver<T>::~G4FSALIntegrationDriver()
0055 {
0056 #ifdef G4VERBOSE
0057   if( fVerboseLevel > 0 )
0058   {
0059     G4cout << "G4FSALIntegration Driver Stats: "
0060            << "#QuickAdvance " << fNoQuickAvanceCalls
0061            << " - #AccurateAdvance " << fNoAccurateAdvanceCalls << G4endl
0062            << "#good steps " << fNoAccurateAdvanceGoodSteps << " "
0063            << "#bad steps " << fNoAccurateAdvanceBadSteps << G4endl;
0064   }
0065 #endif
0066 }
0067
0068 // Runge-Kutta driver with adaptive stepsize control. Integrate starting
0069 // values at y_current over hstep x2 with accuracy eps.
0070 // On output ystart is replaced by values at the end of the integration
0071 // interval. RightHandSide is the right-hand side of ODE system.
0072 // The source is similar to odeint routine from NRC p.721-722 .
0073 //
0074 template <class T>
0075 G4bool G4FSALIntegrationDriver<T>::
0076 AccurateAdvance( G4FieldTrack& track, G4double hstep,
0077                  G4double eps, G4double hinitial )
0078 {
0079   ++fNoAccurateAdvanceCalls;
0080
0081   if (hstep < GetMinimumStep())
0082   {
0083     G4double dchord_step = 0.0, dyerr = 0.0;
0084     G4double dydx[G4FieldTrack::ncompSVEC];
0085     Base::GetDerivatives(track, dydx);
0086     return QuickAdvance(track, dydx, hstep, dchord_step, dyerr);
0087   }
0088
0089   G4bool succeeded = false;
0090
0091   G4double hnext, hdid;
0092
0093   G4double y[G4FieldTrack::ncompSVEC], dydx[G4FieldTrack::ncompSVEC];
0094
0095   track.DumpToArray(y);
0096
0097   // hstep somtimes is too small. No need to add large curveLength.
0098   G4double curveLength = 0.0;
0099   G4double endCurveLength = hstep;
0100
0101   G4double h = hstep;
0102   if (hinitial > CLHEP::perMillion * hstep && hinitial < hstep)
0103   {
0104     h = hinitial;
0105   }
0106
0107   Base::GetStepper()->RightHandSide(y, dydx);
0108
0109   for (G4int iter = 0; iter < Base::GetMaxNoSteps(); ++iter)
0110   {
0111     const G4ThreeVector StartPos =
0112           field_utils::makeVector(y, field_utils::Value3D::Position);
0113
0114     OneGoodStep(y, dydx, curveLength, h, eps, hdid, hnext);
0115
0116     const G4ThreeVector EndPos =
0117           field_utils::makeVector(y, field_utils::Value3D::Position);
0118
0119     CheckStep(EndPos, StartPos, hdid);
0120
0121     G4double restCurveLength = endCurveLength - curveLength;
0122     if (restCurveLength < GetSmallestFraction() * hstep)
0123     {
0124       succeeded = true;
0125       break;
0126     }
0127     h = std::min(hnext, restCurveLength);
0128   }
0129
0130
0131   if (succeeded)
0132   {
0133     track.LoadFromArray(y, Base::GetStepper()->GetNumberOfVariables());
0134     track.SetCurveLength(track.GetCurveLength() + curveLength);
0135   }
0136
0137   return succeeded;
0138 }
0139
0140 // Driver for one Runge-Kutta Step with monitoring of local truncation error
0141 // to ensure accuracy and adjust stepsize. Input are dependent variable
0142 // array y[0,...,5] and its derivative dydx[0,...,5] at the
0143 // starting value of the independent variable x . Also input are stepsize
0144 // to be attempted htry, and the required accuracy eps. On output y and x
0145 // are replaced by their new values, hdid is the stepsize that was actually
0146 // accomplished, and hnext is the estimated next stepsize.
0147 // This is similar to the function rkqs from the book:
0148 // Numerical Recipes in C: The Art of Scientific Computing (NRC), Second
0149 // Edition, by William H. Press, Saul A. Teukolsky, William T.
0150 // Vetterling, and Brian P. Flannery (Cambridge University Press 1992),
0151 // 16.2 Adaptive StepSize Control for Runge-Kutta, p. 719
0152 //
0153 template <class T>
0154 void G4FSALIntegrationDriver<T>::
0155 OneGoodStep(G4double y[],
0156             G4double dydx[],
0157             G4double& curveLength,   // InOut
0158             G4double htry,
0159             G4double eps_rel_max,
0160             G4double& hdid,         // Out
0161             G4double& hnext)        // Out
0162 {
0163   G4double error2 = DBL_MAX;
0164
0165   G4double yError[G4FieldTrack::ncompSVEC],
0166            yOut[G4FieldTrack::ncompSVEC],
0167            dydxOut[G4FieldTrack::ncompSVEC];
0168
0169   // Set stepsize to the initial trial value
0170   G4double hstep = htry;
0171
0172   const G4int max_trials = 100;
0173
0174   for (G4int iter = 0; iter < max_trials; ++iter)
0175   {
0176     Base::GetStepper()->Stepper(y, dydx, hstep, yOut, yError, dydxOut);
0177     error2 = field_utils::relativeError2(y, yError, hstep, eps_rel_max);
0178
0179     // Step succeeded.
0180     if (error2 <= 1) {  break; }
0181
0182     hstep = Base::ShrinkStepSize2(hstep, error2);
0183   }
0184
0185   hnext = Base::GrowStepSize2(hstep, error2);
0186   curveLength += (hdid = hstep);
0187
0188   for(G4int k = 0; k < Base::GetStepper()->GetNumberOfVariables(); ++k)
0189   {
0190     y[k] = yOut[k];
0191     dydx[k] = dydxOut[k];
0192   }
0193 }
0194
0195 template <class T>
0196 G4bool G4FSALIntegrationDriver<T>::
0197 QuickAdvance( G4FieldTrack& fieldTrack, const G4double dydxIn[],
0198               G4double hstep,
0199               G4double& dchord_step, G4double& dyerr )
0200 {
0201   ++fNoQuickAvanceCalls;
0202
0203   if (hstep == 0)
0204   {
0205     std::ostringstream message;
0206     message << "Proposed step is zero; hstep = " << hstep << " !";
0207     G4Exception("G4FSALIntegrationDriver ::QuickAdvance()",
0208                 "GeomField1001", JustWarning, message);
0209     return true;
0210   }
0211   if (hstep < 0)
0212   {
0213     std::ostringstream message;
0214     message << "Invalid run condition." << G4endl
0215             << "Proposed step is negative; hstep = "
0216             << hstep << "." << G4endl
0217             << "Requested step cannot be negative! Aborting event.";
0218     G4Exception("G4FSALIntegrationDriver ::QuickAdvance()",
0219                 "GeomField0003", EventMustBeAborted, message);
0220     return false;
0221   }
0222
0223   G4double yError[G4FieldTrack::ncompSVEC],
0224            yIn[G4FieldTrack::ncompSVEC],
0225            yOut[G4FieldTrack::ncompSVEC],
0226            dydxOut[G4FieldTrack::ncompSVEC];
0227
0228   fieldTrack.DumpToArray(yIn);
0229
0230   Base::GetStepper()->Stepper(yIn, dydxIn, hstep, yOut, yError, dydxOut);
0231   dchord_step = Base::GetStepper()->DistChord();
0232
0233   fieldTrack.LoadFromArray(yOut, Base::GetStepper()->GetNumberOfVariables());
0234   fieldTrack.SetCurveLength(fieldTrack.GetCurveLength() + hstep);
0235
0236   dyerr = field_utils::absoluteError(yOut, yError, hstep);
0237
0238   return true;
0239 }
0240
0241 template <class T>
0242 void G4FSALIntegrationDriver<T>::SetSmallestFraction(G4double newFraction)
0243 {
0244   if( newFraction > 1.e-16 && newFraction < 1e-8 )
0245   {
0246     fSmallestFraction = newFraction;
0247   }
0248   else
0249   {
0250     std::ostringstream message;
0251     message << "Smallest Fraction not changed. " << G4endl
0252             << "  Proposed value was " << newFraction << G4endl
0253             << "  Value must be between 1.e-8 and 1.e-16";
0254     G4Exception("G4FSALIntegrationDriver::SetSmallestFraction()",
0255                 "GeomField1001", JustWarning, message);
0256   }
0257 }
0258
0259 template <class T>
0260 void G4FSALIntegrationDriver<T>::CheckStep(
0261     const G4ThreeVector& posIn, const G4ThreeVector& posOut, G4double hdid)
0262 {
0263   const G4double endPointDist = (posOut - posIn).mag();
0264   if (endPointDist >= hdid * (1. + CLHEP::perMillion))
0265   {
0266     ++fNoAccurateAdvanceBadSteps;
0267 #ifdef G4DEBUG_FIELD
0268     // Issue a warning only for gross differences -
0269     // we understand how small difference occur.
0270     if (endPointDist >= hdid * (1. + perThousand))
0271     {
0272       G4Exception("G4FSALIntegrationDriver::CheckStep()",
0273                   "GeomField1002", JustWarning, "endPointDist >= hdid!");
0274     }
0275 #endif
0276   }
0277   else
0278   {
0279     ++fNoAccurateAdvanceGoodSteps;
0280   }
0281 }
0282
0283 template <class T>
0284 inline G4double G4FSALIntegrationDriver<T>::GetMinimumStep() const
0285 {
0286   return fMinimumStep;
0287 }
0288
0289 template <class T>
0290 void G4FSALIntegrationDriver<T>::SetMinimumStep(G4double minimumStepLength)
0291 {
0292   fMinimumStep = minimumStepLength;
0293 }
0294
0295 template <class T>
0296 G4int G4FSALIntegrationDriver<T>::GetVerboseLevel() const
0297 {
0298   return fVerboseLevel;
0299 }
0300
0301 template <class T>
0302 void G4FSALIntegrationDriver<T>::SetVerboseLevel(G4int newLevel)
0303 {
0304   fVerboseLevel = newLevel;
0305 }
0306
0307 template <class T>
0308 G4double G4FSALIntegrationDriver<T>::GetSmallestFraction() const
0309 {
0310   return fSmallestFraction;
0311 }
0312
0313 template <class T>
0314 void G4FSALIntegrationDriver<T>::IncrementQuickAdvanceCalls()
0315 {
0316   ++fNoQuickAvanceCalls;
0317 }
0318
0319 template <class T>
0320 G4double
0321 G4FSALIntegrationDriver<T>::AdvanceChordLimited(G4FieldTrack& track,
0322                                              G4double hstep,
0323                                              G4double eps,
0324                                              G4double chordDistance)
0325 {
0326   return ChordFinderDelegate::AdvanceChordLimitedImpl(track, hstep,
0327                                                       eps, chordDistance);
0328 }
0329
0330 template <class T>
0331 void G4FSALIntegrationDriver<T>::OnStartTracking()
0332 {
0333   ChordFinderDelegate::ResetStepEstimate();
0334 }
0335
0336 template <class T>
0337 void G4FSALIntegrationDriver<T>::OnComputeStep(const G4FieldTrack*)
0338 {
0339 }
0340
0341 template <class T>
0342 G4bool G4FSALIntegrationDriver<T>::DoesReIntegrate() const
0343 {
0344   return true;
0345 }
0346
0347 template <class T>
0348 void G4FSALIntegrationDriver<T>::StreamInfo( std::ostream& os ) const
0349 {
0350   // Write out the parameters / state of the driver
0351
0352   os << "State of G4IntegrationDriver: " << std::endl;
0353   os << "--Base state (G4RKIntegrationDriver): " << std::endl;
0354   Base::StreamInfo( os );
0355   os << "--Own  state (G4IntegrationDriver<>): " << std::endl;
0356   os << "    fMinimumStep =      " << fMinimumStep  << std::endl;
0357   os << "    Smallest Fraction = " << fSmallestFraction << std::endl;
0358
0359   os << "    verbose level     = " << fVerboseLevel << std::endl;
0360   os << "    Reintegrates      = " << DoesReIntegrate() << std::endl;
0361   os << "--Chord Finder Delegate state: " << std::endl;
0362   ChordFinderDelegate::StreamDelegateInfo( os );
0363 }