Inconsistency of Initialization System

[bradcarman]

Take the following model. Both sys1 and sys2 are equivalent, however sys1 does not give any warnings about the state of the initialization system (even though it should). As can be seen, sys2 simply defines the same initial equations but using u0 at the ODEProblem level, here we get the correct warning about initial conditions underdefinied.

using ModelingToolkit
using ModelingToolkit: D_nounits as D, t_nounits as t


function ClosedVolume(use_initial_eqs::Bool; name, p_int, x_int=0, direction)
    pars = @parameters begin
        rho_0 = 876
        beta = 1.2e9
        A =  0.1 
        L =  0.01 
        p_int = p_int
        x_int = x_int
        direction = direction
    end   
    vars = @variables begin
        p(t), [guess = p_int]
        x(t), [guess = x_int]
        rho(t), [guess= rho_0]
        m(t), [guess=rho*x*A]
    end
    eqs = [
        m ~ rho*(L+x*direction)*A 
        rho ~ rho_0*(1 + p/beta)
        D(m) ~ 0 
    ]

    initialization_eqs = if use_initial_eqs
        [
            p ~ p_int
            x ~ x_int
        ]
    else
        Equation[]
    end

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


function System(use_initial_eqs::Bool; name)
    pars = @parameters begin
        m =  100         
        c =  1000 
        force =  3e5 
        p_int = 10e5
        A =  0.1 
    end   
    vars = @variables begin
        p_1(t), [guess=p_int]
        p_2(t), [guess=p_int]
        x(t), [guess=0]
        dx(t), [guess=1]
        ddx(t), [guess=0]
    end
    systems = @named begin
        v1 = ClosedVolume(use_initial_eqs; p_int, direction=+1)
        v2 = ClosedVolume(use_initial_eqs; p_int, direction=-1)
    end
    eqs = [
        x ~ v1.x
        x ~ v2.x

        D(x) ~ dx
        D(dx) ~ ddx
    
        p_1 ~ v1.p
        p_2 ~ v2.p
    
        m*ddx ~ (p_1 - p_2)*A - c*dx + force*sin(2*pi*t/0.01)
    ]
    return ODESystem(eqs, t, vars, pars; name, systems)
end    

@mtkbuild sys1 = System(true)
prob1 = ODEProblem(sys1, [], (0,0))

#=
julia> initialization_equations(sys1)
4-element Vector{Equation}:
 v1₊p(t) ~ v1₊p_int
 v1₊x(t) ~ v1₊x_int
 v2₊p(t) ~ v2₊p_int
 v2₊x(t) ~ v2₊x_int
=#

# Build the equivalent system as sys1 with same initial equations given with u0...
@mtkbuild sys2 = System(false)

u0 = [
    sys2.v1.p => 10e5
    sys2.v1.x => 0.0
    sys2.v2.p => 10e5
    sys2.v2.x => 0.0
]

prob2 = ODEProblem(sys2, u0, (0,0)) # Warning: Initialization system is underdetermined. 0 equations for 1 unknowns.

Current versions...

  [961ee093] ModelingToolkit v9.60.0
  [1dea7af3] OrdinaryDiffEq v6.90.1

[bradcarman]

Note, for sys1, the fully_determined keyword does work as expected...

julia> prob1 = ODEProblem(sys1, [], (0,0); fully_determined = true)
ERROR: InvalidSystemException: The system is structurally singular! Here are the problematic variables: 
 dx(t)

[AayushSabharwal]

Yeah prob1 should print a singular warning

[AayushSabharwal]

Okay so the same consistency check is not valid after the rest of structural simplify which makes this annoying to handle