sciml / dassl.jl Goto Github PK

Is there a reason why you have the Travis tests using the Clang compiler instead of the normal Travis setup? You don't have any binary dependencies so it seems very unnecessary. Also, release is kept as v0.4 now, should set the environment separately. I could put this to a standard Travis script if you want, but I just want to double check that you aren't using the compiler for something special.

using DASSL
u0 = [1.0, 0, 0]
du0 = [-0.04, 0.04, 0.0]
tspan = (0.0,100000.0)

function resrob(r,yp,y,p,t)
r[1] = -0.04y[1] + 1.0e4y[2]y[3]
r[2] = -r[1] - 3.0e7y[2]*y[2] - yp[2]
r[1] -= yp[1]
r[3] = y[1] + y[2] + y[3] - 1.0
end

prob = DAEProblem(resrob,u0,du0,tspan)
sol = solve(prob, dassl())

ERROR: LinearAlgebra.LAPACKException(1)
Stacktrace:
[1] check_channel_state at .\channels.jl:120 [inlined]
[2] take_unbuffered(::Channel{Any}) at .\channels.jl:318
[3] take! at .\channels.jl:306 [inlined]
[4] iterate(::Channel{Any}, ::Nothing) at .\channels.jl:386
[5] iterate at .\channels.jl:385 [inlined]
[6] #dasslSolve#16(::Array{Float64,1}, ::Base.Iterators.Pairs{Symbol,Real,NTuple{7,Symbol},NamedTuple{(:abstol, :reltol, :maxstep, :minstep, :initstep, :maxorder, :factorize_jacobian),Tuple{Float64,Float64,Float64,Float64,Float64,Int64,Bool}}}, ::Function, ::Function, ::Array{Float64,1}, ::Array{Float64,1}) at C:\Users\hzgzh.julia\packages\DASSL\ZVHVR\src\DASSL.jl:163
[7] (::getfield(DASSL, Symbol("#kw##dasslSolve")))(::NamedTuple{(:abstol, :reltol, :maxstep, :minstep, :initstep, :maxorder, :factorize_jacobian),Tuple{Float64,Float64,Float64,Float64,Float64,Int64,Bool}}, ::typeof(dasslSolve), ::Function, ::Array{Float64,1}, ::Array{Float64,1}) at .\none:0
[8] #solve#2(::Bool, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Nothing, ::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}, ::typeof(solve), ::DAEProblem{Array{Float64,1},Array{Float64,1},Tuple{Float64,Float64},true,Nothing,DAEFunction{true,typeof(resrob),Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Nothing,Nothing}, ::dassl) at C:\Users\hzgzh.julia\packages\DASSL\ZVHVR\src\common.jl:36
[9] solve(::DAEProblem{Array{Float64,1},Array{Float64,1},Tuple{Float64,Float64},true,Nothing,DAEFunction{true,typeof(resrob),Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Nothing,Nothing}, ::dassl)
at C:\Users\hzgzh.julia\packages\DASSL\ZVHVR\src\common.jl:15
[10] top-level scope at none:0

@JuliaRegistrator register()

I think AppVeyor testing should be enabled on all of the JuliaDiffEq repositories in addition to Travis. Windows is a large set of users (myself included) and AppVeyor many times will catch Windows-specific issues (usually dealing with the way integers can be different, but also sometimes with how roundoff is handled or obscure facts that are hard to catch by yourself). This code is simple enough that I don't expect it to be different than the Travis results, but it's a very simple no-cost double-check and so I think it should be enabled.

Since the common API has pretty much settled, we should work towards getting it mapped onto DASSL. It shouldn't be difficult at all to get a pretty basic setup together.

Solving the van der Pol (VDP) equation with stringent error tolerances (reltol=abstol<0.2e-8) leads to a "step size too small" error. However, the tests on http://www.dm.uniba.it/~testset/problems/vdpol.php#stats show that DASSL works for VDP up to 10 significant digits, with DASSL.jl I only get 7 significant digits.

Test case:

# VDPOL with DASSL:  works
using DASSL
using ODE

# ode23s works at least as low as rtol=atol=1e-12 producing a relative error of 2e-9
# DASSL errors with  rtol=atol=1e-0 with too small step size
atol = .2e-8
rtol = .2e-8

dof=2
mu = 1e3
eps = 1/mu^2 # rescaled parameter
function fn(t,y)
    # the ode function dydt = f(t,y)
    out = zeros(Float64,dof)
    out[1] = y[2]
    prod = 
    out[2] = ( (1.0 - y[1]^2)*y[2] - y[1]) / eps
    return out
end
function jac(t,y)
    # the Jacobian of f
    #
    # Set to nothing if it is not known
    dfdy = zeros(Float64,dof,dof)
    dfdy[1,1] = 0.0
    dfdy[1,2] = 1.0
    dfdy[2,1] = (-2.0 * y[1]*y[2]-1.0)/eps
    dfdy[2,2] = (1.0 - y[1]^2)/eps
    return dfdy
end

# implicit eqs:
fni(t,y,dy)  = fn(t,y)-dy
Fy(t,y,dy) = jac(t,y)
Fdy(t,y,dy) = -eye(dof)

y0 = [2., 0]
tspan = [0., 2] # integration interval
(tn,yn_nojac,dyn)=DASSL.dasslSolve(fni, y0, tspan, reltol=atol, abstol=rtol)
(tn,yn,dyn)=DASSL.dasslSolve(fni, y0, tspan, Fy = Fy, Fdy = Fdy, reltol=atol, abstol=rtol)
tout,yout = ode23s(fn, y0, tspan, jacobian=jac, reltol=atol, abstol=rtol)

# reference:
refsol = [0.1706167732170483e1, -0.8928097010247975e0] # reference solution at tspan[2]

relerror(sol, refsol) = norm(abs((sol-refsol)./refsol),Inf)
calc_error_scd(sol, refsol) = -log10(norm(relerror(sol,refsol), Inf))

@show relerror(yn_nojac[end],refsol)
@show relerror(yn[end],refsol)
@show relerror(yout[end],refsol)

@show calc_error_scd(yn_nojac[end],refsol)
@show calc_error_scd(yn[end],refsol)
@show calc_error_scd(yout[end],refsol)

produces:

>> julia vdpol_dassl.jl
relerror(yn_nojac[end],refsol) => 3.714141830332293e-8
relerror(yn[end],refsol) => 4.626878704902238e-8
relerror(yout[end],refsol) => 3.559129266522373e-7
calc_error_scd(yn_nojac[end],refsol) => 7.430141516063499
calc_error_scd(yn[end],refsol) => 7.334711885518931
calc_error_scd(yout[end],refsol) => 6.4486562382625605

Remove the definitions of InterpolateAt, InterpolateDerivativeAt etc. from DASSL.jl and move them to the InterPoly.jl package (or other suitable Julia package).

This issue is being filed by a script, but if you reply, I will see it.

PackageEvaluator.jl is a script that runs nightly. It attempts to load all Julia packages and run their test (if available) on both the stable version of Julia (0.2) and the nightly build of the unstable version (0.3).

The results of this script are used to generate a package listing enhanced with testing results.

The status of this package, DASSL, on...

• Julia 0.2 is 'Tests fail, but package loads.'
• Julia 0.3 is 'Tests fail, but package loads.'

'No tests, but package loads.' can be due to their being no tests (you should write some if you can!) but can also be due to PackageEvaluator not being able to find your tests. Consider adding a test/runtests.jl file.

'Package doesn't load.' is the worst-case scenario. Sometimes this arises because your package doesn't have BinDeps support, or needs something that can't be installed with BinDeps. If this is the case for your package, please file an issue and an exception can be made so your package will not be tested.

This automatically filed issue is a one-off message. Starting soon, issues will only be filed when the testing status of your package changes in a negative direction (gets worse). If you'd like to opt-out of these status-change messages, reply to this message.

It would be nice if a function could be passed in for the Jacobian instead of hard-coded finite differences.

I would like to move this to JuliaDiffEq so it could be updated to v0.5 testing/compatibility and be added as a backend to AlgebraicDiffEq.jl. What do you think about that?

If MAXORDER constant is changed to 3 the BoundsError() shows up.

File src/DASSL.jl line 686 has the function:

# generate a function that computes approximate jacobian using forward
# finite differences
function numerical_jacobian(F,reltol,abstol,weights)
    function numjac(t, y, dy, a)
        ep      = eps(one(t))   # this is the machine epsilon
        h       = 1/a           # h ~ 1/a
        wt      = weights(y,reltol,abstol)
        # delta for approximation of jacobian.  I removed the
        # sign(h_next*dy0) from the definition of delta because it was
        # causing trouble when dy0==0 (which happens for ord==1)
        edelta  = diagm(0=>max.(abs.(y),abs.(h*dy),wt)*sqrt(ep))

        b=dy-a*y
        f(y1) = F(t,y1,a*y1+b)

        n   = length(y)
        jac = Array{eltype(y)}(undef,n,n)
        for i=1:n
            jac[:,i]=(f(y+edelta[:,i])-f(y))/edelta[i,i]
        end

        jac
    end
end

This (central) function should be improved:

1. In the for loop f(y) is called in every iteration, although y does not change. This means that there are n-1 unnecessary calls of the right hand side function. The call of f(y) should be moved out of the for loop and the result stored in a vector.

2. The statement jac = Array{eltype(y)}(undef,n,n) generates the Jacobian matrix whenever this function is called. So, when this function is called 100 times, then memory for 100 Jacobians are allocated. This should be improved: At the start of the simulation, memory for one Jacobian should be allocated that is then updated within this function.

On the mailing list, someone was trying to call DASSL from Python via pyjulia. This should work, but the Python function gets passed as a (callable) ::PyObject argument, not a ::Function.

It would be better to leave the arguments untyped rather than specify ::Function, in order that any callable object can be passed. Note that there is no performance advantage to declaring argument types; the main reason to declare arguments is to avoid method ambiguities.

I tried to run the convergence test:

...v0.3/DASSL/test((HEAD detached at a17ab3f))  >> julia convergence.jl 
Warning: could not import Base.Text into Tk
INFO: Stepsize too small (h=3.637978807091713e-16 at t=0.0.
INFO: Stepsize too small (h=3.637978807091713e-16 at t=0.0.
ERROR: dimensions must match
 in - at array.jl:723
 in dasslTestConvergence at /home/mauro/.julia/v0.3/DASSL/test/convergence.jl:35
 in dasslDrawConvergence at /home/mauro/.julia/v0.3/DASSL/test/convergence.jl:59
 in include at ./boot.jl:245
 in include_from_node1 at loading.jl:128
 in process_options at ./client.jl:285
 in _start at ./client.jl:354
while loading /home/mauro/.julia/v0.3/DASSL/test/convergence.jl, in expression starting on line 114

This issue is used to trigger TagBot; feel free to unsubscribe.

If you haven't already, you should update your TagBot.yml to include issue comment triggers.
Please see this post on Discourse for instructions and more details.

If you'd like for me to do this for you, comment TagBot fix on this issue.
I'll open a PR within a few hours, please be patient!