OrdinaryDiffEq.jl icon indicating copy to clipboard operation
OrdinaryDiffEq.jl copied to clipboard

Potentially missing update of `fsalfirst`

Open baggepinnen opened this issue 2 years ago • 0 comments

The following code makes use of the integrator interface to take a single step of length dt from a state set using set_u!. If the procedure is repeated twice, the results differ.

using OrdinaryDiffEq, Test
f = (u,p,t) -> -u
u0 = [0.0] # Throwaway value
dt = 1
prob = ODEProblem(f, u0, (0, Inf))
int = init(prob, RK4(); dt, save_everystep=false, adaptive=false) # adaptive=false to minimize complexity

u0_2 = [1.0]
OrdinaryDiffEq.SciMLBase.set_ut!(int, u0_2, 0) # Start integrating from u != u0 above
OrdinaryDiffEq.step!(int, dt, true)
u1_1 = copy(int.u)

OrdinaryDiffEq.SciMLBase.set_ut!(int, u0_2, 0)
OrdinaryDiffEq.step!(int, dt, true)
u1_2 = copy(int.u)

@test u1_1 == u1_2

Test Failed at /home/fredrikb/.julia/dev/JuliaSimControls/test/MPC/test_integrators.jl:220
  Expression: u1_1 == u1_2
   Evaluated: [0.41666666666666674] == [0.375]
ERROR: There was an error during testing

If one manually sets u,uprev and inserts a call to initialize!, the results are identical

prob = ODEProblem(f, u0, (0, Inf))
int = init(prob, RK4(); dt, save_everystep=false, adaptive=false) # adaptive=false to minimize complexity

OrdinaryDiffEq.SciMLBase.set_ut!(int, u0_2, 0) # Start integrating from u != u0 above
int.u = int.uprev = u0_2
OrdinaryDiffEq.initialize!(int, int.cache) # Force update of fsalfirst
OrdinaryDiffEq.step!(int, dt, true)
u1_1 = copy(int.u)


OrdinaryDiffEq.SciMLBase.set_ut!(int, u0_2, 0)
OrdinaryDiffEq.step!(int, dt, true)
u1_2 = copy(int.u)

@test u1_1 == u1_2

julia> @test u1_1 == u1_2
Test Passed
  Expression: u1_1 == u1_2
   Evaluated: [0.375] == [0.375]

baggepinnen avatar Mar 29 '22 10:03 baggepinnen