The parameter for the `cvode` function is `f` and so is the function you
want to use. These get confused, and it tries to use the "function" `J`
instead. Changing the parameter name to something other than `{f}` should
work
On Wednesday, December 16, 2015 at 4:26:41 PM UTC+2, Simon Frost wrote:
>
> Dear Julia Users,
>
> I'm trying to speed up some code that employs passing functions as
> arguments. One part of the code solves an ODE; if I use CVODE from
> Sundials, and rewrite the function to accept callable types, I get a
> ReadOnlyMemoryError - as I don't know the Sundials API, can someone help me
> with where I'm going wrong? Code below.
>
> Best
> Simon
>
> **
>
> using Sundials
>
> function cvode{f}(::Type{f}, y0::Vector{Float64}, t::Vector{Float64};
> reltol::Float64=1e-4, abstol::Float64=1e-6)
> neq = length(y0)
> mem = Sundials.CVodeCreate(Sundials.CV_BDF, Sundials.CV_NEWTON)
> flag = Sundials.CVodeInit(mem, cfunction(Sundials.cvodefun, Int32,
> (Sundials.realtype, Sundials.N_Vector, Sundials.N_Vector, Ref{Function})),
> t[1], Sundials.nvector(y0))
> flag = Sundials.CVodeSetUserData(mem, f)
> flag = Sundials.CVodeSStolerances(mem, reltol, abstol)
> flag = Sundials.CVDense(mem, neq)
> yres = zeros(length(t), length(y0))
> yres[1,:] = y0
> y = copy(y0)
> tout = [0.0]
> for k in 2:length(t)
> flag = Sundials.CVode(mem, t[k], y, tout, Sundials.CV_NORMAL)
> yres[k,:] = y
> end
> Sundials.CVodeFree([mem])
> return yres
> end
>
> function f(t, y, ydot)
> ydot[1] = 0.1*(-72-y[1])+0.1*1.4*exp((y[1]+48)/1.4)+10
> ydot[3] = 0.
> ydot[2] = 0.
> end
>
> immutable J; end
> call(::Type{J},t, y, ydot) = f(t, y, ydot)
>
> t = [0.1:0.0001:1]
> res = cvode(J, [-60.0, 0.0, 0.0], t);
>
