Connectors

[YingboMa]

connect(ps..., ::Val{:flow}) = ...
connect(ps...) = connect(ps..., Val(promote_ connection_properties(ps...)))

We also need arbitrary metadata from SymbolicUtils.

[ValentinKaisermayer]

Any ideas for dispatching the connection equation generation?
In the example there is already the problem, that there are essentially two physical domains and hence two types of Ports (electrical and thermal) but the type is ODESystem in both cases. In the example two methods are used connect_pin and connect_heat which is not really nice - I feel.
One could generalize this by using effort and flow variables for the ports (i.e. potential and current, temperature and heat flow rate, pressure and volume or mass flow rate). The connection would than ensure that the effort is equal for all ports and that the flow is conserved.
However, not all physical domain models will just have those two variables for the ports. For instance modelling thermal-fluid systems requires additional variables defined for each port. See the Modelica specification for an example on StreamConnectors and FluidPorts.

[ValentinKaisermayer]

I see this is already being worked on. https://mtk.sciml.ai/dev/basics/ContextualVariables/#Flow-Variables-(TODO)-1

[ChrisRackauckas]

Yes, indeed. It's a fully extendable metadata system:

https://github.com/JuliaSymbolics/Symbolics.jl/pull/117/files#diff-ddf745d665b6822a8aea0de26378b4fba95e0cafaa53a858e3a881976eecc599R17

and the next thing to do is a dispatchable connection generation system on top of that. That all should come fairly shortly. Indeed the reasoning is because thinking about only physical models is too shallow: there's all of physiology, psychietry, etc. to model as well, and thinking they all fit into flow and stream would be short sighted IMO.

[ValentinKaisermayer]

Nice!

In Modelica there is a notion of an insidePort and outsidePort. For instance if one wants to combine multiple components in a new one, also some internal connections will be needed. It would be nice to also automate them. However, e.g. the balance equations for those internal connections will be different than for outside connections.

See the following example of a three way mixing valve:

using ModelingToolkit, Plots, DifferentialEquations

@parameters t

"""Port for hydraulic systems"""
function HydraulicPort(; name)
    @variables begin
        p(t) # [Pa] Pressure
        mdot(t) # [kg/s] Mass flow rate
    end
    return ODESystem(Equation[], t, [p, mdot], []; name=name, default_u0=[p => 1e5, mdot => 0.0])
end

"""Infinite reservoir; sets BC for hydraulic system"""
function Reservoir(; name, p=1e5)
    p_val = p # [Pa] Pressure of the reservoir
    @parameters p
    @named a = HydraulicPort()
    eqs = [
        a.p ~ p
    ]
    return ODESystem(eqs, t, [], [p]; systems=[a], default_p=Dict(p => p_val), name=name)
end

# Smooth approximation of sign(x) * sqrt(abs(x)), see: https://www.maplesoft.com/documentation_center/online_manuals/modelica/Modelica_Fluid_Utilities.html#Modelica.Fluid.Utilities.regRoot
regRoot(x, delta) = x / (x^2 + delta^2)^0.25
@register regRoot(x, delta)

"""Variable flow resistance in a hydraulic network. Input variable: `A`"""
function VariableFlowResistance(; name, A_max=1e-4, A_leakage=1e-9)
    @variables A(t) # [m^2] Effective area of valve
    @named a = HydraulicPort()
    @named b = HydraulicPort()
    
    rho = 998.0 # [kg/m^3] Density of water
    p_atm = 101325.0 # [Pa] Atmospheric pressure
    B_lam = 0.999 # [-] Laminar flow pressure ratio
    
    p_cr = (a.p + b.p + 2 * p_atm) / 2 * (1 - B_lam) # critical pressure
    Δp = b.p - a.p

    eqs = [
        a.mdot + b.mdot ~ 0 # Mass balance: mdot_a + mdot_b = 0 
        b.mdot ~ max(min(A, A_max), A_leakage) * sqrt(2 * rho) * regRoot(Δp, p_cr) # Momentum balance: mdot_b == Cd * A * sqrt(2 * rho) * sqrt(abs(Δp)) * sign(Δp)
    ]

    return ODESystem(eqs, t, [A], []; systems=[a, b], default_u0=Dict(A => A_max / 2), name=name)
