Reinit! doesn't work on MTK integrators

Specific Domains Modelling & Simulations
, ,
1

It seems like reinit! doesn’t work for ModelingToolkit systems. Is this true, what is the workaround and why don’t I get an error or warning? Trying to run this mwe shows that the integrator state is unchanged after reinit!.

using OrdinaryDiffEq
using LinearAlgebra
using ModelingToolkit

# Create a simple ODE system
function simple_system()
    @independent_variables t
    D = Differential(t)
    
    @variables x(t) y(t) z(t)
    @parameters p
    
    eqs = [
        D(x) ~ y + p,
        D(y) ~ z,
        D(z) ~ -x - 0.1*y
    ]
    
    @named sys = ODESystem(eqs, t, [x, y, z], [p])
    sys = structural_simplify(sys)
    return sys
end

function demonstrate_issue()
    sys = simple_system()
    
    # Create two different initial states
    initial_u0 = [1.0, 0.0, 0.0]
    new_u0 = [5.0, 2.0, 3.0]
    
    # Create the problem with initial u0
    prob = ODEProblem(sys, initial_u0, (0.0, 10.0), [0.5])
    
    # Initialize with FBDF solver
    println("Testing with FBDF solver:")
    fbdf_integrator = init(prob, FBDF(); dt=0.1, abstol=1e-6, reltol=1e-6)
    
    println("Initial state: ", fbdf_integrator.u)
    
    # Try to reinitialize with new_u0
    reinit!(fbdf_integrator, new_u0)
    
    println("After reinit!: ", fbdf_integrator.u)
    println("Expected:      ", new_u0)
    println("Are equal? ", fbdf_integrator.u == new_u0)
    
    # Compare with Tsit5 solver
    println("\nTesting with Tsit5 solver:")
    tsit_integrator = init(prob, Tsit5(); dt=0.1, abstol=1e-6, reltol=1e-6)
    
    println("Initial state: ", tsit_integrator.u)
    
    # Try to reinitialize with new_u0
    reinit!(tsit_integrator, new_u0)
    
    println("After reinit!: ", tsit_integrator.u)
    println("Expected:      ", new_u0)
    println("Are equal? ", tsit_integrator.u == new_u0)
    
    # Reset cache and history
    println("\nTesting FBDF with reset cache:")
    reinit!(fbdf_integrator, new_u0; reset_dt=true, reinit_cache=true, erase_sol=true)
    
    println("After reset cache: ", fbdf_integrator.u)
    println("Expected:         ", new_u0)
    println("Are equal? ", fbdf_integrator.u == new_u0)
end

demonstrate_issue()

Gives the output:

julia> include("mwe/reinit_problem.jl")
Initial state: [1.0, 0.0, 0.0]
After reinit!: [1.0, 0.0, 0.0]
Expected:      [5.0, 2.0, 3.0]
Are equal? false

Testing with Tsit5 solver:
Initial state: [1.0, 0.0, 0.0]
After reinit!: [1.0, 0.0, 0.0]
Expected:      [5.0, 2.0, 3.0]
Are equal? false

Testing FBDF with reset cache:
After reset cache: [1.0, 0.0, 0.0]
Expected:         [5.0, 2.0, 3.0]
Are equal? false
2

Using setu works, but still, it would be nice to get some kind of warning when trying to use reinit!.

3

I agree with you, this is a problem.

1 Like
4

handled now.