Skip to content

Commit

Permalink
Merge pull request #559 from Lightup1/baoab-dt
Browse files Browse the repository at this point in the history
fix dt scaling when scale_noise==false
  • Loading branch information
ChrisRackauckas authored Dec 16, 2023
2 parents b6d2ec1 + 24026eb commit 2297905
Show file tree
Hide file tree
Showing 3 changed files with 44 additions and 4 deletions.
4 changes: 2 additions & 2 deletions src/caches/dynamical_caches.jl
Original file line number Diff line number Diff line change
Expand Up @@ -20,7 +20,7 @@ end
function alg_cache(alg::BAOAB,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}
k = zero(rate_prototype.x[1])
c1 = exp(-alg.gamma*dt)
c2 = sqrt(1 - alg.scale_noise*c1^2) # if scale_noise == false, c2 = 1
c2 = alg.scale_noise ? sqrt((1 - c1^2)/abs(dt)) : 1 # if scale_noise == false, c2 = 1
BAOABConstantCache(k, uEltypeNoUnits(1//2), uEltypeNoUnits(c1), uEltypeNoUnits(c2))
end

Expand All @@ -34,7 +34,7 @@ function alg_cache(alg::BAOAB,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototy

half = uEltypeNoUnits(1//2)
c1 = exp(-alg.gamma*dt)
c2 = sqrt(1 - alg.scale_noise*c1^2) # if scale_noise == false, c2 = 1
c2 = alg.scale_noise ? sqrt((1 - c1^2)/abs(dt)) : 1 # if scale_noise == false, c2 = 1

tmp = zero(u)

Expand Down
4 changes: 2 additions & 2 deletions src/perform_step/dynamical.jl
Original file line number Diff line number Diff line change
Expand Up @@ -42,7 +42,7 @@ end
u2 = u1 + half*dt*du2

# O
noise = integrator.g(u2,p,t+dt*half).*W.dW / sqdt
noise = integrator.g(u2,p,t+dt*half).*W.dW
du3 = c1*du2 + c2*noise

# A
Expand All @@ -69,7 +69,7 @@ end

# O
integrator.g(gtmp,utmp,p,t+dt*half)
@.. noise = gtmp*W.dW / sqdt
@.. noise = gtmp*W.dW
@.. dutmp = c1*dutmp + c2*noise

# A
Expand Down
40 changes: 40 additions & 0 deletions test/sde/sde_dynamical.jl
Original file line number Diff line number Diff line change
Expand Up @@ -45,3 +45,43 @@ end
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

@testset "Scalar u, scale_noise=false" begin
u0 = 0
v0 = 1

ff_harmonic = DynamicalSDEFunction(f1_harmonic,f2_harmonic,g)
prob1 = DynamicalSDEProblem(ff_harmonic,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=γ,scale_noise=false),(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, scale_noise=false" 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)

ff_harmonic = DynamicalSDEFunction(f1_harmonic,f2_harmonic,g)
prob1 = DynamicalSDEProblem(ff_harmonic,v0,u0,(0.0,5.0))
sol1 = solve(prob1,BAOAB(gamma=γ,scale_noise=false);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=γ,scale_noise=false);dt=1/10)

@test sol1[:] sol2[:]

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=γ,scale_noise=false),(1/2)^10;trajectories=Int(1e2),use_noise_grid=false)
@test abs(sim1.𝒪est[:weak_final]-1) < 0.3
end

0 comments on commit 2297905

Please sign in to comment.