Add support for Sundials 3.1.0 (also 3.0.0)

SConstruct

@@ -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'))
+    if env['use_lapack'] and LooseVersion(env['sundials_version']) >= LooseVersion('3.0'):
-    linkSharedLibs.extend(('sundials_cvodes', 'sundials_ida', 'sundials_nvecserial'))
+        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'])

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<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;

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

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;

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"
+        #if CT_SUNDIALS_VERSION >= 30
-        #include "cvodes/cvodes_band.h"
+            #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(),
+    double rcond = ct_dgbcon('1', m_n, m_kl, m_ku, ludata.data(),
-                      iwork_.data(), rinfo);
+        ldab, m_ipiv->data.data(), a1norm, work_.data(), iwork_.data(), rinfo);
     if (rinfo != 0) {
         throw CanteraError("BandMatrix::rcond()", "DGBCON returned INFO = {}", rinfo);
     }

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
+#if CT_SUNDIALS_VERSION >= 30
-    #include "cvodes/cvodes_lapack.h"
+    #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"
+    #if CT_SUNDIALS_USE_LAPACK
-    #include "cvodes/cvodes_band.h"
+        #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
+        #if CT_SUNDIALS_VERSION >= 30
-            CVLapackDense(m_cvode_mem, N);
+            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
+        #if CT_SUNDIALS_VERSION >= 30
-            CVLapackBand(m_cvode_mem, N, nu, nl);
+            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",

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"
+#if CT_SUNDIALS_VERSION >= 30
-#include "ida/ida_spgmr.h"
+    #if CT_SUNDIALS_USE_LAPACK
-#include "ida/ida_band.h"
+        #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.");