What is the status of MTK and vectors of symbolic variables? I am working my way through connectors and components, aided by Google AI Studio, but I have a feeling it is hallucinating some proposals…
To be a little bit more concrete, this is what Google AI Studio proposed for one of my components:
vars = @variables begin
m(t)[1:Ns] = mp, [unit = u"kg"] # Ns is the number of substances
...
end
eqs = [
...
# --- Dynamics ---
[D(m[i]) ~ a.md * xa[i] + b.md * xb[i] for i in 1:Ns]...
]
...
Google AI Studio suggests that the problem lies with this use of a vector of symbolic variables. For the composed model, the following statements work fine:
@named __sys_t23 = System(eqs_t23_v2, t)
@named _sys_t23 = compose(__sys_t23, [m_s1, m_s2, p_s1, p_s2, p_s3, p_se,
v_12, v_23, v_3e,
t_1, t_2, t_3])
However, the final simplification step fails:
sys_t23 = mtkcompile(_sys_t23)
The error message is:
ArgumentError: Cannot call `getindex` on non-array symbolics - found array of shape UnitRange{Int64}[].
Google AI Studio has a few suggestions for how to fix the problem, but none of them work. As a last resort, Google AI Studio suggests to hard-code the number of substances in my system. But that leads to components of quite limited value!?
Without a reproducible example I can’t be sure, but you might be missing some collects of your symbolic arrays.
Yeah, I probably miss some collect… I remember that was important in the past, but I don’t recall exactly how. I will try to make a MRE tomorrow. For now, I have another problem – but will discuss that in another thread.
So… here are a port definition and a simple valve model which are hard-coded to 2 substances. I would prefer to make it general for Ns substances. But my attempt doesn’t work…
Hardcoded to 2 substances
Port/connector:
@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
Valve model:
@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"]
md(t), [unit = u"kg/s", description = "Mass flow from a to 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
md ~ K * Dp/(Dp^2 + p_eps)^(1/4),
a.md ~ md,
# 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
Unsuccessful generalization
Fluid port:
@connector function FluidPort(; name, Ns=1)
vars = @variables begin
p(t), [guess = 0.0, description="Pressure, Pa"]
md(t), [connect = Flow, guess = 0.0, description="Mass flow, kg/s"]
x(t)[1:Ns], [connect = Stream, guess = 1/Ns, description="Mass fractions, -"]
end
System(Equation[], t, vars, []; name = name)
end
Valve:
@component function Valve(; name, Ns=1, Kv=1.0)
@named a = FluidPort(Ns=Ns)
@named b = FluidPort(Ns=Ns)
pars = @parameters begin
# Units must be (mass flow) / (pressure^0.5)
K = Kv, [unit = u"kg/(s*Pa^0.5)", description = "Flow coefficient"]
p_eps = 1e-4, [unit = u"Pa", 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)
# Should use 'p_eps' here for stability
a.md ~ K * Dp/(Dp^2 + p_eps)^(1/4)
# 3. Continuity: No storage
a.md + b.md ~ 0,
# 4. Stream Propagation: Same composition on both sides
[a.x[i] ~ instream(b.x[i]) for i in 1:Ns]...,
[b.x[i] ~ instream(a.x[i]) for i in 1:Ns]...
]
compose(System(eqs, t, vars, pars; name = name), a, b)
end
Which MTK version are you using?
Latest? MTK v11.6 or 11.7
Did you carefully read ModelingToolkit.jl/NEWS.md at master · SciML/ModelingToolkit.jl · GitHub ?
I think better documentation for MTK 11 isn’t yet available.
My college wrote:
But a lot better array handling, and no .~ is needed anymore
The biggest change is that now you have to assign the defaults to the entire variable:
defaults = [pos[i] => 0.0 for i in 1:3] is not possible anymore
defaults = [pos => zeros(3)] is how you have to do it now
1 Like
OK - I have tried to read it. I am not sure I got any wiser, but I will play around with it.
1 Like
Be careful with what the AI proposes, it often helps to provide the latest docs to your AI tool so it doesn’t hallucinate. MTK v11 was released in december, with quite some breaking changes, but with big advantages as well. Mainly way faster constructors and mtkcompile because of type stability by using Symbolics@7 and SymbolicUtils@4.
D(x) ~ y should usually work where both x and y are symbolic arrays. If you are getting problems you can scalarize your equations:
eqs = Symbolics.scalarize.(reduce(vcat, Symbolics.scalarize.(eqs)))
Collect is usually not needed. Can you provide a reproducible example of the code that doesn’t work?
1 Like
OK – I have tried to make a very rudimentary “library” for simulating flow of liquid into a tank, and out of the tank via a valve. This example works if I hard-code the model to 2 substances.
Here are my connector + components:
@connector function FluidPort(; name, Ns=1)
vars = @variables begin
p(t), [guess = 0.0, description="Pressure, Pa"]
md(t), [connect = Flow, guess = 0.0, description="Mass flow, kg/s"]
x(t)[1:Ns], [connect = Stream, guess = 1/Ns, description="Mass fractions, -"]
end
System(Equation[], t, vars, []; name = name)
end
@component function PressureBoundary(; name, Ns=1, p_set=101325.0)
@named a = FluidPort(Ns=Ns)
x_p = ones(Ns) ./ Ns
pars = @parameters begin
p = p_set, [unit=u"Pa"]
x[1:Ns] = x_p, [unit = u"kg/kg"]
end
eqs = [
a.p ~ p
a.x ~ x
]
compose(System(eqs, t, [], pars; name = name), a)
end
@component function PressureReference(; name, Ns=1,p_set=101325.0)
@named a = FluidPort(Ns=Ns)
x_p = ones(Ns) ./ Ns
pars = @parameters begin
p = p_set, [unit=u"Pa"]
x[1:Ns] = x_p, [unit = u"kg/kg"]
end
eqs = [
a.p ~ p
a.md ~ 0.0
a.x ~ x
]
compose(System(eqs, t, [], pars; name = name), a)
end
@component function PumpFixed(; name, md_set=1.0, x_set=[1])
Ns = length(x_set)
@named a = FluidPort(Ns=Ns)
pars = @parameters begin
md = md_set, [unit=u"kg/s", description = "Specified mass flow rate"]
x[1:Ns] = x_set,[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.x ~ x
]
compose(System(eqs, t, [], pars; name = name), a)
end
@component function Valve(; name, Ns=1, Kv=1.0)
@named a = FluidPort(Ns=Ns)
@named b = FluidPort(Ns=Ns)
pars = @parameters begin
# Units must be (mass flow) / (pressure^0.5)
K = Kv, [unit = u"kg/(s*Pa^0.5)", description = "Flow coefficient"]
eps = 1e-4, [unit = u"Pa^2", 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. Valve quasi-turbulent flow
a.md ~ K * Dp/(Dp^2 + eps)^(1/4)
# 3. Continuity: No storage
a.md + b.md ~ 0
# 4. Stream Propagation: Same composition on both sides
a.x ~ instream(b.x)
b.x ~ instream(a.x)
]
compose(System(eqs, t, vars, pars; name = name), a, b)
end
@component function Tank(; name, Ap=1.0, rho_p=1000.0, m0=[1])
Ns = length(m0)
m0 = abs.(m0) .+ 1e-9 # Prevent zero initial mass
@named a = FluidPort(Ns=Ns) # Top Port (Headspace)
@named b = FluidPort(Ns=Ns) # Bottom Port (Drain)
# Gravity from Unitful
gp = 9.81 # 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
m(t)[1:Ns] = m0, [unit = u"kg"]
h(t), [unit = u"m"]
x(t)[1:Ns], [unit = u"kg/kg"]
end
xa = [ifelse(a.md > 0, instream(a.x[i]), x[i]) for i in 1:Ns]
xb = [ifelse(b.md > 0, instream(b.x[i]), x[i]) for i in 1:Ns]
eqs = [
# --- Geometry ---
h ~ sum(m) / (rho * A)
# --- Pressure
b.p ~ a.p + rho * g * h
# --- Mass fractions ---
x ~ m ./ sum(m)
# --- Streams ---
a.x ~ x
b.x ~ x
# --- Dynamics ---
D(m) ~ a.md * xa + b.md * xb
]
compose(System(eqs, t, vars, pars; name = name), a, b)
end
Next, to use this “library” [I’m sure something similar exists in ModelingToolkitStandardLibrary.jl, but I am not interested in that at the moment], I want to simulate a simple, open liquid tank with a gas pressure over the liquid volume. I need to specify the pressure in this gas volume (via component PressureReference), and the mass flow of liquid (two substances) into the liquid volume (component PumpFixed). These two components should be connected to port a of the tank (component Tank). Next, in the bottom of the tank, I connect port b to the inlet of a valve (component Valve, port a), and the outlet of the valve to the pressure boundary PressureBoundary).
Here is the parameters I use, with component instantiation, and connecting equations + generation of model:
# Parameters
# - initial masses in tanks
m0_1 = [2,1] #, [description = "Initial masses tank 1, kg"]
Ns = length(m0_1)
#
# - 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"]
# - gas pressure boundary conditions
p_a = 101325.0 #, [description = "Atmospheric pressure
# - tank parameters
A_1 = 0.3 #, [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)"]
#
@named mp_1 = PumpFixed(; md_set=md_in_1, x_set=x_in_1)
@named pr_1 = PressureReference(; Ns=Ns,p_set=p_a)
@named pb_1 = PressureBoundary(; Ns=Ns,p_set=p_a)
@named va_1e = Valve(; Ns=Ns, Kv=K_12)
@named tnk_1 = Tank(; Ap=A_1, rho_p=rho, m0=m0_1)
#
eqs = [
connect(mp_1.a, pr_1.a, tnk_1.a)
connect(tnk_1.b, va_1e.a)
connect(va_1e.b, pb_1.a)
]
@named __sys = System(eqs, t)
@named _sys = compose(__sys, [mp_1, pr_1, pb_1, va_1e, tnk_1])
sys = mtkcompile(_sys)
Function mtkcompile responds with an error message (Julia v1.12.4, MTK v11.7.2, Symbolics v7.8.0, Windows 11 latest version on an Intel based laptop from 2025).
ArgumentError: Cannot call `getindex` on non-array symbolics - found array of shape UnitRange{Int64}[].
Stacktrace:
[1] throw_not_array(x::SymbolicUtils.SmallVec{UnitRange{Int64}, Vector{UnitRange{Int64}}})
I am sure I do something wrong, but I have not seen a thorough discussion of the new possibilities with symbolic vectors in MTKL.
I gave it a try and couldn’t find the solution, scalarizing the equations did not help, neither did rewriting to [a.x[i] ~ instream(b.x[i]) for i in 1:Ns].... Probably worth submitting an issue. @cryptic.ax maybe you have an idea?
1 Like
OK – the code works now! Thanks for quick fixes to MTK.
1 Like