Merge pull request #2126 from SciML/myb/domain

Domain feature unit tests and fixes

Author: YingboMa

src/systems/connectors.jl

@@ -262,6 +262,7 @@ function generate_connection_set(sys::AbstractSystem, find = nothing, replace =
     domain_free_connectionsets = filter(connectionsets) do cset
         !any(s -> is_domain_connector(s.sys.sys), cset.set)
     end
+
     _, domainset = merge(connectionsets, true)
     sys, (merge(domain_free_connectionsets), domainset)
 end
@@ -395,7 +396,7 @@ function Graphs.outneighbors(g::SystemDomainGraph, n::Int)
     for s in g.csets[i].set
         sys = s.sys.sys
         is_domain_connector(s.sys.sys) && continue
-        ts = TearingState(s.sys.namespace)
+        ts = TearingState(expand_connections(s.sys.namespace))
         graph = ts.structure.graph
         mm = linear_subsys_adjmat!(ts)
         lineqs = BitSet(mm.nzrows)

test/stream_connectors.jl

@@ -3,10 +3,34 @@ using ModelingToolkit
 @variables t
 
 @connector function TwoPhaseFluidPort(; name, P = 0.0, m_flow = 0.0, h_outflow = 0.0)
-    vars = @variables h_outflow(t)=h_outflow [connect = Stream] m_flow(t)=m_flow [
-        connect = Flow,
-    ] P(t)=P
-    ODESystem(Equation[], t, vars, []; name = name)
+    pars = @parameters begin
+        rho
+        bulk
+        viscosity
+    end
+
+    vars = @variables begin
+        (h_outflow(t) = h_outflow), [connect = Stream]
+        (m_flow(t) = m_flow), [connect = Flow]
+        P(t) = P
+    end
+
+    ODESystem(Equation[], t, vars, pars; name = name)
+end
+
+@connector function TwoPhaseFluid(; name, R, B, V)
+    pars = @parameters begin
+        rho = R
+        bulk = B
+        viscosity = V
+    end
+
+    vars = @variables begin m_flow(t), [connect = Flow] end
+
+    # equations ---------------------------
+    eqs = Equation[m_flow ~ 0]
+
+    ODESystem(eqs, t, vars, pars; name)
 end
 
 function MassFlowSource_h(; name,
@@ -90,6 +114,7 @@ function N1M1(; name,
     sys = compose(sys, subs)
 end
 
+@named fluid = TwoPhaseFluid(; R = 876, B = 1.2e9, V = 0.034)
 @named n1m1 = N1M1()
 @named pipe = AdiabaticStraightPipe()
 @named sink = MassFlowSource_h(m_flow_in = -0.01, h_in = 400e3)
@@ -98,7 +123,13 @@ eqns = [connect(n1m1.port_a, pipe.port_a)
         connect(pipe.port_b, sink.port)]
 
 @named sys = ODESystem(eqns, t)
-@named n1m1Test = compose(sys, n1m1, pipe, sink)
+
+eqns = [connect(fluid, n1m1.port_a)
+        connect(n1m1.port_a, pipe.port_a)
+        connect(pipe.port_b, sink.port)]
+
+@named n1m1Test = ODESystem(eqns, t, [], []; systems = [fluid, n1m1, pipe, sink])
+
 @test_nowarn structural_simplify(n1m1Test)
 @unpack source, port_a = n1m1
 @test sort(equations(expand_connections(n1m1)), by = string) == [0 ~ port_a.m_flow
@@ -125,6 +156,7 @@ eqns = [connect(n1m1.port_a, pipe.port_a)
        0 ~ n1m1.port_a.m_flow + pipe.port_a.m_flow
        0 ~ n1m1.source.port1.m_flow - n1m1.port_a.m_flow
        0 ~ pipe.port_b.m_flow + sink.port.m_flow
+       fluid.m_flow ~ 0
        n1m1.port_a.P ~ pipe.port_a.P
        n1m1.source.port1.P ~ n1m1.port_a.P
        n1m1.source.port1.P ~ n1m1.source.P
@@ -254,3 +286,194 @@ sys = expand_connections(compose(simple, [vp1, vp2, vp3]))
 end
 
 @test_nowarn @named a = VectorHeatPort()
+
+# --------------------------------------------------
+# Test the new Domain feature
+
+sys_ = expand_connections(n1m1Test)
+sys_defs = ModelingToolkit.defaults(sys_)
+csys = complete(n1m1Test)
+@test Symbol(sys_defs[csys.pipe.port_a.rho]) == Symbol(csys.fluid.rho)
+@test Symbol(sys_defs[csys.pipe.port_b.rho]) == Symbol(csys.fluid.rho)
+
+# Testing the domain feature with non-stream system...
+
+@connector function HydraulicPort(; P, name)
+    pars = @parameters begin
+        p_int = P
+        rho
+        bulk
+        viscosity
+    end
+
+    vars = @variables begin
+        p(t) = p_int
+        dm(t), [connect = Flow]
+    end
+
+    # equations ---------------------------
+    eqs = Equation[]
+
+    ODESystem(eqs, t, vars, pars; name, defaults = [dm => 0])
+end
+
+@connector function Fluid(; name, R, B, V)
+    pars = @parameters begin
+        rho = R
+        bulk = B
+        viscosity = V
+    end
+
+    vars = @variables begin dm(t), [connect = Flow] end
+
+    # equations ---------------------------
+    eqs = [
+        dm ~ 0,
+    ]
+
+    ODESystem(eqs, t, vars, pars; name)
+end
+
+function StepSource(; P, name)
+    pars = @parameters begin p_int = P end
+
+    vars = []
+
+    # nodes -------------------------------
+    systems = @named begin H = HydraulicPort(; P = p_int) end
+
+    # equations ---------------------------
+    eqs = [
+        H.p ~ p_int * (t > 0.01),
+    ]
+
+    ODESystem(eqs, t, vars, pars; name, systems)
+end
+
+function Pipe(; P, R, name)
+    pars = @parameters begin
+        p_int = P
+        resistance = R
+    end
+
+    vars = []
+
+    # nodes -------------------------------
+    systems = @named begin
+        HA = HydraulicPort(; P = p_int)
+        HB = HydraulicPort(; P = p_int)
+    end
+
+    # equations ---------------------------
+    eqs = [HA.p - HB.p ~ HA.dm * resistance / HA.viscosity
+           0 ~ HA.dm + HB.dm]
+
+    ODESystem(eqs, t, vars, pars; name, systems)
+end
+
+function StaticVolume(; P, V, name)
+    D = Differential(t)
+
+    pars = @parameters begin
+        p_int = P
+        vol = V
+    end
+
+    vars = @variables begin
+        p(t) = p_int
+        vrho(t)
+        drho(t) = 0
+    end
+
+    # nodes -------------------------------
+    systems = @named begin H = HydraulicPort(; P = p_int) end
+
+    # fluid props ------------------------
+    rho_0 = H.rho
+
+    # equations ---------------------------
+    eqs = [D(vrho) ~ drho
+           vrho ~ rho_0 * (1 + p / H.bulk)
+           H.p ~ p
+           H.dm ~ drho * V]
+
+    ODESystem(eqs, t, vars, pars; name, systems,
+              defaults = [vrho => rho_0 * (1 + p_int / H.bulk)])
+end
+
+function TwoFluidSystem(; name)
+    pars = []
+    vars = []
+
+    # nodes -------------------------------
+    systems = @named begin
+        fluid_a = Fluid(; R = 876, B = 1.2e9, V = 0.034)
+        source_a = StepSource(; P = 10e5)
+        pipe_a = Pipe(; P = 0, R = 1e6)
+        volume_a = StaticVolume(; P = 0, V = 0.1)
+
+        fluid_b = Fluid(; R = 1000, B = 2.5e9, V = 0.00034)
+        source_b = StepSource(; P = 10e5)
+        pipe_b = Pipe(; P = 0, R = 1e6)
+        volume_b = StaticVolume(; P = 0, V = 0.1)
+    end
+
+    # equations ---------------------------
+    eqs = [connect(fluid_a, source_a.H)
+           connect(source_a.H, pipe_a.HA)
+           connect(pipe_a.HB, volume_a.H)
+           connect(fluid_b, source_b.H)
+           connect(source_b.H, pipe_b.HA)
+           connect(pipe_b.HB, volume_b.H)]
+
+    ODESystem(eqs, t, vars, pars; name, systems)
+end
+
+@named two_fluid_system = TwoFluidSystem()
+sys = expand_connections(two_fluid_system)
+
+sys_defs = ModelingToolkit.defaults(sys)
+csys = complete(two_fluid_system)
+
+@test Symbol(sys_defs[csys.volume_a.H.rho]) == Symbol(csys.fluid_a.rho)
+@test Symbol(sys_defs[csys.volume_b.H.rho]) == Symbol(csys.fluid_b.rho)
+
+@test_nowarn structural_simplify(two_fluid_system)
+
+function OneFluidSystem(; name)
+    pars = []
+    vars = []
+
+    # nodes -------------------------------
+    systems = @named begin
+        fluid = Fluid(; R = 876, B = 1.2e9, V = 0.034)
+
+        source_a = StepSource(; P = 10e5)
+        pipe_a = Pipe(; P = 0, R = 1e6)
+        volume_a = StaticVolume(; P = 0, V = 0.1)
+
+        source_b = StepSource(; P = 20e5)
+        pipe_b = Pipe(; P = 0, R = 1e6)
+        volume_b = StaticVolume(; P = 0, V = 0.1)
+    end
+
+    # equations ---------------------------
+    eqs = [connect(fluid, source_a.H, source_b.H)
+           connect(source_a.H, pipe_a.HA)
+           connect(pipe_a.HB, volume_a.H)
+           connect(source_b.H, pipe_b.HA)
+           connect(pipe_b.HB, volume_b.H)]
+
+    ODESystem(eqs, t, vars, pars; name, systems)
+end
+
+@named one_fluid_system = OneFluidSystem()
+sys = expand_connections(one_fluid_system)
+
+sys_defs = ModelingToolkit.defaults(sys)
+csys = complete(one_fluid_system)
+
+@test Symbol(sys_defs[csys.volume_a.H.rho]) == Symbol(csys.fluid.rho)
+@test Symbol(sys_defs[csys.volume_b.H.rho]) == Symbol(csys.fluid.rho)
+
+@test_nowarn structural_simplify(one_fluid_system)