G4DormandPrince745 Class Reference

#include <G4DormandPrince745.hh>

Inheritance diagram for G4DormandPrince745:

Public Member Functions
	G4DormandPrince745 (G4EquationOfMotion *equation, G4int numberOfVariables=6)
virtual void	Stepper (const G4double yInput[], const G4double dydx[], G4double hstep, G4double yOutput[], G4double yError[]) override
void	Stepper (const G4double yInput[], const G4double dydx[], G4double hstep, G4double yOutput[], G4double yError[], G4double dydxOutput[])
void	SetupInterpolation ()
void	Interpolate (G4double tau, G4double yOut[]) const
virtual G4double	DistChord () const override
virtual G4int	IntegratorOrder () const override
const G4String &	StepperType () const
const G4String &	StepperDescription () const
const field_utils::State &	GetYOut () const
void	Interpolate4thOrder (G4double yOut[], G4double tau) const
void	SetupInterpolation5thOrder ()
void	Interpolate5thOrder (G4double yOut[], G4double tau) const
G4EquationOfMotion *	GetSpecificEquation ()
Public Member Functions inherited from G4MagIntegratorStepper
	G4MagIntegratorStepper (G4EquationOfMotion *Equation, G4int numIntegrationVariables, G4int numStateVariables=12, G4bool isFSAL=false)
virtual	~G4MagIntegratorStepper ()=default
	G4MagIntegratorStepper (const G4MagIntegratorStepper &)=delete
G4MagIntegratorStepper &	operator= (const G4MagIntegratorStepper &)=delete
virtual void	Stepper (const G4double y[], const G4double dydx[], G4double h, G4double yout[], G4double yerr[])=0
virtual G4double	DistChord () const =0
void	NormaliseTangentVector (G4double vec[6])
void	NormalisePolarizationVector (G4double vec[12])
void	RightHandSide (const G4double y[], G4double dydx[]) const
void	RightHandSide (const G4double y[], G4double dydx[], G4double field[]) const
G4int	GetNumberOfVariables () const
G4int	GetNumberOfStateVariables () const
virtual G4int	IntegratorOrder () const =0
G4int	IntegrationOrder ()
G4EquationOfMotion *	GetEquationOfMotion ()
const G4EquationOfMotion *	GetEquationOfMotion () const
void	SetEquationOfMotion (G4EquationOfMotion *newEquation)
unsigned long	GetfNoRHSCalls ()
void	ResetfNORHSCalls ()
G4bool	IsFSAL () const

Additional Inherited Members
Protected Member Functions inherited from G4MagIntegratorStepper
void	SetIntegrationOrder (G4int order)
void	SetFSAL (G4bool flag=true)

Detailed Description

Definition at line 46 of file G4DormandPrince745.hh.

Constructor & Destructor Documentation

G4DormandPrince745()

G4DormandPrince745::G4DormandPrince745	(	G4EquationOfMotion *	equation,
		G4int	numberOfVariables = 6
	)

Definition at line 67 of file G4DormandPrince745.cc.

    : G4MagIntegratorStepper(equation, noIntegrationVariables)
{
}

Member Function Documentation

DistChord()

G4double G4DormandPrince745::DistChord ( ) const

overridevirtual

Implements G4MagIntegratorStepper.

Definition at line 206 of file G4DormandPrince745.cc.

{
    // Coefficients were taken from Some Practical Runge-Kutta Formulas
    // by Lawrence F. Shampine, page 149, c*
    //
    const G4double hf1 = 6025192743.0 / 30085553152.0,
                   hf3 = 51252292925.0 / 65400821598.0,
                   hf4 = - 2691868925.0 / 45128329728.0,
                   hf5 = 187940372067.0 / 1594534317056.0,
                   hf6 = - 1776094331.0 / 19743644256.0,
                   hf7 = 11237099.0 / 235043384.0;
 
    G4ThreeVector mid;
 
    for(G4int i = 0; i < 3; ++i) 
    {
        mid[i] = fyIn[i] + 0.5 * fLastStepLength * (
            hf1 * fdydxIn[i] + hf3 * ak3[i] + 
            hf4 * ak4[i] + hf5 * ak5[i] + hf6 * ak6[i] + hf7 * ak7[i]);
    }
 
    const G4ThreeVector begin = makeVector(fyIn, Value3D::Position);
    const G4ThreeVector end = makeVector(fyOut, Value3D::Position);
 
    return G4LineSection::Distline(mid, begin, end);
}

