Geant4-11
Public Member Functions | Protected Member Functions | Private Attributes | Static Private Attributes
G4BogackiShampine45 Class Reference

#include <G4BogackiShampine45.hh>

Inheritance diagram for G4BogackiShampine45:
G4MagIntegratorStepper

Public Member Functions

G4double DistChord () const
 
 G4BogackiShampine45 (const G4BogackiShampine45 &)=delete
 
 G4BogackiShampine45 (G4EquationOfMotion *EqRhs, G4int numberOfVariables=6, G4bool primary=true)
 
G4EquationOfMotion * GetEquationOfMotion ()
 
const G4EquationOfMotion * GetEquationOfMotion () const
 
unsigned long GetfNoRHSCalls ()
 
void GetLastDydx (G4double dyDxLast[])
 
G4int GetNumberOfStateVariables () const
 
G4int GetNumberOfVariables () const
 
G4int IntegrationOrder ()
 
G4int IntegratorOrder () const
 
void Interpolate (G4double tau, G4double yOut[])
 
void InterpolateHigh (G4double tau, G4double yOut[]) const
 
G4bool IsFSAL () const
 
void NormalisePolarizationVector (G4double vec[12])
 
void NormaliseTangentVector (G4double vec[6])
 
G4BogackiShampine45 & operator= (const G4BogackiShampine45 &)=delete
 
void PrepareConstants ()
 
void ResetfNORHSCalls ()
 
void RightHandSide (const G4double y[], G4double dydx[]) const
 
void RightHandSide (const G4double y[], G4double dydx[], G4double field[]) const
 
void SetEquationOfMotion (G4EquationOfMotion *newEquation)
 
void SetupInterpolation ()
 
void SetupInterpolationHigh ()
 
void Stepper (const G4double y[], const G4double dydx[], G4double h, G4double yout[], G4double yerr[])
 
 ~G4BogackiShampine45 ()
 

Protected Member Functions

void SetFSAL (G4bool flag=true)
 
void SetIntegrationOrder (G4int order)
 

Private Attributes

G4double * ak10
 
G4double * ak11
 
G4double * ak2
 
G4double * ak3
 
G4double * ak4
 
G4double * ak5
 
G4double * ak6
 
G4double * ak7
 
G4double * ak8
 
G4double * ak9
 
G4BogackiShampine45 * fAuxStepper = nullptr
 
G4EquationOfMotion * fEquation_Rhs = nullptr
 
G4int fIntegrationOrder = -1
 
G4bool fIsFSAL = false
 
G4double * fLastDyDx
 
G4double * fLastFinalVector
 
G4double * fLastInitialVector
 
G4double fLastStepLength = -1.0
 
G4double * fMidError
 
G4double * fMidVector
 
const G4int fNoIntegrationVariables = 0
 
unsigned long fNoRHSCalls = 0UL
 
const G4int fNoStateVariables = 0
 
G4bool fPreparedInterpolation = false
 
G4double * p [6]
 
G4double * yIn
 
G4double * yTemp
 

Static Private Attributes

static G4double bi [12][7]
 
static G4bool fPreparedConstants = false
 

Detailed Description

Definition at line 51 of file G4BogackiShampine45.hh.

Constructor & Destructor Documentation

◆ G4BogackiShampine45() [1/2]

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

Definition at line 70 of file G4BogackiShampine45.cc.

73 : G4MagIntegratorStepper(EqRhs, noIntegrationVariables)
74{
75 const G4int numberOfVariables = noIntegrationVariables;
76
77 // New Chunk of memory being created for use by the stepper
78
79 // aki - for storing intermediate RHS
80 ak2 = new G4double[numberOfVariables];
81 ak3 = new G4double[numberOfVariables];
82 ak4 = new G4double[numberOfVariables];
83 ak5 = new G4double[numberOfVariables];
84 ak6 = new G4double[numberOfVariables];
85 ak7 = new G4double[numberOfVariables];
86 ak8 = new G4double[numberOfVariables];
87 ak9 = new G4double[numberOfVariables];
88 ak10 = new G4double[numberOfVariables];
89 ak11 = new G4double[numberOfVariables];
90
91 for (auto i = 0; i < 6; ++i)
92 {
93 p[i]= new G4double[numberOfVariables];
94 }
95
96 assert ( GetNumberOfStateVariables() >= 8 );
97 const G4int numStateVars = std::max(noIntegrationVariables,
98 GetNumberOfStateVariables() );
99
100 // Must ensure space extra 'state' variables exists - i.e. yIn[7]
101 yTemp = new G4double[numStateVars];
102 yIn = new G4double[numStateVars] ;
103
104 fLastInitialVector = new G4double[numStateVars] ;
105 fLastFinalVector = new G4double[numStateVars] ;
106 fLastDyDx = new G4double[numberOfVariables]; // Only derivatives
107
108 fMidVector = new G4double[numberOfVariables];
109 fMidError = new G4double[numberOfVariables];
110
111 if( ! fPreparedConstants )
112 {
113 PrepareConstants();
114 }
115
116 if( primary )
117 {
118 fAuxStepper = new G4BogackiShampine45(EqRhs, numberOfVariables, false);
119 }
120}
G4double
double G4double
Definition: G4Types.hh:83
G4int
int G4int
Definition: G4Types.hh:85
G4BogackiShampine45::G4BogackiShampine45
G4BogackiShampine45(G4EquationOfMotion *EqRhs, G4int numberOfVariables=6, G4bool primary=true)
Definition: G4BogackiShampine45.cc:70
G4BogackiShampine45::fAuxStepper
G4BogackiShampine45 * fAuxStepper
Definition: G4BogackiShampine45.hh:115
G4BogackiShampine45::fPreparedConstants
static G4bool fPreparedConstants
Definition: G4BogackiShampine45.hh:120
G4MagIntegratorStepper::G4MagIntegratorStepper
G4MagIntegratorStepper(G4EquationOfMotion *Equation, G4int numIntegrationVariables, G4int numStateVariables=12, G4bool isFSAL=false)
Definition: G4MagIntegratorStepper.cc:38
G4MagIntegratorStepper::GetNumberOfStateVariables
G4int GetNumberOfStateVariables() const
G4INCL::Math::max
T max(const T t1, const T t2)
brief Return the largest of the two arguments
Definition: G4INCLGlobals.hh:112

