diff --git a/Project.toml b/Project.toml index 2f48cc68fe..1de504cbae 100644 --- a/Project.toml +++ b/Project.toml @@ -123,7 +123,7 @@ Libdl = "1" LinearAlgebra = "1" Logging = "1" MLStyle = "0.4.17" -ModelingToolkitStandardLibrary = "2.19" +ModelingToolkitStandardLibrary = "2.20" Moshi = "0.3" NaNMath = "0.3, 1" NonlinearSolve = "4.3" diff --git a/examples/electrical_components.jl b/examples/electrical_components.jl deleted file mode 100644 index 1f4f151c21..0000000000 --- a/examples/electrical_components.jl +++ /dev/null @@ -1,95 +0,0 @@ -using Test -using ModelingToolkit, OrdinaryDiffEq -using ModelingToolkit: t_nounits as t, D_nounits as D - -@connector function Pin(; name) - sts = @variables v(t) [guess = 1.0] i(t) [guess = 1.0, connect = Flow] - ODESystem(Equation[], t, sts, []; name = name) -end - -@component function Ground(; name) - @named g = Pin() - eqs = [g.v ~ 0] - compose(ODESystem(eqs, t, [], []; name = name), g) -end - -@component function OnePort(; name) - @named p = Pin() - @named n = Pin() - sts = @variables v(t) [guess = 1.0] i(t) [guess = 1.0] - eqs = [v ~ p.v - n.v - 0 ~ p.i + n.i - i ~ p.i] - compose(ODESystem(eqs, t, sts, []; name = name), p, n) -end - -@component function Resistor(; name, R = 1.0) - @named oneport = OnePort() - @unpack v, i = oneport - ps = @parameters R = R - eqs = [ - v ~ i * R - ] - extend(ODESystem(eqs, t, [], ps; name = name), oneport) -end - -@component function Capacitor(; name, C = 1.0) - @named oneport = OnePort() - @unpack v, i = oneport - ps = @parameters C = C - eqs = [ - D(v) ~ i / C - ] - extend(ODESystem(eqs, t, [], ps; name = name), oneport) -end - -@component function ConstantVoltage(; name, V = 1.0) - @named oneport = OnePort() - @unpack v = oneport - ps = @parameters V = V - eqs = [ - V ~ v - ] - extend(ODESystem(eqs, t, [], ps; name = name), oneport) -end - -@component function Inductor(; name, L = 1.0) - @named oneport = OnePort() - @unpack v, i = oneport - ps = @parameters L = L - eqs = [ - D(i) ~ v / L - ] - extend(ODESystem(eqs, t, [], ps; name = name), oneport) -end - -@connector function HeatPort(; name) - @variables T(t) [guess = 293.15] Q_flow(t) [guess = 0.0, connect = Flow] - ODESystem(Equation[], t, [T, Q_flow], [], name = name) -end - -@component function HeatingResistor(; name, R = 1.0, TAmbient = 293.15, alpha = 1.0) - @named p = Pin() - @named n = Pin() - @named h = HeatPort() - @variables v(t) RTherm(t) - @parameters R=R TAmbient=TAmbient alpha=alpha - eqs = [RTherm ~ R * (1 + alpha * (h.T - TAmbient)) - v ~ p.i * RTherm - h.Q_flow ~ -v * p.i # -LossPower - v ~ p.v - n.v - 0 ~ p.i + n.i] - compose(ODESystem(eqs, t, [v, RTherm], [R, TAmbient, alpha], - name = name), p, n, h) -end - -@component function HeatCapacitor(; name, rho = 8050, V = 1, cp = 460, TAmbient = 293.15) - @parameters rho=rho V=V cp=cp - C = rho * V * cp - @named h = HeatPort() - eqs = [ - D(h.T) ~ h.Q_flow / C - ] - compose(ODESystem(eqs, t, [], [rho, V, cp], - name = name), h) -end diff --git a/examples/rc_model.jl b/examples/rc_model.jl deleted file mode 100644 index 158436d980..0000000000 --- a/examples/rc_model.jl +++ /dev/null @@ -1,17 +0,0 @@ -include("electrical_components.jl") - -R = 1.0 -C = 1.0 -V = 1.0 -@named resistor = Resistor(R = R) -@named capacitor = Capacitor(C = C) -@named source = ConstantVoltage(V = V) -@named ground = Ground() - -rc_eqs = [connect(source.p, resistor.p) - connect(resistor.n, capacitor.p) - connect(capacitor.n, source.n) - connect(capacitor.n, ground.g)] - -@named rc_model = ODESystem(rc_eqs, t) -rc_model = compose(rc_model, [resistor, capacitor, source, ground]) diff --git a/examples/serial_inductor.jl b/examples/serial_inductor.jl deleted file mode 100644 index 302df32c17..0000000000 --- a/examples/serial_inductor.jl +++ /dev/null @@ -1,33 +0,0 @@ -include("electrical_components.jl") - -@named source = ConstantVoltage(V = 10.0) -@named resistor = Resistor(R = 1.0) -@named inductor1 = Inductor(L = 1.0e-2) -@named inductor2 = Inductor(L = 2.0e-2) -@named ground = Ground() - -eqs = [connect(source.p, resistor.p) - connect(resistor.n, inductor1.p) - connect(inductor1.n, inductor2.p) - connect(source.n, inductor2.n) - connect(inductor2.n, ground.g)] - -@named ll_model = ODESystem(eqs, t) -ll_model = compose(ll_model, [source, resistor, inductor1, inductor2, ground]) - -@named source = ConstantVoltage(V = 10.0) -@named resistor1 = Resistor(R = 1.0) -@named resistor2 = Resistor(R = 1.0) -@named inductor1 = Inductor(L = 1.0e-2) -@named inductor2 = Inductor(L = 2.0e-2) -@named ground = Ground() - -eqs = [connect(source.p, inductor1.p) - connect(inductor1.n, resistor1.p) - connect(inductor1.n, resistor2.p) - connect(resistor1.n, resistor2.n) - connect(resistor2.n, inductor2.p) - connect(source.n, inductor2.n) - connect(inductor2.n, ground.g)] -@named ll2_model = ODESystem(eqs, t) -ll2_model = compose(ll2_model, [source, resistor1, resistor2, inductor1, inductor2, ground]) diff --git a/test/common/rc_model.jl b/test/common/rc_model.jl new file mode 100644 index 0000000000..8eec8d048f --- /dev/null +++ b/test/common/rc_model.jl @@ -0,0 +1,26 @@ +import ModelingToolkitStandardLibrary.Electrical as El +import ModelingToolkitStandardLibrary.Blocks as Bl + +@mtkmodel RCModel begin + @parameters begin + R = 1.0 + C = 1.0 + V = 1.0 + end + @components begin + resistor = El.Resistor(R = R) + capacitor = El.Capacitor(C = C) + shape = Bl.Constant(k = V) + source = El.Voltage() + ground = El.Ground() + end + @equations begin + connect(shape.output, source.V) + connect(source.p, resistor.p) + connect(resistor.n, capacitor.p) + connect(capacitor.n, source.n) + connect(capacitor.n, ground.g) + end +end + +@named rc_model = RCModel() diff --git a/test/common/serial_inductor.jl b/test/common/serial_inductor.jl new file mode 100644 index 0000000000..1e930cee56 --- /dev/null +++ b/test/common/serial_inductor.jl @@ -0,0 +1,47 @@ +import ModelingToolkitStandardLibrary.Electrical as El +import ModelingToolkitStandardLibrary.Blocks as Bl + +@mtkmodel LLModel begin + @components begin + shape = Bl.Constant(k = 10.0) + source = El.Voltage() + resistor = El.Resistor(R = 1.0) + inductor1 = El.Inductor(L = 1.0e-2) + inductor2 = El.Inductor(L = 2.0e-2) + ground = El.Ground() + end + @equations begin + connect(shape.output, source.V) + connect(source.p, resistor.p) + connect(resistor.n, inductor1.p) + connect(inductor1.n, inductor2.p) + connect(source.n, inductor2.n) + connect(inductor2.n, ground.g) + end +end + +@named ll_model = LLModel() + +@mtkmodel LL2Model begin + @components begin + shape = Bl.Constant(k = 10.0) + source = El.Voltage() + resistor1 = El.Resistor(R = 1.0) + resistor2 = El.Resistor(R = 1.0) + inductor1 = El.Inductor(L = 1.0e-2) + inductor2 = El.Inductor(L = 2.0e-2) + ground = El.Ground() + end + @equations begin + connect(shape.output, source.V) + connect(source.p, inductor1.p) + connect(inductor1.n, resistor1.p) + connect(inductor1.n, resistor2.p) + connect(resistor1.n, resistor2.n) + connect(resistor2.n, inductor2.p) + connect(source.n, inductor2.n) + connect(inductor2.n, ground.g) + end +end + +@named ll2_model = LL2Model() diff --git a/test/components.jl b/test/components.jl index 8ac40f6fbb..43c0d1acf8 100644 --- a/test/components.jl +++ b/test/components.jl @@ -4,211 +4,159 @@ using ModelingToolkit: get_component_type using ModelingToolkit.BipartiteGraphs using ModelingToolkit.StructuralTransformations using ModelingToolkit: t_nounits as t, D_nounits as D -include("../examples/rc_model.jl") - -function check_contract(sys) - state = ModelingToolkit.get_tearing_state(sys) - graph = state.structure.graph - fullvars = state.fullvars - sys = tearing_substitution(sys) - - eqs = equations(sys) - for (i, eq) in enumerate(eqs) - actual = union(ModelingToolkit.vars(eq.lhs), ModelingToolkit.vars(eq.rhs)) - actual = filter(!ModelingToolkit.isparameter, collect(actual)) - current = Set(fullvars[𝑠neighbors(graph, i)]) - @test isempty(setdiff(actual, current)) +using ModelingToolkitStandardLibrary.Electrical +using ModelingToolkitStandardLibrary.Blocks +using LinearAlgebra +using ModelingToolkitStandardLibrary.Thermal +include("common/rc_model.jl") + +@testset "Basics" begin + @unpack resistor, capacitor, source = rc_model + function check_contract(sys) + state = ModelingToolkit.get_tearing_state(sys) + graph = state.structure.graph + fullvars = state.fullvars + sys = tearing_substitution(sys) + + eqs = equations(sys) + for (i, eq) in enumerate(eqs) + actual = union(ModelingToolkit.vars(eq.lhs), ModelingToolkit.vars(eq.rhs)) + actual = filter(!ModelingToolkit.isparameter, collect(actual)) + current = Set(fullvars[𝑠neighbors(graph, i)]) + @test isempty(setdiff(actual, current)) + end end -end -function check_rc_sol(sol) - rpi = sol[rc_model.resistor.p.i] - rpifun = sol.prob.f.observed(rc_model.resistor.p.i) - @test rpifun.(sol.u, (sol.prob.p,), sol.t) == rpi - @test any(!isequal(rpi[1]), rpi) # test that we don't have a constant system - @test sol[rc_model.resistor.p.i] == sol[resistor.p.i] == sol[capacitor.p.i] - @test sol[rc_model.resistor.n.i] == sol[resistor.n.i] == -sol[capacitor.p.i] - @test sol[rc_model.capacitor.n.i] == sol[capacitor.n.i] == -sol[capacitor.p.i] - @test iszero(sol[rc_model.ground.g.i]) - @test iszero(sol[rc_model.ground.g.v]) - @test sol[rc_model.resistor.v] == sol[resistor.v] == - sol[source.p.v] - sol[capacitor.p.v] -end - -@named pin = Pin() -@test get_component_type(pin).name == :Pin -@test get_component_type(rc_model.resistor).name == :Resistor - -completed_rc_model = complete(rc_model) -@test isequal(completed_rc_model.resistor.n.i, resistor.n.i) -@test ModelingToolkit.n_expanded_connection_equations(capacitor) == 2 -@test length(equations(structural_simplify(rc_model, allow_parameter = false))) == 2 -sys = structural_simplify(rc_model) -@test_throws ModelingToolkit.RepeatedStructuralSimplificationError structural_simplify(sys) -@test length(equations(sys)) == 1 -check_contract(sys) -@test !isempty(ModelingToolkit.defaults(sys)) -u0 = [capacitor.v => 0.0] -prob = ODEProblem(sys, u0, (0, 10.0)) -sol = solve(prob, Rodas4()) -check_rc_sol(sol) - -# https://discourse.julialang.org/t/using-optimization-parameters-in-modelingtoolkit/82099 -let - @parameters param_r1 param_c1 - @named resistor = Resistor(R = param_r1) - @named capacitor = Capacitor(C = param_c1) - @named source = ConstantVoltage(V = 1.0) - @named ground = Ground() + function check_rc_sol(sol) + rpi = sol[rc_model.resistor.p.i] + rpifun = sol.prob.f.observed(rc_model.resistor.p.i) + @test rpifun.(sol.u, (sol.prob.p,), sol.t) == rpi + @test any(!isequal(rpi[1]), rpi) # test that we don't have a constant system + @test sol[rc_model.resistor.p.i] == sol[resistor.p.i] == sol[capacitor.p.i] + @test sol[rc_model.resistor.n.i] == sol[resistor.n.i] == -sol[capacitor.p.i] + @test sol[rc_model.capacitor.n.i] == sol[capacitor.n.i] == -sol[capacitor.p.i] + @test iszero(sol[rc_model.ground.g.i]) + @test iszero(sol[rc_model.ground.g.v]) + @test sol[rc_model.resistor.v] == sol[resistor.v] == + sol[source.p.v] - sol[capacitor.p.v] + end - rc_eqs = [connect(source.p, resistor.p) - connect(resistor.n, capacitor.p) - connect(capacitor.n, source.n) - connect(capacitor.n, ground.g)] + @named pin = Pin() + @test get_component_type(pin).name == :Pin + @test get_component_type(rc_model.resistor).name == :Resistor - @named _rc_model = ODESystem(rc_eqs, t) - @named rc_model = compose(_rc_model, - [resistor, capacitor, source, ground]) + completed_rc_model = complete(rc_model) + @test isequal(completed_rc_model.resistor.n.i, resistor.n.i) + @test ModelingToolkit.n_expanded_connection_equations(capacitor) == 2 + @test length(equations(structural_simplify(rc_model, allow_parameter = false))) == 2 sys = structural_simplify(rc_model) - u0 = [ - capacitor.v => 0.0 - ] - - params = [param_r1 => 1.0, param_c1 => 1.0] - tspan = (0.0, 10.0) - - prob = ODEProblem(sys, u0, tspan, params) - @test solve(prob, Tsit5()).retcode == ReturnCode.Success + @test_throws ModelingToolkit.RepeatedStructuralSimplificationError structural_simplify(sys) + @test length(equations(sys)) == 1 + check_contract(sys) + @test !isempty(ModelingToolkit.defaults(sys)) + u0 = [capacitor.v => 0.0] + prob = ODEProblem(sys, u0, (0, 10.0)) + sol = solve(prob, Rodas4()) + check_rc_sol(sol) end -let - # 1478 - @named resistor2 = Resistor(R = R) - rc_eqs2 = [connect(source.p, resistor.p) - connect(resistor.n, resistor2.p) - connect(resistor2.n, capacitor.p) - connect(capacitor.n, source.n) - connect(capacitor.n, ground.g)] - - @named _rc_model2 = ODESystem(rc_eqs2, t) - @named rc_model2 = compose(_rc_model2, - [resistor, resistor2, capacitor, source, ground]) - sys2 = structural_simplify(rc_model2) - prob2 = ODEProblem(sys2, [source.p.i => 0.0], (0, 10.0), guesses = u0) - sol2 = solve(prob2, Rosenbrock23()) - @test sol2[source.p.i] ≈ sol2[rc_model2.source.p.i] ≈ -sol2[capacitor.i] - - prob3 = ODEProblem(sys2, [], (0, 10.0), guesses = u0) - sol3 = solve(prob2, Rosenbrock23()) - @test sol3[unknowns(rc_model2), end] ≈ sol2[unknowns(rc_model2), end] -end +@testset "Outer/inner connections" begin + sys = structural_simplify(rc_model) -# Outer/inner connections -function rc_component(; name, R = 1, C = 1) - @parameters R=R C=C - @named p = Pin() - @named n = Pin() - @named resistor = Resistor(R = R) # test parent scope default of @named - @named capacitor = Capacitor(C = ParentScope(C)) - eqs = [connect(p, resistor.p); - connect(resistor.n, capacitor.p); - connect(capacitor.n, n)] - @named sys = ODESystem(eqs, t, [], [R, C]) - compose(sys, [p, n, resistor, capacitor]; name = name) -end + prob = ODEProblem(sys, [sys.capacitor.v => 0.0], (0.0, 10.0)) + sol = solve(prob, Rodas4()) + function rc_component(; name, R = 1, C = 1) + local sys + @parameters R=R C=C + @named p = Pin() + @named n = Pin() + @named resistor = Resistor(R = R) # test parent scope default of @named + @named capacitor = Capacitor(C = C) + eqs = [connect(p, resistor.p); + connect(resistor.n, capacitor.p); + connect(capacitor.n, n)] + @named sys = ODESystem(eqs, t, [], [R, C]) + compose(sys, [p, n, resistor, capacitor]; name = name) + end -@named ground = Ground() -@named source = ConstantVoltage(V = 1) -@named rc_comp = rc_component() -eqs = [connect(source.p, rc_comp.p) - connect(source.n, rc_comp.n) - connect(source.n, ground.g)] -@named sys′ = ODESystem(eqs, t) -@named sys_inner_outer = compose(sys′, [ground, source, rc_comp]) -@test_nowarn show(IOBuffer(), MIME"text/plain"(), sys_inner_outer) -expand_connections(sys_inner_outer, debug = true) -sys_inner_outer = structural_simplify(sys_inner_outer) -@test !isempty(ModelingToolkit.defaults(sys_inner_outer)) -u0 = [rc_comp.capacitor.v => 0.0] -prob = ODEProblem(sys_inner_outer, u0, (0, 10.0), sparse = true) -sol_inner_outer = solve(prob, Rodas4()) -@test sol[capacitor.v] ≈ sol_inner_outer[rc_comp.capacitor.v] - -u0 = [ - capacitor.v => 0.0 -] -prob = ODEProblem(sys, u0, (0, 10.0)) -sol = solve(prob, Tsit5()) - -@test sol[resistor.p.i] == sol[capacitor.p.i] -@test sol[resistor.n.i] == -sol[capacitor.p.i] -@test sol[capacitor.n.i] == -sol[capacitor.p.i] -@test iszero(sol[ground.g.i]) -@test iszero(sol[ground.g.v]) -@test sol[resistor.v] == sol[source.p.v] - sol[capacitor.p.v] + @named ground = Ground() + @named shape = Constant(k = 1) + @named source = Voltage() + @named rc_comp = rc_component() + eqs = [connect(shape.output, source.V) + connect(source.p, rc_comp.p) + connect(source.n, rc_comp.n) + connect(source.n, ground.g)] + @named sys′ = ODESystem(eqs, t) + @named sys_inner_outer = compose(sys′, [ground, shape, source, rc_comp]) + @test_nowarn show(IOBuffer(), MIME"text/plain"(), sys_inner_outer) + expand_connections(sys_inner_outer, debug = true) + sys_inner_outer = structural_simplify(sys_inner_outer) + @test !isempty(ModelingToolkit.defaults(sys_inner_outer)) + u0 = [rc_comp.capacitor.v => 0.0] + prob = ODEProblem(sys_inner_outer, u0, (0, 10.0), sparse = true) + sol_inner_outer = solve(prob, Rodas4()) + @test sol[sys.capacitor.v] ≈ sol_inner_outer[rc_comp.capacitor.v] + + prob = ODEProblem(sys, [sys.capacitor.v => 0.0], (0, 10.0)) + sol = solve(prob, Tsit5()) + + @test sol[sys.resistor.p.i] == sol[sys.capacitor.p.i] + @test sol[sys.resistor.n.i] == -sol[sys.capacitor.p.i] + @test sol[sys.capacitor.n.i] == -sol[sys.capacitor.p.i] + @test iszero(sol[sys.ground.g.i]) + @test iszero(sol[sys.ground.g.v]) + @test sol[sys.resistor.v] == sol[sys.source.p.v] - sol[sys.capacitor.p.v] +end #using Plots #plot(sol) -include("../examples/serial_inductor.jl") -sys = structural_simplify(ll_model) -@test length(equations(sys)) == 2 -u0 = unknowns(sys) .=> 0 -@test_nowarn ODEProblem( - sys, [], (0, 10.0), guesses = u0, warn_initialize_determined = false) -prob = DAEProblem(sys, D.(unknowns(sys)) .=> 0, [], (0, 0.5), guesses = u0) -sol = solve(prob, DFBDF()) -@test sol.retcode == SciMLBase.ReturnCode.Success - -sys2 = structural_simplify(ll2_model) -@test length(equations(sys2)) == 3 -u0 = unknowns(sys) .=> 0 -prob = ODEProblem(sys, u0, (0, 10.0)) -@test_nowarn sol = solve(prob, FBDF()) - -@variables x1(t) x2(t) x3(t) x4(t) -@named sys1_inner = ODESystem([D(x1) ~ x1], t) -@named sys1_partial = compose(ODESystem([D(x2) ~ x2], t; name = :foo), sys1_inner) -@named sys1 = extend(ODESystem([D(x3) ~ x3], t; name = :foo), sys1_partial) -@named sys2 = compose(ODESystem([D(x4) ~ x4], t; name = :foo), sys1) -@test_nowarn sys2.sys1.sys1_inner.x1 # test the correct nesting - -# compose tests -function record_fun(; name) - pars = @parameters a=10 b=100 - ODESystem(Equation[], t, [], pars; name) +include("common/serial_inductor.jl") +@testset "Serial inductor" begin + sys = structural_simplify(ll_model) + @test length(equations(sys)) == 2 + u0 = unknowns(sys) .=> 0 + @test_nowarn ODEProblem( + sys, [], (0, 10.0), guesses = u0, warn_initialize_determined = false) + prob = DAEProblem(sys, D.(unknowns(sys)) .=> 0, [], (0, 0.5), guesses = u0) + sol = solve(prob, DFBDF()) + @test sol.retcode == SciMLBase.ReturnCode.Success + + sys2 = structural_simplify(ll2_model) + @test length(equations(sys2)) == 3 + u0 = unknowns(sys) .=> 0 + prob = ODEProblem(sys, u0, (0, 10.0)) + @test_nowarn sol = solve(prob, FBDF()) end -function first_model(; name) - @named foo = record_fun() +@testset "Compose/extend" begin + @variables x1(t) x2(t) x3(t) x4(t) + @named sys1_inner = ODESystem([D(x1) ~ x1], t) + @named sys1_partial = compose(ODESystem([D(x2) ~ x2], t; name = :foo), sys1_inner) + @named sys1 = extend(ODESystem([D(x3) ~ x3], t; name = :foo), sys1_partial) + @named sys2 = compose(ODESystem([D(x4) ~ x4], t; name = :foo), sys1) + @test_nowarn sys2.sys1.sys1_inner.x1 # test the correct nesting + + # compose tests + function record_fun(; name) + pars = @parameters a=10 b=100 + ODESystem(Equation[], t, [], pars; name) + end + + function first_model(; name) + @named foo = record_fun() - defs = Dict() - defs[foo.a] = 3 - defs[foo.b] = 300 - pars = @parameters x=2 y=20 - compose(ODESystem(Equation[], t, [], pars; name, defaults = defs), foo) + defs = Dict() + defs[foo.a] = 3 + defs[foo.b] = 300 + pars = @parameters x=2 y=20 + compose(ODESystem(Equation[], t, [], pars; name, defaults = defs), foo) + end + @named goo = first_model() + @unpack foo = goo + @test ModelingToolkit.defaults(goo)[foo.a] == 3 + @test ModelingToolkit.defaults(goo)[foo.b] == 300 end -@named goo = first_model() -@unpack foo = goo -@test ModelingToolkit.defaults(goo)[foo.a] == 3 -@test ModelingToolkit.defaults(goo)[foo.b] == 300 - -#= -model Circuit - Ground ground; - Load load; - Resistor resistor; -equation - connect(load.p , ground.p); - connect(resistor.p, ground.p); -end Circuit; -model Load - extends TwoPin; - Resistor resistor; -equation - connect(p, resistor.p); - connect(resistor.n, n); -end Load; -=# function Load(; name) R = 1 @@ -236,62 +184,68 @@ end @test structural_simplify(foo) isa ModelingToolkit.AbstractSystem # BLT tests -using LinearAlgebra -function parallel_rc_model(i; name, source, ground, R, C) - resistor = HeatingResistor(name = Symbol(:resistor, i), R = R) - capacitor = Capacitor(name = Symbol(:capacitor, i), C = C) - heat_capacitor = HeatCapacitor(name = Symbol(:heat_capacitor, i)) - - rc_eqs = [connect(source.p, resistor.p) - connect(resistor.n, capacitor.p) - connect(capacitor.n, source.n, ground.g) - connect(resistor.h, heat_capacitor.h)] - - compose(ODESystem(rc_eqs, t, name = Symbol(name, i)), - [resistor, capacitor, source, ground, heat_capacitor]) +@testset "BLT ordering" begin + function parallel_rc_model(i; name, shape, source, ground, R, C) + resistor = Resistor(name = Symbol(:resistor, i), R = R, T_dep = true) + capacitor = Capacitor(name = Symbol(:capacitor, i), C = C) + heat_capacitor = HeatCapacitor(name = Symbol(:heat_capacitor, i)) + + rc_eqs = [connect(shape.output, source.V) + connect(source.p, resistor.p) + connect(resistor.n, capacitor.p) + connect(capacitor.n, source.n, ground.g) + connect(resistor.heat_port, heat_capacitor.port)] + + compose(ODESystem(rc_eqs, t, name = Symbol(name, i)), + [resistor, capacitor, source, ground, heat_capacitor]) + end + V = 2.0 + @named shape = Constant(k = V) + @named source = Voltage() + @named ground = Ground() + N = 50 + Rs = 10 .^ range(0, stop = -4, length = N) + Cs = 10 .^ range(-3, stop = 0, length = N) + rc_systems = map(1:N) do i + parallel_rc_model(i; name = :rc, source, ground, shape, R = Rs[i], C = Cs[i]) + end + @variables E(t) = 0.0 + eqs = [ + D(E) ~ sum(((i, sys),) -> getproperty(sys, Symbol(:resistor, i)).heat_port.Q_flow, + enumerate(rc_systems)) + ] + @named _big_rc = ODESystem(eqs, t, [E], []) + @named big_rc = compose(_big_rc, rc_systems) + ts = TearingState(expand_connections(big_rc)) + # this is block upper triangular, so `istriu` needs a little leeway + @test istriu(but_ordered_incidence(ts)[1], -2) end -V = 2.0 -@named source = ConstantVoltage(V = V) -@named ground = Ground() -N = 50 -Rs = 10 .^ range(0, stop = -4, length = N) -Cs = 10 .^ range(-3, stop = 0, length = N) -rc_systems = map(1:N) do i - parallel_rc_model(i; name = :rc, source = source, ground = ground, R = Rs[i], C = Cs[i]) -end; -@variables E(t) = 0.0 -eqs = [ - D(E) ~ sum(((i, sys),) -> getproperty(sys, Symbol(:resistor, i)).h.Q_flow, - enumerate(rc_systems)) -] -@named _big_rc = ODESystem(eqs, t, [E], []) -@named big_rc = compose(_big_rc, rc_systems) -ts = TearingState(expand_connections(big_rc)) -@test istriu(but_ordered_incidence(ts)[1]) # Test using constants inside subsystems -function FixedResistor(; name, R = 1.0) - @named oneport = OnePort() - @unpack v, i = oneport - @constants R = R - eqs = [ - v ~ i * R - ] - extend(ODESystem(eqs, t, [], []; name = name), oneport) +@testset "Constants inside subsystems" begin + function FixedResistor(; name, R = 1.0) + @named oneport = OnePort() + @unpack v, i = oneport + @constants R = R + eqs = [ + v ~ i * R + ] + extend(ODESystem(eqs, t, [], []; name = name), oneport) + end + capacitor = Capacitor(; name = :c1, C = 1.0) + resistor = FixedResistor(; name = :r1) + ground = Ground(; name = :ground) + rc_eqs = [connect(capacitor.n, resistor.p) + connect(resistor.n, capacitor.p) + connect(capacitor.n, ground.g)] + + @named _rc_model = ODESystem(rc_eqs, t) + @named rc_model = compose(_rc_model, + [resistor, capacitor, ground]) + sys = structural_simplify(rc_model) + prob = ODEProblem(sys, [sys.c1.v => 0.0], (0, 10.0)) + sol = solve(prob, Tsit5()) end -capacitor = Capacitor(; name = :c1) -resistor = FixedResistor(; name = :r1) -ground = Ground(; name = :ground) -rc_eqs = [connect(capacitor.n, resistor.p) - connect(resistor.n, capacitor.p) - connect(capacitor.n, ground.g)] - -@named _rc_model = ODESystem(rc_eqs, t) -@named rc_model = compose(_rc_model, - [resistor, capacitor, ground]) -sys = structural_simplify(rc_model) -prob = ODEProblem(sys, u0, (0, 10.0)) -sol = solve(prob, Tsit5()) @testset "docstrings (#1155)" begin """ diff --git a/test/error_handling.jl b/test/error_handling.jl index 59aa6b79a1..d5605e3cbc 100644 --- a/test/error_handling.jl +++ b/test/error_handling.jl @@ -1,8 +1,9 @@ using Test using ModelingToolkit +using ModelingToolkit: t_nounits as t, D_nounits as D import ModelingToolkit: ExtraVariablesSystemException, ExtraEquationsSystemException -include("../examples/electrical_components.jl") +using ModelingToolkitStandardLibrary.Electrical function UnderdefinedConstantVoltage(; name, V = 1.0) val = V diff --git a/test/funcaffect.jl b/test/funcaffect.jl index 3004044d61..05543c9161 100644 --- a/test/funcaffect.jl +++ b/test/funcaffect.jl @@ -1,4 +1,6 @@ using ModelingToolkit, Test, OrdinaryDiffEq +using ModelingToolkitStandardLibrary.Electrical +using ModelingToolkitStandardLibrary.Blocks using ModelingToolkit: t_nounits as t, D_nounits as D @constants h=1 zr=0 @@ -134,19 +136,21 @@ prob = ODEProblem(s2, [resistor.v => 10.0], (0, 2.01)) sol = solve(prob, Tsit5()) @test ctx[1] == 2 -include("../examples/rc_model.jl") +include("common/rc_model.jl") function affect5!(integ, u, p, ctx) @test integ.u[u.capacitor₊v] ≈ 0.3 integ.ps[p.C] *= 200 end -@named rc_model = ODESystem(rc_eqs, t, +@unpack capacitor = rc_model +@named event_sys = ODESystem(Equation[], t; continuous_events = [ [capacitor.v ~ 0.3] => (affect5!, [capacitor.v], [capacitor.C => :C], [capacitor.C], nothing) ]) -rc_model = compose(rc_model, [resistor, capacitor, source, ground]) +rc_model = extend(rc_model, event_sys) +# rc_model = compose(rc_model, [resistor, capacitor, source, ground]) sys = structural_simplify(rc_model) u0 = [capacitor.v => 0.0 @@ -177,15 +181,23 @@ function Capacitor2(; name, C = 1.0) oneport) end -@named capacitor2 = Capacitor2(C = C) +@named begin + capacitor2 = Capacitor2(C = 1.0) + resistor = Resistor(R = 1.0) + capacitor = Capacitor(C = 1.0) + shape = Constant(k = 1.0) + source = Voltage() + ground = Ground() +end -rc_eqs2 = [connect(source.p, resistor.p) +rc_eqs2 = [connect(shape.output, source.V) + connect(source.p, resistor.p) connect(resistor.n, capacitor2.p) connect(capacitor2.n, source.n) connect(capacitor2.n, ground.g)] @named rc_model2 = ODESystem(rc_eqs2, t) -rc_model2 = compose(rc_model2, [resistor, capacitor2, source, ground]) +rc_model2 = compose(rc_model2, [resistor, capacitor2, shape, source, ground]) sys2 = structural_simplify(rc_model2) u0 = [capacitor2.v => 0.0 diff --git a/test/print_tree.jl b/test/print_tree.jl index 7584a22c60..351e95f612 100644 --- a/test/print_tree.jl +++ b/test/print_tree.jl @@ -1,6 +1,6 @@ using ModelingToolkit, AbstractTrees, Test -include("../examples/rc_model.jl") +include("common/rc_model.jl") io = IOBuffer() print_tree(io, rc_model) @@ -12,9 +12,12 @@ str = """rc_model ├─ capacitor │ ├─ p │ └─ n + ├─ shape + │ └─ output ├─ source │ ├─ p - │ └─ n + │ ├─ n + │ └─ V └─ ground └─ g """ diff --git a/test/serialization.jl b/test/serialization.jl index feb3c7e4e7..ee380b6616 100644 --- a/test/serialization.jl +++ b/test/serialization.jl @@ -22,12 +22,14 @@ for prob in [ run(`$(Base.julia_cmd()) -e $(_cmd)`) end -include("../examples/rc_model.jl") +include("common/rc_model.jl") +@unpack capacitor = rc_model io = IOBuffer() -write(io, rc_model) +write(io, expand_connections(rc_model)) str = String(take!(io)) + sys = include_string(@__MODULE__, str) -@test sys == flatten(rc_model) # this actually kind of works, but the variables would have different identities. +@test sys == expand_connections(rc_model) # this actually kind of works, but the variables would have different identities. # check answer ss = structural_simplify(rc_model)