CoolProp  6.6.0
An open-source fluid property and humid air property database
ODEIntegrators.cpp
Go to the documentation of this file.
1 
2 #include "ODEIntegrators.h"
3 #include "Eigen/Core"
4 #include "CPstrings.h"
5 #include "Exceptions.h"
6 #include <algorithm>
7 
8 bool ODEIntegrators::AdaptiveRK54(AbstractODEIntegrator& ode, double tstart, double tend, double hmin, double hmax, double eps_allowed,
9  double step_relax) {
10  // Get the starting array of variables of integration
11  std::vector<double> xold = ode.get_initial_array();
12  const long N = static_cast<long>(xold.size());
13 
14  // Start at an index of 0
15  int Itheta = 0;
16  double t0 = tstart;
17  double h = hmin;
18 
19  // Figure out if t is increasing or decreasing in the integration and set a flag
20  bool forwards_integration = ((tend - tstart) > 0);
21  // If backwards integration, flip the sign of the step
22  if (!forwards_integration) {
23  h *= -1;
24  }
25 
26  double max_error;
27 
28  std::vector<double> xnew1(N), xnew2(N), xnew3(N), xnew4(N), xnew5(N), f1(N), f2(N), f3(N), f4(N), f5(N), f6(N), error(N), xnew(N);
29 
30  // t is the independent variable here, where t takes on values in the bounded range [tmin,tmax]
31  do {
32 
33  // Check for termination
34  bool abort = ode.premature_termination();
35  if (abort) {
36  return abort;
37  }
38 
39  bool stepAccepted = false, disableAdaptive = false;
40 
41  while (!stepAccepted) {
42 
43  // reset the flag
44  disableAdaptive = false;
45 
46  // If the step would go beyond the end of the region of integration,
47  // just take a step to the end of the region of integration
48  if (forwards_integration && (t0 + h > tend)) {
49  disableAdaptive = true;
50  h = tend - t0;
51  }
52  if (!forwards_integration && (t0 + h < tend)) {
53  disableAdaptive = true;
54  h = tend - t0;
55  }
56 
57  ode.pre_step_callback();
58 
59  // We check stepAccepted again because if the derived class
60  // sets the variable stepAccepted, we should not actually do the evaluation
61  if (!stepAccepted) {
62 
63  Eigen::Map<Eigen::VectorXd> xold_w(&(xold[0]), N);
64 
65  if (std::abs(h) < hmin && !disableAdaptive) {
66  // Step is too small, just use the minimum step size
67  h = (forwards_integration) ? hmin : -hmin;
68  disableAdaptive = true;
69  }
70 
71  // Step 1: derivatives evaluated at old values
72  ode.derivs(t0, xold, f1);
73 
74  // Call post derivative callback after the first derivative evaluation (which might cache values)
75  ode.post_deriv_callback();
76 
77  Eigen::Map<Eigen::VectorXd> xnew1_w(&(xnew1[0]), N), f1_w(&(f1[0]), N);
78  xnew1_w = xold_w + h * (1.0 / 5.0) * f1_w;
79 
80  ode.derivs(t0 + 1.0 / 5.0 * h, xnew1, f2);
81  Eigen::Map<Eigen::VectorXd> xnew2_w(&(xnew2[0]), N), f2_w(&(f2[0]), N);
82  xnew2_w = xold_w + h * (+3.0 / 40.0 * f1_w + 9.0 / 40.0 * f2_w);
83 
84  ode.derivs(t0 + 3.0 / 10.0 * h, xnew2, f3);
85  Eigen::Map<Eigen::VectorXd> xnew3_w(&(xnew3[0]), N), f3_w(&(f3[0]), N);
86  xnew3_w = xold_w + h * (3.0 / 10.0 * f1_w - 9.0 / 10.0 * f2_w + 6.0 / 5.0 * f3_w);
87 
88  ode.derivs(t0 + 3.0 / 5.0 * h, xnew3, f4);
89  Eigen::Map<Eigen::VectorXd> xnew4_w(&(xnew4[0]), N), f4_w(&(f4[0]), N);
90  xnew4_w = xold_w + h * (-11.0 / 54.0 * f1_w + 5.0 / 2.0 * f2_w - 70.0 / 27.0 * f3_w + 35.0 / 27.0 * f4_w);
91 
92  ode.derivs(t0 + h, xnew4, f5);
93  Eigen::Map<Eigen::VectorXd> xnew5_w(&(xnew5[0]), N), f5_w(&(f5[0]), N);
94  xnew5_w =
95  xold_w
96  + h * (1631.0 / 55296 * f1_w + 175.0 / 512.0 * f2_w + 575.0 / 13824.0 * f3_w + 44275.0 / 110592.0 * f4_w + 253.0 / 4096.0 * f5_w);
97 
98  // Updated values at the next step using 5-th order
99  ode.derivs(t0 + 7.0 / 8.0 * h, xnew5, f6);
100  Eigen::Map<Eigen::VectorXd> xnew_w(&(xnew[0]), N), f6_w(&(f6[0]), N);
101  xnew_w = xold_w + h * (37.0 / 378.0 * f1_w + 250.0 / 621.0 * f3_w + 125.0 / 594.0 * f4_w + 512.0 / 1771.0 * f6_w);
102 
103  Eigen::Map<Eigen::VectorXd> error_w(&(error[0]), N);
104  error_w =
105  h
106  * (-277.0 / 64512.0 * f1_w + 6925.0 / 370944.0 * f3_w - 6925.0 / 202752.0 * f4_w - 277.0 / 14336.0 * f5_w + 277.0 / 7084.0 * f6_w);
107 
108  max_error = error_w.norm();
109 
110  // If the error is too large, make the step size smaller and try
111  // the step again
112  if (disableAdaptive) {
113  // Accept the step regardless of whether the error
114  // is too large or not
115  stepAccepted = true;
116  } else {
117  if (max_error > eps_allowed) {
118  // Take a smaller step next time, try again on this step
119  // But only if adaptive mode is on
120  // If eps_allowed == max_error (approximately), force the step to change to avoid infinite loop
121  h *= std::min(step_relax * pow(eps_allowed / max_error, 0.3), 0.999);
122  stepAccepted = false;
123  } else {
124  stepAccepted = true;
125  }
126  }
127 
128  } else {
129  std::cout << format("accepted");
130  }
131  }
132 
133  // Step has been accepted, update variables
134  t0 += h;
135  Itheta += 1;
136  xold = xnew;
137 
138  ode.post_step_callback(t0, h, xnew);
139 
140  // The error is already below the threshold
141  if (max_error < eps_allowed && disableAdaptive == false && max_error > 0) {
142  // Take a bigger step next time, since eps_allowed>max_error, but don't
143  // let the steps get much larger too quickly
144  h *= step_relax * pow(eps_allowed / max_error, 0.2);
145  }
146 
147  // Constrain the step to not be too large
148  if (forwards_integration) {
149  h = std::min(h, hmax);
150  } else {
151  h = -std::min(std::abs(h), hmax);
152  }
153 
154  // Overshot the end, oops... That's an error
155  if (forwards_integration && (t0 - tend > +1e-3)) {
156  throw CoolProp::ValueError(format("t0 - tend [%g] > 1e-3", t0 - tend));
157  }
158  if (!forwards_integration && (t0 - tend < -1e-3)) {
159  throw CoolProp::ValueError(format("t0 - tend [%g] < -1e-3", t0 - tend));
160  }
161  } while (((forwards_integration) && t0 < tend - 1e-10) || ((!forwards_integration) && t0 > tend + 1e-10));
162 
163  // No termination was requested
164  return false;
165 }
166 
167 #if defined(ENABLE_CATCH)
168 # include <catch2/catch_all.hpp>
169 
170 TEST_CASE("Integrate y'=y", "[ODEIntegrator]") {
171  class SimpleODEIntegrator : public ODEIntegrators::AbstractODEIntegrator
172  {
173  public:
174  std::vector<double> t, h, y;
175 
176  virtual std::vector<double> get_initial_array() const {
177  return std::vector<double>(1, 1);
178  }
179 
180  virtual void pre_step_callback(){};
181 
182  virtual void post_deriv_callback(){};
183 
184  virtual void post_step_callback(double t, double h, std::vector<double>& y) {
185  this->t.push_back(t);
186  this->h.push_back(h);
187  this->y.push_back(y[0]);
188  };
189 
190  virtual bool premature_termination() {
191  return false;
192  };
193 
194  virtual void derivs(double t, std::vector<double>& y, std::vector<double>& yprime) {
195  yprime[0] = y[0];
196  };
197  };
198 
199  SimpleODEIntegrator simple;
200  ODEIntegrators::AdaptiveRK54(simple, 0, 4, 1e-4, 0.5, 1e-7, 0.9);
201  double yfinal_integration = simple.y[simple.y.size() - 1];
202  double tfinal_integration = simple.t[simple.t.size() - 1];
203 
204  double yfinal_analytic = exp(4.0);
205  double error = yfinal_integration / yfinal_analytic - 1;
206 
207  CAPTURE(yfinal_analytic);
208  CAPTURE(yfinal_integration);
209  CAPTURE(tfinal_integration);
210  CHECK(std::abs(error) < 1e-6);
211  CHECK(std::abs(tfinal_integration - 4) < 1e-10);
212  int rr = 0;
213 }
214 #endif