Cantera  2.1.2
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  * Evaluate the residual function. Called by the integrator.
83  * @param t time. (input)
84  * @param y solution vector. (input)
85  * @param ydot rate of change of solution vector. (input)
86  * @param r residual vector (output)
87  */
88  virtual int eval(const doublereal t, const doublereal* const y,
89  const doublereal* const ydot,
90  doublereal* const r) {
91  throw CanteraError("ResidEval::eval()", "base class called");
92  }
93 
94  virtual int evalSS(const doublereal t, const doublereal* const y,
95  doublereal* const r) {
96  return eval(t, y, 0, r);
97  }
98 
99  virtual int evalSimpleTD(const doublereal t, const doublereal* const y,
100  const doublereal* const yold, doublereal deltaT,
101  doublereal* const r) {
102  int nn = nEquations();
103  vector_fp ydot(nn);
104  for (int i = 0; i < nn; i++) {
105  ydot[i] = (y[i] - yold[i]) / deltaT;
106  }
107  return eval(t, y, DATA_PTR(ydot), r);
108  }
109 
110  //! Fill in the initial conditions
111  /*!
112  * Values for both the solution and the value of ydot may be provided.
113  *
114  * @param[in] t0 Time
115  * @param[out] y Solution vector
116  * @param[out] ydot Rate of change of solution vector.
117  *
118  * @return Returns a flag to indicate that operation is successful.
119  * 1 Means a successful operation
120  * -0 or neg value Means an unsuccessful operation
121  */
122  virtual int getInitialConditions(const doublereal t0, doublereal* const y,
123  doublereal* const ydot) {
124  initSizes();
125  throw CanteraError("ResidEval::GetInitialConditions()", "base class called");
126  return 1;
127  }
128 
129  //! Return the number of equations in the equation system
130  virtual int nEquations() const = 0;
131 
132 
133  //! Write out to a file or to standard output the current solution
134  /*!
135  * ievent is a description of the event that caused this
136  * function to be called.
137  */
138  virtual void writeSolution(int ievent, const double time,
139  const double deltaT,
140  const int time_step_num,
141  const double* y, const double* ydot) {
142  int k;
143  printf("ResidEval::writeSolution\n");
144  printf(" Time = %g, ievent = %d, deltaT = %g\n", time, ievent, deltaT);
145  if (ydot) {
146  printf(" k y[] ydot[]\n");
147  for (k = 0; k < nEquations(); k++) {
148  printf("%d %g %g\n", k, y[k], ydot[k]);
149  }
150  } else {
151  printf(" k y[]\n");
152  for (k = 0; k < nEquations(); k++) {
153  printf("%d %g \n", k, y[k]);
154  }
155  }
156  }
157 
158  //! Return the number of parameters in the calculation
159  /*!
160  * This is the number of parameters in the sensitivity calculation. We have
161  * set this to zero and have included it for later expansion
162  */
163  int nparams() const {
164  return 0;
165  }
166 
167 protected:
168 
169  //! Mapping vector that stores whether a degree of freedom is a DAE or not
170  /*!
171  * The first index is the equation number. The second index is 1 if it is a DAE,
172  * and zero if it is not.
173  */
174  std::vector<int> m_alg;
175  std::map<int, int> m_constrain;
176 
177 private:
178 
179 };
180 
181 }
182 
183 #endif
virtual void initSizes()
Initialization function.
Definition: ResidEval.h:60
This file contains definitions of terms that are used in internal routines and are unlikely to need m...
virtual void writeSolution(int ievent, const double time, const double deltaT, const int time_step_num, const double *y, const double *ydot)
Write out to a file or to standard output the current solution.
Definition: ResidEval.h:138
virtual int getInitialConditions(const doublereal t0, doublereal *const y, doublereal *const ydot)
Fill in the initial conditions.
Definition: ResidEval.h:122
Virtual base class for DAE residual function evaluators.
Definition: ResidEval.h:33
int nparams() const
Return the number of parameters in the calculation.
Definition: ResidEval.h:163
Base class for exceptions thrown by Cantera classes.
Definition: ctexceptions.h:68
virtual int nEquations() const =0
Return the number of equations in the equation system.
std::vector< double > vector_fp
Turn on the use of stl vectors for the basic array type within cantera Vector of doubles.
Definition: ct_defs.h:165
virtual void constrain(const int k, const int flag)
Constrain solution component k.
Definition: ResidEval.h:48
virtual void setAlgebraic(const int k)
Specify that solution component k is purely algebraic - that is, the derivative of this component doe...
Definition: ResidEval.h:70
#define DATA_PTR(vec)
Creates a pointer to the start of the raw data for a vector.
Definition: ct_defs.h:36
virtual int eval(const doublereal t, const doublereal *const y, const doublereal *const ydot, doublereal *const r)
Evaluate the residual function.
Definition: ResidEval.h:88
std::vector< int > m_alg
Mapping vector that stores whether a degree of freedom is a DAE or not.
Definition: ResidEval.h:174
Definitions for the classes that are thrown when Cantera experiences an error condition (also contain...