include/Geant4/G4IntegrationDriver.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 // G4IntegrationDriver inline implementation
0027 //
0028 // Author: Dmitry Sorokin (CERN, Google Summer of Code 2017), 20.10.2017
0029 // Supervision: John Apostolakis (CERN)
0030 // --------------------------------------------------------------------
0031
0032 #include "G4FieldUtils.hh"
0033
0034 #include <CLHEP/Units/SystemOfUnits.h>
0035
0036 template <class T>
0037 G4IntegrationDriver<T>::
0038 G4IntegrationDriver ( G4double hminimum, T* pStepper,
0039                       G4int numComponents, G4int statisticsVerbose )
0040     : G4RKIntegrationDriver<T>(pStepper),
0041       fMinimumStep(hminimum),
0042       fVerboseLevel(statisticsVerbose)
0043 {
0044   if (numComponents != Base::GetStepper()->GetNumberOfVariables())
0045   {
0046     std::ostringstream message;
0047     message << "Driver's number of integrated components "
0048             << numComponents
0049             << " != Stepper's number of components "
0050             << pStepper->GetNumberOfVariables();
0051     G4Exception("G4IntegrationDriver","GeomField0002",
0052                 FatalException, message);
0053   }
0054 }
0055
0056 template <class T>
0057 G4IntegrationDriver<T>::~G4IntegrationDriver()
0058 {
0059 #ifdef G4VERBOSE
0060   if (fVerboseLevel > 0)
0061   {
0062     G4cout << "G4Integration Driver Stats: "
0063            << "#QuickAdvance " << fNoQuickAvanceCalls
0064            << " - #AccurateAdvance " << fNoAccurateAdvanceCalls << " "
0065            << "#good steps " << fNoAccurateAdvanceGoodSteps << " "
0066            << "#bad steps " << fNoAccurateAdvanceBadSteps << G4endl;
0067   }
0068 #endif
0069 }
0070
0071 template <class T>
0072 G4double G4IntegrationDriver<T>::AdvanceChordLimited(G4FieldTrack& track,
0073                                                      G4double stepMax,
0074                                                      G4double epsStep,
0075                                                      G4double chordDistance)
0076 {
0077   return ChordFinderDelegate::AdvanceChordLimitedImpl(track, stepMax, epsStep,
0078                                                       chordDistance);
0079 }
0080
0081 template <class T>
0082 void G4IntegrationDriver<T>::OnStartTracking()
0083 {
0084   ChordFinderDelegate::ResetStepEstimate();
0085 }
0086
0087 template <class T>
0088 void G4IntegrationDriver<T>::OnComputeStep(const G4FieldTrack*)
0089 {
0090 }
0091
0092 template <class T>
0093 G4bool G4IntegrationDriver<T>::DoesReIntegrate() const
0094 {
0095   return true;
0096 }
0097
0098 // Runge-Kutta driver with adaptive stepsize control. Integrate starting
0099 // values at y_current over hstep x2 with accuracy eps.
0100 // On output ystart is replaced by values at the end of the integration
0101 // interval. RightHandSide is the right-hand side of ODE system.
0102 // The source is similar to odeint routine from NRC p.721-722 .
0103 //
0104 template <class T>
0105 G4bool G4IntegrationDriver<T>::
0106 AccurateAdvance(G4FieldTrack& track, G4double hstep,
0107                 G4double eps, G4double hinitial)
0108 {
0109   ++fNoAccurateAdvanceCalls;
0110
0111   if (hstep == 0.0)
0112   {
0113     std::ostringstream message;
0114     message << "Proposed step is zero; hstep = " << hstep << " !";
0115     G4Exception("G4IntegrationDriver::AccurateAdvance()",
0116                 "GeomField1001", JustWarning, message);
0117     return true;
0118   }
0119
0120   if (hstep < 0)
0121   {
0122     std::ostringstream message;
0123     message << "Invalid run condition." << G4endl
0124             << "Proposed step is negative; hstep = " << hstep << "."
0125             << G4endl
0126             << "Requested step cannot be negative! Aborting event.";
0127     G4Exception("G4IntegrationDriver::AccurateAdvance()",
0128                 "GeomField0003", EventMustBeAborted, message);
0129     return false;
0130   }
0131
0132   G4double hnext, hdid;
0133
0134   G4double dydx[G4FieldTrack::ncompSVEC];
0135   G4bool succeeded = true;
0136
0137   G4double y[G4FieldTrack::ncompSVEC];
0138   track.DumpToArray(y);
0139
0140   const G4double startCurveLength = track.GetCurveLength();
0141   const G4double endCurveLength = startCurveLength + hstep;
0142   const G4double hThreshold =
0143         std::min(eps * hstep, fSmallestFraction * startCurveLength);
0144
0145   G4double h = hstep;
0146   if (hinitial > CLHEP::perMillion * hstep && hinitial < hstep)
0147   {
0148     h = hinitial;
0149   }
0150
0151   G4double curveLength = startCurveLength;
0152
0153   for (G4int nstp = 0; nstp < Base::GetMaxNoSteps(); ++nstp)
0154   {
0155     const G4ThreeVector StartPos =
0156           field_utils::makeVector(y, field_utils::Value3D::Position);
0157
0158     Base::GetStepper()->RightHandSide(y, dydx);
0159
0160     if (h > GetMinimumStep())
0161     {
0162       OneGoodStep(y, dydx, curveLength, h, eps, hdid, hnext);
0163     }
0164     else
0165     {
0166       G4FieldTrack yFldTrk('0');
0167       G4double dchord_step, dyerr, dyerr_len;
0168       yFldTrk.LoadFromArray(y, Base::GetStepper()->GetNumberOfVariables());
0169       yFldTrk.SetCurveLength(curveLength);
0170
0171       QuickAdvance(yFldTrk, dydx, h, dchord_step, dyerr_len);
0172
0173       yFldTrk.DumpToArray(y);
0174
0175       if (h == 0.0)
0176       {
0177         G4Exception("G4IntegrationDriver::AccurateAdvance()",
0178                     "GeomField0003", FatalException,
0179                     "Integration Step became Zero!");
0180       }
0181       dyerr = dyerr_len / h;
0182       hdid = h;
0183       curveLength += hdid;
0184       hnext = Base::ComputeNewStepSize(dyerr / eps, h);
0185     }
0186
0187     const G4ThreeVector EndPos =
0188           field_utils::makeVector(y, field_utils::Value3D::Position);
0189
0190     CheckStep(EndPos, StartPos, hdid);
0191
0192     //  Avoid numerous small last steps
0193     if (h < hThreshold || curveLength >= endCurveLength)
0194     {
0195       break;
0196     }
0197
0198     h = std::max(hnext, GetMinimumStep());
0199     if (curveLength + h > endCurveLength)
0200     {
0201       h = endCurveLength - curveLength;
0202     }
0203   }
0204   // Have we reached the end ?
0205   // --> a better test might be x-endCurveLength > an_epsilon
0206   succeeded = (curveLength >= endCurveLength);
0207     // If it was a "forced" last step
0208
0209   track.LoadFromArray(y, Base::GetStepper()->GetNumberOfVariables());
0210   track.SetCurveLength(curveLength);
0211
0212   return succeeded;
0213 }
0214
0215 // Driver for one Runge-Kutta Step with monitoring of local truncation error
0216 // to ensure accuracy and adjust stepsize. Input are dependent variable
0217 // array y[0,...,5] and its derivative dydx[0,...,5] at the
0218 // starting value of the independent variable x . Also input are stepsize
0219 // to be attempted htry, and the required accuracy eps. On output y and x
0220 // are replaced by their new values, hdid is the stepsize that was actually
0221 // accomplished, and hnext is the estimated next stepsize.
0222 // This is similar to the function rkqs from the book:
0223 // Numerical Recipes in C: The Art of Scientific Computing (NRC), Second
0224 // Edition, by William H. Press, Saul A. Teukolsky, William T.
0225 // Vetterling, and Brian P. Flannery (Cambridge University Press 1992),
0226 // 16.2 Adaptive StepSize Control for Runge-Kutta, p. 719
0227 //
0228 template <class T>
0229 void G4IntegrationDriver<T>::OneGoodStep(G4double y[],           // InOut
0230                                          const G4double dydx[],
0231                                          G4double& curveLength,  // InOut
0232                                          G4double htry,
0233                                          G4double eps_rel_max,
0234                                          G4double& hdid,      // Out
0235                                          G4double& hnext)    // Out
0236
0237 {
0238   G4double error2 = DBL_MAX;
0239
0240   G4double yerr[G4FieldTrack::ncompSVEC], ytemp[G4FieldTrack::ncompSVEC];
0241
0242   G4double h = htry;
0243
0244   const G4int max_trials = 100;
0245
0246   for (G4int iter = 0; iter < max_trials; ++iter)
0247   {
0248     Base::GetStepper()->Stepper(y, dydx, h, ytemp, yerr);
0249     error2 = field_utils::relativeError2(y, yerr, std::max(h, fMinimumStep),
0250                                              eps_rel_max);
0251     if (error2 <= 1.0)
0252     {
0253       break;
0254     }
0255
0256     h = Base::ShrinkStepSize2(h, error2);
0257
0258     G4double xnew = curveLength + h;
0259     if(xnew == curveLength)
0260     {
0261       std::ostringstream message;
0262       message << "Stepsize underflow in Stepper !" << G4endl
0263               << "- Step's start x=" << curveLength
0264               << " and end x= " << xnew
0265               << " are equal !! " << G4endl
0266               << "  Due to step-size= " << h
0267               << ". Note that input step was " << htry;
0268       G4Exception("G4IntegrationDriver::OneGoodStep()",
0269                   "GeomField1001", JustWarning, message);
0270       break;
0271     }
0272   }
0273
0274   hnext = Base::GrowStepSize2(h, error2);
0275   curveLength += (hdid = h);
0276
0277   field_utils::copy(y, ytemp, Base::GetStepper()->GetNumberOfVariables());
0278 }
0279
0280 template <class T>
0281 G4bool G4IntegrationDriver<T>::QuickAdvance(G4FieldTrack& track,    // INOUT
0282                                       const G4double dydx[],
0283                                             G4double hstep,
0284                                             G4double& dchord_step,
0285                                             G4double& dyerr)
0286 {
0287   ++fNoQuickAvanceCalls;
0288
0289   G4double yIn[G4FieldTrack::ncompSVEC],
0290            yOut[G4FieldTrack::ncompSVEC],
0291            yError[G4FieldTrack::ncompSVEC];
0292
0293   track.DumpToArray(yIn);
0294
0295   Base::GetStepper()->Stepper(yIn, dydx, hstep, yOut, yError);
0296
0297   dchord_step = Base::GetStepper()->DistChord();
0298   dyerr = field_utils::absoluteError(yOut, yError, hstep);
0299   track.LoadFromArray(yOut, Base::GetStepper()->GetNumberOfVariables());
0300   track.SetCurveLength(track.GetCurveLength() + hstep);
0301
0302   return true;
0303 }
0304
0305 template <class T>
0306 void G4IntegrationDriver<T>::CheckStep(const G4ThreeVector& posIn,
0307                                        const G4ThreeVector& posOut,
0308                                              G4double hdid)
0309 {
0310   const G4double endPointDist = (posOut - posIn).mag();
0311   if (endPointDist >= hdid * (1. + CLHEP::perMillion))
0312   {
0313     ++fNoAccurateAdvanceBadSteps;
0314 #ifdef G4DEBUG_FIELD
0315     // Issue a warning only for gross differences -
0316     // we understand how small difference occur.
0317     if (endPointDist >= hdid * (1. + perThousand))
0318     {
0319        G4Exception("G4IntegrationDriver::CheckStep()",
0320                    "GeomField1002", JustWarning,
0321                    "endPointDist >= hdid!");
0322     }
0323 #endif
0324   }
0325   else
0326   {
0327     ++fNoAccurateAdvanceGoodSteps;
0328   }
0329 }
0330
0331 template <class T>
0332 inline G4double G4IntegrationDriver<T>::GetMinimumStep() const
0333 {
0334   return fMinimumStep;
0335 }
0336
0337 template <class T>
0338 void G4IntegrationDriver<T>::SetMinimumStep(G4double minimumStepLength)
0339 {
0340   fMinimumStep = minimumStepLength;
0341 }
0342
0343 template <class T>
0344 G4int G4IntegrationDriver<T>::GetVerboseLevel() const
0345 {
0346   return fVerboseLevel;
0347 }
0348
0349 template <class T>
0350 void G4IntegrationDriver<T>::SetVerboseLevel(G4int newLevel)
0351 {
0352   fVerboseLevel = newLevel;
0353 }
0354
0355 template <class T>
0356 G4double G4IntegrationDriver<T>::GetSmallestFraction() const
0357 {
0358   return fSmallestFraction;
0359 }
0360
0361 template <class T>
0362 void G4IntegrationDriver<T>::SetSmallestFraction(G4double newFraction)
0363 {
0364   if (newFraction > 1.e-16 && newFraction < 1e-8)
0365   {
0366     fSmallestFraction = newFraction;
0367   }
0368   else
0369   {
0370     std::ostringstream message;
0371     message << "Smallest Fraction not changed. " << G4endl
0372             << "  Proposed value was " << newFraction << G4endl
0373             << "  Value must be between 1.e-8 and 1.e-16";
0374     G4Exception("G4IntegrationDriver::SetSmallestFraction()",
0375                 "GeomField1001", JustWarning, message);
0376   }
0377 }
0378
0379 template <class T>
0380 void G4IntegrationDriver<T>::IncrementQuickAdvanceCalls()
0381 {
0382   ++fNoQuickAvanceCalls;
0383 }
0384
0385 template <class T>
0386 void G4IntegrationDriver<T>::StreamInfo( std::ostream& os ) const
0387 {
0388   // Write out the parameters / state of the driver
0389
0390   os << "State of G4IntegrationDriver: " << std::endl;
0391   os << "--Base state (G4RKIntegrationDriver): " << std::endl;
0392   Base::StreamInfo( os );
0393   os << "--Own  state (G4IntegrationDriver<>): " << std::endl;
0394   os << "    fMinimumStep =      " << fMinimumStep  << std::endl;
0395   os << "    Smallest Fraction = " << fSmallestFraction << std::endl;
0396
0397   os << "    verbose level     = " << fVerboseLevel << std::endl;
0398   os << "    Reintegrates      = " << DoesReIntegrate() << std::endl;
0399   os << "--Chord Finder Delegate state: " << std::endl;
0400   ChordFinderDelegate::StreamDelegateInfo( os );
0401 }