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

Merged Revisision_1.8.0 changes into the trunk.

Browse source
This commit is contained in:
Harry Moffat 2009-09-22 17:26:18 +00:00
commit 15c9ca0294
7 changed files with 532 additions and 403 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 #else
#ifdef SUNDIALS_VERSION_23 #if defined(SUNDIALS_VERSION_23) || defined (SUNDIALS_VERSION_24)
#include <sundials/sundials_types.h> #include <sundials/sundials_types.h>
#include <sundials/sundials_math.h> #include <sundials/sundials_math.h>
#include <sundials/sundials_nvector.h> #include <sundials/sundials_nvector.h>
@ -42,383 +42,471 @@ unsupported sundials version!
#endif #endif
#if defined (SUNDIALS_VERSION_24)
#define CV_SS 1
#define CV_SV 2
#endif
#endif #endif
inline static N_Vector nv(void* x) { inline static N_Vector nv(void* x) {
return reinterpret_cast<N_Vector>(x); return reinterpret_cast<N_Vector>(x);
} }
namespace Cantera { namespace Cantera {
class FuncData { class FuncData {
public: public:
FuncData(FuncEval* f, int npar = 0) { FuncData(FuncEval* f, int npar = 0) {
m_pars.resize(npar, 1.0); m_pars.resize(npar, 1.0);
m_func = f; m_func = f;
} }
virtual ~FuncData() {} virtual ~FuncData() {}
vector_fp m_pars; vector_fp m_pars;
FuncEval* m_func; FuncEval* m_func;
}; };
} }
extern "C" { extern "C" {
/** /**
* Function called by cvodes to evaluate ydot given y. The cvode * Function called by cvodes to evaluate ydot given y. The cvode
* integrator allows passing in a void* pointer to access * integrator allows passing in a void* pointer to access
* external data. This pointer is cast to a pointer to a instance * external data. This pointer is cast to a pointer to a instance
* of class FuncEval. The equations to be integrated should be * of class FuncEval. The equations to be integrated should be
* specified by deriving a class from FuncEval that evaluates the * specified by deriving a class from FuncEval that evaluates the
* desired equations. * desired equations.
* @ingroup odeGroup * @ingroup odeGroup
*/ */
static int cvodes_rhs(realtype t, N_Vector y, N_Vector ydot, static int cvodes_rhs(realtype t, N_Vector y, N_Vector ydot,
void *f_data) { void *f_data) {
double* ydata = NV_DATA_S(y); //N_VDATA(y); double* ydata = NV_DATA_S(y); //N_VDATA(y);
double* ydotdata = NV_DATA_S(ydot); //N_VDATA(ydot); double* ydotdata = NV_DATA_S(ydot); //N_VDATA(ydot);
Cantera::FuncData* d = (Cantera::FuncData*)f_data; Cantera::FuncData* d = (Cantera::FuncData*)f_data;
Cantera::FuncEval* f = d->m_func; Cantera::FuncEval* f = d->m_func;
//try { //try {
if (d->m_pars.size() == 0) if (d->m_pars.size() == 0)
f->eval(t, ydata, ydotdata, NULL); f->eval(t, ydata, ydotdata, NULL);
else else
f->eval(t, ydata, ydotdata, DATA_PTR(d->m_pars)); f->eval(t, ydata, ydotdata, DATA_PTR(d->m_pars));
//} //}
//catch (...) { //catch (...) {
//Cantera::showErrors(); //Cantera::showErrors();
//Cantera::error("Teminating execution"); //Cantera::error("Teminating execution");
//} //}
return 0; return 0;
} }
} }
namespace Cantera { namespace Cantera {
/** /**
* Constructor. Default settings: dense jacobian, no user-supplied * Constructor. Default settings: dense jacobian, no user-supplied
* Jacobian function, Newton iteration. * Jacobian function, Newton iteration.
*/ */
CVodesIntegrator::CVodesIntegrator() : m_neq(0), CVodesIntegrator::CVodesIntegrator() :
m_cvode_mem(0), m_neq(0),
m_t0(0.0), m_cvode_mem(0),
m_y(0), m_t0(0.0),
m_abstol(0), m_y(0),
m_type(DENSE+NOJAC), m_abstol(0),
m_itol(CV_SS), m_type(DENSE+NOJAC),
m_method(CV_BDF), m_itol(CV_SS),
m_iter(CV_NEWTON), m_method(CV_BDF),
m_maxord(0), m_iter(CV_NEWTON),
m_reltol(1.e-9), m_maxord(0),
m_abstols(1.e-15), m_reltol(1.e-9),
m_reltolsens(1.0e-5), m_abstols(1.e-15),
m_abstolsens(1.0e-4), m_reltolsens(1.0e-5),
m_nabs(0), m_abstolsens(1.0e-4),
m_hmax(0.0), m_nabs(0),
m_maxsteps(20000), m_np(0), m_hmax(0.0),
m_mupper(0), m_mlower(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]; //m_ropt.resize(OPT_SIZE,0.0);
//fill(m_iopt, m_iopt+OPT_SIZE,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;
//delete[] m_iopt;
/// 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;
}
double& CVodesIntegrator::solution(int k){ double& CVodesIntegrator::solution(int k){
return NV_Ith_S(nv(m_y),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));
} }
for (int i=0; i<n; i++) {
void CVodesIntegrator::setTolerances(double reltol, int n, double* abstol) { NV_Ith_S(nv(m_abstol), i) = abstol[i];
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;
} }
m_reltol = reltol;
}
void CVodesIntegrator::setTolerances(double reltol, double abstol) { void CVodesIntegrator::setTolerances(double reltol, double abstol) {
m_itol = CV_SS; m_itol = CV_SS;
m_reltol = reltol; m_reltol = reltol;
m_abstols = abstol; m_abstols = abstol;
} }
void CVodesIntegrator::setSensitivityTolerances(double reltol, double abstol) { void CVodesIntegrator::setSensitivityTolerances(double reltol, double abstol) {
m_reltolsens = reltol; m_reltolsens = reltol;
m_abstolsens = abstol; m_abstolsens = abstol;
} }
void CVodesIntegrator::setProblemType(int probtype) { void CVodesIntegrator::setProblemType(int probtype) {
m_type = probtype; m_type = probtype;
} }
void CVodesIntegrator::setMethod(MethodType t) { void CVodesIntegrator::setMethod(MethodType t) {
if (t == BDF_Method) if (t == BDF_Method)
m_method = CV_BDF; m_method = CV_BDF;
else if (t == Adams_Method) else if (t == Adams_Method)
m_method = CV_ADAMS; m_method = CV_ADAMS;
else else
throw CVodesErr("unknown method"); throw CVodesErr("unknown method");
} }
void CVodesIntegrator::setMaxStepSize(doublereal hmax) { void CVodesIntegrator::setMaxStepSize(doublereal hmax) {
m_hmax = hmax; m_hmax = hmax;
if (m_cvode_mem) if (m_cvode_mem)
CVodeSetMaxStep(m_cvode_mem, hmax); CVodeSetMaxStep(m_cvode_mem, hmax);
} }
void CVodesIntegrator::setMinStepSize(doublereal hmin) { void CVodesIntegrator::setMinStepSize(doublereal hmin) {
m_hmin = hmin; m_hmin = hmin;
if (m_cvode_mem) if (m_cvode_mem)
CVodeSetMinStep(m_cvode_mem, hmin); CVodeSetMinStep(m_cvode_mem, hmin);
} }
void CVodesIntegrator::setMaxSteps(int nmax) { void CVodesIntegrator::setMaxSteps(int nmax) {
m_maxsteps = nmax; m_maxsteps = nmax;
if (m_cvode_mem) if (m_cvode_mem)
CVodeSetMaxNumSteps(m_cvode_mem, m_maxsteps); CVodeSetMaxNumSteps(m_cvode_mem, m_maxsteps);
} }
void CVodesIntegrator::setIterator(IterType t) { void CVodesIntegrator::setIterator(IterType t) {
if (t == Newton_Iter) if (t == Newton_Iter)
m_iter = CV_NEWTON; m_iter = CV_NEWTON;
else if (t == Functional_Iter) else if (t == Functional_Iter)
m_iter = CV_FUNCTIONAL; m_iter = CV_FUNCTIONAL;
else else
throw CVodesErr("unknown iterator"); throw CVodesErr("unknown iterator");
} }
void CVodesIntegrator::sensInit(double t0, FuncEval& func) { void CVodesIntegrator::sensInit(double t0, FuncEval& func) {
m_np = func.nparams(); m_np = func.nparams();
long int nv = func.neq(); long int nv = func.neq();
doublereal* data; doublereal* data;
int n, j; int n, j;
N_Vector y; N_Vector y;
y = N_VNew_Serial(nv); y = N_VNew_Serial(nv);
m_yS = N_VCloneVectorArray_Serial(m_np, y); m_yS = N_VCloneVectorArray_Serial(m_np, y);
for (n = 0; n < m_np; n++) { for (n = 0; n < m_np; n++) {
data = NV_DATA_S(m_yS[n]); data = NV_DATA_S(m_yS[n]);
for (j = 0; j < nv; j++) { for (j = 0; j < nv; j++) {
data[j] =0.0; 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;
} }
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) { int flag;
if (k < 0 || k >= m_neq)
throw CVodesErr("sensitivity: k out of range ("+int2str(p)+")"); #if defined(SUNDIALS_VERSION_22) || defined(SUNDIALS_VERSION23)
if (p < 0 || p >= m_np) flag = CVodeSensMalloc(m_cvode_mem, m_np, CV_STAGGERED, m_yS);
throw CVodesErr("sensitivity: p out of range ("+int2str(p)+")"); if (flag != CV_SUCCESS) {
return NV_Ith_S(m_yS[p],k); 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 HAS_SUNDIALS
#undef SUNDIALS_VERSION_22 #undef SUNDIALS_VERSION_22
#undef SUNDIALS_VERSION_23 #undef SUNDIALS_VERSION_23
#undef SUNDIALS_VERSION_24
//-------- LAPACK / BLAS --------- //-------- LAPACK / BLAS ---------

74
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 if test ${use_sundials} = 1; then
AC_DEFINE(HAS_SUNDIALS)
echo "using CVODES from SUNDIALS... Sensitivity analysis enabled." echo "using CVODES from SUNDIALS... Sensitivity analysis enabled."
CVODE_LIBS='-lsundials_cvodes -lsundials_nvecserial' CVODE_LIBS='-lsundials_cvodes -lsundials_nvecserial'
IDA_LIBS='-lsundials_ida -lsundials_nvecserial' IDA_LIBS='-lsundials_ida -lsundials_nvecserial'
if test "$SUNDIALS_VERSION" = "2.2"; then if test "$SUNDIALS_VERSION" = "2.2"; then
AC_DEFINE(SUNDIALS_VERSION_22) AC_DEFINE(HAS_SUNDIALS)
sundials_include='-I'${SUNDIALS_HOME}'/include -I'${SUNDIALS_HOME}'/include/sundials -I'${SUNDIALS_HOME}'/include/cvodes -I'${SUNDIALS_HOME}'/include/ida' AC_DEFINE(SUNDIALS_VERSION_22)
echo "sundials include directory: " ${sundials_include} sundials_include='-I'${SUNDIALS_HOME}'/include -I'${SUNDIALS_HOME}'/include/sundials -I'${SUNDIALS_HOME}'/include/cvodes -I'${SUNDIALS_HOME}'/include/ida'
echo "sundials library directory: " $SUNDIALS_LIB_DIR echo "sundials include directory: " ${sundials_include}
sundials_lib_dir=$SUNDIALS_LIB_DIR echo "sundials library directory: " $SUNDIALS_LIB_DIR
sundials_lib="-L$SUNDIALS_LIB_DIR -lsundials_cvodes -lsundials_ida -lsundials_nvecserial" sundials_lib_dir=$SUNDIALS_LIB_DIR
sundials_lib_dep="$SUNDIALS_LIB_DIR/libsundials_cvodes.a $SUNDIALS_LIB_DIR/libsundials_ida.a $SUNDIALS_LIB_DIR/libsundials_nvecserial.a" 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 else
AC_DEFINE(SUNDIALS_VERSION_23) echo "ERROR: unknown or unsupported sundials version #: $SUNDIALS_VERSION"
sundials_include='-I'${SUNDIALS_INC_DIR} echo " Supported versions are 2.2, 2.3, and 2.4"
echo "sundials include directory: " ${sundials_include} echo " Please fix or turn off the sundials option by setting USE_SUNDIALS to no"
echo "sundials library directory: " $SUNDIALS_LIB_DIR use_sundials=0
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
fi fi
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 # 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. # older version, set SUNDIALS_VERSION to the version you have.
# Acceptable values are '2.2' and '2.3' only; anything else will cause # Acceptable values are '2.2', '2.3', or '2.4' ; anything else will cause
# Cantera to not use sundials. # Cantera to
SUNDIALS_VERSION=${SUNDIALS_VERSION:='2.3'} SUNDIALS_VERSION=${SUNDIALS_VERSION:='2.4'}
#----------------------------------------------------------------- #-----------------------------------------------------------------
# BLAS and LAPACK # BLAS and LAPACK

18
test_problems/surfSolverTest/surfaceSolver2_blessed.out
View file

@ -21,9 +21,9 @@ Number of reactions = 8
3 4.8454e-06 3.0554e-06 6.8484e+04 4.0505e+04 Csoot-* 3 4.8454e-06 3.0554e-06 6.8484e+04 4.0505e+04 Csoot-*
4 2.6364e-05 2.1519e-05 1.1250e+04 5.8320e+03 Csoot-* 4 2.6364e-05 2.1519e-05 1.1250e+04 5.8320e+03 Csoot-*
5 1.3012e-03 1.2749e-03 1.4427e+03 7.3573e+02 Csoot-* 5 1.3012e-03 1.2749e-03 1.4427e+03 7.3573e+02 Csoot-*
6 4.7372e+00 4.7359e+00 5.8161e-01 2.9647e-01 Csoot-* 6 4.7372e+00 4.7359e+00 5.8160e-01 2.9647e-01 Csoot-*
7 6.1771e-08 3.1736e-08 7 6.4112e-08 3.2378e-08
FIN 7 6.1771e-08 2.1591e-11 -- success FIN 7 6.4112e-08 3.6822e-12 -- success
Gas Temperature = 1.4e+03 Gas Temperature = 1.4e+03
Gas Pressure = 1.01e+05 Gas Pressure = 1.01e+05
Gas Phase: gas (0) Gas Phase: gas (0)
@ -81,8 +81,8 @@ Sum of coverages = 1
Iter Time Del_t Damp DelX Resid Name-Time Name-Damp Iter Time Del_t Damp DelX Resid Name-Time Name-Damp
----------------------------------------------------------------------------------- -----------------------------------------------------------------------------------
1 5.3218e+03 2.7005e+03 1 5.3218e+03 2.7005e+03
2 8.9325e-06 4.0777e-06 2 8.8765e-06 4.0574e-06
FIN 2 8.9325e-06 6.2499e-12 -- success FIN 2 8.8765e-06 8.4527e-11 -- success
Gas Temperature = 1.4e+03 Gas Temperature = 1.4e+03
Gas Pressure = 1.01e+05 Gas Pressure = 1.01e+05
Gas Phase: gas (0) Gas Phase: gas (0)
@ -140,8 +140,8 @@ Sum of coverages = 1
Iter Time Del_t Damp DelX Resid Name-Time Name-Damp Iter Time Del_t Damp DelX Resid Name-Time Name-Damp
----------------------------------------------------------------------------------- -----------------------------------------------------------------------------------
1 2.1569e+05 9.5571e+04 1 2.1569e+05 9.5571e+04
2 1.7622e-04 2.0108e-04 2 7.8671e-04 4.6185e-04
FIN 2 1.7622e-04 2.0335e-10 -- success FIN 2 7.8671e-04 1.5306e-10 -- success
Gas Temperature = 1.5e+03 Gas Temperature = 1.5e+03
Gas Pressure = 1.01e+05 Gas Pressure = 1.01e+05
Gas Phase: gas (0) Gas Phase: gas (0)
@ -198,8 +198,8 @@ Sum of coverages = 1
Iter Time Del_t Damp DelX Resid Name-Time Name-Damp Iter Time Del_t Damp DelX Resid Name-Time Name-Damp
----------------------------------------------------------------------------------- -----------------------------------------------------------------------------------
1 2.8875e-10 1.4792e-10 1 3.0635e-10 1.5694e-10
FIN 1 2.8875e-10 1.2324e-10 -- success FIN 1 3.0635e-10 1.2314e-10 -- success
Gas Temperature = 1.5e+03 Gas Temperature = 1.5e+03
Gas Pressure = 1.01e+05 Gas Pressure = 1.01e+05
Gas Phase: gas (0) Gas Phase: gas (0)

18
test_problems/surfSolverTest/surfaceSolver_blessed.out
View file

@ -19,9 +19,9 @@ Number of reactions = 8
3 4.8454e-06 3.0554e-06 6.8484e+04 4.0505e+04 Csoot-* 3 4.8454e-06 3.0554e-06 6.8484e+04 4.0505e+04 Csoot-*
4 2.6364e-05 2.1519e-05 1.1250e+04 5.8320e+03 Csoot-* 4 2.6364e-05 2.1519e-05 1.1250e+04 5.8320e+03 Csoot-*
5 1.3012e-03 1.2749e-03 1.4427e+03 7.3573e+02 Csoot-* 5 1.3012e-03 1.2749e-03 1.4427e+03 7.3573e+02 Csoot-*
6 4.7372e+00 4.7359e+00 5.8161e-01 2.9647e-01 Csoot-* 6 4.7372e+00 4.7359e+00 5.8160e-01 2.9647e-01 Csoot-*
7 6.1771e-08 3.1736e-08 7 6.4112e-08 3.2378e-08
FIN 7 6.1771e-08 2.1591e-11 -- success FIN 7 6.4112e-08 3.6822e-12 -- success
Gas Temperature = 1.4e+03 Gas Temperature = 1.4e+03
Gas Pressure = 1.01e+05 Gas Pressure = 1.01e+05
Gas Phase: gas (0) Gas Phase: gas (0)
@ -72,8 +72,8 @@ Sum of coverages = 1
Iter Time Del_t Damp DelX Resid Name-Time Name-Damp Iter Time Del_t Damp DelX Resid Name-Time Name-Damp
----------------------------------------------------------------------------------- -----------------------------------------------------------------------------------
1 5.3218e+03 2.7005e+03 1 5.3218e+03 2.7005e+03
2 8.9325e-06 4.0777e-06 2 8.8765e-06 4.0574e-06
FIN 2 8.9325e-06 6.2499e-12 -- success FIN 2 8.8765e-06 8.4527e-11 -- success
Gas Temperature = 1.4e+03 Gas Temperature = 1.4e+03
Gas Pressure = 1.01e+05 Gas Pressure = 1.01e+05
Gas Phase: gas (0) Gas Phase: gas (0)
@ -124,8 +124,8 @@ Sum of coverages = 1
Iter Time Del_t Damp DelX Resid Name-Time Name-Damp Iter Time Del_t Damp DelX Resid Name-Time Name-Damp
----------------------------------------------------------------------------------- -----------------------------------------------------------------------------------
1 2.1569e+05 9.5571e+04 1 2.1569e+05 9.5571e+04
2 1.7622e-04 2.0108e-04 2 7.8671e-04 4.6185e-04
FIN 2 1.7622e-04 2.0335e-10 -- success FIN 2 7.8671e-04 1.5306e-10 -- success
Gas Temperature = 1.5e+03 Gas Temperature = 1.5e+03
Gas Pressure = 1.01e+05 Gas Pressure = 1.01e+05
Gas Phase: gas (0) Gas Phase: gas (0)
@ -175,8 +175,8 @@ Sum of coverages = 1
Iter Time Del_t Damp DelX Resid Name-Time Name-Damp Iter Time Del_t Damp DelX Resid Name-Time Name-Damp
----------------------------------------------------------------------------------- -----------------------------------------------------------------------------------
1 2.8875e-10 1.4792e-10 1 3.0635e-10 1.5694e-10
FIN 1 2.8875e-10 1.2324e-10 -- success FIN 1 3.0635e-10 1.2314e-10 -- success
Gas Temperature = 1.5e+03 Gas Temperature = 1.5e+03
Gas Pressure = 1.01e+05 Gas Pressure = 1.01e+05
Gas Phase: gas (0) Gas Phase: gas (0)