Fix for Sundials 6.6 behavior changes

.github/workflows/main.yml

@@ -468,6 +468,8 @@ jobs:
           fmt-ver: 9.1
         - sundials-ver: 6.4.1
           fmt-ver: 10
+        - sundials-ver: 6.6
+          fmt-ver: 10
       fail-fast: false
     steps:
       - uses: actions/checkout@v3

SConstruct

@@ -1512,7 +1512,7 @@ if env['system_sundials'] == 'y':
     if sundials_ver < parse_version("3.0") or sundials_ver >= parse_version("7.0"):
         logger.error(f"Sundials version {env['sundials_version']!r} is not supported.")
         sys.exit(1)
-    elif sundials_ver > parse_version("6.4.1"):
+    elif sundials_ver > parse_version("6.6.0"):
         logger.warning(f"Sundials version {env['sundials_version']!r} has not been tested.")
 
     logger.info(f"Using system installation of Sundials version {sundials_version!r}.")

include/cantera/kinetics/ImplicitSurfChem.h

@@ -150,7 +150,7 @@ class ImplicitSurfChem : public FuncEval
     // overloaded methods of class FuncEval
 
     //! Return the number of equations
-    virtual size_t neq() {
+    virtual size_t neq() const {
         return m_nv;
     }
 

include/cantera/numerics/CVodesIntegrator.h

@@ -90,15 +90,27 @@ class CVodesIntegrator : public Integrator
 private:
     void sensInit(double t0, FuncEval& func);
 
+    //! Check whether a CVODES method indicated an error. If so, throw an exception
+    //! containing the method name and the error code stashed by the cvodes_err() function.
+    void checkError(long flag, const string& ctMethod, const string& cvodesMethod) const;
+
     size_t m_neq = 0;
     void* m_cvode_mem = nullptr;
     SundialsContext m_sundials_ctx; //!< SUNContext object for Sundials>=6.0
     void* m_linsol = nullptr; //!< Sundials linear solver object
     void* m_linsol_matrix = nullptr; //!< matrix used by Sundials
     FuncEval* m_func = nullptr;
     double m_t0 = 0.0;
-    double m_time; //!< The current integrator time
+
+    //! The current system time, corresponding to #m_y
+    double m_time;
+
+    //! The latest time reached by the integrator. May be greater than #m_time.
+    double m_tInteg;
+
+    //! The system state at #m_time
     N_Vector m_y = nullptr;
+
     N_Vector m_abstol = nullptr;
     N_Vector m_dky = nullptr;
     string m_type = "DENSE";

include/cantera/numerics/FuncEval.h

@@ -150,10 +150,10 @@ class FuncEval
     }
 
     //! Number of equations.
-    virtual size_t neq()=0;
+    virtual size_t neq() const = 0;
 
     //! Number of sensitivity parameters.
-    virtual size_t nparams() {
+    virtual size_t nparams() const {
         return m_sens_params.size();
     }
 

include/cantera/numerics/IdasIntegrator.h

@@ -107,7 +107,11 @@ class IdasIntegrator : public Integrator
     FuncEval* m_func = nullptr;
 
     double m_t0 = 0.0; //!< The start time for the integrator
-    double m_time; //!< The current integrator time
+    double m_time; //!< The current integrator time, corresponding to #m_y
+
+    //! The latest time reached by the integrator. May be greater than #m_time
+    double m_tInteg;
+
     N_Vector m_y = nullptr; //!< The current system state
     N_Vector m_ydot = nullptr; //!< The time derivatives of the system state
     N_Vector m_abstol = nullptr; //!< Absolute tolerances for each state variable

include/cantera/zeroD/Reactor.h

@@ -157,7 +157,7 @@ class Reactor : public ReactorBase
 
     //! Number of sensitivity parameters associated with this reactor
     //! (including walls)
-    virtual size_t nSensParams();
+    virtual size_t nSensParams() const;
 
     //! Add a sensitivity parameter associated with the reaction number *rxn*
     //! (in the homogeneous phase).
@@ -183,14 +183,14 @@ class Reactor : public ReactorBase
 
     //! Check whether Reactor object uses advance limits
     //! @returns           True if at least one limit is set, False otherwise
-    bool hasAdvanceLimits() {
+    bool hasAdvanceLimits() const {
         return !m_advancelimits.empty();
     }
 
     //! Retrieve absolute step size limits during advance
     //! @param[out] limits array of step size limits with length neq
     //! @returns           True if at least one limit is set, False otherwise
-    bool getAdvanceLimits(double* limits);
+    bool getAdvanceLimits(double* limits) const;
 
     //! Set individual step size limit for component name *nm*
     //! @param nm component name

include/cantera/zeroD/ReactorNet.h

@@ -60,7 +60,7 @@ class ReactorNet : public FuncEval
     }
 
     //! Get the maximum integrator step.
-    double maxTimeStep() {
+    double maxTimeStep() const {
         return m_maxstep;
     }
 
@@ -205,11 +205,11 @@ class ReactorNet : public FuncEval
                       double* ydot, double* p, Array2D* j);
 
     // overloaded methods of class FuncEval
-    virtual size_t neq() {
+    virtual size_t neq() const {
         return m_nv;
     }
 
-    size_t nReactors() {
+    size_t nReactors() const {
         return m_reactors.size();
     }
 
@@ -227,7 +227,7 @@ class ReactorNet : public FuncEval
 
     virtual void getConstraints(double* constraints);
 
-    virtual size_t nparams() {
+    virtual size_t nparams() const {
         return m_sens_params.size();
     }
 
@@ -254,7 +254,7 @@ class ReactorNet : public FuncEval
                                         double scale);
 
     //! The name of the p-th sensitivity parameter added to this ReactorNet.
-    const std::string& sensitivityParameterName(size_t p) {
+    const string& sensitivityParameterName(size_t p) const {
         return m_paramNames.at(p);
     }
 
@@ -284,10 +284,10 @@ class ReactorNet : public FuncEval
     void setAdvanceLimits(const double* limits);
 
     //! Check whether ReactorNet object uses advance limits
-    bool hasAdvanceLimits();
+    bool hasAdvanceLimits() const;
 
     //! Retrieve absolute step size limits during advance
-    bool getAdvanceLimits(double* limits);
+    bool getAdvanceLimits(double* limits) const;
 
     virtual void preconditionerSetup(double t, double* y, double gamma);
 
@@ -302,7 +302,7 @@ class ReactorNet : public FuncEval
 
 protected:
     //! Check that preconditioning is supported by all reactors in the network
-    virtual void checkPreconditionerSupported();
+    virtual void checkPreconditionerSupported() const;
 
     //! Update the preconditioner based on the already computed jacobian values
     virtual void updatePreconditioner(double gamma);
@@ -315,7 +315,7 @@ class ReactorNet : public FuncEval
     //! Returns the order used for last solution step of the ODE integrator
     //! The function is intended for internal use by ReactorNet::advance
     //! and deliberately not exposed in external interfaces.
-    virtual int lastOrder();
+    virtual int lastOrder() const;
 
     std::vector<Reactor*> m_reactors;
     std::unique_ptr<Integrator> m_integ;

samples/cxx/custom/custom.cpp

@@ -113,7 +113,7 @@ class ReactorODEs : public FuncEval {
      * Number of equations in the ODE system.
      *   - overridden from FuncEval, called by the integrator during initialization.
      */
-    size_t neq() {
+    size_t neq() const {
         return m_nEqs;
     }