G4FSALBogackiShampine45 Class Reference

G4FSALBogackiShampine45 is an integrator of particle's equation of motion based on the Bogacki-Shampine - 8 - 5(4) FSAL implementation. More...

#include <G4FSALBogackiShampine45.hh>

Inheritance diagram for G4FSALBogackiShampine45:

Public Member Functions
	G4FSALBogackiShampine45 (G4EquationOfMotion *EqRhs, G4int numberOfVariables=6, G4bool primary=true)
	~G4FSALBogackiShampine45 () override
	G4FSALBogackiShampine45 (const G4FSALBogackiShampine45 &)=delete
G4FSALBogackiShampine45 &	operator= (const G4FSALBogackiShampine45 &)=delete
void	Stepper (const G4double y[], const G4double dydx[], G4double h, G4double yout[], G4double yerr[], G4double nextDydx[]) override
void	interpolate (const G4double yInput[], const G4double dydx[], G4double yOut[], G4double Step, G4double tau)
G4double	DistChord () const override
G4int	IntegratorOrder () const override
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
void	NormaliseTangentVector (G4double vec[6])
void	NormalisePolarizationVector (G4double vec[12])
void	RightHandSide (const double y[], double dydx[])
G4int	GetNumberOfVariables () const
G4int	GetNumberOfStateVariables () const
G4EquationOfMotion *	GetEquationOfMotion ()
void	SetEquationOfMotion (G4EquationOfMotion *newEquation)
G4int	GetfNoRHSCalls ()
void	increasefNORHSCalls ()
void	ResetfNORHSCalls ()

Detailed Description

G4FSALBogackiShampine45 is an integrator of particle's equation of motion based on the Bogacki-Shampine - 8 - 5(4) FSAL implementation.

Definition at line 45 of file G4FSALBogackiShampine45.hh.

Constructor & Destructor Documentation

G4FSALBogackiShampine45() [1/2]

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

Constructor for G4FSALBogackiShampine45.

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 59 of file G4FSALBogackiShampine45.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];
    ak8 = new G4double[numberOfVariables];
 
    ak9 = new G4double[numberOfVariables];
    ak10 = new G4double[numberOfVariables];
    ak11 = new G4double[numberOfVariables];
    DyDx = new G4double[numberOfVariables];
    
    assert ( GetNumberOfStateVariables() >= 8 );
    const G4int numStateVars = std::max(noIntegrationVariables,
                                        GetNumberOfStateVariables() );  
 
    // Must ensure space extra 'state' variables exists - i.e. yIn[7]
    //
    yTemp = new G4double[numStateVars];
    yIn = new G4double[numStateVars] ;
    
    fLastInitialVector = new G4double[numStateVars] ;
    fLastFinalVector = new G4double[numStateVars] ;
    fLastDyDx = new G4double[numberOfVariables];  // Only derivatives
    
    fMidVector = new G4double[numStateVars];
    fMidError =  new G4double[numStateVars];
    
    pseudoDydx_for_DistChord = new G4double[numberOfVariables];
    
    fMidVector = new G4double[numberOfVariables];
    fMidError =  new G4double[numberOfVariables];
    if( primary )
    {
      fAuxStepper = new G4FSALBogackiShampine45(EqRhs, numberOfVariables,
                                                !primary);
    }
    if( !fPreparedConstants )
    {
       PrepareConstants();
    }
}

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

~G4FSALBogackiShampine45()

G4FSALBogackiShampine45::~G4FSALBogackiShampine45 ( )

override

Destructor.

Definition at line 116 of file G4FSALBogackiShampine45.cc.

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

G4FSALBogackiShampine45() [2/2]

G4FSALBogackiShampine45::G4FSALBogackiShampine45 ( const G4FSALBogackiShampine45 & )

delete

Copy constructor and assignment operator not allowed.

Member Function Documentation

DistChord()

G4double G4FSALBogackiShampine45::DistChord ( ) const

overridevirtual

Returns the distance from chord line.

Implements G4VFSALIntegrationStepper.

Definition at line 293 of file G4FSALBogackiShampine45.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 G4FSALBogackiShampine45::IntegratorOrder ( ) const

inlineoverridevirtual

Returns the order, 4, of integration.

Implements G4VFSALIntegrationStepper.

