fix out of place init dt for arrays

SciML · Jan 16, 2017 · 8769c76 · 8769c76
1 parent e6a1cca
commit 8769c76
Show file tree

Hide file tree

Showing 3 changed files with 10 additions and 6 deletions.
diff --git a/src/initdt.jl b/src/initdt.jl
@@ -1,5 +1,6 @@
 
-function ode_determine_initdt{uType,tType,uEltypeNoUnits}(u0::uType,t::tType,tdir,dtmax,abstol,reltol::uEltypeNoUnits,internalnorm,f,order)
+function ode_determine_initdt{tType,uType,F}(u0,t::tType,tdir,dtmax,abstol,reltol,internalnorm,prob::AbstractODEProblem{uType,tType,true,F},order)
+  f = prob.f
   f₀ = similar(u0./t); f₁ = similar(u0./t); u₁ = similar(u0)
   sk = abstol+abs.(u0)*reltol
   d₀ = internalnorm(u0./sk)
@@ -28,11 +29,12 @@ function ode_determine_initdt{uType,tType,uEltypeNoUnits}(u0::uType,t::tType,tdi
   dt = tdir*min(100dt₀,dt₁)
 end
 
-function ode_determine_initdt{uType<:Number,tType,uEltypeNoUnits}(u0::uType,t::tType,tdir,dtmax,abstol,reltol::uEltypeNoUnits,internalnorm,f,order)
+function ode_determine_initdt{uType,tType,F}(u0::uType,t,tdir,dtmax,abstol,reltol,internalnorm,prob::AbstractODEProblem{uType,tType,false,F},order)
+  f = prob.f
   sk = abstol+abs(u0)*reltol
-  d₀ = abs(u0/sk)
+  d₀ = internalnorm(u0/sk)
   f₀ = f(t,u0)
-  d₁ = abs(f₀/sk)
+  d₁ = internalnorm(f₀/sk)
   T0 = typeof(d₀)
   T1 = typeof(d₁)
   if d₀ < T0(1//10^(5)) || d₁ < T1(1//10^(5))
@@ -43,7 +45,7 @@ function ode_determine_initdt{uType<:Number,tType,uEltypeNoUnits}(u0::uType,t::t
   dt₀ = min(dt₀,tdir*dtmax)
   u₁ = u0 + tdir*dt₀*f₀
   f₁ = f(t+tdir*dt₀,u₁)
-  d₂ = abs((f₁-f₀)./(abstol+abs(u0)*reltol))/dt₀*tType(1)
+  d₂ = internalnorm((f₁-f₀)./(abstol+abs(u0)*reltol))/dt₀*tType(1)
   if max(d₁,d₂) <= T1(1//10^(15))
     dt₁ = max(tType(1//10^(6)),dt₀*1//10^(3))
   else

diff --git a/src/solve.jl b/src/solve.jl
@@ -85,7 +85,7 @@ function init{uType,tType,isinplace,algType<:OrdinaryDiffEqAlgorithm,F,recompile
   tTypeNoUnits   = typeof(recursive_one(t))
 
   if dt == zero(dt) && adaptive
-    dt = tType(ode_determine_initdt(u0,t,tdir,dtmax,uEltype(abstol),uEltypeNoUnits(reltol),internalnorm,f,order))
+    dt = tType(ode_determine_initdt(u0,t,tdir,dtmax,uEltype(abstol),uEltypeNoUnits(reltol),internalnorm,prob,order))
   end
 
   if sign(dt)!=tdir && dt!=tType(0)

diff --git a/test/ode/ode_inplace_tests.jl b/test/ode/ode_inplace_tests.jl
@@ -28,3 +28,5 @@ alloc1 = @allocated sol = solve(prob_ode_large2Dlinear,Euler(),dt=1//2^(8),save_
 alloc2 = @allocated sol2 = solve(prob_ode_large2Dlinear,Euler(),sol[:],sol.t,sol.k;dt=1//2^(8),save_timeseries=true)
 
 @test alloc2 <= alloc1
+
+sol =solve(prob,Tsit5())