-
-
Notifications
You must be signed in to change notification settings - Fork 234
Open
Labels
bugSomething isn't workingSomething isn't working
Description
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.1Metadata
Metadata
Assignees
Labels
bugSomething isn't workingSomething isn't working