SDE reversal

Author: frankschae

src/integrators/integrator_utils.jl

@@ -41,7 +41,7 @@ end
       modify_dtnew_for_tstops!(integrator)
       reject_step!(integrator)
       integrator.dt = integrator.dtnew
-      integrator.sqdt = sqrt(abs(integrator.dt))
+      integrator.sqdt = integrator.tdir * sqrt(abs(integrator.dt))
     end
   end
 
@@ -292,7 +292,7 @@ end
 
   # Allow RSWM1 on Wiener Process to change dt
   !isnothing(integrator.W) && (integrator.dt = integrator.W.dt)
-  integrator.sqdt = @fastmath sqrt(abs(integrator.dt)) # It can change dt, like in RSwM1
+  integrator.sqdt = @fastmath integrator.tdir*sqrt(abs(integrator.dt)) # It can change dt, like in RSwM1
 end
 
 @inline function handle_tstop!(integrator)

src/perform_step/low_order.jl

@@ -174,7 +174,7 @@ end
   if alg_interpretation(integrator.alg) == :Ito
     utilde =  K + L*integrator.sqdt
     ggprime = (integrator.g(utilde,p,t).-L)./(integrator.sqdt)
-    mil_correction = ggprime.*(W.dW.^2 .- dt)./2
+    mil_correction = ggprime.*(W.dW.^2 .- abs(dt))./2
   elseif alg_interpretation(integrator.alg) == :Stratonovich
     utilde = uprev + L*integrator.sqdt
     ggprime = (integrator.g(utilde,p,t).-L)./(integrator.sqdt)
@@ -183,6 +183,7 @@ end
     error("Alg interpretation invalid. Use either :Ito or :Stratonovich")
   end
   u = K+L.*W.dW+mil_correction
+
   if integrator.opts.adaptive
     du2 = integrator.f(K,p,t+dt)
     Ed = dt*(du2 - du1)/2
@@ -196,31 +197,6 @@ end
   integrator.u = u
 end
 
-#=
-@muladd function perform_step!(integrator,cache::RKMilCache)
-  @unpack du1,du2,K,tmp,L = cache
-  @unpack t,dt,uprev,u,W,p,f = integrator
-  integrator.f(du1,uprev,p,t)
-  integrator.g(L,uprev,p,t)
-  @.. K = uprev + dt * du1
-  @.. tmp = K + integrator.sqdt * L
-  integrator.g(du2,tmp,p,t)
-  if alg_interpretation(integrator.alg) == :Ito
-    @.. tmp = (du2-L)/(2integrator.sqdt)*(W.dW^2 - dt)
-  elseif alg_interpretation(integrator.alg) == :Stratonovich
-    @.. tmp = (du2-L)/(2integrator.sqdt)*(W.dW^2)
-  else
-    error("Alg interpretation invalid. Use either :Ito or :Stratonovich")
-  end
-  @.. u = K+L*W.dW + tmp
-  if integrator.opts.adaptive
-    @.. tmp = (tmp)/(integrator.opts.abstol + max(integrator.opts.internalnorm(uprev,t),integrator.opts.internalnorm(u,t))*integrator.opts.reltol)
-    integrator.EEst = integrator.opts.internalnorm(tmp,t)
-  end
-  integrator.u = u
-end
-=#
-
 @muladd function perform_step!(integrator,cache::RKMilCache)
   @unpack du1,du2,K,tmp,L = cache
   @unpack t,dt,uprev,u,W,p,f = integrator
@@ -231,7 +207,7 @@ end
   if alg_interpretation(integrator.alg) == :Ito
     @.. tmp = K + integrator.sqdt * L
     integrator.g(du2,tmp,p,t)
-    @.. tmp = (du2-L)/(2integrator.sqdt)*(W.dW.^2 - dt)
+    @.. tmp = (du2-L)/(2integrator.sqdt)*(W.dW.^2 - abs(dt))
   elseif alg_interpretation(integrator.alg) == :Stratonovich
     @.. tmp = uprev + integrator.sqdt * L
     integrator.g(du2,tmp,p,t)