`step!` fails to build a solution sometimes (i.e. integrator and solution time can fall out of sync)

[bradcarman]

The below MWE shows that running 24 steps of the integrator using `step!(integrator, dt, true) gives a solution with 25 rows and the solution end time matches the integrator time. But when running by 25 steps, the solution does not build for the final step, leaving the integrator and solution time out of sync.

using Test
using OrdinaryDiffEq

dt = 4e-4
prob = ODEProblem((u, p, t) -> -u, 1.0, nothing)
integrator = init(prob, ImplicitEuler(); tspan=(0., 1.0), adaptive=false, dt=dt,  initializealg=NoInit(), verbose=true, saveat=dt)
[step!(integrator, dt, true) for i=1:24];
@test integrator.sol.t[end] ≈ integrator.t

SciMLBase.set_ut!(integrator, copy(integrator.sol(0)), 0)
[step!(integrator, dt, true) for i=1:25];
@test integrator.sol.t[end] ≈ integrator.t

[YingboMa]

The problem is caused by

julia> dt = 4e-4;

julia> sdt = 0.0; for i in 1:25
           sdt += dt
           (0:dt:1)[i+1] <= sdt || println("$i * dt != sdt")
       end
25 * dt != sdt

just use

integrator = init(prob, ImplicitEuler(); tspan=(0., 1.0), adaptive=false, dt=dt,  initializealg=NoInit(), verbose=true)
step!(integrator, dt)

because adaptive = false already forces the step size to be dt.

[bradcarman]

The bug is not about the resolution of the time. I understand that 25*dt != sdt, but what is true is 25*dt ≈ sdt.

What I'm seeing is the solution is not updated, so there should be 26 rows, but there are only 25!

using Test
using OrdinaryDiffEq

dt = 4e-4
prob = ODEProblem((u, p, t) -> -u, 1.0, nothing)
integrator = init(prob, ImplicitEuler(); tspan=(0., 1.0), adaptive=false, dt=dt,  initializealg=NoInit(), verbose=true, saveat=dt)
[step!(integrator, dt) for i=1:24];
@test size(integrator.sol,1) == 24+1

SciMLBase.set_ut!(integrator, copy(integrator.sol(0)), 0)
[step!(integrator, dt) for i=1:25];
@test size(integrator.sol,1) == 25+1

This will pass the first test, but fail the second

Test Failed 
  Expression: size(integrator.sol, 1) == 25 + 1
   Evaluated: 25 == 26
ERROR: There was an error during testing

[YingboMa]

You have saveat=dt, it rejects to save at sdt because it's not at saveat

[bradcarman]

saveat=dt is a requirement because I can't reset the integrator without it.

julia> SciMLBase.set_ut!(integrator, copy(integrator.sol(0)), 0)
ERROR: Integrator state cannot be reset unless it is initialized with save_everystep=false

[bradcarman]

For my real system I'm attempting to run a step and extract the observables, therefore I need a solution object to do that. For the MWE of course the solution object is not necessary because there is only one variable. I need a way to step the integrator and get the solution for every step.

[YingboMa]

You could use the iterator interface to save the solution manually

julia> integrator = init(prob, ImplicitEuler(); tspan=(0., 1.0), save_on=false, save_everystep=false, verbose=true);

julia> for (u,t) in tuples(integrator)
           # do something
       end

julia> SciMLBase.set_ut!(integrator, copy(prob.u0), 0)
retcode: Default
Interpolation: 1st order linear
t: 1-element Vector{Float64}:
 0.0
u: 1-element Vector{Float64}:
 1.0

julia> for (u,t) in tuples(integrator)
           # do something
       end

Then you can use prob.f.observed to get the observed value.

[bradcarman]

OK, that works, I'm able to do the following which allows me to get the solved variables (integrator.u) and observables from integrator.f.observed() like so...

step!(integrator, dt, true)    
integrator.f.observed(num, integrator.u, integrator.p, integrator.t)