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" 22 #include "sunlinsol/sunlinsol_dense.h" 23 #include "sunlinsol/sunlinsol_band.h" 25 #include "sunlinsol/sunlinsol_spgmr.h" 26 #include "cvodes/cvodes_direct.h" 27 #include "cvodes/cvodes_diag.h" 28 #include "cvodes/cvodes_spils.h" 30 #if CT_SUNDIALS_USE_LAPACK 31 #include "cvodes/cvodes_lapack.h" 33 #include "cvodes/cvodes_dense.h" 34 #include "cvodes/cvodes_band.h" 36 #include "cvodes/cvodes_diag.h" 37 #include "cvodes/cvodes_spgmr.h" 43 #if CT_SUNDIALS_VERSION < 25 44 typedef int sd_size_t;
46 typedef long int sd_size_t;
61 static int cvodes_rhs(realtype t, N_Vector y, N_Vector ydot,
void* f_data)
64 return f->
eval_nothrow(t, NV_DATA_S(y), NV_DATA_S(ydot));
70 static void cvodes_err(
int error_code,
const char* module,
71 const char*
function,
char* msg,
void* eh_data)
79 CVodesIntegrator::CVodesIntegrator() :
101 m_maxErrTestFails(0),
104 m_mupper(0), m_mlower(0),
109 CVodesIntegrator::~CVodesIntegrator()
113 CVodeSensFree(m_cvode_mem);
115 CVodeFree(&m_cvode_mem);
118 #if CT_SUNDIALS_VERSION >= 30 119 SUNLinSolFree((SUNLinearSolver)
m_linsol);
124 N_VDestroy_Serial(m_y);
127 N_VDestroy_Serial(m_abstol);
130 N_VDestroyVectorArray_Serial(m_yS, static_cast<sd_size_t>(m_np));
136 return NV_Ith_S(m_y, k);
141 return NV_DATA_S(m_y);
150 N_VDestroy_Serial(m_abstol);
152 m_abstol = N_VNew_Serial(static_cast<sd_size_t>(n));
154 for (
size_t i=0; i<n; i++) {
155 NV_Ith_S(m_abstol, i) = abstol[i];
169 m_reltolsens = reltol;
170 m_abstolsens = abstol;
185 throw CanteraError(
"CVodesIntegrator::setMethod",
"unknown method");
193 CVodeSetMaxStep(m_cvode_mem, hmax);
201 CVodeSetMinStep(m_cvode_mem, hmin);
205 void CVodesIntegrator::setMaxSteps(
int nmax)
209 CVodeSetMaxNumSteps(m_cvode_mem, m_maxsteps);
215 m_maxErrTestFails = n;
217 CVodeSetMaxErrTestFails(m_cvode_mem, n);
226 m_iter = CV_FUNCTIONAL;
228 throw CanteraError(
"CVodesIntegrator::setIterator",
"unknown iterator");
232 void CVodesIntegrator::sensInit(
double t0,
FuncEval& func)
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]);
242 N_VDestroy_Serial(y);
244 int flag = CVodeSensInit(m_cvode_mem, static_cast<sd_size_t>(m_np),
245 CV_STAGGERED, CVSensRhsFn(0), m_yS);
247 if (flag != CV_SUCCESS) {
248 throw CanteraError(
"CVodesIntegrator::sensInit",
"Error in CVodeSensInit");
251 for (
size_t n = 0; n < m_np; n++) {
256 flag = CVodeSensSStolerances(m_cvode_mem, m_reltolsens, atol.data());
268 N_VDestroy_Serial(m_y);
270 m_y = N_VNew_Serial(static_cast<sd_size_t>(m_neq));
273 if (m_itol == CV_SV && m_nabs < m_neq) {
275 "not enough absolute tolerance values specified.");
281 CVodeFree(&m_cvode_mem);
287 m_cvode_mem = CVodeCreate(m_method, m_iter);
290 "CVodeCreate failed.");
293 int flag = CVodeInit(m_cvode_mem,
cvodes_rhs, m_t0, m_y);
294 if (flag != CV_SUCCESS) {
295 if (flag == CV_MEM_FAIL) {
297 "Memory allocation failed.");
298 }
else if (flag == CV_ILL_INPUT) {
300 "Illegal value for CVodeInit input argument.");
303 "CVodeInit failed.");
306 CVodeSetErrHandlerFn(m_cvode_mem, &
cvodes_err,
this);
308 if (m_itol == CV_SV) {
309 flag = CVodeSVtolerances(m_cvode_mem, m_reltol, m_abstol);
311 flag = CVodeSStolerances(m_cvode_mem, m_reltol, m_abstols);
313 if (flag != CV_SUCCESS) {
314 if (flag == CV_MEM_FAIL) {
316 "Memory allocation failed.");
317 }
else if (flag == CV_ILL_INPUT) {
319 "Illegal value for CVodeInit input argument.");
322 "CVodeInit failed.");
326 flag = CVodeSetUserData(m_cvode_mem, &func);
327 if (flag != CV_SUCCESS) {
329 "CVodeSetUserData failed.");
333 flag = CVodeSetSensParams(m_cvode_mem, func.
m_sens_params.data(),
335 if (flag != CV_SUCCESS) {
337 "CVodeSetSensParams failed.");
343 void CVodesIntegrator::reinitialize(
double t0,
FuncEval& func)
351 int result = CVodeReInit(m_cvode_mem, m_t0, m_y);
352 if (result != CV_SUCCESS) {
354 "CVodeReInit failed. result = {}", result);
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);
367 #if CT_SUNDIALS_USE_LAPACK 372 CVDlsSetLinearSolver(m_cvode_mem, (SUNLinearSolver)
m_linsol,
375 #if CT_SUNDIALS_USE_LAPACK 376 CVLapackDense(m_cvode_mem, N);
378 CVDense(m_cvode_mem, N);
381 }
else if (m_type == DIAG) {
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);
388 CVSpgmr(m_cvode_mem, PREC_NONE, 0);
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);
398 #if CT_SUNDIALS_USE_LAPACK 403 CVDlsSetLinearSolver(m_cvode_mem, (SUNLinearSolver)
m_linsol,
406 #if CT_SUNDIALS_USE_LAPACK 407 CVLapackBand(m_cvode_mem, N, nu, nl);
409 CVBand(m_cvode_mem, N, nu, nl);
414 "unsupported option");
418 CVodeSetMaxOrd(m_cvode_mem, m_maxord);
420 if (m_maxsteps > 0) {
421 CVodeSetMaxNumSteps(m_cvode_mem, m_maxsteps);
424 CVodeSetMaxStep(m_cvode_mem, m_hmax);
427 CVodeSetMinStep(m_cvode_mem, m_hmin);
429 if (m_maxErrTestFails > 0) {
430 CVodeSetMaxErrTestFails(m_cvode_mem, m_maxErrTestFails);
439 int flag = CVode(m_cvode_mem, tout, m_y, &
m_time, CV_NORMAL);
440 if (flag != CV_SUCCESS) {
442 if (!f_errs.empty()) {
443 f_errs =
"Exceptions caught during RHS evaluation:\n" + f_errs;
446 "CVodes error encountered. Error code: {}\n{}\n" 448 "Components with largest weighted error estimates:\n{}",
456 int flag = CVode(m_cvode_mem, tout, m_y, &
m_time, CV_ONE_STEP);
457 if (flag != CV_SUCCESS) {
459 if (!f_errs.empty()) {
460 f_errs =
"Exceptions caught during RHS evaluation:\n" + f_errs;
463 "CVodes error encountered. Error code: {}\n{}\n" 465 "Components with largest weighted error estimates:\n{}",
476 CVodeGetNumRhsEvals(m_cvode_mem, &ne);
480 double CVodesIntegrator::sensitivity(
size_t k,
size_t p)
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);
496 throw CanteraError(
"CVodesIntegrator::sensitivity",
497 "sensitivity: k out of range ({})", k);
500 throw CanteraError(
"CVodesIntegrator::sensitivity",
501 "sensitivity: p out of range ({})", p);
503 return NV_Ith_S(m_yS[p],k);
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);
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);
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]));
Backward Differentiation.
virtual void setProblemType(int probtype)
Set the problem type.
void applyOptions()
Applies user-specified options to the underlying CVODES solver.
Wrapper class for 'cvodes' integrator from LLNL.
virtual int nEvals() const
The number of function evaluations.
vector_fp m_paramScales
Scaling factors for each sensitivity parameter.
void clearErrors()
Clear any previously-stored suppressed errors.
std::string getErrors() const
Return a string containing the text of any suppressed errors.
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.
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.
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.
vector_fp m_sens_params
Values for the problem parameters for which sensitivities are computed This is the array which is per...
virtual size_t nparams()
Number of sensitivity parameters.
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.
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.
std::vector< double > vector_fp
Turn on the use of stl vectors for the basic array type within cantera Vector of doubles.
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.
Namespace for the Cantera kernel.
MethodType
Specifies the method used to integrate the system of equations.
double m_time
The current integrator time.
virtual void initialize(double t0, FuncEval &func)
Initialize the integrator for a new problem.