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, DISCONTINUOUS };
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>
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>
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  while (this->sys->numEvents() > 0)
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 (
110  // End condition when z is continuous
111  ((mode != DISCONTINUOUS) && (fabs(z[1][i]) <= err_tol)) ||
112  // End condition when z is discontinuous
113  ((mode == DISCONTINUOUS) && (h <= err_tol))
114  )
115  {
116  events[i] = found_event = true;
117  }
118  // There is an event in (0,h)
119  else
120  {
121  if (mode == INTERPOLATE)
122  {
123  double tcandidate = z[0][i]*h/(z[0][i]-z[1][i]);
124  // Don't let the step size go to zero
125  if (tcandidate < h/4.0) tcandidate = h/4.0;
126  if (tcandidate < tguess) tguess = tcandidate;
127  }
128  event_in_interval = true;
129  }
130  }
131  }
132  // Guess at a new h and calculate qend for that time
133  if (event_in_interval)
134  {
135  if (mode == BISECTION || mode == DISCONTINUOUS) h /= 2.0;
136  else h = tguess;
137  for (int i = 0; i < this->sys->numVars(); i++)
138  qend[i] = qstart[i];
139  solver->advance(qend,h);
140  }
141  else return found_event;
142  }
143  // Will never reach this line
144  return false;
145 }
146 
150 template <typename X>
152  public event_locator_impl<X>
153 {
154  public:
155  bisection_event_locator(ode_system<X>* sys, double err_tol):
157  sys,err_tol,event_locator_impl<X>::BISECTION){}
158 };
159 
163 template <typename X>
165  public event_locator_impl<X>
166 {
167  public:
168  linear_event_locator(ode_system<X>* sys, double err_tol):
171 };
172 
176 template <typename X>
178  public event_locator_impl<X>
179 {
180  public:
181  discontinuous_event_locator(ode_system<X>* sys, double err_tol):
184 };
185 
190 template <typename X> class null_event_locator:
191  public event_locator<X>
192 {
193  public:
195  ~null_event_locator(){}
196  bool find_events(bool*, const double*, double*, ode_solver<X>*, double&)
197  {
198  return false;
199  }
200 };
201 
202 } // end of namespace
203 
204 #endif
int numEvents() const
Get the number of state events.
Definition: adevs_hybrid.h:54
Mode
Definition: adevs_event_locators.h:54
bool find_events(bool *, const double *, double *, ode_solver< X > *, double &)
Definition: adevs_event_locators.h:196
Definition: adevs_event_locators.h:164
Definition: adevs_event_locators.h:44
Definition: adevs_event_locators.h:190
Definition: adevs_hybrid.h:406
Definition: adevs_event_locators.h:151
int numVars() const
Get the number of state variables.
Definition: adevs_hybrid.h:52
Definition: adevs_fmi.h:57
virtual void advance(double *q, double h)=0
Definition: adevs_event_locators.h:177
bool find_events(bool *events, const double *qstart, double *qend, ode_solver< X > *solver, double &h)
Definition: adevs_event_locators.h:90
virtual void state_event_func(const double *q, double *z)=0
Compute the state event functions for state q and put them in z.
Definition: adevs_hybrid.h:45
Definition: adevs_hybrid.h:377