integrator_cvode.h

Go to the documentation of this file.

#pragma once
 
#include <vector>
 
#include <cvode/cvode.h>               // prototypes for CVODE fcts., consts.
#include <nvector/nvector_serial.h>    // access to serial N_Vector
#include <sunmatrix/sunmatrix_dense.h> // access to dense SUNMatrix
#include <sunlinsol/sunlinsol_dense.h> // access to dense SUNLinearSolver
 
 
class integrator_cvode {
 
public: 
 
 
    SUNContext      sun;          
    void           *cmem;         
    N_Vector        vars;         
    N_Vector        atol;         
    sunrealtype        rtol;         
    unsigned        nvar;         
    SUNMatrix       J;            
    SUNLinearSolver LS;           
    int             rv;           
 
 
 
integrator_cvode(
                 int (*Func)(sunrealtype, N_Vector, N_Vector, void*),
                 void *                _user_data,
                 const int             _nvar, 
                 const double          _rtol,
                 const std::vector<double> &_atol) :
    nvar(_nvar),
    rtol(_rtol) {
 
    rv = SUNContext_Create(SUN_COMM_NULL, &sun);
 
    vars = N_VNew_Serial(nvar, sun);
    atol = N_VNew_Serial(nvar, sun);
    for(int k=0; k<nvar; ++k)
        NV_Ith_S(atol, k) = _atol[k];
 
    cmem = CVodeCreate(CV_BDF, sun);
    rv   = CVodeSetUserData(cmem, _user_data);
    rv   = CVodeInit(cmem, Func, 0.0, vars);
    rv   = CVodeSVtolerances(cmem, rtol, atol);
    rv   = CVodeSetMaxNumSteps(cmem, 5000);
 
    J    = SUNDenseMatrix(nvar, nvar, sun);   // linear solver matrix J
    LS   = SUNLinSol_Dense(vars, J, sun);          // set linear solver
    rv   = CVodeSetLinearSolver(cmem, LS, J);  // associate matrix J and solver LS
}
 
//--------------
 
 
int integrate(std::vector<double> &y, const sunrealtype dt) {
 
    for(int k=0; k<nvar; ++k)
        NV_Ith_S(vars, k) = y[k];
 
    sunrealtype t;
 
    rv = CVodeReInit(cmem, 0.0, vars);
    rv = CVode(cmem, dt, vars, &t, CV_NORMAL);
 
    for(int k=0; k<nvar; ++k)
        y[k] = NV_Ith_S(vars,k);
 
    return rv;
}
 
 
~integrator_cvode() {
 
    N_VDestroy(vars);
    N_VDestroy(atol);
    CVodeFree(&cmem);
    SUNLinSolFree(LS);
    SUNMatDestroy(J);
    SUNContext_Free(&sun);
}
};