G4FSALDormandPrince745 Class Reference

#include <G4FSALDormandPrince745.hh>

Inheritance diagram for G4FSALDormandPrince745:

Public Member Functions
	G4FSALDormandPrince745 (G4EquationOfMotion *EqRhs, G4int numberOfVariables=6, G4bool primary=true)
	~G4FSALDormandPrince745 ()
	G4FSALDormandPrince745 (const G4FSALDormandPrince745 &)=delete
G4FSALDormandPrince745 &	operator= (const G4FSALDormandPrince745 &)=delete
void	Stepper (const G4double y[], const G4double dydx[], G4double h, G4double yout[], G4double yerr[], G4double nextDydx[])
void	interpolate (const G4double yInput[], const G4double dydx[], G4double yOut[], G4double Step, G4double tau)
void	SetupInterpolate (const G4double yInput[], const G4double dydx[], const G4double Step)
void	Interpolate (const G4double yInput[], const G4double dydx[], const G4double Step, G4double yOut[], G4double tau)
G4double	DistChord () const
G4int	IntegratorOrder () const
G4bool	isFSAL () const
Public Member Functions inherited from G4VFSALIntegrationStepper
	G4VFSALIntegrationStepper (G4EquationOfMotion *Equation, G4int numIntegrationVariables, G4int numStateVariables=12)
virtual	~G4VFSALIntegrationStepper ()=default
	G4VFSALIntegrationStepper (const G4VFSALIntegrationStepper &)=delete
G4VFSALIntegrationStepper &	operator= (const G4VFSALIntegrationStepper &)=delete
virtual void	Stepper (const G4double y[], const G4double dydx[], G4double h, G4double yout[], G4double yerr[], G4double lastDydx[])=0
virtual G4double	DistChord () const =0
void	NormaliseTangentVector (G4double vec[6])
void	NormalisePolarizationVector (G4double vec[12])
void	RightHandSide (const double y[], double dydx[])
G4int	GetNumberOfVariables () const
G4int	GetNumberOfStateVariables () const
virtual G4int	IntegratorOrder () const =0
G4EquationOfMotion *	GetEquationOfMotion ()
void	SetEquationOfMotion (G4EquationOfMotion *newEquation)
G4int	GetfNoRHSCalls ()
void	increasefNORHSCalls ()
void	ResetfNORHSCalls ()

Detailed Description

Definition at line 40 of file G4FSALDormandPrince745.hh.

Constructor & Destructor Documentation

G4FSALDormandPrince745() [1/2]

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

Definition at line 51 of file G4FSALDormandPrince745.cc.

  : G4VFSALIntegrationStepper(EqRhs, noIntegrationVariables)
{
  const G4int numberOfVariables = noIntegrationVariables;
    
  // New Chunk of memory being created for use by the stepper
    
  // aki - for storing intermediate RHS
  //
  ak2 = new G4double[numberOfVariables];
  ak3 = new G4double[numberOfVariables];
  ak4 = new G4double[numberOfVariables];
  ak5 = new G4double[numberOfVariables];
  ak6 = new G4double[numberOfVariables];
  ak7 = new G4double[numberOfVariables];
 
  // Also always allocate arrays for interpolation stages    
  //
  ak8 = new G4double[numberOfVariables];
  ak9 = new G4double[numberOfVariables];
    
  yTemp = new G4double[numberOfVariables] ;
  yIn = new G4double[numberOfVariables] ;
    
  pseudoDydx_for_DistChord = new G4double[numberOfVariables];
 
  fInitialDyDx = new G4double[numberOfVariables];    
  fLastInitialVector = new G4double[numberOfVariables] ;
  fLastFinalVector = new G4double[numberOfVariables] ;
  fLastDyDx = new G4double[numberOfVariables];
    
  fMidVector = new G4double[numberOfVariables];
  fMidError =  new G4double[numberOfVariables];
    
  if( primary )
  {
    fAuxStepper = new G4FSALDormandPrince745(EqRhs,numberOfVariables,!primary);
  }
}