Definition at line 111 of file G4FSALBogackiShampine45.hh.

{ return 4; }

interpolate()

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

Calculates the output at the tau fraction of step.

Parameters

[in]	yInput	Starting values array of integration variables.
[in]	dydx	Derivatives array.
[out]	yOut	Interpolation output.
[in]	Step	The given step size.
[in]	tau	The tau fraction of the step.

Definition at line 410 of file G4FSALBogackiShampine45.cc.

{
    const G4double a91 = 455.0/6144.0 ,
                   a92 = 0.0 ,
                   a93 = 10256301.0/35409920.0 ,
                   a94 = 2307361.0/17971200.0 ,
                   a95 = -387.0/102400.0 ,
                   a96 = 73.0/5130.0 ,
                   a97 = -7267.0/215040.0 ,
                   a98 = 1.0/32.0 ,
 
                   a101 = -837888343715.0/13176988637184.0 ,
                   a102 = 30409415.0/52955362.0 ,
                   a103 = -48321525963.0/759168069632.0 ,
                   a104 = 8530738453321.0/197654829557760.0 ,
                   a105 = 1361640523001.0/1626788720640.0 ,
                   a106 = -13143060689.0/38604458898.0 ,
                   a107 = 18700221969.0/379584034816.0 ,
                   a108 = -5831595.0/847285792.0 ,
                   a109 = -5183640.0/26477681.0 ,
 
                   a111 = 98719073263.0/1551965184000.0 ,
                   a112 = 1307.0/123552.0 ,
                   a113 = 4632066559387.0/70181753241600.0 ,
                   a114 = 7828594302389.0/382182512025600.0 ,
                   a115 = 40763687.0/11070259200.0 ,
                   a116 = 34872732407.0/224610586200.0 ,
                   a117 = -2561897.0/30105600.0 ,
                   a118 =  1.0/10.0 ,
                   a119 = -1.0/10.0 ,
                   a1110 = -1403317093.0/11371610250.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];
    }
    
    // The number of variables to be integrated over
    //
    yOut[7] = yTemp[7]  = yIn[7];
 
    // Calculating extra stages
    //
    for(G4int i=0; i<numberOfVariables; ++i)
    {
        yTemp[i] = yIn[i] + Step*(a91*dydx[i] + a92*ak2[i] + a93*ak3[i] +
                                  a94*ak4[i] + a95*ak5[i] + a96*ak6[i] +
                                  a97*ak7[i] + a98*ak8[i] );
    }
    
    RightHandSide(yTemp, ak9);
    
    for(G4int i=0; i<numberOfVariables; ++i)
    {
        yTemp[i] = yIn[i] + Step*(a101*dydx[i] + a102*ak2[i] + a103*ak3[i] +
                                  a104*ak4[i] + a105*ak5[i] + a106*ak6[i] +
                                  a107*ak7[i] + a108*ak8[i] + a109*ak9[i] );
    }
    
    RightHandSide(yTemp, ak10);
 
    for(G4int i=0; i<numberOfVariables; ++i)
    {
        yTemp[i] = yIn[i] + Step*(a111*dydx[i] + a112*ak2[i] + a113*ak3[i] +
                                  a114*ak4[i] + a115*ak5[i] + a116*ak6[i] +
                                  a117*ak7[i] + a118*ak8[i] + a119*ak9[i] +
                                  a1110*ak10[i] );
    }
    
    RightHandSide(yTemp, ak11);
    
    G4double tau0 = tau;
 
    // Calculating the polynomials
    //
    for(auto i=1; i<=11; ++i) // i is NOT the coordinate no., it's stage no.
    {
        b[i] = 0.0;
        tau = tau0;
        for(auto j=1; j<=6; ++j)
        {
            b[i] += bi[i][j]*tau;
            tau*=tau0;
        }
    }
    
    for(G4int i=0; i<numberOfVariables; ++i)
    {
        yOut[i] = yIn[i] + Step*(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] +
                                 b[10]*ak10[i] + b[11]*ak11[i] );
    }  
}

operator=()

G4FSALBogackiShampine45 & G4FSALBogackiShampine45::operator= ( const G4FSALBogackiShampine45 & )

delete

Stepper()

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

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.
[out]	nextDydx	Last derivatives array for the next step.

