Merge remote-tracking branch 'origin/master'

src/StochasticDiffEq.jl

@@ -43,6 +43,7 @@ module StochasticDiffEq
   include("algorithms.jl")
   include("options_type.jl")
   include("derivative_wrappers.jl")
+  include("derivative_utils.jl")
   include("interp_func.jl")
   include("caches/cache_types.jl")
   include("caches/basic_method_caches.jl")

src/derivative_utils.jl

@@ -0,0 +1,47 @@
+function calc_J!(integrator, cache, is_compos)
+    @unpack t,dt,uprev,u,f,p = integrator
+    @unpack du1,uf,J,jac_config = cache
+    if has_jac(f)
+      f(Val{:jac}, J, uprev, p, t)
+    else
+      uf.t = t
+      uf.p = p
+      jacobian!(J, uf, uprev, du1, integrator, jac_config)
+      if is_compos
+        integrator.eigen_est = norm(J, Inf)
+      end
+    end
+end
+
+function calc_W!(integrator, cache, γdt, repeat_step)
+  @inbounds begin
+    @unpack t,dt,uprev,u,f,p = integrator
+    @unpack J,W,jac_config = cache
+    is_compos = typeof(integrator.alg) <: StochasticCompositeAlgorithm
+    mass_matrix = integrator.sol.prob.mass_matrix
+
+    new_W = true
+    if has_invW(f)
+      # skip calculation of inv(W) if step is repeated
+      !repeat_step && f(Val{:invW},W,uprev,p,γdt,t) # W == inverse W
+      is_compos && calc_J!(integrator, cache, true)
+    else
+      # skip calculation of J if step is repeated
+      if repeat_step || (!integrator.last_stepfail && cache.newton_iters == 1 && cache.ηold < integrator.alg.new_jac_conv_bound)
+        new_jac = false
+      else # Compute a new Jacobian
+        new_jac = true
+        calc_J!(integrator, cache, is_compos)
+      end
+      # skip calculation of W if step is repeated
+      if !repeat_step && (integrator.iter < 1 || new_jac || abs(dt - (t-integrator.tprev)) > 100eps(typeof(integrator.t)))
+        for j in 1:length(u), i in 1:length(u)
+            @inbounds W[i,j] = mass_matrix[i,j]-γdt*J[i,j]
+        end
+      else
+        new_W = false
+      end
+    end
+    return new_W
+  end
+end

src/integrators/type.jl

@@ -1,4 +1,4 @@
-mutable struct SDEIntegrator{T1,uType,uEltype,tType,P,tTypeNoUnits,uEltypeNoUnits,randType,rateType,solType,cacheType,progType,F4,F5,OType,noiseType} <: AbstractSDEIntegrator
+mutable struct SDEIntegrator{T1,uType,uEltype,tType,P,eigenType,tTypeNoUnits,uEltypeNoUnits,randType,rateType,solType,cacheType,progType,F4,F5,OType,noiseType} <: AbstractSDEIntegrator
   f::F4
   g::F5
   noise::noiseType
