Cantera  2.4.0
CVodesIntegrator.cpp
Go to the documentation of this file.
1 //! @file CVodesIntegrator.cpp
2 
3 // This file is part of Cantera. See License.txt in the top-level directory or
4 // at http://www.cantera.org/license.txt for license and copyright information.
5 
8 
9 #include <iostream>
10 using namespace std;
11 
12 #include "sundials/sundials_types.h"
13 #include "sundials/sundials_math.h"
14 #include "sundials/sundials_nvector.h"
15 #include "nvector/nvector_serial.h"
16 #include "cvodes/cvodes.h"
17 #if CT_SUNDIALS_VERSION >= 30
18  #if CT_SUNDIALS_USE_LAPACK
19  #include "sunlinsol/sunlinsol_lapackdense.h"
20  #include "sunlinsol/sunlinsol_lapackband.h"
21  #else
22  #include "sunlinsol/sunlinsol_dense.h"
23  #include "sunlinsol/sunlinsol_band.h"
24  #endif
25  #include "sunlinsol/sunlinsol_spgmr.h"
26  #include "cvodes/cvodes_direct.h"
27  #include "cvodes/cvodes_diag.h"
28  #include "cvodes/cvodes_spils.h"
29 #else
30  #if CT_SUNDIALS_USE_LAPACK
31  #include "cvodes/cvodes_lapack.h"
32  #else
33  #include "cvodes/cvodes_dense.h"
34  #include "cvodes/cvodes_band.h"
35  #endif
36  #include "cvodes/cvodes_diag.h"
37  #include "cvodes/cvodes_spgmr.h"
38 #endif
39 
40 #define CV_SS 1
41 #define CV_SV 2
42 
43 #if CT_SUNDIALS_VERSION < 25
44 typedef int sd_size_t;
45 #else
46 typedef long int sd_size_t;
47 #endif
48 
49 namespace Cantera
50 {
51 
52 extern "C" {
53  /**
54  * Function called by cvodes to evaluate ydot given y. The CVODE integrator
55  * allows passing in a void* pointer to access external data. This pointer
56  * is cast to a pointer to a instance of class FuncEval. The equations to be
57  * integrated should be specified by deriving a class from FuncEval that
58  * evaluates the desired equations.
59  * @ingroup odeGroup
60  */
61  static int cvodes_rhs(realtype t, N_Vector y, N_Vector ydot, void* f_data)
62  {
63  FuncEval* f = (FuncEval*) f_data;
64  return f->eval_nothrow(t, NV_DATA_S(y), NV_DATA_S(ydot));
65  }
66 
67  //! Function called by CVodes when an error is encountered instead of
68  //! writing to stdout. Here, save the error message provided by CVodes so
69  //! that it can be included in the subsequently raised CanteraError.
70  static void cvodes_err(int error_code, const char* module,
71  const char* function, char* msg, void* eh_data)
72  {
73  CVodesIntegrator* integrator = (CVodesIntegrator*) eh_data;
74  integrator->m_error_message = msg;
75  integrator->m_error_message += "\n";
76  }
77 }
78 
79 CVodesIntegrator::CVodesIntegrator() :
80  m_neq(0),
81  m_cvode_mem(0),
82  m_linsol(0),
83  m_linsol_matrix(0),
84  m_func(0),
85  m_t0(0.0),
86  m_y(0),
87  m_abstol(0),
88  m_type(DENSE+NOJAC),
89  m_itol(CV_SS),
90  m_method(CV_BDF),
91  m_iter(CV_NEWTON),
92  m_maxord(0),
93  m_reltol(1.e-9),
94  m_abstols(1.e-15),
95  m_reltolsens(1.0e-5),
96  m_abstolsens(1.0e-4),
97  m_nabs(0),
98  m_hmax(0.0),
99  m_hmin(0.0),
100  m_maxsteps(20000),
101  m_maxErrTestFails(0),
102  m_yS(nullptr),
103  m_np(0),
104  m_mupper(0), m_mlower(0),
105  m_sens_ok(false)
106 {
107 }
108 
109 CVodesIntegrator::~CVodesIntegrator()
110 {
111  if (m_cvode_mem) {
112  if (m_np > 0) {
113  CVodeSensFree(m_cvode_mem);
114  }
115  CVodeFree(&m_cvode_mem);
116  }
117 
118  #if CT_SUNDIALS_VERSION >= 30
119  SUNLinSolFree((SUNLinearSolver) m_linsol);
120  SUNMatDestroy((SUNMatrix) m_linsol_matrix);
121  #endif
122 
123  if (m_y) {
124  N_VDestroy_Serial(m_y);
125  }
126  if (m_abstol) {
127  N_VDestroy_Serial(m_abstol);
128  }
129  if (m_yS) {
130  N_VDestroyVectorArray_Serial(m_yS, static_cast<sd_size_t>(m_np));
131  }
132 }
133 
134 double& CVodesIntegrator::solution(size_t k)
135 {
136  return NV_Ith_S(m_y, k);
137 }
138 
140 {
141  return NV_DATA_S(m_y);
142 }
143 
144 void CVodesIntegrator::setTolerances(double reltol, size_t n, double* abstol)
145 {
146  m_itol = CV_SV;
147  m_nabs = n;
148  if (n != m_neq) {
149  if (m_abstol) {
150  N_VDestroy_Serial(m_abstol);
151  }
152  m_abstol = N_VNew_Serial(static_cast<sd_size_t>(n));
153  }
154  for (size_t i=0; i<n; i++) {
155  NV_Ith_S(m_abstol, i) = abstol[i];
156  }
157  m_reltol = reltol;
158 }
159 
160 void CVodesIntegrator::setTolerances(double reltol, double abstol)
161 {
162  m_itol = CV_SS;
163  m_reltol = reltol;
164  m_abstols = abstol;
165 }
166 
167 void CVodesIntegrator::setSensitivityTolerances(double reltol, double abstol)
168 {
169  m_reltolsens = reltol;
170  m_abstolsens = abstol;
171 }
172 
174 {
175  m_type = probtype;
176 }
177 
179 {
180  if (t == BDF_Method) {
181  m_method = CV_BDF;
182  } else if (t == Adams_Method) {
183  m_method = CV_ADAMS;
184  } else {
185  throw CanteraError("CVodesIntegrator::setMethod", "unknown method");
186  }
187 }
188 
189 void CVodesIntegrator::setMaxStepSize(doublereal hmax)
190 {
191  m_hmax = hmax;
192  if (m_cvode_mem) {
193  CVodeSetMaxStep(m_cvode_mem, hmax);
194  }
195 }
196 
197 void CVodesIntegrator::setMinStepSize(doublereal hmin)
198 {
199  m_hmin = hmin;
200  if (m_cvode_mem) {
201  CVodeSetMinStep(m_cvode_mem, hmin);
202  }
203 }
204 
205 void CVodesIntegrator::setMaxSteps(int nmax)
206 {
207  m_maxsteps = nmax;
208  if (m_cvode_mem) {
209  CVodeSetMaxNumSteps(m_cvode_mem, m_maxsteps);
210  }
211 }
212 
214 {
215  m_maxErrTestFails = n;
216  if (m_cvode_mem) {
217  CVodeSetMaxErrTestFails(m_cvode_mem, n);
218  }
219 }
220 
222 {
223  if (t == Newton_Iter) {
224  m_iter = CV_NEWTON;
225  } else if (t == Functional_Iter) {
226  m_iter = CV_FUNCTIONAL;
227  } else {
228  throw CanteraError("CVodesIntegrator::setIterator", "unknown iterator");
229  }
230 }
231 
232 void CVodesIntegrator::sensInit(double t0, FuncEval& func)
233 {
234  m_np = func.nparams();
235  m_sens_ok = false;
236 
237  N_Vector y = N_VNew_Serial(static_cast<sd_size_t>(func.neq()));
238  m_yS = N_VCloneVectorArray_Serial(static_cast<sd_size_t>(m_np), y);
239  for (size_t n = 0; n < m_np; n++) {
240  N_VConst(0.0, m_yS[n]);
241  }
242  N_VDestroy_Serial(y);
243 
244  int flag = CVodeSensInit(m_cvode_mem, static_cast<sd_size_t>(m_np),
245  CV_STAGGERED, CVSensRhsFn(0), m_yS);
246 
247  if (flag != CV_SUCCESS) {
248  throw CanteraError("CVodesIntegrator::sensInit", "Error in CVodeSensInit");
249  }
250  vector_fp atol(m_np);
251  for (size_t n = 0; n < m_np; n++) {
252  // This scaling factor is tuned so that reaction and species enthalpy
253  // sensitivities can be computed simultaneously with the same abstol.
254  atol[n] = m_abstolsens / func.m_paramScales[n];
255  }
256  flag = CVodeSensSStolerances(m_cvode_mem, m_reltolsens, atol.data());
257 }
258 
260 {
261  m_neq = func.neq();
262  m_t0 = t0;
263  m_time = t0;
264  m_func = &func;
265  func.clearErrors();
266 
267  if (m_y) {
268  N_VDestroy_Serial(m_y); // free solution vector if already allocated
269  }
270  m_y = N_VNew_Serial(static_cast<sd_size_t>(m_neq)); // allocate solution vector
271  N_VConst(0.0, m_y);
272  // check abs tolerance array size
273  if (m_itol == CV_SV && m_nabs < m_neq) {
274  throw CanteraError("CVodesIntegrator::initialize",
275  "not enough absolute tolerance values specified.");
276  }
277 
278  func.getState(NV_DATA_S(m_y));
279 
280  if (m_cvode_mem) {
281  CVodeFree(&m_cvode_mem);
282  }
283 
284  //! Specify the method and the iteration type. Cantera Defaults:
285  //! CV_BDF - Use BDF methods
286  //! CV_NEWTON - use Newton's method
287  m_cvode_mem = CVodeCreate(m_method, m_iter);
288  if (!m_cvode_mem) {
289  throw CanteraError("CVodesIntegrator::initialize",
290  "CVodeCreate failed.");
291  }
292 
293  int flag = CVodeInit(m_cvode_mem, cvodes_rhs, m_t0, m_y);
294  if (flag != CV_SUCCESS) {
295  if (flag == CV_MEM_FAIL) {
296  throw CanteraError("CVodesIntegrator::initialize",
297  "Memory allocation failed.");
298  } else if (flag == CV_ILL_INPUT) {
299  throw CanteraError("CVodesIntegrator::initialize",
300  "Illegal value for CVodeInit input argument.");
301  } else {
302  throw CanteraError("CVodesIntegrator::initialize",
303  "CVodeInit failed.");
304  }
305  }
306  CVodeSetErrHandlerFn(m_cvode_mem, &cvodes_err, this);
307 
308  if (m_itol == CV_SV) {
309  flag = CVodeSVtolerances(m_cvode_mem, m_reltol, m_abstol);
310  } else {
311  flag = CVodeSStolerances(m_cvode_mem, m_reltol, m_abstols);
312  }
313  if (flag != CV_SUCCESS) {
314  if (flag == CV_MEM_FAIL) {
315  throw CanteraError("CVodesIntegrator::initialize",
316  "Memory allocation failed.");
317  } else if (flag == CV_ILL_INPUT) {
318  throw CanteraError("CVodesIntegrator::initialize",
319  "Illegal value for CVodeInit input argument.");
320  } else {
321  throw CanteraError("CVodesIntegrator::initialize",
322  "CVodeInit failed.");
323  }
324  }
325 
326  flag = CVodeSetUserData(m_cvode_mem, &func);
327  if (flag != CV_SUCCESS) {
328  throw CanteraError("CVodesIntegrator::initialize",
329  "CVodeSetUserData failed.");
330  }
331  if (func.nparams() > 0) {
332  sensInit(t0, func);
333  flag = CVodeSetSensParams(m_cvode_mem, func.m_sens_params.data(),
334  func.m_paramScales.data(), NULL);
335  if (flag != CV_SUCCESS) {
336  throw CanteraError("CVodesIntegrator::initialize",
337  "CVodeSetSensParams failed.");
338  }
339  }
340  applyOptions();
341 }
342 
343 void CVodesIntegrator::reinitialize(double t0, FuncEval& func)
344 {
345  m_t0 = t0;
346  m_time = t0;
347  func.getState(NV_DATA_S(m_y));
348  m_func = &func;
349  func.clearErrors();
350 
351  int result = CVodeReInit(m_cvode_mem, m_t0, m_y);
352  if (result != CV_SUCCESS) {
353  throw CanteraError("CVodesIntegrator::reinitialize",
354  "CVodeReInit failed. result = {}", result);
355  }
356  applyOptions();
357 }
358 
360 {
361  if (m_type == DENSE + NOJAC) {
362  sd_size_t N = static_cast<sd_size_t>(m_neq);
363  #if CT_SUNDIALS_VERSION >= 30
364  SUNLinSolFree((SUNLinearSolver) m_linsol);
365  SUNMatDestroy((SUNMatrix) m_linsol_matrix);
366  m_linsol_matrix = SUNDenseMatrix(N, N);
367  #if CT_SUNDIALS_USE_LAPACK
368  m_linsol = SUNLapackDense(m_y, (SUNMatrix) m_linsol_matrix);
369  #else
370  m_linsol = SUNDenseLinearSolver(m_y, (SUNMatrix) m_linsol_matrix);
371  #endif
372  CVDlsSetLinearSolver(m_cvode_mem, (SUNLinearSolver) m_linsol,
373  (SUNMatrix) m_linsol_matrix);
374  #else
375  #if CT_SUNDIALS_USE_LAPACK
376  CVLapackDense(m_cvode_mem, N);
377  #else
378  CVDense(m_cvode_mem, N);
379  #endif
380  #endif
381  } else if (m_type == DIAG) {
382  CVDiag(m_cvode_mem);
383  } else if (m_type == GMRES) {
384  #if CT_SUNDIALS_VERSION >= 30
385  m_linsol = SUNSPGMR(m_y, PREC_NONE, 0);
386  CVSpilsSetLinearSolver(m_cvode_mem, (SUNLinearSolver) m_linsol);
387  #else
388  CVSpgmr(m_cvode_mem, PREC_NONE, 0);
389  #endif
390  } else if (m_type == BAND + NOJAC) {
391  sd_size_t N = static_cast<sd_size_t>(m_neq);
392  long int nu = m_mupper;
393  long int nl = m_mlower;
394  #if CT_SUNDIALS_VERSION >= 30
395  SUNLinSolFree((SUNLinearSolver) m_linsol);
396  SUNMatDestroy((SUNMatrix) m_linsol_matrix);
397  m_linsol_matrix = SUNBandMatrix(N, nu, nl, nu+nl);
398  #if CT_SUNDIALS_USE_LAPACK
399  m_linsol = SUNLapackBand(m_y, (SUNMatrix) m_linsol_matrix);
400  #else
401  m_linsol = SUNBandLinearSolver(m_y, (SUNMatrix) m_linsol_matrix);
402  #endif
403  CVDlsSetLinearSolver(m_cvode_mem, (SUNLinearSolver) m_linsol,
404  (SUNMatrix) m_linsol_matrix);
405  #else
406  #if CT_SUNDIALS_USE_LAPACK
407  CVLapackBand(m_cvode_mem, N, nu, nl);
408  #else
409  CVBand(m_cvode_mem, N, nu, nl);
410  #endif
411  #endif
412  } else {
413  throw CanteraError("CVodesIntegrator::applyOptions",
414  "unsupported option");
415  }
416 
417  if (m_maxord > 0) {
418  CVodeSetMaxOrd(m_cvode_mem, m_maxord);
419  }
420  if (m_maxsteps > 0) {
421  CVodeSetMaxNumSteps(m_cvode_mem, m_maxsteps);
422  }
423  if (m_hmax > 0) {
424  CVodeSetMaxStep(m_cvode_mem, m_hmax);
425  }
426  if (m_hmin > 0) {
427  CVodeSetMinStep(m_cvode_mem, m_hmin);
428  }
429  if (m_maxErrTestFails > 0) {
430  CVodeSetMaxErrTestFails(m_cvode_mem, m_maxErrTestFails);
431  }
432 }
433 
435 {
436  if (tout == m_time) {
437  return;
438  }
439  int flag = CVode(m_cvode_mem, tout, m_y, &m_time, CV_NORMAL);
440  if (flag != CV_SUCCESS) {
441  string f_errs = m_func->getErrors();
442  if (!f_errs.empty()) {
443  f_errs = "Exceptions caught during RHS evaluation:\n" + f_errs;
444  }
445  throw CanteraError("CVodesIntegrator::integrate",
446  "CVodes error encountered. Error code: {}\n{}\n"
447  "{}"
448  "Components with largest weighted error estimates:\n{}",
449  flag, m_error_message, f_errs, getErrorInfo(10));
450  }
451  m_sens_ok = false;
452 }
453 
454 double CVodesIntegrator::step(double tout)
455 {
456  int flag = CVode(m_cvode_mem, tout, m_y, &m_time, CV_ONE_STEP);
457  if (flag != CV_SUCCESS) {
458  string f_errs = m_func->getErrors();
459  if (!f_errs.empty()) {
460  f_errs = "Exceptions caught during RHS evaluation:\n" + f_errs;
461  }
462  throw CanteraError("CVodesIntegrator::step",
463  "CVodes error encountered. Error code: {}\n{}\n"
464  "{}"
465  "Components with largest weighted error estimates:\n{}",
466  flag, f_errs, m_error_message, getErrorInfo(10));
467 
468  }
469  m_sens_ok = false;
470  return m_time;
471 }
472 
474 {
475  long int ne;
476  CVodeGetNumRhsEvals(m_cvode_mem, &ne);
477  return ne;
478 }
479 
480 double CVodesIntegrator::sensitivity(size_t k, size_t p)
481 {
482  if (m_time == m_t0) {
483  // calls to CVodeGetSens are only allowed after a successful time step.
484  return 0.0;
485  }
486  if (!m_sens_ok && m_np) {
487  int flag = CVodeGetSens(m_cvode_mem, &m_time, m_yS);
488  if (flag != CV_SUCCESS) {
489  throw CanteraError("CVodesIntegrator::sensitivity",
490  "CVodeGetSens failed. Error code: {}", flag);
491  }
492  m_sens_ok = true;
493  }
494 
495  if (k >= m_neq) {
496  throw CanteraError("CVodesIntegrator::sensitivity",
497  "sensitivity: k out of range ({})", k);
498  }
499  if (p >= m_np) {
500  throw CanteraError("CVodesIntegrator::sensitivity",
501  "sensitivity: p out of range ({})", p);
502  }
503  return NV_Ith_S(m_yS[p],k);
504 }
505 
507 {
508  N_Vector errs = N_VNew_Serial(static_cast<sd_size_t>(m_neq));
509  N_Vector errw = N_VNew_Serial(static_cast<sd_size_t>(m_neq));
510  CVodeGetErrWeights(m_cvode_mem, errw);
511  CVodeGetEstLocalErrors(m_cvode_mem, errs);
512 
513  vector<tuple<double, double, size_t> > weightedErrors;
514  for (size_t i=0; i<m_neq; i++) {
515  double err = NV_Ith_S(errs, i) * NV_Ith_S(errw, i);
516  weightedErrors.emplace_back(-abs(err), err, i);
517  }
518  N_VDestroy(errs);
519  N_VDestroy(errw);
520 
521  N = std::min(N, static_cast<int>(m_neq));
522  sort(weightedErrors.begin(), weightedErrors.end());
523  fmt::memory_buffer s;
524  for (int i=0; i<N; i++) {
525  format_to(s, "{}: {}\n",
526  get<2>(weightedErrors[i]), get<1>(weightedErrors[i]));
527  }
528  return to_string(s);
529 }
530 
531 }
Backward Differentiation.
Definition: Integrator.h:33
virtual void setProblemType(int probtype)
Set the problem type.
void applyOptions()
Applies user-specified options to the underlying CVODES solver.
Wrapper class for &#39;cvodes&#39; integrator from LLNL.
virtual int nEvals() const
The number of function evaluations.
vector_fp m_paramScales
Scaling factors for each sensitivity parameter.
Definition: FuncEval.h:89
void clearErrors()
Clear any previously-stored suppressed errors.
Definition: FuncEval.h:79
std::string getErrors() const
Return a string containing the text of any suppressed errors.
Definition: FuncEval.cpp:45
STL namespace.
void * m_linsol
Sundials linear solver object.
virtual void setMaxStepSize(double hmax)
Set the maximum step size.
virtual void setMethod(MethodType t)
Set the solution method.
virtual void getState(double *y)
Fill in the vector y with the current state of the system.
Definition: FuncEval.h:53
virtual void setMinStepSize(double hmin)
Set the minimum step size.
virtual void setIterator(IterType t)
Set the linear iterator.
virtual double * solution()
The current value of the solution of the system of equations.
Functional Iteration.
Definition: Integrator.h:45
virtual size_t neq()=0
Number of equations.
virtual void integrate(double tout)
Integrate the system of equations.
void * m_linsol_matrix
matrix used by Sundials
virtual void setSensitivityTolerances(double reltol, double abstol)
Set the sensitivity error tolerances.
virtual std::string getErrorInfo(int N)
Returns a string listing the weighted error estimates associated with each solution component...
Base class for exceptions thrown by Cantera classes.
Definition: ctexceptions.h:65
vector_fp m_sens_params
Values for the problem parameters for which sensitivities are computed This is the array which is per...
Definition: FuncEval.h:81
virtual size_t nparams()
Number of sensitivity parameters.
Definition: FuncEval.h:61
static void cvodes_err(int error_code, const char *module, const char *function, char *msg, void *eh_data)
Function called by CVodes when an error is encountered instead of writing to stdout.
virtual void setTolerances(double reltol, size_t n, double *abstol)
Set error tolerances.
virtual void setMaxErrTestFails(int n)
Set the maximum permissible number of error test failures.
int eval_nothrow(double t, double *y, double *ydot)
Evaluate the right-hand side using return code to indicate status.
Definition: FuncEval.cpp:12
bool m_sens_ok
Indicates whether the sensitivities stored in m_yS have been updated for at the current integrator ti...
static int cvodes_rhs(realtype t, N_Vector y, N_Vector ydot, void *f_data)
Function called by cvodes to evaluate ydot given y.
IterType
Specifies the method used for iteration.
Definition: Integrator.h:41
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:157
Newton Iteration.
Definition: Integrator.h:43
virtual doublereal step(double tout)
Integrate the system of equations.
Contains declarations for string manipulation functions within Cantera.
std::string m_error_message
Error message information provide by CVodes.
Virtual base class for ODE right-hand-side function evaluators.
Definition: FuncEval.h:26
Namespace for the Cantera kernel.
Definition: AnyMap.cpp:8
MethodType
Specifies the method used to integrate the system of equations.
Definition: Integrator.h:32
double m_time
The current integrator time.
virtual void initialize(double t0, FuncEval &func)
Initialize the integrator for a new problem.