GetSpecificEquation()

G4EquationOfMotion * G4DormandPrince745::GetSpecificEquation ( )

inline

Definition at line 88 of file G4DormandPrince745.hh.

{ return GetEquationOfMotion(); }

GetYOut()

const field_utils::State & G4DormandPrince745::GetYOut ( ) const

inline

Definition at line 81 of file G4DormandPrince745.hh.

{ return fyOut; }

IntegratorOrder()

virtual G4int G4DormandPrince745::IntegratorOrder ( ) const

inlineoverridevirtual

Implements G4MagIntegratorStepper.

Definition at line 76 of file G4DormandPrince745.hh.

{ return 4; }

Interpolate()

void G4DormandPrince745::Interpolate ( G4double  tau, G4double  yOut[]  ) const

inline

Definition at line 68 of file G4DormandPrince745.hh.

    {
        Interpolate4thOrder(yOut, tau);
    }

Interpolate4thOrder()

void G4DormandPrince745::Interpolate4thOrder	(	G4double	yOut[],
		G4double	tau
	)		const

Definition at line 238 of file G4DormandPrince745.cc.

{
    const G4int numberOfVariables = GetNumberOfVariables();
    
    const G4double tau2 = tau * tau,
                   tau3 = tau * tau2,
                   tau4 = tau2 * tau2;
 
    const G4double bf1 = 1.0 / 11282082432.0 * (
        157015080.0 * tau4 - 13107642775.0 * tau3 + 34969693132.0 * tau2 - 
        32272833064.0 * tau + 11282082432.0);
 
    const G4double bf3 = - 100.0 / 32700410799.0 * tau * (
        15701508.0 * tau3 - 914128567.0 * tau2 + 2074956840.0 * tau - 
        1323431896.0);
    
    const G4double bf4 = 25.0 / 5641041216.0 * tau * (
        94209048.0 * tau3 - 1518414297.0 * tau2 + 2460397220.0 * tau - 
        889289856.0);
    
    const G4double bf5 = - 2187.0 / 199316789632.0 * tau * (
        52338360.0 * tau3 - 451824525.0 * tau2 + 687873124.0 * tau - 
        259006536.0);
 
    const G4double bf6 = 11.0 / 2467955532.0 * tau * (
        106151040.0 * tau3 - 661884105.0 * tau2 + 
        946554244.0 * tau - 361440756.0);
    
    const G4double bf7 = 1.0 / 29380423.0 * tau * (1.0 - tau) * (
        8293050.0 * tau2 - 82437520.0 * tau + 44764047.0);
 
    for(G4int i = 0; i < numberOfVariables; ++i)
    {
        yOut[i] = fyIn[i] + fLastStepLength * tau * (
            bf1 * fdydxIn[i] + bf3 * ak3[i] + bf4 * ak4[i] + 
            bf5 * ak5[i] + bf6 * ak6[i] + bf7 * ak7[i]);
    }
}

Referenced by Interpolate().

Interpolate5thOrder()

void G4DormandPrince745::Interpolate5thOrder	(	G4double	yOut[],
		G4double	tau
	)		const

Definition at line 334 of file G4DormandPrince745.cc.

