Dag Sverre Seljebotn <da...@...> writes:
> The challenge here though is about function pointers.
>
> Basically you need to define
> cdef extern from ...:
> ctypedef double realtype
> ctypedef struct N_Vector:
> ...more or less copy definition from nvector header files...
> ctypedef int (*CVRhsFn)(realtype t, N_Vector y, N_Vector ydot, void
> *user_data)
> int CVodeInit(void* mem, CVRhsFn f, realtype t0, N_Vector y0)
>
> Note that here CVRhsFn is a /type/, it is a pointer to a function, and
> CVodeInit is a function which takes such a pointer (a "callback"). So
> you'd do
>
> cdef int your_rhs(realtype t, N_Vector y, N_vector ydot, void* user_data):
> ...
>
> CVodeInit(..., your_rhs, t0, y0)
>
> and now Sundials will know that your_rhs should be called.
I'm pleased to report that Sundials is finally doing my bidding from within
Cython [1]! As I become more familiar with the cdef extern syntax, it's fairly
easy to import new functions from Sundials/CVODE. I'll come back to efficient
indexing [2] later.
If I may make a wish, I would like Cython to find more declarations
automatically. For instance, having
cdef extern * from "cvode/cvode.h"
or, if that is impractical,
cdef extern from "cvode/cvode.h":
CVodeCreate, CVRhsFn, CVodeMalloc, CVodeFree, CVode, CV_ADAMS, CV_BDF
instead of:
cdef extern from "cvode/cvode.h":
void* CVodeCreate(int lmm, int iter)
ctypedef int (*CVRhsFn)(realtype t, N_Vector y, N_Vector ydot, void *f_data)
int CVodeMalloc(void *cvode_mem, CVRhsFn f, realtype t0, N_Vector y0,
int itol, realtype reltol, void *abstol)
void CVodeFree(void **cvode_mem)
int CVode(void *cvode_mem, realtype tout, N_Vector yout, realtype *tret,
int itask)
# constants
## /* lmm */
DEF CV_ADAMS = 1
DEF CV_BDF = 2
...
I haven't benchmarked anything yet, but the relative ease of adapting the C
tutorial for CVODE [3] is uplifting. As a last question, what is the cleanest
way to transfer arrays of double between realtype* and Numpy vectors (one-
dimensional arrays)?
Thank you very much to all who helped me!
Best regards,
Jon Olav
[1]
== setup.py ==
"""Run with e.g. python setup.py build_ext --inplace
Or, to remove intermediate files from previous builds:
rm -rf build hello.so hello.c && python setup.py build_ext --inplace && python -
c "import hello"
"""
from distutils.core import setup
from distutils.extension import Extension
from Cython.Distutils import build_ext
import numpy as np
setup(
cmdclass = {'build_ext': build_ext},
ext_modules = [
Extension("hello",
["hello.pyx"],
include_dirs=['/site/VERSIONS/sundials-2.3.0/include',
np.get_include ()],
library_dirs=['/site/VERSIONS/sundials-2.3.0/lib'],
libraries=['sundials_cvode', 'sundials_nvecserial']),
]
)
== hello.pyx ==
"""
Cythonizing simple usage of CVODE
Build and test with
rm -rf build hello.so hello.c && python setup.py build_ext --inplace && python -
c "import hello"
"""
cdef extern from "sundials/sundials_types.h":
ctypedef double realtype
cdef extern from "sundials/sundials_nvector.h":
cdef struct _generic_N_Vector:
void* content
ctypedef _generic_N_Vector* N_Vector
N_Vector N_VNew_Serial(long int vec_length)
void N_VDestroy_Serial(N_Vector v)
void N_VPrint_Serial(N_Vector v)
cdef extern from "nvector/nvector_serial.h":
cdef struct _N_VectorContent_Serial:
long int length
realtype* data
ctypedef _N_VectorContent_Serial* N_VectorContent_Serial
cdef extern from "cvode/cvode.h":
void* CVodeCreate(int lmm, int iter)
ctypedef int (*CVRhsFn)(realtype t, N_Vector y, N_Vector ydot, void *f_data)
int CVodeMalloc(void *cvode_mem, CVRhsFn f, realtype t0, N_Vector y0,
int itol, realtype reltol, void *abstol)
void CVodeFree(void **cvode_mem)
int CVode(void *cvode_mem, realtype tout, N_Vector yout, realtype *tret,
int itask)
# constants
## /* lmm */
DEF CV_ADAMS = 1
DEF CV_BDF = 2
## /* iter */
DEF CV_FUNCTIONAL = 1
DEF CV_NEWTON = 2
## /* itol */
DEF CV_SS = 1
DEF CV_SV = 2
DEF CV_WF = 3
## /* itask */
DEF CV_NORMAL = 1
DEF CV_ONE_STEP = 2
DEF CV_NORMAL_TSTOP = 3
DEF CV_ONE_STEP_TSTOP = 4
## /* CVODE return flags */
DEF CV_SUCCESS = 0
DEF CV_TSTOP_RETURN = 1
DEF CV_ROOT_RETURN = 2
# how to find this stuff:
# grep "CVDense(" /xanadu/site/common/VERSIONS/sundials-2.3.0/include/*/*
# /xanadu/site/common/VERSIONS/sundials-2.3.0/include/
# ... cvode/cvode_dense.h:int CVDense(void *cvode_mem, long int N);
cdef extern from "cvode/cvode_dense.h":
int CVDense(void *cvode_mem, long int N)
import numpy as np
cdef nv2arr(N_Vector v):
"""Create new numpy array from N_Vector"""
cdef long int n = (<N_VectorContent_Serial>v.content).length
cdef realtype* v_data = (<N_VectorContent_Serial>v.content).data
cdef long int i
x = np.empty(n)
for i in range(n):
x[i] = v_data[i]
return x
cimport numpy as np
# cdef N_Vector arr2nv(np.ndarray[np.int_t] x):
cdef N_Vector arr2nv(x):
"""Create new N_Vector from numpy array"""
cdef long int n = len(x)
cdef N_Vector v = N_VNew_Serial(n)
cdef realtype* v_data = (<N_VectorContent_Serial>v.content).data
cdef long int i
for i in range(n):
v_data[i] = x[i]
return v
# following https://computation.llnl.gov/casc/sundials/documentation/
# ... cv_guide/node5.html#SECTION00540000000000000000
# set problem dimensions
N = 2
# set vector of initial values
cdef N_Vector y0 = arr2nv(np.arange(N)) # N_VMake_Serial(N, ydata)
# create CVODE object
cdef void* cvode_mem = CVodeCreate(CV_BDF, CV_NEWTON)
## ODE right-hand side
cdef int vanderpol (realtype t, N_Vector y_nv, N_Vector ydot_nv, void *f_data):
"""van der Pol equation adapted from example_ode.py"""
cdef realtype* y = (<N_VectorContent_Serial>y_nv.content).data
cdef realtype* ydot = (<N_VectorContent_Serial>ydot_nv.content).data
ydot[0] = y[1]
ydot[1] = (1 - y[0] * y[0]) * y[1] - y[0] # eps[0] * ...
return 0
### testing RHS
print "Current state:"
N_VPrint_Serial(y0) # 0, 1
cdef N_Vector ydot = arr2nv(np.zeros(N)) # N_VMake_Serial(N, ydata)
vanderpol(0.0, y0, ydot, <void*>None)
print "Rates of change:"
N_VPrint_Serial(ydot) # 1, 1
# initialize CVODE solver
# specify integration tolerances
# (note: Sundials 2.3 uses CVodeMalloc instead of CVodeInit + CVodeSStolerances)
cdef realtype abstol = 1e-4
flag = CVodeMalloc(cvode_mem, vanderpol, 0.0, y0, CV_SS, 1e-4, &abstol)
# if flag == CV_SUCCESS: print "CVodeMalloc returned CV_SUCCESS"
# set optional inputs
# attach linear solver module
flag = CVDense(cvode_mem, N)
# set linear solver optional inputs
# specify rootfinding problem
# advance solution in time
tmax = 150
y = np.zeros((tmax, N))
t = np.zeros((tmax,))
cdef realtype tout = 2.0
cdef realtype tret
cdef int i, j
cdef realtype* v_data = (<N_VectorContent_Serial>y0.content).data
for i in range(tmax):
flag = CVode(cvode_mem, tout, y0, &tret, CV_ONE_STEP)
if flag != CV_SUCCESS:
break
t[i] = tret
for j in range(N):
y[i,j] = v_data[j]
# get optional outputs
# deallocate memory for solution vector
N_VDestroy_Serial(ydot)
N_VDestroy_Serial(y0)
# free solver memory
CVodeFree(&cvode_mem)
# display results
import matplotlib as mpl
mpl.use("svg")
import pylab
pylab.plot(t, y)
pylab.savefig("cvode.svgz")
[2] http://docs.cython.org/docs/numpy_tutorial.html#efficient-indexing
[3] https://computation.llnl.gov/casc/sundials/documentation/cv_guide/
node5.html#SECTION00540000000000000000
_______________________________________________
Cython-dev mailing list
[email protected]
http://codespeak.net/mailman/listinfo/cython-dev
