diff --git a/src/lib/sundials_includes.pxd b/src/lib/sundials_includes.pxd index 671bc3db..9b4394f0 100644 --- a/src/lib/sundials_includes.pxd +++ b/src/lib/sundials_includes.pxd @@ -41,6 +41,7 @@ IF SUNDIALS_VERSION >= (6,0,0): int SUNContext_Create(SUNComm comm, SUNContext* ctx) noexcept ELSE: int SUNContext_Create(void* comm, SUNContext* ctx) noexcept + int SUNContext_Free(SUNContext* ctx) noexcept IF SUNDIALS_VERSION >= (7,0,0): cdef extern from "sundials/sundials_context.h": @@ -294,6 +295,9 @@ ELSE: cdef inline int with_superlu() noexcept: return 0 IF SUNDIALS_VERSION >= (4,0,0): + cdef extern from "sundials/sundials_nonlinearsolver.h": + int SUNNonlinSolFree(SUNNonlinearSolver NLS) noexcept + cdef extern from "cvodes/cvodes.h": IF SUNDIALS_VERSION >= (6,0,0): void* CVodeCreate(int lmm, SUNContext ctx) noexcept diff --git a/src/ode_event_locator.c b/src/ode_event_locator.c index d9025721..1300911f 100644 --- a/src/ode_event_locator.c +++ b/src/ode_event_locator.c @@ -85,7 +85,8 @@ int f_event_locator(int n_y, int n_g, double TOL_inp, double t_low, double *t_hi imax = 0; for (i = 0; i < n_g; i++){ if ((g_low[i] > 0) != (g_high[i] > 0)){ - gfrac = _abs(g_high[i]/(g_low[i] - g_high[i])); + double denom = g_low[i] - g_high[i]; + gfrac = (denom != 0.0) ? _abs(g_high[i]/denom) : 0.0; if (gfrac >= maxfrac){ maxfrac = gfrac; imax = i; @@ -97,7 +98,12 @@ int f_event_locator(int n_y, int n_g, double TOL_inp, double t_low, double *t_hi if ((g_high[imax] == 0) || (g_low[imax] == 0)){ t_mid = (t_low + *t_high)/2.; }else{ - t_mid = *t_high - (*t_high - t_low)*g_high[imax]/(g_high[imax] - alpha*g_low[imax]); + double denom = g_high[imax] - alpha*g_low[imax]; + if (denom != 0.0){ + t_mid = *t_high - (*t_high - t_low)*g_high[imax]/denom; + } else { + t_mid = (t_low + *t_high)/2.; + } } /* Check if t_mid is to close to current brackets and adjust inwards if so is the case. */ diff --git a/src/solvers/kinsol.pyx b/src/solvers/kinsol.pyx index 341f75a1..d4111fbb 100644 --- a/src/solvers/kinsol.pyx +++ b/src/solvers/kinsol.pyx @@ -107,11 +107,15 @@ cdef class KINSOL(Algebraic): def __dealloc__(self): if self.y_temp != NULL: - #Deallocate N_Vector N_VDestroy(self.y_temp) + if self.y_scale != NULL: + N_VDestroy(self.y_scale) + + if self.f_scale != NULL: + N_VDestroy(self.f_scale) + if self.kinsol_mem != NULL: - #Free Memory SUNDIALS.KINFree(&self.kinsol_mem) IF SUNDIALS_VERSION >= (3,0,0): diff --git a/src/solvers/sundials.pyx b/src/solvers/sundials.pyx index e0843152..48cbf165 100644 --- a/src/solvers/sundials.pyx +++ b/src/solvers/sundials.pyx @@ -173,18 +173,31 @@ cdef class IDA(Implicit_ODE): def __dealloc__(self): if self.yTemp != NULL: - #Deallocate N_Vector N_VDestroy(self.yTemp) if self.ydTemp != NULL: - #Deallocate N_Vector N_VDestroy(self.ydTemp) if self.nv_atol != NULL: N_VDestroy(self.nv_atol) + if self.pData.dimSens > 0: + if self.pData.p != NULL: + free(self.pData.p) + if self.pData.pbar != NULL: + free(self.pData.pbar) + if self.ySO != NULL: + for i in range(self.pData.dimSens): + if self.ySO[i] != NULL: + N_VDestroy(self.ySO[i]) + free(self.ySO) + if self.ydSO != NULL: + for i in range(self.pData.dimSens): + if self.ydSO[i] != NULL: + N_VDestroy(self.ydSO[i]) + free(self.ydSO) + if self.ida_mem != NULL: - #Free Memory SUNDIALS.IDAFree(&self.ida_mem) IF SUNDIALS_VERSION >= (3,0,0): @@ -546,6 +559,8 @@ cdef class IDA(Implicit_ODE): #Set stop time flag = SUNDIALS.IDASetStopTime(self.ida_mem, tf) if flag < 0: + N_VDestroy(yout) + N_VDestroy(ydout) raise IDAError(flag, t) if opts["report_continuously"] or opts["output_list"] is None: @@ -596,6 +611,8 @@ cdef class IDA(Implicit_ODE): #Integration loop flag = SUNDIALS.IDASolve(self.ida_mem,tout,&tret,yout,ydout,IDA_NORMAL) if flag < 0: + N_VDestroy(yout) + N_VDestroy(ydout) raise IDAError(flag, tret) #Store results @@ -650,11 +667,15 @@ cdef class IDA(Implicit_ODE): #Set stop time flag = SUNDIALS.IDASetStopTime(self.ida_mem, tf) if flag < 0: + N_VDestroy(yout) + N_VDestroy(ydout) raise IDAError(flag, t) #Integration loop flag = SUNDIALS.IDASolve(self.ida_mem,tf,&tret,yout,ydout,IDA_ONE_STEP) if flag < 0: + N_VDestroy(yout) + N_VDestroy(ydout) raise IDAError(flag, tret) #Store results @@ -761,9 +782,17 @@ cdef class IDA(Implicit_ODE): flag = SUNDIALS.IDAGetErrWeights(self.ida_mem, eweight) if flag < 0: + N_VDestroy(ele) + N_VDestroy(eweight) + IF SUNDIALS_VERSION >= (6,0,0): + SUNDIALS.SUNContext_Free(&ctx) raise IDAError(flag) flag = SUNDIALS.IDAGetEstLocalErrors(self.ida_mem, ele) if flag < 0: + N_VDestroy(ele) + N_VDestroy(eweight) + IF SUNDIALS_VERSION >= (6,0,0): + SUNDIALS.SUNContext_Free(&ctx) raise IDAError(flag) pyweight = nv2arr(eweight) @@ -771,17 +800,14 @@ cdef class IDA(Implicit_ODE): err = pyweight*pyele - N_VDestroy(ele) #Deallocate - N_VDestroy(eweight) #Deallocate + N_VDestroy(ele) + N_VDestroy(eweight) + IF SUNDIALS_VERSION >= (6,0,0): + SUNDIALS.SUNContext_Free(&ctx) return err cpdef np.ndarray interpolate(self,double t,int k = 0): - """ - Calls the internal IDAGetDky for the interpolated values at time t. - t must be within the last internal step. k is the derivative of y which - can be from zero to the current order. - """ cdef flag cdef np.ndarray res IF SUNDIALS_VERSION >= (6,0,0): @@ -798,36 +824,20 @@ cdef class IDA(Implicit_ODE): flag = SUNDIALS.IDAGetDky(self.ida_mem, t, k, dky) if flag < 0: + N_VDestroy(dky) + IF SUNDIALS_VERSION >= (6,0,0): + SUNDIALS.SUNContext_Free(&ctx) raise IDAError(flag, t) res = nv2arr(dky) - N_VDestroy(dky) #Deallocate + N_VDestroy(dky) + IF SUNDIALS_VERSION >= (6,0,0): + SUNDIALS.SUNContext_Free(&ctx) return res cpdef interpolate_sensitivity(self,double t, int k = 0, int i=-1): - """ - This method calls the internal method IDAGetSensDky which computes the k-th derivatives - of the interpolating polynomials for the sensitivity variables at time t. - - Parameters:: - - t - - Specifies the time at which sensitivity information is requested. The time - t must fall within the interval defined by the last successful step taken - by IDAS. - - k - - The order of derivatives. - - i - - Specifies the sensitivity derivative vector to be returned (0<=i<=Ns) - - Return:: - - A matrix containing the Ns vectors or a vector if i is specified. - """ IF SUNDIALS_VERSION >= (6,0,0): cdef SUNDIALS.SUNContext ctx = NULL IF SUNDIALS_VERSION >= (7,0,0): @@ -849,22 +859,32 @@ cdef class IDA(Implicit_ODE): flag = SUNDIALS.IDAGetSensDky1(self.ida_mem, t, k, x, dkyS) if flag<0: + N_VDestroy(dkyS) + IF SUNDIALS_VERSION >= (6,0,0): + SUNDIALS.SUNContext_Free(&ctx) raise IDAError(flag, t) matrix += [nv2arr(dkyS)] N_VDestroy(dkyS) + IF SUNDIALS_VERSION >= (6,0,0): + SUNDIALS.SUNContext_Free(&ctx) return np.array(matrix) else: flag = SUNDIALS.IDAGetSensDky1(self.ida_mem, t, k, i, dkyS) if flag <0: + N_VDestroy(dkyS) + IF SUNDIALS_VERSION >= (6,0,0): + SUNDIALS.SUNContext_Free(&ctx) raise IDAError(flag, t) res = nv2arr(dkyS) N_VDestroy(dkyS) + IF SUNDIALS_VERSION >= (6,0,0): + SUNDIALS.SUNContext_Free(&ctx) return res @@ -1606,7 +1626,6 @@ cdef class CVode(Explicit_ODE): def __dealloc__(self): if self.yTemp != NULL: - #Deallocate N_Vector N_VDestroy(self.yTemp) if self.nv_atol != NULL: @@ -1615,8 +1634,24 @@ cdef class CVode(Explicit_ODE): if self.nv_rtol != NULL: N_VDestroy(self.nv_rtol) + if self.pData.dimSens > 0: + if self.pData.p != NULL: + free(self.pData.p) + if self.pData.pbar != NULL: + free(self.pData.pbar) + if self.ySO != NULL: + for i in range(self.pData.dimSens): + if self.ySO[i] != NULL: + N_VDestroy(self.ySO[i]) + free(self.ySO) + + IF SUNDIALS_VERSION >= (4,0,0): + if self.sun_nonlinearsolver != NULL: + SUNDIALS.SUNNonlinSolFree(self.sun_nonlinearsolver) + if self.sun_nonlinearsolver_sens != NULL: + SUNDIALS.SUNNonlinSolFree(self.sun_nonlinearsolver_sens) + if self.cvode_mem != NULL: - #Free Memory SUNDIALS.CVodeFree(&self.cvode_mem) IF SUNDIALS_VERSION >= (3,0,0): @@ -1627,9 +1662,6 @@ cdef class CVode(Explicit_ODE): SUNDIALS.SUNLinSolFree(self.sun_linearsolver) cpdef get_local_errors(self): - """ - Returns the vector of estimated local errors at the current step. - """ cdef int flag if self.cvode_mem == NULL: @@ -1644,16 +1676,20 @@ cdef class CVode(Explicit_ODE): SUNDIALS.SUNContext_Create(comm, &ctx) cdef N_Vector ele=N_VNew_Serial(self.pData.dim, ctx) ELSE: - cdef N_Vector ele=N_VNew_Serial(self.pData.dim) #Allocates a new N_Vector + cdef N_Vector ele=N_VNew_Serial(self.pData.dim) flag = SUNDIALS.CVodeGetEstLocalErrors(self.cvode_mem, ele) if flag < 0: + N_VDestroy(ele) + IF SUNDIALS_VERSION >= (6,0,0): + SUNDIALS.SUNContext_Free(&ctx) raise CVodeError(flag, self.t) ele_py = nv2arr(ele) - #Deallocate N_Vector N_VDestroy(ele) + IF SUNDIALS_VERSION >= (6,0,0): + SUNDIALS.SUNContext_Free(&ctx) return ele_py @@ -1722,9 +1758,6 @@ cdef class CVode(Explicit_ODE): return qcur cpdef get_error_weights(self): - """ - Returns the solution error weights at the current step. - """ cdef int flag if self.cvode_mem == NULL: @@ -1739,16 +1772,20 @@ cdef class CVode(Explicit_ODE): SUNDIALS.SUNContext_Create(comm, &ctx) cdef N_Vector eweight=N_VNew_Serial(self.pData.dim, ctx) ELSE: - cdef N_Vector eweight=N_VNew_Serial(self.pData.dim) #Allocates a new N_Vector + cdef N_Vector eweight=N_VNew_Serial(self.pData.dim) flag = SUNDIALS.CVodeGetErrWeights(self.cvode_mem, eweight) if flag < 0: + N_VDestroy(eweight) + IF SUNDIALS_VERSION >= (6,0,0): + SUNDIALS.SUNContext_Free(&ctx) raise CVodeError(flag, self.t) eweight_py = nv2arr(eweight) - #Deallocate N_Vector N_VDestroy(eweight) + IF SUNDIALS_VERSION >= (6,0,0): + SUNDIALS.SUNContext_Free(&ctx) return eweight_py @@ -1963,11 +2000,6 @@ cdef class CVode(Explicit_ODE): cpdef np.ndarray interpolate(self,double t,int k = 0): - """ - Calls the internal CVodeGetDky for the interpolated values at time t. - t must be within the last internal step. k is the derivative of y which - can be from zero to the current order. - """ cdef flag cdef np.ndarray res IF SUNDIALS_VERSION >= (6,0,0): @@ -1979,42 +2011,25 @@ cdef class CVode(Explicit_ODE): SUNDIALS.SUNContext_Create(comm, &ctx) cdef N_Vector dky=N_VNew_Serial(self.pData.dim, ctx) ELSE: - cdef N_Vector dky=N_VNew_Serial(self.pData.dim) #Allocates a new N_Vector + cdef N_Vector dky=N_VNew_Serial(self.pData.dim) flag = SUNDIALS.CVodeGetDky(self.cvode_mem, t, k, dky) if flag < 0: + N_VDestroy(dky) + IF SUNDIALS_VERSION >= (6,0,0): + SUNDIALS.SUNContext_Free(&ctx) raise CVodeError(flag, t) res = nv2arr(dky) - #Deallocate N_Vector N_VDestroy(dky) + IF SUNDIALS_VERSION >= (6,0,0): + SUNDIALS.SUNContext_Free(&ctx) return res cpdef np.ndarray interpolate_sensitivity(self, realtype t, int k = 0, int i=-1): - """ - This method calls the internal method CVodeGetSensDky which computes the k-th derivatives - of the interpolating polynomials for the sensitivity variables at time t. - - Parameters:: - - t - - Specifies the time at which sensitivity information is requested. The time - t must fall within the interval defined by the last successful step taken - by CVodeS. - - k - - The order of derivatives. - - i - - Specifies the sensitivity derivative vector to be returned (0<=i<=Ns) - - Return:: - - A matrix containing the Ns vectors or a vector if i is specified. - """ IF SUNDIALS_VERSION >= (6,0,0): cdef SUNDIALS.SUNContext ctx = NULL IF SUNDIALS_VERSION >= (7,0,0): @@ -2035,21 +2050,31 @@ cdef class CVode(Explicit_ODE): for x in range(self.pData.dimSens): flag = SUNDIALS.CVodeGetSensDky1(self.cvode_mem, t, k, x, dkyS) if flag<0: + N_VDestroy(dkyS) + IF SUNDIALS_VERSION >= (6,0,0): + SUNDIALS.SUNContext_Free(&ctx) raise CVodeError(flag, t) matrix += [nv2arr(dkyS)] N_VDestroy(dkyS) + IF SUNDIALS_VERSION >= (6,0,0): + SUNDIALS.SUNContext_Free(&ctx) return np.array(matrix) else: flag = SUNDIALS.CVodeGetSensDky1(self.cvode_mem, t, k, i, dkyS) if flag <0: + N_VDestroy(dkyS) + IF SUNDIALS_VERSION >= (6,0,0): + SUNDIALS.SUNContext_Free(&ctx) raise CVodeError(flag, t) res = nv2arr(dkyS) N_VDestroy(dkyS) + IF SUNDIALS_VERSION >= (6,0,0): + SUNDIALS.SUNContext_Free(&ctx) return res @@ -2085,11 +2110,13 @@ cdef class CVode(Explicit_ODE): #Set stop time flag = SUNDIALS.CVodeSetStopTime(self.cvode_mem, tf) if flag < 0: + N_VDestroy(yout) raise CVodeError(flag, t) #Integration loop flag = SUNDIALS.CVode(self.cvode_mem,tf,yout,&tret,CV_ONE_STEP) if flag < 0: + N_VDestroy(yout) raise CVodeError(flag, tret) #Store results