Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Option stop_at_tdt not correctly documented #307

Open
ufechner7 opened this issue May 24, 2021 · 18 comments
Open

Option stop_at_tdt not correctly documented #307

ufechner7 opened this issue May 24, 2021 · 18 comments

Comments

@ufechner7
Copy link

ufechner7 commented May 24, 2021

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
Copy link
Author

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
Copy link
Author

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 ufechner7 changed the title Option stop_at_tdt = true not respected by function step! Option stop_at_tdt not correctly documented May 24, 2021
@ufechner7
Copy link
Author

ufechner7 commented May 24, 2021

.

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
Copy link
Member

What do you mean that doesn't work?

@ufechner7
Copy link
Author

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

@ChrisRackauckas
Copy link
Member

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
Copy link
Member

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
Copy link
Author

ufechner7 commented May 24, 2021

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
Copy link
Member

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
Copy link
Member

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
Copy link
Member

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
Copy link
Author

ufechner7 commented May 24, 2021

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
Copy link
Member

Agreed.

@ChrisRackauckas
Copy link
Member

TimeChoiceIterator might be a little more natural though?

@ChrisRackauckas
Copy link
Member

for (u,t) in TimeChoiceIterator(integrator,0:dt:tend)
  @show u,t
end

might be the nicest way.

@ufechner7
Copy link
Author

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

@ChrisRackauckas
Copy link
Member

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
Copy link
Author

Great!

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants