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

Add support for Sundials 3.1.0 (also 3.0.0)

Browse source
This commit is contained in:
Ray Speth 2018-02-10 23:09:53 -05:00
parent 50a3fecb23
commit 9c50f7528c
7 changed files with 208 additions and 52 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

17
SConstruct
View file

@ -938,6 +938,9 @@ def get_expression_value(includes, expression):
s.extend(('#define Q(x) #x',
'#define QUOTE(x) Q(x)',
'#include <iostream>',
'#ifndef SUNDIALS_PACKAGE_VERSION', # name change in Sundials >= 3.0
'#define SUNDIALS_PACKAGE_VERSION SUNDIALS_VERSION',
'#endif',
'int main(int argc, char** argv) {',
' std::cout << %s << std::endl;' % expression,
' return 0;',
@ -1012,7 +1015,7 @@ if env['system_sundials'] == 'y':
# Ignore the minor version, e.g. 2.4.x -> 2.4
env['sundials_version'] = '.'.join(sundials_version.split('.')[:2])
if env['sundials_version'] not in ('2.4','2.5','2.6','2.7'):
if env['sundials_version'] not in ('2.4','2.5','2.6','2.7','3.0','3.1'):
print("""ERROR: Sundials version %r is not supported.""" % env['sundials_version'])
sys.exit(1)
print("""INFO: Using system installation of Sundials version %s.""" % sundials_version)
@ -1642,11 +1645,19 @@ linkSharedLibs = ['cantera_shared']
if env['system_sundials'] == 'y':
env['sundials_libs'] = ['sundials_cvodes', 'sundials_ida', 'sundials_nvecserial']
linkLibs.extend(('sundials_cvodes', 'sundials_ida', 'sundials_nvecserial'))
linkSharedLibs.extend(('sundials_cvodes', 'sundials_ida', 'sundials_nvecserial'))
if env['use_lapack'] and LooseVersion(env['sundials_version']) >= LooseVersion('3.0'):
if env.get('has_sundials_lapack'):
env['sundials_libs'].extend(('sundials_sunlinsollapackdense',
'sundials_sunlinsollapackband'))
else:
env['sundials_libs'].extend(('sundials_sunlinsoldense',
'sundials_sunlinsolband'))
else:
env['sundials_libs'] = []
linkLibs.extend(env['sundials_libs'])
linkSharedLibs.extend(env['sundials_libs'])
# Add LAPACK and BLAS to the link line
if env['blas_lapack_libs']:
linkLibs.extend(env['blas_lapack_libs'])

15
include/cantera/numerics/BandMatrix.h
View file

@ -16,13 +16,6 @@
namespace Cantera
{
#if CT_USE_LAPACK
typedef vector_int pivot_vector_t;
#else
typedef std::vector<long int> pivot_vector_t;
#endif
//! A class for banded matrices, involving matrix inversion processes.
//! The class is based upon the LAPACK banded storage matrix format.
/*!
@ -46,6 +39,7 @@ public:
* Create an \c 0 by \c 0 matrix, and initialize all elements to \c 0.
*/
BandMatrix();
~BandMatrix();
//! Creates a banded matrix and sets all elements to zero
/*!
@ -135,9 +129,6 @@ public:
//! Return the number of rows of storage needed for the band storage
size_t ldim() const;
//! Return a reference to the pivot vector
pivot_vector_t& ipiv();
//! Multiply A*b and write result to \c prod.
virtual void mult(const doublereal* b, doublereal* prod) const;
virtual void leftMult(const doublereal* const b, doublereal* const prod) const;
@ -288,8 +279,10 @@ protected:
//! value of zero
doublereal m_zero;
class PivData; // pImpl wrapper class
//! Pivot vector
pivot_vector_t m_ipiv;
std::unique_ptr<PivData> m_ipiv;
//! Vector of column pointers
std::vector<doublereal*> m_colPtrs;

2
include/cantera/numerics/CVodesIntegrator.h
View file

@ -82,6 +82,8 @@ private:
size_t m_neq;
void* m_cvode_mem;
void* m_linsol; //!< Sundials linear solver object
void* m_linsol_matrix; //!< matrix used by Sundials
FuncEval* m_func;
double m_t0;
double m_time; //!< The current integrator time

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

@ -223,6 +223,8 @@ public:
protected:
//! Pointer to the IDA memory for the problem
void* m_ida_mem;
void* m_linsol; //!< Sundials linear solver object
void* m_linsol_matrix; //!< matrix used by Sundials
//! Initial value of the time
doublereal m_t0;

59
src/numerics/BandMatrix.cpp
View file

@ -11,10 +11,18 @@
#include "cantera/numerics/ctlapack.h"
#else
#if CT_SUNDIALS_USE_LAPACK
#include "cvodes/cvodes_lapack.h"
#if CT_SUNDIALS_VERSION >= 30
#include "sunlinsol/sunlinsol_lapackband.h"
#else
#include "cvodes/cvodes_lapack.h"
#endif
#else
#include "cvodes/cvodes_dense.h"
#include "cvodes/cvodes_band.h"
#if CT_SUNDIALS_VERSION >= 30
#include "sunlinsol/sunlinsol_band.h"
#else
#include "cvodes/cvodes_dense.h"
#include "cvodes/cvodes_band.h"
#endif
#endif
#endif
@ -26,27 +34,46 @@ using namespace std;
namespace Cantera
{
// pImpl wrapper class for vector of Sundials index types to avoid needing to
// include Sundials headers in BandMatrix.h
struct BandMatrix::PivData {
#if CT_USE_LAPACK
vector_int data;
#elif CT_SUNDIALS_VERSION >= 30
std::vector<sunindextype> data;
#else
std::vector<long int> data;
#endif
};
BandMatrix::BandMatrix() :
m_n(0),
m_kl(0),
m_ku(0),
m_zero(0.0),
m_ipiv{new PivData()},
m_info(0)
{
}
BandMatrix::~BandMatrix()
{
// Needs to be defined here so m_ipiv can be deleted
}
BandMatrix::BandMatrix(size_t n, size_t kl, size_t ku, doublereal v) :
m_n(n),
m_kl(kl),
m_ku(ku),
m_zero(0.0),
m_ipiv{new PivData()},
m_info(0)
{
data.resize(n*(2*kl + ku + 1));
ludata.resize(n*(2*kl + ku + 1));
fill(data.begin(), data.end(), v);
fill(ludata.begin(), ludata.end(), 0.0);
m_ipiv.resize(m_n);
m_ipiv->data.resize(m_n);
m_colPtrs.resize(n);
m_lu_col_ptrs.resize(n);
size_t ldab = (2*kl + ku + 1);
@ -62,6 +89,7 @@ BandMatrix::BandMatrix(const BandMatrix& y) :
m_kl(0),
m_ku(0),
m_zero(0.0),
m_ipiv{new PivData()},
m_info(y.m_info)
{
m_n = y.m_n;
@ -69,7 +97,7 @@ BandMatrix::BandMatrix(const BandMatrix& y) :
m_ku = y.m_ku;
data = y.data;
ludata = y.ludata;
m_ipiv = y.m_ipiv;
m_ipiv->data = y.m_ipiv->data;
m_colPtrs.resize(m_n);
m_lu_col_ptrs.resize(m_n);
size_t ldab = (2 *m_kl + m_ku + 1);
@ -88,7 +116,7 @@ BandMatrix& BandMatrix::operator=(const BandMatrix& y)
m_n = y.m_n;
m_kl = y.m_kl;
m_ku = y.m_ku;
m_ipiv = y.m_ipiv;
m_ipiv->data = y.m_ipiv->data;
data = y.data;
ludata = y.ludata;
m_colPtrs.resize(m_n);
@ -109,7 +137,7 @@ void BandMatrix::resize(size_t n, size_t kl, size_t ku, doublereal v)
m_ku = ku;
data.resize(n*(2*kl + ku + 1));
ludata.resize(n*(2*kl + ku + 1));
m_ipiv.resize(m_n);
m_ipiv->data.resize(m_n);
fill(data.begin(), data.end(), v);
m_colPtrs.resize(m_n);
m_lu_col_ptrs.resize(m_n);
@ -195,11 +223,6 @@ size_t BandMatrix::ldim() const
return 2*m_kl + m_ku + 1;
}
pivot_vector_t& BandMatrix::ipiv()
{
return m_ipiv;
}
void BandMatrix::mult(const doublereal* b, doublereal* prod) const
{
for (size_t m = 0; m < m_n; m++) {
@ -231,13 +254,13 @@ int BandMatrix::factor()
ludata = data;
#if CT_USE_LAPACK
ct_dgbtrf(nRows(), nColumns(), nSubDiagonals(), nSuperDiagonals(),
ludata.data(), ldim(), ipiv().data(), m_info);
ludata.data(), ldim(), m_ipiv->data.data(), m_info);
#else
long int nu = static_cast<long int>(nSuperDiagonals());
long int nl = static_cast<long int>(nSubDiagonals());
long int smu = nu + nl;
m_info = bandGBTRF(m_lu_col_ptrs.data(), static_cast<long int>(nColumns()),
nu, nl, smu, m_ipiv.data());
nu, nl, smu, m_ipiv->data.data());
#endif
if (m_info != 0) {
throw Cantera::CanteraError("BandMatrix::factor",
@ -264,13 +287,13 @@ int BandMatrix::solve(doublereal* b, size_t nrhs, size_t ldb)
#if CT_USE_LAPACK
ct_dgbtrs(ctlapack::NoTranspose, nColumns(), nSubDiagonals(),
nSuperDiagonals(), nrhs, ludata.data(), ldim(),
ipiv().data(), b, ldb, m_info);
m_ipiv->data.data(), b, ldb, m_info);
#else
long int nu = static_cast<long int>(nSuperDiagonals());
long int nl = static_cast<long int>(nSubDiagonals());
long int smu = nu + nl;
double** a = m_lu_col_ptrs.data();
bandGBTRS(a, static_cast<long int>(nColumns()), smu, nl, m_ipiv.data(), b);
bandGBTRS(a, static_cast<long int>(nColumns()), smu, nl, m_ipiv->data.data(), b);
m_info = 0;
#endif
@ -327,8 +350,8 @@ doublereal BandMatrix::rcond(doublereal a1norm)
#if CT_USE_LAPACK
size_t ldab = (2 *m_kl + m_ku + 1);
int rinfo = 0;
double rcond = ct_dgbcon('1', m_n, m_kl, m_ku, ludata.data(), ldab, m_ipiv.data(), a1norm, work_.data(),
iwork_.data(), rinfo);
double rcond = ct_dgbcon('1', m_n, m_kl, m_ku, ludata.data(),
ldab, m_ipiv->data.data(), a1norm, work_.data(), iwork_.data(), rinfo);
if (rinfo != 0) {
throw CanteraError("BandMatrix::rcond()", "DGBCON returned INFO = {}", rinfo);
}

79
src/numerics/CVodesIntegrator.cpp
View file

@ -14,14 +14,28 @@ using namespace std;
#include "sundials/sundials_nvector.h"
#include "nvector/nvector_serial.h"
#include "cvodes/cvodes.h"
#if CT_SUNDIALS_USE_LAPACK
#include "cvodes/cvodes_lapack.h"
#if CT_SUNDIALS_VERSION >= 30
#if CT_SUNDIALS_USE_LAPACK
#include "sunlinsol/sunlinsol_lapackdense.h"
#include "sunlinsol/sunlinsol_lapackband.h"
#else
#include "sunlinsol/sunlinsol_dense.h"
#include "sunlinsol/sunlinsol_band.h"
#endif
#include "sunlinsol/sunlinsol_spgmr.h"
#include "cvodes/cvodes_direct.h"
#include "cvodes/cvodes_diag.h"
#include "cvodes/cvodes_spils.h"
#else
#include "cvodes/cvodes_dense.h"
#include "cvodes/cvodes_band.h"
#if CT_SUNDIALS_USE_LAPACK
#include "cvodes/cvodes_lapack.h"
#else
#include "cvodes/cvodes_dense.h"
#include "cvodes/cvodes_band.h"
#endif
#include "cvodes/cvodes_diag.h"
#include "cvodes/cvodes_spgmr.h"
#endif
#include "cvodes/cvodes_diag.h"
#include "cvodes/cvodes_spgmr.h"
#define CV_SS 1
#define CV_SV 2
@ -65,6 +79,8 @@ extern "C" {
CVodesIntegrator::CVodesIntegrator() :
m_neq(0),
m_cvode_mem(0),
m_linsol(0),
m_linsol_matrix(0),
m_func(0),
m_t0(0.0),
m_y(0),
@ -98,6 +114,12 @@ CVodesIntegrator::~CVodesIntegrator()
}
CVodeFree(&m_cvode_mem);
}
#if CT_SUNDIALS_VERSION >= 30
SUNLinSolFree((SUNLinearSolver) m_linsol);
SUNMatDestroy((SUNMatrix) m_linsol_matrix);
#endif
if (m_y) {
N_VDestroy_Serial(m_y);
}
@ -338,23 +360,54 @@ void CVodesIntegrator::applyOptions()
{
if (m_type == DENSE + NOJAC) {
sd_size_t N = static_cast<sd_size_t>(m_neq);
#if CT_SUNDIALS_USE_LAPACK
CVLapackDense(m_cvode_mem, N);
#if CT_SUNDIALS_VERSION >= 30
SUNLinSolFree((SUNLinearSolver) m_linsol);
SUNMatDestroy((SUNMatrix) m_linsol_matrix);
m_linsol_matrix = SUNDenseMatrix(N, N);
#if CT_SUNDIALS_USE_LAPACK
m_linsol = SUNLapackDense(m_y, (SUNMatrix) m_linsol_matrix);
#else
m_linsol = SUNDenseLinearSolver(m_y, (SUNMatrix) m_linsol_matrix);
#endif
CVDlsSetLinearSolver(m_cvode_mem, (SUNLinearSolver) m_linsol,
(SUNMatrix) m_linsol_matrix);
#else
CVDense(m_cvode_mem, N);
#if CT_SUNDIALS_USE_LAPACK
CVLapackDense(m_cvode_mem, N);
#else
CVDense(m_cvode_mem, N);
#endif
#endif
} else if (m_type == DIAG) {
CVDiag(m_cvode_mem);
} else if (m_type == GMRES) {
CVSpgmr(m_cvode_mem, PREC_NONE, 0);
#if CT_SUNDIALS_VERSION >= 30
m_linsol = SUNSPGMR(m_y, PREC_NONE, 0);
CVSpilsSetLinearSolver(m_cvode_mem, (SUNLinearSolver) m_linsol);
#else
CVSpgmr(m_cvode_mem, PREC_NONE, 0);
#endif
} else if (m_type == BAND + NOJAC) {
sd_size_t N = static_cast<sd_size_t>(m_neq);
long int nu = m_mupper;
long int nl = m_mlower;
#if CT_SUNDIALS_USE_LAPACK
CVLapackBand(m_cvode_mem, N, nu, nl);
#if CT_SUNDIALS_VERSION >= 30
SUNLinSolFree((SUNLinearSolver) m_linsol);
SUNMatDestroy((SUNMatrix) m_linsol_matrix);
m_linsol_matrix = SUNBandMatrix(N, nu, nl, nu+nl);
#if CT_SUNDIALS_USE_LAPACK
m_linsol = SUNLapackBand(m_y, (SUNMatrix) m_linsol_matrix);
#else
m_linsol = SUNBandLinearSolver(m_y, (SUNMatrix) m_linsol_matrix);
#endif
CVDlsSetLinearSolver(m_cvode_mem, (SUNLinearSolver) m_linsol,
(SUNMatrix) m_linsol_matrix);
#else
CVBand(m_cvode_mem, N, nu, nl);
#if CT_SUNDIALS_USE_LAPACK
CVLapackBand(m_cvode_mem, N, nu, nl);
#else
CVBand(m_cvode_mem, N, nu, nl);
#endif
#endif
} else {
throw CanteraError("CVodesIntegrator::applyOptions",

86
src/numerics/IDA_Solver.cpp
View file

@ -9,9 +9,22 @@
#include "sundials/sundials_types.h"
#include "sundials/sundials_math.h"
#include "ida/ida.h"
#include "ida/ida_dense.h"
#include "ida/ida_spgmr.h"
#include "ida/ida_band.h"
#if CT_SUNDIALS_VERSION >= 30
#if CT_SUNDIALS_USE_LAPACK
#include "sunlinsol/sunlinsol_lapackdense.h"
#include "sunlinsol/sunlinsol_lapackband.h"
#else
#include "sunlinsol/sunlinsol_dense.h"
#include "sunlinsol/sunlinsol_band.h"
#endif
#include "sunlinsol/sunlinsol_spgmr.h"
#include "ida/ida_direct.h"
#include "ida/ida_spils.h"
#else
#include "ida/ida_dense.h"
#include "ida/ida_spgmr.h"
#include "ida/ida_band.h"
#endif
#include "nvector/nvector_serial.h"
using namespace std;
@ -94,6 +107,28 @@ extern "C" {
* In the case of a recoverable error return, the integrator will attempt to
* recover by reducing the stepsize (which changes cj).
*/
#if CT_SUNDIALS_VERSION >= 30
static int ida_jacobian(realtype t, realtype c_j, N_Vector y, N_Vector yp,
N_Vector r, SUNMatrix Jac, void *f_data,
N_Vector tmp1, N_Vector tmp2, N_Vector tmp3)
{
Cantera::ResidData* d = (Cantera::ResidData*) f_data;
Cantera::ResidJacEval* f = d->m_func;
Cantera::IDA_Solver* s = d->m_solver;
double delta_t = s->getCurrentStepFromIDA();
double** cols;
if (SUNMatGetID(Jac) == SUNMATRIX_DENSE) {
cols = SM_COLS_D(Jac);
} else if (SUNMatGetID(Jac) == SUNMATRIX_BAND) {
cols = SM_COLS_B(Jac);
} else {
throw Cantera::CanteraError("ida_jacobian", "Unknown SUNMatrix type");
}
f->evalJacobianDP(t, delta_t, c_j, NV_DATA_S(y), NV_DATA_S(yp),
cols, NV_DATA_S(r));
return 0;
}
#else
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)
{
@ -105,6 +140,8 @@ extern "C" {
Jac->cols, NV_DATA_S(r));
return 0;
}
#endif
}
namespace Cantera
@ -265,7 +302,11 @@ void IDA_Solver::setJacobianType(int formJac)
{
m_formJac = formJac;
if (m_ida_mem && m_formJac == 1) {
int flag = IDADlsSetDenseJacFn(m_ida_mem, ida_jacobian);
#if CT_SUNDIALS_VERSION >= 30
int flag = IDADlsSetJacFn(m_ida_mem, ida_jacobian);
#else
int flag = IDADlsSetDenseJacFn(m_ida_mem, ida_jacobian);
#endif
if (flag != IDA_SUCCESS) {
throw CanteraError("IDA_Solver::setJacobianType",
"IDADlsSetDenseJacFn failed.");
@ -376,7 +417,21 @@ void IDA_Solver::init(doublereal t0)
// set the linear solver type
if (m_type == 1 || m_type == 0) {
long int N = m_neq;
flag = IDADense(m_ida_mem, N);
int flag;
#if CT_SUNDIALS_VERSION >= 30
SUNLinSolFree((SUNLinearSolver) m_linsol);
SUNMatDestroy((SUNMatrix) m_linsol_matrix);
m_linsol_matrix = SUNDenseMatrix(N, N);
#if CT_SUNDIALS_USE_LAPACK
m_linsol = SUNLapackDense(m_y, (SUNMatrix) m_linsol_matrix);
#else
m_linsol = SUNDenseLinearSolver(m_y, (SUNMatrix) m_linsol_matrix);
#endif
flag = IDADlsSetLinearSolver(m_ida_mem, (SUNLinearSolver) m_linsol,
(SUNMatrix) m_linsol_matrix);
#else
flag = IDADense(m_ida_mem, N);
#endif
if (flag) {
throw CanteraError("IDA_Solver::init", "IDADense failed");
}
@ -384,14 +439,31 @@ void IDA_Solver::init(doublereal t0)
long int N = m_neq;
long int nu = m_mupper;
long int nl = m_mlower;
IDABand(m_ida_mem, N, nu, nl);
#if CT_SUNDIALS_VERSION >= 30
SUNLinSolFree((SUNLinearSolver) m_linsol);
SUNMatDestroy((SUNMatrix) m_linsol_matrix);
m_linsol_matrix = SUNBandMatrix(N, nu, nl, nu+nl);
#if CT_SUNDIALS_USE_LAPACK
m_linsol = SUNLapackBand(m_y, (SUNMatrix) m_linsol_matrix);
#else
m_linsol = SUNBandLinearSolver(m_y, (SUNMatrix) m_linsol_matrix);
#endif
IDADlsSetLinearSolver(m_ida_mem, (SUNLinearSolver) m_linsol,
(SUNMatrix) m_linsol_matrix);
#else
IDABand(m_ida_mem, N, nu, nl);
#endif
} else {
throw CanteraError("IDA_Solver::init",
"unsupported linear solver type");
}
if (m_formJac == 1) {
flag = IDADlsSetDenseJacFn(m_ida_mem, ida_jacobian);
#if CT_SUNDIALS_VERSION >= 30
flag = IDADlsSetJacFn(m_ida_mem, ida_jacobian);
#else
flag = IDADlsSetDenseJacFn(m_ida_mem, ida_jacobian);
#endif
if (flag != IDA_SUCCESS) {
throw CanteraError("IDA_Solver::init",
"IDADlsSetDenseJacFn failed.");