References ak10, ak11, ak2, ak3, ak4, ak5, ak6, ak7, ak8, ak9, fAuxStepper, fLastDyDx, fLastFinalVector, fLastInitialVector, fMidError, fMidVector, fPreparedConstants, G4BogackiShampine45(), G4MagIntegratorStepper::GetNumberOfStateVariables(), G4INCL::Math::max(), p, PrepareConstants(), yIn, and yTemp.

Referenced by G4BogackiShampine45().

◆ ~G4BogackiShampine45()

G4BogackiShampine45::~G4BogackiShampine45 ( )

Definition at line 124 of file G4BogackiShampine45.cc.

125{
126 // Clear all previously allocated memory for stepper and DistChord
127 //
128 delete [] ak2;
129 delete [] ak3;
130 delete [] ak4;
131 delete [] ak5;
132 delete [] ak6;
133 delete [] ak7;
134 delete [] ak8;
135 delete [] ak9;
136 delete [] ak10;
137 delete [] ak11;
138
139 for (auto i = 0; i < 6; ++i)
140 {
141 delete [] p[i];
142 }
143
144 delete [] yTemp;
145 delete [] yIn;
146
147 delete [] fLastInitialVector;
148 delete [] fLastFinalVector;
149 delete [] fLastDyDx;
150 delete [] fMidVector;
151 delete [] fMidError;
152
153 delete fAuxStepper;
154}

References ak10, ak11, ak2, ak3, ak4, ak5, ak6, ak7, ak8, ak9, fAuxStepper, fLastDyDx, fLastFinalVector, fLastInitialVector, fMidError, fMidVector, p, yIn, and yTemp.

◆ G4BogackiShampine45() [2/2]

G4BogackiShampine45::G4BogackiShampine45 ( const G4BogackiShampine45 &  )
delete

Member Function Documentation

◆ DistChord()

G4double G4BogackiShampine45::DistChord ( ) const
virtual

Implements G4MagIntegratorStepper.

Definition at line 310 of file G4BogackiShampine45.cc.

311{
312 G4double distLine, distChord;
313 G4ThreeVector initialPoint, finalPoint, midPoint;
314
315 // Store last initial and final points
316 // (they will be overwritten in self-Stepper call!)
317 //
318 initialPoint = G4ThreeVector(fLastInitialVector[0],
319 fLastInitialVector[1], fLastInitialVector[2]);
320 finalPoint = G4ThreeVector(fLastFinalVector[0],
321 fLastFinalVector[1], fLastFinalVector[2]);
322
323#if 1
324 // Old method -- Do half a step using StepNoErr
325 //
326 fAuxStepper->Stepper( fLastInitialVector, fLastDyDx, 0.5*fLastStepLength,
327 fMidVector, fMidError);
328#else
329 // New method -- Using interpolation,
330 // requires only 3 extra stages (ie 3 extra field evaluations )
331
332 // Use Interpolation, instead of auxiliary stepper to evaluate midpoint
333 //
334 if( ! fPreparedInterpolation )
335 {
336 G4BogackiShampine45* cThis = const_cast<G4BogackiShampine45 *>(this);
337 cThis-> SetupInterpolationHigh();
338 }
339 // For calculating the output at the tau fraction of Step
340 //
341 G4double tau = 0.5;
342 InterpolateHigh( tau, fMidVector );
343#endif
344
345 midPoint = G4ThreeVector( fMidVector[0], fMidVector[1], fMidVector[2]);
346
347 // Use stored values of Initial and Endpoint + new Midpoint to evaluate
348 // distance of Chord
349
350 if (initialPoint != finalPoint)
351 {
352 distLine = G4LineSection::Distline( midPoint,initialPoint,finalPoint );
353 distChord = distLine;
354 }
355 else
356 {
357 distChord = (midPoint-initialPoint).mag();
358 }
359 return distChord;
360}
G4ThreeVector
CLHEP::Hep3Vector G4ThreeVector
Definition: G4ThreeVector.hh:36
G4BogackiShampine45::Stepper
void Stepper(const G4double y[], const G4double dydx[], G4double h, G4double yout[], G4double yerr[])
Definition: G4BogackiShampine45.cc:171
G4BogackiShampine45::InterpolateHigh
void InterpolateHigh(G4double tau, G4double yOut[]) const
Definition: G4BogackiShampine45.cc:563
G4LineSection::Distline
static G4double Distline(const G4ThreeVector &OtherPnt, const G4ThreeVector &LinePntA, const G4ThreeVector &LinePntB)
Definition: G4LineSection.hh:72

References G4LineSection::Distline(), fAuxStepper, fLastDyDx, fLastFinalVector, fLastInitialVector, fLastStepLength, fMidError, fMidVector, fPreparedInterpolation, InterpolateHigh(), SetupInterpolationHigh(), and Stepper().

◆ GetEquationOfMotion() [1/2]

G4EquationOfMotion * G4MagIntegratorStepper::GetEquationOfMotion ( )
inlineinherited

Referenced by G4MagInt_Driver::GetEquationOfMotion(), G4OldMagIntDriver::GetEquationOfMotion(), G4NystromRK4::GetField(), G4DormandPrince745::GetSpecificEquation(), G4NystromRK4::SetDistanceForConstantField(), G4RKG3_Stepper::StepNoErr(), and G4RKG3_Stepper::Stepper().