end

"""Three-way mixing valve"""
function MixingValve(; name, A_max=1e-4, A_leakage=1e-9, y=0.5)
    y_val = y
    @parameters y # opening of the valve [0,1]
    @named valve_1 = VariableFlowResistance(; A_max, A_leakage) 
    @named valve_2 = VariableFlowResistance(; A_max, A_leakage)
    @named a = HydraulicPort()
    @named b = HydraulicPort()
    @named c = HydraulicPort()
    eqs = [
        valve_1.A ~ y * A_max
        valve_2.A ~ (1 - y) * A_max
        valve_1.a.mdot - a.mdot ~ 0 # Inflowing mass flow rates are always positive
        a.p ~ valve_1.a.p
        valve_2.a.mdot - b.mdot ~ 0 # Inflowing mass flow rates are always positive
        b.p ~ valve_2.a.p
        valve_1.b.mdot + valve_2.b.mdot - c.mdot ~ 0 # Inflowing mass flow rates are always positive
        c.p ~ valve_1.b.p
        c.p ~ valve_2.b.p
    ]
    return ODESystem(eqs, t, [], [y]; systems=[a, b, c, valve_1, valve_2], default_p=Dict(y => y_val), name=name)
end

"""Generates the connection equations between `HydraulicPorts`"""
function connect(ports...)
    # Mass balance
    eqs = [
        sum(port -> port.mdot, ports) ~ 0.0, # sum(mdot) = 0
    ]

    # Pressure balance
    for port in ports[2:end]
        push!(eqs, ports[1].p ~ port.p) # p_1 = p_2, p_1 = p_3,...
    end

    return eqs
end

@named r1 = Reservoir(; p=2e5)
@named r2 = Reservoir(; p=2e5)
@named mixing_valve = MixingValve(; A_max=0.7 * 1e-4, y=1/3)
@named r3 = Reservoir(; p=1e5)
connections = [
    connect(r1.a, mixing_valve.a)
    connect(r2.a, mixing_valve.b)
    connect(mixing_valve.c, r3.a)
]
@named model = ODESystem(connections, t; systems=[r1, r2, mixing_valve, r3])
sys = structural_simplify(model)
tspan = (0.0, 10.0)
prob = ODEProblem(sys, Pair[], tspan)

plot()
plot!(0:10, t->prob.f.observed(mixing_valve.a.mdot, prob.u0, prob.p, t), label="mixing_valve.a.mdot")
plot!(0:10, t->prob.f.observed(mixing_valve.b.mdot, prob.u0, prob.p, t), label="mixing_valve.b.mdot")
plot!(0:10, t->prob.f.observed(mixing_valve.c.mdot, prob.u0, prob.p, t), label="mixing_valve.c.mdot")

In this example I tried to re-use the VariableFlowResistance component in the MixingValve component.
However, since inflowing mass flow rates are by definition positive the mass balance equations inside the mixing valve are different than the mass balance equations outside the mixing valve.
This is also specified in the Modelica specification I posted earlier:

// Standard connection equation for flow variables
0 = sum(m_j.c.m_flow for j in 1:N) + sum(-ck.m_flow for k in 1:M);

where the mass flow rates from outside connectors (ports) are counted negative and hence connect(ports...) can not be used directly.
So metadata (val, connect, unit) could be not enough.
Or one has to define struct FlowExternal <: AbsstractFlow end and struct FlowInternal <: AbsstractFlow end and dispatch the flow balance equation generation based on that.

[crlaugh]

I'm not familiar with insidePort or outsidePort in Modelica. Can you share the equivalent Modelica code for this model?

[ValentinKaisermayer]

As far as I know there are no separate types but it is inferred during compile time of the model. See: https://specification.modelica.org/master/connectors-and-connections.html

[ValentinKaisermayer]

You can see this also in the implementation of TinyModia.jl https://github.com/ModiaSim/TinyModia.jl/blob/666a92cf91e93e93391f64e76a09fb3ca3577406/src/NamedTupleModels.jl#L101
Otherwise, compositions in submodels are not possible.

