Merge branch 'master' into myb/sde

Author: YingboMa

.github/workflows/Downstream.yml

@@ -26,6 +26,9 @@ jobs:
           - {user: SciML, repo: StructuralIdentifiability.jl, group: All}
           - {user: SciML, repo: ModelingToolkitStandardLibrary.jl}
           - {user: SciML, repo: ModelOrderReduction.jl, group: All}
+          - {user: SciML, repo: MethodOfLines.jl, group: Interface}
+          - {user: SciML, repo: MethodOfLines.jl, group: 2D_Diffusion}
+          - {user: SciML, repo: MethodOfLines.jl, group: DAE}
           - {user: ai4energy, repo: Ai4EComponentLib.jl, group: Downstream}
     steps:
       - uses: actions/checkout@v2

src/ModelingToolkit.jl

@@ -163,8 +163,9 @@ end
 
 export AbstractTimeDependentSystem, AbstractTimeIndependentSystem,
        AbstractMultivariateSystem
-export ODESystem, ODEFunction, ODEFunctionExpr, ODEProblemExpr, convert_system
-export System
+
+export ODESystem, ODEFunction, ODEFunctionExpr, ODEProblemExpr, convert_system,
+       add_accumulations, System
 export DAEFunctionExpr, DAEProblemExpr
 export SDESystem, SDEFunction, SDEFunctionExpr, SDEProblemExpr
 export SystemStructure

src/systems/connectors.jl