◆ GetEquationOfMotion() [2/2]

const G4EquationOfMotion * G4MagIntegratorStepper::GetEquationOfMotion ( ) const
inlineinherited

◆ GetfNoRHSCalls()

unsigned long G4MagIntegratorStepper::GetfNoRHSCalls ( )
inlineinherited

◆ GetLastDydx()

void G4BogackiShampine45::GetLastDydx ( G4double  dyDxLast[])

Definition at line 156 of file G4BogackiShampine45.cc.

157{
158 const G4int numberOfVariables = GetNumberOfVariables();
159
160 for(G4int i=0; i < numberOfVariables; ++i )
161 {
162 dyDxLast[i] = ak9[i];
163 }
164}
G4MagIntegratorStepper::GetNumberOfVariables
G4int GetNumberOfVariables() const

References ak9, and G4MagIntegratorStepper::GetNumberOfVariables().

◆ GetNumberOfStateVariables()

G4int G4MagIntegratorStepper::GetNumberOfStateVariables ( ) const
inlineinherited

Referenced by G4BogackiShampine45(), G4CashKarpRKF45::G4CashKarpRKF45(), G4SimpleRunge::G4SimpleRunge(), G4TSimpleRunge< T_Equation, N >::G4TSimpleRunge(), G4TsitourasRK45::G4TsitourasRK45(), G4BogackiShampine23::makeStep(), G4RK547FEq1::makeStep(), G4RK547FEq2::makeStep(), G4RK547FEq3::makeStep(), G4ConstRK4::Stepper(), G4MagErrorStepper::Stepper(), and G4TMagErrorStepper< T_Stepper, T_Equation, N >::Stepper().

◆ GetNumberOfVariables()

G4int G4MagIntegratorStepper::GetNumberOfVariables ( ) const
inlineinherited

Referenced by G4ExplicitEuler::DumbStepper(), G4ImplicitEuler::DumbStepper(), G4ClassicalRK4::DumbStepper(), G4ChordFinder::G4ChordFinder(), GetLastDydx(), G4DoLoMcPriRK34::Interpolate(), G4TsitourasRK45::Interpolate(), G4DormandPrince745::Interpolate4thOrder(), G4DormandPrince745::Interpolate5thOrder(), G4DormandPrinceRK56::Interpolate_high(), G4DormandPrinceRK56::Interpolate_low(), InterpolateHigh(), G4BogackiShampine23::makeStep(), G4RK547FEq1::makeStep(), G4RK547FEq2::makeStep(), G4RK547FEq3::makeStep(), G4DormandPrinceRK56::SetupInterpolate_high(), G4DormandPrinceRK56::SetupInterpolate_low(), G4DormandPrince745::SetupInterpolation5thOrder(), SetupInterpolationHigh(), Stepper(), G4CashKarpRKF45::Stepper(), G4DoLoMcPriRK34::Stepper(), G4DormandPrinceRK56::Stepper(), G4DormandPrinceRK78::Stepper(), G4MagErrorStepper::Stepper(), G4TsitourasRK45::Stepper(), and G4DormandPrince745::Stepper().

◆ IntegrationOrder()

G4int G4MagIntegratorStepper::IntegrationOrder ( )
inlineinherited

◆ IntegratorOrder()

G4int G4BogackiShampine45::IntegratorOrder ( ) const
inlinevirtual

Implements G4MagIntegratorStepper.

Definition at line 97 of file G4BogackiShampine45.hh.

97{ return 4; }

◆ Interpolate()

void G4BogackiShampine45::Interpolate ( G4double  tau,
G4double  yOut[] 
)
inline

Definition at line 82 of file G4BogackiShampine45.hh.

84 {
85 InterpolateHigh( tau, yOut);
86 // InterpolateHigh( yInput, dydx, Step, yOut, tau);
87 }

References InterpolateHigh().

◆ InterpolateHigh()

void G4BogackiShampine45::InterpolateHigh ( G4double  tau,
G4double  yOut[] 
) const

Definition at line 563 of file G4BogackiShampine45.cc.

564{
565 G4int numberOfVariables = GetNumberOfVariables();
566
567 G4Exception("G4BogackiShampine45::InterpolateHigh()", "GeomField0001",
568 FatalException, "Method is not yet validated.");
569
570 // const G4double *yIn= fLastInitialVector;
571 // const G4double *dydx= fLastDyDx;
572 const G4double Step = fLastStepLength;
573
574#if 1
575 G4int nwant = numberOfVariables;
576 const G4int norder= 6;
577 G4int l, k;
578
579 for (l = 0; l < nwant; ++l)
580 {
581 yOut[l] = p[norder-1][l] * tau;
582 }
583 for (k = norder - 2; k >= 1; --k)
584 {
585 for (l = 0; l < nwant; ++l)
586 {
587 yOut[l] = ( yOut[l] + p[k][l] ) * tau;
588 }
589 }
590 for (l = 0; l < nwant; ++l)
591 {
592 yOut[l] = ( yOut[l] + Step * ak8[l] ) * tau + yIn[l];
593 }
594 // The derivative at the end-point is nextDydx[i] = ak8[i];
595#else
596 // The scheme tries to do the same as the DormandPrince745 routine,
597 // but fails
598
599 G4double b[12];
600 const G4double* dydx = fLastDyDx;
601
602 G4double tau0 = tau;
603
604 for(G4int iStage=1; iStage<=11; ++iStage) // iStage = stage number
605 {
606 b[iStage] = 0.0;
607 tau = tau0;
608 for(G4int j=6; j>=1; --j) // j reversed
609 {
610 b[iStage] += bi[iStage][j] * tau;
611 tau *= tau0;
612 }
613 }
614
615 for(G4int i=0; i<numberOfVariables; ++i)
616 {
617 yOut[i] = yIn[i] + Step*(b[1]*dydx[i] + b[2]*ak2[i] + b[3]*ak3[i] +
618 b[4]*ak4[i] + b[5]*ak5[i] + b[6]*ak6[i] +
619 b[7]*ak7[i] + b[8]*ak8[i] + b[9]*ak9[i] +
620 b[10]*ak10[i] + b[11]*ak11[i] );
621 }
622#endif
623}
FatalException
@ FatalException
Definition: G4ExceptionSeverity.hh:61
G4Exception
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Definition: G4Exception.cc:35
G4BogackiShampine45::bi
static G4double bi[12][7]
Definition: G4BogackiShampine45.hh:121

