research/cantera
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

Sundials 2.4 compatibility.

Browse source 
Added compatibility for sundials 2.4 along the 1.8.0 branch. I'll
add it to the main trunk as well shortly. This hasn't been fully tested
except for the cxx_examples problem.
This commit is contained in:
Harry Moffat 2009-09-22 16:57:14 +00:00
parent f28e01aa23
commit 55bd4a766a
5 changed files with 4924 additions and 5451 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

768
Cantera/src/numerics/CVodesIntegrator.cpp
View file

@ -25,7 +25,7 @@ using namespace std;
#else
#ifdef SUNDIALS_VERSION_23
#if defined(SUNDIALS_VERSION_23) || defined (SUNDIALS_VERSION_24)
#include <sundials/sundials_types.h>
#include <sundials/sundials_math.h>
#include <sundials/sundials_nvector.h>
@ -42,383 +42,471 @@ unsupported sundials version!
#endif
#if defined (SUNDIALS_VERSION_24)
#define CV_SS 1
#define CV_SV 2
#endif
#endif
inline static N_Vector nv(void* x) {
return reinterpret_cast<N_Vector>(x);
return reinterpret_cast<N_Vector>(x);
}
namespace Cantera {
class FuncData {
public:
FuncData(FuncEval* f, int npar = 0) {
m_pars.resize(npar, 1.0);
m_func = f;
}
virtual ~FuncData() {}
vector_fp m_pars;
FuncEval* m_func;
};
class FuncData {
public:
FuncData(FuncEval* f, int npar = 0) {
m_pars.resize(npar, 1.0);
m_func = f;
}
virtual ~FuncData() {}
vector_fp m_pars;
FuncEval* m_func;
};
}
extern "C" {
/**
* Function called by cvodes to evaluate ydot given y. The cvode
* integrator allows passing in a void* pointer to access
* external data. This pointer is cast to a pointer to a instance
* of class FuncEval. The equations to be integrated should be
* specified by deriving a class from FuncEval that evaluates the
* desired equations.
* @ingroup odeGroup
*/
static int cvodes_rhs(realtype t, N_Vector y, N_Vector ydot,
void *f_data) {
double* ydata = NV_DATA_S(y); //N_VDATA(y);
double* ydotdata = NV_DATA_S(ydot); //N_VDATA(ydot);
Cantera::FuncData* d = (Cantera::FuncData*)f_data;
Cantera::FuncEval* f = d->m_func;
//try {
if (d->m_pars.size() == 0)
f->eval(t, ydata, ydotdata, NULL);
else
f->eval(t, ydata, ydotdata, DATA_PTR(d->m_pars));
//}
//catch (...) {
//Cantera::showErrors();
//Cantera::error("Teminating execution");
//}
return 0;
}
/**
* Function called by cvodes to evaluate ydot given y. The cvode
* integrator allows passing in a void* pointer to access
* external data. This pointer is cast to a pointer to a instance
* of class FuncEval. The equations to be integrated should be
* specified by deriving a class from FuncEval that evaluates the
* desired equations.
* @ingroup odeGroup
*/
static int cvodes_rhs(realtype t, N_Vector y, N_Vector ydot,
void *f_data) {
double* ydata = NV_DATA_S(y); //N_VDATA(y);
double* ydotdata = NV_DATA_S(ydot); //N_VDATA(ydot);
Cantera::FuncData* d = (Cantera::FuncData*)f_data;
Cantera::FuncEval* f = d->m_func;
//try {
if (d->m_pars.size() == 0)
f->eval(t, ydata, ydotdata, NULL);
else
f->eval(t, ydata, ydotdata, DATA_PTR(d->m_pars));
//}
//catch (...) {
//Cantera::showErrors();
//Cantera::error("Teminating execution");
//}
return 0;
}
}
namespace Cantera {
/**
* Constructor. Default settings: dense jacobian, no user-supplied
* Jacobian function, Newton iteration.
*/
CVodesIntegrator::CVodesIntegrator() : m_neq(0),
m_cvode_mem(0),
m_t0(0.0),
m_y(0),
m_abstol(0),
m_type(DENSE+NOJAC),
m_itol(CV_SS),
m_method(CV_BDF),
m_iter(CV_NEWTON),
m_maxord(0),
m_reltol(1.e-9),
m_abstols(1.e-15),
m_reltolsens(1.0e-5),
m_abstolsens(1.0e-4),
m_nabs(0),
m_hmax(0.0),
m_maxsteps(20000), m_np(0),
m_mupper(0), m_mlower(0)
{
//m_ropt.resize(OPT_SIZE,0.0);
//m_iopt = new long[OPT_SIZE];
//fill(m_iopt, m_iopt+OPT_SIZE,0);
/**
* Constructor. Default settings: dense jacobian, no user-supplied
* Jacobian function, Newton iteration.
*/
CVodesIntegrator::CVodesIntegrator() :
m_neq(0),
m_cvode_mem(0),
m_t0(0.0),
m_y(0),
m_abstol(0),
m_type(DENSE+NOJAC),
m_itol(CV_SS),
m_method(CV_BDF),
m_iter(CV_NEWTON),
m_maxord(0),
m_reltol(1.e-9),
m_abstols(1.e-15),
m_reltolsens(1.0e-5),
m_abstolsens(1.0e-4),
m_nabs(0),
m_hmax(0.0),
m_maxsteps(20000), m_np(0),
m_mupper(0), m_mlower(0)
{
//m_ropt.resize(OPT_SIZE,0.0);
//m_iopt = new long[OPT_SIZE];
//fill(m_iopt, m_iopt+OPT_SIZE,0);
}
/// Destructor.
CVodesIntegrator::~CVodesIntegrator()
{
if (m_cvode_mem) {
if (m_np > 0)
CVodeSensFree(m_cvode_mem);
CVodeFree(&m_cvode_mem);
}
if (m_y) N_VDestroy_Serial(nv(m_y));
if (m_abstol) N_VDestroy_Serial(nv(m_abstol));
delete m_fdata;
/// Destructor.
CVodesIntegrator::~CVodesIntegrator()
{
if (m_cvode_mem) {
if (m_np > 0)
CVodeSensFree(m_cvode_mem);
CVodeFree(&m_cvode_mem);
}
if (m_y) N_VDestroy_Serial(nv(m_y));
if (m_abstol) N_VDestroy_Serial(nv(m_abstol));
delete m_fdata;
//delete[] m_iopt;
}
//delete[] m_iopt;
}
double& CVodesIntegrator::solution(int k){
return NV_Ith_S(nv(m_y),k);
}
double& CVodesIntegrator::solution(int k){
return NV_Ith_S(nv(m_y),k);
}
double* CVodesIntegrator::solution(){ return NV_DATA_S(nv(m_y));
double* CVodesIntegrator::solution(){ return NV_DATA_S(nv(m_y));
}
void CVodesIntegrator::setTolerances(double reltol, int n, double* abstol) {
m_itol = CV_SV;
m_nabs = n;
if (n != m_neq) {
if (m_abstol) N_VDestroy_Serial(nv(m_abstol));
m_abstol = reinterpret_cast<void*>(N_VNew_Serial(n));
}
void CVodesIntegrator::setTolerances(double reltol, int n, double* abstol) {
m_itol = CV_SV;
m_nabs = n;
if (n != m_neq) {
if (m_abstol) N_VDestroy_Serial(nv(m_abstol));
m_abstol = reinterpret_cast<void*>(N_VNew_Serial(n));
}
for (int i=0; i<n; i++) {
NV_Ith_S(nv(m_abstol), i) = abstol[i];
}
m_reltol = reltol;
for (int i=0; i<n; i++) {
NV_Ith_S(nv(m_abstol), i) = abstol[i];
}
m_reltol = reltol;
}
void CVodesIntegrator::setTolerances(double reltol, double abstol) {
m_itol = CV_SS;
m_reltol = reltol;
m_abstols = abstol;
}
void CVodesIntegrator::setTolerances(double reltol, double abstol) {
m_itol = CV_SS;
m_reltol = reltol;
m_abstols = abstol;
}
void CVodesIntegrator::setSensitivityTolerances(double reltol, double abstol) {
m_reltolsens = reltol;
m_abstolsens = abstol;
}
void CVodesIntegrator::setSensitivityTolerances(double reltol, double abstol) {
m_reltolsens = reltol;
m_abstolsens = abstol;
}
void CVodesIntegrator::setProblemType(int probtype) {
m_type = probtype;
}
void CVodesIntegrator::setProblemType(int probtype) {
m_type = probtype;
}
void CVodesIntegrator::setMethod(MethodType t) {
if (t == BDF_Method)
m_method = CV_BDF;
else if (t == Adams_Method)
m_method = CV_ADAMS;
else
throw CVodesErr("unknown method");
}
void CVodesIntegrator::setMethod(MethodType t) {
if (t == BDF_Method)
m_method = CV_BDF;
else if (t == Adams_Method)
m_method = CV_ADAMS;
else
throw CVodesErr("unknown method");
}
void CVodesIntegrator::setMaxStepSize(doublereal hmax) {
m_hmax = hmax;
if (m_cvode_mem)
CVodeSetMaxStep(m_cvode_mem, hmax);
}
void CVodesIntegrator::setMaxStepSize(doublereal hmax) {
m_hmax = hmax;
if (m_cvode_mem)
CVodeSetMaxStep(m_cvode_mem, hmax);
}
void CVodesIntegrator::setMinStepSize(doublereal hmin) {
m_hmin = hmin;
if (m_cvode_mem)
CVodeSetMinStep(m_cvode_mem, hmin);
}
void CVodesIntegrator::setMinStepSize(doublereal hmin) {
m_hmin = hmin;
if (m_cvode_mem)
CVodeSetMinStep(m_cvode_mem, hmin);
}
void CVodesIntegrator::setMaxSteps(int nmax) {
m_maxsteps = nmax;
if (m_cvode_mem)
CVodeSetMaxNumSteps(m_cvode_mem, m_maxsteps);
}
void CVodesIntegrator::setMaxSteps(int nmax) {
m_maxsteps = nmax;
if (m_cvode_mem)
CVodeSetMaxNumSteps(m_cvode_mem, m_maxsteps);
}
void CVodesIntegrator::setIterator(IterType t) {
if (t == Newton_Iter)
m_iter = CV_NEWTON;
else if (t == Functional_Iter)
m_iter = CV_FUNCTIONAL;
else
throw CVodesErr("unknown iterator");
}
void CVodesIntegrator::setIterator(IterType t) {
if (t == Newton_Iter)
m_iter = CV_NEWTON;
else if (t == Functional_Iter)
m_iter = CV_FUNCTIONAL;
else
throw CVodesErr("unknown iterator");
}
void CVodesIntegrator::sensInit(double t0, FuncEval& func) {
m_np = func.nparams();
long int nv = func.neq();
void CVodesIntegrator::sensInit(double t0, FuncEval& func) {
m_np = func.nparams();
long int nv = func.neq();
doublereal* data;
int n, j;
N_Vector y;
y = N_VNew_Serial(nv);
m_yS = N_VCloneVectorArray_Serial(m_np, y);
for (n = 0; n < m_np; n++) {
data = NV_DATA_S(m_yS[n]);
for (j = 0; j < nv; j++) {
data[j] =0.0;
}
}
int flag;
flag = CVodeSensMalloc(m_cvode_mem, m_np, CV_STAGGERED, m_yS);
if (flag != CV_SUCCESS)
throw CVodesErr("Error in CVodeSensMalloc");
vector_fp atol(m_np, m_abstolsens);
double rtol = m_reltolsens;
flag = CVodeSetSensTolerances(m_cvode_mem, CV_SS, rtol, DATA_PTR(atol));
}
void CVodesIntegrator::initialize(double t0, FuncEval& func)
{
m_neq = func.neq();
m_t0 = t0;
if (m_y) {
N_VDestroy_Serial(nv(m_y)); // free solution vector if already allocated
}
m_y = reinterpret_cast<void*>(N_VNew_Serial(m_neq)); // allocate solution vector
for (int i=0; i<m_neq; i++) {
NV_Ith_S(nv(m_y), i) = 0.0;
}
// check abs tolerance array size
if (m_itol == CV_SV && m_nabs < m_neq)
throw CVodesErr("not enough absolute tolerance values specified.");
func.getInitialConditions(m_t0, m_neq, NV_DATA_S(nv(m_y)));
if (m_cvode_mem) CVodeFree(&m_cvode_mem);
m_cvode_mem = CVodeCreate(m_method, m_iter);
if (!m_cvode_mem) throw CVodesErr("CVodeCreate failed.");
int flag = 0;
if (m_itol == CV_SV) {
// vector atol
flag = CVodeMalloc(m_cvode_mem, cvodes_rhs, m_t0, nv(m_y), m_itol,
m_reltol, nv(m_abstol));
}
else {
// scalar atol
flag = CVodeMalloc(m_cvode_mem, cvodes_rhs, m_t0, nv(m_y), m_itol,
m_reltol, &m_abstols);
}
if (flag != CV_SUCCESS) {
if (flag == CV_MEM_FAIL) {
throw CVodesErr("Memory allocation failed.");
}
else if (flag == CV_ILL_INPUT) {
throw CVodesErr("Illegal value for CVodeMalloc input argument.");
}
else
throw CVodesErr("CVodeMalloc failed.");
}
if (m_type == DENSE + NOJAC) {
long int N = m_neq;
CVDense(m_cvode_mem, N);
}
else if (m_type == DIAG) {
CVDiag(m_cvode_mem);
}
else if (m_type == GMRES) {
CVSpgmr(m_cvode_mem, PREC_NONE, 0);
}
else if (m_type == BAND + NOJAC) {
long int N = m_neq;
long int nu = m_mupper;
long int nl = m_mlower;
CVBand(m_cvode_mem, N, nu, nl);
}
else {
throw CVodesErr("unsupported option");
}
// pass a pointer to func in m_data
m_fdata = new FuncData(&func, func.nparams());
//m_data = (void*)&func;
flag = CVodeSetFdata(m_cvode_mem, (void*)m_fdata);
if (flag != CV_SUCCESS)
throw CVodesErr("CVodeSetFdata failed.");
if (func.nparams() > 0) {
sensInit(t0, func);
flag = CVodeSetSensParams(m_cvode_mem, DATA_PTR(m_fdata->m_pars),
NULL, NULL);
}
// set options
if (m_maxord > 0)
flag = CVodeSetMaxOrd(m_cvode_mem, m_maxord);
if (m_maxsteps > 0)
flag = CVodeSetMaxNumSteps(m_cvode_mem, m_maxsteps);
if (m_hmax > 0)
flag = CVodeSetMaxStep(m_cvode_mem, m_hmax);
}
void CVodesIntegrator::reinitialize(double t0, FuncEval& func)
{
m_t0 = t0;
//try {
func.getInitialConditions(m_t0, m_neq, NV_DATA_S(nv(m_y)));
//}
//catch (CanteraError) {
//showErrors();
//error("Teminating execution");
//}
int result, flag;
if (m_itol == CV_SV) {
result = CVodeReInit(m_cvode_mem, cvodes_rhs, m_t0, nv(m_y),
m_itol, m_reltol,
nv(m_abstol));
}
else {
result = CVodeReInit(m_cvode_mem, cvodes_rhs, m_t0, nv(m_y),
m_itol, m_reltol,
&m_abstols);
}
if (result != CV_SUCCESS)
throw CVodesErr("CVReInit failed. result = "+int2str(result));
if (m_type == DENSE + NOJAC) {
long int N = m_neq;
CVDense(m_cvode_mem, N);
}
else if (m_type == DIAG) {
CVDiag(m_cvode_mem);
}
else if (m_type == BAND + NOJAC) {
long int N = m_neq;
long int nu = m_mupper;
long int nl = m_mlower;
CVBand(m_cvode_mem, N, nu, nl);
}
else if (m_type == GMRES) {
CVSpgmr(m_cvode_mem, PREC_NONE, 0);
}
else {
throw CVodesErr("unsupported option");
}
// set options
if (m_maxord > 0)
flag = CVodeSetMaxOrd(m_cvode_mem, m_maxord);
if (m_maxsteps > 0)
flag = CVodeSetMaxNumSteps(m_cvode_mem, m_maxsteps);
if (m_hmax > 0)
flag = CVodeSetMaxStep(m_cvode_mem, m_hmax);
}
void CVodesIntegrator::integrate(double tout)
{
double t;
int flag;
flag = CVode(m_cvode_mem, tout, nv(m_y), &t, CV_NORMAL);
if (flag != CV_SUCCESS)
throw CVodesErr(" CVodes error encountered.");
if (m_np > 0) {
CVodeGetSens(m_cvode_mem, tout, m_yS);
}
}
double CVodesIntegrator::step(double tout)
{
double t;
int flag;
flag = CVode(m_cvode_mem, tout, nv(m_y), &t, CV_ONE_STEP);
if (flag != CV_SUCCESS)
throw CVodesErr(" CVodes error encountered.");
return t;
doublereal* data;
int n, j;
N_Vector y;
y = N_VNew_Serial(nv);
m_yS = N_VCloneVectorArray_Serial(m_np, y);
for (n = 0; n < m_np; n++) {
data = NV_DATA_S(m_yS[n]);
for (j = 0; j < nv; j++) {
data[j] =0.0;
}
int CVodesIntegrator::nEvals() const {
long int ne;
CVodeGetNumRhsEvals(m_cvode_mem, &ne);
return ne;
//return m_iopt[NFE];
}
double CVodesIntegrator::sensitivity(int k, int p) {
if (k < 0 || k >= m_neq)
throw CVodesErr("sensitivity: k out of range ("+int2str(p)+")");
if (p < 0 || p >= m_np)
throw CVodesErr("sensitivity: p out of range ("+int2str(p)+")");
return NV_Ith_S(m_yS[p],k);
int flag;
#if defined(SUNDIALS_VERSION_22) || defined(SUNDIALS_VERSION23)
flag = CVodeSensMalloc(m_cvode_mem, m_np, CV_STAGGERED, m_yS);
if (flag != CV_SUCCESS) {
throw CVodesErr("Error in CVodeSensMalloc");
}
vector_fp atol(m_np, m_abstolsens);
double rtol = m_reltolsens;
flag = CVodeSetSensTolerances(m_cvode_mem, CV_SS, rtol, DATA_PTR(atol));
#elif defined(SUNDIALS_VERSION_24)
flag = CVodeSensInit(m_cvode_mem, m_np, CV_STAGGERED,
CVSensRhsFn (0), m_yS);
if (flag != CV_SUCCESS) {
throw CVodesErr("Error in CVodeSensMalloc");
}
vector_fp atol(m_np, m_abstolsens);
double rtol = m_reltolsens;
flag = CVodeSensSStolerances(m_cvode_mem, rtol, DATA_PTR(atol));
#endif
}
void CVodesIntegrator::initialize(double t0, FuncEval& func)
{
m_neq = func.neq();
m_t0 = t0;
if (m_y) {
N_VDestroy_Serial(nv(m_y)); // free solution vector if already allocated
}
m_y = reinterpret_cast<void*>(N_VNew_Serial(m_neq)); // allocate solution vector
for (int i=0; i<m_neq; i++) {
NV_Ith_S(nv(m_y), i) = 0.0;
}
// check abs tolerance array size
if (m_itol == CV_SV && m_nabs < m_neq)
throw CVodesErr("not enough absolute tolerance values specified.");
func.getInitialConditions(m_t0, m_neq, NV_DATA_S(nv(m_y)));
if (m_cvode_mem) CVodeFree(&m_cvode_mem);
/*
* Specify the method and the iteration type:
* Cantera Defaults:
* CV_BDF - Use BDF methods
* CV_NEWTON - use newton's method
*/
m_cvode_mem = CVodeCreate(m_method, m_iter);
if (!m_cvode_mem) throw CVodesErr("CVodeCreate failed.");
int flag = 0;
#if defined(SUNDIALS_VERSION_22) || defined(SUNDIALS_VERSION23)
if (m_itol == CV_SV) {
// vector atol
flag = CVodeMalloc(m_cvode_mem, cvodes_rhs, m_t0, nv(m_y), m_itol,
m_reltol, nv(m_abstol));
}
else {
// scalar atol
flag = CVodeMalloc(m_cvode_mem, cvodes_rhs, m_t0, nv(m_y), m_itol,
m_reltol, &m_abstols);
}
if (flag != CV_SUCCESS) {
if (flag == CV_MEM_FAIL) {
throw CVodesErr("Memory allocation failed.");
}
else if (flag == CV_ILL_INPUT) {
throw CVodesErr("Illegal value for CVodeMalloc input argument.");
}
else
throw CVodesErr("CVodeMalloc failed.");
}
#elif defined(SUNDIALS_VERSION_24)
flag = CVodeInit(m_cvode_mem, cvodes_rhs, m_t0, nv(m_y));
if (flag != CV_SUCCESS) {
if (flag == CV_MEM_FAIL) {
throw CVodesErr("Memory allocation failed.");
} else if (flag == CV_ILL_INPUT) {
throw CVodesErr("Illegal value for CVodeInit input argument.");
} else {
throw CVodesErr("CVodeInit failed.");
}
}
if (m_itol == CV_SV) {
flag = CVodeSVtolerances(m_cvode_mem, m_reltol, nv(m_abstol));
} else {
flag = CVodeSStolerances(m_cvode_mem, m_reltol, m_abstols);
}
if (flag != CV_SUCCESS) {
if (flag == CV_MEM_FAIL) {
throw CVodesErr("Memory allocation failed.");
} else if (flag == CV_ILL_INPUT) {
throw CVodesErr("Illegal value for CVodeInit input argument.");
} else {
throw CVodesErr("CVodeInit failed.");
}
}
#else
printf("unknown sundials verson\n");
exit(-1);
#endif
if (m_type == DENSE + NOJAC) {
long int N = m_neq;
CVDense(m_cvode_mem, N);
}
else if (m_type == DIAG) {
CVDiag(m_cvode_mem);
}
else if (m_type == GMRES) {
CVSpgmr(m_cvode_mem, PREC_NONE, 0);
}
else if (m_type == BAND + NOJAC) {
long int N = m_neq;
long int nu = m_mupper;
long int nl = m_mlower;
CVBand(m_cvode_mem, N, nu, nl);
}
else {
throw CVodesErr("unsupported option");
}
// pass a pointer to func in m_data
m_fdata = new FuncData(&func, func.nparams());
//m_data = (void*)&func;
#if defined(SUNDIALS_VERSION_22) || defined(SUNDIALS_VERSION23)
flag = CVodeSetFdata(m_cvode_mem, (void*)m_fdata);
if (flag != CV_SUCCESS) {
throw CVodesErr("CVodeSetFdata failed.");
}
#elif defined(SUNDIALS_VERSION_24)
flag = CVodeSetUserData(m_cvode_mem, (void*)m_fdata);
if (flag != CV_SUCCESS) {
throw CVodesErr("CVodeSetUserData failed.");
}
#endif
if (func.nparams() > 0) {
sensInit(t0, func);
flag = CVodeSetSensParams(m_cvode_mem, DATA_PTR(m_fdata->m_pars),
NULL, NULL);
}
// set options
if (m_maxord > 0)
flag = CVodeSetMaxOrd(m_cvode_mem, m_maxord);
if (m_maxsteps > 0)
flag = CVodeSetMaxNumSteps(m_cvode_mem, m_maxsteps);
if (m_hmax > 0)
flag = CVodeSetMaxStep(m_cvode_mem, m_hmax);
}
void CVodesIntegrator::reinitialize(double t0, FuncEval& func)
{
m_t0 = t0;
//try {
func.getInitialConditions(m_t0, m_neq, NV_DATA_S(nv(m_y)));
//}
//catch (CanteraError) {
//showErrors();
//error("Teminating execution");
//}
int result, flag;
#if defined(SUNDIALS_VERSION_22) || defined(SUNDIALS_VERSION23)
if (m_itol == CV_SV) {
result = CVodeReInit(m_cvode_mem, cvodes_rhs, m_t0, nv(m_y),
m_itol, m_reltol,
nv(m_abstol));
}
else {
result = CVodeReInit(m_cvode_mem, cvodes_rhs, m_t0, nv(m_y),
m_itol, m_reltol,
&m_abstols);
}
if (result != CV_SUCCESS) {
throw CVodesErr("CVodeReInit failed. result = "+int2str(result));
}
#elif defined(SUNDIALS_VERSION_24)
result = CVodeReInit(m_cvode_mem, m_t0, nv(m_y));
if (result != CV_SUCCESS) {
throw CVodesErr("CVodeReInit failed. result = "+int2str(result));
}
#endif
if (m_type == DENSE + NOJAC) {
long int N = m_neq;
CVDense(m_cvode_mem, N);
}
else if (m_type == DIAG) {
CVDiag(m_cvode_mem);
}
else if (m_type == BAND + NOJAC) {
long int N = m_neq;
long int nu = m_mupper;
long int nl = m_mlower;
CVBand(m_cvode_mem, N, nu, nl);
}
else if (m_type == GMRES) {
CVSpgmr(m_cvode_mem, PREC_NONE, 0);
}
else {
throw CVodesErr("unsupported option");
}
// set options
if (m_maxord > 0)
flag = CVodeSetMaxOrd(m_cvode_mem, m_maxord);
if (m_maxsteps > 0)
flag = CVodeSetMaxNumSteps(m_cvode_mem, m_maxsteps);
if (m_hmax > 0)
flag = CVodeSetMaxStep(m_cvode_mem, m_hmax);
}
void CVodesIntegrator::integrate(double tout)
{
double t;
int flag;
double tretn;
flag = CVode(m_cvode_mem, tout, nv(m_y), &t, CV_NORMAL);
if (flag != CV_SUCCESS)
throw CVodesErr(" CVodes error encountered.");
#if defined(SUNDIALS_VERSION_22) || defined(SUNDIALS_VERSION23)
if (m_np > 0) {
CVodeGetSens(m_cvode_mem, tout, m_yS);
}
#elif defined(SUNDIALS_VERSION_24)
if (m_np > 0) {
CVodeGetSens(m_cvode_mem, &tretn, m_yS);
if (fabs(tretn - tout) > 1.0E-5) {
throw CVodesErr("Time of Sensitivities different than time of tout");
}
}
#endif
}
double CVodesIntegrator::step(double tout)
{
double t;
int flag;
flag = CVode(m_cvode_mem, tout, nv(m_y), &t, CV_ONE_STEP);
if (flag != CV_SUCCESS)
throw CVodesErr(" CVodes error encountered.");
return t;
}
int CVodesIntegrator::nEvals() const {
long int ne;
CVodeGetNumRhsEvals(m_cvode_mem, &ne);
return ne;
//return m_iopt[NFE];
}
double CVodesIntegrator::sensitivity(int k, int p) {
if (k < 0 || k >= m_neq)
throw CVodesErr("sensitivity: k out of range ("+int2str(p)+")");
if (p < 0 || p >= m_np)
throw CVodesErr("sensitivity: p out of range ("+int2str(p)+")");
return NV_Ith_S(m_yS[p],k);
}
}

1
config.h.in
View file

@ -47,6 +47,7 @@ typedef int ftnlen; // Fortran hidden string length type
#undef HAS_SUNDIALS
#undef SUNDIALS_VERSION_22
#undef SUNDIALS_VERSION_23
#undef SUNDIALS_VERSION_24
//-------- LAPACK / BLAS ---------

9550
configure vendored
View file

File diff suppressed because one or more lines are too long

48
configure.in
View file

@ -365,29 +365,45 @@ sundials_lib_dep=
if test ${use_sundials} = 1; then
AC_DEFINE(HAS_SUNDIALS)
echo "using CVODES from SUNDIALS... Sensitivity analysis enabled."
CVODE_LIBS='-lsundials_cvodes -lsundials_nvecserial'
IDA_LIBS='-lsundials_ida -lsundials_nvecserial'
if test "$SUNDIALS_VERSION" = "2.2"; then
AC_DEFINE(SUNDIALS_VERSION_22)
sundials_include='-I'${SUNDIALS_HOME}'/include -I'${SUNDIALS_HOME}'/include/sundials -I'${SUNDIALS_HOME}'/include/cvodes -I'${SUNDIALS_HOME}'/include/ida'
echo "sundials include directory: " ${sundials_include}
echo "sundials library directory: " $SUNDIALS_LIB_DIR
sundials_lib_dir=$SUNDIALS_LIB_DIR
sundials_lib="-L$SUNDIALS_LIB_DIR -lsundials_cvodes -lsundials_ida -lsundials_nvecserial"
sundials_lib_dep="$SUNDIALS_LIB_DIR/libsundials_cvodes.a $SUNDIALS_LIB_DIR/libsundials_ida.a $SUNDIALS_LIB_DIR/libsundials_nvecserial.a"
AC_DEFINE(HAS_SUNDIALS)
AC_DEFINE(SUNDIALS_VERSION_22)
sundials_include='-I'${SUNDIALS_HOME}'/include -I'${SUNDIALS_HOME}'/include/sundials -I'${SUNDIALS_HOME}'/include/cvodes -I'${SUNDIALS_HOME}'/include/ida'
echo "sundials include directory: " ${sundials_include}
echo "sundials library directory: " $SUNDIALS_LIB_DIR
sundials_lib_dir=$SUNDIALS_LIB_DIR
sundials_lib="-L$SUNDIALS_LIB_DIR -lsundials_cvodes -lsundials_ida -lsundials_nvecserial"
sundials_lib_dep="$SUNDIALS_LIB_DIR/libsundials_cvodes.a $SUNDIALS_LIB_DIR/libsundials_ida.a $SUNDIALS_LIB_DIR/libsundials_nvecserial.a"
elif test "$SUNDIALS_VERSION" = "2.3"; then
AC_DEFINE(HAS_SUNDIALS)
AC_DEFINE(SUNDIALS_VERSION_23)
sundials_include='-I'${SUNDIALS_INC_DIR}
echo "sundials include directory: " ${sundials_include}
echo "sundials library directory: " $SUNDIALS_LIB_DIR
sundials_lib_dir=$SUNDIALS_LIB_DIR
sundials_lib="-L$SUNDIALS_LIB_DIR -lsundials_cvodes -lsundials_ida -lsundials_nvecserial"
sundials_lib_dep="$SUNDIALS_LIB_DIR/libsundials_cvodes.a $SUNDIALS_LIB_DIR/libsundials_ida.a $SUNDIALS_LIB_DIR/libsundials_nvecserial.a"
# python tools/src/sundials_version.py $SUNDIALS_HOME
elif test "$SUNDIALS_VERSION" = "2.4"; then
AC_DEFINE(HAS_SUNDIALS)
AC_DEFINE(SUNDIALS_VERSION_24)
sundials_include='-I'${SUNDIALS_INC_DIR}
echo "sundials include directory: " ${sundials_include}
echo "sundials library directory: " $SUNDIALS_LIB_DIR
sundials_lib_dir=$SUNDIALS_LIB_DIR
sundials_lib="-L$SUNDIALS_LIB_DIR -lsundials_cvodes -lsundials_ida -lsundials_nvecserial"
sundials_lib_dep="$SUNDIALS_LIB_DIR/libsundials_cvodes.a $SUNDIALS_LIB_DIR/libsundials_ida.a $SUNDIALS_LIB_DIR/libsundials_nvecserial.a"
# python tools/src/sundials_version.py $SUNDIALS_HOME
else
AC_DEFINE(SUNDIALS_VERSION_23)
sundials_include='-I'${SUNDIALS_INC_DIR}
echo "sundials include directory: " ${sundials_include}
echo "sundials library directory: " $SUNDIALS_LIB_DIR
sundials_lib_dir=$SUNDIALS_LIB_DIR
sundials_lib="-L$SUNDIALS_LIB_DIR -lsundials_cvodes -lsundials_ida -lsundials_nvecserial"
sundials_lib_dep="$SUNDIALS_LIB_DIR/libsundials_cvodes.a $SUNDIALS_LIB_DIR/libsundials_ida.a $SUNDIALS_LIB_DIR/libsundials_nvecserial.a"
# python tools/src/sundials_version.py $SUNDIALS_HOME
echo "ERROR: unknown or unsupported sundials version #: $SUNDIALS_VERSION"
echo " Supported versions are 2.2, 2.3, and 2.4"
echo " Please fix or turn off the sundials option by setting USE_SUNDIALS to no"
use_sundials=0
fi
fi

8
preconfig
View file

@ -315,11 +315,11 @@ USE_SUNDIALS=${USE_SUNDIALS:='default'}
# It is recommended that you install the newest release of sundials
# (currently 2.3.0) before building Cantera. But if you want to use an
# (currently 2.4.0) before building Cantera. But if you want to use an
# older version, set SUNDIALS_VERSION to the version you have.
# Acceptable values are '2.2' and '2.3' only; anything else will cause
# Cantera to not use sundials.
SUNDIALS_VERSION=${SUNDIALS_VERSION:='2.3'}
# Acceptable values are '2.2', '2.3', or '2.4' ; anything else will cause
# Cantera to
SUNDIALS_VERSION=${SUNDIALS_VERSION:='2.4'}
#-----------------------------------------------------------------
# BLAS and LAPACK