[ValentinKaisermayer]

My take on this:

struct Effort end
struct Flow end

function Pin(; name)
    @parameters t
    @variables begin
        v(t), [connect=Effort] # [V] Potential
        i(t), [connect=Flow] # [A] Current
    end
    return ODESystem(Equation[], t, [v, i], []; name=name, defaults=[v => 0.0, i => 0.0])
end

function get_vars_of_type(sys, T)
    sts = ModelingToolkit.get_states(sys)
    vars = filter(s->ModelingToolkit.getmetadata(s, ModelingToolkit.VariableConnectType) == T, sts)
    return getproperty.(Ref(sys), ModelingToolkit.getname.(vars))
end

function connect(connectors...;
    is_inside=fill(true, length(connectors))) # Which connector are inside connectors?

    eqs = [ # Flow balance
        (0 .~ sum([(i_ ? 1 : -1) * get_vars_of_type(c, Flow) for (c, i_) in zip(connectors, is_inside)]))...
    ]

    for (i, c) in enumerate(connectors[2:end]) # Effort equality
        push!(eqs, (get_vars_of_type(connectors[1], Effort) .~ get_vars_of_type(c, Effort))...)
    end
    
    return eqs
end

Which can handle composite components:

function RealCapacitor(; name, C=1.0, R_esr=0.1, v0=0.0)
    C_val, R_esr_val = C, R_esr

    @parameters t
    @named p = Pin()
    @named n = Pin()
    @named R_esr = Resistor(; R=R_esr_val)
    @named C = Capacitor(; C=C_val, v0=v0)
    connections = [
        connect(p, R_esr.p, is_inside=[false, true])
        connect(R_esr.n, C.p, is_inside=[true, true])
        connect(C.n, n, is_inside=[false, true])
    ]
    return ODESystem(connections, t, [], []; systems=[p, n, R_esr, C], name=name)
end

@named R = Resistor(; R=1)
@named C = RealCapacitor(; C=1, R_esr=0.1, v0=0.0)
@named source = ConstantVoltage(; V=1)
@named ground = Ground()
connections = [
    connect(source.p, R.p)
    connect(R.n, C.p)
    connect(C.n, source.n, ground.p)
]
@parameters t
@named model = ODESystem(connections, t; systems=[R, C, source, ground])
sys = structural_simplify(model)
prob = ODAEProblem(sys, Pair[], (0, 10.0))
sol = solve(prob, Tsit5())
plot(sol; vars=[R.v, C.R_esr.v])

A bit wonky and the user has to know exactly what he/she does but it works.
It would be nice to know that e.g. for this connection

connect(p, R_esr.p)

the first system is a top-level component and has no parent and the second has a parent. Maybe this could be used to infer if the connector is an inside or outside connector.

However, there are other difficulties.
For a Stream connector one would have to know to which other connectors it is connected.
In Modelica inStream(c) and actualStream(c) can be used inside components.
The actualStream(c) is easy, it just depends on the specific connector.
But the inStream(c) depends on all connected connectors.
This could be resolved by simply adding the corresponding in stream variable to the connector itself.

"""Port for a thermo-fluid system"""
function ThermoFluidPort(; name)
    @parameters t
    @variables begin
        p(t), [connect=Effort] # [Pa] Pressure at the port
        mdot(t), [connect=Flow] # [kg/s] Mass flow rate at the port. mdot > 0 means inward flow into the component
        T_outflow(t), [connect=Stream] # [K] T if mdot < 0 i.e. the stream flows out of the component (Modelica uses sp. enthalpy)
        T_in_stream(t), [connect=InStream] # [K] T in stream direction (Not part of the Modelica specification of the stream connector)
    end
    return ODESystem(Equation[], t, [p, mdot, T_outflow, T_in_stream], []; 
        name=name, 
        defaults=[p => 1e5, mdot => 0, h_outflow => 273.15 + 20, T_in_stream => 273.15 + 20])
end

When the connection equations are generated one would know all other connectors and can build the equation for T_in_stream according to: https://specification.modelica.org/master/stream-connectors.html