References ak10, ak11, ak2, ak3, ak4, ak5, ak6, ak7, ak8, ak9, bi, FatalException, fLastDyDx, fLastStepLength, G4Exception(), G4MagIntegratorStepper::GetNumberOfVariables(), p, and yIn.

Referenced by DistChord(), and Interpolate().

◆ IsFSAL()

G4bool G4MagIntegratorStepper::IsFSAL ( ) const
inlineinherited

◆ NormalisePolarizationVector()

void G4MagIntegratorStepper::NormalisePolarizationVector ( G4double  vec[12])
inlineinherited

Referenced by G4SimpleHeum::DumbStepper(), and G4ClassicalRK4::DumbStepper().

◆ NormaliseTangentVector()

void G4MagIntegratorStepper::NormaliseTangentVector ( G4double  vec[6])
inlineinherited

◆ operator=()

G4BogackiShampine45 & G4BogackiShampine45::operator= ( const G4BogackiShampine45 &  )
delete

◆ PrepareConstants()

void G4BogackiShampine45::PrepareConstants ( )

Definition at line 487 of file G4BogackiShampine45.cc.

488{
489 for(auto i=1; i<= 11; ++i)
490 {
491 bi[i][1] = 0.0 ;
492 }
493
494 for(auto i=1; i<=6; ++i)
495 {
496 bi[2][i] = 0.0 ;
497 }
498
499 bi[1][6] = -12134338393.0 / 1050809760.0 ,
500 bi[1][5] = -1620741229.0 / 50038560.0 ,
501 bi[1][4] = -2048058893.0 / 59875200.0 ,
502 bi[1][3] = -87098480009.0 / 5254048800.0 ,
503 bi[1][2] = -11513270273.0 / 3502699200.0 ,
504 //
505 bi[3][6] = -33197340367.0 / 1218433216.0 ,
506 bi[3][5] = -539868024987.0 / 6092166080.0 ,
507 bi[3][4] = -39991188681.0 / 374902528.0 ,
508 bi[3][3] = -69509738227.0 / 1218433216.0 ,
509 bi[3][2] = -29327744613.0 / 2436866432.0 ,
510 //
511 bi[4][6] = -284800997201.0 / 19905339168.0 ,
512 bi[4][5] = -7896875450471.0 / 165877826400.0 ,
513 bi[4][4] = -333945812879.0 / 5671036800.0 ,
514 bi[4][3] = -16209923456237.0 / 497633479200.0 ,
515 bi[4][2] = -2382590741699.0 / 331755652800.0 ,
516 //
517 bi[5][6] = -540919.0 / 741312.0 ,
518 bi[5][5] = -103626067.0 / 43243200.0 ,
519 bi[5][4] = -633779.0 / 211200.0 ,
520 bi[5][3] = -32406787.0 / 18532800.0 ,
521 bi[5][2] = -36591193.0 / 86486400.0 ,
522 //
523 bi[6][6] = 7157998304.0 / 374350977.0 ,
524 bi[6][5] = 30405842464.0 / 623918295.0 ,
525 bi[6][4] = 183022264.0 / 5332635.0 ,
526 bi[6][3] = -3357024032.0 / 1871754885.0 ,
527 bi[6][2] = -611586736.0 / 89131185.0 ,
528 //
529 bi[7][6] = -138073.0 / 9408.0 ,
530 bi[7][5] = -719433.0 / 15680.0 ,
531 bi[7][4] = -1620541.0 / 31360.0 ,
532 bi[7][3] = -385151.0 / 15680.0 ,
533 bi[7][2] = -65403.0 / 15680.0 ,
534 //
535 bi[8][6] = 1245.0 / 64.0 ,
536 bi[8][5] = 3991.0 / 64.0 ,
537 bi[8][4] = 4715.0 / 64.0 ,
538 bi[8][3] = 2501.0 / 64.0 ,
539 bi[8][2] = 149.0 / 16.0 ,
540 bi[8][1] = 1.0 ,
541 //
542 bi[9][6] = 55.0 / 3.0 ,
543 bi[9][5] = 71.0 ,
544 bi[9][4] = 103.0 ,
545 bi[9][3] = 199.0 / 3.0 ,
546 bi[9][2] = 16.0 ,
547 //
548 bi[10][6] = -1774004627.0 / 75810735.0 ,
549 bi[10][5] = -1774004627.0 / 25270245.0 ,
550 bi[10][4] = -26477681.0 / 359975.0 ,
551 bi[10][3] = -11411880511.0 / 379053675.0 ,
552 bi[10][2] = -423642896.0 / 126351225.0 ,
553 //
554 bi[11][6] = 35.0 ,
555 bi[11][5] = 105.0 ,
556 bi[11][4] = 117.0 ,
557 bi[11][3] = 59.0 ,
558 bi[11][2] = 12.0 ;
559
560 fPreparedConstants = true;
561}

References bi, and fPreparedConstants.

Referenced by G4BogackiShampine45().

◆ ResetfNORHSCalls()

void G4MagIntegratorStepper::ResetfNORHSCalls ( )
inlineinherited

◆ RightHandSide() [1/2]

