G4CashKarpRKF45 Class Reference

G4CashKarpRKF45 implements the Cash-Karp Runge-Kutta-Fehlberg 4/5 method, an embedded fourth order method (giving fifth-order accuracy) for the solution of an ODE. Two different fourth order estimates are calculated; their difference gives an error estimate. It is used to integrate the equations of the motion of a particle in a magnetic field. More...

#include <G4CashKarpRKF45.hh>

Inheritance diagram for G4CashKarpRKF45:

Public Member Functions
	G4CashKarpRKF45 (G4EquationOfMotion *EqRhs, G4int numberOfVariables=6, G4bool primary=true)
	~G4CashKarpRKF45 () override
	G4CashKarpRKF45 (const G4CashKarpRKF45 &)=delete
G4CashKarpRKF45 &	operator= (const G4CashKarpRKF45 &)=delete
void	Stepper (const G4double y[], const G4double dydx[], G4double h, G4double yout[], G4double yerr[]) override
G4double	DistChord () const override
G4int	IntegratorOrder () const override
G4StepperType	StepperType () const override
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
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
G4int	IntegrationOrder ()
G4EquationOfMotion *	GetEquationOfMotion ()
const G4EquationOfMotion *	GetEquationOfMotion () const
void	SetEquationOfMotion (G4EquationOfMotion *newEquation)
unsigned long	GetfNoRHSCalls ()
void	ResetfNORHSCalls ()
G4bool	IsFSAL () const
G4bool	isQSS () const
void	SetIsQSS (G4bool val)

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

Detailed Description

G4CashKarpRKF45 implements the Cash-Karp Runge-Kutta-Fehlberg 4/5 method, an embedded fourth order method (giving fifth-order accuracy) for the solution of an ODE. Two different fourth order estimates are calculated; their difference gives an error estimate. It is used to integrate the equations of the motion of a particle in a magnetic field.

Definition at line 53 of file G4CashKarpRKF45.hh.

Constructor & Destructor Documentation

G4CashKarpRKF45() [1/2]

G4CashKarpRKF45::G4CashKarpRKF45	(	G4EquationOfMotion *	EqRhs,
		G4int	numberOfVariables = 6,
		G4bool	primary = true )

Constructor for G4CashKarpRKF45.

Parameters

[in]	EqRhs	Pointer to the provided equation of motion.
[in]	numberOfVariables	The number of integration variables.
[in]	primary	Flag for initialisation of the auxiliary stepper.

Definition at line 47 of file G4CashKarpRKF45.cc.

  : G4MagIntegratorStepper(EqRhs, noIntegrationVariables)
{
  const G4int numberOfVariables =
      std::max( noIntegrationVariables,
               ( ( (noIntegrationVariables-1)/4 + 1 ) * 4 ) );
  // For better alignment with cache-line
  
  ak2 = new G4double[numberOfVariables] ;
  ak3 = new G4double[numberOfVariables] ;
  ak4 = new G4double[numberOfVariables] ;
  ak5 = new G4double[numberOfVariables] ;
  ak6 = new G4double[numberOfVariables] ;
 
  // Must ensure space extra 'state' variables exists - i.e. yIn[7]
  const G4int numStateMax  = std::max(GetNumberOfStateVariables(), 8);  
  const G4int numStateVars = std::max(noIntegrationVariables,
                                      numStateMax );
                                   // GetNumberOfStateVariables() ); 
                                      
  yTemp = new G4double[numStateVars] ;
  yIn = new G4double[numStateVars] ;
 
  fLastInitialVector = new G4double[numStateVars] ;
  fLastFinalVector = new G4double[numStateVars] ;
  fLastDyDx = new G4double[numberOfVariables];
 
  fMidVector = new G4double[numStateVars];
  fMidError =  new G4double[numStateVars];
  if( primary )
  { 
    fAuxStepper = new G4CashKarpRKF45(EqRhs, numberOfVariables, !primary);
  }
}

Referenced by G4CashKarpRKF45(), G4CashKarpRKF45(), and operator=().

