# MTK -- \`instream\` problem

**URL:** <https://discourse.julialang.org/t/mtk-instream-problem/135171>\
**Category:** Modelling & Simulations\
**Created:** [January 20, 2026, 8:51pm UTC](https://discourse.julialang.org/t/mtk-instream-problem/135171 "2026-01-20T20:51:43Z")\
**Posts on this page:** 2\
**Page:** 1

<div class="post-metadata">

**Author:** ![BLI](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bli/32/37206_2.png) [@BLI](https://discourse.julialang.org/u/BLI)\
**Post date:** [January 20, 2026, 8:51pm UTC](https://discourse.julialang.org/t/mtk-instream-problem/135171/1 "2026-01-20T20:51:43Z")

</div>

I’m trying to understand ModelingToolkit (MTK) connectors, components, and now `Stream` variables in connector.

As an _exercise_ for myself, I am looking at 2 substance, ideal solution liquid tank systems:

 ![image](https://global.discourse-cdn.com/julialang/original/3X/6/0/60c619dc376d11837ba470fbaf476addfdd06fe6.png)

I consider _3 cases_:

1. A single tank case, blue ellipse
2. A two tank case, the green ellipse
3. A three tank case; all tanks

[m is mass, \dot{m} is mass flow rate — `md` in code, h is level, p is pressure, x is mass fraction, A is cross sectional area, \rho is density — `rho` in code.]

I have simulated a “flat” model without using components for (i) the single tank case, and (ii) the three tank case. There is no problem solving them.

I have then created a connector (`FluidPort`) and components and simulated Case 1 (dashed lines below), which exactly matches the flat model of Case 1 (pale, solid lines below).

 ![image](https://global.discourse-cdn.com/julialang/original/3X/a/6/a6dc49b9e5794cf3b92f3788fc849a60de5cbf71.png)

I have also used components to simulate Case 2; no problems.

Finally, I have used components to simulate Case 3 (which works with the flat model), but for Case 3 with components, I get an error message:

```julia
MethodError: no method matching instream(::Float64)
The function `instream` exists, but no method is defined for this combination of argument types.

```

**Question** : Why do I get an error message in the 3 tank case? What am I doing wrong?

* * *

### Building blocks

#### Basic code: package import + components

_Package import_:

```julia
using Unitful
using ModelingToolkit
using ModelingToolkit: t_nounits as t, D_nounits as D
using OrdinaryDiffEq

```

**NB** : the components are not meant to be production quality!! This is a self-learning project. I have specifically hard-coded the components to be valid for 2 substances.

_`FluidPort` connector_:

```julia
@connector function FluidPort(; name)
    vars = @variables begin 
        p(t), [guess = 0.0, description="Pressure, Pa"]
        md(t), [connect = Flow, guess = 0.0, description="Mass flow, kg/s"]
        x1(t), [connect = Stream, guess = 1/2, description="Mass fraction 1, -"]
        x2(t), [connect = Stream, guess = 1/2, description="Mass fraction 2, -"] 
    end
    System(Equation[], t, vars, []; name = name)
end

```

#### One port sources/sinks

_`PressureBoundary`_: specifying pressure and mass fractions – may supply/draw mass, e.g., through a valve.

```julia
@component function PressureBoundary(; name, p_set=101325.0, x_set=[0.5, 0.5])
    @named a = FluidPort()
    x_s = sum(x_set)
    x_s > 0 || error("PressureBoundary $name: Sum of x_set must be greater than zero")
    x_p = x_set ./ x_s  
    pars = @parameters begin
        p = p_set, [unit = u"Pa"]
        x1 = x_p[1], [unit = u"kg/kg"]
        x2 = x_p[2], [unit = u"kg/kg"]
    end
    eqs = [
        a.p ~ p
        a.x1 ~ x1
        a.x2 ~ x2
    ]
    compose(System(eqs, t, [], pars; name = name), a)
end

```

_`PressureReference`_: same as `PressureBoundary`, but with mass flow forced to zero – used to specify pressure.

```julia
@component function PressureReference(; name, p_set=101325.0, x_set=[0.5, 0.5])
    @named a = FluidPort()
    x_s = sum(x_set)
    x_s > 0 || error("PressureBoundary $name: Sum of x_set must be greater than zero")
    x_p = x_set ./ x_s  
    pars = @parameters begin
        p = p_set, [unit = u"Pa"]
        x1 = x_p[1], [unit = u"kg/kg"]
        x2 = x_p[2], [unit = u"kg/kg"]
    end
    eqs = [
        a.p ~ p,
        a.md ~ 0, # THE KEY DIFFERENCE from PressureBoundary: Blocks mass flow
        a.x1 ~ x1,
        a.x2 ~ x2
    ]
    compose(System(eqs, t, [], pars; name = name), a)
end

```

_`MassFlowPump`_: an ideal pump (no pressure dependence) to provide specified total mass flow rate and specified mass fractions.

```julia
@component function MassFlowPump(; name, md_set=1.0, x_set=[0.5, 0.5])
    @named a = FluidPort()
    x_s = sum(x_set)
    x_s > 0 || error("MassFlowPump $name: Sum of x_set must be greater than zero")
    x_p = x_set ./ x_s
    pars = @parameters begin
        md = md_set, [unit=u"kg/s", description = "Specified mass flow rate"]
        x1 = x_p[1], [unit = u"kg/kg", description = "Specified mass fractions"]
        x2 = x_p[2], [unit = u"kg/kg", description = "Specified mass fractions"]
    end
    eqs = [
        # Inflow to system means outflow from component (negative md)
        a.md ~ -md,
        # Define the stream variables for the port
        a.x1 ~ x1,
        a.x2 ~ x2
    ]
    compose(System(eqs, t, [], pars; name = name), a)
end

```

#### Two port components

_`Valve`_: essentially a turbulent flow valve, with flow proportional to the square root of the pressure difference.

```julia
@component function Valve(; name, Kv=1.0)
    @named a = FluidPort()
    @named b = FluidPort()
    pars = @parameters begin
        # Units must be (mass flow) / (pressure^0.5)
        K = Kv, [description = "Flow coefficient"]
        p_eps = 1e-4, [description = "Regularization factor"]
    end
    vars = @variables begin
        Dp(t), [unit = u"Pa", description = "Pressure drop a-b"]
    end
    eqs = [
        # 1. Pressure drop
        Dp ~ a.p - b.p,
        # 2. Flow Equation: Square Root Law (turbulent flow)
        # - use 'p_eps' here for stability
        a.md ~ K * sign(Dp) * sqrt(abs(Dp) + p_eps),
        # 3. Continuity: No storage
        a.md + b.md ~ 0,
        # 4. Stream Propagation: same composition on both sides
        a.x1 ~ instream(b.x1),
        a.x2 ~ instream(b.x2),
        b.x1 ~ instream(a.x1),
        b.x2 ~ instream(a.x2)
    ]
    compose(System(eqs, t, vars, pars; name = name), a, b)
end

```

_`Tank`_: simple liquid tank with specified cross sectional area, and unlimited possible liquid level. The pressure of gas above liquid is specified (`PressureReference`) where also liquid may be pumped in (port `a`) and the pressure at the bottom (port `b`) is this gas level + the hydrostatic pressure of the liquid.

```julia
@component function Tank(; name, Ap=1.0, rho_p=1000.0, m0=[1,1])
    @named a = FluidPort() # Top Port (Headspace)
    @named b = FluidPort() # Bottom Port (Drain)
    # Gravity from Unitful
    gp = ustrip(u"m/s^2", Unitful.gn)
    # Symbolic parameters   
    pars = @parameters begin
        A = Ap, [unit = u"m^2"]
        rho = rho_p, [unit = u"kg/m^3"]
        g = gp, [unit = u"m/s^2"]
    end
    # Symbolic variables
    vars = @variables begin
        m1(t) = m0[1], [unit = u"kg"]
        m2(t) = m0[2], [unit = u"kg"]
        h(t), [unit = u"m"]
        x1(t), [unit = u"kg/kg"]
        x2(t), [unit = u"kg/kg"]
        x1a(t), [unit = u"kg/kg"]
        x2a(t), [unit = u"kg/kg"]
        x1b(t), [unit = u"kg/kg"]
        x2b(t), [unit = u"kg/kg"]
    end
    eqs = [
        #
        # --- Level, assuming incompressible, ideal solution liquid
        h ~ (m1 + m2) / (rho * A),
        # --- Pressure Logic
        # ... pressure at bottom = pressure at the top + liquid head.
        b.p ~ a.p + rho * g * h,
        # --- Mass fraction
        x1 ~ m1 / (m1 + m2),
        x2 ~ m2 / (m1 + m2),
        # --- Actual mass fractions at ports
        x1a ~ ifelse(a.md > 0, instream(a.x1), x1),
        x2a ~ ifelse(a.md > 0, instream(a.x2), x2),
        x1b ~ ifelse(b.md > 0, instream(b.x1), x1),
        x2b ~ ifelse(b.md > 0, instream(b.x2), x2),
        # --- Setting stream values ---
        a.x1 ~ x1,
        a.x2 ~ x2,
        b.x1 ~ x1,
        b.x2 ~ x2,
        # --- Dynamics ---
        D(m1) ~ a.md * x1a + b.md * x1b,
        D(m2) ~ a.md * x2a + b.md * x2b
    ]
    compose(System(eqs, t, vars, pars; name = name), a, b)
end

```

### Case 1: single tank

Component model:

```julia
# Parameters
# - initial masses in tanks
m0 = [2,1] #, [description = "Initial masses tank 1, kg"]
#
# - inlet flow rates and compositions
md_in = 2.0 #, [description = "Inlet mass flow tank 1, kg/s"]
x_in = [0.8,0.2] #, [description = "Inlet mass fractions tank 1, kg/kg"]
# - gas pressure boundary conditions
p_a = 101325.0 #, [description = "Atmospheric pressure
# - tank parameters
A = 0.3 #, [description = "Tank cross section area, m^2"]
rho = 1000.0 #, [description = "Fluid density, kg/m^3"]
# - valve flow coefficients
K = 0.3 #, [description = "Flow coefficient tank 1 to tank 2, kg/(s*Pa^0.5)"]
#
# Building the system
#
@named md_s = MassFlowPump(; md_set=md_in, x_set=x_in)
@named p_si = PressureReference(; p_set=p_a)
@named p_se = PressureBoundary(; p_set=p_a)
@named v_e = Valve(; Kv=K)
@named tnk = Tank(; Ap=A, rho_p=rho, m0=m0)

eqs_t1_2 = [
    connect(md_s.a, p_si.a, tnk.a)
    connect(tnk.b, v_e.a)
    connect(v_e.b, p_se.a)
]

@named __sys_t1_2 = System(eqs_t1_2, t)
@named _sys_t1_2 = compose(__sys_t1_2, [md_s,p_si, p_se,v_e,tnk])  

```

Compare this to Case 1 (blue ellipse). The pumped mass flow (`md_s` has a single port, `a`) is connected with the pressure reference (`p_si`; single port, `a`) at the top (port `a`) of the tank `tnk`.

At the bottom of the tank `tnk` (port `b`), the tank is connected to the inlet of the effluent valve (`v_e`). Finally, the outlet of the valve `v_e` is connected to the pressure boundary `p_se` (which may receive mass flow).

The symbolic model is simplified via `mtkcompile`, a numeric problem is generated via `ODEProblem`, and the solution is found:

```julia
sys_t1_2 = mtkcompile(_sys_t1_2)

tspan = (0.0, 10)
prob_t1_2 = ODEProblem(sys_t1_2, [], tspan)

sol_t1_2 = solve(prob_t1_2)

```

_Note_: the `mtkcompile(_sys_t1_2)` statement produce warnings:

```julia
Warning: Internal error: Variable tnk₊a₊md(t) was marked as being in 0 ~ -tnk₊x1a(t) + ifelse(tnk₊a₊md(t) > 0, ModelingToolkitBase.instream_rt(Val{2}(), Val{0}(), md_s₊a₊md(t), p_si₊a₊md(t), md_s₊a₊x1(t), p_si₊a₊x1(t)), tnk₊x1(t)), but was actually zero
└ @ ModelingToolkitTearing ...

```

I didn’t worry about these, but it would be nice if there was a way to get rid of them.

### Case 2: two tanks

Component model (green circle):

```julia
# Parameters
# - initial masses in tanks
m0_1 = [2,1] #, [description = "Initial masses tank 1, kg"]
m0_2 = [5,8] #, [description = "Initial masses tank 2, kg"]
#
# - inlet flow rates and compositions
md_in_1 = 2.0 #, [description = "Inlet mass flow tank 1, kg/s"]
x_in_1 = [0.8,0.2] #, [description = "Inlet mass fractions tank 1, kg/kg"]
md_in_2 = 2.5 #, [description = "Inlet mass flow tank 2, kg/s"]
x_in_2 = [0.4,0.6] #, [description = "Inlet mass fractions tank 2, kg/kg"]
# - gas pressure boundary/reference conditions
p_a = 101325.0 #, [description = "Atmospheric pressure
# - tank parameters
A_1 = 0.3 #, [description = "Tank cross section area, m^2"]
A_2 = 0.6 #, [description = "Tank cross section area, m^2"]
rho = 1000.0 #, [description = "Fluid density, kg/m^3"]
# - valve flow coefficients
K_12 = 0.1 #, [description = "Flow coefficient tank 1 to tank 2, kg/(s*Pa^0.5)"]
K_e = 0.4 #, [description = "Flow coefficient tank 2 to tank 3, kg/(s*Pa^0.5)"]

#
@named m_s1 = MassFlowPump(; md_set=md_in_1, x_set=x_in_1)
@named m_s2 = MassFlowPump(; md_set=md_in_2, x_set=x_in_2)
@named p_s1 = PressureReference(; p_set=p_a)
@named p_s2 = PressureReference(; p_set=p_a)
@named p_e = PressureBoundary(; p_set=p_a)
@named v_12 = Valve(; Kv=K_12)
@named v_e = Valve(; Kv=K_e)
@named t_1 = Tank(; Ap=A_1, rho_p=rho, m0=m0_1)
@named t_2 = Tank(; Ap=A_2, rho_p=rho, m0=m0_2)

eqs_t2 = [
    connect(m_s1.a, p_s1.a, t_1.a)
    connect(m_s2.a, p_s2.a, t_2.a)
    connect(t_1.b, v_12.a)
    connect(v_12.b, t_2.b,v_e.a)
    connect(v_e.b, p_e.a)
]

@named __sys_t2 = System(eqs_t2, t)
@named _sys_t2 = compose(__sys_t2, [m_s1, m_s2, p_s1, p_s2, p_e, 
                                        v_12, v_e, t_1, t_2])
sys_t2 = mtkcompile(_sys_t2)

tspan = (0.0, 10)
prob_t2 = ODEProblem(sys_t2, [], tspan)

sol_t2 = solve(prob_t2)

```

Again, same warning after `mtkcompile`. But the solution is found, and makes perfect sense:

 ![image](https://global.discourse-cdn.com/julialang/original/3X/0/8/082a7e3b38af63f83725a24ec3aa5717db097356.png)

Because of a change in which tank has the highest level (highest bottom pressure), the flow direction \dot{m}\_{1,2} changes direction, with consequences for which mass fractions are transferred:

 ![image](https://global.discourse-cdn.com/julialang/original/3X/4/e/4eb8ffc96c0ffcd24c0a451ada60ee5035379a72.png)

### Case 3: three tanks

Component model:

```julia
# Parameters
# - initial masses in tanks
m0_1 = [2,1] #, [description = "Initial masses tank 1, kg"]
m0_2 = [5,8] #, [description = "Initial masses tank 2, kg"]
m0_3 = [3,5] #, [description = "Initial masses tank 3, kg"]
#
# - inlet flow rates and compositions
md_in_1 = 2.0 #, [description = "Inlet mass flow tank 1, kg/s"]
x_in_1 = [0.8,0.2] #, [description = "Inlet mass fractions tank 1, kg/kg"]
md_in_2 = 2.5 #, [description = "Inlet mass flow tank 2, kg/s"]
x_in_2 = [0.4,0.6] #, [description = "Inlet mass fractions tank 2, kg/kg"]
# - gas pressure boundary conditions
p_a = 101325.0 #, [description = "Atmospheric pressure
# - tank parameters
A_1 = 0.3 #, [description = "Tank cross section area, m^2"]
A_2 = 0.6 #, [description = "Tank cross section area, m^2"]
A_3 = 0.4 #, [description = "Tank cross section area, m^2"]
rho = 1000.0 #, [description = "Fluid density, kg/m^3"]
# - valve flow coefficients
K_12 = 0.1 #, [description = "Flow coefficient tank 1 to tank 2, kg/(s*Pa^0.5)"]
K_23 = 0.4 #, [description = "Flow coefficient tank 2 to tank 3, kg/(s*Pa^0.5)"]
K_3e = 0.3 #, [description = "Flow coefficient effluent tank 3, kg/(s*Pa^0.5)"]

#
@named mp_s1 = MassFlowPump(; md_set=md_in_1, x_set=x_in_1)
@named mp_s2 = MassFlowPump(; md_set=md_in_2, x_set=x_in_2)
@named p_r1 = PressureReference(; p_set=p_a)
@named p_r2 = PressureReference(; p_set=p_a)
@named p_r3 = PressureReference(; p_set=p_a)
@named p_be = PressureBoundary(; p_set=p_a)
@named va_12 = Valve(; Kv=K_12)
@named va_23 = Valve(; Kv=K_23)
@named va_3e = Valve(; Kv=K_3e)
@named tnk_1 = Tank(; Ap=A_1, rho_p=rho, m0=m0_1)
@named tnk_2 = Tank(; Ap=A_2, rho_p=rho, m0=m0_2)
@named tnk_3 = Tank(; Ap=A_3, rho_p=rho, m0=m0_3)

eqs_t3 = [
    connect(mp_s1.a, p_r1.a, tnk_1.a)
    connect(mp_s2.a, p_r2.a, tnk_2.a)
    connect(tnk_1.b, va_12.a)
    connect(va_12.b, tnk_2.b, va_23.a)
    connect(va_23.b, tnk_3.a, p_r3.a)
    connect(tnk_2.b, va_3e.a)
    connect(va_3e.b, p_be.a)
]

@named __sys_t3 = System(eqs_t3, t)
@named _sys_t3 = compose(__sys_t3, [mp_s1, mp_s2, p_r1, p_r2, p_r3, p_be,
                                       va_12, va_23, va_3e,
                                       tnk_1, tnk_2, tnk_3])

sys_t3 = mtkcompile(_sys_t3)

```

Same warning as before, but the simplified symbolic model is created. Next, creating numeric problem:

```julia
tspan = (0.0, 10)
prob_t3 = ODEProblem(sys_t3, [], tspan)

```

This works, and creates the expected output:

```julia
ODEProblem with uType Vector{Float64} and tType Float64. In-place: true
Initialization status: FULLY_DETERMINED
Non-trivial mass matrix: false
timespan: (0.0, 10.0)
u0: 6-element Vector{Float64}:
 5.0
 3.0
 8.0
 5.0
 1.0
 2.0

```

Finally, attempting to solve the problem:

```julia
sol_t3 = solve(prob_t3)

```

crashes with the following error message:

```julia
MethodError: no method matching instream(::Float64)
The function `instream` exists, but no method is defined for this combination of argument types.

Closest candidates are:
  instream(!Matched::SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymReal})
   @ ModelingToolkitBase C:\Users\Be...

```

---

<div class="post-metadata">

**Author:** ![BLI](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bli/32/37206_2.png) [@BLI](https://discourse.julialang.org/u/BLI)\
**Post date:** [January 21, 2026, 3:08pm UTC](https://discourse.julialang.org/t/mtk-instream-problem/135171/2 "2026-01-21T15:08:25Z")

</div>

Hm. I found a tiny, but critical bug in my connections. I’ll check that it works before updating.

* * *

OK here was my problem:

 ![image](https://global.discourse-cdn.com/julialang/original/3X/2/0/20e3c23fc4f2982f3e912deaddc9bfbf6ab0305c.png)

I had connected my valve for flow \dot{m}\_{1,2} _simultaneously_ to the bottom port of Tank 2 and to the flow \dot{m}\_{2,3}. This works if the flow \dot{m}\_{1,2} is negative, i.e., level in tank 2 is higher than the level in tank 1: then both flows \dot{m}\_{1,2} and \dot{m}\_{2,3} _leave_ tank 2, and they take the mass fractions from tank 2.

However, when \dot{m}\_{1,2} changes direction and becomes positive, this creates a connection where mass is added via \dot{m}\_{1,2} and tank 2, and leaves via \dot{m}\_{2,3}. This implies that the mass fractions leaving with \dot{m}\_{2,3} are computed as a weighted mixture of the mass fraction from tank 1 (via \dot{m}\_{1,2} and the mass fraction from tank 2. But it should be the mass fraction from tank 2.

Because of this, it is necessary to have _two_ ports at the bottom of tank 2, and connect \dot{m}\_{1,2} to one of those ports, and \dot{m}\_{2,3}$ to the other port. That makes the code compute the correct results.