Cantera  2.0
ResidEval.h
Go to the documentation of this file.
1 /**
2  * @file ResidEval.h
3  */
4
5 // Copyright 2006 California Institute of Technology
6
7 #ifndef CT_RESIDEVAL_H
8 #define CT_RESIDEVAL_H
9
10 #include "cantera/base/ct_defs.h"
12
13 #include <cstdio>
14
15 namespace Cantera
16 {
17
18 const int c_NONE = 0;
19 const int c_GE_ZERO = 1;
20 const int c_GT_ZERO = 2;
21 const int c_LE_ZERO = -1;
22 const int c_LT_ZERO = -2;
23
24 /**
25  * Virtual base class for DAE residual function evaluators.
26  * Classes derived from ResidEval evaluate the residual function
27  * \f[
28  * \vec{F}(t,\vec{y}, \vec{y^\prime})
29  * \f]
30  * The DAE solver attempts to find a solution y(t) such that F = 0.
31  * @ingroup DAE_Group
32  */
33 class ResidEval
34 {
35 public:
36  ResidEval() {}
37  virtual ~ResidEval() {}
38
39  /**
40  * Constrain solution component k. Possible values for
41  * 'flag' are:
42  * - c_NONE no constraint
43  * - c_GE_ZERO >= 0
44  * - c_GT_ZERO > 0
45  * - c_LE_ZERO <= 0
46  * - c_LT_ZERO < 0
47  */
48  virtual void constrain(const int k, const int flag) {
49  m_constrain[k] = flag;
50  }
51  int constraint(const int k) const {
52  std::map<int,int>::const_iterator i = m_constrain.find(k);
53  if (i != m_constrain.end()) {
54  return i->second;
55  }
56  return c_NONE;
57  }
58
59  //! Initialization function
60  virtual void initSizes() {
61  int neq = nEquations();
62  m_alg.resize(neq, 0);
63  }
64
65  /**
66  * Specify that solution component k is purely algebraic -
67  * that is, the derivative of this component does not appear
68  * in the residual function.
69  */
70  virtual void setAlgebraic(const int k) {
71  if ((int) m_alg.size() < (k+1)) {
72  initSizes();
73  }
74  m_alg[k] = 1;
75  }
76
77  virtual bool isAlgebraic(const int k) {
78  return (m_alg[k] == 1);
79  }
80
81
82  /**
83  * Evaluate the residual function. Called by the
84  * integrator.
85  * @param t time. (input)
86  * @param y solution vector. (input)
87  * @param ydot rate of change of solution vector. (input)
88  * @param r residual vector (output)
89  */
90  virtual int eval(const doublereal t, const doublereal* const y,
91  const doublereal* const ydot,
92  doublereal* const r) {
93  throw CanteraError("ResidEval::eval()", "base class called");
94  }
95
96  virtual int evalSS(const doublereal t, const doublereal* const y,
97  doublereal* const r) {
98  return eval(t, y, 0, r);
99  }
100
101  virtual int evalSimpleTD(const doublereal t, const doublereal* const y,
102  const doublereal* const yold, doublereal deltaT,
103  doublereal* const r) {
104  int nn = nEquations();
105  vector_fp ydot(nn);
106  for (int i = 0; i < nn; i++) {
107  ydot[i] = (y[i] - yold[i]) / deltaT;
108  }
109  return eval(t, y, DATA_PTR(ydot), r);
110  }
111
112  /**
113  * Fill the solution and derivative vectors with the initial
114  * conditions at initial time t0.
115  * @return 1 Everything is fine
116  * 0 or neg Something went wrong
117  */
118  virtual int getInitialConditions(const doublereal t0, doublereal* const y,
119  doublereal* const ydot) {
120  initSizes();
121  throw CanteraError("ResidEval::GetInitialConditions()", "base class called");
122  return 1;
123  }
124
125  //! Return the number of equations in the equation system
126  virtual int nEquations() const = 0;
127
128
129  //! Write out to a file or to standard output the current solution
130  /*!
131  * ievent is a description of the event that caused this
132  * function to be called.
133  */
134  virtual void writeSolution(int ievent, const double time,
135  const double deltaT,
136  const int time_step_num,
137  const double* y, const double* ydot) {
138  int k;
139  printf("ResidEval::writeSolution\n");
140  printf(" Time = %g, ievent = %d, deltaT = %g\n", time, ievent, deltaT);
141  if (ydot) {
142  printf(" k y[] ydot[]\n");
143  for (k = 0; k < nEquations(); k++) {
144  printf("%d %g %g\n", k, y[k], ydot[k]);
145  }
146  } else {
147  printf(" k y[]\n");
148  for (k = 0; k < nEquations(); k++) {
149  printf("%d %g \n", k, y[k]);
150  }
151  }
152  }
153
154  //! Return the number of parameters in the calculation
155  /*!
156  * This is the number of parameters in the sensitivity calculation. We have
157  * set this to zero and have included it for later expansion
158  */
159  int nparams() const {
160  return 0;
161  }
162
163 protected:
164
165  //! Mapping vector that stores whether a degree of freedom is a DAE or not
166  /*!
167  * The first index is the equation number. The second index is 1 if it is a DAE,
168  * and zero if it is not.
169  */
170  std::vector<int> m_alg;
171  std::map<int, int> m_constrain;
172
173 private:
174
175 };
176
177 }
178
179 #endif