adevs
adevs_event_locators.h
1 
31 #ifndef _adevs_event_locators_h_
32 #define _adevs_event_locators_h_
33 #include "adevs_hybrid.h"
34 #include <cmath>
35 #include <iostream>
36 
37 namespace adevs
38 {
39 
44 template <typename X> class event_locator_impl:
45  public event_locator<X>
46 {
47  public:
54  enum Mode { INTERPOLATE, BISECTION };
55  event_locator_impl(ode_system<X>* sys, double err_tol, Mode mode);
57  bool find_events(bool* events, const double* qstart, double* qend,
58  ode_solver<X>* solver, double& h);
59  private:
60  double *z[2]; // State events at the start and end of [0,h]
61  const double err_tol; // Error tolerance
62  Mode mode;
63 
64  int sign(double x) const
65  {
66  if (x < 0.0) return -1;
67  else if (x > 0.0) return 1;
68  else return 0;
69  }
70 };
71 
72 template <typename X>
73 event_locator_impl<X>::event_locator_impl(ode_system<X>* sys,
74  double err_tol, Mode mode):
75  event_locator<X>(sys),
76  err_tol(err_tol),
77  mode(mode)
78 {
79  z[0] = new double[sys->numEvents()];
80  z[1] = new double[sys->numEvents()];
81 }
82 
83 template <typename X>
84 event_locator_impl<X>::~event_locator_impl()
85 {
86  delete [] z[0]; delete [] z[1];
87 }
88 
89 template <typename X>
91  const double* qstart, double* qend, ode_solver<X>* solver, double& h)
92 {
93  // Calculate the state event functions at the start
94  // of the interval
95  this->sys->state_event_func(qstart,z[0]);
96  // Look for the first event inside of the interval [0,h]
97  for (;;)
98  {
99  double tguess = h;
100  bool event_in_interval = false, found_event = false;
101  this->sys->state_event_func(qend,z[1]);
102  // Do any of the z functions change sign? Have we found an event?
103  for (int i = 0; i < this->sys->numEvents(); i++)
104  {
105  events[i] = false;
106  if (sign(z[1][i]) != sign(z[0][i]))
107  {
108  // Event at h > 0
109  if (fabs(z[1][i]) <= err_tol)
110  {
111  events[i] = found_event = true;
112  }
113  // There is an event in (0,h)
114  else
115  {
116  if (mode == INTERPOLATE)
117  {
118  double tcandidate = z[0][i]*h/(z[0][i]-z[1][i]);
119  // Don't let the step size go to zero
120  if (tcandidate < h/4.0) tcandidate = h/4.0;
121  if (tcandidate < tguess) tguess = tcandidate;
122  }
123  event_in_interval = true;
124  }
125  }
126  }
127  // Guess at a new h and calculate qend for that time
128  if (event_in_interval)
129  {
130  if (mode == BISECTION) h /= 2.0;
131  else h = tguess;
132  for (int i = 0; i < this->sys->numVars(); i++)
133  qend[i] = qstart[i];
134  solver->advance(qend,h);
135  }
136  else return found_event;
137  }
138  // Will never reach this line
139  return false;
140 }
141 
145 template <typename X>
147  public event_locator_impl<X>
148 {
149  public:
150  bisection_event_locator(ode_system<X>* sys, double err_tol):
152  sys,err_tol,event_locator_impl<X>::BISECTION){}
153 };
154 
158 template <typename X>
160  public event_locator_impl<X>
161 {
162  public:
163  linear_event_locator(ode_system<X>* sys, double err_tol):
166 };
167 
168 } // end of namespace
169 
170 #endif
Mode
Definition: adevs_event_locators.h:54
Definition: adevs_event_locators.h:159
Definition: adevs_event_locators.h:44
Definition: adevs_hybrid.h:396
Definition: adevs_event_locators.h:146
virtual void advance(double *q, double h)=0
bool find_events(bool *events, const double *qstart, double *qend, ode_solver< X > *solver, double &h)
Definition: adevs_event_locators.h:90
Definition: adevs_hybrid.h:45
Definition: adevs_hybrid.h:367