Skip to content

Commit e226e9a

Browse files
devmotiondevmotion
committed
Merge master and remove calls of ode_addsteps!
2 parents 73ee9f5 + 9223d96 commit e226e9a

File tree

11 files changed

+637
-126
lines changed

11 files changed

+637
-126
lines changed

‎src/DelayDiffEq.jl‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -18,6 +18,7 @@ import OrdinaryDiffEq: Rosenbrock23Cache, Rosenbrock32Cache, ImplicitEulerCache,
1818
TrapezoidCache
1919

2020
include("integrator_type.jl")
21+
include("integrator_utils.jl")
2122
include("integrator_interface.jl")
2223
include("history_function.jl")
2324
include("algorithms.jl")

‎src/integrator_interface.jl‎

Lines changed: 17 additions & 44 deletions
Original file line numberDiff line numberDiff line change
@@ -4,29 +4,17 @@
44
Update solution of `integrator`, if necessary or forced by `force_save`.
55
"""
66
function savevalues!(integrator::DDEIntegrator, force_save=false)
7-
# update ODE integrator to interval [tprev, t] with corresponding values
8-
# u(tprev) and u(t), and interpolation data k of this interval
9-
integrator.integrator.t = integrator.t
10-
integrator.integrator.tprev = integrator.tprev
11-
12-
# copy u(tprev) since it is overwritten by integrator at the end of apply_step!
13-
if typeof(integrator.u) <: AbstractArray
14-
recursivecopy!(integrator.integrator.u, integrator.u)
15-
recursivecopy!(integrator.integrator.uprev, integrator.uprev)
16-
else
17-
integrator.integrator.u = integrator.u
18-
integrator.integrator.uprev = integrator.uprev
19-
end
20-
21-
# copy interpolation data (fsalfirst overwritten at the end of apply_step!, which also
22-
# updates k[1] when using chaches for which k[1] points to fsalfirst)
23-
recursivecopy!(integrator.integrator.k, integrator.k)
24-
25-
# add steps for interpolation to ODE integrator when needed
26-
OrdinaryDiffEq.ode_addsteps!(integrator.integrator, integrator.f)
7+
# update ODE integrator
8+
update_ode_integrator!(integrator)
279

2810
# update solution of ODE integrator
2911
savevalues!(integrator.integrator, force_save)
12+
13+
# delete part of ODE solution that is not required for DDE solution
14+
reduce_solution!(integrator,
15+
# function values at later time points might be necessary for
16+
# calculation of next step, thus keep those interpolation data
17+
integrator.integrator.tprev - maximum(integrator.prob.lags))
3018
end
3119

3220
"""
@@ -36,19 +24,13 @@ Clean up solution of `integrator`.
3624
"""
3725
function postamble!(integrator::DDEIntegrator)
3826
# update ODE integrator
39-
integrator.integrator.t = integrator.t
40-
if typeof(integrator.u) <: AbstractArray
41-
recursivecopy!(integrator.integrator.u, integrator.u)
42-
else
43-
integrator.integrator.u = integrator.u
44-
end
45-
recursivecopy!(integrator.integrator.k, integrator.k)
46-
47-
# add steps for interpolation to ODE integrator when needed
48-
OrdinaryDiffEq.ode_addsteps!(integrator.integrator, integrator.f)
27+
update_ode_integrator!(integrator)
4928

5029
# clean up solution of ODE integrator
5130
OrdinaryDiffEq.postamble!(integrator.integrator)
31+
32+
# delete part of ODE solution that is not required for DDE solution
33+
reduce_solution!(integrator, integrator.integrator.sol.t[end])
5234
end
5335

5436
"""
@@ -64,15 +46,9 @@ function perform_step!(integrator::DDEIntegrator)
6446
recursivecopy!(integrator.k_cache, integrator.k)
6547
integrator.integrator.EEst = integrator.EEst
6648

67-
# add steps to interpolation data ODE integrator if necessary
68-
OrdinaryDiffEq.ode_addsteps!(integrator.integrator, integrator.f)
69-
7049
# perform always at least one calculation
7150
perform_step!(integrator, integrator.cache)
7251

73-
# shrink interpolation data of ODE problem
74-
resize!(integrator.integrator.k, integrator.integrator.kshortsize)
75-
7652
# if dt is greater than the minimal lag, then use a fixed-point iteration
7753
if integrator.dt > minimum(integrator.prob.lags) && isfinite(integrator.EEst)
7854

@@ -92,6 +68,7 @@ function perform_step!(integrator::DDEIntegrator)
9268
# in the next iterations when evaluating the history function
9369
integrator.integrator.t = integrator.t + integrator.dt
9470
integrator.integrator.tprev = integrator.t
71+
integrator.integrator.dt = integrator.dt
9572
if typeof(integrator.integrator.uprev) <: AbstractArray
9673
recursivecopy!(integrator.integrator.uprev, integrator.uprev)
9774
else
@@ -110,17 +87,12 @@ function perform_step!(integrator::DDEIntegrator)
11087
end
11188
recursivecopy!(integrator.integrator.k, integrator.k)
11289

113-
# add steps to interpolation data of ODE integrator if necessary
114-
OrdinaryDiffEq.ode_addsteps!(integrator.integrator, integrator.f)
115-
11690
# calculate next step
11791
perform_step!(integrator, integrator.cache)
11892

