G4TCashKarpRKF45< T_Equation, N > Class Template Reference

G4TCashKarpRKF45 is a templated version of Cash-Karp 4th/5th order embedded stepper. More...

#include <G4TCashKarpRKF45.hh>

Inheritance diagram for G4TCashKarpRKF45< T_Equation, N >:

Public Member Functions
	G4TCashKarpRKF45 (T_Equation *EqRhs, G4bool primary=true)
virtual	~G4TCashKarpRKF45 ()
	G4TCashKarpRKF45 (const G4TCashKarpRKF45 &)=delete
G4TCashKarpRKF45 &	operator= (const G4TCashKarpRKF45 &)=delete
void	StepWithError (const G4double yInput[], const G4double dydx[], G4double Step, G4double yOut[], G4double yErr[])
virtual void	Stepper (const G4double yInput[], const G4double dydx[], G4double hstep, G4double yOutput[], G4double yError[]) override final
void	RightHandSideInl (const G4double y[], G4double dydx[])
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

template<class T_Equation, unsigned int N = 6>
class G4TCashKarpRKF45< T_Equation, N >

G4TCashKarpRKF45 is a templated version of Cash-Karp 4th/5th order embedded stepper.

Definition at line 67 of file G4TCashKarpRKF45.hh.

Constructor & Destructor Documentation

G4TCashKarpRKF45() [1/2]

template<class T_Equation, unsigned int N>

G4TCashKarpRKF45< T_Equation, N >::G4TCashKarpRKF45	(	T_Equation *	EqRhs,
		G4bool	primary = true )

Definition at line 128 of file G4TCashKarpRKF45.hh.

    : G4MagIntegratorStepper(dynamic_cast<G4EquationOfMotion*>(EqRhs), N )
    , fEquation_Rhs(EqRhs)
{
  if( dynamic_cast<G4EquationOfMotion*>(EqRhs) == nullptr )
  {
    G4Exception("G4TCashKarpRKF45: constructor", "GeomField0001",
                FatalException, "Equation is not an G4EquationOfMotion.");      
  }
    
  fLastInitialVector = new G4double[N];
  fLastFinalVector   = new G4double[N];
  fLastDyDx          = new G4double[N];
  
  fMidVector = new G4double[N];
  fMidError  = new G4double[N];
  
  if(primary)
  {
     fAuxStepper = new G4TCashKarpRKF45<T_Equation, N> (EqRhs, !primary);
  }
}

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

~G4TCashKarpRKF45()

template<class T_Equation, unsigned int N>

G4TCashKarpRKF45< T_Equation, N >::~G4TCashKarpRKF45 ( )

virtual

Definition at line 153 of file G4TCashKarpRKF45.hh.

{
  delete[] fLastInitialVector;
  delete[] fLastFinalVector;
  delete[] fLastDyDx;
  delete[] fMidVector;
  delete[] fMidError;
  
  delete fAuxStepper;
}

G4TCashKarpRKF45() [2/2]

template<class T_Equation, unsigned int N = 6>

G4TCashKarpRKF45< T_Equation, N >::G4TCashKarpRKF45 ( const G4TCashKarpRKF45< T_Equation, N > & )

delete

Member Function Documentation

DistChord()

template<class T_Equation, unsigned int N>

G4double G4TCashKarpRKF45< T_Equation, N >::DistChord ( ) const

inlineoverridevirtual

Estimates the maximum distance of a chord from the true path over the segment last integrated.

Implements G4MagIntegratorStepper.

Definition at line 275 of file G4TCashKarpRKF45.hh.