void G4MagIntegratorStepper::RightHandSide ( const G4double  y[],
G4double  dydx[] 
) const
inlineinherited

Referenced by G4MagInt_Driver::AccurateAdvance(), G4OldMagIntDriver::AccurateAdvance(), G4ImplicitEuler::DumbStepper(), G4SimpleHeum::DumbStepper(), G4SimpleRunge::DumbStepper(), G4ClassicalRK4::DumbStepper(), G4MagInt_Driver::GetDerivatives(), G4OldMagIntDriver::GetDerivatives(), G4BogackiShampine23::makeStep(), G4RK547FEq1::makeStep(), G4RK547FEq2::makeStep(), G4RK547FEq3::makeStep(), G4DormandPrinceRK56::SetupInterpolate_high(), G4DormandPrinceRK56::SetupInterpolate_low(), G4DormandPrince745::SetupInterpolation5thOrder(), SetupInterpolationHigh(), Stepper(), G4CashKarpRKF45::Stepper(), G4DoLoMcPriRK34::Stepper(), G4DormandPrinceRK56::Stepper(), G4DormandPrinceRK78::Stepper(), G4MagErrorStepper::Stepper(), G4TsitourasRK45::Stepper(), and G4DormandPrince745::Stepper().

◆ RightHandSide() [2/2]

void G4MagIntegratorStepper::RightHandSide ( const G4double  y[],
G4double  dydx[],
G4double  field[] 
) const
inlineinherited

◆ SetEquationOfMotion()

void G4MagIntegratorStepper::SetEquationOfMotion ( G4EquationOfMotion *  newEquation)
inlineinherited

Referenced by G4MagInt_Driver::SetEquationOfMotion(), and G4OldMagIntDriver::SetEquationOfMotion().

◆ SetFSAL()

void G4MagIntegratorStepper::SetFSAL ( G4bool  flag = true)
inlineprotectedinherited

Referenced by G4BogackiShampine23::G4BogackiShampine23().

◆ SetIntegrationOrder()

void G4MagIntegratorStepper::SetIntegrationOrder ( G4int  order)
inlineprotectedinherited

Referenced by G4BogackiShampine23::G4BogackiShampine23().

◆ SetupInterpolation()

void G4BogackiShampine45::SetupInterpolation ( )
inline

Definition at line 75 of file G4BogackiShampine45.hh.

76 {
77 SetupInterpolationHigh(); // ( yInput, dydx, Step);
78 }

References SetupInterpolationHigh().

◆ SetupInterpolationHigh()

void G4BogackiShampine45::SetupInterpolationHigh ( )

Definition at line 362 of file G4BogackiShampine45.cc.

