path: root/sundials6.patch

diff --git a/.github/workflows/build.yml b/.github/workflows/build.yml
index a394b94..ec7cbe4 100644
--- a/.github/workflows/build.yml
+++ b/.github/workflows/build.yml
@@ -10,7 +10,7 @@ jobs:
     runs-on: ubuntu-latest
     strategy:
       matrix:
-        sundials_version: [2.7.0, 3.2.0]
+        sundials_version: [2.7.0, 3.2.0, 4.1.0, 5.8.0, 6.7.0, 7.0.0]
     steps:
     - uses: actions/checkout@v3
     - name: Install system
diff --git a/src/lib/sundials_callbacks_ida_cvode.pxi b/src/lib/sundials_callbacks_ida_cvode.pxi
index 366c3e2..38b98d4 100644
--- a/src/lib/sundials_callbacks_ida_cvode.pxi
+++ b/src/lib/sundials_callbacks_ida_cvode.pxi
@@ -742,32 +742,64 @@ cdef int ida_jacv(realtype t, N_Vector yy, N_Vector yp, N_Vector rr, N_Vector vv
 
 # Error handling callback functions
 # =================================
+IF SUNDIALS_VERSION >= (7,0,0):
+    cdef extern from "sundials/sundials_context.h":
+        ctypedef _SUNContext * SUNContext
+        cdef struct _SUNContext:
+            pass
+        ctypedef int SUNErrCode
+    cdef void cv_err(int line, const char* func, const char* file, const char* msg, SUNErrCode error_code, void* problem_data, SUNContext sunctx):
+        """
+        This method overrides the default handling of error messages.
+        """
+        cdef ProblemData pData = <ProblemData>problem_data
+        
+        if error_code > 0 and pData.verbose > 0: #Warning
+            print '[CVode Warning]', msg
+        
+        if pData.verbose > 2: #Verbosity is greater than NORMAL, print warnings and errors
+            if error_code < 0: #Error
+                print '[CVode Error]', msg
+ELSE:
+    cdef void cv_err(int error_code, const char *module, const char *function, char *msg, void *problem_data):
+        """
+        This method overrides the default handling of error messages.
+        """
+        cdef ProblemData pData = <ProblemData>problem_data
+        
+        if error_code > 0 and pData.verbose > 0: #Warning
+            print '[CVode Warning]', msg
+        
+        if pData.verbose > 2: #Verbosity is greater than NORMAL, print warnings and errors
+            if error_code < 0: #Error
+                print '[CVode Error]', msg
 
-cdef void cv_err(int error_code, const char *module, const char *function, char *msg, void *problem_data):
-    """
-    This method overrides the default handling of error messages.
-    """
-    cdef ProblemData pData = <ProblemData>problem_data
-    
-    if error_code > 0 and pData.verbose > 0: #Warning
-        print '[CVode Warning]', msg
-    
-    if pData.verbose > 2: #Verbosity is greater than NORMAL, print warnings and errors
-        if error_code < 0: #Error
-            print '[CVode Error]', msg
-            
-cdef void ida_err(int error_code, const char *module, const char *function, char *msg, void *problem_data):
-    """
-    This method overrides the default handling of error messages.
-    """
-    cdef ProblemData pData = <ProblemData>problem_data
-    
-    if error_code > 0 and pData.verbose > 0: #Warning
-        print '[IDA Warning]', msg
-    
-    if pData.verbose > 2: #Verbosity is greater than NORMAL, print warnings and errors
-        if error_code < 0: #Error
-            print '[IDA Error]', msg
+IF SUNDIALS_VERSION >= (7,0,0):
+    cdef void ida_err(int line, const char* func, const char* file, const char* msg, SUNErrCode error_code, void* problem_data, SUNContext sunctx):
+        """
+        This method overrides the default handling of error messages.
+        """
+        cdef ProblemData pData = <ProblemData>problem_data
+        
+        if error_code > 0 and pData.verbose > 0: #Warning
+            print '[IDA Warning]', msg
+        
+        if pData.verbose > 2: #Verbosity is greater than NORMAL, print warnings and errors
+            if error_code < 0: #Error
+                print '[IDA Error]', msg
+ELSE:
+    cdef void ida_err(int error_code, const char *module, const char *function, char *msg, void *problem_data):
+        """
+        This method overrides the default handling of error messages.
+        """
+        cdef ProblemData pData = <ProblemData>problem_data
+        
+        if error_code > 0 and pData.verbose > 0: #Warning
+            print '[IDA Warning]', msg
+        
+        if pData.verbose > 2: #Verbosity is greater than NORMAL, print warnings and errors
+            if error_code < 0: #Error
+                print '[IDA Error]', msg
 
 
 cdef class ProblemData:
diff --git a/src/lib/sundials_callbacks_kinsol.pxi b/src/lib/sundials_callbacks_kinsol.pxi
index fc73035..513077e 100644
--- a/src/lib/sundials_callbacks_kinsol.pxi
+++ b/src/lib/sundials_callbacks_kinsol.pxi
@@ -65,8 +65,13 @@ ELSE:
             return KINDLS_SUCCESS
         except Exception:
             return KINDLS_JACFUNC_RECVR #Recoverable Error (See Sundials description)
-            
-cdef int kin_jacv(N_Vector vv, N_Vector Jv, N_Vector vx, int* new_u,
+
+IF SUNDIALS_VERSION >= (6,0,0):
+    ctypedef bint kin_jacv_bool
+ELSE:
+    ctypedef int kin_jacv_bool
+
+cdef int kin_jacv(N_Vector vv, N_Vector Jv, N_Vector vx, kin_jacv_bool* new_u,
             void *problem_data):
     cdef ProblemDataEquationSolver pData = <ProblemDataEquationSolver>problem_data
     cdef N.ndarray x  = nv2arr(vx)
@@ -217,21 +222,40 @@ ELSE:
             return KIN_SYSFUNC_FAIL
         
         return KIN_SUCCESS
-        
 
-cdef void kin_err(int err_code, const char *module, const char *function, char *msg, void *eh_data):
-    cdef ProblemDataEquationSolver pData = <ProblemDataEquationSolver>eh_data
-    
-    if err_code > 0: #Warning
-        category = 1
-    elif err_code < 0: #Error
-        category = -1
-    else:
-        category = 0
-    
-    print "Error occured in <function: %s>."%function
-    print "<message: %s>"%msg
-    #print "<functionNorm: %g, scaledStepLength: %g, tolerance: %g>"%(fnorm, snorm, pData.TOL)
+IF SUNDIALS_VERSION >= (7,0,0):
+    cdef extern from "sundials/sundials_context.h":
+        ctypedef _SUNContext * SUNContext
+        cdef struct _SUNContext:
+            pass
+        ctypedef int SUNErrCode
+    cdef void kin_err(int line, const char* function, const char* file, const char* msg, SUNErrCode err_code, void* eh_data, SUNContext sunctx):
+        cdef ProblemDataEquationSolver pData = <ProblemDataEquationSolver>eh_data
+        
+        if err_code > 0: #Warning
+            category = 1
+        elif err_code < 0: #Error
+            category = -1
+        else:
+            category = 0
+        
+        print "Error occured in <function: %s>."%function
+        print "<message: %s>"%msg
+        #print "<functionNorm: %g, scaledStepLength: %g, tolerance: %g>"%(fnorm, snorm, pData.TOL)
+ELSE:
+    cdef void kin_err(int err_code, const char *module, const char *function, char *msg, void *eh_data):
+        cdef ProblemDataEquationSolver pData = <ProblemDataEquationSolver>eh_data
+        
+        if err_code > 0: #Warning
+            category = 1
+        elif err_code < 0: #Error
+            category = -1
+        else:
+            category = 0
+        
+        print "Error occured in <function: %s>."%function
+        print "<message: %s>"%msg
+        #print "<functionNorm: %g, scaledStepLength: %g, tolerance: %g>"%(fnorm, snorm, pData.TOL)
 
 
 cdef void kin_info(const char *module, const char *function, char *msg, void *eh_data):
diff --git a/src/lib/sundials_includes.pxd b/src/lib/sundials_includes.pxd
index f1e43ac..ea4c7c1 100644
--- a/src/lib/sundials_includes.pxd
+++ b/src/lib/sundials_includes.pxd
@@ -40,9 +40,22 @@ IF SUNDIALS_VERSION >= (6,0,0):
             pass
         int SUNContext_Create(void* comm, SUNContext* ctx)
 
-cdef extern from "sundials/sundials_types.h":
+IF SUNDIALS_VERSION >= (7,0,0):
+    cdef extern from "sundials/sundials_context.h":
+        ctypedef int SUNErrCode
+        ctypedef void (*SUNErrHandlerFn)(int line, const char* func, const char* file, const char* msg, SUNErrCode err_code, void* err_user_data, SUNContext sunctx)
+        SUNErrCode SUNContext_PushErrHandler(SUNContext sunctx, SUNErrHandlerFn err_fn, void* err_user_data)
+
+IF SUNDIALS_VERSION >= (6,0,0):
+    cdef extern from "sundials/sundials_types.h":
+        ctypedef double sunrealtype
+        ctypedef bint sunbooleantype
     ctypedef double realtype
-    ctypedef bint booleantype # should be bool instead of bint, but there is a bug in Cython
+    ctypedef bint booleantype
+ELSE:
+    cdef extern from "sundials/sundials_types.h":
+        ctypedef double realtype
+        ctypedef bint booleantype # should be bool instead of bint, but there is a bug in Cython
 
 cdef extern from "sundials/sundials_nvector.h":
     ctypedef _generic_N_Vector* N_Vector
@@ -333,12 +346,7 @@ cdef extern from "cvodes/cvodes.h":
     
     #Functions for retrieving results
     int CVodeGetDky(void *cvode_mem, realtype t, int k, N_Vector dky)
-    
-    #Functions for error handling
-    ctypedef void (*CVErrHandlerFn)(int error_code, const char *module, const char *function, char *msg,
-                                   void *eh_data)
-    int CVodeSetErrHandlerFn(void *cvode_mem, CVErrHandlerFn ehfun, void* eh_data)
-    
+
     #Functions for discontinuity handling
     ctypedef int (*CVRootFn)(realtype tt, N_Vector yy, realtype *gout, void *user_data)
     int CVodeRootDirection(void *cvode_mem, int *rootdir)
@@ -408,11 +416,31 @@ cdef extern from "cvodes/cvodes.h":
     int CVodeGetStgrSensNumNonlinSolvIters(void *cvode_mem, long int *nSTGR1niters)
     int CVodeGetStgrSensNumNonlinSolvConvFails(void *cvode_mem, long int *nSTGR1ncfails)
 
-cdef extern from "cvodes/cvodes_spils.h":
-    ctypedef int (*CVSpilsJacTimesVecFn)(N_Vector v, N_Vector Jv, realtype t,
-            N_Vector y, N_Vector fy, void *user_data, N_Vector tmp)
+IF SUNDIALS_VERSION >= (6,0,0):
+    cdef extern from "cvodes/cvodes_ls.h":
+        ctypedef int (*CVLsJacFn)(sunrealtype t, N_Vector y, N_Vector fy, SUNMatrix Jac, void* user_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3);
+        ctypedef int (*CVLsPrecSolveFn)(sunrealtype t, N_Vector y, N_Vector fy, N_Vector r, N_Vector z, sunrealtype gamma, sunrealtype delta, int lr, void* user_data);
+        ctypedef int (*CVLsJacTimesSetupFn)(realtype t, N_Vector y, N_Vector fy, void *user_data)
+        ctypedef int (*CVLsJacTimesVecFn)(N_Vector v, N_Vector Jv, realtype t, N_Vector y, N_Vector fy, void *user_data, N_Vector tmp)
+ELSE:
+    cdef extern from "cvodes/cvodes_spils.h":
+        ctypedef int (*CVSpilsJacTimesVecFn)(N_Vector v, N_Vector Jv, realtype t, N_Vector y, N_Vector fy, void *user_data, N_Vector tmp)
 
-IF SUNDIALS_VERSION >= (3,0,0):
+
+IF SUNDIALS_VERSION >= (6,0,0):
+    cdef extern from "cvode/cvode_ls.h":
+       int CVodeSetJacTimes(void *cvode_mem, CVLsJacTimesSetupFn jtsetup, CVLsJacTimesVecFn jtimes)
+       int CVodeSetLinearSolver(void* cvode_mem, SUNLinearSolver LS, SUNMatrix A)
+       int CVodeSetJacFn(void* cvode_mem, CVLsJacFn jac)
+    IF SUNDIALS_WITH_SUPERLU:
+        cdef extern from "sunlinsol/sunlinsol_superlumt.h":
+            SUNLinearSolver SUNLinSol_SuperLUMT(N_Vector y, SUNMatrix A, int num_threads, SUNContext ctx)
+    ELSE:
+        cdef inline SUNLinearSolver SUNLinSol_SuperLUMT(N_Vector y, SUNMatrix A, int num_threads, SUNContext ctx): return NULL
+
+    cdef inline int cv_spils_jtsetup_dummy(realtype t, N_Vector y, N_Vector fy, void *user_data): return 0
+    cdef inline tuple version(): return (6,0,0)
+ELIF SUNDIALS_VERSION >= (3,0,0):
     cdef extern from "cvodes/cvodes_direct.h":
         ctypedef int (*CVDlsDenseJacFn)(realtype t, N_Vector y, N_Vector fy, 
                        SUNMatrix Jac, void *user_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3)
@@ -435,21 +463,13 @@ IF SUNDIALS_VERSION >= (3,0,0):
 				  N_Vector r, N_Vector z,
 				  realtype gamma, realtype delta, int lr, void *user_data)
 
-
     IF SUNDIALS_WITH_SUPERLU:
-        IF SUNDIALS_VERSION >= (6,0,0):
-            cdef extern from "sunlinsol/sunlinsol_superlumt.h":
-                SUNLinearSolver SUNLinSol_SuperLUMT(N_Vector y, SUNMatrix A, int num_threads, SUNContext ctx)
-        ELSE:
-            cdef extern from "sunlinsol/sunlinsol_superlumt.h":
-                SUNLinearSolver SUNSuperLUMT(N_Vector y, SUNMatrix A, int num_threads)
+        cdef extern from "sunlinsol/sunlinsol_superlumt.h":
+            SUNLinearSolver SUNSuperLUMT(N_Vector y, SUNMatrix A, int num_threads)
     ELSE:
-        IF SUNDIALS_VERSION >= (6,0,0):
-            cdef inline SUNLinearSolver SUNLinSol_SuperLUMT(N_Vector y, SUNMatrix A, int num_threads, SUNContext ctx): return NULL
-        ELSE:
-            cdef inline SUNLinearSolver SUNSuperLUMT(N_Vector y, SUNMatrix A, int num_threads): return NULL
+         cdef inline SUNLinearSolver SUNSuperLUMT(N_Vector y, SUNMatrix A, int num_threads): return NULL
     
-    cdef inline int cv_spils_jtsetup_dummy(realtype t, N_Vector y, N_Vector fy, void *user_data): return 0    
+    cdef inline int cv_spils_jtsetup_dummy(realtype t, N_Vector y, N_Vector fy, void *user_data): return 0
     cdef inline tuple version(): return (3,0,0)
 ELSE:
     cdef extern from "cvodes/cvodes_dense.h":
@@ -494,20 +514,29 @@ ELSE:
         cdef inline int CVSlsSetSparseJacFn(void *cvode_mem, CVSlsSparseJacFn jac): return -1
         cdef inline int CVSlsGetNumJacEvals(void *cvode_mem, long int *njevals): return -1
         cdef inline tuple version(): return (2,5,0)
-    
-cdef extern from "cvodes/cvodes_spils.h":
-    IF SUNDIALS_VERSION >= (4,0,0):
-        int CVodeSetPreconditioner(void *cvode_mem, CVSpilsPrecSetupFn psetup, CVSpilsPrecSolveFn psolve)
+
+IF SUNDIALS_VERSION >= (6,0,0):
+    cdef extern from "cvode/cvode_ls.h":
+        ctypedef int (*CVLsPrecSetupFn)(sunrealtype t, N_Vector y, N_Vector fy, sunbooleantype jok, sunbooleantype* jcurPtr, sunrealtype gamma, void* user_data)
+        int CVodeSetPreconditioner(void* cvode_mem, CVLsPrecSetupFn pset, CVLsPrecSolveFn psolve)
         int CVodeGetNumJtimesEvals(void *cvode_mem, long int *njvevals) #Number of jac*vector evals
         int CVodeGetNumRhsEvals(void *cvode_mem, long int *nfevalsLS) #Number of res evals due to jacÄvector evals
         int CVodeGetNumPrecEvals(void *cvode_mem, long int *npevals)
         int CVodeGetNumPrecSolves(void *cvode_mem, long int *npsolves)
-    ELSE:
-        int CVSpilsSetPreconditioner(void *cvode_mem, CVSpilsPrecSetupFn psetup, CVSpilsPrecSolveFn psolve)
-        int CVSpilsGetNumJtimesEvals(void *cvode_mem, long int *njvevals) #Number of jac*vector evals
-        int CVSpilsGetNumRhsEvals(void *cvode_mem, long int *nfevalsLS) #Number of res evals due to jacÄvector evals
-        int CVSpilsGetNumPrecEvals(void *cvode_mem, long int *npevals)
-        int CVSpilsGetNumPrecSolves(void *cvode_mem, long int *npsolves)
+ELSE:
+    cdef extern from "cvodes/cvodes_spils.h":
+        IF SUNDIALS_VERSION >= (4,0,0):
+            int CVodeSetPreconditioner(void *cvode_mem, CVSpilsPrecSetupFn psetup, CVSpilsPrecSolveFn psolve)
+            int CVodeGetNumJtimesEvals(void *cvode_mem, long int *njvevals) #Number of jac*vector evals
+            int CVodeGetNumRhsEvals(void *cvode_mem, long int *nfevalsLS) #Number of res evals due to jacÄvector evals
+            int CVodeGetNumPrecEvals(void *cvode_mem, long int *npevals)
+            int CVodeGetNumPrecSolves(void *cvode_mem, long int *npsolves)
+        ELSE:
+            int CVSpilsSetPreconditioner(void *cvode_mem, CVSpilsPrecSetupFn psetup, CVSpilsPrecSolveFn psolve)
+            int CVSpilsGetNumJtimesEvals(void *cvode_mem, long int *njvevals) #Number of jac*vector evals
+            int CVSpilsGetNumRhsEvals(void *cvode_mem, long int *nfevalsLS) #Number of res evals due to jacÄvector evals
+            int CVSpilsGetNumPrecEvals(void *cvode_mem, long int *npevals)
+            int CVSpilsGetNumPrecSolves(void *cvode_mem, long int *npsolves)
 
 cdef extern from "idas/idas.h":
     ctypedef int (*IDAResFn)(realtype tt, N_Vector yy, N_Vector yp, N_Vector rr, void *user_data)
@@ -536,13 +565,7 @@ cdef extern from "idas/idas.h":
     
     #Functions for retrieving results
     int IDAGetDky(void *ida_mem, realtype t, int k, N_Vector dky)
-    
-    #Functions for error handling
-    ctypedef void (*IDAErrHandlerFn)(int error_code, const char *module, const char *function, char *msg,
-                                    void *eh_data)
-    int IDASetErrHandlerFn(void *ida_mem,IDAErrHandlerFn ehfun, void* eh_data)
-    
-    
+
     #Functions for discontinuity handling
     ctypedef int (*IDARootFn)(realtype tt, N_Vector yy, N_Vector yp, realtype *gout, void *user_data)
     int IDASetRootDirection(void *ida_mem, int *rootdir)
@@ -615,12 +638,25 @@ cdef extern from "idas/idas.h":
     
     #End Sensitivities
     #=================
-    
-cdef extern from "idas/idas_spils.h":
-    ctypedef int (*IDASpilsJacTimesVecFn)(realtype tt, N_Vector yy, N_Vector yp, N_Vector rr, 
-            N_Vector v, N_Vector Jv, realtype cj, void *user_data,N_Vector tmp1, N_Vector tmp2)
-    
-IF SUNDIALS_VERSION >= (3,0,0):
+
+IF SUNDIALS_VERSION >= (6,0,0):
+    cdef extern from "ida/ida_ls.h":
+        ctypedef int (*IDALsJacFn)(sunrealtype t, sunrealtype c_j, N_Vector y, N_Vector yp, N_Vector r, SUNMatrix Jac, void* user_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3)
+        ctypedef int (*IDALsJacTimesSetupFn)(sunrealtype tt, N_Vector yy, N_Vector yp, N_Vector rr, sunrealtype c_j, void* user_data)
+        ctypedef int (*IDALsJacTimesVecFn)(sunrealtype tt, N_Vector yy, N_Vector yp, N_Vector rr, N_Vector v, N_Vector Jv, sunrealtype c_j, void* user_data, N_Vector tmp1, N_Vector tmp2)
+ELSE:
+    cdef extern from "idas/idas_spils.h":
+        ctypedef int (*IDASpilsJacTimesVecFn)(realtype tt, N_Vector yy, N_Vector yp, N_Vector rr, 
+                N_Vector v, N_Vector Jv, realtype cj, void *user_data,N_Vector tmp1, N_Vector tmp2)
+
+IF SUNDIALS_VERSION >= (6,0,0):
+    cdef extern from "ida/ida_ls.h":
+        int IDASetJacFn(void* ida_mem, IDALsJacFn jac)
+        int IDASetLinearSolver(void* ida_mem, SUNLinearSolver LS, SUNMatrix A)
+        int IDASetJacTimes(void* ida_mem, IDALsJacTimesSetupFn jtsetup, IDALsJacTimesVecFn jtimes)
+
+    cdef inline int ida_spils_jtsetup_dummy(realtype tt, N_Vector yy, N_Vector yp, N_Vector rr, realtype c_j, void *user_data): return 0
+ELIF SUNDIALS_VERSION >= (3,0,0):
     cdef extern from "idas/idas_direct.h":
         ctypedef int (*IDADlsDenseJacFn)(realtype tt, realtype cj, N_Vector yy, 
                        N_Vector yp, N_Vector rr, SUNMatrix Jac, void *user_data, 
@@ -660,13 +696,18 @@ ELSE:
     cdef extern from "idas/idas_spils.h":
         int IDASpilsSetJacTimesVecFn(void *ida_mem, IDASpilsJacTimesVecFn ida_jacv)
 
-cdef extern from "idas/idas_spils.h":
-    IF SUNDIALS_VERSION >= (4,0,0):
+IF SUNDIALS_VERSION >= (6,0,0):
+    cdef extern from "ida/ida_ls.h": 
         int IDAGetNumJtimesEvals(void *ida_mem, long int *njvevals) #Number of jac*vector
         int IDAGetNumResEvals(void *ida_mem, long int *nfevalsLS) #Number of rhs due to jac*vector
-    ELSE:
-        int IDASpilsGetNumJtimesEvals(void *ida_mem, long int *njvevals) #Number of jac*vector
-        int IDASpilsGetNumResEvals(void *ida_mem, long int *nfevalsLS) #Number of rhs due to jac*vector
+ELSE:
+    cdef extern from "idas/idas_spils.h":
+        IF SUNDIALS_VERSION >= (4,0,0):
+            int IDAGetNumJtimesEvals(void *ida_mem, long int *njvevals) #Number of jac*vector
+            int IDAGetNumResEvals(void *ida_mem, long int *nfevalsLS) #Number of rhs due to jac*vector
+        ELSE:
+            int IDASpilsGetNumJtimesEvals(void *ida_mem, long int *njvevals) #Number of jac*vector
+            int IDASpilsGetNumResEvals(void *ida_mem, long int *nfevalsLS) #Number of rhs due to jac*vector
 
 
 ####################
@@ -678,7 +719,6 @@ cdef extern from "idas/idas_spils.h":
 cdef extern from "kinsol/kinsol.h":
     # user defined functions
     ctypedef int (*KINSysFn)(N_Vector uu, N_Vector fval, void *user_data )
-    ctypedef void (*KINErrHandlerFn)(int error_code, char *module, char *function, char *msg, void *user_data)
     ctypedef void (*KINInfoHandlerFn)(const char *module, const char *function, char *msg, void *user_data)
     # initialization routines
     IF SUNDIALS_VERSION >= (6,0,0):
@@ -689,7 +729,6 @@ cdef extern from "kinsol/kinsol.h":
 
     # optional input spec. functions,
     # for specificationsdocumentation cf. kinsol.h line 218-449
-    int KINSetErrHandlerFn(void *kinmem, KINErrHandlerFn ehfun, void *eh_data)
     int KINSetInfoHandlerFn(void *kinmem, KINInfoHandlerFn ihfun, void *ih_data)
     int KINSetUserData(void *kinmem, void *user_data)
     int KINSetPrintLevel(void *kinmemm, int printfl)
@@ -729,8 +768,24 @@ cdef extern from "kinsol/kinsol.h":
     # fuction used to deallocate memory used by KINSOL
     void KINFree(void **kinmem)
 
+# Functions for error handling
+IF SUNDIALS_VERSION < (7,0,0):
+    cdef extern from "kinsol/kinsol.h":
+        ctypedef void (*KINErrHandlerFn)(int error_code, char *module, char *function, char *msg, void *user_data)
+        int KINSetErrHandlerFn(void *kinmem, KINErrHandlerFn ehfun, void *eh_data)
+    cdef extern from "cvodes/cvodes.h":
+        ctypedef void (*CVErrHandlerFn)(int error_code, const char *module, const char *function, char *msg, void *eh_data)
+        int CVodeSetErrHandlerFn(void *cvode_mem, CVErrHandlerFn ehfun, void* eh_data)
+    cdef extern from "idas/idas.h":
+        ctypedef void (*IDAErrHandlerFn)(int error_code, const char *module, const char *function, char *msg, void *eh_data) 
+        int IDASetErrHandlerFn(void *ida_mem,IDAErrHandlerFn ehfun, void* eh_data)
 
-IF SUNDIALS_VERSION >= (3,0,0):
+IF SUNDIALS_VERSION >= (6,0,0):
+    cdef extern from "kinsol/kinsol_ls.h":
+        ctypedef int (*KINLsJacFn)(N_Vector u, N_Vector fu, SUNMatrix J, void* user_data, N_Vector tmp1, N_Vector tmp2)
+        int KINSetLinearSolver(void* kinmem, SUNLinearSolver LS, SUNMatrix A)
+        int KINSetJacFn(void* kinmem, KINLsJacFn jac)
+ELIF SUNDIALS_VERSION >= (3,0,0):
     cdef extern from "kinsol/kinsol_direct.h":
         ctypedef int (*KINDlsDenseJacFn)(N_Vector u, N_Vector fu, SUNMatrix J, void *user_data, N_Vector tmp1, N_Vector tmp2)
         IF SUNDIALS_VERSION < (4,0,0):
@@ -767,40 +822,57 @@ ELSE:
         ctypedef int (*KINSpilsPrecSetupFn)(N_Vector u, N_Vector uscale,
                     N_Vector fval, N_Vector fscale, void *problem_data, N_Vector tmp1, N_Vector tmp2)
 
-cdef extern from "kinsol/kinsol_direct.h":
-    # optional output fcts for linear direct solver
-    int KINDlsGetWorkSpace(void *kinmem, long int *lenrwB, long int *leniwB)
-    IF SUNDIALS_VERSION >= (4,0,0):
+IF SUNDIALS_VERSION >= (6,0,0):
+    cdef extern from "kinsol/kinsol_ls.h":
+        ctypedef int (*KINLsPrecSetupFn)(N_Vector uu, N_Vector uscale, N_Vector fval, N_Vector fscale, void* user_data);
+        ctypedef int (*KINLsPrecSolveFn)(N_Vector uu, N_Vector uscale, N_Vector fval, N_Vector fscale, N_Vector vv, void* user_data)
+        ctypedef int (*KINLsJacTimesVecFn)(N_Vector v, N_Vector Jv, N_Vector uu, sunbooleantype* new_uu, void* J_data)
         int KINGetLastLinFlag(void *kinmem, long int *flag)
         int KINGetNumJacEvals(void *kinmem, long int *njevalsB)
         int KINGetNumFuncEvals(void *kinmem, long int *nfevalsB)
-    ELSE:
-        int KINDlsGetLastFlag(void *kinmem, long int *flag)
-        int KINDlsGetNumJacEvals(void *kinmem, long int *njevalsB)
-        int KINDlsGetNumFuncEvals(void *kinmem, long int *nfevalsB)
-    char *KINDlsGetReturnFlagName(int flag)
-
-cdef extern from "kinsol/kinsol_spils.h":
-    ctypedef int (*KINSpilsJacTimesVecFn)(N_Vector vv, N_Vector Jv, N_Vector vx, int* new_u,
-                void *problem_data)
-    IF SUNDIALS_VERSION >= (4,0,0):
-        int KINSetJacTimesVecFn(void *kinmem, KINSpilsJacTimesVecFn jacv)
-        int KINSetPreconditioner(void *kinmem, KINSpilsPrecSetupFn psetup, KINSpilsPrecSolveFn psolve)
+        int KINSetJacTimesVecFn(void* kinmem, KINLsJacTimesVecFn jtv)
+        int KINSetPreconditioner(void* kinmem, KINLsPrecSetupFn psetup, KINLsPrecSolveFn psolve)
         int KINGetNumLinIters(void *kinmem, long int *nliters)
         int KINGetNumLinConvFails(void *kinmem, long int *nlcfails)
         int KINGetNumPrecEvals(void *kinmem, long int *npevals)
         int KINGetNumPrecSolves(void *kinmem, long int *npsolves)
         int KINGetNumJtimesEvals(void *kinmem, long int *njevals)
         int KINGetNumFuncEvals(void *kinmem, long int *nfevalsLS)
-    ELSE:
-        int KINSpilsSetJacTimesVecFn(void *kinmem, KINSpilsJacTimesVecFn jacv)
-        int KINSpilsSetPreconditioner(void *kinmem, KINSpilsPrecSetupFn psetup, KINSpilsPrecSolveFn psolve)
-        int KINSpilsGetNumLinIters(void *kinmem, long int *nliters)
-        int KINSpilsGetNumConvFails(void *kinmem, long int *nlcfails)
-        int KINSpilsGetNumPrecEvals(void *kinmem, long int *npevals)
-        int KINSpilsGetNumPrecSolves(void *kinmem, long int *npsolves)
-        int KINSpilsGetNumJtimesEvals(void *kinmem, long int *njevals)
-        int KINSpilsGetNumFuncEvals(void *kinmem, long int *nfevalsLS)
+ELSE:
+    cdef extern from "kinsol/kinsol_direct.h":
+        # optional output fcts for linear direct solver
+        int KINDlsGetWorkSpace(void *kinmem, long int *lenrwB, long int *leniwB)
+        IF SUNDIALS_VERSION >= (4,0,0):
+            int KINGetLastLinFlag(void *kinmem, long int *flag)
+            int KINGetNumJacEvals(void *kinmem, long int *njevalsB)
+            int KINGetNumFuncEvals(void *kinmem, long int *nfevalsB)
+        ELSE:
+            int KINDlsGetLastFlag(void *kinmem, long int *flag)
+            int KINDlsGetNumJacEvals(void *kinmem, long int *njevalsB)
+            int KINDlsGetNumFuncEvals(void *kinmem, long int *nfevalsB)
+        char *KINDlsGetReturnFlagName(int flag)
+
+    cdef extern from "kinsol/kinsol_spils.h":
+        ctypedef int (*KINSpilsJacTimesVecFn)(N_Vector vv, N_Vector Jv, N_Vector vx, int* new_u,
+                    void *problem_data)
+        IF SUNDIALS_VERSION >= (4,0,0):
+            int KINSetJacTimesVecFn(void *kinmem, KINSpilsJacTimesVecFn jacv)
+            int KINSetPreconditioner(void *kinmem, KINSpilsPrecSetupFn psetup, KINSpilsPrecSolveFn psolve)
+            int KINGetNumLinIters(void *kinmem, long int *nliters)
+            int KINGetNumLinConvFails(void *kinmem, long int *nlcfails)
+            int KINGetNumPrecEvals(void *kinmem, long int *npevals)
+            int KINGetNumPrecSolves(void *kinmem, long int *npsolves)
+            int KINGetNumJtimesEvals(void *kinmem, long int *njevals)
+            int KINGetNumFuncEvals(void *kinmem, long int *nfevalsLS)
+        ELSE:
+            int KINSpilsSetJacTimesVecFn(void *kinmem, KINSpilsJacTimesVecFn jacv)
+            int KINSpilsSetPreconditioner(void *kinmem, KINSpilsPrecSetupFn psetup, KINSpilsPrecSolveFn psolve)
+            int KINSpilsGetNumLinIters(void *kinmem, long int *nliters)
+            int KINSpilsGetNumConvFails(void *kinmem, long int *nlcfails)
+            int KINSpilsGetNumPrecEvals(void *kinmem, long int *npevals)
+            int KINSpilsGetNumPrecSolves(void *kinmem, long int *npsolves)
+            int KINSpilsGetNumJtimesEvals(void *kinmem, long int *njevals)
+            int KINSpilsGetNumFuncEvals(void *kinmem, long int *nfevalsLS)
 
 #=========================
 # END SUNDIALS DEFINITIONS
diff --git a/src/solvers/kinsol.pyx b/src/solvers/kinsol.pyx
index c0dbaac..32ed0f0 100644
--- a/src/solvers/kinsol.pyx
+++ b/src/solvers/kinsol.pyx
@@ -210,7 +210,10 @@ cdef class KINSOL(Algebraic):
                 raise KINSOLError(flag)
             
             #Specify the error handling
-            flag = SUNDIALS.KINSetErrHandlerFn(self.kinsol_mem, kin_err, <void*>self.pData)
+            IF SUNDIALS_VERSION >= (7,0,0):
+                flag = SUNDIALS.SUNContext_PushErrHandler(ctx, kin_err, <void*>self.pData)
+            ELSE:
+                flag = SUNDIALS.KINSetErrHandlerFn(self.kinsol_mem, kin_err, <void*>self.pData)
             if flag < 0:
                 raise KINSOLError(flag) 
                 
diff --git a/src/solvers/sundials.pyx b/src/solvers/sundials.pyx
index 85aa459..20b96be 100644
--- a/src/solvers/sundials.pyx
+++ b/src/solvers/sundials.pyx
@@ -342,7 +342,10 @@ cdef class IDA(Implicit_ODE):
                     raise IDAError(flag,self.t)
             
             #Specify the error handling
-            flag = SUNDIALS.IDASetErrHandlerFn(self.ida_mem, ida_err, <void*>self.pData)
+            IF SUNDIALS_VERSION >= (7,0,0):
+                flag = SUNDIALS.SUNContext_PushErrHandler(ctx, ida_err, <void*>self.pData)
+            ELSE:
+                flag = SUNDIALS.IDASetErrHandlerFn(self.ida_mem, ida_err, <void*>self.pData)
             if flag < 0:
                 raise IDAError(flag, self.t)
                 
@@ -1843,7 +1846,10 @@ cdef class CVode(Explicit_ODE):
                     raise CVodeError(flag, self.t)
                     
             #Specify the error handling
-            flag = SUNDIALS.CVodeSetErrHandlerFn(self.cvode_mem, cv_err, <void*>self.pData)
+            IF SUNDIALS_VERSION >= (7,0,0):
+                flag = SUNDIALS.SUNContext_PushErrHandler(ctx, cv_err, <void*>self.pData)
+            ELSE:
+                flag = SUNDIALS.CVodeSetErrHandlerFn(self.cvode_mem, cv_err, <void*>self.pData)
             if flag < 0:
                 raise CVodeError(flag, self.t)