Skip to content

Commit b7579ea

Browse files
ChrisRackauckasChrisRackauckas
authored
Merge pull request #67 from JuliaDiffEq/double_implicit
Noise implicit methods
2 parents 1ef7813 + 68433d7 commit b7579ea

10 files changed

+620
-7
lines changed

‎src/StochasticDiffEq.jl

Lines changed: 4 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -49,6 +49,7 @@ module StochasticDiffEq
4949
include("caches/lamba_caches.jl")
5050
include("caches/iif_caches.jl")
5151
include("caches/sdirk_caches.jl")
52+
include("caches/split_step_caches.jl")
5253
include("caches/sra_caches.jl")
5354
include("caches/rossler_caches.jl")
5455
include("caches/kencarp_caches.jl")
@@ -68,6 +69,7 @@ module StochasticDiffEq
6869
include("perform_step/sri.jl")
6970
include("perform_step/sra.jl")
7071
include("perform_step/sdirk.jl")
72+
include("perform_step/split_step.jl")
7173
include("perform_step/kencarp.jl")
7274
include("perform_step/predcorr.jl")
7375
include("perform_step/split.jl")
@@ -86,7 +88,8 @@ module StochasticDiffEq
8688

8789
export SplitEM, IIF1M, IIF2M, IIF1Mil
8890

89-
export ImplicitEM, ImplicitEulerHeun, ImplicitRKMil
91+
export ImplicitEM, ImplicitEulerHeun, ISSEM, ISSEulerHeun,
92+
ImplicitRKMil
9093

9194
export StochasticDiffEqRODEAlgorithm, StochasticDiffEqRODEAdaptiveAlgorithm,
9295
StochasticDiffEqRODECompositeAlgorithm

‎src/alg_utils.jl

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -10,6 +10,8 @@ alg_order(alg::LambaEM) = 1//2
1010
alg_order(alg::ImplicitEM) = 1//2
1111
alg_order(alg::ImplicitEulerHeun) = 1//2
1212
alg_order(alg::ImplicitRKMil) = 1//1
13+
alg_order(alg::ISSEM) = 1//2
14+
alg_order(alg::ISSEulerHeun) = 1//2
1315
alg_order(alg::SplitEM) = 1//2
1416
alg_order(alg::PCEuler) = 1//2
1517
alg_order(alg::IIF1M) = 1//2
@@ -69,6 +71,8 @@ alg_compatible(prob,alg::SplitEM) = true
6971
alg_compatible(prob,alg::PCEuler) = true
7072
alg_compatible(prob,alg::ImplicitEM) = true
7173
alg_compatible(prob,alg::ImplicitEulerHeun) = true
74+
alg_compatible(prob,alg::ISSEM) = true
75+
alg_compatible(prob,alg::ISSEulerHeun) = true
7276
alg_compatible(prob,alg::RKMil) = is_diagonal_noise(prob)
7377
alg_compatible(prob,alg::ImplicitRKMil) = is_diagonal_noise(prob)
7478
alg_compatible(prob,alg::RKMilCommute) = true # No good check for commutative noise

‎src/algorithms.jl