363{
364 // Coefficients for the additional stages
365 //
366 const G4double
367 a91 = 455.0/6144.0 ,
368 a92 = 0.0 ,
369 a93 = 10256301.0/35409920.0 ,
370 a94 = 2307361.0/17971200.0 ,
371 a95 = -387.0/102400.0 ,
372 a96 = 73.0/5130.0 ,
373 a97 = -7267.0/215040.0 ,
374 a98 = 1.0/32.0 ,
375
376 a101 = -837888343715.0/13176988637184.0 ,
377 a102 = 30409415.0/52955362.0 ,
378 a103 = -48321525963.0/759168069632.0 ,
379 a104 = 8530738453321.0/197654829557760.0 ,
380 a105 = 1361640523001.0/1626788720640.0 ,
381 a106 = -13143060689.0/38604458898.0 ,
382 a107 = 18700221969.0/379584034816.0 ,
383 a108 = -5831595.0/847285792.0 ,
384 a109 = -5183640.0/26477681.0 ,
385
386 a111 = 98719073263.0/1551965184000.0 ,
387 a112 = 1307.0/123552.0 ,
388 a113 = 4632066559387.0/70181753241600.0 ,
389 a114 = 7828594302389.0/382182512025600.0 ,
390 a115 = 40763687.0/11070259200.0 ,
391 a116 = 34872732407.0/224610586200.0 ,
392 a117 = -2561897.0/30105600.0 ,
393 a118 = 1.0/10.0 ,
394 a119 = -1.0/10.0 ,
395 a1110 = -1403317093.0/11371610250.0 ;
396
397 const G4int numberOfVariables= this->GetNumberOfVariables();
398 const G4double* dydx= fLastDyDx;
399 const G4double Step = fLastStepLength;
400
401 yTemp[7] = yIn[7];
402
403 // Evaluate the extra stages
404 //
405 for(G4int i=0; i<numberOfVariables; ++i)
406 {
407 yTemp[i] = yIn[i] + Step*(a91*dydx[i] + a92*ak2[i] + a93*ak3[i] +
408 a94*ak4[i] + a95*ak5[i] + a96*ak6[i] +
409 a97*ak7[i] + a98*ak8[i] );
410 }
411
412 RightHandSide(yTemp, ak9); // 9th stage
413
414 for(G4int i=0; i<numberOfVariables; ++i)
415 {
416 yTemp[i] = yIn[i] + Step*(a101*dydx[i] + a102*ak2[i] + a103*ak3[i] +
417 a104*ak4[i] + a105*ak5[i] + a106*ak6[i] +
418 a107*ak7[i] + a108*ak8[i] + a109*ak9[i] );
419 }
420
421 RightHandSide(yTemp, ak10); // 10th stage
422
423 for(G4int i=0; i<numberOfVariables; ++i)
424 {
425 yTemp[i] = yIn[i] + Step*(a111*dydx[i] + a112*ak2[i] + a113*ak3[i] +
426 a114*ak4[i] + a115*ak5[i] + a116*ak6[i] +
427 a117*ak7[i] + a118*ak8[i] + a119*ak9[i] +
428 a1110*ak10[i] );
429 }
430 RightHandSide(yTemp, ak11); // 11th stage
431
432 // In future we can restrict the number of variables interpolated
433 //
434 G4int nwant = numberOfVariables;
435
436 // Form the coefficients of the interpolating polynomial in its shifted
437 // and scaled form. The terms are grouped to minimize the errors
438 // of the transformation, to cope with ill-conditioning. ( From RKSUITE )
439 //
440 for (G4int l = 0; l < nwant; ++l)
441 {
442 // Coefficient of tau^6
443 p[5][l] = bi[5][6]*ak5[l] +
444 ((bi[10][6]*ak10[l] + bi[8][6]*ak8[l]) +
445 (bi[7][6]*ak7[l] + bi[6][6]*ak6[l])) +
446 ((bi[4][6]*ak4[l] + bi[9][6]*ak9[l]) +
447 (bi[3][6]*ak3[l] + bi[11][6]*ak11[l]) +
448 bi[1][6]*dydx[l]);
449 // Coefficient of tau^5
450 p[4][l] = (bi[10][5]*ak10[l] + bi[9][5]*ak9[l]) +
451 ((bi[7][5]*ak7[l] + bi[6][5]*ak6[l]) +
452 bi[5][5]*ak5[l]) + ((bi[4][5]*ak4[l] +
453 bi[8][5]*ak8[l]) + (bi[3][5]*ak3[l] +
454 bi[11][5]*ak11[l]) + bi[1][5]*dydx[l]);
455 // Coefficient of tau^4
456 p[3][l] = ((bi[4][4]*ak4[l] + bi[8][4]*ak8[l]) +
457 (bi[7][4]*ak7[l] + bi[6][4]*ak6[l]) +
458 bi[5][4]*ak5[l]) + ((bi[10][4]*ak10[l] +
459 bi[9][4]*ak9[l]) + (bi[3][4]*ak3[l] +
460 bi[11][4]*ak11[l]) + bi[1][4]*dydx[l]);
461 // Coefficient of tau^3
462 p[2][l] = bi[5][3]*ak5[l] + bi[6][3]*ak6[l] +
463 ((bi[3][3]*ak3[l] + bi[9][3]*ak9[l]) +
464 (bi[10][3]*ak10[l]+ bi[8][3]*ak8[l]) + bi[1][3]*dydx[l]) +
465 ((bi[4][3]*ak4[l] + bi[11][3]*ak11[l]) + bi[7][3]*ak7[l]);
466 // Coefficient of tau^2
467 p[1][l] = bi[5][2]*ak5[l] + ((bi[6][2]*ak6[l] +
468 bi[8][2]*ak8[l]) + bi[1][2]*dydx[l]) +
469 ((bi[3][2]*ak3[l] + bi[9][2]*ak9[l]) +
470 bi[10][2]*ak10[l])+ ((bi[4][2]*ak4[l] +
471 bi[11][2]*ak2[l]) + bi[7][2]*ak7[l]);
472 }
473
474 // Scale all the coefficients by the step size.
475 //
476 for (G4int i = 0; i < 6; ++i)
477 {
478 for (G4int l = 0; l < nwant; ++l)
479 {
480 p[i][l] *= Step;
481 }
482 }
483
484 fPreparedInterpolation = true;
485}
G4MagIntegratorStepper::RightHandSide
void RightHandSide(const G4double y[], G4double dydx[]) const

References ak10, ak11, ak2, ak3, ak4, ak5, ak6, ak7, ak8, ak9, bi, fLastDyDx, fLastStepLength, fPreparedInterpolation, G4MagIntegratorStepper::GetNumberOfVariables(), p, G4MagIntegratorStepper::RightHandSide(), yIn, and yTemp.

Referenced by DistChord(), and SetupInterpolation().

◆ Stepper()

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

Implements G4MagIntegratorStepper.

Definition at line 171 of file G4BogackiShampine45.cc.

