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