{
    // Define the coefficients for the polynomials
    //
    G4double bi[10][5];
    
    //  COEFFICIENTS OF   bi[1]
    bi[1][0] = 1.0,
    bi[1][1] = -38039.0 / 7040.0,
    bi[1][2] = 125923.0 / 10560.0,
    bi[1][3] = -19683.0 / 1760.0,
    bi[1][4] = 3303.0 / 880.0,
    //  --------------------------------------------------------
    //
    //  COEFFICIENTS OF  bi[2]
    bi[2][0] = 0.0,
    bi[2][1] = 0.0,
    bi[2][2] = 0.0,
    bi[2][3] = 0.0,
    bi[2][4] = 0.0,
    //  --------------------------------------------------------
    //
    //  COEFFICIENTS OF  bi[3]
    bi[3][0] = 0.0,
    bi[3][1] = -12500.0 / 4081.0,
    bi[3][2] = 205000.0 / 12243.0,
    bi[3][3] = -90000.0 / 4081.0,
    bi[3][4] = 36000.0 / 4081.0,
    //  --------------------------------------------------------
    //
    //  COEFFICIENTS OF  bi[4]
    bi[4][0] = 0.0,
    bi[4][1] = -3125.0 / 704.0,
    bi[4][2] = 25625.0 / 1056.0,
    bi[4][3] = -5625.0 / 176.0,
    bi[4][4] = 1125.0 / 88.0,
    //  --------------------------------------------------------
    //
    //  COEFFICIENTS OF  bi[5]
    bi[5][0] = 0.0,
    bi[5][1] = 164025.0 / 74624.0,
    bi[5][2] = -448335.0 / 37312.0,
    bi[5][3] = 295245.0 / 18656.0,
    bi[5][4] = -59049.0 / 9328.0,
    //  --------------------------------------------------------
    //
    //  COEFFICIENTS OF  bi[6]
    bi[6][0] = 0.0,
    bi[6][1] = -25.0 / 28.0,
    bi[6][2] = 205.0 / 42.0,
    bi[6][3] = -45.0 / 7.0,
    bi[6][4] = 18.0 / 7.0,
    //  --------------------------------------------------------
    //
    //  COEFFICIENTS OF  bi[7]
    bi[7][0] = 0.0,
    bi[7][1] = -2.0 / 11.0,
    bi[7][2] = 73.0 / 55.0,
    bi[7][3] = -171.0 / 55.0,
    bi[7][4] = 108.0 / 55.0,
    //  --------------------------------------------------------
    //
    //  COEFFICIENTS OF  bi[8]
    bi[8][0] = 0.0,
    bi[8][1] = 189.0 / 22.0,
    bi[8][2] = -1593.0 / 55.0,
    bi[8][3] = 3537.0 / 110.0,
    bi[8][4] = -648.0 / 55.0,
    //  --------------------------------------------------------
    //
    //  COEFFICIENTS OF  bi[9]
    bi[9][0] = 0.0,
    bi[9][1] = 351.0 / 110.0,
    bi[9][2] = -999.0 / 55.0,
    bi[9][3] = 2943.0 / 110.0,
    bi[9][4] = -648.0 / 55.0;
    //  --------------------------------------------------------
        
    // Calculating the polynomials
 
    G4double b[10];
    std::memset(b, 0.0, sizeof(b));
 
    G4double tauPower = 1.0;
    for(G4int j = 0; j <= 4; ++j)
    {       
       for(G4int iStage = 1; iStage <= 9; ++iStage)
       {
          b[iStage] += bi[iStage][j] * tauPower;
       }
       tauPower *= tau;       
    }
    
    const G4int numberOfVariables = GetNumberOfVariables();
    const G4double stepLen = fLastStepLength * tau;
    for(G4int i = 0; i < numberOfVariables; ++i)
    {
        yOut[i] = fyIn[i] + stepLen * (
            b[1] * fdydxIn[i] + b[2] * ak2[i] + b[3] * ak3[i] +
            b[4] * ak4[i] + b[5] * ak5[i] + b[6] * ak6[i] +
            b[7] * ak7[i] + b[8] * ak8[i] + b[9] * ak9[i]
        );
    }
}

SetupInterpolation()

void G4DormandPrince745::SetupInterpolation ( )

inline

Definition at line 66 of file G4DormandPrince745.hh.

{}

