full scalar noise

Author: ChrisRackauckas

src/caches/rossler_caches.jl

@@ -20,7 +20,11 @@ du_cache(c::SRA1Cache) = (c.chi2,c.E₁,c.E₂,c.gt,c.k₁,c.k₂,c.gpdt)
 user_cache(c::SRA1Cache) = (c.u,c.uprev,c.tmp,c.tmp1)
 
 function alg_cache(alg::SRA1,prob,u,ΔW,ΔZ,rate_prototype,noise_rate_prototype,uEltypeNoUnits,tTypeNoUnits,uprev,f,t,::Type{Val{true}})
-  chi2 = similar(ΔW)
+  if typeof(ΔW) <: Union{SArray,Number}
+    chi2 = copy(ΔW)
+  else
+    chi2 = similar(ΔW)
+  end
   tmp1 = zeros(u)
   E₁ = zeros(rate_prototype); gt = zeros(noise_rate_prototype); gpdt = zeros(noise_rate_prototype)
   E₂ = zeros(rate_prototype); k₁ = zeros(rate_prototype); k₂ = zeros(rate_prototype)
@@ -81,7 +85,12 @@ function alg_cache(alg::SRA,prob,u,ΔW,ΔZ,rate_prototype,noise_rate_prototype,u
     push!(H0,zeros(u))
   end
   A0temp = zeros(rate_prototype); B0temp = zeros(rate_prototype)
-  ftmp = zeros(rate_prototype); gtmp = zeros(noise_rate_prototype); chi2 = similar(ΔW)
+  ftmp = zeros(rate_prototype); gtmp = zeros(noise_rate_prototype);
+  if typeof(ΔW) <: Union{SArray,Number}
+    chi2 = copy(ΔW)
+  else
+    chi2 = similar(ΔW)
+  end
   atemp = zeros(rate_prototype); btemp = zeros(rate_prototype); E₂ = zeros(rate_prototype); E₁temp = zeros(rate_prototype)
   E₁ = zeros(rate_prototype)
   tmp = zeros(u)

src/integrators/sra.jl

@@ -68,7 +68,11 @@ end
   integrator.g(t,uprev,gt)
   integrator.g(t+dt,uprev,gpdt)
   integrator.f(t,uprev,k₁); k₁*=dt
-  @. chi2 = (W.dW + W.dZ/sqrt(3))/2 #I_(1,0)/h
+  if typeof(W.dW) <: Union{SArray,Number}
+    chi2 = @. (W.dW + W.dZ/sqrt(3))/2 #I_(1,0)/h
+  else
+    @. chi2 = (W.dW + W.dZ/sqrt(3))/2 #I_(1,0)/h
+  end
 
   if is_diagonal_noise(integrator.sol.prob)
     @. E₁ = chi2*gpdt
@@ -155,7 +159,13 @@ end
   @unpack t,dt,uprev,u,W = integrator
   @unpack H0,A0temp,B0temp,ftmp,gtmp,chi2,atemp,btemp,E₁,E₁temp,E₂,tmp = cache
   @unpack c₀,c₁,A₀,B₀,α,β₁,β₂,stages = cache.tab
-  @. chi2 = (W.dW + W.dZ/sqrt(3))/2 #I_(1,0)/h
+
+  if typeof(W.dW) <: Union{SArray,Number}
+    chi2 = @. (W.dW + W.dZ/sqrt(3))/2 #I_(1,0)/h
+  else
+    @. chi2 = (W.dW + W.dZ/sqrt(3))/2 #I_(1,0)/h
+  end
+  
   for i in 1:stages
     fill!(H0[i],zero(eltype(integrator.u)))
   end

src/integrators/sri.jl

@@ -56,11 +56,17 @@ end
   @unpack c₀,c₁,A₀,A₁,B₀,B₁,α,β₁,β₂,β₃,β₄,stages,error_terms = cache.tab
   @unpack H0,H1,A0temp,A1temp,B0temp,B1temp,A0temp2,A1temp2,B0temp2,B1temp2,atemp,btemp,E₁,E₂,E₁temp,ftemp,gtemp,chi1,chi2,chi3,tmp = cache
   @unpack t,dt,uprev,u,W = integrator
-  @tight_loop_macros for i in eachindex(u)
-    @inbounds chi1[i] = .5*(W.dW[i].^2 - dt)/integrator.sqdt #I_(1,1)/sqrt(h)
-    @inbounds chi2[i] = .5*(W.dW[i] + W.dZ[i]/sqrt(3)) #I_(1,0)/h
-    @inbounds chi3[i] = 1/6 * (W.dW[i].^3 - 3*W.dW[i]*dt)/dt #I_(1,1,1)/h
+
+  if typeof(W.dW) <: Union{SArray,Number}
+    chi1 = @. (W.dW.^2 - dt)/2integrator.sqdt #I_(1,1)/sqrt(h)
+    chi2 = @. (W.dW + W.dZ/sqrt(3))/2 #I_(1,0)/h
+    chi3 = @. (W.dW.^3 - 3W.dW*dt)/6dt #I_(1,1,1)/h
+  else
+    @. chi1 = (W.dW.^2 - dt)/2integrator.sqdt #I_(1,1)/sqrt(h)
+    @. chi2 = (W.dW + W.dZ/sqrt(3))/2 #I_(1,0)/h
+    @. chi3 = (W.dW.^3 - 3W.dW*dt)/6dt #I_(1,1,1)/h
   end
+
   for i=1:stages
     fill!(H0[i],zero(eltype(integrator.u)))
     fill!(H1[i],zero(eltype(integrator.u)))

test/scalar_noise.jl

@@ -3,8 +3,21 @@ using StochasticDiffEq, DiffEqNoiseProcess
 f(t,u,du) = (du .= u)
 g(t,u,du) = (du .= u)
 u0 = rand(4,2)
+
 W = WienerProcess(0.0,0.0,0.0)
 prob = SDEProblem(f,g,u0,(0.0,1.0),noise=W)
 sol = solve(prob,SRIW1())
 
-similar(0.0)
+W = WienerProcess(0.0,0.0,0.0)
+prob = SDEProblem(f,g,u0,(0.0,1.0),noise=W)
+sol = solve(prob,SRI())
+
+g(t,u,du) = (du .= 1)
+
+W = WienerProcess(0.0,0.0,0.0)
+prob = SDEProblem(f,g,u0,(0.0,1.0),noise=W)
+sol = solve(prob,SRA1())
+
+W = WienerProcess(0.0,0.0,0.0)
+prob = SDEProblem(f,g,u0,(0.0,1.0),noise=W)
+sol = solve(prob,SRA())