{
  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->G4TCashKarpRKF45::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()

template<class T_Equation, unsigned int N = 6>

G4int G4TCashKarpRKF45< T_Equation, N >::IntegratorOrder ( ) const

inlineoverridevirtual

Returns the order of the integrator, i.e. its error behaviour is of the order O(h^order).

Implements G4MagIntegratorStepper.

Definition at line 101 of file G4TCashKarpRKF45.hh.

{ return 4; }

operator=()

template<class T_Equation, unsigned int N = 6>

G4TCashKarpRKF45 & G4TCashKarpRKF45< T_Equation, N >::operator= ( const G4TCashKarpRKF45< T_Equation, N > & )

delete

RightHandSideInl()

template<class T_Equation, unsigned int N = 6>

void G4TCashKarpRKF45< T_Equation, N >::RightHandSideInl ( const G4double y[], G4double dydx[] )

inline

Definition at line 93 of file G4TCashKarpRKF45.hh.

    {
      fEquation_Rhs->T_Equation::RightHandSide(y, dydx);
    }

Referenced by StepWithError().

Stepper()

template<class T_Equation, unsigned int N>

void G4TCashKarpRKF45< T_Equation, N >::Stepper ( const G4double y[], const G4double dydx[], G4double h, G4double yout[], G4double yerr[] )

inlinefinaloverridevirtual

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 312 of file G4TCashKarpRKF45.hh.

{
  assert( yOutput != yInput );
  assert( yError  != yInput );
  
  StepWithError( yInput, dydx, Step, yOutput, yError);
}

StepperType()

template<class T_Equation, unsigned int N = 6>

G4StepperType G4TCashKarpRKF45< T_Equation, N >::StepperType ( ) const

inlineoverridevirtual

Returns the stepper type ID ('kUserStepper'). This function should be overriden in derived classes.

Reimplemented from G4MagIntegratorStepper.

Definition at line 103 of file G4TCashKarpRKF45.hh.

{ return kTCashKarpRKF45; }

StepWithError()

template<class T_Equation, unsigned int N = 6>

void G4TCashKarpRKF45< T_Equation, N >::StepWithError ( const G4double yInput[], const G4double dydx[], G4double Step, G4double yOut[], G4double yErr[] )

inline

Definition at line 176 of file G4TCashKarpRKF45.hh.

{
  // const G4double a2 = 0.2 , a3 = 0.3 , a4 = 0.6 , a5 = 1.0 , a6 = 0.875;
 
  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];
 
  //  Saving yInput because yInput and yOut can be aliases for same array
  for(unsigned int i = 0; i < N; ++i)
  {
    yIn[i] = yInput[i];
  }
  // RightHandSideInl(yIn, dydx) ;              // 1st Step
 
  for(unsigned int i = 0; i < N; ++i)
  {
    yTemp[i] = yIn[i] + b21 * Step * dydx[i];
  }
  this->RightHandSideInl(yTemp, ak2);  // 2nd Step
 
  for(unsigned int i = 0; i < N; ++i)
  {
    yTemp[i] = yIn[i] + Step * (b31 * dydx[i] + b32 * ak2[i]);
  }
  this->RightHandSideInl(yTemp, ak3);  // 3rd Step
  
  for(unsigned int i = 0; i < N; ++i)
  {
    yTemp[i] = yIn[i] + Step * (b41 * dydx[i] + b42 * ak2[i] + b43 * ak3[i]);
  }
  this->RightHandSideInl(yTemp, ak4);  // 4th Step
  
  for(unsigned int i = 0; i < N; ++i)
  {
    yTemp[i] = yIn[i] + Step * (b51 * dydx[i] + b52 * ak2[i] + b53 * ak3[i] +
                                b54 * ak4[i]);
  }
  this->RightHandSideInl(yTemp, ak5);  // 5th Step
  
  for(unsigned int i = 0; i < N; ++i)
  {
    yTemp[i] = yIn[i] + Step * (b61 * dydx[i] + b62 * ak2[i] + b63 * ak3[i] +
                                b64 * ak4[i] + b65 * ak5[i]);
  }
  this->RightHandSideInl(yTemp, ak6);  // 6th Step
 
  for(unsigned int i = 0; i < N; ++i)
  {
    // Accumulate increments with proper weights
     
    yOut[i] = yIn[i] +
       Step * (c1 * dydx[i] + c3 * ak3[i] + c4 * ak4[i] + c6 * ak6[i]);
  }
  for(unsigned int i = 0; i < N; ++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]);
  }
  for(unsigned int i = 0; i < N; ++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;
}

Referenced by Stepper().

---

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

• G4TCashKarpRKF45.hh