Cantera  2.1.2
CVodeInt.cpp
Go to the documentation of this file.
1 /**
2  * @file CVodeInt.cpp
3  */
4 
5 // Copyright 2001 California Institute of Technology
6 
7 #include "CVodeInt.h"
8 using namespace std;
9 
10 // cvode includes
11 #include "../../ext/cvode/include/llnltyps.h"
12 #include "../../ext/cvode/include/llnlmath.h"
13 #include "../../ext/cvode/include/cvode.h"
14 #include "../../ext/cvode/include/cvdense.h"
15 #include "../../ext/cvode/include/cvdiag.h"
16 #include "../../ext/cvode/include/cvspgmr.h"
17 #include "../../ext/cvode/include/cvode.h"
18 
20 
21 extern "C" {
22 
23  /**
24  * Function called by cvode to evaluate ydot given y. The cvode
25  * integrator allows passing in a void* pointer to access
26  * external data. This pointer is cast to a pointer to a instance
27  * of class FuncEval. The equations to be integrated should be
28  * specified by deriving a class from FuncEval that evaluates the
29  * desired equations.
30  * @ingroup odeGroup
31  */
32  static void cvode_rhs(integer N, real t, N_Vector y, N_Vector ydot,
33  void* f_data)
34  {
35  double* ydata = N_VDATA(y);
36  double* ydotdata = N_VDATA(ydot);
38  f->eval(t, ydata, ydotdata, NULL);
39  }
40 
41  /**
42  * Function called by cvode to evaluate the Jacobian matrix.
43  * (temporary)
44  * @ingroup odeGroup
45  */
46  static void cvode_jac(integer N, DenseMat J, RhsFn f, void* f_data,
47  real t, N_Vector y, N_Vector fy, N_Vector ewt, real h, real uround,
48  void* jac_data, long int* nfePtr, N_Vector vtemp1, N_Vector vtemp2,
49  N_Vector vtemp3)
50  {
51  // get pointers to start of data
52  double* ydata = N_VDATA(y);
53  double* fydata = N_VDATA(fy);
54  double* ewtdata = N_VDATA(ewt);
55  double* ydot = N_VDATA(vtemp1);
56 
57  Cantera::FuncEval* func = (Cantera::FuncEval*)f_data;
58 
59  int i,j;
60  double* col_j;
61  double ysave, dy;
62  for (j=0; j < N; j++) {
63  col_j = (J->data)[j];
64  ysave = ydata[j];
65  dy = 1.0/ewtdata[j];
66  ydata[j] = ysave + dy;
67  dy = ydata[j] - ysave;
68  func->eval(t, ydata, ydot, NULL);
69  for (i=0; i < N; i++) {
70  col_j[i] = (ydot[i] - fydata[i])/dy;
71  }
72  ydata[j] = ysave;
73  }
74  }
75 }
76 
77 namespace Cantera
78 {
79 CVodeInt::CVodeInt() : m_neq(0),
80  m_cvode_mem(0),
81  m_t0(0.0),
82  m_y(0),
83  m_abstol(0),
84  m_type(DENSE+NOJAC),
85  m_itol(0),
86  m_method(BDF),
87  m_iter(NEWTON),
88  m_maxord(0),
89  m_reltol(1.e-9),
90  m_abstols(1.e-15),
91  m_nabs(0),
92  m_hmax(0.0),
93  m_maxsteps(20000)
94 {
95  m_ropt.resize(OPT_SIZE,0.0);
96  m_iopt = new long[OPT_SIZE];
97  fill(m_iopt, m_iopt+OPT_SIZE,0);
98 }
99 
100 CVodeInt::~CVodeInt()
101 {
102  if (m_cvode_mem) {
103  CVodeFree(m_cvode_mem);
104  }
105  if (m_y) {
106  N_VFree(m_y);
107  }
108  if (m_abstol) {
109  N_VFree(m_abstol);
110  }
111  delete[] m_iopt;
112 }
113 
114 double& CVodeInt::solution(size_t k)
115 {
116  return N_VIth(m_y, int(k));
117 }
119 {
120  return N_VDATA(m_y);
121 }
122 
123 void CVodeInt::setTolerances(double reltol, size_t n, double* abstol)
124 {
125  m_itol = 1;
126  m_nabs = int(n);
127  if (m_nabs != m_neq) {
128  if (m_abstol) {
129  N_VFree(m_abstol);
130  }
131  m_abstol = N_VNew(m_nabs, 0);
132  }
133  for (int i=0; i<m_nabs; i++) {
134  N_VIth(m_abstol, i) = abstol[i];
135  }
136  m_reltol = reltol;
137 }
138 
139 void CVodeInt::setTolerances(double reltol, double abstol)
140 {
141  m_itol = 0;
142  m_reltol = reltol;
143  m_abstols = abstol;
144 }
145 
146 void CVodeInt::setProblemType(int probtype)
147 {
148  m_type = probtype;
149 }
150 
152 {
153  if (t == BDF_Method) {
154  m_method = BDF;
155  } else if (t == Adams_Method) {
156  m_method = ADAMS;
157  } else {
158  throw CVodeErr("unknown method");
159  }
160 }
161 
162 void CVodeInt::setMaxStepSize(doublereal hmax)
163 {
164  m_hmax = hmax;
165  m_ropt[HMAX] = hmax;
166 }
167 
168 void CVodeInt::setMinStepSize(doublereal hmin)
169 {
170  m_hmin = hmin;
171  m_ropt[HMIN] = hmin;
172 }
173 
174 void CVodeInt::setMaxSteps(int nmax)
175 {
176  m_maxsteps = nmax;
177  m_iopt[MXSTEP] = m_maxsteps;
178 }
179 
181 {
182  if (t == Newton_Iter) {
183  m_iter = NEWTON;
184  } else if (t == Functional_Iter) {
185  m_iter = FUNCTIONAL;
186  } else {
187  throw CVodeErr("unknown iterator");
188  }
189 }
190 
191 void CVodeInt::initialize(double t0, FuncEval& func)
192 {
193  m_neq = int(func.neq());
194  m_t0 = t0;
195 
196  if (m_y) {
197  N_VFree(m_y); // free solution vector if already allocated
198  }
199  m_y = N_VNew(m_neq, 0); // allocate solution vector
200  // check abs tolerance array size
201  if (m_itol == 1 && m_nabs < m_neq) {
202  throw CVodeErr("not enough absolute tolerance values specified.");
203  }
204  func.getInitialConditions(m_t0, m_neq, N_VDATA(m_y));
205 
206  // set options
207  m_iopt[MXSTEP] = m_maxsteps;
208  m_iopt[MAXORD] = m_maxord;
209  m_ropt[HMAX] = m_hmax;
210 
211  if (m_cvode_mem) {
212  CVodeFree(m_cvode_mem);
213  }
214 
215  // pass a pointer to func in m_data
216  m_data = (void*)&func;
217 
218  if (m_itol) {
219  m_cvode_mem = CVodeMalloc(m_neq, cvode_rhs, m_t0, m_y, m_method,
220  m_iter, m_itol, &m_reltol,
221  m_abstol, m_data, NULL, 1, m_iopt,
222  DATA_PTR(m_ropt), NULL);
223  } else {
224  m_cvode_mem = CVodeMalloc(m_neq, cvode_rhs, m_t0, m_y, m_method,
225  m_iter, m_itol, &m_reltol,
226  &m_abstols, m_data, NULL, 1, m_iopt,
227  DATA_PTR(m_ropt), NULL);
228  }
229 
230  if (!m_cvode_mem) {
231  throw CVodeErr("CVodeMalloc failed.");
232  }
233 
234  if (m_type == DENSE + NOJAC) {
235  CVDense(m_cvode_mem, NULL, NULL);
236  } else if (m_type == DENSE + JAC) {
237  CVDense(m_cvode_mem, cvode_jac, NULL);
238  } else if (m_type == DIAG) {
239  CVDiag(m_cvode_mem);
240  } else if (m_type == GMRES) {
241  CVSpgmr(m_cvode_mem, NONE, MODIFIED_GS, 0, 0.0,
242  NULL, NULL, NULL);
243  } else {
244  throw CVodeErr("unsupported option");
245  }
246 }
247 
248 void CVodeInt::reinitialize(double t0, FuncEval& func)
249 {
250  m_t0 = t0;
251  func.getInitialConditions(m_t0, m_neq, N_VDATA(m_y));
252 
253  // set options
254  m_iopt[MXSTEP] = m_maxsteps;
255  m_iopt[MAXORD] = m_maxord;
256  m_ropt[HMAX] = m_hmax;
257 
258  //if (m_cvode_mem) CVodeFree(m_cvode_mem);
259 
260  // pass a pointer to func in m_data
261  m_data = (void*)&func;
262  int result;
263  if (m_itol) {
264  result = CVReInit(m_cvode_mem, cvode_rhs, m_t0, m_y, m_method,
265  m_iter, m_itol, &m_reltol,
266  m_abstol, m_data, NULL, 1, m_iopt,
267  DATA_PTR(m_ropt), NULL);
268  } else {
269  result = CVReInit(m_cvode_mem, cvode_rhs, m_t0, m_y, m_method,
270  m_iter, m_itol, &m_reltol,
271  &m_abstols, m_data, NULL, 1, m_iopt,
272  DATA_PTR(m_ropt), NULL);
273  }
274 
275  if (result != 0) {
276  throw CVodeErr("CVReInit failed.");
277  }
278 
279  if (m_type == DENSE + NOJAC) {
280  CVDense(m_cvode_mem, NULL, NULL);
281  } else if (m_type == DENSE + JAC) {
282  CVDense(m_cvode_mem, cvode_jac, NULL);
283  } else if (m_type == DIAG) {
284  CVDiag(m_cvode_mem);
285  } else if (m_type == GMRES) {
286  CVSpgmr(m_cvode_mem, NONE, MODIFIED_GS, 0, 0.0,
287  NULL, NULL, NULL);
288  } else {
289  throw CVodeErr("unsupported option");
290  }
291 }
292 
293 void CVodeInt::integrate(double tout)
294 {
295  double t;
296  int flag;
297  flag = CVode(m_cvode_mem, tout, m_y, &t, NORMAL);
298  if (flag != SUCCESS) {
299  throw CVodeErr(" CVode error encountered. Error code: " + int2str(flag));
300  }
301 }
302 
303 double CVodeInt::step(double tout)
304 {
305  double t;
306  int flag;
307  flag = CVode(m_cvode_mem, tout, m_y, &t, ONE_STEP);
308  if (flag != SUCCESS) {
309  throw CVodeErr(" CVode error encountered. Error code: " + int2str(flag));
310  }
311  return t;
312 }
313 
314 int CVodeInt::nEvals() const
315 {
316  return m_iopt[NFE];
317 }
318 }
Backward Differentiation.
Definition: Integrator.h:33
Exception thrown when a CVODE error is encountered.
Definition: CVodeInt.h:21
std::string int2str(const int n, const std::string &fmt)
Convert an int to a string using a format converter.
Definition: stringUtils.cpp:40
static void cvode_jac(integer N, DenseMat J, RhsFn f, void *f_data, real t, N_Vector y, N_Vector fy, N_Vector ewt, real h, real uround, void *jac_data, long int *nfePtr, N_Vector vtemp1, N_Vector vtemp2, N_Vector vtemp3)
Function called by cvode to evaluate the Jacobian matrix.
Definition: CVodeInt.cpp:46
virtual void eval(double t, double *y, double *ydot, double *p)=0
Evaluate the right-hand-side function.
virtual double * solution()
The current value of the solution of the system of equations.
Definition: CVodeInt.cpp:118
virtual void setMethod(MethodType t)
Set the solution method.
Definition: CVodeInt.cpp:151
virtual void integrate(double tout)
Integrate the system of equations.
Definition: CVodeInt.cpp:293
Functional Iteration.
Definition: Integrator.h:45
virtual size_t neq()=0
Number of equations.
virtual void setIterator(IterType t)
Set the linear iterator.
Definition: CVodeInt.cpp:180
virtual int nEvals() const
The number of function evaluations.
Definition: CVodeInt.cpp:314
virtual void setProblemType(int probtype)
Set the problem type.
Definition: CVodeInt.cpp:146
static void cvode_rhs(integer N, real t, N_Vector y, N_Vector ydot, void *f_data)
Function called by cvode to evaluate ydot given y.
Definition: CVodeInt.cpp:32
virtual void setTolerances(double reltol, size_t n, double *abstol)
Set or reset the number of equations.
Definition: CVodeInt.cpp:123
virtual void initialize(double t0, FuncEval &func)
Initialize the integrator for a new problem.
Definition: CVodeInt.cpp:191
virtual void setMaxStepSize(double hmax)
Set the maximum step size.
Definition: CVodeInt.cpp:162
IterType
Specifies the method used for iteration.
Definition: Integrator.h:41
virtual void setMinStepSize(double hmin)
Set the minimum step size.
Definition: CVodeInt.cpp:168
Newton Iteration.
Definition: Integrator.h:43
Contains declarations for string manipulation functions within Cantera.
#define DATA_PTR(vec)
Creates a pointer to the start of the raw data for a vector.
Definition: ct_defs.h:36
Virtual base class for ODE right-hand-side function evaluators.
Definition: FuncEval.h:23
virtual doublereal step(double tout)
Integrate the system of equations.
Definition: CVodeInt.cpp:303
MethodType
Specifies the method used to integrate the system of equations.
Definition: Integrator.h:32
virtual void getInitialConditions(double t0, size_t leny, double *y)=0
Fill the solution vector with the initial conditions at initial time t0.