Geant4 11.2.2
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
G4BogackiShampine23.cc
Go to the documentation of this file.
1//
2// ********************************************************************
3// * License and Disclaimer *
4// * *
5// * The Geant4 software is copyright of the Copyright Holders of *
6// * the Geant4 Collaboration. It is provided under the terms and *
7// * conditions of the Geant4 Software License, included in the file *
8// * LICENSE and available at http://cern.ch/geant4/license . These *
9// * include a list of copyright holders. *
10// * *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work make any representation or warranty, express or implied, *
14// * regarding this software system or assume any liability for its *
15// * use. Please see the license in the file LICENSE and URL above *
16// * for the full disclaimer and the limitation of liability. *
17// * *
18// * This code implementation is the result of the scientific and *
19// * technical work of the GEANT4 collaboration. *
20// * By using, copying, modifying or distributing the software (or *
21// * any work based on the software) you agree to acknowledge its *
22// * use in resulting scientific publications, and indicate your *
23// * acceptance of all terms of the Geant4 Software license. *
24// ********************************************************************
25//
26// G4BogackiShampine23 implementation
27//
28// Bogacki-Shampine - 4 - 3(2) non-FSAL implementation
29//
30// Implementation of the method proposed in the publication
31// "A 3(2) pair of Runge - Kutta formulas"
32// by P. Bogacki and L. F. Shampine,
33// Appl. Math. Lett., vol. 2, no. 4, pp. 321-325, Jan. 1989.
34//
35// The Bogacki shampine method has the following Butcher's tableau
36//
37// 0 |
38// 1/2|1/2
39// 3/4|0 3/4
40// 1 |2/9 1/3 4/9
41// -------------------
42// |2/9 1/3 4/9 0
43// |7/24 1/4 1/3 1/8
44//
45// Created: Somnath Banerjee, Google Summer of Code 2015, 20 May 2015
46// Supervision: John Apostolakis, CERN
47// --------------------------------------------------------------------
48
50#include "G4LineSection.hh"
51#include "G4FieldUtils.hh"
52
53using namespace field_utils;
54
56 G4int integrationVariables)
57 : G4MagIntegratorStepper(EqRhs, integrationVariables)
58{
60 SetFSAL(true);
61}
62
63void G4BogackiShampine23::makeStep(const G4double yInput[],
64 const G4double dydx[],
65 const G4double hstep,
66 G4double yOutput[],
67 G4double* dydxOutput,
68 G4double* yError) const
69{
70
73 {
74 yOutput[i] = yTemp[i] = yInput[i];
75 }
76
79
80 const G4double b21 = 0.5 ,
81 b31 = 0., b32 = 3.0 / 4.0,
82 b41 = 2.0 / 9.0, b42 = 1.0 / 3.0, b43 = 4.0 / 9.0;
83
84 const G4double dc1 = b41 - 7.0 / 24.0, dc2 = b42 - 1.0 / 4.0,
85 dc3 = b43 - 1.0 / 3.0, dc4 = - 1.0 / 8.0;
86
87 // RightHandSide(yInput, dydx);
88 for(G4int i = 0; i < GetNumberOfVariables(); ++i)
89 {
90 yTemp[i] = yInput[i] + b21 * hstep * dydx[i];
91 }
92
93 RightHandSide(yTemp, ak2);
94 for(G4int i = 0; i < GetNumberOfVariables(); ++i)
95 {
96 yTemp[i] = yInput[i] + hstep * (b31 * dydx[i] + b32 * ak2[i]);
97 }
98
99 RightHandSide(yTemp, ak3);
100 for(G4int i = 0; i < GetNumberOfVariables(); ++i)
101 {
102 yOutput[i] = yInput[i] + hstep * (b41*dydx[i] + b42*ak2[i] + b43*ak3[i]);
103 }
104
105 if ((dydxOutput != nullptr) && (yError != nullptr))
106 {
107 RightHandSide(yOutput, dydxOutput);
108 for(G4int i = 0; i < GetNumberOfVariables(); ++i)
109 {
110 yError[i] = hstep * (dc1 * dydx[i] + dc2 * ak2[i] +
111 dc3 * ak3[i] + dc4 * dydxOutput[i]);
112 }
113 }
114}
115
117 const G4double dydx[],
118 G4double hstep,
119 G4double yOutput[],
120 G4double yError[])
121{
122 copy(fyIn, yInput);
123 copy(fdydx, dydx);
124 fhstep = hstep;
125
126 makeStep(fyIn, fdydx, fhstep, fyOut, fdydxOut, yError);
127
128 copy(yOutput, fyOut);
129}
130
132 const G4double dydx[],
133 G4double hstep,
134 G4double yOutput[],
135 G4double yError[],
136 G4double dydxOutput[])
137{
138 copy(fyIn, yInput);
139 copy(fdydx, dydx);
140 fhstep = hstep;
141
142 makeStep(fyIn, fdydx, fhstep, fyOut, fdydxOut, yError);
143
144 copy(yOutput, fyOut);
145 copy(dydxOutput, fdydxOut);
146}
147
149{
151 makeStep(fyIn, fdydx, fhstep / 2., yMid);
152
153 const G4ThreeVector begin = makeVector(fyIn, Value3D::Position);
154 const G4ThreeVector mid = makeVector(yMid, Value3D::Position);
155 const G4ThreeVector end = makeVector(fyOut, Value3D::Position);
156
157 return G4LineSection::Distline(mid, begin, end);
158}
double G4double
Definition G4Types.hh:83
int G4int
Definition G4Types.hh:85
G4double DistChord() const override
void Stepper(const G4double yInput[], const G4double dydx[], G4double hstep, G4double yOutput[], G4double yError[]) override
G4BogackiShampine23(G4EquationOfMotion *EqRhs, G4int numberOfVariables=6)
static G4double Distline(const G4ThreeVector &OtherPnt, const G4ThreeVector &LinePntA, const G4ThreeVector &LinePntB)
G4int GetNumberOfVariables() const
void RightHandSide(const G4double y[], G4double dydx[]) const
void SetIntegrationOrder(G4int order)
void SetFSAL(G4bool flag=true)
G4int GetNumberOfStateVariables() const
G4ThreeVector makeVector(const ArrayType &array, Value3D value)
void copy(G4double dst[], const G4double src[], std::size_t size=G4FieldTrack::ncompSVEC)