~G4FSALDormandPrince745()

G4FSALDormandPrince745::~G4FSALDormandPrince745

(

)

Definition at line 95 of file G4FSALDormandPrince745.cc.

{
    // Clear all previously allocated memory for stepper and DistChord
 
    delete [] ak2;  ak2 = nullptr;
    delete [] ak3;  ak3 = nullptr;
    delete [] ak4;  ak4 = nullptr;
    delete [] ak5;  ak5 = nullptr;
    delete [] ak6;  ak6 = nullptr;
    delete [] ak7;  ak7 = nullptr;
    delete [] ak8;  ak8 = nullptr;
    delete [] ak9;  ak9 = nullptr;
    
    delete [] yTemp; yTemp = nullptr;
    delete [] yIn;   yIn = nullptr;
 
    delete [] pseudoDydx_for_DistChord;  pseudoDydx_for_DistChord = nullptr;
    delete [] fInitialDyDx;              fInitialDyDx = nullptr;
    
    delete [] fLastInitialVector;    fLastInitialVector = nullptr;
    delete [] fLastFinalVector;      fLastFinalVector  = nullptr;
    delete [] fLastDyDx;             fLastDyDx  = nullptr;
    delete [] fMidVector;            fMidVector = nullptr;
    delete [] fMidError;             fMidError  = nullptr;
    
    delete fAuxStepper; fAuxStepper = nullptr;
}

G4FSALDormandPrince745() [2/2]

G4FSALDormandPrince745::G4FSALDormandPrince745 ( const G4FSALDormandPrince745 &  )

delete

Member Function Documentation

DistChord()

G4double G4FSALDormandPrince745::DistChord ( ) const

virtual

Implements G4VFSALIntegrationStepper.

Definition at line 243 of file G4FSALDormandPrince745.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, pseudoDydx_for_DistChord );
    
    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 G4FSALDormandPrince745::IntegratorOrder ( ) const

inlinevirtual

Implements G4VFSALIntegrationStepper.

Definition at line 78 of file G4FSALDormandPrince745.hh.

{return 4; }

Interpolate()

void G4FSALDormandPrince745::Interpolate	(	const G4double	yInput[],
		const G4double	dydx[],
		const G4double	Step,
		G4double	yOut[],
		G4double	tau
	)

Definition at line 384 of file G4FSALDormandPrince745.cc.

{
    // Define the coefficients for the polynomials
 
    G4double bi[10][5], b[10];
    G4int numberOfVariables = GetNumberOfVariables();
    
    //  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 ;
    //  --------------------------------------------------------
 
    for(G4int i = 0; i< numberOfVariables; ++i)
    {
        yIn[i] = yInput[i];
    }
    
    G4double tau0 = tau;
 
    // Calculating the polynomials
    //    
    for(auto i=1; i<=9; ++i)  // i is NOT the coordinate no., it's stage no.
    {
        b[i] = 0;
        tau = 1.0;
        for(auto j=0; j<=4; ++j)
        {
            b[i] += bi[i][j]*tau;
            tau*=tau0;
        }
    }
    
    for(G4int i=0; i<numberOfVariables; ++i) // Here i IS the coordinate no.
    {
        yOut[i] = yIn[i] + Step*tau0*(b[1]*dydx[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] );
    }
}

interpolate()

void G4FSALDormandPrince745::interpolate	(	const G4double	yInput[],
		const G4double	dydx[],
		G4double	yOut[],
		G4double	Step,
		G4double	tau
	)

Definition at line 280 of file G4FSALDormandPrince745.cc.

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

isFSAL()

G4bool G4FSALDormandPrince745::isFSAL ( ) const

inline

Definition at line 79 of file G4FSALDormandPrince745.hh.

{ return true; }

operator=()

G4FSALDormandPrince745 & G4FSALDormandPrince745::operator= ( const G4FSALDormandPrince745 &  )

