From bcd8ce21268a8cb584102d73cff286c5bb6d818c Mon Sep 17 00:00:00 2001 From: Frank Schaefer Date: Mon, 11 Jan 2021 14:23:37 +0100 Subject: [PATCH 01/14] RKMil constant cache perform step and tests --- src/perform_step/low_order.jl | 62 +++++++++++++++++++++- test/commutative_tests.jl | 35 ++++++++++++ test/sde/sde_convergence_tests.jl | 6 +++ test/sde/sde_linear_tests.jl | 5 ++ test/sde/sde_twodimlinear_tests.jl | 4 ++ test/weak_convergence/iip_weak.jl | 29 +++++++--- test/weak_convergence/multidim_iip_weak.jl | 6 +++ test/weak_convergence/oop_weak.jl | 11 ++-- 8 files changed, 145 insertions(+), 13 deletions(-) diff --git a/src/perform_step/low_order.jl b/src/perform_step/low_order.jl index 168805d6d..ce4755dac 100644 --- a/src/perform_step/low_order.jl +++ b/src/perform_step/low_order.jl @@ -246,12 +246,70 @@ end end end +@muladd function perform_step!(integrator,cache::RKMilCommuteConstantCache,f=integrator.f) + @unpack t,dt,uprev,u,W,p = integrator + dW = W.dW; sqdt = integrator.sqdt + + ggprime_norm = 0.0 + + mil_correction = zero(u) + if alg_interpretation(integrator.alg) == :Ito + if typeof(dW) <: Number + WikJ = 0.5 .* ((dW) .* (dW)' .- dt) + else + WikJ = 0.5 .* (vec(dW) .* vec(dW)' .- dt .* Eye{eltype(W.dW)}(length(W.dW))) + end + else + if typeof(dW) <: Number + WikJ = 0.5 * dW .* dW' + else + WikJ = 0.5 .* vec(dW) .* vec(dW)' + end + end + + du1 = integrator.f(uprev,p,t) + L = integrator.g(uprev,p,t) + + K = uprev + dt*du1 + for j = 1:length(dW) + if typeof(dW) <: Number + Kj = K + sqdt*L + else + Kj = K + sqdt*@view(L[:,j]) + end + gtmp = integrator.g(Kj,p,t) + Dgj = (gtmp - L)/sqdt + if integrator.opts.adaptive + ggprime_norm += integrator.opts.internalnorm(Dgj,t) + end + if typeof(dW) <: Number + tmp = Dgj*WikJ + else + tmp = Dgj*@view(WikJ[:,j]) + end + mil_correction += tmp + end + tmp = L*dW + u = uprev + dt*du1 + tmp + mil_correction + + if integrator.opts.adaptive + En = integrator.opts.internalnorm(dW,t)^3*ggprime_norm^2 / 6 + du2 = integrator.f(K,p,t+dt) + tmp = integrator.opts.internalnorm(integrator.opts.delta * dt * (du2 - du1) / 2,t) + En + + tmp = calculate_residuals(tmp, uprev, u, integrator.opts.abstol, + integrator.opts.reltol, integrator.opts.internalnorm, t) + integrator.EEst = integrator.opts.internalnorm(tmp,t) + + end + integrator.u = u +end + @muladd function perform_step!(integrator,cache::RKMilCommuteCache,f=integrator.f) @unpack du1,du2,K,gtmp,L = cache @unpack t,dt,uprev,u,W,p = integrator @unpack WikJ,mil_correction,Kj,Dgj,tmp = cache dW = W.dW; sqdt = integrator.sqdt - f = integrator.f; g = integrator.g ggprime_norm = 0.0 @@ -269,7 +327,7 @@ end for j = 1:length(dW) @.. Kj = K + sqdt*@view(L[:,j]) # This works too #Kj .= uprev .+ sqdt*L[:,j] - g(gtmp,Kj,p,t) + integrator.g(gtmp,Kj,p,t) @.. Dgj = (gtmp - L)/sqdt if integrator.opts.adaptive ggprime_norm += integrator.opts.internalnorm(Dgj,t) diff --git a/test/commutative_tests.jl b/test/commutative_tests.jl index b49bb54b1..8e8560889 100644 --- a/test/commutative_tests.jl +++ b/test/commutative_tests.jl @@ -14,6 +14,11 @@ function f_commute(du,u,p,t) du .+= 1.01u end +function f_commute_oop(u,p,t) + du = A*u + du += 1.01u +end + function f_commute_analytic(u0,p,t,W) tmp = (A+1.01I-2*(B^2))*t + B*sum(W) exp(tmp)*u0 @@ -30,6 +35,13 @@ function σ(du,u,p,t) du[2,4] = σ_const*u[2] end +function σ_oop(u,p,t) + σ_const*[ + u[1] u[1] u[1] u[1] + u[2] u[2] u[2] u[2] + ] +end + ff_commute = SDEFunction(f_commute,σ,analytic=f_commute_analytic) prob = SDEProblem(ff_commute,σ,u0,(0.0,1.0),noise_rate_prototype=rand(2,4)) @@ -51,3 +63,26 @@ sim2 = test_convergence(dts,prob,RKMilGeneral(ii_approx=IICommutative()),traject sim2 = test_convergence(dts,prob,RKMilGeneral(p=2),trajectories=Int(2e2)) @test abs(sim2.𝒪est[:final] - 1.0) < 0.2 + + +ff_commute_oop = SDEFunction(f_commute_oop,σ_oop,analytic=f_commute_analytic) + +proboop = SDEProblem(ff_commute_oop,σ_oop,u0,(0.0,1.0),noise_rate_prototype=rand(2,4)) + +sol = solve(proboop,RKMilCommute(),dt=1/2^(8)) +sol = solve(proboop,RKMilGeneral(ii_approx=IICommutative()),dt=1/2^(8)) +sol = solve(proboop,RKMilGeneral(p=2),dt=1/2^(10)) +sol = solve(proboop,EM(),dt=1/2^(10)) + +dts = (1/2) .^ (10:-1:3) #14->7 good plot +sim2 = test_convergence(dts,proboop,EM(),trajectories=Int(1e2)) +@test abs(sim2.𝒪est[:final] - 0.5) < 0.2 + +sim2 = test_convergence(dts,proboop,RKMilCommute(),trajectories=Int(2e2)) +@test abs(sim2.𝒪est[:final] - 1.0) < 0.2 + +sim2 = test_convergence(dts,proboop,RKMilGeneral(ii_approx=IICommutative()),trajectories=Int(2e2)) +@test abs(sim2.𝒪est[:final] - 1.0) < 0.2 + +sim2 = test_convergence(dts,proboop,RKMilGeneral(p=2),trajectories=Int(2e2)) +@test abs(sim2.𝒪est[:final] - 1.0) < 0.2 diff --git a/test/sde/sde_convergence_tests.jl b/test/sde/sde_convergence_tests.jl index a76046ee2..f39c06520 100644 --- a/test/sde/sde_convergence_tests.jl +++ b/test/sde/sde_convergence_tests.jl @@ -23,6 +23,8 @@ sim = test_convergence(dts,prob,LambaEM(),trajectories=Int(1e2)) @test abs(sim.𝒪est[:l2]-.5) < 0.2 sim2 = test_convergence(dts,prob,RKMil(),trajectories=Int(2e2)) @test abs(sim2.𝒪est[:l∞]-1) < 0.2 +sim2 = test_convergence(dts,prob,RKMilCommute(),trajectories=Int(2e2)) +@test abs(sim2.𝒪est[:l∞]-1) < 0.2 sim2 = test_convergence(dts,prob,RKMilGeneral(),trajectories=Int(2e2)) @test abs(sim2.𝒪est[:l∞]-1) < 0.2 print(".") @@ -101,6 +103,8 @@ sim = test_convergence(dts,prob,ImplicitRKMil(),trajectories=Int(1e2)) @test abs(sim.𝒪est[:l2]-1) < 0.2 sim2 = test_convergence(dts,prob,RKMil(),trajectories=Int(1e2)) @test abs(sim2.𝒪est[:l∞]-1) < 0.22 +sim2 = test_convergence(dts,prob,RKMilCommute(),trajectories=Int(1e2)) +@test abs(sim2.𝒪est[:l∞]-1) < 0.22 sim2 = test_convergence(dts,prob,RKMilGeneral(),trajectories=Int(1e2)) @test abs(sim2.𝒪est[:l∞]-1) < 0.22 print(".") @@ -169,6 +173,8 @@ sim = test_convergence(dts,prob,ImplicitRKMil(),trajectories=Int(1e2)) @test abs(sim.𝒪est[:l2]-1) < 0.2 sim2 = test_convergence(dts,prob,RKMil(),trajectories=Int(1e1)) @test abs(sim2.𝒪est[:l∞]-1) < 0.2 +sim2 = test_convergence(dts,prob,RKMilCommute(),trajectories=Int(1e1)) +@test abs(sim2.𝒪est[:l∞]-1) < 0.2 sim2 = test_convergence(dts,prob,RKMilGeneral(),trajectories=Int(1e1)) @test abs(sim2.𝒪est[:l∞]-1) < 0.2 print(".") diff --git a/test/sde/sde_linear_tests.jl b/test/sde/sde_linear_tests.jl index a92516738..e99dcd394 100644 --- a/test/sde/sde_linear_tests.jl +++ b/test/sde/sde_linear_tests.jl @@ -8,6 +8,7 @@ prob = prob_sde_linear println("Solve and Plot") sol = solve(prob,EM(),dt=1//2^(4)) sol = solve(prob,RKMil(),dt=1//2^(4)) +sol = solve(prob,RKMilCommute(),dt=1//2^(4)) sol = solve(prob,RKMilGeneral(),dt=1//2^(4)) sol = solve(prob,SRI(),dt=1//2^(4)) sol = solve(prob,SRIW1(),dt=1//2^(4)) @@ -22,6 +23,8 @@ sim2 = test_convergence(dts,prob,RKMil(),trajectories=trajectories) sim21 = test_convergence(dts,prob,RKMilGeneral(),trajectories=trajectories) +sim22 = test_convergence(dts,prob,RKMilCommute(),trajectories=trajectories) + sim3 = test_convergence(dts,prob,SRI(),trajectories=trajectories) #TEST_PLOT && plot(plot(sim),plot(sim2),plot(sim3),layout=@layout([a b c]),size=(1200,600)) @@ -30,6 +33,8 @@ sim3 = test_convergence(dts,prob,SRI(),trajectories=trajectories) @test abs(sim.𝒪est[:l2]-.5) + abs(sim21.𝒪est[:l∞]-1) + abs(sim3.𝒪est[:final]-1.5)<.441 #High tolerance since low dts for testing! +@test abs(sim22.𝒪est[:l∞]-1)<.3 + # test reinit integrator = init(prob,EM(),dt=1//2^(4)) solve!(integrator) diff --git a/test/sde/sde_twodimlinear_tests.jl b/test/sde/sde_twodimlinear_tests.jl index bca3ae99b..0121894db 100644 --- a/test/sde/sde_twodimlinear_tests.jl +++ b/test/sde/sde_twodimlinear_tests.jl @@ -10,6 +10,7 @@ sol = solve(prob,EM(),dt=1/2^(3)) sol = solve(prob,LambaEM(),dt=1/2^(3)) sol = solve(prob,LambaEulerHeun(),dt=1/2^(3)) sol = solve(prob,RKMil(),dt=1/2^(3)) +sol = solve(prob,RKMilCommute(),dt=1/2^(3)) sol = solve(prob,RKMilGeneral(),dt=1/2^(3)) sol = solve(prob,SRI(),dt=1/2^(3)) sol = solve(prob,SRIW1(),dt=1/2^(3)) @@ -78,6 +79,9 @@ sim = test_convergence(dts,prob,ImplicitRKMil(nlsolve=StochasticDiffEq.NLFunctio sim2 = test_convergence(dts,prob,RKMil(),trajectories=1000) @test abs(sim2.𝒪est[:l∞]-1) < 0.2 +sim2 = test_convergence(dts,prob,RKMilCommute(),trajectories=1000) +@test abs(sim2.𝒪est[:l∞]-1) < 0.2 + sim2 = test_convergence(dts,prob,RKMilGeneral(),trajectories=1000) @test abs(sim2.𝒪est[:l∞]-1) < 0.2 diff --git a/test/weak_convergence/iip_weak.jl b/test/weak_convergence/iip_weak.jl index 7a8774df9..989f5d7f3 100644 --- a/test/weak_convergence/iip_weak.jl +++ b/test/weak_convergence/iip_weak.jl @@ -1,12 +1,18 @@ using StochasticDiffEq, DiffEqDevTools, Test using Random -using DiffEqProblemLibrary.SDEProblemLibrary: importsdeproblems -importsdeproblems() -using DiffEqProblemLibrary.SDEProblemLibrary: prob_sde_linear + +f_linear_iip(du,u,p,t) = @.(du = 1.01*u) +σ_linear_iip(du,u,p,t) = @.(du = 0.87*u) + +linear_analytic(u0,p,t,W) = @.(u0*exp(0.63155t+0.87W)) + +prob_sde_linear_iip = SDEProblem(SDEFunction(f_linear_iip,σ_linear_iip, + analytic=linear_analytic),σ_linear_iip,[1/2],(0.0,1.0)) + Random.seed!(100) dts = 1 .//2 .^(7:-1:3) #14->7 good plot -prob = prob_sde_linear +prob = prob_sde_linear_iip println("EM") sim = test_convergence(dts,prob,EM(),save_everystep=false,trajectories=Int(1e4), weak_timeseries_errors=false) @@ -15,7 +21,7 @@ sim = test_convergence(dts,prob,EM(),save_everystep=false,trajectories=Int(1e4) #@test abs(sim.𝒪est[:weak_l∞]-1) < 0.3 println("SimplifiedEM") dts = 1 .//2 .^(8:-1:2) -sim = test_convergence(dts,prob,SimplifiedEM(),save_everystep=false,trajectories=Int(5e4), +sim = test_convergence(dts,prob,SimplifiedEM(),save_everystep=false,trajectories=Int(1e5), weak_timeseries_errors=false) @test abs(sim.𝒪est[:weak_final]-1) < 0.3 #@test abs(sim.𝒪est[:weak_l2]-1) < 0.3 @@ -27,6 +33,12 @@ sim = test_convergence(dts,prob,RKMil(),save_everystep=false,trajectories=Int(1e @test abs(sim.𝒪est[:weak_final]-1) < 0.3 #@test abs(sim.𝒪est[:weak_l2]-1) < 0.3 #@test abs(sim.𝒪est[:weak_l∞]-1) < 0.3 +println("RKMilCommute") +sim = test_convergence(dts,prob,RKMilCommute(),save_everystep=false,trajectories=Int(1e4), + weak_timeseries_errors=false) +@test abs(sim.𝒪est[:weak_final]-1) < 0.3 +#@test abs(sim.𝒪est[:weak_l2]-1) < 0.3 +#@test abs(sim.𝒪est[:weak_l∞]-1) < 0.3 println("RKMilGeneral") sim = test_convergence(dts,prob,RKMilGeneral(),save_everystep=false,trajectories=Int(1e4), weak_timeseries_errors=false) @@ -40,9 +52,9 @@ sim = test_convergence(dts,prob,SROCK1(),save_everystep=false,trajectories=Int(4 #@test abs(sim.𝒪est[:weak_l2]-1) < 0.3 #@test abs(sim.𝒪est[:weak_l∞]-1) < 0.3 println("SROCK2") -sim = test_convergence(dts,prob,SROCK2(),save_everystep=false,trajectories=Int(4e4), +sim = test_convergence(dts,prob,SROCK2(),save_everystep=false,trajectories=Int(1e5), weak_timeseries_errors=false) -@test abs(sim.𝒪est[:weak_final]-1.5) < 0.3 +@test abs(sim.𝒪est[:weak_final]-1.5) < 0.31 #@test abs(sim.𝒪est[:weak_l2]-2) < 0.3 #@test abs(sim.𝒪est[:weak_l∞]-2) < 0.3 println("SROCKEM") @@ -64,7 +76,8 @@ sim = test_convergence(dts,prob,SKSROCK(),save_everystep=false,trajectories=Int( #@test abs(sim.𝒪est[:weak_l2]-1) < 0.3 #@test abs(sim.𝒪est[:weak_l∞]-1) < 0.3 println("SROCKC2") -sim = test_convergence(dts,prob,SROCKC2(),save_everystep=false,trajectories=Int(6e4), +#Random.seed!(1) +sim = test_convergence(dts,prob,SROCKC2(),save_everystep=false,trajectories=Int(5e5), weak_timeseries_errors=false) @test abs(sim.𝒪est[:weak_final]-2) < 0.35 #@test abs(sim.𝒪est[:weak_l2]-2) < 0.3 diff --git a/test/weak_convergence/multidim_iip_weak.jl b/test/weak_convergence/multidim_iip_weak.jl index 0a7ced3d0..d57a96eae 100644 --- a/test/weak_convergence/multidim_iip_weak.jl +++ b/test/weak_convergence/multidim_iip_weak.jl @@ -25,6 +25,12 @@ sim = test_convergence(dts,prob,RKMil(),save_everystep=false,trajectories=Int(1e @test abs(sim.𝒪est[:weak_final]-1) < 0.3 #@test abs(sim.𝒪est[:weak_l2]-1) < 0.3 #@test abs(sim.𝒪est[:weak_l∞]-1) < 0.3 +println("RKMilCommute") +sim = test_convergence(dts,prob,RKMilCommute(),save_everystep=false,trajectories=Int(1e4), + weak_timeseries_errors=false) +@test abs(sim.𝒪est[:weak_final]-1) < 0.3 +#@test abs(sim.𝒪est[:weak_l2]-1) < 0.3 +#@test abs(sim.𝒪est[:weak_l∞]-1) < 0.3 println("RKMilGeneral") sim = test_convergence(dts,prob,RKMilGeneral(),save_everystep=false,trajectories=Int(1e4), weak_timeseries_errors=false) diff --git a/test/weak_convergence/oop_weak.jl b/test/weak_convergence/oop_weak.jl index a5abb0131..4bc4b94ef 100644 --- a/test/weak_convergence/oop_weak.jl +++ b/test/weak_convergence/oop_weak.jl @@ -2,7 +2,7 @@ using StochasticDiffEq, DiffEqDevTools, Test using Random using DiffEqProblemLibrary.SDEProblemLibrary: importsdeproblems importsdeproblems() -using DiffEqProblemLibrary.SDEProblemLibrary: prob_sde_linear, prob_sde_2Dlinear +using DiffEqProblemLibrary.SDEProblemLibrary: prob_sde_linear Random.seed!(100) dts = 1 .//2 .^(7:-1:3) #14->7 good plot @@ -24,17 +24,22 @@ sim = test_convergence(dts,prob,RKMil(),save_everystep=false,trajectories=Int(1e @test abs(sim.𝒪est[:weak_final]-1) < 0.3 #@test abs(sim.𝒪est[:weak_l2]-1) < 0.3 #@test abs(sim.𝒪est[:weak_l∞]-1) < 0.3 +println("RKMilCommute") +sim = test_convergence(dts,prob,RKMilCommute(),save_everystep=false,trajectories=Int(1e4)) +@test abs(sim.𝒪est[:weak_final]-1) < 0.3 +#@test abs(sim.𝒪est[:weak_l2]-1) < 0.3 +#@test abs(sim.𝒪est[:weak_l∞]-1) < 0.3 println("RKMilGeneral") sim = test_convergence(dts,prob,RKMilGeneral(),save_everystep=false,trajectories=Int(1e4)) @test abs(sim.𝒪est[:weak_final]-1) < 0.3 #@test abs(sim.𝒪est[:weak_l2]-1) < 0.3 #@test abs(sim.𝒪est[:weak_l∞]-1) < 0.3 -println("RKMilGeneral") +println("SROCK1") sim = test_convergence(dts,prob,SROCK1(),save_everystep=false,trajectories=Int(1e4)) @test abs(sim.𝒪est[:weak_final]-1) < 0.3 #@test abs(sim.𝒪est[:weak_l2]-1) < 0.3 #@test abs(sim.𝒪est[:weak_l∞]-1) < 0.3 -println("RKMilGeneral") +println("SROCK2") dts = 1 .//2 .^(10:-1:1) #14->7 good plot sim = test_convergence(dts,prob,SROCK2(),save_everystep=false,trajectories=Int(4e4)) @test abs(sim.𝒪est[:weak_final]-1.5) < 0.3 From 3021adf3f4e296a49ba195d31cc276aeb87341cd Mon Sep 17 00:00:00 2001 From: Frank Schaefer Date: Mon, 11 Jan 2021 14:24:04 +0100 Subject: [PATCH 02/14] RKMil general small bug fix --- src/iterated_integrals.jl | 2 +- src/solve.jl | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/iterated_integrals.jl b/src/iterated_integrals.jl index 5837ccbca..462f65dd3 100644 --- a/src/iterated_integrals.jl +++ b/src/iterated_integrals.jl @@ -226,7 +226,7 @@ function get_WikJ(ΔW,prob,alg) elseif alg.ii_approx isa IICommutative return WikJCommute_oop() else - return KPWJ_oop #WikJGeneral_oop(ΔW) + return KPWJ_oop() #WikJGeneral_oop(ΔW) end end end diff --git a/src/solve.jl b/src/solve.jl index 7b731ddd9..77089137e 100644 --- a/src/solve.jl +++ b/src/solve.jl @@ -362,7 +362,7 @@ function DiffEqBase.__init( W = WienerProcess(t,rand_prototype,rand_prototype2, save_everystep=save_noise, rng = Xorshifts.Xoroshiro128Plus(_seed)) - elseif alg===RKMilGeneral() + elseif typeof(alg) <: RKMilGeneral m = length(rand_prototype) if typeof(rand_prototype) <: Number || alg.p == nothing rand_prototype2 = nothing From 1044aabbf3f2ad3bfd91e692a915ed945f927d61 Mon Sep 17 00:00:00 2001 From: Frank Schaefer Date: Mon, 11 Jan 2021 18:31:00 +0100 Subject: [PATCH 03/14] handle iterated integral and diagonal case --- src/caches/basic_method_caches.jl | 14 +++-- src/iterated_integrals.jl | 16 +++++ src/perform_step/low_order.jl | 99 +++++++++++++++++++------------ 3 files changed, 84 insertions(+), 45 deletions(-) diff --git a/src/caches/basic_method_caches.jl b/src/caches/basic_method_caches.jl index dcca27c40..0bc22bc28 100644 --- a/src/caches/basic_method_caches.jl +++ b/src/caches/basic_method_caches.jl @@ -115,7 +115,9 @@ function alg_cache(alg::RKMil,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototy RKMilCache(u,uprev,du1,du2,K,tmp,L) end -struct RKMilCommuteConstantCache <: StochasticDiffEqConstantCache end +struct RKMilCommuteConstantCache{WikType} <: StochasticDiffEqConstantCache + WikJ::WikType +end @cache struct RKMilCommuteCache{uType,rateType,rateNoiseType,WikType} <: StochasticDiffEqMutableCache u::uType uprev::uType @@ -125,24 +127,24 @@ struct RKMilCommuteConstantCache <: StochasticDiffEqConstantCache end gtmp::rateNoiseType L::rateNoiseType WikJ::WikType - Dg::WikType mil_correction::rateType Kj::uType Dgj::rateNoiseType tmp::uType end -alg_cache(alg::RKMilCommute,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,jump_rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{false}}) = RKMilCommuteConstantCache() +function alg_cache(alg::RKMilCommute,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,jump_rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{false}}) WikJ = WikJ = get_WikJ(ΔW,prob,alg) + RKMilCommuteConstantCache{typeof(WikJ)}(WikJ) +end function alg_cache(alg::RKMilCommute,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,jump_rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{true}}) du1 = zero(rate_prototype); du2 = zero(rate_prototype) K = zero(rate_prototype); gtmp = zero(noise_rate_prototype); L = zero(noise_rate_prototype); tmp = zero(rate_prototype) - WikJ = false .* vec(ΔW) .* vec(ΔW)' - Dg = false .* vec(ΔW) .* vec(ΔW)' + WikJ = get_WikJ(ΔW,prob,alg) mil_correction = zero(rate_prototype) Kj = zero(u); Dgj = zero(noise_rate_prototype) - RKMilCommuteCache(u,uprev,du1,du2,K,gtmp,L,WikJ,Dg,mil_correction,Kj,Dgj,tmp) + RKMilCommuteCache(u,uprev,du1,du2,K,gtmp,L,WikJ,mil_correction,Kj,Dgj,tmp) end struct RKMilGeneralConstantCache{WikType} <: StochasticDiffEqConstantCache diff --git a/src/iterated_integrals.jl b/src/iterated_integrals.jl index 462f65dd3..3fa2bb956 100644 --- a/src/iterated_integrals.jl +++ b/src/iterated_integrals.jl @@ -231,6 +231,22 @@ function get_WikJ(ΔW,prob,alg) end end +function get_WikJ(ΔW,prob,alg::RKMilCommute) + if isinplace(prob) + if typeof(ΔW) <: Number || is_diagonal_noise(prob) + return WikJDiagonal_iip(ΔW) + else + return WikJCommute_iip(ΔW) + end + else + if typeof(ΔW) <: Number || is_diagonal_noise(prob) + return WikJDiagonal_oop() + else + return WikJCommute_oop() + end + end +end + function get_iterated_I(dt, dW, dZ, Wik::WikJDiagonal_oop, p=nothing, C=1, γ=1//1) WikJ = 1//2 .* dW .* dW WikJ diff --git a/src/perform_step/low_order.jl b/src/perform_step/low_order.jl index ce4755dac..7ad86136d 100644 --- a/src/perform_step/low_order.jl +++ b/src/perform_step/low_order.jl @@ -249,21 +249,18 @@ end @muladd function perform_step!(integrator,cache::RKMilCommuteConstantCache,f=integrator.f) @unpack t,dt,uprev,u,W,p = integrator dW = W.dW; sqdt = integrator.sqdt + Wik = cache.WikJ ggprime_norm = 0.0 + WikJ = get_iterated_I(dt, dW, W.dZ, Wik) + mil_correction = zero(u) if alg_interpretation(integrator.alg) == :Ito - if typeof(dW) <: Number - WikJ = 0.5 .* ((dW) .* (dW)' .- dt) - else - WikJ = 0.5 .* (vec(dW) .* vec(dW)' .- dt .* Eye{eltype(W.dW)}(length(W.dW))) - end - else - if typeof(dW) <: Number - WikJ = 0.5 * dW .* dW' + if typeof(dW) <: Number || is_diagonal_noise(integrator.sol.prob) + WikJ = WikJ .- 1//2 .* dt else - WikJ = 0.5 .* vec(dW) .* vec(dW)' + WikJ -= 1//2 .* UniformScaling(dt) end end @@ -271,27 +268,35 @@ end L = integrator.g(uprev,p,t) K = uprev + dt*du1 - for j = 1:length(dW) - if typeof(dW) <: Number - Kj = K + sqdt*L - else - Kj = K + sqdt*@view(L[:,j]) - end - gtmp = integrator.g(Kj,p,t) + + if is_diagonal_noise(integrator.sol.prob) + tmp = (alg_interpretation(integrator.alg) == :Ito ? K : uprev) .+ integrator.sqdt .* L + gtmp = integrator.g(tmp,p,t) Dgj = (gtmp - L)/sqdt - if integrator.opts.adaptive + ggprime_norm = integrator.opts.internalnorm(Dgj,t) + u = @.. K + L*dW + Dgj*WikJ + else + for j = 1:length(dW) + if typeof(dW) <: Number + Kj = K + sqdt*L + else + Kj = K + sqdt*@view(L[:,j]) + end + gtmp = integrator.g(Kj,p,t) + Dgj = (gtmp - L)/sqdt + if integrator.opts.adaptive ggprime_norm += integrator.opts.internalnorm(Dgj,t) + end + if typeof(dW) <: Number + tmp = Dgj*WikJ + else + tmp = Dgj*@view(WikJ[:,j]) + end + mil_correction += tmp end - if typeof(dW) <: Number - tmp = Dgj*WikJ - else - tmp = Dgj*@view(WikJ[:,j]) - end - mil_correction += tmp + tmp = L*dW + u = uprev + dt*du1 + tmp + mil_correction end - tmp = L*dW - u = uprev + dt*du1 + tmp + mil_correction - if integrator.opts.adaptive En = integrator.opts.internalnorm(dW,t)^3*ggprime_norm^2 / 6 du2 = integrator.f(K,p,t+dt) @@ -313,30 +318,46 @@ end ggprime_norm = 0.0 + Wik = cache.WikJ + + get_iterated_I!(dt, dW, W.dZ, Wik) + WikJ = Wik.WikJ + @.. mil_correction = zero(u) if alg_interpretation(integrator.alg) == :Ito - WikJ .= 0.5 .* (vec(dW) .* vec(dW)' .- dt .* Eye{eltype(W.dW)}(length(W.dW))) - else - WikJ .= 0.5 .* vec(dW) .* vec(dW)' + if typeof(dW) <: Number || is_diagonal_noise(integrator.sol.prob) + @.. WikJ -= 1//2*dt + else + WikJ -= 1//2 .* UniformScaling(dt) + end end integrator.f(du1,uprev,p,t) integrator.g(L,uprev,p,t) @.. K = uprev + dt*du1 - for j = 1:length(dW) - @.. Kj = K + sqdt*@view(L[:,j]) # This works too - #Kj .= uprev .+ sqdt*L[:,j] - integrator.g(gtmp,Kj,p,t) + + if is_diagonal_noise(integrator.sol.prob) + tmp .= (alg_interpretation(integrator.alg) == :Ito ? K : uprev) .+ integrator.sqdt .* L + integrator.g(gtmp,tmp,p,t) @.. Dgj = (gtmp - L)/sqdt - if integrator.opts.adaptive - ggprime_norm += integrator.opts.internalnorm(Dgj,t) + ggprime_norm = integrator.opts.internalnorm(Dgj,t) + @.. u = K + L*dW + Dgj*WikJ + else + for j = 1:length(dW) + @.. Kj = K + sqdt*@view(L[:,j]) # This works too + #Kj .= uprev .+ sqdt*L[:,j] + integrator.g(gtmp,Kj,p,t) + @.. Dgj = (gtmp - L)/sqdt + if integrator.opts.adaptive + ggprime_norm += integrator.opts.internalnorm(Dgj,t) + end + mul!(tmp,Dgj,@view(WikJ[:,j])) + mil_correction .+= tmp end - mul!(tmp,Dgj,@view(WikJ[:,j])) - mil_correction .+= tmp + mul!(tmp,L,dW) + @.. u .= uprev + dt*du1 + tmp + mil_correction end - mul!(tmp,L,dW) - @.. u .= uprev + dt*du1 + tmp + mil_correction if integrator.opts.adaptive En = integrator.opts.internalnorm(W.dW,t)^3*ggprime_norm^2 / 6 From 4fd706a1279c1a0d7a2d3577e289ee9e5d0196b4 Mon Sep 17 00:00:00 2001 From: Frank Schaefer Date: Mon, 11 Jan 2021 21:56:51 +0100 Subject: [PATCH 04/14] add local (test error looks random) --- test/sde/sde_twodimlinear_tests.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/test/sde/sde_twodimlinear_tests.jl b/test/sde/sde_twodimlinear_tests.jl index 0121894db..9bafc8aa1 100644 --- a/test/sde/sde_twodimlinear_tests.jl +++ b/test/sde/sde_twodimlinear_tests.jl @@ -109,6 +109,7 @@ print(".") eigen_est = (integrator) -> integrator.eigen_est = 10.0 for Alg in [SROCK1, SROCK2], alg in [Alg(), Alg(eigen_est=eigen_est)] + local sim2 sim2 = test_convergence(dts,prob,alg,trajectories=100) @test abs(sim2.𝒪est[:l∞]-1) < 0.2 end From 613e90fd5ff796d22ca8325b5625165fafc577a2 Mon Sep 17 00:00:00 2001 From: Frank Schaefer Date: Mon, 11 Jan 2021 22:50:26 +0100 Subject: [PATCH 05/14] try more trajectories --- test/sde/sde_twodimlinear_tests.jl | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/test/sde/sde_twodimlinear_tests.jl b/test/sde/sde_twodimlinear_tests.jl index 9bafc8aa1..1edf109da 100644 --- a/test/sde/sde_twodimlinear_tests.jl +++ b/test/sde/sde_twodimlinear_tests.jl @@ -107,11 +107,13 @@ sim2 = test_convergence(dts,prob,WangLi3SMil_F(),trajectories=100) print(".") +Random.seed!(100) eigen_est = (integrator) -> integrator.eigen_est = 10.0 for Alg in [SROCK1, SROCK2], alg in [Alg(), Alg(eigen_est=eigen_est)] local sim2 - sim2 = test_convergence(dts,prob,alg,trajectories=100) + sim2 = test_convergence(dts,prob,alg,trajectories=150) @test abs(sim2.𝒪est[:l∞]-1) < 0.2 + @show sim2.𝒪est[:l∞] end sim2 = test_convergence(dts,prob,SROCKEM(strong_order_1=false),trajectories=100) From b2090ec3d766825188612b191dfd619c4b97b8df Mon Sep 17 00:00:00 2001 From: Frank Schaefer Date: Tue, 12 Jan 2021 20:06:50 +0100 Subject: [PATCH 06/14] fix week iip tests --- test/gpu/sde_weak_brusselator_adaptive.jl | 2 +- test/weak_convergence/iip_weak.jl | 5 +++-- 2 files changed, 4 insertions(+), 3 deletions(-) diff --git a/test/gpu/sde_weak_brusselator_adaptive.jl b/test/gpu/sde_weak_brusselator_adaptive.jl index d195a31b9..34dfa93e0 100644 --- a/test/gpu/sde_weak_brusselator_adaptive.jl +++ b/test/gpu/sde_weak_brusselator_adaptive.jl @@ -45,7 +45,7 @@ ensembleprob = EnsembleProblem(prob, prob_func = prob_func) #Performance check with nvvp # CUDAnative.CUDAdrv.@profile # check either on CPU with EnsembleCPUArray() or on GPU with EnsembleGPUArray() -sol = @time solve(ensembleprob,DRI1(),EnsembleCPUArray(),trajectories=numtraj) +@test_nowarn sol = @time solve(ensembleprob,DRI1(),EnsembleCPUArray(),trajectories=numtraj) #sol = @time solve(ensembleprob,DRI1(),EnsembleGPUArray(),trajectories=numtraj) diff --git a/test/weak_convergence/iip_weak.jl b/test/weak_convergence/iip_weak.jl index 989f5d7f3..32fa22388 100644 --- a/test/weak_convergence/iip_weak.jl +++ b/test/weak_convergence/iip_weak.jl @@ -52,9 +52,10 @@ sim = test_convergence(dts,prob,SROCK1(),save_everystep=false,trajectories=Int(4 #@test abs(sim.𝒪est[:weak_l2]-1) < 0.3 #@test abs(sim.𝒪est[:weak_l∞]-1) < 0.3 println("SROCK2") +Random.seed!(2) sim = test_convergence(dts,prob,SROCK2(),save_everystep=false,trajectories=Int(1e5), weak_timeseries_errors=false) -@test abs(sim.𝒪est[:weak_final]-1.5) < 0.31 +@test abs(sim.𝒪est[:weak_final]-1.5) < 0.3 #@test abs(sim.𝒪est[:weak_l2]-2) < 0.3 #@test abs(sim.𝒪est[:weak_l∞]-2) < 0.3 println("SROCKEM") @@ -76,7 +77,7 @@ sim = test_convergence(dts,prob,SKSROCK(),save_everystep=false,trajectories=Int( #@test abs(sim.𝒪est[:weak_l2]-1) < 0.3 #@test abs(sim.𝒪est[:weak_l∞]-1) < 0.3 println("SROCKC2") -#Random.seed!(1) +Random.seed!(3) sim = test_convergence(dts,prob,SROCKC2(),save_everystep=false,trajectories=Int(5e5), weak_timeseries_errors=false) @test abs(sim.𝒪est[:weak_final]-2) < 0.35 From fc110003b069100ce19a25cc1ad55e7d82dc88d0 Mon Sep 17 00:00:00 2001 From: Frank Schaefer Date: Tue, 12 Jan 2021 23:08:14 +0100 Subject: [PATCH 07/14] double weak adaptive time out --- .buildkite/pipeline.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.buildkite/pipeline.yml b/.buildkite/pipeline.yml index c7a460aaf..acd6ebb4b 100644 --- a/.buildkite/pipeline.yml +++ b/.buildkite/pipeline.yml @@ -29,7 +29,7 @@ steps: fastcpu: "true" env: GROUP: 'WeakAdaptive' - timeout_in_minutes: 120 + timeout_in_minutes: 240 # Don't run Buildkite if the commit message includes the text [skip tests] if: build.message !~ /\[skip tests\]/ From f25420f231bad91abfdf8979df3cefe45955ad9e Mon Sep 17 00:00:00 2001 From: Frank Schaefer Date: Wed, 13 Jan 2021 11:01:51 +0100 Subject: [PATCH 08/14] move WeakConvergence1 to buildkite pipeline with larger time window --- .buildkite/pipeline.yml | 2 +- .github/workflows/CI.yml | 1 - test/runtests.jl | 19 ++++++++++--------- 3 files changed, 11 insertions(+), 11 deletions(-) diff --git a/.buildkite/pipeline.yml b/.buildkite/pipeline.yml index acd6ebb4b..889210e06 100644 --- a/.buildkite/pipeline.yml +++ b/.buildkite/pipeline.yml @@ -12,7 +12,7 @@ steps: fastcpu: "true" env: GROUP: 'WeakConvergence' - timeout_in_minutes: 120 + timeout_in_minutes: 240 # Don't run Buildkite if the commit message includes the text [skip tests] if: build.message !~ /\[skip tests\]/ diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index afa7f512a..39c0e36e4 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -22,7 +22,6 @@ jobs: - WeakConvergence2 - WeakConvergence3 - WeakConvergence4 - - WeakConvergence5 steps: - uses: actions/checkout@v2 - uses: julia-actions/setup-julia@v1 diff --git a/test/runtests.jl b/test/runtests.jl index 8d9683991..6c10fd63c 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -74,30 +74,31 @@ const is_APPVEYOR = Sys.iswindows() && haskey(ENV,"APPVEYOR") end if !is_APPVEYOR && (GROUP == "All" || GROUP == "WeakConvergence1") - @time @safetestset "OOP Weak Convergence Tests" begin include("weak_convergence/oop_weak.jl") end - @time @safetestset "Additive Weak Convergence Tests" begin include("weak_convergence/additive_weak.jl") end - @time @safetestset "IIP Weak Convergence Tests" begin include("weak_convergence/iip_weak.jl") end - end - - if !is_APPVEYOR && (GROUP == "All" || GROUP == "WeakConvergence2") @time @safetestset "Multidimensional IIP Weak Convergence Tests" begin include("weak_convergence/multidim_iip_weak.jl") end @time @safetestset "Platen's PL1WM weak second order" begin include("weak_convergence/PL1WM.jl") end end - if !is_APPVEYOR && (GROUP == "All" || GROUP == "WeakConvergence3") + if !is_APPVEYOR && (GROUP == "All" || GROUP == "WeakConvergence2") @time @safetestset "Roessler weak SRK Tests" begin include("weak_convergence/srk_weak_final.jl") end end - if !is_APPVEYOR && (GROUP == "All" || GROUP == "WeakConvergence4") + if !is_APPVEYOR && (GROUP == "All" || GROUP == "WeakConvergence3") @time @safetestset "Roessler weak SRK (non-diagonal) Tests" begin include("weak_convergence/srk_weak_final_non_diagonal.jl") end @time @safetestset "Weak Stratonovich Tests" begin include("weak_convergence/weak_strat.jl") end end - if !is_APPVEYOR && (GROUP == "All" || GROUP == "WeakConvergence5") + if !is_APPVEYOR && (GROUP == "All" || GROUP == "WeakConvergence4") @time @safetestset "Weak Stratonovich (non-diagonal) Tests" begin include("weak_convergence/weak_strat_non_diagonal.jl") end @time @safetestset "SIE SME weak Tests" begin include("weak_convergence/SIE_SME.jl") end end + if !is_APPVEYOR && GROUP == "WeakConvergence" + #activate_gpu_env() + @time @safetestset "OOP Weak Convergence Tests" begin include("weak_convergence/oop_weak.jl") end + @time @safetestset "Additive Weak Convergence Tests" begin include("weak_convergence/additive_weak.jl") end + @time @safetestset "IIP Weak Convergence Tests" begin include("weak_convergence/iip_weak.jl") end + end + if !is_APPVEYOR && GROUP == "WeakAdaptive" activate_gpu_env() @time @safetestset "Weak adaptive step size Brusselator " begin include("gpu/sde_weak_brusselator_adaptive.jl") end From 8b1cf1693f9f8731b27e2b95aa4e0feb72a76215 Mon Sep 17 00:00:00 2001 From: Frank Schaefer Date: Wed, 13 Jan 2021 14:14:06 +0100 Subject: [PATCH 09/14] use in oop_weak same seed as in iip_weak --- test/weak_convergence/oop_weak.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/weak_convergence/oop_weak.jl b/test/weak_convergence/oop_weak.jl index 4bc4b94ef..2fdd6ee7c 100644 --- a/test/weak_convergence/oop_weak.jl +++ b/test/weak_convergence/oop_weak.jl @@ -62,8 +62,8 @@ sim = test_convergence(dts,prob,SKSROCK(),save_everystep=false,trajectories=Int( #@test abs(sim.𝒪est[:weak_l2]-1) < 0.3 #@test abs(sim.𝒪est[:weak_l∞]-1) < 0.3 println("SROCKC2") -dts = 1 .//2 .^(7:-1:2) -sim = test_convergence(dts,prob,SROCKC2(),save_everystep=false,trajectories=Int(1e4)) +Random.seed!(3) +sim = test_convergence(dts,prob,SROCKC2(),save_everystep=false,trajectories=Int(5e5)) @test abs(sim.𝒪est[:weak_final]-2) < 0.3 #@test abs(sim.𝒪est[:weak_l2]-2) < 0.3 #@test abs(sim.𝒪est[:weak_l∞]-2) < 0.3 From 15db8715dd139ba5c55cce503497ee23a931aba3 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Wed, 13 Jan 2021 11:06:52 -0500 Subject: [PATCH 10/14] Update pipeline.yml --- .buildkite/pipeline.yml | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/.buildkite/pipeline.yml b/.buildkite/pipeline.yml index 889210e06..1d5daf8c3 100644 --- a/.buildkite/pipeline.yml +++ b/.buildkite/pipeline.yml @@ -7,9 +7,7 @@ steps: # - JuliaCI/julia-coverage#v0.3: # codecov: true agents: - queue: "juliagpu" - cuda: "*" - fastcpu: "true" + queue: "julia" env: GROUP: 'WeakConvergence' timeout_in_minutes: 240 From 710f624b591387613bf57869bb5aea66417a19ef Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Wed, 13 Jan 2021 11:23:08 -0500 Subject: [PATCH 11/14] Update pipeline.yml --- .buildkite/pipeline.yml | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/.buildkite/pipeline.yml b/.buildkite/pipeline.yml index 1d5daf8c3..b009e443f 100644 --- a/.buildkite/pipeline.yml +++ b/.buildkite/pipeline.yml @@ -7,7 +7,8 @@ steps: # - JuliaCI/julia-coverage#v0.3: # codecov: true agents: - queue: "julia" + queue: "juliagpu" + fastcpu: "true" env: GROUP: 'WeakConvergence' timeout_in_minutes: 240 From 0e935059e7bfeaaa3e4dc7fce50a4bb481db6d3e Mon Sep 17 00:00:00 2001 From: Frank Schaefer Date: Wed, 13 Jan 2021 23:44:21 +0100 Subject: [PATCH 12/14] next try for SROCK methods --- test/weak_convergence/iip_weak.jl | 12 +++++------- test/weak_convergence/oop_weak.jl | 10 ++++------ 2 files changed, 9 insertions(+), 13 deletions(-) diff --git a/test/weak_convergence/iip_weak.jl b/test/weak_convergence/iip_weak.jl index 32fa22388..8e9a77f4c 100644 --- a/test/weak_convergence/iip_weak.jl +++ b/test/weak_convergence/iip_weak.jl @@ -20,14 +20,12 @@ sim = test_convergence(dts,prob,EM(),save_everystep=false,trajectories=Int(1e4) #@test abs(sim.𝒪est[:weak_l2]-1) < 0.3 #@test abs(sim.𝒪est[:weak_l∞]-1) < 0.3 println("SimplifiedEM") -dts = 1 .//2 .^(8:-1:2) -sim = test_convergence(dts,prob,SimplifiedEM(),save_everystep=false,trajectories=Int(1e5), +sim = test_convergence(dts,prob,SimplifiedEM(),save_everystep=false,trajectories=Int(5e4), weak_timeseries_errors=false) @test abs(sim.𝒪est[:weak_final]-1) < 0.3 #@test abs(sim.𝒪est[:weak_l2]-1) < 0.3 #@test abs(sim.𝒪est[:weak_l∞]-1) < 0.35 println("RKMil") -dts = 1 .//2 .^(7:-1:3) sim = test_convergence(dts,prob,RKMil(),save_everystep=false,trajectories=Int(1e4), weak_timeseries_errors=false) @test abs(sim.𝒪est[:weak_final]-1) < 0.3 @@ -52,12 +50,13 @@ sim = test_convergence(dts,prob,SROCK1(),save_everystep=false,trajectories=Int(4 #@test abs(sim.𝒪est[:weak_l2]-1) < 0.3 #@test abs(sim.𝒪est[:weak_l∞]-1) < 0.3 println("SROCK2") -Random.seed!(2) -sim = test_convergence(dts,prob,SROCK2(),save_everystep=false,trajectories=Int(1e5), +dts = 1 .//2 .^(8:-1:1) #14->7 good plot +sim = test_convergence(dts,prob,SROCK2(),save_everystep=false,trajectories=Int(1e4), weak_timeseries_errors=false) @test abs(sim.𝒪est[:weak_final]-1.5) < 0.3 #@test abs(sim.𝒪est[:weak_l2]-2) < 0.3 #@test abs(sim.𝒪est[:weak_l∞]-2) < 0.3 +dts = 1 .//2 .^(7:-1:2) #14->7 good plot println("SROCKEM") sim = test_convergence(dts,prob,SROCKEM(strong_order_1=false),save_everystep=false,trajectories=Int(1e4), weak_timeseries_errors=false) @@ -77,8 +76,7 @@ sim = test_convergence(dts,prob,SKSROCK(),save_everystep=false,trajectories=Int( #@test abs(sim.𝒪est[:weak_l2]-1) < 0.3 #@test abs(sim.𝒪est[:weak_l∞]-1) < 0.3 println("SROCKC2") -Random.seed!(3) -sim = test_convergence(dts,prob,SROCKC2(),save_everystep=false,trajectories=Int(5e5), +sim = test_convergence(dts,prob,SROCKC2(),save_everystep=false,trajectories=Int(1e6), weak_timeseries_errors=false) @test abs(sim.𝒪est[:weak_final]-2) < 0.35 #@test abs(sim.𝒪est[:weak_l2]-2) < 0.3 diff --git a/test/weak_convergence/oop_weak.jl b/test/weak_convergence/oop_weak.jl index 2fdd6ee7c..485d35227 100644 --- a/test/weak_convergence/oop_weak.jl +++ b/test/weak_convergence/oop_weak.jl @@ -13,13 +13,11 @@ println("EM") #@test abs(sim.𝒪est[:weak_l2]-1) < 0.3 #@test abs(sim.𝒪est[:weak_l∞]-1) < 0.3 println("SimplifiedEM") -dts = 1 .//2 .^(8:-1:2) @time sim = test_convergence(dts,prob,SimplifiedEM(),save_everystep=false,trajectories=Int(5e4)) @test abs(sim.𝒪est[:weak_final]-1) < 0.3 #@test abs(sim.𝒪est[:weak_l2]-1) < 0.3 #@test abs(sim.𝒪est[:weak_l∞]-1) < 0.35 println("RKMil") -dts = 1 .//2 .^(7:-1:3) sim = test_convergence(dts,prob,RKMil(),save_everystep=false,trajectories=Int(1e4)) @test abs(sim.𝒪est[:weak_final]-1) < 0.3 #@test abs(sim.𝒪est[:weak_l2]-1) < 0.3 @@ -40,8 +38,8 @@ sim = test_convergence(dts,prob,SROCK1(),save_everystep=false,trajectories=Int(1 #@test abs(sim.𝒪est[:weak_l2]-1) < 0.3 #@test abs(sim.𝒪est[:weak_l∞]-1) < 0.3 println("SROCK2") -dts = 1 .//2 .^(10:-1:1) #14->7 good plot -sim = test_convergence(dts,prob,SROCK2(),save_everystep=false,trajectories=Int(4e4)) +dts = 1 .//2 .^(8:-1:1) #14->7 good plot +sim = test_convergence(dts,prob,SROCK2(),save_everystep=false,trajectories=Int(1e4)) @test abs(sim.𝒪est[:weak_final]-1.5) < 0.3 #@test abs(sim.𝒪est[:weak_l2]-2) < 0.3 #@test abs(sim.𝒪est[:weak_l∞]-2) < 0.3 @@ -62,8 +60,8 @@ sim = test_convergence(dts,prob,SKSROCK(),save_everystep=false,trajectories=Int( #@test abs(sim.𝒪est[:weak_l2]-1) < 0.3 #@test abs(sim.𝒪est[:weak_l∞]-1) < 0.3 println("SROCKC2") -Random.seed!(3) -sim = test_convergence(dts,prob,SROCKC2(),save_everystep=false,trajectories=Int(5e5)) +sim = test_convergence(dts,prob,SROCKC2(),save_everystep=false,trajectories=Int(1e6)) +@show sim.𝒪est[:weak_final] @test abs(sim.𝒪est[:weak_final]-2) < 0.3 #@test abs(sim.𝒪est[:weak_l2]-2) < 0.3 #@test abs(sim.𝒪est[:weak_l∞]-2) < 0.3 From be59cbceecd086dbf5586a1dde382522d74816d7 Mon Sep 17 00:00:00 2001 From: Frank Schaefer Date: Thu, 14 Jan 2021 11:01:18 +0100 Subject: [PATCH 13/14] update pipeline --- .buildkite/pipeline.yml | 2 +- test/weak_convergence/iip_weak.jl | 3 ++- test/weak_convergence/oop_weak.jl | 2 +- 3 files changed, 4 insertions(+), 3 deletions(-) diff --git a/.buildkite/pipeline.yml b/.buildkite/pipeline.yml index b009e443f..7202337c8 100644 --- a/.buildkite/pipeline.yml +++ b/.buildkite/pipeline.yml @@ -11,7 +11,7 @@ steps: fastcpu: "true" env: GROUP: 'WeakConvergence' - timeout_in_minutes: 240 + timeout_in_minutes: 480 # Don't run Buildkite if the commit message includes the text [skip tests] if: build.message !~ /\[skip tests\]/ diff --git a/test/weak_convergence/iip_weak.jl b/test/weak_convergence/iip_weak.jl index 8e9a77f4c..ae847b922 100644 --- a/test/weak_convergence/iip_weak.jl +++ b/test/weak_convergence/iip_weak.jl @@ -76,8 +76,9 @@ sim = test_convergence(dts,prob,SKSROCK(),save_everystep=false,trajectories=Int( #@test abs(sim.𝒪est[:weak_l2]-1) < 0.3 #@test abs(sim.𝒪est[:weak_l∞]-1) < 0.3 println("SROCKC2") -sim = test_convergence(dts,prob,SROCKC2(),save_everystep=false,trajectories=Int(1e6), +@time sim = test_convergence(dts,prob,SROCKC2(),save_everystep=false,trajectories=Int(1e6), weak_timeseries_errors=false) +@show sim.𝒪est[:weak_final] @test abs(sim.𝒪est[:weak_final]-2) < 0.35 #@test abs(sim.𝒪est[:weak_l2]-2) < 0.3 #@test abs(sim.𝒪est[:weak_l∞]-2) < 0.3 diff --git a/test/weak_convergence/oop_weak.jl b/test/weak_convergence/oop_weak.jl index 485d35227..028da7df6 100644 --- a/test/weak_convergence/oop_weak.jl +++ b/test/weak_convergence/oop_weak.jl @@ -60,7 +60,7 @@ sim = test_convergence(dts,prob,SKSROCK(),save_everystep=false,trajectories=Int( #@test abs(sim.𝒪est[:weak_l2]-1) < 0.3 #@test abs(sim.𝒪est[:weak_l∞]-1) < 0.3 println("SROCKC2") -sim = test_convergence(dts,prob,SROCKC2(),save_everystep=false,trajectories=Int(1e6)) +@time sim = test_convergence(dts,prob,SROCKC2(),save_everystep=false,trajectories=Int(1e6)) @show sim.𝒪est[:weak_final] @test abs(sim.𝒪est[:weak_final]-2) < 0.3 #@test abs(sim.𝒪est[:weak_l2]-2) < 0.3 From 82a42a3dfa4d8e8447315026d7e58872863050da Mon Sep 17 00:00:00 2001 From: Frank Schaefer Date: Thu, 14 Jan 2021 16:45:06 +0100 Subject: [PATCH 14/14] increase trajectories --- test/weak_convergence/iip_weak.jl | 4 ++-- test/weak_convergence/oop_weak.jl | 6 +++--- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/test/weak_convergence/iip_weak.jl b/test/weak_convergence/iip_weak.jl index ae847b922..0be0d96eb 100644 --- a/test/weak_convergence/iip_weak.jl +++ b/test/weak_convergence/iip_weak.jl @@ -20,7 +20,7 @@ sim = test_convergence(dts,prob,EM(),save_everystep=false,trajectories=Int(1e4) #@test abs(sim.𝒪est[:weak_l2]-1) < 0.3 #@test abs(sim.𝒪est[:weak_l∞]-1) < 0.3 println("SimplifiedEM") -sim = test_convergence(dts,prob,SimplifiedEM(),save_everystep=false,trajectories=Int(5e4), +sim = test_convergence(dts,prob,SimplifiedEM(),save_everystep=false,trajectories=Int(1e5), weak_timeseries_errors=false) @test abs(sim.𝒪est[:weak_final]-1) < 0.3 #@test abs(sim.𝒪est[:weak_l2]-1) < 0.3 @@ -51,7 +51,7 @@ sim = test_convergence(dts,prob,SROCK1(),save_everystep=false,trajectories=Int(4 #@test abs(sim.𝒪est[:weak_l∞]-1) < 0.3 println("SROCK2") dts = 1 .//2 .^(8:-1:1) #14->7 good plot -sim = test_convergence(dts,prob,SROCK2(),save_everystep=false,trajectories=Int(1e4), +sim = test_convergence(dts,prob,SROCK2(),save_everystep=false,trajectories=Int(5e4), weak_timeseries_errors=false) @test abs(sim.𝒪est[:weak_final]-1.5) < 0.3 #@test abs(sim.𝒪est[:weak_l2]-2) < 0.3 diff --git a/test/weak_convergence/oop_weak.jl b/test/weak_convergence/oop_weak.jl index 028da7df6..af23f27ee 100644 --- a/test/weak_convergence/oop_weak.jl +++ b/test/weak_convergence/oop_weak.jl @@ -13,7 +13,7 @@ println("EM") #@test abs(sim.𝒪est[:weak_l2]-1) < 0.3 #@test abs(sim.𝒪est[:weak_l∞]-1) < 0.3 println("SimplifiedEM") -@time sim = test_convergence(dts,prob,SimplifiedEM(),save_everystep=false,trajectories=Int(5e4)) +@time sim = test_convergence(dts,prob,SimplifiedEM(),save_everystep=false,trajectories=Int(1e5)) @test abs(sim.𝒪est[:weak_final]-1) < 0.3 #@test abs(sim.𝒪est[:weak_l2]-1) < 0.3 #@test abs(sim.𝒪est[:weak_l∞]-1) < 0.35 @@ -39,7 +39,7 @@ sim = test_convergence(dts,prob,SROCK1(),save_everystep=false,trajectories=Int(1 #@test abs(sim.𝒪est[:weak_l∞]-1) < 0.3 println("SROCK2") dts = 1 .//2 .^(8:-1:1) #14->7 good plot -sim = test_convergence(dts,prob,SROCK2(),save_everystep=false,trajectories=Int(1e4)) +sim = test_convergence(dts,prob,SROCK2(),save_everystep=false,trajectories=Int(5e4)) @test abs(sim.𝒪est[:weak_final]-1.5) < 0.3 #@test abs(sim.𝒪est[:weak_l2]-2) < 0.3 #@test abs(sim.𝒪est[:weak_l∞]-2) < 0.3 @@ -62,7 +62,7 @@ sim = test_convergence(dts,prob,SKSROCK(),save_everystep=false,trajectories=Int( println("SROCKC2") @time sim = test_convergence(dts,prob,SROCKC2(),save_everystep=false,trajectories=Int(1e6)) @show sim.𝒪est[:weak_final] -@test abs(sim.𝒪est[:weak_final]-2) < 0.3 +@test abs(sim.𝒪est[:weak_final]-2) < 0.35 #@test abs(sim.𝒪est[:weak_l2]-2) < 0.3 #@test abs(sim.𝒪est[:weak_l∞]-2) < 0.3