119-
# shrink interpolation data of ODE integrator
120-
resize!(integrator.integrator.k, integrator.integrator.kshortsize)
121-
12293
# calculate residuals of fixed-point iteration
123-
if typeof(integrator.resid) <: AbstractArray
94+
# can be fixed with new @muladd: https://github.com/JuliaDiffEq/DiffEqBase.jl/pull/57
95+
if typeof(integrator.u) <: AbstractArray
12496
@tight_loop_macros @inbounds for (i, atol, rtol) in
12597
zip(eachindex(integrator.u),
12698
Iterators.cycle(integrator.fixedpoint_abstol),
@@ -140,7 +112,7 @@ function perform_step!(integrator::DDEIntegrator)
140112

141113
# stop fixed-point iteration when error estimate of integrator or error estimate
142114
# of fixed-point iteration are not finite
143-
if !(isfinite(fixedpointEEst)) || !(isfinite(integrator.EEst))
115+
if !isfinite(fixedpointEEst) || !isfinite(integrator.EEst)
144116
# assure that integrator is reset to cached values
145117
integrator.EEst = max(fixedpointEEst, integrator.EEst)
146118
break
@@ -167,6 +139,7 @@ function perform_step!(integrator::DDEIntegrator)
167139
# u(tprev) and u(t), and interpolation data k of this interval
168140
integrator.integrator.t = integrator.t
169141
integrator.integrator.tprev = integrator.tprev
142+
integrator.integrator.dt = integrator.integrator.t - integrator.integrator.tprev
170143
if typeof(integrator.u) <: AbstractArray
171144
recursivecopy!(integrator.integrator.u, integrator.uprev)
172145
recursivecopy!(integrator.integrator.uprev, integrator.uprev_cache)

‎src/integrator_type.jl‎

Lines changed: 15 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,6 @@
11
mutable struct DDEIntegrator{algType<:OrdinaryDiffEqAlgorithm,uType,tType,absType,relType,
22
residType,tTypeNoUnits,tdirType,ksEltype,SolType,rateType,F,
3-
ProgressType,CacheType,IType,ProbType,NType,O} <:
3+
ProgressType,CacheType,IType,ProbType,NType,O,tstopsType} <:
44
AbstractDDEIntegrator
55

66
sol::SolType
@@ -20,6 +20,7 @@ mutable struct DDEIntegrator{algType<:OrdinaryDiffEqAlgorithm,uType,tType,absTyp
2020
resid::residType # This would have to resize for resizing DDE to work
2121
fixedpoint_norm::NType
2222
max_fixedpoint_iters::Int
23+
minimal_solution::Bool
2324
alg::algType
2425
rate_prototype::rateType
2526
notsaveat_idxs::Vector{Int}
@@ -43,29 +44,31 @@ mutable struct DDEIntegrator{algType<:OrdinaryDiffEqAlgorithm,uType,tType,absTyp
4344
u_modified::Bool
4445
opts::O
4546
integrator::IType
47+
saveat::tstopsType
4648
fsalfirst::rateType
4749
fsallast::rateType
4850

4951
# incomplete initialization without fsalfirst and fsallast
5052
function DDEIntegrator{algType,uType,tType,absType,relType,residType,tTypeNoUnits,
5153
tdirType,ksEltype,SolType,rateType,F,ProgressType,CacheType,
52-
IType,ProbType,NType,O}(
54+
IType,ProbType,NType,O,tstopsType}(
5355
sol,prob,u,k,t,dt,f,uprev,tprev,uprev_cache,k_cache,
54-
k_integrator_cache,fixedpoint_abstol,fixedpoint_reltol,resid,
55-
fixedpoint_norm,max_fixedpoint_iters,alg,rate_prototype,
56-
notsaveat_idxs,dtcache,dtchangeable,dtpropose,tdir,EEst,qold,
57-
q11,iter,saveiter,saveiter_dense,prog,cache,kshortsize,
58-
just_hit_tstop,accept_step,isout,reeval_fsal,u_modified,opts,
59-
integrator) where
56+
k_integrator_cache,fixedpoint_abstol,fixedpoint_reltol,resid,
57+
fixedpoint_norm,max_fixedpoint_iters,minimal_solution,alg,
58+
rate_prototype,notsaveat_idxs,dtcache,dtchangeable,dtpropose,
59+
tdir,EEst,qold,q11,iter,saveiter,saveiter_dense,prog,cache,
60+
kshortsize,just_hit_tstop,accept_step,isout,reeval_fsal,
61+
u_modified,opts,integrator,saveat) where
6062
{algType<:OrdinaryDiffEqAlgorithm,uType,tType,absType,relType,residType,
6163
tTypeNoUnits,tdirType,ksEltype,SolType,rateType,F,ProgressType,CacheType,IType,
62-
ProbType,NType,O}
64+
ProbType,NType,O,tstopsType}
6365

6466
new(sol,prob,u,k,t,dt,f,uprev,tprev,uprev_cache,k_cache,k_integrator_cache,
6567
fixedpoint_abstol,fixedpoint_reltol,resid,fixedpoint_norm,max_fixedpoint_iters,
66-
alg,rate_prototype,notsaveat_idxs,dtcache,dtchangeable,dtpropose,tdir,EEst,qold,
67-
q11,iter,saveiter,saveiter_dense,prog,cache,kshortsize,just_hit_tstop,
68-
accept_step,isout,reeval_fsal,u_modified,opts,integrator)
68+
minimal_solution,alg,rate_prototype,notsaveat_idxs,dtcache,dtchangeable,
69+
dtpropose,tdir,EEst,qold,q11,iter,saveiter,saveiter_dense,prog,cache,
70+
kshortsize,just_hit_tstop,accept_step,isout,reeval_fsal,u_modified,opts,
71+
integrator,saveat)
6972
end
7073
end
7174

0 commit comments

Comments
 (0)