@@ -31,6 +31,7 @@ mutable struct SDEIntegrator{T1,uType,uEltype,tType,P,tTypeNoUnits,uEltypeNoUnit
   iter::Int
   success_iter::Int
   prog::progType
+  eigen_est::eigenType
   EEst::tTypeNoUnits
   q::tTypeNoUnits
   qold::tTypeNoUnits

src/perform_step/implicit_split_step.jl

@@ -130,6 +130,8 @@ end
   mass_matrix = integrator.sol.prob.mass_matrix
   theta = integrator.alg.theta
 
+  repeat_step = false
+
   if integrator.success_iter > 0 && !integrator.u_modified && integrator.alg.extrapolant == :interpolant
     current_extrapolant!(u,t+dt,integrator)
   elseif integrator.alg.extrapolant == :linear
@@ -138,30 +140,7 @@ end
     copy!(u,uprev)
   end
 
-  if has_invW(f)
-    f(Val{:invW},W,uprev,p,dt,t) # W == inverse W
-  else
-    if !integrator.last_stepfail && cache.newton_iters == 1 && cache.ηold < integrator.alg.new_jac_conv_bound
-      new_jac = false
-    else # Compute a new Jacobian
-      new_jac = true
-      if has_jac(f)
-        f(Val{:jac},J,uprev,p,t)
-      else
-        uf.t = t
-        jacobian!(J, uf, uprev, du1, integrator, jac_config)
-      end
-    end
-    if integrator.iter < 1 || new_jac || abs(dt - (t-integrator.tprev)) > 100eps()
-      new_W = true
-      thetadt = dt*theta
-      for j in 1:length(u), i in 1:length(u)
-          @inbounds W[i,j] = mass_matrix[i,j]-thetadt*J[i,j]
-      end
-    else
-      new_W = false
-    end
-  end
+  new_W = calc_W!(integrator, cache, dt, repeat_step)
 
   integrator.f(tmp,uprev,p,t)
 

src/perform_step/kencarp.jl

@@ -265,33 +265,7 @@ end
 
   γdt = γ*dt
 
-  if has_invW(f)
-    # skip calculation of inv(W) if step is repeated
-    !repeat_step && f(Val{:invW},W,uprev,p,γdt,t) # W == inverse W
-  else
-    # skip calculation of J if step is repeated
-    if repeat_step || (!integrator.last_stepfail && cache.newton_iters == 1 && cache.ηold < integrator.alg.new_jac_conv_bound)
-      new_jac = false
-    else # Compute a new Jacobian
-      new_jac = true
-      if has_jac(f)
-        f(Val{:jac},J,uprev,p,t)
-      else
-        uf.t = t
-        jacobian!(J, uf, uprev, du1, integrator, jac_config)
-      end
-    end
-    # skip calculation of W if step is repeated
-    if !repeat_step && (integrator.iter < 1 || new_jac || abs(dt - (t-integrator.tprev)) > 100eps(typeof(integrator.t)))
-      new_W = true
-      mass_matrix = integrator.sol.prob.mass_matrix
-      for j in 1:length(u), i in 1:length(u)
-          @inbounds W[i,j] = mass_matrix[i,j]-γdt*J[i,j]
-      end
-    else
-      new_W = false
-    end
-  end
+  new_W = calc_W!(integrator, cache, γdt, repeat_step)
 
   if !repeat_step && !integrator.last_stepfail
     f(z₁, integrator.uprev, p, integrator.t)

src/perform_step/sdirk.jl

@@ -141,6 +141,8 @@ end
   mass_matrix = integrator.sol.prob.mass_matrix
   theta = integrator.alg.theta
 
+  repeat_step = false
+
   if integrator.success_iter > 0 && !integrator.u_modified && integrator.alg.extrapolant == :interpolant
     current_extrapolant!(u,t+dt,integrator)
   elseif integrator.alg.extrapolant == :linear
@@ -149,30 +151,7 @@ end
     copy!(u,uprev)
   end
 
-  if has_invW(f)
-    f(Val{:invW},W,uprev,p,dt,t) # W == inverse W
-  else
-    if !integrator.last_stepfail && cache.newton_iters == 1 && cache.ηold < integrator.alg.new_jac_conv_bound
-      new_jac = false
-    else # Compute a new Jacobian
-      new_jac = true
-      if has_jac(f)
-        f(Val{:jac},J,uprev,p,t)
-      else
-        uf.t = t
-        jacobian!(J, uf, uprev, du1, integrator, jac_config)
-      end
-    end
-    if integrator.iter < 1 || new_jac || abs(dt - (t-integrator.tprev)) > 100eps()
-      new_W = true
-      thetadt = dt*theta
-      for j in 1:length(u), i in 1:length(u)
-          @inbounds W[i,j] = mass_matrix[i,j]-thetadt*J[i,j]
-      end
-    else
-      new_W = false
-    end
-  end
+  new_W = calc_W!(integrator, cache, dt, repeat_step)
 
   ##############################################################################
 

src/solve.jl

@@ -283,6 +283,7 @@ function init(
     end
   end
 
+  eigen_est = 1/oneunit(tType)
   EEst = tTypeNoUnits(1)
   q = tTypeNoUnits(1)
   just_hit_tstop = false
@@ -320,7 +321,7 @@ function init(
   end
 
   integrator =    SDEIntegrator{typeof(alg),uType,uBottomEltype,tType,typeof(p),
-                  tTypeNoUnits,
+                  typeof(eigen_est),tTypeNoUnits,
                   uEltypeNoUnits,typeof(W),rateType,typeof(sol),typeof(cache),
                   typeof(prog),FType,GType,typeof(opts),typeof(noise)}(
                   f,g,noise,uprev,tprev,t,u,p,tType(dt),tType(dt),tType(dt),dtcache,T,tdir,
@@ -330,7 +331,7 @@ function init(
                   saveiter,
                   alg,sol,
                   cache,tType(dt),W,
-                  opts,iter,success_iter,prog,EEst,q,
+                  opts,iter,success_iter,prog,eigen_est,EEst,q,
                   tTypeNoUnits(qoldinit),q11)
 
   if initialize_integrator