| Geant4 Cross Reference |
1 // ********************************************************************
2 // * License and Disclaimer *
3 // * *
4 // * The Geant4 software is copyright of the Copyright Holders of *
5 // * the Geant4 Collaboration. It is provided under the terms and *
6 // * conditions of the Geant4 Software License, included in the file *
7 // * LICENSE and available at http://cern.ch/geant4/license . These *
8 // * include a list of copyright holders. *
9 // * *
10 // * Neither the authors of this software system, nor their employing *
11 // * institutes,nor the agencies providing financial support for this *
12 // * work make any representation or warranty, express or implied, *
13 // * regarding this software system or assume any liability for its *
14 // * use. Please see the license in the file LICENSE and URL above *
15 // * for the full disclaimer and the limitation of liability. *
16 // * *
17 // * This code implementation is the result of the scientific and *
18 // * technical work of the GEANT4 collaboration. *
19 // * By using, copying, modifying or distributing the software (or *
20 // * any work based on the software) you agree to acknowledge its *
21 // * use in resulting scientific publications, and indicate your *
22 // * acceptance of all terms of the Geant4 Software license. *
23 // ********************************************************************
24 //
25 // G4BulirschStoer
26 //
27 // Class description:
28 //
29 // The Bulirsch-Stoer is a controlled driver that adjusts both step size
30 // and order of the method. The algorithm uses the modified midpoint and
31 // a polynomial extrapolation computes the solution.
32
33 // Author: Dmitry Sorokin, Google Summer of Code 2016
34 // Supervision: John Apostolakis, CERN
35 // --------------------------------------------------------------------
36 #ifndef G4BULIRSCH_STOER_HH
37 #define G4BULIRSCH_STOER_HH
38
39 #include "G4ModifiedMidpoint.hh"
40
41 #include "G4FieldTrack.hh"
42
43 class G4BulirschStoer
44 {
45 public:
46
47 enum class step_result
48 {
49 success,
50 fail
51 };
52
53 G4BulirschStoer( G4EquationOfMotion* equation, G4int nvar,
54 G4double eps_rel, G4double max_dt = DBL_MAX);
55
56 inline void set_max_dt(G4double max_dt);
57 inline void set_max_relative_error(G4double eps_rel);
58
59 // Stepper method
60 //
61 step_result try_step(const G4double in[], const G4double dxdt[],
62 G4double& t, G4double out[], G4double& dt);
63
64 // Reset the internal state of the stepper
65 //
66 void reset();
67
68 inline void SetEquationOfMotion(G4EquationOfMotion* equation);
69 inline G4EquationOfMotion* GetEquationOfMotion();
70
71 inline G4int GetNumberOfVariables() const;
72
73 private:
74
75 const static G4int m_k_max = 8;
76
77 void extrapolate(std::size_t k, G4double xest[]);
78 G4double calc_h_opt(G4double h, G4double error, std::size_t k) const;
79
80 G4bool set_k_opt(std::size_t k, G4double& dt);
81 G4bool in_convergence_window(G4int k) const;
82 G4bool should_reject(G4double error, G4int k) const;
83
84 // Number of vars to be integrated
85 G4int fnvar;
86
87 // Relative tolerance
88 G4double m_eps_rel;
89
90 // Modified midpoint algorithm
91 G4ModifiedMidpoint m_midpoint;
92
93 G4bool m_last_step_rejected{false};
94 G4bool m_first{true};
95
96 G4double m_dt_last{0.0};
97 // G4double m_t_last;
98
99 // Max allowed time step
100 G4double m_max_dt;
101
102 G4int m_current_k_opt;
103
104 // G4double m_xnew[G4FieldTrack::ncompSVEC];
105 G4double m_err[G4FieldTrack::ncompSVEC];
106 // G4double m_dxdt[G4FieldTrack::ncompSVEC];
107
108 // Stores the successive interval counts
109 G4int m_interval_sequence[m_k_max+1];
110
111 // Extrapolation coeffs (Neville’s algorithm)
112 G4double m_coeff[m_k_max+1][m_k_max];
113
114 // Costs for interval count
115 G4int m_cost[m_k_max+1];
116
117 // Sequence of states for extrapolation
118 G4double m_table[m_k_max][G4FieldTrack::ncompSVEC];
119
120 // Optimal step size
121 G4double h_opt[m_k_max+1];
122
123 // Work per unit step
124 G4double work[m_k_max+1];
125 };
126
127 #include "G4BulirschStoer.icc"
128
129 #endif
130