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
#include <sunmatrix/sunmatrix_band.h>  // access to dense SUNMatrix
#include <sunlinsol/sunlinsol_band.h>  // access to dense SUNLinearSolver
 
 
class integrator_cvode {
 
public: 
 
 
    SUNContext      sun;           
    void           *cmem;          
    N_Vector        vars;          
    N_Vector        atol;          
    realtype        rtol;          
    unsigned        nvar;          
    SUNMatrix       J;             
    SUNLinearSolver LS;            
    int             rv;            
 
 
 
integrator_cvode(
                 int (*Func)(realtype, N_Vector, N_Vector, void*),
                 void *                _user_data,
                 const int             _nvar, 
                 const double          _rtol,
                 const std::vector<double> &_atol,
                 const int             mu,
                 const int             ml,
                 std::vector<double>  &y) :
                     nvar(_nvar),
                     rtol(_rtol) {
 
    rv = SUNContext_Create(NULL, &sun);
 
    vars = N_VNew_Serial(nvar, sun);
    atol = N_VNew_Serial(nvar, sun);
    for(int k=0; k<nvar; ++k) {
        NV_Ith_S(vars, k) = y[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, 20000); //5000);
 
    J    = SUNBandMatrix(nvar, mu, ml, sun);   // linear solver matrix J
    LS   = SUNLinSol_Band(vars, J, sun);          // set linear solver
    rv   = CVodeSetLinearSolver(cmem, LS, J);  // associate matrix J and solver LS
}
 
 
int integrate(std::vector<double> &y, const realtype dt) {
 
    for(int k=0; k<nvar; ++k)
        NV_Ith_S(vars, k) = y[k];
 
    realtype 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);
}
};