Merge branch 'master' into Assimulo-3.2.x

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

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

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
 # "<project> v<release> 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

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",

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 <http://www.gnu.org/licenses/>.
+
+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()

src/lib/sundials_callbacks.pxi

@@ -55,13 +55,13 @@ cdef inline N.ndarray nv2arr(N_Vector v):
     cdef long int n = (<N_VectorContent_Serial>v.content).length
     cdef realtype* v_data = (<N_VectorContent_Serial>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 = (<N_VectorContent_Serial>v.content).length
     cdef realtype* v_data = (<N_VectorContent_Serial>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

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):

thirdparty/hairer/radau5_c_py.pyx

@@ -53,15 +53,15 @@ 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)
 cdef void c2py_mat_F(np.ndarray[double, ndim=2, mode='fortran'] dest, double* source, int dim):
     """
     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, <void*>fcn_PY, &x, &y_vec[0], &xend,
-                         &h__, &rtol_vec[0], &rtol_vec[0], &itol, callback_jac, <void*> jac_PY,
+                         &h__, &rtol_vec[0], &atol_vec[0], &itol, callback_jac, <void*> jac_PY,
                          &ijac, &mljac, &mujac, callback_mas, <void*> mas_PY, &imas, &mlmas, &mumas,
                          callback_solout, <void*>solout_PY, &iout, &work_vec[0], &lwork, &iwork_vec[0], &liwork, &rpar,
                          &ipar, &idid)