Lines changed: 65 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -115,8 +115,10 @@ Base.@pure ImplicitEM(;chunk_size=0,autodiff=true,diff_type=Val{:central},
115115
theta = 1/2,symplectic=false,
116116
max_newton_iter=7,new_jac_conv_bound = 1e-3,
117117
controller = :Predictive) =
118-
ImplicitEM{chunk_size,autodiff,typeof(linsolve),typeof(diff_type),
119-
typeof(κ),typeof(tol),typeof(new_jac_conv_bound),controller}(
118+
ImplicitEM{chunk_size,autodiff,
119+
typeof(linsolve),typeof(diff_type),
120+
typeof(κ),typeof(tol),
121+
typeof(new_jac_conv_bound),controller}(
120122
linsolve,diff_type,κ,tol,
121123
symplectic ? 1/2 : theta,
122124
extrapolant,
@@ -141,8 +143,10 @@ Base.@pure ImplicitEulerHeun(;chunk_size=0,autodiff=true,diff_type=Val{:central}
141143
theta = 1/2,symplectic = false,
142144
max_newton_iter=7,new_jac_conv_bound = 1e-3,
143145
controller = :Predictive) =
144-
ImplicitEulerHeun{chunk_size,autodiff,typeof(linsolve),typeof(diff_type),
145-
typeof(κ),typeof(tol),typeof(new_jac_conv_bound),controller}(
146+
ImplicitEulerHeun{chunk_size,autodiff,
147+
typeof(linsolve),typeof(diff_type),
148+
typeof(κ),typeof(tol),
149+
typeof(new_jac_conv_bound),controller}(
146150
linsolve,diff_type,κ,tol,
147151
symplectic ? 1/2 : theta,
148152
extrapolant,min_newton_iter,
@@ -166,14 +170,70 @@ Base.@pure ImplicitRKMil(;chunk_size=0,autodiff=true,diff_type=Val{:central},
166170
theta = 1/2,symplectic = false,
167171
max_newton_iter=7,new_jac_conv_bound = 1e-3,
168172
controller = :Predictive,interpretation=:Ito) =
169-
ImplicitRKMil{chunk_size,autodiff,typeof(linsolve),typeof(diff_type),
173+
ImplicitRKMil{chunk_size,autodiff,
174+
typeof(linsolve),typeof(diff_type),
170175
typeof(κ),typeof(tol),typeof(new_jac_conv_bound),
171176
controller,interpretation}(
172177
linsolve,diff_type,κ,tol,
173178
symplectic ? 1/2 : theta,
174179
extrapolant,min_newton_iter,
175180
max_newton_iter,new_jac_conv_bound,symplectic)
176181

182+
struct ISSEM{CS,AD,F,S,K,T,T2,Controller} <: StochasticDiffEqNewtonAdaptiveAlgorithm{CS,AD,Controller}
183+
linsolve::F
184+
diff_type::S
185+
κ::K
186+
tol::T
187+
theta::T2
188+
extrapolant::Symbol
189+
min_newton_iter::Int
190+
max_newton_iter::Int
191+
new_jac_conv_bound::T2
192+
symplectic::Bool
193+
end
194+
Base.@pure ISSEM(;chunk_size=0,autodiff=true,diff_type=Val{:central},
195+
linsolve=DEFAULT_LINSOLVE,κ=nothing,tol=nothing,
196+
extrapolant=:constant,min_newton_iter=1,
197+
theta = 1,symplectic=false,
198+
max_newton_iter=7,new_jac_conv_bound = 1e-3,
199+
controller = :Predictive) =
200+
ISSEM{chunk_size,autodiff,
201+
typeof(linsolve),typeof(diff_type),
202+
typeof(κ),typeof(tol),
203+
typeof(new_jac_conv_bound),controller}(
204+
linsolve,diff_type,κ,tol,
205+
symplectic ? 1/2 : theta,
206+
extrapolant,
207+
min_newton_iter,
208+
max_newton_iter,new_jac_conv_bound,symplectic)
209+
210+
struct ISSEulerHeun{CS,AD,F,S,K,T,T2,Controller} <: StochasticDiffEqNewtonAdaptiveAlgorithm{CS,AD,Controller}
211+
linsolve::F
212+
diff_type::S
213+
κ::K
214+
tol::T
215+
theta::T2
216+
extrapolant::Symbol
217+
min_newton_iter::Int
218+
max_newton_iter::Int
219+
new_jac_conv_bound::T2
220+
symplectic::Bool
221+
end
222+
Base.@pure ISSEulerHeun(;chunk_size=0,autodiff=true,diff_type=Val{:central},
223+
linsolve=DEFAULT_LINSOLVE,κ=nothing,tol=nothing,
224+
extrapolant=:constant,min_newton_iter=1,
225+
theta = 1,symplectic=false,
226+
max_newton_iter=7,new_jac_conv_bound = 1e-3,
227+
controller = :Predictive) =
228+
ISSEulerHeun{chunk_size,autodiff,
229+
typeof(linsolve),typeof(diff_type),
230+
typeof(κ),typeof(tol),typeof(new_jac_conv_bound),controller}(
231+
linsolve,diff_type,κ,tol,
232+
symplectic ? 1/2 : theta,
233+
extrapolant,
234+
min_newton_iter,
235+
max_newton_iter,new_jac_conv_bound,symplectic)
236+
177237
struct SKenCarp{CS,AD,F,FDT,K,T,T2,Controller} <: StochasticDiffEqNewtonAdaptiveAlgorithm{CS,AD,Controller}
178238
linsolve::F
179239
diff_type::FDT

‎src/caches/split_step_caches.jl

Lines changed: 182 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,182 @@
1+
mutable struct ISSEMCache{uType,rateType,J,JC,UF,uEltypeNoUnits,noiseRateType,F} <: StochasticDiffEqMutableCache
2+
u::uType
3+
uprev::uType
4+
du1::rateType
5+
fsalfirst::rateType
6+
k::rateType
7+
z::uType
8+
dz::uType
9+
tmp::uType
10+
gtmp::noiseRateType
11+
gtmp2::rateType
12+
J::J
13+
W::J
14+
jac_config::JC
15+
linsolve::F
16+
uf::UF
17+
ηold::uEltypeNoUnits
18+
κ::uEltypeNoUnits
19+
tol::uEltypeNoUnits
20+
newton_iters::Int
21+
end
22+
23+
u_cache(c::ISSEMCache) = (c.uprev2,c.z,c.dz)
24+
du_cache(c::ISSEMCache) = (c.k,c.fsalfirst)
25+
26+
function alg_cache(alg::ISSEM,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,
27+
uEltypeNoUnits,uBottomEltype,tTypeNoUnits,uprev,f,t,::Type{Val{true}})
28+
du1 = zeros(rate_prototype)
29+
J = zeros(uEltypeNoUnits,length(u),length(u)) # uEltype?
30+
W = similar(J)
31+
z = similar(u)
32+
dz = similar(u); tmp = similar(u); gtmp = similar(noise_rate_prototype)
33+
fsalfirst = zeros(rate_prototype)
34+
k = zeros(rate_prototype)
35+
36+
uf = DiffEqDiffTools.UJacobianWrapper(f,t,p)
37+
linsolve = alg.linsolve(Val{:init},uf,u)
38+
jac_config = build_jac_config(alg,f,uf,du1,uprev,u,tmp,dz)
39+
ηold = one(uEltypeNoUnits)
40+
41+
if alg.κ != nothing
42+
κ = alg.κ
43+
else
44+
κ = uEltypeNoUnits(1//100)
45+
end
46+
if alg.tol != nothing
47+
tol = alg.tol
48+
else
49+
reltol = 1e-1 # TODO: generalize
50+
tol = min(0.03,first(reltol)^(0.5))
51+
end
52+
53+
if is_diagonal_noise(prob)
54+
gtmp2 = gtmp
55+
else
56+
gtmp2 = similar(rate_prototype)
57+
end
58+
59+
ISSEMCache(u,uprev,du1,fsalfirst,k,z,dz,tmp,gtmp,gtmp2,J,W,jac_config,linsolve,uf,
60+
ηold,κ,tol,10000)
61+
end
62+
63+
mutable struct ISSEMConstantCache{F,uEltypeNoUnits} <: StochasticDiffEqConstantCache
64+
uf::F
65+
ηold::uEltypeNoUnits
66+
κ::uEltypeNoUnits
67+
tol::uEltypeNoUnits
68+
newton_iters::Int
69+
end
70+
71+
function alg_cache(alg::ISSEM,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,
72+
uEltypeNoUnits,uBottomEltype,tTypeNoUnits,uprev,f,t,::Type{Val{false}})
73+
uf = DiffEqDiffTools.UDerivativeWrapper(f,t,p)
74+
ηold = one(uEltypeNoUnits)
75+
76+
if alg.κ != nothing
77+
κ = alg.κ
78+
else
79+
κ = uEltypeNoUnits(1//100)
80+
end
81+
if alg.tol != nothing
82+
tol = alg.tol
83+
else
84+
reltol = 1e-1 # TODO: generalize
85+
tol = min(0.03,first(reltol)^(0.5))
86+
end
87+
88+
ISSEMConstantCache(uf,ηold,κ,tol,100000)
89+
end
90+
91+
mutable struct ISSEulerHeunCache{uType,rateType,J,JC,UF,uEltypeNoUnits,noiseRateType,F} <: StochasticDiffEqMutableCache
92+
u::uType
93+
uprev::uType
94+
du1::rateType
95+
fsalfirst::rateType
96+
k::rateType
97+
z::uType
98+
dz::uType
99+
tmp::uType
100+
gtmp::noiseRateType
101+
gtmp2::rateType
102+
gtmp3::noiseRateType
103+
J::J
104+
W::J
105+
jac_config::JC
106+
linsolve::F
107+
uf::UF
108+
ηold::uEltypeNoUnits
109+
κ::uEltypeNoUnits
110+
tol::uEltypeNoUnits
111+
newton_iters::Int
112+
end
113+
114+
u_cache(c::ISSEulerHeunCache) = (c.uprev2,c.z,c.dz)
115+
du_cache(c::ISSEulerHeunCache) = (c.k,c.fsalfirst)
116+
117+
function alg_cache(alg::ISSEulerHeun,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,
118+
uEltypeNoUnits,uBottomEltype,tTypeNoUnits,uprev,f,t,::Type{Val{true}})
119+
du1 = zeros(rate_prototype)
120+
J = zeros(uEltypeNoUnits,length(u),length(u)) # uEltype?
121+
W = similar(J)
122+
z = similar(u)
123+
dz = similar(u); tmp = similar(u); gtmp = similar(noise_rate_prototype)
124+
fsalfirst = zeros(rate_prototype)
125+
k = zeros(rate_prototype)
126+
127+
uf = DiffEqDiffTools.UJacobianWrapper(f,t,p)
128+
linsolve = alg.linsolve(Val{:init},uf,u)
129+
jac_config = build_jac_config(alg,f,uf,du1,uprev,u,tmp,dz)
130+
ηold = one(uEltypeNoUnits)
131+
132+
if alg.κ != nothing
133+
κ = alg.κ
134+
else
135+
κ = uEltypeNoUnits(1//100)
136+
end
137+
if alg.tol != nothing
138+
tol = alg.tol
139+
else
140+
reltol = 1e-1 # TODO: generalize
141+
tol = min(0.03,first(reltol)^(0.5))
142+
end
143+
144+
gtmp2 = similar(rate_prototype)
145+
146+
if is_diagonal_noise(prob)
147+
gtmp3 = gtmp2
148+
else
149+
gtmp3 = similar(noise_rate_prototype)
150+
end
151+
152+
ISSEulerHeunCache(u,uprev,du1,fsalfirst,k,z,dz,tmp,gtmp,gtmp2,gtmp3,
153+
J,W,jac_config,linsolve,uf,ηold,κ,tol,10000)
154+
end
155+
156+
mutable struct ISSEulerHeunConstantCache{F,uEltypeNoUnits} <: StochasticDiffEqConstantCache
157+
uf::F
158+
ηold::uEltypeNoUnits
159+
κ::uEltypeNoUnits
160+
tol::uEltypeNoUnits
161+
newton_iters::Int
162+
end
163+
164+
function alg_cache(alg::ISSEulerHeun,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,
165+
uEltypeNoUnits,uBottomEltype,tTypeNoUnits,uprev,f,t,::Type{Val{false}})
166+
uf = DiffEqDiffTools.UDerivativeWrapper(f,t,p)
167+
ηold = one(uEltypeNoUnits)
168+
169+
if alg.κ != nothing
170+
κ = alg.κ
171+
else
172+
κ = uEltypeNoUnits(1//100)
173+
end
174+
if alg.tol != nothing
175+
tol = alg.tol
176+
else
177+
reltol = 1e-1 # TODO: generalize
178+
tol = min(0.03,first(reltol)^(0.5))
179+
end
180+
181+
ISSEulerHeunConstantCache(uf,ηold,κ,tol,100000)
182+
end

‎src/perform_step/sdirk.jl

Lines changed: 0 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -111,7 +111,6 @@
111111

112112
cache.ηold = η
113113
cache.newton_iters = iter
114-
u = uprev + dt*(1-theta)*ftmp + theta*z + gtmp
115114

116115
if integrator.opts.adaptive
117116

0 commit comments

Comments
 (0)