diff --git a/CHANGELOG b/CHANGELOG index 3b02302f..b9f61ad3 100644 --- a/CHANGELOG +++ b/CHANGELOG @@ -1,5 +1,8 @@ --- CHANGELOG --- +--- Assimulo-3.2.9 --- + * Bugfix for C version of Radau5 solver + --- Assimulo-3.2.8 --- * Sundials 5.x port * Improved asserts in tests diff --git a/doc/sphinx/source/changelog.rst b/doc/sphinx/source/changelog.rst index fa7927e3..d500d8c9 100644 --- a/doc/sphinx/source/changelog.rst +++ b/doc/sphinx/source/changelog.rst @@ -3,6 +3,9 @@ Changelog ========== +--- Assimulo-3.2.9 --- + * Bugfix for C version of Radau5 solver + --- Assimulo-3.2.8 --- * Sundials 5.x port * Improved asserts in tests diff --git a/doc/sphinx/source/conf.py b/doc/sphinx/source/conf.py index 3d25dce6..5b91d8da 100644 --- a/doc/sphinx/source/conf.py +++ b/doc/sphinx/source/conf.py @@ -50,9 +50,9 @@ # built documents. # # The short X.Y version. -version = '3.2.8' +version = '3.2.9' # The full version, including alpha/beta/rc tags. -release = '3.2.8' +release = '3.2.9' # The language for content autogenerated by Sphinx. Refer to documentation # for a list of supported languages. @@ -109,7 +109,7 @@ # The name for this set of Sphinx documents. If None, it defaults to # " v documentation". -html_title = 'Assimulo 3.2.8 documentation' +html_title = 'Assimulo 3.2.9 documentation' # A shorter title for the navigation bar. Default is the same as html_title. #html_short_title = None diff --git a/examples/__init__.py b/examples/__init__.py index 7d9a2f98..2e318065 100644 --- a/examples/__init__.py +++ b/examples/__init__.py @@ -18,7 +18,7 @@ __all__ = ["cvode_gyro", "cvode_basic", "cvode_with_disc", "cvode_with_initial_sensitivity", "cvode_with_jac", "cvode_with_jac_spgmr", "cvode_with_parameters", "euler_basic", "euler_with_disc", "rungekutta4_basic", "rungekutta34_basic", "rungekutta34_with_disc", "ida_with_disc", "ida_with_initial_sensitivity", - "ida_with_jac", "ida_with_parameters", "radau5ode_vanderpol", "radau5ode_with_disc" , "ida_with_jac_spgmr", + "ida_with_jac", "ida_with_parameters", "radau5ode_vanderpol", "radau5ode_with_disc" , "ida_with_jac_spgmr", "radau5ode_with_jac_sparse", "radau5dae_vanderpol", "dopri5_basic", "dopri5_with_disc", "rodasode_vanderpol", "glimda_vanderpol", "lsodar_vanderpol", "lsodar_with_disc", "mech_system_pendulum", "euler_vanderpol", "cvode_with_parameters_modified", diff --git a/examples/radau5ode_with_jac_sparse.py b/examples/radau5ode_with_jac_sparse.py new file mode 100644 index 00000000..e76784df --- /dev/null +++ b/examples/radau5ode_with_jac_sparse.py @@ -0,0 +1,99 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- + +# Copyright (C) 2021 Modelon AB +# +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU Lesser General Public License as published by +# the Free Software Foundation, version 3 of the License. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU Lesser General Public License for more details. +# +# You should have received a copy of the GNU Lesser General Public License +# along with this program. If not, see . + +import numpy as N +import scipy.sparse as SP +import nose +from assimulo.solvers import Radau5ODE +from assimulo.problem import Explicit_Problem + + +def run_example(with_plots=True, solver = 'c'): + r""" + Example for demonstrating the use of a user supplied Jacobian (sparse). + Based on the SUNDIALS example cvRoberts_sps.c + + ODE: + + .. math:: + + \dot y_1 &= -0.04y_1 + 1e4 y_2 y_3 \\ + \dot y_2 &= - \dot y_1 - \dot y_3 \\ + \dot y_3 &= 3e7 y_2^2 + + + on return: + + - :dfn:`exp_mod` problem instance + + - :dfn:`exp_sim` solver instance + + """ + + #Defines the rhs + def f(t,y): + yd_0 = -0.04*y[0] + 1e4*y[1]*y[2] + yd_2 = 3e7*y[1]*y[1] + yd_1 = -yd_0 - yd_2 + return N.array([yd_0,yd_1,yd_2]) + + #Defines the Jacobian + def jac(t,y): + + colptrs = [0,3,6,9] + rowvals = [0, 1, 2, 0, 1, 2, 0, 1, 2] + data = [-0.04, 0.04, 0.0, 1e4*y[2], -1e4*y[2]-6e7*y[1], 6e7*y[1], 1e4*y[1], -1e4*y[1], 0.0] + + J = SP.csc_matrix((data, rowvals, colptrs)) + return J + + #Defines an Assimulo explicit problem + y0 = [1.0, 0.0, 0.0] #Initial conditions + + exp_mod = Explicit_Problem(f, y0, name = 'Example using analytic (sparse) Jacobian') + + exp_mod.jac = jac #Sets the Jacobian + exp_mod.jac_nnz = 9 + + + exp_sim = Radau5ODE(exp_mod) #Create a Radau5 solver + + #Set the parameters + exp_sim.atol = [1e-8,1e-14,1e-6] #Default 1e-6 + exp_sim.rtol = 1e-4 #Default 1e-6 + exp_sim.solver = solver + + #Simulate + t, y = exp_sim.simulate(0.4) #Simulate 0.4 seconds + + #Plot + if with_plots: + import pylab as P + P.plot(t,y[:,1],linestyle="dashed",marker="o") #Plot the solution + P.xlabel('Time') + P.ylabel('State') + P.title(exp_mod.name) + P.show() + + #Basic tests + nose.tools.assert_almost_equal(y[-1][0],0.9851,3) + + return exp_mod, exp_sim + + +if __name__=='__main__': + mod,sim = run_example() diff --git a/src/lib/sundials_callbacks.pxi b/src/lib/sundials_callbacks.pxi index d8e3a920..0a8db24f 100644 --- a/src/lib/sundials_callbacks.pxi +++ b/src/lib/sundials_callbacks.pxi @@ -55,13 +55,13 @@ cdef inline N.ndarray nv2arr(N_Vector v): cdef long int n = (v.content).length cdef realtype* v_data = (v.content).data cdef N.ndarray[realtype, ndim=1, mode='c'] x=N.empty(n) - memcpy(x.data, v_data, n*sizeof(realtype)) + memcpy(PyArray_DATA(x), v_data, n*sizeof(realtype)) return x cdef inline void nv2arr_inplace(N_Vector v, N.ndarray o): cdef long int n = (v.content).length cdef realtype* v_data = (v.content).data - memcpy(o.data, v_data, n*sizeof(realtype)) + memcpy(PyArray_DATA(o), v_data, n*sizeof(realtype)) cdef inline void nv2mat_inplace(int Ns, N_Vector *v, N.ndarray o): cdef long int i,j, Nf @@ -73,5 +73,5 @@ cdef inline void nv2mat_inplace(int Ns, N_Vector *v, N.ndarray o): cdef inline realtype2arr(realtype *data, int n): """Create new numpy array from realtype*""" cdef N.ndarray[realtype, ndim=1, mode='c'] x=N.empty(n) - memcpy(x.data, data, n*sizeof(realtype)) + memcpy(PyArray_DATA(x), data, n*sizeof(realtype)) return x diff --git a/tests/test_examples.py b/tests/test_examples.py index 724382f3..9088eec2 100644 --- a/tests/test_examples.py +++ b/tests/test_examples.py @@ -154,6 +154,14 @@ def test_radau5ode_with_disc_c(self): @testattr(stddist = True) def test_radau5ode_with_disc_f(self): radau5ode_with_disc.run_example(with_plots=False,solver='c') + + @testattr(stddist = True) + def test_radau5ode_with_jac_sparse_c(self): + radau5ode_with_jac_sparse.run_example(with_plots=False,solver='c') + + @testattr(stddist = True) + def test_radau5ode_with_jac_sparse_f(self): + radau5ode_with_jac_sparse.run_example(with_plots=False,solver='f') @testattr(stddist = True) def test_radau5dae_vanderpol_c(self): diff --git a/thirdparty/hairer/radau5_c_py.pyx b/thirdparty/hairer/radau5_c_py.pyx index 25d0e66a..210093bf 100644 --- a/thirdparty/hairer/radau5_c_py.pyx +++ b/thirdparty/hairer/radau5_c_py.pyx @@ -53,7 +53,7 @@ cdef void c2py(np.ndarray[double, ndim=1, mode='c'] dest, double* source, int di """ Copy (double *) C vector to 1D numpy array """ - memcpy(dest.data, source, dim*sizeof(double)) + memcpy(PyArray_DATA(dest), source, dim*sizeof(double)) @cython.boundscheck(False) @cython.wraparound(False) @@ -61,7 +61,7 @@ cdef void c2py_mat_F(np.ndarray[double, ndim=2, mode='fortran'] dest, double* so """ Copy (double *) C matrix (Fotran-style column major ordering) to 2D numpy array """ - memcpy(dest.data, source, dim*sizeof(double)) + memcpy(PyArray_DATA(dest), source, dim*sizeof(double)) cdef int callback_fcn(integer* n, doublereal* x, doublereal* y_in, doublereal* y_out, doublereal* rpar, integer* ipar, void* fcn_PY): @@ -219,7 +219,7 @@ cpdef radau5(fcn_PY, doublereal x, np.ndarray y, cdef np.ndarray[integer, mode="c", ndim=1] iwork_vec = iwork_in radau5_c_py.radau5_c(&n, callback_fcn, fcn_PY, &x, &y_vec[0], &xend, - &h__, &rtol_vec[0], &rtol_vec[0], &itol, callback_jac, jac_PY, + &h__, &rtol_vec[0], &atol_vec[0], &itol, callback_jac, jac_PY, &ijac, &mljac, &mujac, callback_mas, mas_PY, &imas, &mlmas, &mumas, callback_solout, solout_PY, &iout, &work_vec[0], &lwork, &iwork_vec[0], &liwork, &rpar, &ipar, &idid)