Merge branch '07'

Author: ChrisRackauckas

README.md

@@ -3,7 +3,6 @@ DASSL.jl
 
 [![Build Status](https://travis-ci.org/JuliaDiffEq/DASSL.jl.svg?branch=master)](https://travis-ci.org/JuliaDiffEq/DASSL.jl)
 [![Coverage Status](https://img.shields.io/coveralls/pwl/DASSL.jl.svg)](https://coveralls.io/r/pwl/DASSL.jl)
-[![DASSL](http://pkg.julialang.org/badges/DASSL_0.5.svg)](http://pkg.julialang.org/?pkg=DASSL)
 
 This is an implementation of DASSL algorithm for solving algebraic
 differential equations.  To inastall a stable version run
@@ -34,7 +33,7 @@ prob = DAEProblem(resrob,u0,du0,tspan)
 sol = solve(prob, dassl())
 ```
 
-For more details on using this interface, [see the ODE tutorial](http://docs.juliadiffeq.org/latest/tutorials/ode_example.html). 
+For more details on using this interface, [see the ODE tutorial](http://docs.juliadiffeq.org/latest/tutorials/ode_example.html).
 
 Examples
 --------

src/DASSL.jl

@@ -6,6 +6,8 @@ export dasslIterator, dasslSolve
 
 using Reexport
 @reexport using DiffEqBase
+using LinearAlgebra
+
 import DiffEqBase: solve
 
 export dassl
@@ -20,7 +22,8 @@ mutable struct JacData
     jac # Jacobian matrix for the newton solver
 end
 
-function dasslStep(F,
+function dasslStep(channel,
+                   F,
                    y0::AbstractVector{T},
                    tstart::Real;
                    reltol   = 1e-3,
@@ -45,9 +48,9 @@ function dasslStep(F,
     ord      = 1                # initial method order
     tout     = [tstart]         # initial time
     h        = initstep         # current step size
-    yout     = Array(typeof(y0),1)
+    yout     = Array{typeof(y0)}(undef,1)
     yout[1]  = y0
-    dyout    = Array(typeof(y0),1)
+    dyout    = Array{typeof(y0)}(undef,1)
     dyout[1] = dy0              # initial guess for dy0
     num_rejected = 0            # number of rejected steps
     num_fail = 0                # number of consecutive failures
@@ -128,12 +131,12 @@ function dasslStep(F,
 
             # remove old results
             if length(tout) > ord+3
-                shift!(tout)
-                shift!(yout)
-                shift!(dyout)
+                popfirst!(tout)
+                popfirst!(yout)
+                popfirst!(dyout)
             end
 
-            produce(tout[end],yout[end],dyout[end])
+            push!(channel,(tout[end],yout[end],dyout[end]))
 
             # determine the new step size and order, including the current step
             (r,ord) = newStepOrder(tout,yout,normy,err,ord,num_fail,maxorder)
@@ -147,14 +150,13 @@ end
 
 
 # the iterator version of the dasslSolve
-dasslIterator(F, y0, t0; args...) = Task(()->dasslStep(F, y0, t0; args...))
-
+dasslIterator(F, y0, t0; args...) = Channel((channel)->dasslStep(channel,F, y0, t0; args...))
 
 # solves the equation F with initial data y0 over for times t in tspan=[t0,t1]
 function dasslSolve(F, y0::AbstractVector, tspan; dy0 = zero(y0), args...)
-    tout  = Array(typeof(tspan[1]),1)
-    yout  = Array(typeof(y0),1)
-    dyout = Array(typeof(y0),1)
+    tout  = Array{typeof(tspan[1])}(undef,1)
+    yout  = Array{typeof(y0)}(undef,1)
+    dyout = Array{typeof(y0)}(undef,1)
     tout[1]  = tspan[1]
     yout[1]  = y0
     dyout[1] = dy0
@@ -451,7 +453,7 @@ function stepper(ord::Integer,
                              f_newton, # we want to find zeroes of this function
                              normy)     # the norm used to estimate error needs weights
 
-    alpha = Array(eltype(t),ord+1)
+    alpha = Array{eltype(t)}(undef,ord+1)
 
     for i = 1:ord
         alpha[i] = h_next/(t_next-t[end-i+1])
@@ -470,7 +472,7 @@ function stepper(ord::Integer,
     alpha[ord+1] = h_next/(t_next-t0)
 
     alpha0 = -sum(alpha[1:ord])
-    M      =  max(alpha[ord+1],abs(alpha[ord+1]+alphas-alpha0))
+    M      =  max(alpha[ord+1],abs.(alpha[ord+1]+alphas-alpha0))
     err::eltype(t) =  normy((yc-y0))*M
 
 
@@ -540,7 +542,7 @@ function newton_iteration(f,
     # after the first iteration the normy turned out to be very small,
     # terminate and return the first correction step
 
-    ep    = eps(eltype(abs(y0))) # this is the epsilon for type y0
+    ep    = eps(eltype(abs.(y0))) # this is the epsilon for type y0
 
     if norm1 < 100*ep*normy(y0)
         status=0
@@ -586,7 +588,7 @@ function dassl_norm(v, wt)
 end
 
 function dassl_weights(y,reltol,abstol)
-    reltol*abs(y).+abstol
+    @. reltol*abs(y)+abstol
 end
 
 # returns the value of the interpolation polynomial at the point x0
@@ -689,13 +691,13 @@ function numerical_jacobian(F,reltol,abstol,weights)
         # 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(max(abs(y),abs(h*dy),wt)*sqrt(ep))
+        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),n,n)
+        jac = Array{eltype(y)}(undef,n,n)
         for i=1:n
             jac[:,i]=(f(y+edelta[:,i])-f(y))/edelta[i,i]
         end

src/common.jl

@@ -1,4 +1,4 @@
-abstract type DASSLDAEAlgorithm <: AbstractODEAlgorithm end
+abstract type DASSLDAEAlgorithm <: DiffEqBase.AbstractDAEAlgorithm end
 struct dassl <: DASSLDAEAlgorithm
   maxorder
   factorize_jacobian
@@ -7,10 +7,12 @@ end
 dassl(;maxorder = 6,factorize_jacobian = true) = dassl(maxorder,factorize_jacobian)
 
 function solve(
-    prob::AbstractDAEProblem{uType,duType,tType,isinplace},
+    prob::DiffEqBase.AbstractDAEProblem{uType,duType,tupType,isinplace},
     alg::DASSLDAEAlgorithm,args...;timeseries_errors=true,
     abstol=1e-5,reltol=1e-3,dt = 1e-4, dtmin = 0.0, dtmax = Inf,
-    callback=nothing,kwargs...) where {uType,duType,tType,isinplace}
+    callback=nothing,kwargs...) where {uType,duType,tupType,isinplace}
+
+    tType = eltype(tupType)
 
     if callback != nothing || prob.callback != nothing
         error("DASSL is not compatible with callbacks.")
@@ -52,6 +54,6 @@ function solve(
     end
     =#
 
-    build_solution(prob,alg,ts,timeseries,
+    DiffEqBase.build_solution(prob,alg,ts,timeseries,
                       timeseries_errors = timeseries_errors)
 end

test/common.jl

@@ -1,4 +1,6 @@
 using DiffEqProblemLibrary, DiffEqBase, DASSL
+using DiffEqProblemLibrary.DAEProblemLibrary: importdaeproblems; importdaeproblems()
+using DiffEqProblemLibrary.DAEProblemLibrary: prob_dae_resrob
 
 prob = prob_dae_resrob
 

test/runtests.jl

@@ -1,14 +1,12 @@
-using DASSL
-using Base.Test
-using FactCheck
+using DASSL, Test
 
-facts("Testing maxorder") do
+@testset "Testing maxorder" begin
 
     F(t,y,dy)=(dy+y.^2)
     Fy(t,y,dy)=diagm(2y)
     Fdy(t,y,dy)=eye(length(y))
 
-    sol(t)=1./(1 .+ t)
+    sol(t)=1.0./(1 .+ t)
     tspan=[0.0, 1.0]
 
     atol = 1.0e-5
@@ -19,53 +17,30 @@ facts("Testing maxorder") do
     for order=1:6
         # scalar version
         (tn,yn,dyn)=DASSL.dasslSolve(F, sol(0.0), tspan, maxorder = order)
-        aerror = maximum(abs(yn-sol(tn)))
-        rerror = maximum(abs(yn-sol(tn))./abs(sol(tn)))
+        aerror = maximum(abs.(yn-sol(tn)))
+        rerror = maximum(abs.(yn-sol(tn))./abs.(sol(tn)))
         nsteps = length(tn)
 
-        @fact aerror --> less_than(2*nsteps*atol)
-        @fact rerror --> less_than(2*nsteps*rtol)
+        @test aerror < (2*nsteps*atol)
+        @test rerror < (2*nsteps*rtol)
 
         # vector version
         (tnV,ynV,dynV)=DASSL.dasslSolve(F,[sol(0.0)], tspan, maxorder = order)
 
-        @fact  vcat(ynV...) --> yn
-        @fact vcat(dynV...) --> dyn
+        @test  vcat(ynV...) == yn
+        @test  vcat(dynV...) == dyn
 
         # analytical jacobian version (vector)
         (tna,yna,dyna)=dasslSolve(F, [sol(0.0)], tspan, maxorder = order, Fy = Fy, Fdy = Fdy)
-        aerror = maximum(abs(map(first,yn)-sol(tn)))
-        rerror = maximum(abs(map(first,yn)-sol(tn))./abs(sol(tn)))
+        aerror = maximum(abs.(map(first,yn)-sol(tn)))
+        rerror = maximum(abs.(map(first,yn)-sol(tn))./abs.(sol(tn)))
         nsteps = length(tn)
 
-        @fact aerror --> less_than(2*nsteps*atol)
-        @fact rerror --> less_than(2*nsteps*rtol)
+        @test aerror < (2*nsteps*atol)
+        @test rerror < (2*nsteps*rtol)
     end
 end
 
-facts("Testing minimal error tolerances") do
-    eps=1e-6
-    # van der Pol equation
-    Fvdp(t,y,dy)=([dy[1]+y[2],
-                   eps*dy[2]-y[1]+y[2]*(1.0-y[1]^2)])
-    y0 = [2., 0]
-    tn = (0.,[2.,0.],[0.,-2/eps])
-
-    tol = 1e0
-    while true
-        tol/=10
-        sol=dasslIterator(Fvdp, y0, 1.0; reltol=tol, abstol=tol)
-
-        try
-            (tn,yn,dyn)=consume(sol)
-        catch
-            break
-        end
-
-    end
-    @fact tol --> roughly(1e-15)
-end
-
-facts("Testing common interface") do
+@testset "Testing common interface" begin
   include("common.jl")
 end