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 // Created: D.Sorokin, 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
0102     G4double h = hstep;
0103     if (hinitial > CLHEP::perMillion * hstep && hinitial < hstep)
0104     {
0105         h = hinitial;
0106     }
0107
0108     Base::GetStepper()->RightHandSide(y, dydx);
0109
0110     for (G4int iter = 0; iter < Base::GetMaxNoSteps(); ++iter)
0111     {
0112         const G4ThreeVector StartPos =
0113             field_utils::makeVector(y, field_utils::Value3D::Position);
0114
0115         OneGoodStep(y, dydx, curveLength, h, eps, hdid, hnext);
0116
0117         const G4ThreeVector EndPos =
0118             field_utils::makeVector(y, field_utils::Value3D::Position);
0119
0120         CheckStep(EndPos, StartPos, hdid);
0121
0122         G4double restCurveLength = endCurveLength - curveLength;
0123         if (restCurveLength < GetSmallestFraction() * hstep)
0124         {
0125             succeeded = true;
0126             break;
0127         }
0128         h = std::min(hnext, restCurveLength);
0129     }
0130
0131
0132     if (succeeded)
0133     {
0134         track.LoadFromArray(y, Base::GetStepper()->GetNumberOfVariables());
0135         track.SetCurveLength(track.GetCurveLength() + curveLength);
0136     }
0137
0138     return succeeded;
0139 }
0140
0141 // Driver for one Runge-Kutta Step with monitoring of local truncation error
0142 // to ensure accuracy and adjust stepsize. Input are dependent variable
0143 // array y[0,...,5] and its derivative dydx[0,...,5] at the
0144 // starting value of the independent variable x . Also input are stepsize
0145 // to be attempted htry, and the required accuracy eps. On output y and x
0146 // are replaced by their new values, hdid is the stepsize that was actually
0147 // accomplished, and hnext is the estimated next stepsize.
0148 // This is similar to the function rkqs from the book:
0149 // Numerical Recipes in C: The Art of Scientific Computing (NRC), Second
0150 // Edition, by William H. Press, Saul A. Teukolsky, William T.
0151 // Vetterling, and Brian P. Flannery (Cambridge University Press 1992),
0152 // 16.2 Adaptive StepSize Control for Runge-Kutta, p. 719
0153 //
0154 template <class T>
0155 void G4FSALIntegrationDriver<T>::
0156 OneGoodStep(G4double y[],
0157             G4double dydx[],
0158             G4double& curveLength,   // InOut
0159             G4double htry,
0160             G4double eps_rel_max,
0161             G4double& hdid,         // Out
0162             G4double& hnext)        // Out
0163 {
0164     G4double error2 = DBL_MAX;
0165
0166     G4double yError[G4FieldTrack::ncompSVEC],
0167              yOut[G4FieldTrack::ncompSVEC],
0168              dydxOut[G4FieldTrack::ncompSVEC];
0169
0170     // Set stepsize to the initial trial value
0171     G4double hstep = htry;
0172
0173     const G4int max_trials = 100;
0174
0175     for (G4int iter = 0; iter < max_trials; ++iter)
0176     {
0177         Base::GetStepper()->Stepper(y, dydx, hstep, yOut, yError, dydxOut);
0178         error2 = field_utils::relativeError2(y, yError, hstep, eps_rel_max);
0179
0180         // Step succeeded.
0181         if (error2 <= 1) {  break; }
0182
0183         hstep = Base::ShrinkStepSize2(hstep, error2);
0184     }
0185
0186     hnext = Base::GrowStepSize2(hstep, error2);
0187     curveLength += (hdid = hstep);
0188
0189     for(G4int k = 0; k < Base::GetStepper()->GetNumberOfVariables(); ++k)
0190     {
0191         y[k] = yOut[k];
0192         dydx[k] = dydxOut[k];
0193     }
0194 }
0195
0196 template <class T>
0197 G4bool G4FSALIntegrationDriver<T>::
0198 QuickAdvance( G4FieldTrack& fieldTrack, const G4double dydxIn[],
0199               G4double hstep,
0200               G4double& dchord_step, G4double& dyerr )
0201 {
0202     ++fNoQuickAvanceCalls;
0203
0204     if (hstep == 0)
0205     {
0206         std::ostringstream message;
0207         message << "Proposed step is zero; hstep = " << hstep << " !";
0208         G4Exception("G4FSALIntegrationDriver ::QuickAdvance()",
0209                   "GeomField1001", JustWarning, message);
0210         return true;
0211     }
0212     if (hstep < 0)
0213     {
0214         std::ostringstream message;
0215         message << "Invalid run condition." << G4endl
0216                 << "Proposed step is negative; hstep = "
0217                 << hstep << "." << G4endl
0218                 << "Requested step cannot be negative! Aborting event.";
0219         G4Exception("G4FSALIntegrationDriver ::QuickAdvance()",
0220                     "GeomField0003", EventMustBeAborted, message);
0221         return false;
0222     }
0223
0224     G4double yError[G4FieldTrack::ncompSVEC],
0225              yIn[G4FieldTrack::ncompSVEC],
0226              yOut[G4FieldTrack::ncompSVEC],
0227              dydxOut[G4FieldTrack::ncompSVEC];
0228
0229     fieldTrack.DumpToArray(yIn);
0230
0231     Base::GetStepper()->Stepper(yIn, dydxIn, hstep, yOut, yError, dydxOut);
0232     dchord_step = Base::GetStepper()->DistChord();
0233
0234     fieldTrack.LoadFromArray(yOut, Base::GetStepper()->GetNumberOfVariables());
0235     fieldTrack.SetCurveLength(fieldTrack.GetCurveLength() + hstep);
0236
0237     dyerr = field_utils::absoluteError(yOut, yError, hstep);
0238
0239     return true;
0240 }
0241
0242 template <class T>
0243 void G4FSALIntegrationDriver<T>::SetSmallestFraction(G4double newFraction)
0244 {
0245     if( newFraction > 1.e-16 && newFraction < 1e-8 )
0246     {
0247         fSmallestFraction = newFraction;
0248     }
0249     else
0250     {
0251         std::ostringstream message;
0252         message << "Smallest Fraction not changed. " << G4endl
0253                 << "  Proposed value was " << newFraction << G4endl
0254                 << "  Value must be between 1.e-8 and 1.e-16";
0255         G4Exception("G4FSALIntegrationDriver::SetSmallestFraction()",
0256                     "GeomField1001", JustWarning, message);
0257     }
0258 }
0259
0260 template <class T>
0261 void G4FSALIntegrationDriver<T>::CheckStep(
0262     const G4ThreeVector& posIn, const G4ThreeVector& posOut, G4double hdid)
0263 {
0264     const G4double endPointDist = (posOut - posIn).mag();
0265     if (endPointDist >= hdid * (1. + CLHEP::perMillion))
0266     {
0267         ++fNoAccurateAdvanceBadSteps;
0268 #ifdef G4DEBUG_FIELD
0269         // Issue a warning only for gross differences -
0270         // we understand how small difference occur.
0271         if (endPointDist >= hdid * (1. + perThousand))
0272         {
0273            G4Exception("G4FSALIntegrationDriver::CheckStep()",
0274                        "GeomField1002", JustWarning,
0275                        "endPointDist >= hdid!");
0276         }
0277 #endif
0278     }
0279     else
0280     {
0281        ++fNoAccurateAdvanceGoodSteps;
0282     }
0283 }
0284
0285 template <class T>
0286 inline G4double G4FSALIntegrationDriver<T>::GetMinimumStep() const
0287 {
0288     return fMinimumStep;
0289 }
0290
0291 template <class T>
0292 void G4FSALIntegrationDriver<T>::SetMinimumStep(G4double minimumStepLength)
0293 {
0294     fMinimumStep = minimumStepLength;
0295 }
0296
0297 template <class T>
0298 G4int G4FSALIntegrationDriver<T>::GetVerboseLevel() const
0299 {
0300     return fVerboseLevel;
0301 }
0302
0303 template <class T>
0304 void G4FSALIntegrationDriver<T>::SetVerboseLevel(G4int newLevel)
0305 {
0306     fVerboseLevel = newLevel;
0307 }
0308
0309 template <class T>
0310 G4double G4FSALIntegrationDriver<T>::GetSmallestFraction() const
0311 {
0312     return fSmallestFraction;
0313 }
0314
0315 template <class T>
0316 void G4FSALIntegrationDriver<T>::IncrementQuickAdvanceCalls()
0317 {
0318     ++fNoQuickAvanceCalls;
0319 }
0320
0321 template <class T>
0322 G4double
0323 G4FSALIntegrationDriver<T>::AdvanceChordLimited(G4FieldTrack& track,
0324                                              G4double hstep,
0325                                              G4double eps,
0326                                              G4double chordDistance)
0327 {
0328    return ChordFinderDelegate::AdvanceChordLimitedImpl(track, hstep,
0329                                                        eps, chordDistance);
0330 }
0331
0332 template <class T>
0333 void G4FSALIntegrationDriver<T>::StreamInfo( std::ostream& os ) const
0334 {
0335 // Write out the parameters / state of the driver
0336   os << "State of G4IntegrationDriver: " << std::endl;
0337   os << "--Base state (G4RKIntegrationDriver): " << std::endl;
0338   Base::StreamInfo( os );
0339   os << "--Own  state (G4IntegrationDriver<>): " << std::endl;
0340   os << "    fMinimumStep =      " << fMinimumStep  << std::endl;
0341   os << "    Smallest Fraction = " << fSmallestFraction << std::endl;
0342
0343   os << "    verbose level     = " << fVerboseLevel << std::endl;
0344   os << "    Reintegrates      = " << DoesReIntegrate() << std::endl;
0345   os << "--Chord Finder Delegate state: " << std::endl;
0346   ChordFinderDelegate::StreamDelegateInfo( os );
0347 }