include/Geant4/G4TMagErrorStepper.hh

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 // G4TMagErrorStepper
0027 //
0028 // Class description:
0029 //
0030 // Templated version of G4MagErrorStepper.
0031 // Adapted from G4G4TMagErrorStepper class.
0032
0033 // Author: Josh Xie (CERN, Google Summer of Code 2014), June 2014
0034 // Supervisors:  Sandro Wenzel, John Apostolakis (CERN)
0035 // --------------------------------------------------------------------
0036 #ifndef G4TMAG_ERROR_STEPPER_HH
0037 #define G4TMAG_ERROR_STEPPER_HH
0038
0039 #include "G4Types.hh"
0040 #include "G4MagIntegratorStepper.hh"
0041 #include "G4ThreeVector.hh"
0042 #include "G4LineSection.hh"
0043
0044 /**
0045  * @brief G4TMagErrorStepper is a templated version of G4MagErrorStepper.
0046  */
0047
0048 template <class T_Stepper, class T_Equation, unsigned int N>
0049 class G4TMagErrorStepper : public G4MagIntegratorStepper
0050 {
0051   public:
0052
0053     G4TMagErrorStepper(T_Equation* EqRhs, G4int numberOfVariables,
0054                        G4int numStateVariables = 12)
0055       : G4MagIntegratorStepper(EqRhs, numberOfVariables, numStateVariables)
0056       , fEquation_Rhs(EqRhs) { ; }
0057
0058     virtual ~G4TMagErrorStepper() = default;
0059
0060     G4TMagErrorStepper(const G4TMagErrorStepper&) = delete;
0061     G4TMagErrorStepper& operator=(const G4TMagErrorStepper&) = delete;
0062
0063     inline void RightHandSide(G4double y[], G4double dydx[])
0064     {
0065       fEquation_Rhs->T_Equation::RightHandSide(y, dydx);
0066     }
0067
0068     inline void Stepper(const G4double yInput[], const G4double dydx[],
0069                         G4double hstep, G4double yOutput[], G4double yError[]) override final;
0070
0071     inline G4double DistChord() const override final;
0072     G4StepperType StepperType() const override { return kTMagErrorStepper; }
0073
0074   private:
0075
0076     // STATE
0077     G4ThreeVector fInitialPoint, fMidPoint, fFinalPoint;
0078     // Data stored in order to find the chord
0079
0080     // Dependent Objects, owned --- part of the STATE
0081     G4double yInitial[N < 8 ? 8 : N];
0082     G4double yMiddle[N < 8 ? 8 : N];
0083     G4double dydxMid[N < 8 ? 8 : N];
0084     G4double yOneStep[N < 8 ? 8 : N];
0085     // The following arrays are used only for temporary storage
0086     // they are allocated at the class level only for efficiency -
0087     // so that calls to new and delete are not made in Stepper().
0088
0089     T_Equation* fEquation_Rhs;
0090 };
0091
0092 // ------------   Implementation -----------------------
0093
0094 template <class T_Stepper, class T_Equation, unsigned int N >
0095 void G4TMagErrorStepper<T_Stepper,T_Equation,N>::
0096 Stepper(const G4double yInput[],
0097         const G4double dydx[],
0098               G4double hstep,
0099               G4double yOutput[],
0100               G4double yError[])
0101 // The stepper for the Runge Kutta integration. The stepsize
0102 // is fixed, with the Step size given by hstep.
0103 // Integrates ODE starting values y[0 to N].
0104 // Outputs yout[] and its estimated error yerr[].
0105 {
0106   const unsigned int maxvar = GetNumberOfStateVariables();
0107
0108   //  Saving yInput because yInput and yOutput can be aliases for same array
0109   for(unsigned int i = 0; i < N; ++i)
0110      yInitial[i] = yInput[i];
0111   yInitial[7] =
0112      yInput[7];  // Copy the time in case ... even if not really needed
0113   yMiddle[7]  = yInput[7];  // Copy the time from initial value
0114   yOneStep[7] = yInput[7];  // As it contributes to final value of yOutput ?
0115   // yOutput[7] = yInput[7];  // -> dumb stepper does it too for RK4
0116   for(unsigned int i = N; i < maxvar; ++i)
0117      yOutput[i] = yInput[i];
0118
0119   G4double halfStep = hstep * 0.5;
0120
0121   // Do two half steps
0122
0123   static_cast<T_Stepper*>(this)->DumbStepper(yInitial, dydx, halfStep,
0124                                                yMiddle);
0125   this->RightHandSide(yMiddle, dydxMid);
0126   static_cast<T_Stepper*>(this)->DumbStepper(yMiddle, dydxMid, halfStep,
0127                                              yOutput);
0128
0129   // Store midpoint, chord calculation
0130
0131   fMidPoint = G4ThreeVector(yMiddle[0], yMiddle[1], yMiddle[2]);
0132
0133   // Do a full Step
0134   static_cast<T_Stepper*>(this)->DumbStepper(yInitial, dydx, hstep, yOneStep);
0135   for(unsigned int i = 0; i < N; ++i)
0136   {
0137     yError[i] = yOutput[i] - yOneStep[i];
0138     yOutput[i] +=
0139        yError[i] *
0140        T_Stepper::IntegratorCorrection;  // Provides accuracy increased
0141     // by 1 order via the
0142     // Richardson Extrapolation
0143   }
0144
0145   fInitialPoint = G4ThreeVector(yInitial[0], yInitial[1], yInitial[2]);
0146   fFinalPoint   = G4ThreeVector(yOutput[0], yOutput[1], yOutput[2]);
0147
0148   return;
0149 }
0150
0151
0152 template <class T_Stepper, class T_Equation, unsigned int N >
0153 inline G4double
0154 G4TMagErrorStepper<T_Stepper,T_Equation,N>::DistChord() const
0155 {
0156   // Estimate the maximum distance from the curve to the chord
0157   //
0158   //  We estimate this using the distance of the midpoint to
0159   //  chord (the line between
0160   //
0161   //  Method below is good only for angle deviations < 2 pi,
0162   //   This restriction should not a problem for the Runge cutta methods,
0163   //   which generally cannot integrate accurately for large angle deviations.
0164   G4double distLine, distChord;
0165
0166   if(fInitialPoint != fFinalPoint)
0167   {
0168     distLine = G4LineSection::Distline(fMidPoint, fInitialPoint, fFinalPoint);
0169     // This is a class method that gives distance of Mid
0170     //  from the Chord between the Initial and Final points.
0171
0172     distChord = distLine;
0173   }
0174   else
0175   {
0176      distChord = (fMidPoint - fInitialPoint).mag();
0177   }
0178
0179   return distChord;
0180 }
0181
0182 #endif