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

Added support for Sundials 2.5

Browse source
This commit is contained in:
Ray Speth 2012-05-25 23:13:51 +00:00
parent 396d34e1ab
commit 31ffb5e400
6 changed files with 42 additions and 48 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

14
SConstruct
View file

@ -864,14 +864,14 @@ elif env['use_sundials'] == 'y' and not env['HAS_SUNDIALS']:
print "ERROR: Unable to find Sundials headers and / or libraries."
print "See config.log for details."
sys.exit(1)
elif env['use_sundials'] == 'y' and env['sundials_version'] not in ('2.2','2.3','2.4'):
print """ERROR: Sundials version %r is not supported."""
elif env['use_sundials'] == 'y' and env['sundials_version'] not in ('2.2','2.3','2.4','2.5'):
print """ERROR: Sundials version %r is not supported.""" % env['sundials_version']
sys.exit(1)
# Deprecation warnings for old Sundials versions
if env.get('sundials_version') in ('2.2', '2.3'):
print 'WARNING: Support for Sundials %s is deprecated and will be removed.' % env['sundials_version']
print 'WARNING: Upgrading to Sundials 2.4 is strongly recommended.'
print 'WARNING: Upgrading to Sundials 2.5 is strongly recommended.'
# Deprecation warnings for numarray and numeric
if env.get('python_array') in ('numarray', 'numeric'):
@ -961,10 +961,10 @@ cdefine('HAS_NUMERIC', 'python_array', 'numeric')
cdefine('HAS_NO_PYTHON', 'python_package', 'none')
cdefine('HAS_SUNDIALS', 'use_sundials', 'y')
if env['use_sundials']:
cdefine('SUNDIALS_VERSION_22', 'sundials_version', '2.2')
cdefine('SUNDIALS_VERSION_23', 'sundials_version', '2.3')
cdefine('SUNDIALS_VERSION_24', 'sundials_version', '2.4')
if env['use_sundials'] == 'y':
configh['SUNDIALS_VERSION'] = env['sundials_version'].replace('.','')
else:
configh['SUNDIALS_VERSION'] = 0
cdefine('WITH_ELECTROLYTES', 'with_electrolytes')
cdefine('WITH_IDEAL_SOLUTIONS', 'with_ideal_solutions')

4
include/cantera/base/config.h.in
View file

@ -50,9 +50,7 @@ typedef int ftnlen; // Fortran hidden string length type
%(HAS_SUNDIALS)s
%(SUNDIALS_VERSION_22)s
%(SUNDIALS_VERSION_23)s
%(SUNDIALS_VERSION_24)s
%(SUNDIALS_VERSION)s
//-------- LAPACK / BLAS ---------

6
include/cantera/numerics/IDA_Solver.h
View file

@ -13,9 +13,9 @@
#include "DAE_Solver.h"
#include "cantera/base/ctexceptions.h"
#ifdef HAS_SUNDIALS
#if HAS_SUNDIALS
#ifdef SUNDIALS_VERSION_22
#if SUNDIALS_VERSION == 22
#include "nvector_serial.h"
#else
#include "sundials/sundials_nvector.h"
@ -27,7 +27,7 @@
#define IDA_WF 3
#endif
#if defined(SUNDIALS_VERSION_24)
#if SUNDIALS_VERSION >= 24
#define REAL_WORKSPACE_SIZE 0
#endif

37
src/numerics/CVodesIntegrator.cpp
View file

