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