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>
5 namespace 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 
32  const RKIntegrator::RKData::Data & s,
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 (!(((float) s.time+h - (float) s.time) > 0) ) {
70  std::cerr << "Warning, RK Integrator step underflow" << std::endl;
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  }
80  else {
81  hnext = Rmax*h;
82  }
83  }
84  break;
85  }
86  stepsize=hnext;
87  return;
88  }
89 
90 
91 
93  delete eeStepper;
94  }
95 
97  return new AdaptiveRKStepper(*this);
98  }
99 
101  }
102 
104  return T;
105  }
106 
107  const double & AdaptiveRKStepper::tolerance() const{
108  return T;
109  }
110 
112  return sStepsize;
113  }
114  const double & AdaptiveRKStepper::startingStepsize() const{
115  return sStepsize;
116  }
117 
119  return S;
120  }
121 
122  const double & AdaptiveRKStepper::safetyFactor() const{
123  return S;
124  }
125 
127  return Rmin;
128  }
129  const double & AdaptiveRKStepper::rmin() const{
130  return Rmin;
131  }
132 
134  return Rmax;
135  }
136  const double & AdaptiveRKStepper::rmax() const{
137  return Rmax;
138  }
139 
140 }