@ -12,7 +12,7 @@
#include <iostream>
using namespace std;
#ifdef SUNDIALS_VERSION_22
#if SUNDIALS_VERSION == 22
#include "sundials_types.h"
#include "sundials_math.h"
@ -25,7 +25,7 @@ using namespace std;
#else
#if defined(SUNDIALS_VERSION_23) || defined (SUNDIALS_VERSION_24)
#if SUNDIALS_VERSION >= 23
#include "sundials/sundials_types.h"
#include "sundials/sundials_math.h"
#include "sundials/sundials_nvector.h"
@ -37,15 +37,12 @@ using namespace std;
#include "cvodes/cvodes_band.h"
#else
unsupported sundials version!
#error unsupported Sundials version!
#endif
#if defined (SUNDIALS_VERSION_24)
#if SUNDIALS_VERSION >= 24
#define CV_SS 1
#define CV_SV 2
#endif
#endif
@ -270,7 +267,7 @@ void CVodesIntegrator::sensInit(double t0, FuncEval& func)
int flag;
#if defined(SUNDIALS_VERSION_22) || defined(SUNDIALS_VERSION_23)
#if SUNDIALS_VERSION <= 23
flag = CVodeSensMalloc(m_cvode_mem, m_np, CV_STAGGERED, m_yS);
if (flag != CV_SUCCESS) {
throw CVodesErr("Error in CVodeSensMalloc");
@ -278,7 +275,7 @@ void CVodesIntegrator::sensInit(double t0, FuncEval& func)
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)
#elif SUNDIALS_VERSION >= 24
flag = CVodeSensInit(m_cvode_mem, m_np, CV_STAGGERED,
CVSensRhsFn(0), m_yS);
@ -288,7 +285,6 @@ void CVodesIntegrator::sensInit(double t0, FuncEval& func)
vector_fp atol(m_np, m_abstolsens);
double rtol = m_reltolsens;
flag = CVodeSensSStolerances(m_cvode_mem, rtol, DATA_PTR(atol));
#endif
}
@ -328,7 +324,7 @@ void CVodesIntegrator::initialize(double t0, FuncEval& func)
}
int flag = 0;
#if defined(SUNDIALS_VERSION_22) || defined(SUNDIALS_VERSION_23)
#if SUNDIALS_VERSION <= 23
if (m_itol == CV_SV) {
// vector atol
flag = CVodeMalloc(m_cvode_mem, cvodes_rhs, m_t0, nv(m_y), m_itol,
@ -347,7 +343,7 @@ void CVodesIntegrator::initialize(double t0, FuncEval& func)
throw CVodesErr("CVodeMalloc failed.");
}
}
#elif defined(SUNDIALS_VERSION_24)
#elif SUNDIALS_VERSION >= 24
flag = CVodeInit(m_cvode_mem, cvodes_rhs, m_t0, nv(m_y));
if (flag != CV_SUCCESS) {
@ -374,13 +370,8 @@ void CVodesIntegrator::initialize(double t0, FuncEval& func)
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);
@ -401,12 +392,12 @@ void CVodesIntegrator::initialize(double t0, FuncEval& func)
m_fdata = new FuncData(&func, func.nparams());
//m_data = (void*)&func;
#if defined(SUNDIALS_VERSION_22) || defined(SUNDIALS_VERSION_23)
#if SUNDIALS_VERSION <= 23
flag = CVodeSetFdata(m_cvode_mem, (void*)m_fdata);
if (flag != CV_SUCCESS) {
throw CVodesErr("CVodeSetFdata failed.");
}
#elif defined(SUNDIALS_VERSION_24)
#elif SUNDIALS_VERSION >= 24
flag = CVodeSetUserData(m_cvode_mem, (void*)m_fdata);
if (flag != CV_SUCCESS) {
throw CVodesErr("CVodeSetUserData failed.");
@ -444,7 +435,7 @@ void CVodesIntegrator::reinitialize(double t0, FuncEval& func)
int result;
#if defined(SUNDIALS_VERSION_22) || defined(SUNDIALS_VERSION_23)
#if SUNDIALS_VERSION <= 23
if (m_itol == CV_SV) {
result = CVodeReInit(m_cvode_mem, cvodes_rhs, m_t0, nv(m_y),
m_itol, m_reltol,
@ -457,7 +448,7 @@ void CVodesIntegrator::reinitialize(double t0, FuncEval& func)
if (result != CV_SUCCESS) {
throw CVodesErr("CVodeReInit failed. result = "+int2str(result));
}
#elif defined(SUNDIALS_VERSION_24)
#elif SUNDIALS_VERSION >= 24
result = CVodeReInit(m_cvode_mem, m_t0, nv(m_y));
if (result != CV_SUCCESS) {
throw CVodesErr("CVodeReInit failed. result = "+int2str(result));
@ -501,11 +492,11 @@ void CVodesIntegrator::integrate(double tout)
if (flag != CV_SUCCESS) {
throw CVodesErr(" CVodes error encountered. Error code: " + int2str(flag));
}
#if defined(SUNDIALS_VERSION_22) || defined(SUNDIALS_VERSION_23)
#if SUNDIALS_VERSION <= 23
if (m_np > 0) {
CVodeGetSens(m_cvode_mem, tout, m_yS);
}
#elif defined(SUNDIALS_VERSION_24)
#elif SUNDIALS_VERSION >= 24
double tretn;
if (m_np > 0) {
CVodeGetSens(m_cvode_mem, &tretn, m_yS);

2
src/numerics/CVodesIntegrator.h
View file

@ -13,7 +13,7 @@
#include "cantera/base/ctexceptions.h"
#include "cantera/base/ct_defs.h"
#ifdef SUNDIALS_VERSION_22
#if SUNDIALS_VERSION == 22
#include "nvector_serial.h"
#else
#include "sundials/sundials_nvector.h"

27
src/numerics/IDA_Solver.cpp
View file

@ -10,7 +10,7 @@
#include <iostream>
#ifdef SUNDIALS_VERSION_24
#if SUNDIALS_VERSION >= 24
#include "sundials/sundials_types.h"
#include "sundials/sundials_math.h"
#include "ida/ida.h"
@ -21,6 +21,11 @@
using namespace std;
#if SUNDIALS_VERSION < 25
typedef int sd_size_t;
#else
typedef long int sd_size_t;
#endif
inline static N_Vector nv(void* x)
@ -93,7 +98,7 @@ extern "C" {
/*!
*
*
* typedef int (*IDADlsDenseJacFn)(int N, realtype t, realtype c_j,
* typedef int (*IDADlsDenseJacFn)(sd_size_t N, realtype t, realtype c_j,
* N_Vector y, N_Vector yp, N_Vector r,
* DlsMat Jac, void *user_data,
* N_Vector tmp1, N_Vector tmp2, N_Vector tmp3);
@ -105,7 +110,7 @@ extern "C" {
* In the case of a recoverable error return, the integrator will
* attempt to recover by reducing the stepsize (which changes cj).
*/
static int ida_jacobian(int nrows, realtype t, realtype c_j, N_Vector y, N_Vector ydot, N_Vector r,
static int ida_jacobian(sd_size_t nrows, realtype t, realtype c_j, N_Vector y, N_Vector ydot, N_Vector r,
DlsMat Jac, void* f_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3)
{
doublereal* ydata = NV_DATA_S(y);
@ -367,7 +372,7 @@ void IDA_Solver::init(doublereal t0)
if (m_itol == IDA_SV) {
#if defined(SUNDIALS_VERSION_22) || defined(SUNDIALS_VERSION_23)
#if SUNDIALS_VERSION <= 23
// vector atol
flag = IDAMalloc(m_ida_mem, ida_resid, m_t0, nv(m_y), nv(m_ydot),
m_itol, m_reltol, nv(m_abstol));
@ -381,7 +386,7 @@ void IDA_Solver::init(doublereal t0)
}
}
#elif defined(SUNDIALS_VERSION_24)
#elif SUNDIALS_VERSION >= 24
flag = IDAInit(m_ida_mem, ida_resid, m_t0, nv(m_y), nv(m_ydot));
if (flag != IDA_SUCCESS) {
if (flag == IDA_MEM_FAIL) {
@ -398,7 +403,7 @@ void IDA_Solver::init(doublereal t0)
}
#endif
} else {
#if defined(SUNDIALS_VERSION_22) || defined(SUNDIALS_VERSION_23)
#if SUNDIALS_VERSION <= 23
// scalar atol
flag = IDAMalloc(m_ida_mem, ida_resid, m_t0, nv(m_y), nv(m_ydot),
m_itol, m_reltol, &m_abstols);
@ -412,7 +417,7 @@ void IDA_Solver::init(doublereal t0)
}
}
#elif defined(SUNDIALS_VERSION_24)
#elif SUNDIALS_VERSION >= 24
flag = IDAInit(m_ida_mem, ida_resid, m_t0, nv(m_y), nv(m_ydot));
if (flag != IDA_SUCCESS) {
if (flag == IDA_MEM_FAIL) {
@ -458,12 +463,12 @@ void IDA_Solver::init(doublereal t0)
// pass a pointer to func in m_data
m_fdata = new ResidData(&m_resid, this, m_resid.nparams());
#if defined(SUNDIALS_VERSION_22) || defined(SUNDIALS_VERSION_23)
#if SUNDIALS_VERSION <= 23
flag = IDASetRdata(m_ida_mem, (void*)m_fdata);
if (flag != IDA_SUCCESS) {
throw IDA_Err("IDASetRdata failed.");
}
#elif defined(SUNDIALS_VERSION_24)
#elif SUNDIALS_VERSION >= 24
flag = IDASetUserData(m_ida_mem, (void*)m_fdata);
if (flag != IDA_SUCCESS) {
throw IDA_Err("IDASetUserData failed.");
@ -549,7 +554,7 @@ void IDA_Solver::correctInitial_Y_given_Yp(doublereal* y, doublereal* yp, doubl
}
flag = IDAGetSolution(m_ida_mem, tout1, nv(m_y), nv(m_ydot));
flag = IDAGetConsistentIC(m_ida_mem, nv(m_y), nv(m_ydot));
if (flag != IDA_SUCCESS) {
throw IDA_Err("IDAGetSolution failed: error = " + int2str(flag));
}
@ -595,7 +600,7 @@ void IDA_Solver::correctInitial_YaYp_given_Yd(doublereal* y, doublereal* yp, dou
}
flag = IDAGetSolution(m_ida_mem, tout1, nv(m_y), nv(m_ydot));
flag = IDAGetConsistentIC(m_ida_mem, nv(m_y), nv(m_ydot));
if (flag != IDA_SUCCESS) {
throw IDA_Err("IDAGetSolution failed: error = " + int2str(flag));
}