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 (!(((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 
103  }
104 
106  return T;
107  }
108 
109  const double & AdaptiveRKStepper::tolerance() const{
110  return T;
111  }
112 
114  return sStepsize;
115  }
116  const double & AdaptiveRKStepper::startingStepsize() const{
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 unsigned int order() const =0
virtual void step(const RKIntegrator::RKData *data, const RKIntegrator::RKData::Data &sdata, RKIntegrator::RKData::Data &ddata, std::vector< double > &errors) const =0
AdaptiveRKStepper(const EEStepper *eeStepper=NULL)
virtual void step(const RKIntegrator::RKData *data, const RKIntegrator::RKData::Data &sdata, RKIntegrator::RKData::Data &ddata, double timeLimit) const
virtual AdaptiveRKStepper * clone() const