Implements G4VFSALIntegrationStepper.

Definition at line 150 of file G4FSALBogackiShampine45.cc.

{
    G4int i;
    
    // The various constants defined on the basis of butcher tableu
 
    const G4double b21 = 1.0/6.0 ,
                   b31 = 2.0/27.0 , b32 = 4.0/27.0,
    
                   b41 = 183.0/1372.0 , b42 = -162.0/343.0, b43 = 1053.0/1372.0,
    
                   b51 = 68.0/297.0, b52 = -4.0/11.0,
                   b53 = 42.0/143.0, b54 = 1960.0/3861.0,
    
                   b61 = 597.0/22528.0, b62 = 81.0/352.0,
                   b63 = 63099.0/585728.0, b64 = 58653.0/366080.0,
                   b65 = 4617.0/20480.0,
    
                   b71 = 174197.0/959244.0, b72 = -30942.0/79937.0,
                   b73 = 8152137.0/19744439.0, b74 = 666106.0/1039181.0,
                   b75 = -29421.0/29068.0,  b76 = 482048.0/414219.0,
    
                   b81 = 587.0/8064.0,  b82 = 0.0,
                   b83 = 4440339.0/15491840.0,  b84 = 24353.0/124800.0,
                   b85 = 387.0/44800.0, b86 = 2152.0/5985.0,
                   b87 = 7267.0/94080.0,
    
    
             //    c1 = 2479.0/34992.0,
             //    c2 = 0.0,
             //    c3 = 123.0/416.0,
             //    c4 = 612941.0/3411720.0,
             //    c5 = 43.0/1440.0,
             //    c6 = 2272.0/6561.0,
             //    c7 = 79937.0/1113912.0,
             //    c8 = 3293.0/556956.0,
    
    // For the embedded higher order method only the difference of values
    // taken and is used directly later instead of defining the last row
    // of butcher table in a separate set of variables and taking the
    // difference there
    
                   dc1 = b81 - 2479.0/34992.0 ,
                   dc2 = 0.0,
                   dc3 = b83 - 123.0/416.0 ,
                   dc4 = b84 - 612941.0/3411720.0,
                   dc5 = b85 - 43.0/1440.0,
                   dc6 = b86 - 2272.0/6561.0,
                   dc7 = b87 - 79937.0/1113912.0,
                   dc8 = -3293.0/556956.0;   // end of declaration
 
    const G4int numberOfVariables = GetNumberOfVariables();
    
    // The number of variables to be integrated over
    //
    yOut[7] = yTemp[7]  = yIn[7];
 
    //  Saving yInput because yInput and yOut can be aliases for same array
    //
    for(i=0; i<numberOfVariables; ++i)
    {
        yIn[i]=yInput[i];
        DyDx[i] = dydx[i];
    }
    // RightHandSide(yIn, dydx) ;   // 1st Step - Not doing, getting passed
    
    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)
    {
        yTemp[i] = yIn[i] + Step*(b71*DyDx[i] + b72*ak2[i] + b73*ak3[i] +
                                  b74*ak4[i] + b75*ak5[i] + b76*ak6[i]);
    }
    RightHandSide(yTemp, ak7);               // 7th Step
    
    for(i=0; i<numberOfVariables; ++i)
    {
        yOut[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(yOut, ak8);                // 8th Step - Final one Using FSAL
 
    
    for(i=0; i<numberOfVariables; ++i)
    {
        
        yErr[i] = Step*(dc1*DyDx[i] + dc2*ak2[i] + dc3*ak3[i] + dc4*ak4[i] +
                        dc5*ak5[i] + dc6*ak6[i] + dc7*ak7[i] + dc8*ak8[i]) ;
        
        
        // FSAL stepper : Must pass the last DyDx for the next step, here ak8
        //
        nextDydx[i] = ak8[i];
        
        // Store Input and Final values, for possible use in calculating chord
        //
        fLastInitialVector[i] = yIn[i] ;
        fLastFinalVector[i]   = yOut[i];
        fLastDyDx[i]          = DyDx[i];
    }
    fLastStepLength = Step;
    
    return;
}

---

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

• G4FSALBogackiShampine45.hh
• G4FSALBogackiShampine45.cc