Fix BAOAB cache to work with scalar u

Author: jamesgardner1421

src/caches/dynamical_caches.jl

@@ -1,5 +1,5 @@
 
-struct BAOABConstantCache{uType,uEltypeNoUnits} <: StochasticDiffEqConstantCache
+mutable struct BAOABConstantCache{uType,uEltypeNoUnits} <: StochasticDiffEqConstantCache
   k::uType
   half::uEltypeNoUnits
   c1::uEltypeNoUnits

src/perform_step/dynamical.jl

@@ -17,7 +17,7 @@ function initialize!(integrator, cache::BAOABConstantCache)
   u1 = integrator.uprev.x[2]
 
   verify_f2(integrator.f.f2, du1, u1, p, t, integrator, cache)
-  cache.k .= integrator.f.f1(du1,u1,p,t)
+  cache.k = integrator.f.f1(du1,u1,p,t)
 end
 
 function initialize!(integrator, cache::BAOABCache)
@@ -31,12 +31,12 @@ end
 
 @muladd function perform_step!(integrator,cache::BAOABConstantCache)
   @unpack t,dt,sqdt,uprev,u,p,W,f = integrator
-  @unpack k, half, c1, c2 = cache
+  @unpack half, c1, c2 = cache
   du1 = uprev.x[1]
   u1 = uprev.x[2]
 
   # B
-  du2 = du1 + half*dt*k
+  du2 = du1 + half*dt*cache.k
 
   # A
   u2 = u1 + half*dt*du2
@@ -49,8 +49,8 @@ end
   u = u2 + half*dt*du3
 
   # B
-  k .= f.f1(du3,u,p,t+dt)
-  du = du3 + half*dt*k
+  cache.k = f.f1(du3,u,p,t+dt)
+  du = du3 + half*dt*cache.k
 
   integrator.u = ArrayPartition((du, u))
 end

test/sde/sde_dynamical.jl

@@ -1,29 +1,47 @@
 using StochasticDiffEq, DiffEqNoiseProcess, Test, DiffEqDevTools, Random
 Random.seed!(1)
 
-u0 = zeros(2)
-v0 = ones(2)
-γ = 1
-
 f1_harmonic(v,u,p,t) = -u
 f2_harmonic(v,u,p,t) = v
 g(u,p,t) = 1
+γ = 1
+
+@testset "Scalar u" begin
+    u0 = 0
+    v0 = 1
+
+    ff_harmonic = DynamicalSDEFunction(f1_harmonic,f2_harmonic,g)
+    prob1 = DynamicalSDEProblem(ff_harmonic,g,v0,u0,(0.0,5.0))
+
+    dts = (1/2) .^ (8:-1:4)
+
+    # Can't use NoiseGrid as noise is not generated with the correct size in convergence.jl. We require noise with shape of v.
+    sim1  = analyticless_test_convergence(dts,prob1,BAOAB(gamma=γ),(1/2)^10;trajectories=Int(2e2),use_noise_grid=false)
+    display(sim1.𝒪est)
+    @test abs(sim1.𝒪est[:weak_final]-1) < 0.3
+end
+
+@testset "Vector u" begin
+
+    u0 = zeros(2)
+    v0 = ones(2)
 
-f1_harmonic_iip(dv,v,u,p,t) = dv .= f1_harmonic(v,u,p,t)
-f2_harmonic_iip(du,v,u,p,t) = du .= f2_harmonic(v,u,p,t)
-g_iip(du,u,p,t) = du .= g(u,p,t)
+    f1_harmonic_iip(dv,v,u,p,t) = dv .= f1_harmonic(v,u,p,t)
+    f2_harmonic_iip(du,v,u,p,t) = du .= f2_harmonic(v,u,p,t)
+    g_iip(du,u,p,t) = du .= g(u,p,t)
 
-ff_harmonic = DynamicalSDEFunction(f1_harmonic,f2_harmonic,g)
-prob1 = DynamicalSDEProblem(ff_harmonic,g,v0,u0,(0.0,5.0))
-sol1 = solve(prob1,BAOAB(gamma=γ);dt=1/10,save_noise=true)
+    ff_harmonic = DynamicalSDEFunction(f1_harmonic,f2_harmonic,g)
+    prob1 = DynamicalSDEProblem(ff_harmonic,g,v0,u0,(0.0,5.0))
+    sol1 = solve(prob1,BAOAB(gamma=γ);dt=1/10,save_noise=true)
 
-prob2 = DynamicalSDEProblem(f1_harmonic_iip,f2_harmonic_iip,g_iip,v0,u0,(0.0,5.0); noise=NoiseWrapper(sol1.W))
-sol2 = solve(prob2,BAOAB(gamma=γ);dt=1/10)
+    prob2 = DynamicalSDEProblem(f1_harmonic_iip,f2_harmonic_iip,g_iip,v0,u0,(0.0,5.0); noise=NoiseWrapper(sol1.W))
+    sol2 = solve(prob2,BAOAB(gamma=γ);dt=1/10)
 
-@test sol1[:] ≈ sol2[:]
+    @test sol1[:] ≈ sol2[:]
 
-dts = (1/2) .^ (8:-1:4)
+    dts = (1/2) .^ (8:-1:4)
 
-# Can't use NoiseGrid as noise is not generated with the correct size in convergence.jl. We require noise with shape of v.
-sim1  = analyticless_test_convergence(dts,prob1,BAOAB(gamma=γ),(1/2)^10;trajectories=Int(1e2),use_noise_grid=false)
-@test abs(sim1.𝒪est[:weak_final]-1) < 0.3
+    # Can't use NoiseGrid as noise is not generated with the correct size in convergence.jl. We require noise with shape of v.
+    sim1  = analyticless_test_convergence(dts,prob1,BAOAB(gamma=γ),(1/2)^10;trajectories=Int(1e2),use_noise_grid=false)
+    @test abs(sim1.𝒪est[:weak_final]-1) < 0.3
+end