~G4CashKarpRKF45()

G4CashKarpRKF45::~G4CashKarpRKF45 ( )

override

Destructor.

Definition at line 88 of file G4CashKarpRKF45.cc.

{
  delete [] ak2;
  delete [] ak3;
  delete [] ak4;
  delete [] ak5;
  delete [] ak6;
 
  delete [] yTemp;
  delete [] yIn;
 
  delete [] fLastInitialVector;
  delete [] fLastFinalVector;
  delete [] fLastDyDx;
  delete [] fMidVector;
  delete [] fMidError; 
 
  delete fAuxStepper;
}

G4CashKarpRKF45() [2/2]

G4CashKarpRKF45::G4CashKarpRKF45 ( const G4CashKarpRKF45 & )

delete

Copy constructor and assignment operator not allowed.

Member Function Documentation

DistChord()

G4double G4CashKarpRKF45::DistChord ( ) const

overridevirtual

Returns the distance from chord line.

Implements G4MagIntegratorStepper.

Definition at line 220 of file G4CashKarpRKF45.cc.

{
  G4double distLine, distChord; 
  G4ThreeVector initialPoint, finalPoint, midPoint;
 
  // Store last initial and final points
  // (they will be overwritten in self-Stepper call!)
  //
  initialPoint = G4ThreeVector( fLastInitialVector[0], 
                                fLastInitialVector[1], fLastInitialVector[2]); 
  finalPoint   = G4ThreeVector( fLastFinalVector[0],  
                                fLastFinalVector[1],  fLastFinalVector[2]); 
 
  // Do half a step using StepNoErr
  //
  fAuxStepper->Stepper( fLastInitialVector, fLastDyDx,
                        0.5 * fLastStepLength, fMidVector, fMidError );
 
  midPoint = G4ThreeVector( fMidVector[0], fMidVector[1], fMidVector[2]);       
 
  // Use stored values of Initial and Endpoint + new Midpoint to evaluate
  // distance of Chord
  //
  if (initialPoint != finalPoint) 
  {
     distLine  = G4LineSection::Distline( midPoint, initialPoint, finalPoint );
     distChord = distLine;
  }
  else
  {
     distChord = (midPoint-initialPoint).mag();
  }
  return distChord;
}

IntegratorOrder()

G4int G4CashKarpRKF45::IntegratorOrder ( ) const

inlineoverridevirtual

Returns the order, 4, of integration.

Implements G4MagIntegratorStepper.

Definition at line 103 of file G4CashKarpRKF45.hh.

{ return 4; }

operator=()

G4CashKarpRKF45 & G4CashKarpRKF45::operator= ( const G4CashKarpRKF45 & )

delete

Stepper()

void G4CashKarpRKF45::Stepper ( const G4double y[], const G4double dydx[], G4double h, G4double yout[], G4double yerr[] )

overridevirtual

The stepper for the Runge Kutta integration. The stepsize is fixed, with the step size given by 'h'. Integrates ODE starting values y[0 to 6]. Outputs yout[] and its estimated error yerr[].

Parameters

[in]	y	Starting values array of integration variables.
[in]	dydx	Derivatives array.
[in]	h	The given step size.
[out]	yout	Integration output.
[out]	yerr	The estimated error.

Implements G4MagIntegratorStepper.

Definition at line 119 of file G4CashKarpRKF45.cc.

