Cantera  2.4.0
IDA_Solver.cpp
Go to the documentation of this file.
1 //! @file IDA_Solver.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 "sundials/sundials_types.h"
10 #include "sundials/sundials_math.h"
11 #include "ida/ida.h"
12 #if CT_SUNDIALS_VERSION >= 30
13  #if CT_SUNDIALS_USE_LAPACK
14  #include "sunlinsol/sunlinsol_lapackdense.h"
15  #include "sunlinsol/sunlinsol_lapackband.h"
16  #else
17  #include "sunlinsol/sunlinsol_dense.h"
18  #include "sunlinsol/sunlinsol_band.h"
19  #endif
20  #include "sunlinsol/sunlinsol_spgmr.h"
21  #include "ida/ida_direct.h"
22  #include "ida/ida_spils.h"
23 #else
24  #include "ida/ida_dense.h"
25  #include "ida/ida_spgmr.h"
26  #include "ida/ida_band.h"
27 #endif
28 #include "nvector/nvector_serial.h"
29 
30 using namespace std;
31 
32 #if CT_SUNDIALS_VERSION < 25
33 typedef int sd_size_t;
34 #else
35 typedef long int sd_size_t;
36 #endif
37 
38 namespace Cantera
39 {
40 
41 //! A simple class to hold an array of parameter values and a pointer to an
42 //! instance of a subclass of ResidEval.
43 class ResidData
44 {
45 public:
46  ResidData(ResidJacEval* f, IDA_Solver* s, int npar = 0) {
47  m_func = f;
48  m_solver = s;
49  }
50 
51  virtual ~ResidData() {
52  }
53 
54  ResidJacEval* m_func;
55  IDA_Solver* m_solver;
56 };
57 }
58 
59 extern "C" {
60  //! Function called by IDA to evaluate the residual, given y and ydot.
61  /*!
62  * IDA allows passing in a void* pointer to access external data. Instead of
63  * requiring the user to provide a residual function directly to IDA (which
64  * would require using the sundials data types N_Vector, etc.), we define
65  * this function as the single function that IDA always calls. The real
66  * evaluation of the residual is done by an instance of a subclass of
67  * ResidEval, passed in to this function as a pointer in the parameters.
68  *
69  * FROM IDA WRITEUP -> What the IDA solver expects as a return flag from its
70  * residual routines:
71  *
72  * A IDAResFn res should return a value of 0 if successful, a positive value
73  * if a recoverable error occured (e.g. yy has an illegal value), or a
74  * negative value if a nonrecoverable error occured. In the latter case, the
75  * program halts. If a recoverable error occured, the integrator will
76  * attempt to correct and retry.
77  */
78  static int ida_resid(realtype t, N_Vector y, N_Vector ydot, N_Vector r, void* f_data)
79  {
81  Cantera::ResidJacEval* f = d->m_func;
82  Cantera::IDA_Solver* s = d->m_solver;
83  double delta_t = s->getCurrentStepFromIDA();
84  // TODO evaluate evalType. Assumed to be Base_ResidEval
85  int flag = f->evalResidNJ(t, delta_t, NV_DATA_S(y), NV_DATA_S(ydot),
86  NV_DATA_S(r));
87  if (flag < 0) {
88  // This signals to IDA that a nonrecoverable error has occurred.
89  return flag;
90  } else {
91  return 0;
92  }
93  }
94 
95  //! Function called by by IDA to evaluate the Jacobian, given y and ydot.
96  /*!
97  * typedef int (*IDADlsDenseJacFn)(sd_size_t N, realtype t, realtype c_j,
98  * N_Vector y, N_Vector yp, N_Vector r,
99  * DlsMat Jac, void *user_data,
100  * N_Vector tmp1, N_Vector tmp2, N_Vector tmp3);
101  *
102  * A IDADlsDenseJacFn should return
103  * - 0 if successful,
104  * - a positive int if a recoverable error occurred, or
105  * - a negative int if a nonrecoverable error occurred.
106  *
107  * In the case of a recoverable error return, the integrator will attempt to
108  * recover by reducing the stepsize (which changes cj).
109  */
110 #if CT_SUNDIALS_VERSION >= 30
111  static int ida_jacobian(realtype t, realtype c_j, N_Vector y, N_Vector yp,
112  N_Vector r, SUNMatrix Jac, void *f_data,
113  N_Vector tmp1, N_Vector tmp2, N_Vector tmp3)
114  {
115  Cantera::ResidData* d = (Cantera::ResidData*) f_data;
116  Cantera::ResidJacEval* f = d->m_func;
117  Cantera::IDA_Solver* s = d->m_solver;
118  double delta_t = s->getCurrentStepFromIDA();
119  double** cols;
120  if (SUNMatGetID(Jac) == SUNMATRIX_DENSE) {
121  cols = SM_COLS_D(Jac);
122  } else if (SUNMatGetID(Jac) == SUNMATRIX_BAND) {
123  cols = SM_COLS_B(Jac);
124  } else {
125  throw Cantera::CanteraError("ida_jacobian", "Unknown SUNMatrix type");
126  }
127  f->evalJacobianDP(t, delta_t, c_j, NV_DATA_S(y), NV_DATA_S(yp),
128  cols, NV_DATA_S(r));
129  return 0;
130  }
131 #else
132  static int ida_jacobian(sd_size_t nrows, realtype t, realtype c_j, N_Vector y, N_Vector ydot, N_Vector r,
133  DlsMat Jac, void* f_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3)
134  {
135  Cantera::ResidData* d = (Cantera::ResidData*) f_data;
136  Cantera::ResidJacEval* f = d->m_func;
137  Cantera::IDA_Solver* s = d->m_solver;
138  double delta_t = s->getCurrentStepFromIDA();
139  f->evalJacobianDP(t, delta_t, c_j, NV_DATA_S(y), NV_DATA_S(ydot),
140  Jac->cols, NV_DATA_S(r));
141  return 0;
142  }
143 #endif
144 
145 }
146 
147 namespace Cantera
148 {
149 
150 IDA_Solver::IDA_Solver(ResidJacEval& f) :
151  DAE_Solver(f),
152  m_ida_mem(0),
153  m_t0(0.0),
154  m_y(0),
155  m_ydot(0),
156  m_id(0),
157  m_constraints(0),
158  m_abstol(0),
159  m_type(0),
160  m_itol(IDA_SS),
161  m_iter(0),
162  m_reltol(1.e-9),
163  m_abstols(1.e-15),
164  m_nabs(0),
165  m_hmax(0.0),
166  m_hmin(0.0),
167  m_h0(0.0),
168  m_maxsteps(20000),
169  m_maxord(0),
170  m_formJac(0),
171  m_tstop(0.0),
172  m_told_old(0.0),
173  m_told(0.0),
174  m_tcurrent(0.0),
175  m_deltat(0.0),
176  m_maxErrTestFails(-1),
177  m_maxNonlinIters(0),
178  m_maxNonlinConvFails(-1),
179  m_setSuppressAlg(0),
180  m_mupper(0),
181  m_mlower(0)
182 {
183 }
184 
185 IDA_Solver::~IDA_Solver()
186 {
187  if (m_ida_mem) {
188  IDAFree(&m_ida_mem);
189  }
190  if (m_y) {
191  N_VDestroy_Serial(m_y);
192  }
193  if (m_ydot) {
194  N_VDestroy_Serial(m_ydot);
195  }
196  if (m_abstol) {
197  N_VDestroy_Serial(m_abstol);
198  }
199  if (m_constraints) {
200  N_VDestroy_Serial(m_constraints);
201  }
202 }
203 
204 doublereal IDA_Solver::solution(int k) const
205 {
206  return NV_Ith_S(m_y,k);
207 }
208 
209 const doublereal* IDA_Solver::solutionVector() const
210 {
211  return NV_DATA_S(m_y);
212 }
213 
214 doublereal IDA_Solver::derivative(int k) const
215 {
216  return NV_Ith_S(m_ydot,k);
217 }
218 
219 const doublereal* IDA_Solver::derivativeVector() const
220 {
221  return NV_DATA_S(m_ydot);
222 }
223 
224 void IDA_Solver::setTolerances(double reltol, double* abstol)
225 {
226  m_itol = IDA_SV;
227  if (!m_abstol) {
228  m_abstol = N_VNew_Serial(m_neq);
229  }
230  for (int i = 0; i < m_neq; i++) {
231  NV_Ith_S(m_abstol, i) = abstol[i];
232  }
233  m_reltol = reltol;
234  if (m_ida_mem) {
235  int flag = IDASVtolerances(m_ida_mem, m_reltol, m_abstol);
236  if (flag != IDA_SUCCESS) {
237  throw CanteraError("IDA_Solver::setTolerances",
238  "Memory allocation failed.");
239  }
240  }
241 }
242 
243 void IDA_Solver::setTolerances(doublereal reltol, doublereal abstol)
244 {
245  m_itol = IDA_SS;
246  m_reltol = reltol;
247  m_abstols = abstol;
248  if (m_ida_mem) {
249  int flag = IDASStolerances(m_ida_mem, m_reltol, m_abstols);
250  if (flag != IDA_SUCCESS) {
251  throw CanteraError("IDA_Solver::setTolerances",
252  "Memory allocation failed.");
253  }
254  }
255 }
256 
257 void IDA_Solver::setLinearSolverType(int solverType)
258 {
259  m_type = solverType;
260 }
261 
263 {
264  setLinearSolverType(0);
265 }
266 
267 void IDA_Solver::setBandedLinearSolver(int m_upper, int m_lower)
268 {
269  m_type = 2;
270  m_upper = m_mupper;
271  m_mlower = m_lower;
272 }
273 
274 void IDA_Solver::setMaxOrder(int n)
275 {
276  m_maxord = n;
277 }
278 
280 {
281  m_maxsteps = n;
282 }
283 
285 {
286  m_h0 = h0;
287 }
288 
289 void IDA_Solver::setStopTime(doublereal tstop)
290 {
291  m_tstop = tstop;
292 }
293 
295 {
296  doublereal hcur;
297  IDAGetCurrentStep(m_ida_mem, &hcur);
298  return hcur;
299 }
300 
302 {
303  m_formJac = formJac;
304  if (m_ida_mem && m_formJac == 1) {
305  #if CT_SUNDIALS_VERSION >= 30
306  int flag = IDADlsSetJacFn(m_ida_mem, ida_jacobian);
307  #else
308  int flag = IDADlsSetDenseJacFn(m_ida_mem, ida_jacobian);
309  #endif
310  if (flag != IDA_SUCCESS) {
311  throw CanteraError("IDA_Solver::setJacobianType",
312  "IDADlsSetDenseJacFn failed.");
313  }
314  }
315 }
316 
317 void IDA_Solver::setMaxErrTestFailures(int maxErrTestFails)
318 {
319  m_maxErrTestFails = maxErrTestFails;
320 }
321 
323 {
324  m_maxNonlinIters = n;
325 }
326 
328 {
330 }
331 
332 void IDA_Solver::inclAlgebraicInErrorTest(bool yesno)
333 {
334  if (yesno) {
335  m_setSuppressAlg = 0;
336  } else {
337  m_setSuppressAlg = 1;
338  }
339 }
340 
341 void IDA_Solver::init(doublereal t0)
342 {
343  m_t0 = t0;
344  m_told = t0;
345  m_told_old = t0;
346  m_tcurrent = t0;
347  if (m_y) {
348  N_VDestroy_Serial(m_y);
349  }
350  if (m_ydot) {
351  N_VDestroy_Serial(m_ydot);
352  }
353  if (m_id) {
354  N_VDestroy_Serial(m_id);
355  }
356  if (m_constraints) {
357  N_VDestroy_Serial(m_constraints);
358  }
359 
360  m_y = N_VNew_Serial(m_neq);
361  m_ydot = N_VNew_Serial(m_neq);
362  m_constraints = N_VNew_Serial(m_neq);
363 
364  for (int i=0; i<m_neq; i++) {
365  NV_Ith_S(m_y, i) = 0.0;
366  NV_Ith_S(m_ydot, i) = 0.0;
367  NV_Ith_S(m_constraints, i) = 0.0;
368  }
369 
370  // get the initial conditions
371  m_resid.getInitialConditions(m_t0, NV_DATA_S(m_y), NV_DATA_S(m_ydot));
372 
373  if (m_ida_mem) {
374  IDAFree(&m_ida_mem);
375  }
376 
377  /* Call IDACreate */
378  m_ida_mem = IDACreate();
379 
380  int flag = 0;
381  if (m_itol == IDA_SV) {
382  flag = IDAInit(m_ida_mem, ida_resid, m_t0, m_y, m_ydot);
383  if (flag != IDA_SUCCESS) {
384  if (flag == IDA_MEM_FAIL) {
385  throw CanteraError("IDA_Solver::init",
386  "Memory allocation failed.");
387  } else if (flag == IDA_ILL_INPUT) {
388  throw CanteraError("IDA_Solver::init",
389  "Illegal value for IDAInit input argument.");
390  } else {
391  throw CanteraError("IDA_Solver::init", "IDAInit failed.");
392  }
393  }
394  flag = IDASVtolerances(m_ida_mem, m_reltol, m_abstol);
395  if (flag != IDA_SUCCESS) {
396  throw CanteraError("IDA_Solver::init", "Memory allocation failed.");
397  }
398  } else {
399  flag = IDAInit(m_ida_mem, ida_resid, m_t0, m_y, m_ydot);
400  if (flag != IDA_SUCCESS) {
401  if (flag == IDA_MEM_FAIL) {
402  throw CanteraError("IDA_Solver::init",
403  "Memory allocation failed.");
404  } else if (flag == IDA_ILL_INPUT) {
405  throw CanteraError("IDA_Solver::init",
406  "Illegal value for IDAInit input argument.");
407  } else {
408  throw CanteraError("IDA_Solver::init", "IDAInit failed.");
409  }
410  }
411  flag = IDASStolerances(m_ida_mem, m_reltol, m_abstols);
412  if (flag != IDA_SUCCESS) {
413  throw CanteraError("IDA_Solver::init", "Memory allocation failed.");
414  }
415  }
416 
417  // set the linear solver type
418  if (m_type == 1 || m_type == 0) {
419  long int N = m_neq;
420  int flag;
421  #if CT_SUNDIALS_VERSION >= 30
422  SUNLinSolFree((SUNLinearSolver) m_linsol);
423  SUNMatDestroy((SUNMatrix) m_linsol_matrix);
424  m_linsol_matrix = SUNDenseMatrix(N, N);
425  #if CT_SUNDIALS_USE_LAPACK
426  m_linsol = SUNLapackDense(m_y, (SUNMatrix) m_linsol_matrix);
427  #else
428  m_linsol = SUNDenseLinearSolver(m_y, (SUNMatrix) m_linsol_matrix);
429  #endif
430  flag = IDADlsSetLinearSolver(m_ida_mem, (SUNLinearSolver) m_linsol,
431  (SUNMatrix) m_linsol_matrix);
432  #else
433  flag = IDADense(m_ida_mem, N);
434  #endif
435  if (flag) {
436  throw CanteraError("IDA_Solver::init", "IDADense failed");
437  }
438  } else if (m_type == 2) {
439  long int N = m_neq;
440  long int nu = m_mupper;
441  long int nl = m_mlower;
442  #if CT_SUNDIALS_VERSION >= 30
443  SUNLinSolFree((SUNLinearSolver) m_linsol);
444  SUNMatDestroy((SUNMatrix) m_linsol_matrix);
445  m_linsol_matrix = SUNBandMatrix(N, nu, nl, nu+nl);
446  #if CT_SUNDIALS_USE_LAPACK
447  m_linsol = SUNLapackBand(m_y, (SUNMatrix) m_linsol_matrix);
448  #else
449  m_linsol = SUNBandLinearSolver(m_y, (SUNMatrix) m_linsol_matrix);
450  #endif
451  IDADlsSetLinearSolver(m_ida_mem, (SUNLinearSolver) m_linsol,
452  (SUNMatrix) m_linsol_matrix);
453  #else
454  IDABand(m_ida_mem, N, nu, nl);
455  #endif
456  } else {
457  throw CanteraError("IDA_Solver::init",
458  "unsupported linear solver type");
459  }
460 
461  if (m_formJac == 1) {
462  #if CT_SUNDIALS_VERSION >= 30
463  flag = IDADlsSetJacFn(m_ida_mem, ida_jacobian);
464  #else
465  flag = IDADlsSetDenseJacFn(m_ida_mem, ida_jacobian);
466  #endif
467  if (flag != IDA_SUCCESS) {
468  throw CanteraError("IDA_Solver::init",
469  "IDADlsSetDenseJacFn failed.");
470  }
471  }
472 
473  // pass a pointer to func in m_data
474  m_fdata.reset(new ResidData(&m_resid, this, m_resid.nparams()));
475  flag = IDASetUserData(m_ida_mem, m_fdata.get());
476  if (flag != IDA_SUCCESS) {
477  throw CanteraError("IDA_Solver::init", "IDASetUserData failed.");
478  }
479 
480  // set options
481  if (m_maxord > 0) {
482  flag = IDASetMaxOrd(m_ida_mem, m_maxord);
483  if (flag != IDA_SUCCESS) {
484  throw CanteraError("IDA_Solver::init", "IDASetMaxOrd failed.");
485  }
486  }
487  if (m_maxsteps > 0) {
488  flag = IDASetMaxNumSteps(m_ida_mem, m_maxsteps);
489  if (flag != IDA_SUCCESS) {
490  throw CanteraError("IDA_Solver::init", "IDASetMaxNumSteps failed.");
491  }
492  }
493  if (m_h0 > 0.0) {
494  flag = IDASetInitStep(m_ida_mem, m_h0);
495  if (flag != IDA_SUCCESS) {
496  throw CanteraError("IDA_Solver::init", "IDASetInitStep failed.");
497  }
498  }
499  if (m_tstop > 0.0) {
500  flag = IDASetStopTime(m_ida_mem, m_tstop);
501  if (flag != IDA_SUCCESS) {
502  throw CanteraError("IDA_Solver::init", "IDASetStopTime failed.");
503  }
504  }
505  if (m_maxErrTestFails >= 0) {
506  flag = IDASetMaxErrTestFails(m_ida_mem, m_maxErrTestFails);
507  if (flag != IDA_SUCCESS) {
508  throw CanteraError("IDA_Solver::init",
509  "IDASetMaxErrTestFails failed.");
510  }
511  }
512  if (m_maxNonlinIters > 0) {
513  flag = IDASetMaxNonlinIters(m_ida_mem, m_maxNonlinIters);
514  if (flag != IDA_SUCCESS) {
515  throw CanteraError("IDA_Solver::init",
516  "IDASetmaxNonlinIters failed.");
517  }
518  }
519  if (m_maxNonlinConvFails >= 0) {
520  flag = IDASetMaxConvFails(m_ida_mem, m_maxNonlinConvFails);
521  if (flag != IDA_SUCCESS) {
522  throw CanteraError("IDA_Solver::init",
523  "IDASetMaxConvFails failed.");
524  }
525  }
526  if (m_setSuppressAlg != 0) {
527  flag = IDASetSuppressAlg(m_ida_mem, m_setSuppressAlg);
528  if (flag != IDA_SUCCESS) {
529  throw CanteraError("IDA_Solver::init", "IDASetSuppressAlg failed.");
530  }
531  }
532 }
533 
534 void IDA_Solver::correctInitial_Y_given_Yp(doublereal* y, doublereal* yp, doublereal tout)
535 {
536  doublereal tout1 = tout;
537  if (tout == 0.0) {
538  double h0 = 1.0E-5;
539  if (m_h0 > 0.0) {
540  h0 = m_h0;
541  }
542  tout1 = m_t0 + h0;
543  }
544 
545  int flag = IDACalcIC(m_ida_mem, IDA_Y_INIT, tout1);
546  if (flag != IDA_SUCCESS) {
547  throw CanteraError("IDA_Solver::correctInitial_Y_given_Yp",
548  "IDACalcIC failed: error = {}", flag);
549  }
550 
551  flag = IDAGetConsistentIC(m_ida_mem, m_y, m_ydot);
552  if (flag != IDA_SUCCESS) {
553  throw CanteraError("IDA_Solver::correctInitial_Y_given_Yp",
554  "IDAGetSolution failed: error = {}", flag);
555  }
556  for (int i = 0; i < m_neq; i++) {
557  y[i] = NV_Ith_S(m_y, i);
558  yp[i] = NV_Ith_S(m_ydot, i);
559  }
560 }
561 
562 void IDA_Solver::correctInitial_YaYp_given_Yd(doublereal* y, doublereal* yp, doublereal tout)
563 {
564  int icopt = IDA_YA_YDP_INIT;
565  doublereal tout1 = tout;
566  if (tout == 0.0) {
567  double h0 = 1.0E-5;
568  if (m_h0 > 0.0) {
569  h0 = m_h0;
570  }
571  tout1 = m_t0 + h0;
572  }
573 
574  int flag = IDACalcIC(m_ida_mem, icopt, tout1);
575  if (flag != IDA_SUCCESS) {
576  throw CanteraError("IDA_Solver::correctInitial_YaYp_given_Yd",
577  "IDACalcIC failed: error = {}", flag);
578  }
579 
580  flag = IDAGetConsistentIC(m_ida_mem, m_y, m_ydot);
581  if (flag != IDA_SUCCESS) {
582  throw CanteraError("IDA_Solver::correctInitial_YaYp_given_Yd",
583  "IDAGetSolution failed: error = {}", flag);
584  }
585  for (int i = 0; i < m_neq; i++) {
586  y[i] = NV_Ith_S(m_y, i);
587  yp[i] = NV_Ith_S(m_ydot, i);
588  }
589 }
590 
591 int IDA_Solver::solve(double tout)
592 {
593  double tretn = tout - 1000;
594  int flag;
595  flag = IDASetStopTime(m_ida_mem, tout);
596  if (flag != IDA_SUCCESS) {
597  throw CanteraError("IDA_Solver::solve", "IDA error encountered.");
598  }
599  while (tretn < tout) {
600  if (tout <= m_tcurrent) {
601  throw CanteraError("IDA_Solver::solve", "tout <= tcurrent");
602  }
603  m_told_old = m_told;
604  m_told = m_tcurrent;
605  flag = IDASolve(m_ida_mem, tout, &tretn, m_y, m_ydot, IDA_ONE_STEP);
606  if (flag < 0) {
607  throw CanteraError("IDA_Solver::solve", "IDA error encountered.");
608  } else if (flag == IDA_TSTOP_RETURN) {
609  // we've reached our goal, and have actually integrated past it
610  } else if (flag == IDA_ROOT_RETURN) {
611  // not sure what to do with this yet
612  } else if (flag == IDA_WARNING) {
613  throw CanteraError("IDA_Solver::solve", "IDA Warning encountered.");
614  }
615  m_tcurrent = tretn;
617  };
618 
619  if (flag != IDA_SUCCESS && flag != IDA_TSTOP_RETURN) {
620  throw CanteraError("IDA_Solver::solve", "IDA error encountered.");
621  }
622  return flag;
623 }
624 
625 double IDA_Solver::step(double tout)
626 {
627  double t;
628  int flag;
629  if (tout <= m_tcurrent) {
630  throw CanteraError("IDA_Solver::step", "tout <= tcurrent");
631  }
632  m_told_old = m_told;
633  m_told = m_tcurrent;
634  flag = IDASolve(m_ida_mem, tout, &t, m_y, m_ydot, IDA_ONE_STEP);
635  if (flag < 0) {
636  throw CanteraError("IDA_Solver::step", "IDA error encountered.");
637  } else if (flag == IDA_TSTOP_RETURN) {
638  // we've reached our goal, and have actually integrated past it
639  } else if (flag == IDA_ROOT_RETURN) {
640  // not sure what to do with this yet
641  } else if (flag == IDA_WARNING) {
642  throw CanteraError("IDA_Solver::step", "IDA Warning encountered.");
643  }
644  m_tcurrent = t;
646  return t;
647 }
648 
649 doublereal IDA_Solver::getOutputParameter(int flag) const
650 {
651  long int lenrw, leniw;
652  switch (flag) {
653  case REAL_WORKSPACE_SIZE:
654  flag = IDAGetWorkSpace(m_ida_mem, &lenrw, &leniw);
655  return doublereal(lenrw);
656  break;
657  }
658  return 0.0;
659 }
660 
661 }
void * m_linsol
Sundials linear solver object.
Definition: IDA_Solver.h:226
virtual void setStopTime(doublereal tstop)
Set the stop time.
Definition: IDA_Solver.cpp:289
virtual int evalJacobianDP(const doublereal t, const doublereal delta_t, doublereal cj, const doublereal *const y, const doublereal *const ydot, doublereal *const *jacobianColPts, doublereal *const resid)
Calculate an analytical Jacobian and the residual at the current time and values. ...
A simple class to hold an array of parameter values and a pointer to an instance of a subclass of Res...
Definition: IDA_Solver.cpp:43
integer m_neq
Number of total equations in the system.
Definition: DAE_Solver.h:249
Wrapper for Sundials IDA solver.
Definition: IDA_Solver.h:40
int m_maxNonlinConvFails
Maximum number of nonlinear convergence failures.
Definition: IDA_Solver.h:297
int nparams() const
Return the number of parameters in the calculation.
Definition: ResidEval.h:157
virtual void setBandedLinearSolver(int m_upper, int m_lower)
Set up the problem to use a band solver.
Definition: IDA_Solver.cpp:267
doublereal m_told
Value of the previous time.
Definition: IDA_Solver.h:279
int m_maxsteps
Maximum number of time steps allowed.
Definition: IDA_Solver.h:258
virtual double getCurrentStepFromIDA()
Get the current step size from IDA via a call.
Definition: IDA_Solver.cpp:294
virtual void setInitialStepSize(doublereal h0)
Set the initial step size.
Definition: IDA_Solver.cpp:284
virtual void setDenseLinearSolver()
Set up the problem to use a dense linear direct solver.
Definition: IDA_Solver.cpp:262
virtual int evalResidNJ(const doublereal t, const doublereal delta_t, const doublereal *const y, const doublereal *const ydot, doublereal *const resid, const ResidEval_Type_Enum evalType=Base_ResidEval, const int id_x=-1, const doublereal delta_x=0.0)
Evaluate the residual function.
int m_maxErrTestFails
maximum number of error test failures
Definition: IDA_Solver.h:288
virtual void correctInitial_YaYp_given_Yd(doublereal *y, doublereal *yp, doublereal tout)
Calculate consistent value of the algebraic constraints and derivatives at the start of the problem...
Definition: IDA_Solver.cpp:562
doublereal m_tstop
maximum time
Definition: IDA_Solver.h:273
STL namespace.
doublereal m_tcurrent
Value of the current time.
Definition: IDA_Solver.h:282
virtual int getInitialConditions(const doublereal t0, doublereal *const y, doublereal *const ydot)
Fill in the initial conditions.
virtual void setMaxNonlinConvFailures(int n)
Set the maximum number of nonlinear solver convergence failures.
Definition: IDA_Solver.cpp:327
doublereal m_told_old
Value of the previous, previous time.
Definition: IDA_Solver.h:276
static int ida_resid(realtype t, N_Vector y, N_Vector ydot, N_Vector r, void *f_data)
Function called by IDA to evaluate the residual, given y and ydot.
Definition: IDA_Solver.cpp:78
Wrappers for the function evaluators for Nonlinear solvers and Time steppers.
Definition: ResidJacEval.h:55
void * m_ida_mem
Pointer to the IDA memory for the problem.
Definition: IDA_Solver.h:225
virtual void setMaxNumSteps(int n)
Set the maximum number of time steps.
Definition: IDA_Solver.cpp:279
void * m_linsol_matrix
matrix used by Sundials
Definition: IDA_Solver.h:227
N_Vector m_y
Current value of the solution vector.
Definition: IDA_Solver.h:233
doublereal m_h0
Value of the initial time step.
Definition: IDA_Solver.h:255
virtual int solve(doublereal tout)
Step the system to a final value of the time.
Definition: IDA_Solver.cpp:591
virtual void correctInitial_Y_given_Yp(doublereal *y, doublereal *yp, doublereal tout)
Calculate consistent value of the starting solution given the starting solution derivatives.
Definition: IDA_Solver.cpp:534
doublereal m_deltat
Value of deltaT for the current step.
Definition: IDA_Solver.h:285
virtual doublereal solution(int k) const
the current value of solution component k.
Definition: IDA_Solver.cpp:204
Base class for exceptions thrown by Cantera classes.
Definition: ctexceptions.h:65
virtual doublereal derivative(int k) const
the current value of the derivative of solution component k.
Definition: IDA_Solver.cpp:214
virtual void setMaxNonlinIterations(int n)
Set the maximum number of nonlinear iterations on a timestep.
Definition: IDA_Solver.cpp:322
virtual void setJacobianType(int formJac)
Set the form of the Jacobian.
Definition: IDA_Solver.cpp:301
virtual doublereal getOutputParameter(int flag) const
Get the value of a solver-specific output parameter.
Definition: IDA_Solver.cpp:649
virtual void init(doublereal t0)
initialize.
Definition: IDA_Solver.cpp:341
Contains declarations for string manipulation functions within Cantera.
int m_maxNonlinIters
Maximum number of nonlinear solver iterations at one solution.
Definition: IDA_Solver.h:294
doublereal m_t0
Initial value of the time.
Definition: IDA_Solver.h:230
N_Vector m_ydot
Current value of the derivative of the solution vector.
Definition: IDA_Solver.h:236
Header file for class IDA_Solver.
virtual doublereal step(doublereal tout)
Take one internal step.
Definition: IDA_Solver.cpp:625
Namespace for the Cantera kernel.
Definition: AnyMap.cpp:8
Wrapper for DAE solvers.
Definition: DAE_Solver.h:75
int m_maxord
maximum time step order of the method
Definition: IDA_Solver.h:261
virtual void setTolerances(doublereal reltol, doublereal *abstol)
Set error tolerances.
Definition: IDA_Solver.cpp:224
static int ida_jacobian(sd_size_t nrows, realtype t, realtype c_j, N_Vector y, N_Vector ydot, N_Vector r, DlsMat Jac, void *f_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3)
Function called by by IDA to evaluate the Jacobian, given y and ydot.
Definition: IDA_Solver.cpp:132
int m_formJac
Form of the Jacobian.
Definition: IDA_Solver.h:270
int m_setSuppressAlg
If true, the algebraic variables don&#39;t contribute to error tolerances.
Definition: IDA_Solver.h:300