176{
177 G4int i;
178
179 // Constants from the Butcher tableu
180 //
181 const G4double
182 b21 = 1.0/6.0 ,
183 b31 = 2.0/27.0 , b32 = 4.0/27.0,
184
185 b41 = 183.0/1372.0 , b42 = -162.0/343.0, b43 = 1053.0/1372.0,
186
187 b51 = 68.0/297.0, b52 = -4.0/11.0,
188 b53 = 42.0/143.0, b54 = 1960.0/3861.0,
189
190 b61 = 597.0/22528.0, b62 = 81.0/352.0,
191 b63 = 63099.0/585728.0, b64 = 58653.0/366080.0,
192 b65 = 4617.0/20480.0,
193
194 b71 = 174197.0/959244.0, b72 = -30942.0/79937.0,
195 b73 = 8152137.0/19744439.0, b74 = 666106.0/1039181.0,
196 b75 = -29421.0/29068.0, b76 = 482048.0/414219.0,
197
198 b81 = 587.0/8064.0, b82 = 0.0,
199 b83 = 4440339.0/15491840.0, b84 = 24353.0/124800.0,
200 b85 = 387.0/44800.0, b86 = 2152.0/5985.0,
201 b87 = 7267.0/94080.0;
202
203// c1 = 2479.0/34992.0,
204// c2 = 0.0,
205// c3 = 123.0/416.0,
206// c4 = 612941.0/3411720.0,
207// c5 = 43.0/1440.0,
208// c6 = 2272.0/6561.0,
209// c7 = 79937.0/1113912.0,
210// c8 = 3293.0/556956.0,
211
212 // For the embedded higher order method only the difference of values
213 // taken and is used directly later (instead of defining the last row
214 // of Butcher table in separate constants and taking the
215 // difference)
216 //
217 const G4double
218 dc1 = b81 - 2479.0 / 34992.0 ,
219 dc2 = 0.0,
220 dc3 = b83 - 123.0 / 416.0 ,
221 dc4 = b84 - 612941.0 / 3411720.0,
222 dc5 = b85 - 43.0 / 1440.0,
223 dc6 = b86 - 2272.0 / 6561.0,
224 dc7 = b87 - 79937.0 / 1113912.0,
225 dc8 = - 3293.0 / 556956.0;
226
227 const G4int numberOfVariables = GetNumberOfVariables();
228
229 // The number of variables to be integrated over
230 //
231 yOut[7] = yTemp[7] = yIn[7] = yInput[7];
232
233 // Saving yInput because yInput and yOut can be aliases for same array
234 //
235 for(i=0; i<numberOfVariables; ++i)
236 {
237 yIn[i]=yInput[i];
238 }
239
240 // RightHandSide(yIn, dydx) ;
241 // 1st Step - Not doing, getting passed
242 //
243 for(i=0; i<numberOfVariables; ++i)
244 {
245 yTemp[i] = yIn[i] + b21*Step*DyDx[i] ;
246 }
247 RightHandSide(yTemp, ak2) ; // 2nd Step
248
249 for(i=0; i<numberOfVariables; ++i)
250 {
251 yTemp[i] = yIn[i] + Step*(b31*DyDx[i] + b32*ak2[i]) ;
252 }
253 RightHandSide(yTemp, ak3) ; // 3rd Step
254
255 for(i=0; i<numberOfVariables; ++i)
256 {
257 yTemp[i] = yIn[i] + Step*(b41*DyDx[i] + b42*ak2[i] + b43*ak3[i]) ;
258 }
259 RightHandSide(yTemp, ak4) ; // 4th Step
260
261 for(i=0; i<numberOfVariables; ++i)
262 {
263 yTemp[i] = yIn[i] + Step*(b51*DyDx[i] + b52*ak2[i] + b53*ak3[i] +
264 b54*ak4[i]) ;
265 }
266 RightHandSide(yTemp, ak5) ; // 5th Step
267
268 for(i=0; i<numberOfVariables; ++i)
269 {
270 yTemp[i] = yIn[i] + Step*(b61*DyDx[i] + b62*ak2[i] + b63*ak3[i] +
271 b64*ak4[i] + b65*ak5[i]) ;
272 }
273 RightHandSide(yTemp, ak6) ; // 6th Step
274
275 for(i=0; i<numberOfVariables; ++i)
276 {
277 yTemp[i] = yIn[i] + Step*(b71*DyDx[i] + b72*ak2[i] + b73*ak3[i] +
278 b74*ak4[i] + b75*ak5[i] + b76*ak6[i]);
279 }
280 RightHandSide(yTemp, ak7); // 7th Step
281
282 for(i=0; i<numberOfVariables; ++i)
283 {
284 yOut[i] = yIn[i] + Step*(b81*DyDx[i] + b82*ak2[i] + b83*ak3[i] +
285 b84*ak4[i] + b85*ak5[i] + b86*ak6[i] +
286 b87*ak7[i]);
287 }
288 RightHandSide(yOut, ak8); // 8th Step - Final one Using FSAL
289
290 for(i=0; i<numberOfVariables; ++i)
291 {
292 yErr[i] = Step*(dc1*DyDx[i] + dc2*ak2[i] + dc3*ak3[i] + dc4*ak4[i] +
293 dc5*ak5[i] + dc6*ak6[i] + dc7*ak7[i] + dc8*ak8[i]) ;
294
295 // Store Input and Final values, for possible use in calculating chord
296 //
297 fLastInitialVector[i] = yIn[i] ;
298 fLastFinalVector[i] = yOut[i];
299 fLastDyDx[i] = DyDx[i];
300 }
301
302 fLastStepLength = Step;
303 fPreparedInterpolation= false;
304
305 return ;
306}

References ak2, ak3, ak4, ak5, ak6, ak7, ak8, fLastDyDx, fLastFinalVector, fLastInitialVector, fLastStepLength, fPreparedInterpolation, G4MagIntegratorStepper::GetNumberOfVariables(), G4MagIntegratorStepper::RightHandSide(), yIn, and yTemp.

Referenced by DistChord().

Field Documentation

◆ ak10

G4double * G4BogackiShampine45::ak10
private

Definition at line 106 of file G4BogackiShampine45.hh.

Referenced by G4BogackiShampine45(), InterpolateHigh(), SetupInterpolationHigh(), and ~G4BogackiShampine45().

◆ ak11

G4double * G4BogackiShampine45::ak11
private

Definition at line 106 of file G4BogackiShampine45.hh.

Referenced by G4BogackiShampine45(), InterpolateHigh(), SetupInterpolationHigh(), and ~G4BogackiShampine45().

◆ ak2

G4double* G4BogackiShampine45::ak2
private

Definition at line 105 of file G4BogackiShampine45.hh.

Referenced by G4BogackiShampine45(), InterpolateHigh(), SetupInterpolationHigh(), Stepper(), and ~G4BogackiShampine45().

◆ ak3

G4double * G4BogackiShampine45::ak3
private

Definition at line 105 of file G4BogackiShampine45.hh.

Referenced by G4BogackiShampine45(), InterpolateHigh(), SetupInterpolationHigh(), Stepper(), and ~G4BogackiShampine45().

◆ ak4

G4double * G4BogackiShampine45::ak4
private

Definition at line 105 of file G4BogackiShampine45.hh.

Referenced by G4BogackiShampine45(), InterpolateHigh(), SetupInterpolationHigh(), Stepper(), and ~G4BogackiShampine45().

◆ ak5

G4double * G4BogackiShampine45::ak5
private

Definition at line 105 of file G4BogackiShampine45.hh.

Referenced by G4BogackiShampine45(), InterpolateHigh(), SetupInterpolationHigh(), Stepper(), and ~G4BogackiShampine45().

