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

Supplying u0 to ODEProblem does not work in some cases #3268

Open
bradcarman opened this issue Dec 11, 2024 · 0 comments
Open

Supplying u0 to ODEProblem does not work in some cases #3268

bradcarman opened this issue Dec 11, 2024 · 0 comments
Labels
bug Something isn't working

Comments

@bradcarman
Copy link
Contributor

If I use InitializationProblem to pre-calculate a u0 for a simple problem, this works without incident...

@mtkmodel Simple begin
    @variables begin
        x(t) = 0
    end
    @equations begin
        D(x) ~ 1
    end
end

@mtkbuild sys = Simple()
initprob = ModelingToolkit.InitializationProblem(sys, 0.0)
initsol = solve(initprob)
u0 = initsol[unknowns(sys)]
prob = ODEProblem(sys, u0, (0, 0))  #OK, no warning

However, for a more complex problem, for some reason I do not get the same behavior...

function System()
    pars = @parameters begin
        rho_0 = 1000
        beta = 2e9
        A = 0.1
        m = 100
        L = 1
        p_s = 100e5
        p_r = 10e5
        C = 2.7*0.5
        c = 1000
        A_p = 0.00094
    end   
    vars = @variables begin
        p_1(t) = p_s
        p_2(t) = p_r
        x(t)=0
        dx(t)=0
        ddx(t), [guess=0]
        rho_1(t), [guess=rho_0]
        rho_2(t), [guess=rho_0]
        drho_1(t), [guess=0]
        drho_2(t), [guess=0]
        dm_1(t), [guess=0]
        dm_2(t), [guess=0]
    end
    
    # let -----
    u_1 = dm_1/(rho_0*A_p)
    u_2 = dm_2/(rho_0*A_p)

    eqs = [
        D(x) ~ dx
        D(dx) ~ ddx
        D(rho_1) ~ drho_1
        D(rho_2) ~ drho_2
        +dm_1 ~ drho_1*(L+x)*A + rho_1*dx*A
        -dm_2 ~ drho_2*(L-x)*A - rho_2*dx*A
        rho_1 ~ rho_0*(1 + p_1/beta)
        rho_2 ~ rho_0*(1 + p_2/beta)
        m*ddx ~ (p_1 - p_2)*A - c*dx
        (p_s - p_1) ~ C*rho_0*(u_1)*abs(u_1)
        (p_2 - p_r) ~ C*rho_0*(u_2)*abs(u_2)
    ]

    return ODESystem(eqs, t, vars, pars; name)
end

@mtkbuild sys = System()
initprob = ModelingToolkit.InitializationProblem(sys, 0.0)
initsol = solve(initprob)
u0 = initsol[unknowns(sys)]
prob = ODEProblem(sys, u0, (0, 0)) # Warning: Initialization system is overdetermined. 4 equations for 0 unknowns.

Note: my system is not overdetermined! I can see this because for one, the InitializationProblem build and solves just fine. Also I can simply remove u0 without incident...

julia> prob = ODEProblem(sys, [], (0, 0))
ODEProblem with uType Vector{Float64} and tType Int64. In-place: true
timespan: (0, 0)
u0: 6-element Vector{Float64}:
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
@bradcarman bradcarman added the bug Something isn't working label Dec 11, 2024
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
bug Something isn't working
Projects
None yet
Development

No branches or pull requests

1 participant