SciML / Sundials.jl

Julia interface to Sundials, including a nonlinear solver (KINSOL), ODE's (CVODE and ARKODE), and DAE's (IDA) in a SciML scientific machine learning enabled manner
https://diffeq.sciml.ai
BSD 2-Clause "Simplified" License
209 stars 78 forks source link

Option stop_at_tdt not correctly documented #307

Open ufechner7 opened 3 years ago

ufechner7 commented 3 years ago

The following code gives an unexpected output:

differential_vars =  ones(Bool, 36)
solver = IDA(linear_solver=:Dense)
dt = 0.5
t_start = 0.0
t_end   = 5*dt
tspan = (t_start, t_end) 

prob = DAEProblem(residual!, yd0, y0, tspan, differential_vars=differential_vars)
integrator = init(prob, solver, abstol=0.000001, reltol=0.001)

for i in 1:5
    step!(integrator, dt, true)
    for (u,t) in tuples(integrator)
        @show u[18], t
    end
end

Output:

(u[18], t) = (368.74700170881215, 0.6359999999999999)
(u[18], t) = (368.7470017992227, 0.7719999999999998)
(u[18], t) = (368.7470020404195, 1.0439999999999996)
(u[18], t) = (368.74700225092147, 1.3159999999999994)
(u[18], t) = (368.74700290175184, 1.859999999999999)
(u[18], t) = (368.7470033622797, 2.4039999999999986)
(u[18], t) = (368.74700352032835, 2.5)

I would expect the results at 0.5, 1, 1.5, 2.0 and 2.5 seconds simulation time.

Full example: https://github.com/ufechner7/KiteViewer/blob/sim/src/RTSim_bug.jl

ufechner7 commented 3 years ago

This code is working:

for i in 1:round(t_end/dt)
    step!(integrator, dt, true)
    u = integrator.u
    t = integrator.t
    @show u[18], t
end
ufechner7 commented 3 years ago

Lets say only the documentation is wrong. It says:

step!(integrator,dt[,stop_at_tdt=false])

at https://diffeq.sciml.ai/stable/basics/integrator/#SciMLBase.step!

ufechner7 commented 3 years ago

.

Perhaps it would also be good to add the working example above to the documentation. Because if it is not added the only thing you see is:

for (u,t) in tuples(integrator)
  @show u,t
end

Which does not work together with the option stop_at_tdt=true.

ChrisRackauckas commented 3 years ago

What do you mean that doesn't work?

ufechner7 commented 3 years ago

Look at the initial example. It gives wrong output, it does not respect dt.

ChrisRackauckas commented 3 years ago

In

for i in 1:5
    step!(integrator, dt, true)
    for (u,t) in tuples(integrator)
        @show u[18], t
    end
end

You're explicitly telling it to only use dt in the first step, and then after you've solved to the end of the interval using adaptive time stepping, use dt 4 more times.

ChrisRackauckas commented 3 years ago
using OrdinaryDiffEq
prob = ODEProblem((u,p,t)->1.01u,1.01,(0.0,1.0))
integrator = init(prob,Tsit5())
for (u,t) in tuples(integrator)
  @show u,t
end
(u, t) = (1.1169134682731792, 0.09962249330449434)
(u, t) = (1.4319565770427753, 0.34563506464944166)
(u, t) = (2.002298156149145, 0.6775695999548759)
(u, t) = (2.7730568905500066, 1.0)

It iterates as documented, and because it does what's documented it causes your "issue".

ufechner7 commented 3 years ago

That is very confusing. Why should the command:

tuples(integrator)

actually cause the integration to continue? My understanding of this command is that it just collects the result of the last integration.

ChrisRackauckas commented 3 years ago

It's an iterator that iterates by stepping.

using OrdinaryDiffEq
prob = ODEProblem((u,p,t)->1.01u,1.01,(0.0,1.0))
integrator = init(prob,Tsit5())
for integ in integrator
  @show integ.t
end
integ.t = 0.09962249330449434
integ.t = 0.34563506464944166
integ.t = 0.6775695999548759
integ.t = 1.0

It's the standard Julia construct for stepping. I don't understand why you were using stepping to view terms when you didn't want stepping.

ChrisRackauckas commented 3 years ago

My understand of this command is that it just collects the result of the last integration.

No, it's an iterator, defined in the iteration section as an iterator.

ChrisRackauckas commented 3 years ago
This type also implements an iterator interface, so one can step n times (or to the last tstop) using the take iterator:

for i in take(integrator,n) end
One can loop to the end by using solve!(integrator) or using the iterator interface:

for i in integrator end
In addition, some helper iterators are provided to help monitor the solution. For example, the tuples iterator lets you view the values:

for (u,t) in tuples(integrator)
  @show u,t
end
and the intervals iterator lets you view the full interval

etc. It's all about different iterators to control and view stepping behavior in nice ways.

ufechner7 commented 3 years ago

Anyway, I think it would be helpful to add an example like this:

dt = 0.1
t_en = 10.0
for i in 1:round(t_end/dt)
    step!(integrator, dt, true)
    u = integrator.u
    t = integrator.t
    @show u, t
end

to the documentation.

ChrisRackauckas commented 3 years ago

Agreed.

ChrisRackauckas commented 3 years ago

TimeChoiceIterator might be a little more natural though?

ChrisRackauckas commented 3 years ago
for (u,t) in TimeChoiceIterator(integrator,0:dt:tend)
  @show u,t
end

might be the nicest way.

ufechner7 commented 3 years ago

But does this do the integration step-by-step with the possibility to update parameters after each step?

ChrisRackauckas commented 3 years ago

Yes, it iterates to the next t and gives you control in the middle of the loop.

for (u,t) in TimeChoiceIterator(integrator,0:dt:tend)
  @show u,t
  integrator.p = ... # modify the parameter before continuing
end
ufechner7 commented 3 years ago

Great!