SetupInterpolation5thOrder()

void G4DormandPrince745::SetupInterpolation5thOrder

(

)

Definition at line 285 of file G4DormandPrince745.cc.

{
    // Coefficients for the additional stages
    //
    const G4double b81 = 6245.0 / 62208.0,
                   b82 = 0.0,
                   b83 = 8875.0 / 103032.0,
                   b84 = -125.0 / 1728.0,
                   b85 = 801.0 / 13568.0,
                   b86 = -13519.0 / 368064.0,
                   b87 = 11105.0 / 368064.0,
        
                   b91 = 632855.0 / 4478976.0,
                   b92 = 0.0,
                   b93 = 4146875.0 / 6491016.0,
                   b94 = 5490625.0 /14183424.0,
                   b95 = -15975.0 / 108544.0,
                   b96 = 8295925.0 / 220286304.0,
                   b97 = -1779595.0 / 62938944.0,
                   b98 = -805.0 / 4104.0;
    
    const G4int numberOfVariables = GetNumberOfVariables();
    State yTemp;
 
    // Evaluate the extra stages
    //
    for(G4int i = 0; i < numberOfVariables; ++i)
    {
        yTemp[i] = fyIn[i] + fLastStepLength * (
            b81 * fdydxIn[i] + b82 * ak2[i] + b83 * ak3[i] +
            b84 * ak4[i] + b85 * ak5[i] + b86 * ak6[i] +
            b87 * ak7[i]
        );
    }
    RightHandSide(yTemp, ak8);          // 8th Stage
    
    for(G4int i = 0; i < numberOfVariables; ++i)
    {
        yTemp[i] = fyIn[i] + fLastStepLength * (
            b91 * fdydxIn[i] + b92 * ak2[i] + b93 * ak3[i] +
            b94 * ak4[i] + b95 * ak5[i] + b96 * ak6[i] +
            b97 * ak7[i] + b98 * ak8[i]
        );
    }
    RightHandSide(yTemp, ak9);          // 9th Stage
}

Stepper() [1/2]

void G4DormandPrince745::Stepper ( const G4double  yInput[], const G4double  dydx[], G4double  hstep, G4double  yOutput[], G4double  yError[]  )

overridevirtual

Implements G4MagIntegratorStepper.

Definition at line 89 of file G4DormandPrince745.cc.

