include/Geant4/G4RKIntegrationDriver.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 // G4RKIntegrationDriver inline implementation
0027 //
0028 // Created: D.Sorokin, 2017
0029 // --------------------------------------------------------------------
0030
0031 #include "globals.hh"
0032
0033 template <class T>
0034 G4RKIntegrationDriver<T>::G4RKIntegrationDriver(T* pStepper)
0035 {
0036     RenewStepperAndAdjustImpl(pStepper);
0037     fMaxStepBase = 250;
0038     fMaxNoSteps = fMaxStepBase / pIntStepper->IntegratorOrder();
0039 }
0040
0041 // Step failed; compute the size of retrial Step.
0042 //
0043 template <class T> G4double G4RKIntegrationDriver<T>::
0044 ShrinkStepSize(G4double h, G4double error) const
0045 {
0046     if (error > errorConstraintShrink)
0047     {
0048         return max_stepping_decrease * h;
0049     }
0050     return GetSafety() * h * std::pow(error, GetPshrnk());
0051 }
0052
0053 template <class T> G4double G4RKIntegrationDriver<T>::
0054 ShrinkStepSize2(G4double h, G4double error2) const
0055 {
0056     if (error2 > errorConstraintShrink * errorConstraintShrink)
0057     {
0058         return max_stepping_decrease * h;
0059     }
0060     return GetSafety() * h * std::pow(error2, 0.5 * GetPshrnk());
0061 }
0062
0063 // Compute size of next Step
0064 //
0065 template<class T> G4double G4RKIntegrationDriver<T>::
0066 GrowStepSize(G4double h, G4double error) const
0067 {
0068     if (error < errorConstraintGrow)
0069     {
0070         return max_stepping_increase * h;
0071     }
0072     return GetSafety() * h * std::pow(error, GetPgrow());
0073 }
0074
0075 template<class T> G4double G4RKIntegrationDriver<T>::
0076 GrowStepSize2(G4double h, G4double error2) const
0077 {
0078     if (error2 < errorConstraintGrow * errorConstraintGrow)
0079     {
0080         return max_stepping_increase * h;
0081     }
0082     return GetSafety() * h * std::pow(error2, 0.5 * GetPgrow());
0083 }
0084
0085 template <class T>
0086 G4double G4RKIntegrationDriver<T>::
0087 ComputeNewStepSize( G4double  errMaxNorm,    // max error  (normalised)
0088                     G4double  hstepCurrent ) // current step size
0089 {
0090     if (errMaxNorm > 1)
0091     {
0092         return ShrinkStepSize(hstepCurrent, errMaxNorm);
0093     }
0094     if (errMaxNorm >= 0)
0095     {
0096         return GrowStepSize(hstepCurrent, errMaxNorm);
0097     }
0098
0099     G4Exception("G4RKIntegrationDriver::ComputeNewStepSize", "GeomField0003",
0100                 FatalException, "Error is negative!");
0101
0102     return max_stepping_increase * hstepCurrent;
0103 }
0104
0105 template <class T>
0106 void G4RKIntegrationDriver<T>::
0107 GetDerivatives( const G4FieldTrack& track, G4double dydx[] ) const
0108 {
0109     G4double y[G4FieldTrack::ncompSVEC];
0110     track.DumpToArray(y);
0111     pIntStepper->RightHandSide(y, dydx);
0112 }
0113
0114 template <class T>
0115 void G4RKIntegrationDriver<T>::GetDerivatives(const G4FieldTrack& track,
0116                                               G4double dydx[],
0117                                               G4double field[]) const
0118 {
0119     G4double y[G4FieldTrack::ncompSVEC];
0120     track.DumpToArray(y);
0121     pIntStepper->RightHandSide(y, dydx, field);
0122 }
0123
0124 template <class T>
0125 void G4RKIntegrationDriver<T>::UpdateErrorConstraints()
0126 {
0127     errorConstraintShrink = std::pow(
0128             max_stepping_decrease / GetSafety(), 1. / GetPshrnk());
0129
0130     errorConstraintGrow = std::pow(
0131             max_stepping_increase / GetSafety(), 1. / GetPgrow());
0132 }
0133
0134 template <class T>
0135 inline G4double G4RKIntegrationDriver<T>::GetSafety() const
0136 {
0137     return safety;
0138 }
0139
0140 template <class T>
0141 inline G4double G4RKIntegrationDriver<T>::GetPshrnk() const
0142 {
0143     return pshrnk;
0144 }
0145
0146 template <class T>
0147 G4double G4RKIntegrationDriver<T>::GetPgrow() const
0148 {
0149     return pgrow;
0150 }
0151
0152 template <class T>
0153 void G4RKIntegrationDriver<T>::ReSetParameters(G4double new_safety)
0154 {
0155     safety = new_safety;
0156     pshrnk = -1.0 / pIntStepper->IntegratorOrder();
0157     pgrow  = -1.0 / (1.0 + pIntStepper->IntegratorOrder());
0158     UpdateErrorConstraints();
0159 }
0160
0161 template <class T>
0162 void G4RKIntegrationDriver<T>::SetSafety(G4double val)
0163 {
0164     safety = val;
0165     UpdateErrorConstraints();
0166 }
0167
0168 template <class T> void G4RKIntegrationDriver<T>::
0169 RenewStepperAndAdjust(G4MagIntegratorStepper* stepper)
0170 {
0171    T* ourStepper = dynamic_cast<T*>(stepper);
0172    if ( ourStepper )
0173    {
0174       RenewStepperAndAdjustImpl( ourStepper );
0175    }
0176    else
0177    {
0178       G4Exception("G4RKIntegrationDriver::RenewStepperAndAdjust()",
0179                   "GeomField0002", FatalException,
0180      "The type of the stepper provided is incorrect for this templated driver");
0181    }
0182 }
0183
0184 template <class T>
0185 void G4RKIntegrationDriver<T>::RenewStepperAndAdjustImpl(T* stepper)
0186 {
0187     pIntStepper = stepper;
0188     ReSetParameters();
0189 }
0190
0191 template <class T>
0192 const T* G4RKIntegrationDriver<T>::GetStepper() const
0193 {
0194     return pIntStepper;
0195 }
0196
0197 template <class T>
0198 T* G4RKIntegrationDriver<T>::GetStepper()
0199 {
0200     return pIntStepper;
0201 }
0202
0203 template <class T>
0204 G4int G4RKIntegrationDriver<T>::GetMaxNoSteps() const
0205 {
0206     return fMaxNoSteps;
0207 }
0208
0209 template <class T>
0210 void G4RKIntegrationDriver<T>::SetMaxNoSteps(G4int val)
0211 {
0212     fMaxNoSteps = val;
0213 }
0214
0215 template <class T>
0216 G4EquationOfMotion* G4RKIntegrationDriver<T>::GetEquationOfMotion()
0217 {
0218     return pIntStepper->GetEquationOfMotion();
0219 }
0220
0221 template <class T>
0222 void G4RKIntegrationDriver<T>::SetEquationOfMotion(G4EquationOfMotion* equation)
0223 {
0224     pIntStepper->SetEquationOfMotion(equation);
0225 }
0226
0227 template <class T>
0228 void G4RKIntegrationDriver<T>::StreamInfo( std::ostream& os ) const
0229 {
0230     os << "State of G4RKIntegrationDriver: " << std::endl;
0231     os << "  Max number of Steps = " << fMaxNoSteps << std::endl;
0232     os << "  Safety factor       = " << safety  << std::endl;
0233     os << "  Power - shrink      = " << pshrnk << std::endl;
0234     os << "  Power - grow        = " << pgrow << std::endl;
0235     os << "  threshold - shrink  = " << errorConstraintShrink << std::endl;
0236     os << "  threshold - grow    = " << errorConstraintGrow   << std::endl;
0237 }