Skip to content

Commit 29c2b54

Browse files
ChrisRackauckasChrisRackauckas
committed
Complete hierarchical
1 parent 452f8fe commit 29c2b54

File tree

4 files changed

+167
-3
lines changed

4 files changed

+167
-3
lines changed

‎src/systems/abstractsystem.jl

Lines changed: 7 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -889,6 +889,12 @@ function namespace_equations(sys::AbstractSystem, ivs = independent_variables(sy
889889
map(eq -> namespace_equation(eq, sys; ivs), eqs)
890890
end
891891

892+
function namespace_initialization_equations(sys::AbstractSystem, ivs = independent_variables(sys))
893+
eqs = initialization_equations(sys)
894+
isempty(eqs) && return Equation[]
895+
map(eq -> namespace_equation(eq, sys; ivs), eqs)
896+
end
897+
892898
function namespace_equation(eq::Equation,
893899
sys,
894900
n = nameof(sys);
@@ -1095,7 +1101,7 @@ function initialization_equations(sys::AbstractSystem)
10951101
else
10961102
eqs = Equation[eqs;
10971103
reduce(vcat,
1098-
namespace_equations.(get_systems(sys));
1104+
namespace_initialization_equations.(get_systems(sys));
10991105
init = Equation[])]
11001106
return eqs
11011107
end

‎src/systems/nonlinear/initializesystem.jl

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -23,6 +23,7 @@ function generate_initializesystem(sys::ODESystem;
2323
diffmap = Dict(getfield.(eqs_diff, :lhs) .=> getfield.(eqs_diff, :rhs))
2424
observed_diffmap = Dict(Differential(get_iv(sys)).(getfield.((observed(sys)), :lhs)) .=>
2525
Differential(get_iv(sys)).(getfield.((observed(sys)), :rhs)))
26+
full_diffmap = merge(diffmap, observed_diffmap)
2627

2728
full_states = unique([sts; getfield.((observed(sys)), :lhs)])
2829
set_full_states = Set(full_states)
@@ -39,8 +40,7 @@ function generate_initializesystem(sys::ODESystem;
3940
filtered_u0 = Pair[]
4041
for x in u0map
4142
y = get(schedule.dummy_sub, x[1], x[1])
42-
y = ModelingToolkit.fixpoint_sub(y, observed_diffmap)
43-
y = get(diffmap, y, y)
43+
y = ModelingToolkit.fixpoint_sub(y, full_diffmap)
4444

4545
if y isa Symbolics.Arr
4646
_y = collect(y)

‎test/hierarchical_initialization_eqs.jl

Lines changed: 157 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,157 @@
1+
using ModelingToolkit, OrdinaryDiffEq
2+
3+
t = only(@variables(t))
4+
D = Differential(t)
5+
"""
6+
A simple linear resistor model
7+
8+
![Resistor](https://encrypted-tbn0.gstatic.com/images?q=tbn:ANd9GcTpJkiEyqh-BRx27pvVH0GLZ4MP_D1oriBwJhnZdgIq7m17z9VKUWaW9MeNQAz1rTML2ho&usqp=CAU)
9+
"""
10+
@component function Resistor(; name, R=1.0)
11+
systems = @named begin
12+
p = Pin()
13+
n = Pin()
14+
end
15+
vars = @variables begin
16+
v(t), [guess=0.0]
17+
i(t), [guess=0.0]
18+
end
19+
params = @parameters begin
20+
R = R, [description = "Resistance of this Resistor"]
21+
end
22+
eqs = [
23+
v ~ p.v - n.v
24+
i ~ p.i
25+
p.i + n.i ~ 0
26+
# Ohm's Law
27+
v ~ i * R
28+
]
29+
return ODESystem(eqs, t, vars, params; systems, name)
30+
end
31+
@connector Pin begin
32+
v(t)
33+
i(t), [connect = Flow]
34+
end
35+
@component function ConstantVoltage(; name, V=1.0)
36+
systems = @named begin
37+
p = Pin()
38+
n = Pin()
39+
end
40+
vars = @variables begin
41+
v(t), [guess=0.0]
42+
i(t), [guess=0.0]
43+
end
44+
params = @parameters begin
45+
V = 10
46+
end
47+
eqs = [
48+
v ~ p.v - n.v
49+
i ~ p.i
50+
p.i + n.i ~ 0
51+
v ~ V
52+
]
53+
return ODESystem(eqs, t, vars, params; systems, name)
54+
end
55+
56+
@component function Capacitor(; name, C=1.0)
57+
systems = @named begin
58+
p = Pin()
59+
n = Pin()
60+
end
61+
vars = @variables begin
62+
v(t), [guess=0.0]
63+
i(t), [guess=0.0]
64+
end
65+
params = @parameters begin
66+
C = C
67+
end
68+
initialization_eqs = [
69+
v ~ 0
70+
]
71+
eqs = [
72+
v ~ p.v - n.v
73+
i ~ p.i
74+
p.i + n.i ~ 0
75+
C * D(v) ~ i
76+
]
77+
return ODESystem(eqs, t, vars, params; systems, name, initialization_eqs)
78+
end
79+
80+
@component function Ground(; name)
81+
systems = @named begin
82+
g = Pin()
83+
end
84+
eqs = [
85+
g.v ~ 0
86+
]
87+
return ODESystem(eqs, t, [], []; systems, name)
88+
end
89+
90+
@component function Inductor(; name, L=1.0)
91+
systems = @named begin
92+
p = Pin()
93+
n = Pin()
94+
end
95+
vars = @variables begin
96+
v(t), [guess = 0.0]
97+
i(t), [guess = 0.0]
98+
end
99+
params = @parameters begin
100+
(L = L)
101+
end
102+
eqs = [
103+
v ~ p.v - n.v
104+
i ~ p.i
105+
p.i + n.i ~ 0
106+
L * D(i) ~ v
107+
]
108+
return ODESystem(eqs, t, vars, params; systems, name)
109+
end
110+
111+
"""
112+
This is an RLC model. This should support markdown. That includes
113+
HTML as well.
114+
"""
115+
@component function RLCModel(; name)
116+
systems = @named begin
117+
resistor = Resistor(R=100)
118+
capacitor = Capacitor(C=0.001)
119+
inductor = Inductor(L=1)
120+
source = ConstantVoltage(V=30)
121+
ground = Ground()
122+
end
123+
initialization_eqs = [
124+
inductor.i ~ 0
125+
]
126+
eqs = [
127+
connect(source.p, inductor.n)
128+
connect(inductor.p, resistor.p, capacitor.p)
129+
connect(resistor.n, ground.g, capacitor.n, source.n)
130+
]
131+
return ODESystem(eqs, t, [], []; systems, name, initialization_eqs)
132+
end
133+
"""Run model RLCModel from 0 to 10"""
134+
function simple()
135+
@mtkbuild model = RLCModel()
136+
u0 = []
137+
prob = ODEProblem(model, u0, (0.0, 10.0))
138+
sol = solve(prob)
139+
end
140+
@test SciMLBase.successful_retcode(simple())
141+
142+
@named model = RLCModel()
143+
@test length(ModelingToolkit.get_initialization_eqs(model)) == 1
144+
syslist = ModelingToolkit.get_systems(model)
145+
@test length(ModelingToolkit.get_initialization_eqs(syslist[1])) == 0
146+
@test length(ModelingToolkit.get_initialization_eqs(syslist[2])) == 1
147+
@test length(ModelingToolkit.get_initialization_eqs(syslist[3])) == 0
148+
@test length(ModelingToolkit.get_initialization_eqs(syslist[4])) == 0
149+
@test length(ModelingToolkit.get_initialization_eqs(syslist[5])) == 0
150+
@test length(ModelingToolkit.initialization_equations(model)) == 2
151+
152+
u0 = []
153+
prob = ODEProblem(structural_simplify(model), u0, (0.0, 10.0))
154+
sol = solve(prob, Rodas5P())
155+
@test length(sol[end]) == 2
156+
@test length(equations(prob.f.initializeprob.f.sys)) == 0
157+
@test length(unknowns(prob.f.initializeprob.f.sys)) == 0

‎test/runtests.jl

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -38,6 +38,7 @@ end
3838
@safetestset "DDESystem Test" include("dde.jl")
3939
@safetestset "NonlinearSystem Test" include("nonlinearsystem.jl")
4040
@safetestset "InitializationSystem Test" include("initializationsystem.jl")
41+
@safetestset "Hierarchical Initialization Equations" include("hierarchical_initialization_eqs.jl")
4142
@safetestset "PDE Construction Test" include("pde.jl")
4243
@safetestset "JumpSystem Test" include("jumpsystem.jl")
4344
@safetestset "Constraints Test" include("constraints.jl")

0 commit comments

Comments
 (0)