From 9c50f7528c71ef470d73939f580d03ea342d85dd Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Sat, 10 Feb 2018 23:09:53 -0500 Subject: [PATCH] Add support for Sundials 3.1.0 (also 3.0.0) --- SConstruct | 17 +++- include/cantera/numerics/BandMatrix.h | 15 +--- include/cantera/numerics/CVodesIntegrator.h | 2 + include/cantera/numerics/IDA_Solver.h | 2 + src/numerics/BandMatrix.cpp | 59 +++++++++----- src/numerics/CVodesIntegrator.cpp | 79 +++++++++++++++---- src/numerics/IDA_Solver.cpp | 86 +++++++++++++++++++-- 7 files changed, 208 insertions(+), 52 deletions(-) diff --git a/SConstruct b/SConstruct index e07be4e8d..76fb3253f 100644 --- a/SConstruct +++ b/SConstruct @@ -938,6 +938,9 @@ def get_expression_value(includes, expression): s.extend(('#define Q(x) #x', '#define QUOTE(x) Q(x)', '#include ', + '#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']) diff --git a/include/cantera/numerics/BandMatrix.h b/include/cantera/numerics/BandMatrix.h index 0e7118049..012061f85 100644 --- a/include/cantera/numerics/BandMatrix.h +++ b/include/cantera/numerics/BandMatrix.h @@ -16,13 +16,6 @@ namespace Cantera { -#if CT_USE_LAPACK -typedef vector_int pivot_vector_t; -#else -typedef std::vector 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 m_ipiv; //! Vector of column pointers std::vector m_colPtrs; diff --git a/include/cantera/numerics/CVodesIntegrator.h b/include/cantera/numerics/CVodesIntegrator.h index 2479b8c4c..290054cc9 100644 --- a/include/cantera/numerics/CVodesIntegrator.h +++ b/include/cantera/numerics/CVodesIntegrator.h @@ -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 diff --git a/include/cantera/numerics/IDA_Solver.h b/include/cantera/numerics/IDA_Solver.h index 3b39db7e3..482e8e895 100644 --- a/include/cantera/numerics/IDA_Solver.h +++ b/include/cantera/numerics/IDA_Solver.h @@ -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; diff --git a/src/numerics/BandMatrix.cpp b/src/numerics/BandMatrix.cpp index 034b8e415..c9b85e553 100644 --- a/src/numerics/BandMatrix.cpp +++ b/src/numerics/BandMatrix.cpp @@ -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 data; +#else + std::vector 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(nSuperDiagonals()); long int nl = static_cast(nSubDiagonals()); long int smu = nu + nl; m_info = bandGBTRF(m_lu_col_ptrs.data(), static_cast(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(nSuperDiagonals()); long int nl = static_cast(nSubDiagonals()); long int smu = nu + nl; double** a = m_lu_col_ptrs.data(); - bandGBTRS(a, static_cast(nColumns()), smu, nl, m_ipiv.data(), b); + bandGBTRS(a, static_cast(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); } diff --git a/src/numerics/CVodesIntegrator.cpp b/src/numerics/CVodesIntegrator.cpp index 2a412d0a7..4ceea9e67 100644 --- a/src/numerics/CVodesIntegrator.cpp +++ b/src/numerics/CVodesIntegrator.cpp @@ -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(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(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", diff --git a/src/numerics/IDA_Solver.cpp b/src/numerics/IDA_Solver.cpp index 6993a88d1..5cd54dea6 100644 --- a/src/numerics/IDA_Solver.cpp +++ b/src/numerics/IDA_Solver.cpp @@ -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.");