CLHEP VERSION Reference Documentation
   
CLHEP Home Page     CLHEP Documentation     CLHEP Bug Reports

AdaptiveRKStepper.cc
Go to the documentation of this file.
3#include <cmath>
4#include <stdexcept>
5namespace Genfun {
6
8 eeStepper(stepper ? stepper->clone():new EmbeddedRKStepper()),
9 T(1.0E-6),
10 sStepsize(0.01),
11 S(0.9),
12 Rmin(0.0),
13 Rmax(5.0),
14 stepsize(sStepsize)
15 {
16 }
17
19 RKStepper(right),
20 eeStepper(right.eeStepper->clone()),
21 T(right.T),
22 sStepsize(right.sStepsize),
23 S(right.S),
24 Rmin(right.Rmin),
25 Rmax(right.Rmax),
26 stepsize(right.sStepsize)
27 {
28 }
29
30
34 double timeLimit) const {
35 //
36 // Adaptive stepsize control
37 //
38 if (s.time==0.0) {
39 stepsize=sStepsize;
40 }
41 const unsigned int p = eeStepper->order(); // Order of the stepper
42 const double deltaMax = T*std::pow(S/Rmax, (int)(p+1)); // Maximum error 4 adjustment.
43 const double TINY = 1.0E-30; // Denominator regularization
44 double hnext;
45 //
46 // Time limited step ?
47 //
48 d.time= timeLimit==0? s.time+stepsize : timeLimit;
49
50 //--------------------------------------//
51 // Take one step, from s to d: //
52 //--------------------------------------//
53 double h = d.time-s.time;
54 while (1) {
55 std::vector<double> errors;
56 eeStepper->step(data, s, d, errors);
57 if (timeLimit!=0.0) return;
58
59 // Take absolute value:
60 for (size_t e=0;e<errors.size();e++) errors[e] = fabs(errors[e]);
61
62 // Select the largest:
63 double delta = (*std::max_element(errors.begin(),errors.end()));
64 if (delta > T) {
65 //
66 // Bail out and try a smaller step.
67 //
68 h = std::max(S*h*std::pow(T/(delta + TINY), 1.0/(p+1)),Rmin*h);
69 if (!(((double) (s.time+h) - (double) s.time) > 0) ) {
70 throw std::runtime_error("Warning, RK Integrator step underflow");
71 }
72 d.time = s.time+h;
73 hnext=h;
74 continue;
75 }
76 else {
77 if (delta < deltaMax) {
78 hnext = S*h*std::pow(T/(delta + TINY),1.0/(p+1));
79 // stepsize is supposed to increase;
80 if (hnext<h) hnext=h;
81 }
82 else {
83 hnext = Rmax*h;
84 }
85 }
86 break;
87 }
88 stepsize=hnext;
89 return;
90 }
91
92
93
95 delete eeStepper;
96 }
97
99 return new AdaptiveRKStepper(*this);
100 }
101
104
106 return T;
107 }
108
109 const double & AdaptiveRKStepper::tolerance() const{
110 return T;
111 }
112
114 return sStepsize;
115 }
117 return sStepsize;
118 }
119
121 return S;
122 }
123
124 const double & AdaptiveRKStepper::safetyFactor() const{
125 return S;
126 }
127
129 return Rmin;
130 }
131 const double & AdaptiveRKStepper::rmin() const{
132 return Rmin;
133 }
134
136 return Rmax;
137 }
138 const double & AdaptiveRKStepper::rmax() const{
139 return Rmax;
140 }
141
142}
virtual void step(const RKIntegrator::RKData *data, const RKIntegrator::RKData::Data &sdata, RKIntegrator::RKData::Data &ddata, std::vector< double > &errors) const =0
virtual unsigned int order() const =0
virtual void step(const RKIntegrator::RKData *data, const RKIntegrator::RKData::Data &sdata, RKIntegrator::RKData::Data &ddata, double timeLimit) const
AdaptiveRKStepper(const EEStepper *eeStepper=NULL)
virtual AdaptiveRKStepper * clone() const