◆ ak6

G4double * G4BogackiShampine45::ak6
private

Definition at line 105 of file G4BogackiShampine45.hh.

Referenced by G4BogackiShampine45(), InterpolateHigh(), SetupInterpolationHigh(), Stepper(), and ~G4BogackiShampine45().

◆ ak7

G4double * G4BogackiShampine45::ak7
private

Definition at line 105 of file G4BogackiShampine45.hh.

Referenced by G4BogackiShampine45(), InterpolateHigh(), SetupInterpolationHigh(), Stepper(), and ~G4BogackiShampine45().

◆ ak8

G4double * G4BogackiShampine45::ak8
private

Definition at line 105 of file G4BogackiShampine45.hh.

Referenced by G4BogackiShampine45(), InterpolateHigh(), SetupInterpolationHigh(), Stepper(), and ~G4BogackiShampine45().

◆ ak9

G4double * G4BogackiShampine45::ak9
private

Definition at line 106 of file G4BogackiShampine45.hh.

Referenced by G4BogackiShampine45(), GetLastDydx(), InterpolateHigh(), SetupInterpolationHigh(), and ~G4BogackiShampine45().

◆ bi

G4double G4BogackiShampine45::bi
staticprivate

Definition at line 121 of file G4BogackiShampine45.hh.

Referenced by InterpolateHigh(), PrepareConstants(), and SetupInterpolationHigh().

◆ fAuxStepper

G4BogackiShampine45* G4BogackiShampine45::fAuxStepper = nullptr
private

Definition at line 115 of file G4BogackiShampine45.hh.

Referenced by DistChord(), G4BogackiShampine45(), and ~G4BogackiShampine45().

◆ fEquation_Rhs

G4EquationOfMotion* G4MagIntegratorStepper::fEquation_Rhs = nullptr
privateinherited

Definition at line 124 of file G4MagIntegratorStepper.hh.

◆ fIntegrationOrder

G4int G4MagIntegratorStepper::fIntegrationOrder = -1
privateinherited

Definition at line 134 of file G4MagIntegratorStepper.hh.

◆ fIsFSAL

G4bool G4MagIntegratorStepper::fIsFSAL = false
privateinherited

Definition at line 136 of file G4MagIntegratorStepper.hh.

◆ fLastDyDx

G4double * G4BogackiShampine45::fLastDyDx
private

Definition at line 111 of file G4BogackiShampine45.hh.

Referenced by DistChord(), G4BogackiShampine45(), InterpolateHigh(), SetupInterpolationHigh(), Stepper(), and ~G4BogackiShampine45().

◆ fLastFinalVector

G4double * G4BogackiShampine45::fLastFinalVector
private

Definition at line 111 of file G4BogackiShampine45.hh.

Referenced by DistChord(), G4BogackiShampine45(), Stepper(), and ~G4BogackiShampine45().

◆ fLastInitialVector

G4double* G4BogackiShampine45::fLastInitialVector
private

Definition at line 111 of file G4BogackiShampine45.hh.

Referenced by DistChord(), G4BogackiShampine45(), Stepper(), and ~G4BogackiShampine45().

◆ fLastStepLength

G4double G4BogackiShampine45::fLastStepLength = -1.0
private

Definition at line 110 of file G4BogackiShampine45.hh.

Referenced by DistChord(), InterpolateHigh(), SetupInterpolationHigh(), and Stepper().

◆ fMidError

G4double * G4BogackiShampine45::fMidError
private

Definition at line 112 of file G4BogackiShampine45.hh.

Referenced by DistChord(), G4BogackiShampine45(), and ~G4BogackiShampine45().

◆ fMidVector

G4double * G4BogackiShampine45::fMidVector
private

Definition at line 112 of file G4BogackiShampine45.hh.

Referenced by DistChord(), G4BogackiShampine45(), and ~G4BogackiShampine45().

◆ fNoIntegrationVariables

const G4int G4MagIntegratorStepper::fNoIntegrationVariables = 0
privateinherited

Definition at line 125 of file G4MagIntegratorStepper.hh.

◆ fNoRHSCalls

unsigned long G4MagIntegratorStepper::fNoRHSCalls = 0UL
mutableprivateinherited

Definition at line 128 of file G4MagIntegratorStepper.hh.

◆ fNoStateVariables

const G4int G4MagIntegratorStepper::fNoStateVariables = 0
privateinherited

Definition at line 126 of file G4MagIntegratorStepper.hh.

◆ fPreparedConstants

G4bool G4BogackiShampine45::fPreparedConstants = false
staticprivate

Definition at line 120 of file G4BogackiShampine45.hh.

Referenced by G4BogackiShampine45(), and PrepareConstants().

◆ fPreparedInterpolation

G4bool G4BogackiShampine45::fPreparedInterpolation = false
private

Definition at line 117 of file G4BogackiShampine45.hh.

Referenced by DistChord(), SetupInterpolationHigh(), and Stepper().

◆ p

G4double* G4BogackiShampine45::p[6]
private

Definition at line 108 of file G4BogackiShampine45.hh.

Referenced by G4BogackiShampine45(), InterpolateHigh(), SetupInterpolationHigh(), and ~G4BogackiShampine45().

◆ yIn

G4double * G4BogackiShampine45::yIn
private

Definition at line 106 of file G4BogackiShampine45.hh.

Referenced by G4BogackiShampine45(), InterpolateHigh(), SetupInterpolationHigh(), Stepper(), and ~G4BogackiShampine45().

◆ yTemp

G4double * G4BogackiShampine45::yTemp
private

Definition at line 106 of file G4BogackiShampine45.hh.

Referenced by G4BogackiShampine45(), SetupInterpolationHigh(), Stepper(), and ~G4BogackiShampine45().

The documentation for this class was generated from the following files:
Generated by doxygen 1.9.3