Merge pull request #483 from SciML/ISSEM

adaptive stepping `ISSEM` fixes

Author: ChrisRackauckas

src/caches/implicit_split_step_caches.jl

@@ -6,33 +6,38 @@
   gtmp2::rateType
   nlsolver::N
   dW_cache::randType
+  k::uType
+  dz::uType
 end
 
-function alg_cache(alg::ISSEM,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,jump_rate_prototype,
-                   ::Type{uEltypeNoUnits},::Type{uBottomEltypeNoUnits},::Type{tTypeNoUnits},uprev,f,t,dt,::Type{Val{true}}) where {uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits}
-  γ, c = alg.theta,zero(t)
-  nlsolver = OrdinaryDiffEq.build_nlsolver(alg,u,uprev,p,t,dt,f,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,γ,c,Val(true))
+function alg_cache(alg::ISSEM, prob, u, ΔW, ΔZ, p, rate_prototype, noise_rate_prototype, jump_rate_prototype,
+  ::Type{uEltypeNoUnits}, ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, f, t, dt, ::Type{Val{true}}) where {uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits}
+  γ, c = alg.theta, zero(t)
+  nlsolver = OrdinaryDiffEq.build_nlsolver(alg, u, uprev, p, t, dt, f, rate_prototype, uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits, γ, c, Val(true))
   fsalfirst = zero(rate_prototype)
   gtmp = zero(noise_rate_prototype)
   if is_diagonal_noise(prob)
-    gtmp2 = gtmp
+    gtmp2 = copy(gtmp)
     dW_cache = nothing
   else
     gtmp2 = zero(rate_prototype)
     dW_cache = zero(ΔW)
   end
 
-  ISSEMCache(u,uprev,fsalfirst,gtmp,gtmp2,nlsolver,dW_cache)
+  k = zero(u)
+  dz = zero(u)
+
+  ISSEMCache(u, uprev, fsalfirst, gtmp, gtmp2, nlsolver, dW_cache, k, dz)
 end
 
 mutable struct ISSEMConstantCache{N} <: StochasticDiffEqConstantCache
   nlsolver::N
 end
 
-function alg_cache(alg::ISSEM,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,jump_rate_prototype,
-                   ::Type{uEltypeNoUnits},::Type{uBottomEltypeNoUnits},::Type{tTypeNoUnits},uprev,f,t,dt,::Type{Val{false}}) where {uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits}
-  γ, c = alg.theta,zero(t)
-  nlsolver = OrdinaryDiffEq.build_nlsolver(alg,u,uprev,p,t,dt,f,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,γ,c,Val(false))
+function alg_cache(alg::ISSEM, prob, u, ΔW, ΔZ, p, rate_prototype, noise_rate_prototype, jump_rate_prototype,
+  ::Type{uEltypeNoUnits}, ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, f, t, dt, ::Type{Val{false}}) where {uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits}
+  γ, c = alg.theta, zero(t)
+  nlsolver = OrdinaryDiffEq.build_nlsolver(alg, u, uprev, p, t, dt, f, rate_prototype, uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits, γ, c, Val(false))
   ISSEMConstantCache(nlsolver)
 end
 
@@ -45,35 +50,40 @@ end
   gtmp3::noiseRateType
   nlsolver::N
   dW_cache::randType
+  k::uType
+  dz::uType
 end
 
-function alg_cache(alg::ISSEulerHeun,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,jump_rate_prototype,
-                   ::Type{uEltypeNoUnits},::Type{uBottomEltypeNoUnits},::Type{tTypeNoUnits},uprev,f,t,dt,::Type{Val{true}}) where {uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits}
-  γ, c = alg.theta,zero(t)
-  nlsolver = OrdinaryDiffEq.build_nlsolver(alg,u,uprev,p,t,dt,f,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,γ,c,Val(true))
+function alg_cache(alg::ISSEulerHeun, prob, u, ΔW, ΔZ, p, rate_prototype, noise_rate_prototype, jump_rate_prototype,
+  ::Type{uEltypeNoUnits}, ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, f, t, dt, ::Type{Val{true}}) where {uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits}
+  γ, c = alg.theta, zero(t)
+  nlsolver = OrdinaryDiffEq.build_nlsolver(alg, u, uprev, p, t, dt, f, rate_prototype, uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits, γ, c, Val(true))
   fsalfirst = zero(rate_prototype)
 
   gtmp = zero(noise_rate_prototype)
   gtmp2 = zero(rate_prototype)
 
   if is_diagonal_noise(prob)
-      gtmp3 = gtmp2
-      dW_cache = nothing
+    gtmp3 = copy(gtmp2)
+    dW_cache = nothing
   else
-      gtmp3 = zero(noise_rate_prototype)
-      dW_cache = zero(ΔW)
+    gtmp3 = zero(noise_rate_prototype)
+    dW_cache = zero(ΔW)
   end
 
-  ISSEulerHeunCache(u,uprev,fsalfirst,gtmp,gtmp2,gtmp3,nlsolver,dW_cache)
+  k = zero(u)
+  dz = zero(u)
+
+  ISSEulerHeunCache(u, uprev, fsalfirst, gtmp, gtmp2, gtmp3, nlsolver, dW_cache, k, dz)
 end
 
 mutable struct ISSEulerHeunConstantCache{N} <: StochasticDiffEqConstantCache
   nlsolver::N
 end
 
-function alg_cache(alg::ISSEulerHeun,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,jump_rate_prototype,
-                   ::Type{uEltypeNoUnits},::Type{uBottomEltypeNoUnits},::Type{tTypeNoUnits},uprev,f,t,dt,::Type{Val{false}}) where {uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits}
-  γ, c = alg.theta,zero(t)
-  nlsolver = OrdinaryDiffEq.build_nlsolver(alg,u,uprev,p,t,dt,f,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,γ,c,Val(false))
+function alg_cache(alg::ISSEulerHeun, prob, u, ΔW, ΔZ, p, rate_prototype, noise_rate_prototype, jump_rate_prototype,
+  ::Type{uEltypeNoUnits}, ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, f, t, dt, ::Type{Val{false}}) where {uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits}
+  γ, c = alg.theta, zero(t)
+  nlsolver = OrdinaryDiffEq.build_nlsolver(alg, u, uprev, p, t, dt, f, rate_prototype, uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits, γ, c, Val(false))
   ISSEulerHeunConstantCache(nlsolver)
 end