delete

SetupInterpolate()

void G4FSALDormandPrince745::SetupInterpolate	(	const G4double	yInput[],
		const G4double	dydx[],
		const G4double	Step
	)

Definition at line 329 of file G4FSALDormandPrince745.cc.

{    
    // Coefficients for the additional stages
    //
    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();
    
    // Saving yInput because yInput and yOut can be aliases for same array
    //
    for(G4int i=0; i<numberOfVariables; ++i)
    {
      yIn[i] = yInput[i];
    }
    
    yTemp[7] = yIn[7];
    
    // Evaluate the extra stages
    //
    for(G4int i=0; i<numberOfVariables; ++i)
    {
        yTemp[i] = yIn[i] + Step*( b81*dydx[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] = yIn[i] + Step * ( b91*dydx[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()

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

virtual

Implements G4VFSALIntegrationStepper.

Definition at line 128 of file G4FSALDormandPrince745.cc.

{
    G4int i;
    
    // 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,
    
             //    c1 = 35.0/384.0, c2 = .0,
             //    c3 = 500.0/1113.0, c4 = 125.0/192.0,
             //    c5 = -2187.0/6784.0, c6 = 11.0/84.0,
             //    c7 = 0,
    
                   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 ;  //end of declaration
 
    const G4int numberOfVariables = 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];
        fInitialDyDx[i] = dydx[i];
    }
    // Ensure that time is initialised - in case it is not integrated
    //
    yOut[7] = yTemp[7]  = yInput[7];
    // RightHandSide(yIn, DyDx) ;   // 1st Step - Not doing, getting passed
    
    for(i=0; i<numberOfVariables; ++i)
    {
        yTemp[i] = yIn[i] + b21*Step*fInitialDyDx[i] ;
    }
    RightHandSide(yTemp, ak2) ;              // 2nd Step
    
    for(i=0; i<numberOfVariables; ++i)
    {
        yTemp[i] = yIn[i] + Step*(b31*fInitialDyDx[i] + b32*ak2[i]) ;
    }
    RightHandSide(yTemp, ak3) ;              // 3rd Step
    
    for(i=0; i<numberOfVariables; ++i)
    {
        yTemp[i] = yIn[i] + Step*(b41*fInitialDyDx[i]
                                + b42*ak2[i] + b43*ak3[i]) ;
    }
    RightHandSide(yTemp, ak4) ;              // 4th Step
    
    for(i=0; i<numberOfVariables; ++i)
    {
        yTemp[i] = yIn[i] + Step*(b51*fInitialDyDx[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*fInitialDyDx[i] + b62*ak2[i]
                                + b63*ak3[i] + b64*ak4[i] + b65*ak5[i]) ;
    }
    RightHandSide(yTemp, ak6) ;              // 6th Step
    
    for(i=0; i<numberOfVariables; ++i)
    {
        yOut[i] = yIn[i] + Step*(b71*fInitialDyDx[i] + b72*ak2[i] + b73*ak3[i]
                               + b74*ak4[i] + b75*ak5[i] + b76*ak6[i]);
    }
    RightHandSide(yOut, ak7);               //7th and Final step
    
    for(i=0; i<numberOfVariables; ++i)
    {
        
        yErr[i] = Step*(dc1*fInitialDyDx[i] + dc2*ak2[i] + dc3*ak3[i]
                      + dc4*ak4[i] + dc5*ak5[i] + dc6*ak6[i] + dc7*ak7[i] ) ;
 
        // Store Input and Final values, for possible use in calculating chord
        //
        fLastInitialVector[i] = yIn[i] ;
        fLastFinalVector[i]   = yOut[i];
        fLastDyDx[i]          = fInitialDyDx[i];
        nextDydx[i] = ak7[i];
    }
    fLastStepLength = Step;
    
    return ;
}

Referenced by DistChord().

---

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

• G4FSALDormandPrince745.hh
• G4FSALDormandPrince745.cc