{
  // const G4int nvar = 6 ;
  // const G4double a2 = 0.2 , a3 = 0.3 , a4 = 0.6 , a5 = 1.0 , a6 = 0.875;
  G4int i;
 
  const G4double  b21 = 0.2 ,
                  b31 = 3.0/40.0 , b32 = 9.0/40.0 ,
                  b41 = 0.3 , b42 = -0.9 , b43 = 1.2 ,
 
                  b51 = -11.0/54.0 , b52 = 2.5 , b53 = -70.0/27.0 ,
                  b54 = 35.0/27.0 ,
 
                  b61 = 1631.0/55296.0 , b62 =   175.0/512.0 ,
                  b63 =  575.0/13824.0 , b64 = 44275.0/110592.0 ,
                  b65 =  253.0/4096.0 ,
 
                  c1 = 37.0/378.0 , c3 = 250.0/621.0 , c4 = 125.0/594.0 ,
                  c6 = 512.0/1771.0 , dc5 = -277.0/14336.0 ;
 
  const G4double dc1 = c1 - 2825.0/27648.0 ,  dc3 = c3 - 18575.0/48384.0 ,
                 dc4 = c4 - 13525.0/55296.0 , dc6 = c6 - 0.25 ;
 
  // Initialise time to t0, needed when it is not updated by the integration.
  //        [ Note: Only for time dependent fields (usually electric) 
  //                  is it neccessary to integrate the time.] 
  yOut[7] = yTemp[7] = yIn[7] = yInput[7]; 
 
  const G4int numberOfVariables= this->GetNumberOfVariables(); 
    // The number of variables to be integrated over
 
  //  Saving yInput because yInput and yOut can be aliases for same array
 
  for(i=0; i<numberOfVariables; ++i) 
  {
    yIn[i]=yInput[i];
  }
  // RightHandSide(yIn, dydx) ;            // 1st Step
 
  for(i=0; i<numberOfVariables; ++i) 
  {
    yTemp[i] = yIn[i] + b21*Step*dydx[i] ;
  }
  RightHandSide(yTemp, ak2) ;              // 2nd Step
 
  for(i=0; i<numberOfVariables; ++i)
  {
    yTemp[i] = yIn[i] + Step*(b31*dydx[i] + b32*ak2[i]) ;
  }
  RightHandSide(yTemp, ak3) ;              // 3rd Step
 
  for(i=0; i<numberOfVariables; ++i)
  {
    yTemp[i] = yIn[i] + Step*(b41*dydx[i] + b42*ak2[i] + b43*ak3[i]) ;
  }
  RightHandSide(yTemp, ak4) ;              // 4th Step
 
  for(i=0; i<numberOfVariables; ++i)
  {
    yTemp[i] = yIn[i] + Step*(b51*dydx[i]
                      + b52*ak2[i] + b53*ak3[i] + b54*ak4[i]) ;
  }
  RightHandSide(yTemp, ak5) ;              // 5th Step
 
  for(i=0; i<numberOfVariables; ++i)
  {
    yTemp[i] = yIn[i] + Step*(b61*dydx[i]
                      + b62*ak2[i] + b63*ak3[i] + b64*ak4[i] + b65*ak5[i]) ;
  }
  RightHandSide(yTemp, ak6) ;              // 6th Step
 
  for(i=0; i<numberOfVariables; ++i)
  {
    // Accumulate increments with proper weights
    //
    yOut[i] = yIn[i] + Step*(c1*dydx[i] + c3*ak3[i] + c4*ak4[i] + c6*ak6[i]) ;
 
    // Estimate error as difference between 4th and 5th order methods
    //
    yErr[i] = Step*(dc1*dydx[i]
            + dc3*ak3[i] + dc4*ak4[i] + dc5*ak5[i] + dc6*ak6[i]) ;
 
    // Store Input and Final values, for possible use in calculating chord
    //
    fLastInitialVector[i] = yIn[i] ;
    fLastFinalVector[i]   = yOut[i];
    fLastDyDx[i]          = dydx[i];
  }
  // NormaliseTangentVector( yOut ); // Not wanted
 
  fLastStepLength = Step;
 
  return;
} 

StepperType()

G4StepperType G4CashKarpRKF45::StepperType ( ) const

inlineoverridevirtual

Returns the stepper type-ID, "kCashKarpRKF45".

Reimplemented from G4MagIntegratorStepper.

Definition at line 108 of file G4CashKarpRKF45.hh.

{ return kCashKarpRKF45; }

---

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

• G4CashKarpRKF45.hh
• G4CashKarpRKF45.cc