Cantera  2.5.1
ResidEval.h
Go to the documentation of this file.
1 //! @file ResidEval.h
2 
3 // This file is part of Cantera. See License.txt in the top-level directory or
4 // at https://cantera.org/license.txt for license and copyright information.
5 
6 #ifndef CT_RESIDEVAL_H
7 #define CT_RESIDEVAL_H
8 
9 #include "cantera/base/ct_defs.h"
11 #include "cantera/base/utilities.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  return getValue(m_constrain, k, c_NONE);
53  }
54 
55  //! Initialization function
56  virtual void initSizes() {
57  int neq = nEquations();
58  m_alg.resize(neq, 0);
59  }
60 
61  /**
62  * Specify that solution component k is purely algebraic - that is, the
63  * derivative of this component does not appear in the residual function.
64  */
65  virtual void setAlgebraic(const int k) {
66  if ((int) m_alg.size() < (k+1)) {
67  initSizes();
68  }
69  m_alg[k] = 1;
70  }
71 
72  virtual bool isAlgebraic(const int k) {
73  return (m_alg[k] == 1);
74  }
75 
76  /**
77  * Evaluate the residual function. Called by the integrator.
78  * @param t time. (input)
79  * @param y solution vector. (input)
80  * @param ydot rate of change of solution vector. (input)
81  * @param r residual vector (output)
82  */
83  virtual int eval(const doublereal t, const doublereal* const y,
84  const doublereal* const ydot,
85  doublereal* const r) {
86  throw NotImplementedError("ResidEval::eval");
87  }
88 
89  virtual int evalSS(const doublereal t, const doublereal* const y,
90  doublereal* const r) {
91  return eval(t, y, 0, r);
92  }
93 
94  virtual int evalSimpleTD(const doublereal t, const doublereal* const y,
95  const doublereal* const yold, doublereal deltaT,
96  doublereal* const r) {
97  int nn = nEquations();
98  vector_fp ydot(nn);
99  for (int i = 0; i < nn; i++) {
100  ydot[i] = (y[i] - yold[i]) / deltaT;
101  }
102  return eval(t, y, ydot.data(), r);
103  }
104 
105  //! Fill in the initial conditions
106  /*!
107  * Values for both the solution and the value of ydot may be provided.
108  *
109  * @param[in] t0 Time
110  * @param[out] y Solution vector
111  * @param[out] ydot Rate of change of solution vector.
112  *
113  * @returns a flag to indicate that operation is successful.
114  * 1 Means a successful operation
115  * -0 or neg value Means an unsuccessful operation
116  */
117  virtual int getInitialConditions(const doublereal t0, doublereal* const y,
118  doublereal* const ydot) {
119  initSizes();
120  throw NotImplementedError("ResidEval::GetInitialConditions");
121  return 1;
122  }
123 
124  //! Return the number of equations in the equation system
125  virtual int nEquations() const = 0;
126 
127  //! Write out to a file or to standard output the current solution
128  /*!
129  * ievent is a description of the event that caused this function to be
130  * called.
131  */
132  virtual void writeSolution(int ievent, const double time,
133  const double deltaT,
134  const int time_step_num,
135  const double* y, const double* ydot) {
136  int k;
137  writelog("ResidEval::writeSolution\n");
138  writelogf(" Time = %g, ievent = %d, deltaT = %g\n", time, ievent, deltaT);
139  if (ydot) {
140  writelogf(" k y[] ydot[]\n");
141  for (k = 0; k < nEquations(); k++) {
142  writelogf("%d %g %g\n", k, y[k], ydot[k]);
143  }
144  } else {
145  writelogf(" k y[]\n");
146  for (k = 0; k < nEquations(); k++) {
147  writelogf("%d %g \n", k, y[k]);
148  }
149  }
150  }
151 
152  //! Return the number of parameters in the calculation
153  /*!
154  * This is the number of parameters in the sensitivity calculation. We have
155  * set this to zero and have included it for later expansion
156  */
157  int nparams() const {
158  return 0;
159  }
160 
161 protected:
162  //! Mapping vector that stores whether a degree of freedom is a DAE or not
163  /*!
164  * The first index is the equation number. The second index is 1 if it is a
165  * DAE, and zero if it is not.
166  */
168  std::map<int, int> m_constrain;
169 };
170 
171 }
172 
173 #endif
An error indicating that an unimplemented function has been called.
Definition: ctexceptions.h:187
Virtual base class for DAE residual function evaluators.
Definition: ResidEval.h:34
virtual int eval(const doublereal t, const doublereal *const y, const doublereal *const ydot, doublereal *const r)
Evaluate the residual function.
Definition: ResidEval.h:83
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:132
virtual int getInitialConditions(const doublereal t0, doublereal *const y, doublereal *const ydot)
Fill in the initial conditions.
Definition: ResidEval.h:117
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:65
virtual int nEquations() const =0
Return the number of equations in the equation system.
int nparams() const
Return the number of parameters in the calculation.
Definition: ResidEval.h:157
vector_int m_alg
Mapping vector that stores whether a degree of freedom is a DAE or not.
Definition: ResidEval.h:167
virtual void constrain(const int k, const int flag)
Constrain solution component k.
Definition: ResidEval.h:48
virtual void initSizes()
Initialization function.
Definition: ResidEval.h:56
This file contains definitions of terms that are used in internal routines and are unlikely to need m...
Definitions for the classes that are thrown when Cantera experiences an error condition (also contain...
std::vector< int > vector_int
Vector of ints.
Definition: ct_defs.h:182
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:180
void writelog(const std::string &fmt, const Args &... args)
Write a formatted message to the screen.
Definition: global.h:158
void writelogf(const char *fmt, const Args &... args)
Write a formatted message to the screen.
Definition: global.h:179
Namespace for the Cantera kernel.
Definition: AnyMap.cpp:264
const U & getValue(const std::map< T, U > &m, const T &key, const U &default_val)
Const accessor for a value in a std::map.
Definition: utilities.h:528
Various templated functions that carry out common vector operations (see Templated Utility Functions)...