@@ -378,8 +378,26 @@ function expand_instream(csets::AbstractVector{<:ConnectionSet}, sys::AbstractSy
         end
     end
 
+    if !isempty(instream_exprs)
+        # map from a namespaced stream variable to a ConnectionSet
+        expr_cset = Dict()
+        for cset in csets
+            crep = first(cset.set)
+            current = namespace == crep.sys.namespace
+            for v in cset.set
+                if (current || !v.isouter)
+                    expr_cset[namespaced_var(v)] = cset.set
+                end
+            end
+        end
+    end
+
     for ex in instream_exprs
-        cset, idx_in_set, sv = get_cset_sv(namespace, ex, csets)
+        ns_sv = only(arguments(ex))
+        full_name_sv = renamespace(namespace, ns_sv)
+        cset = get(expr_cset, full_name_sv, nothing)
+        cset === nothing && error("$ns_sv is not a variable inside stream connectors")
+        idx_in_set, sv = get_cset_sv(full_name_sv, cset)
 
         n_inners = n_outers = 0
         for (i, e) in enumerate(cset)
@@ -505,17 +523,10 @@ function get_current_var(namespace, cele, sv)
     states(renamespace(unnamespace(namespace, cele.sys.namespace), cele.sys.sys), sv)
 end
 
-function get_cset_sv(namespace, ex, csets)
-    ns_sv = only(arguments(ex))
-    full_name_sv = renamespace(namespace, ns_sv)
-
-    for c in csets
-        crep = first(c.set)
-        current = namespace == crep.sys.namespace
-        for (idx_in_set, v) in enumerate(c.set)
-            if (current || !v.isouter) && isequal(namespaced_var(v), full_name_sv)
-                return c.set, idx_in_set, v.v
-            end
+function get_cset_sv(full_name_sv, cset)
+    for (idx_in_set, v) in enumerate(cset)
+        if isequal(namespaced_var(v), full_name_sv)
+            return idx_in_set, v.v
         end
     end
     error("$ns_sv is not a variable inside stream connectors")

src/systems/diffeqs/odesystem.jl

@@ -436,3 +436,21 @@ function convert_system(::Type{<:ODESystem}, sys, t; name = nameof(sys))
     return ODESystem(neweqs, t, newsts, parameters(sys); defaults = defs, name = name,
                      checks = false)
 end
+
+"""
+$(SIGNATURES)
+
+Add accumulation variables for `vars`.
+"""
+function add_accumulations(sys::ODESystem, vars = states(sys))
+    eqs = get_eqs(sys)
+    accs = filter(x -> startswith(string(x), "accumulation_"), states(sys))
+    if !isempty(accs)
+        error("$accs variable names start with \"accumulation_\"")
+    end
+    avars = [rename(v, Symbol(:accumulation_, getname(v))) for v in vars]
+    D = Differential(get_iv(sys))
+    @set! sys.eqs = [eqs; Equation[D(a) ~ v for (a, v) in zip(avars, vars)]]
+    @set! sys.states = [get_states(sys); avars]
+    @set! sys.defaults = merge(get_defaults(sys), Dict(a => 0.0 for a in avars))
+end

test/odesystem.jl

@@ -331,6 +331,15 @@ D = Differential(t)
 @variables x(t) y(t) z(t)
 D = Differential(t)
 @named sys = ODESystem([D(x) ~ y, 0 ~ x + z, 0 ~ x - y], t, [z, y, x], [])
+asys = add_accumulations(sys)
+@variables accumulation_x(t) accumulation_y(t) accumulation_z(t)
+eqs = [0 ~ x + z
+       0 ~ x - y
+       D(accumulation_x) ~ x
+       D(accumulation_y) ~ y
+       D(accumulation_z) ~ z
+       D(x) ~ y]
+@test sort(equations(asys), by = string) == eqs
 sys2 = ode_order_lowering(sys)
 M = ModelingToolkit.calculate_massmatrix(sys2)
 @test M == Diagonal([1, 0, 0])

test/stream_connectors.jl

@@ -45,13 +45,6 @@ function AdiabaticStraightPipe(; name,
 
     eqns = Equation[]
 
-    #=
-    push!(eqns, port_a.P ~ port_b.P)
-    push!(eqns, 0 ~ port_a.m_flow + port_b.m_flow)
-    push!(eqns, port_b.h_outflow ~ instream(port_a.h_outflow))
-    push!(eqns, port_a.h_outflow ~ instream(port_b.h_outflow))
-    =#
-
     push!(eqns, connect(port_a, port_b))
     sys = ODESystem(eqns, t, vars, pars; name = name)
     sys = compose(sys, subs)
@@ -101,15 +94,49 @@ end
 @named pipe = AdiabaticStraightPipe()
 @named sink = MassFlowSource_h(m_flow_in = -0.01, h_in = 400e3)
 
-streams_a = [n1m1.port_a, pipe.port_a]
-streams_b = [pipe.port_b, sink.port]
-
 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)
 @test_nowarn structural_simplify(n1m1Test)
+@unpack source, port_a = n1m1
+@test sort(equations(expand_connections(n1m1)), by = string) == [0 ~ port_a.m_flow
+       0 ~ source.port1.m_flow - port_a.m_flow
+       source.port1.P ~ port_a.P
+       source.port1.P ~ source.P
+       source.port1.h_outflow ~ port_a.h_outflow
+       source.port1.h_outflow ~ source.h]
+@unpack port_a, port_b = pipe
+@test sort(equations(expand_connections(pipe)), by = string) ==
+      [0 ~ -port_a.m_flow - port_b.m_flow
+       0 ~ port_a.m_flow
+       0 ~ port_b.m_flow
+       port_a.P ~ port_b.P
+       port_a.h_outflow ~ instream(port_b.h_outflow)
+       port_b.h_outflow ~ instream(port_a.h_outflow)]
+@test sort(equations(expand_connections(sys)), by = string) ==
+      [0 ~ n1m1.port_a.m_flow + pipe.port_a.m_flow
+       0 ~ pipe.port_b.m_flow + sink.port.m_flow
+       n1m1.port_a.P ~ pipe.port_a.P
+       pipe.port_b.P ~ sink.port.P]
+@test sort(equations(expand_connections(n1m1Test)), by = string) ==
+      [0 ~ -pipe.port_a.m_flow - pipe.port_b.m_flow
+       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
+       n1m1.port_a.P ~ pipe.port_a.P
+       n1m1.source.port1.P ~ n1m1.port_a.P
+       n1m1.source.port1.P ~ n1m1.source.P
+       n1m1.source.port1.h_outflow ~ n1m1.port_a.h_outflow
+       n1m1.source.port1.h_outflow ~ n1m1.source.h
+       pipe.port_a.P ~ pipe.port_b.P
+       pipe.port_a.h_outflow ~ sink.port.h_outflow
+       pipe.port_b.P ~ sink.port.P
+       pipe.port_b.h_outflow ~ n1m1.port_a.h_outflow
+       sink.port.P ~ sink.P
+       sink.port.h_outflow ~ sink.h_in
+       sink.port.m_flow ~ -sink.m_flow_in]
 
 # N1M2 model and test code.
 function N1M2(; name,
@@ -165,9 +192,6 @@ function N2M2(; name,
     @named port_b = TwoPhaseFluidPort()
     @named pipe = AdiabaticStraightPipe()
 
-    streams_a = [port_a, pipe.port_a]
-    streams_b = [pipe.port_b, port_b]
-
     subs = [port_a; port_b; pipe]
 
     eqns = Equation[]
@@ -190,6 +214,18 @@ eqns = [connect(source.port, n2m2.port_a)
 @named n2m2Test = compose(sys, n2m2, source, sink)
 @test_nowarn structural_simplify(n2m2Test)
 
+# stream var
+@named sp1 = TwoPhaseFluidPort()
+@named sp2 = TwoPhaseFluidPort()
+@named sys = ODESystem([connect(sp1, sp2)], t)
+sys_exp = expand_connections(compose(sys, [sp1, sp2]))
+@test sort(equations(sys_exp), by = string) == [0 ~ -sp1.m_flow - sp2.m_flow
+       0 ~ sp1.m_flow
+       0 ~ sp2.m_flow
+       sp1.P ~ sp2.P
+       sp1.h_outflow ~ ModelingToolkit.instream(sp2.h_outflow)
+       sp2.h_outflow ~ ModelingToolkit.instream(sp1.h_outflow)]
+
 # array var
 @connector function VecPin(; name)
     sts = @variables v(t)[1:2]=[1.0, 0.0] i(t)[1:2]=1.0 [connect = Flow]