{
    // The various constants defined on the basis of butcher tableu
    //
    const G4double b21 = 0.2,
                   b31 = 3.0 / 40.0, b32 = 9.0 / 40.0,
                   b41 = 44.0 / 45.0, b42 = -56.0 / 15.0, b43 = 32.0/9.0,
 
        b51 = 19372.0 / 6561.0, b52 = -25360.0 / 2187.0, b53 = 64448.0 / 6561.0,
        b54 = -212.0 / 729.0,
        
        b61 = 9017.0 / 3168.0 , b62 =   -355.0 / 33.0,
        b63 = 46732.0 / 5247.0, b64 = 49.0 / 176.0,
        b65 = -5103.0 / 18656.0,
        
        b71 = 35.0 / 384.0, b72 = 0.,
        b73 = 500.0 / 1113.0, b74 = 125.0 / 192.0,
        b75 = -2187.0 / 6784.0, b76 = 11.0 / 84.0,
    
    //Sum of columns, sum(bij) = ei
    //    e1 = 0. ,
    //    e2 = 1.0/5.0 ,
    //    e3 = 3.0/10.0 ,
    //    e4 = 4.0/5.0 ,
    //    e5 = 8.0/9.0 ,
    //    e6 = 1.0 ,
    //    e7 = 1.0 ,
    
    // Difference between the higher and the lower order method coeff. :
    // b7j are the coefficients of higher order
    
        dc1 = -(b71 - 5179.0 / 57600.0),
        dc2 = -(b72 - .0),
        dc3 = -(b73 - 7571.0 / 16695.0),
        dc4 = -(b74 - 393.0 / 640.0),
        dc5 = -(b75 + 92097.0 / 339200.0),
        dc6 = -(b76 - 187.0 / 2100.0),
        dc7 = -(- 1.0 / 40.0);
    
    const G4int numberOfVariables = GetNumberOfVariables();
    State yTemp;
    
    // The number of variables to be integrated over
    //
    yOut[7] = yTemp[7]  = yInput[7];
 
    //  Saving yInput because yInput and yOut can be aliases for same array
    //
    for(G4int i = 0; i < numberOfVariables; ++i)
    {
        fyIn[i] = yInput[i];
    }
    // RightHandSide(yIn, dydx);    // Not done! 1st stage
    
    for(G4int i = 0; i < numberOfVariables; ++i)
    {
        yTemp[i] = fyIn[i] + b21 * hstep * dydx[i];
    }
    RightHandSide(yTemp, ak2);              // 2nd stage
    
    for(G4int i = 0; i < numberOfVariables; ++i)
    {
        yTemp[i] = fyIn[i] + hstep * (b31 * dydx[i] + b32 * ak2[i]);
    }
    RightHandSide(yTemp, ak3);              // 3rd stage
    
    for(G4int i = 0; i < numberOfVariables; ++i)
    {
        yTemp[i] = fyIn[i] + hstep * (
            b41 * dydx[i] + b42 * ak2[i] + b43 * ak3[i]);
    }
    RightHandSide(yTemp, ak4);              // 4th stage
    
    for(G4int i = 0; i < numberOfVariables; ++i)
    {
        yTemp[i] = fyIn[i] + hstep * (
            b51 * dydx[i] + b52 * ak2[i] + b53 * ak3[i] + b54 * ak4[i]);
    }
    RightHandSide(yTemp, ak5);              // 5th stage
    
    for(G4int i = 0; i < numberOfVariables; ++i)
    {
        yTemp[i] = fyIn[i] + hstep * (
            b61 * dydx[i] + b62 * ak2[i] + 
            b63 * ak3[i] + b64 * ak4[i] + b65 * ak5[i]);
    }
    RightHandSide(yTemp, ak6);              // 6th stage
    
    for(G4int i = 0; i < numberOfVariables; ++i)
    {
        yOut[i] = fyIn[i] + hstep * (
            b71 * dydx[i] + b72 * ak2[i] + b73 * ak3[i] +
            b74 * ak4[i] + b75 * ak5[i] + b76 * ak6[i]);
    }
    RightHandSide(yOut, ak7);               // 7th and Final stage
    
    for(G4int i = 0; i < numberOfVariables; ++i)
    {
        yErr[i] = hstep * (
            dc1 * dydx[i] + dc2 * ak2[i] + 
            dc3 * ak3[i] + dc4 * ak4[i] +
            dc5 * ak5[i] + dc6 * ak6[i] + dc7 * ak7[i]
        ) + 1.5e-18;
 
        // Store Input and Final values, for possible use in calculating chord
        //
        fyOut[i] = yOut[i];
        fdydxIn[i] = dydx[i];        
    }
    
    fLastStepLength = hstep;
}

Referenced by Stepper().

Stepper() [2/2]

void G4DormandPrince745::Stepper	(	const G4double	yInput[],
		const G4double	dydx[],
		G4double	hstep,
		G4double	yOutput[],
		G4double	yError[],
		G4double	dydxOutput[]
	)

Definition at line 73 of file G4DormandPrince745.cc.

{
  Stepper(yInput, dydx, hstep, yOutput, yError);
  field_utils::copy(dydxOutput, ak7);
}

StepperDescription()

const G4String & G4DormandPrince745::StepperDescription ( ) const

inline

Definition at line 79 of file G4DormandPrince745.hh.

{ return gStepperDescription; }

StepperType()

const G4String & G4DormandPrince745::StepperType ( ) const

inline

Definition at line 78 of file G4DormandPrince745.hh.

{ return gStepperType; }

---

The documentation for this class was generated from the following files:

• G4DormandPrince745.hh
• G4DormandPrince745.cc