# MTK & vector equations -- get error message on "axes of broadcast not known"

OK… I’m finally trying to learn MTK by re-implementing old models from Modelica (i.e., the models work).

Now, I’m doing a manual discretization of an advection PDE/simple gas pipe. It is not super realistic, but that is not the point.

• I want to find the pressure in each volume
• I assume that inlet pressure at `x=0` is known, and I “stupidly” assume that the fluid velocity is the same throughout the pipe.

Here is my Julia/MTK code attempt:

``````# Number of volumes
N = 10
#
@variables t m(t)[1:N] n(t)[1:N] md_i(t)[1:N] md_e(t)[1:N] nd_i(t)[1:N] nd_e(t)[1:N] p(t)[1:N] p_i(t)[1:N] Vd(t) p0(t) v(t)
@parameters R M A L T dV dx
Dt = Differential(t)

@register_symbolic u_p(t)
@register_symbolic u_v(t)

eqs = [ Dt(m) ~ md_i - md_e,  # mass balance
m ~ n*M,              # mass vs. mole
p*dV ~ n*R*T,         # ideal gas law, constant temperature
p_i ~ p0,          # boundary condition for pressure
p_i[2:N] ~ p[1:N-1],  # inlet pressure for volume 2 equals outlet pressure from vol 1, etc.
Vd ~ v*A,             # volumetric flow rate vs. velocity
nd_i ~ Vd*p_i/(R*T),  # flow version of ideal gas law
nd_e ~ Vd*p/(R*T),    # -- " --
md_i ~ nd_i*M,        # mass flow vs. molar flow
md_e ~ nd_e*M,        # -- " --
p0 ~ u_p(t),          # setting inlet pressure equal to a user defined time function
v ~ u_v(t)]           # setting velocity equal to a user defined time function

@named pipe = ODESystem(eqs)
``````

This results in an error message:

``````xes of (broadcast(~, Differential(t)((m(t))[1:10]), broadcast(-, md_i(t), md_e(t)))) not known

Stacktrace:
 error(s::String)
@ Base .\error.jl:35
 axes(A::Symbolics.Arr{Any, 1})
...
``````

*Question: What is my error (or list of errors :-o )?

I believe the problem is that all variables in this equation are arrays, and `Dt` and `~` are scalar operators. If you broadcast these operators, it might work

``````Dt.(m) .~ md_i .- md_e
``````

Normally, I `collect(m)` all array variables immediately after creating them since there are a lot of issues with them. The call to `collect` will transform the array variable to an array of variables.

This executes so far at least

``````N = 10
#
@variables t m(t)[1:N] n(t)[1:N] md_i(t)[1:N] md_e(t)[1:N] nd_i(t)[1:N] nd_e(t)[1:N] p(t)[1:N] p_i(t)[1:N] Vd(t) p0(t) v(t)

m, n, md_i, md_e, nd_i, nd_e, p, p_i = collect.((m, n, md_i, md_e, nd_i, nd_e, p, p_i))
@parameters R M A L T dV dx
Dt = Differential(t)

@register_symbolic u_p(t)
@register_symbolic u_v(t)

eqs = [ Dt.(m) .~ md_i - md_e;  # mass balance
m .~ n*M;              # mass vs. mole
p*dV .~ n*R*T;         # ideal gas law; constant temperature
p_i ~ p0;          # boundary condition for pressure
p_i[2:N] .~ p[1:N-1];  # inlet pressure for volume 2 equals outlet pressure from vol 1; etc.
Vd ~ v*A;             # volumetric flow rate vs. velocity
nd_i .~ Vd*p_i/(R*T);  # flow version of ideal gas law
nd_e .~ Vd*p/(R*T);    # -- " --
md_i .~ nd_i*M;        # mass flow vs. molar flow
md_e .~ nd_e*M;        # -- " --
p0 ~ u_p(t);          # setting inlet pressure equal to a user defined time function
v ~ u_v(t)]
``````
1 Like

Still no luck. HOWEVER: if I changed my comma-separation of elements in `eqs` to semi-colon separation, it works.

• Question: why does “;” work, while “,” does not work?

Yes, that’s why my example had semi colons `Dt.(m) .~ md_i - md_e` creates an array, and `;` is the array-concatenation operator. The comma is used to create a vector, so using this would create a vector of vectors instead of one long vector.

AHA! I see Hm…

``````pipe_simp = structural_simplify(pipe)
``````

works fine.

``````p_a = 1.03e5 # Pa

p0_const(t) = 10*pa
p0_step(t) = t<25 ? 10*pa : 12*pa

v_const(t) = 0.2 # m/s
v_step(t) = t<25 ? 0.2 : 0.3
``````

works (obviously).

``````# Parameters
par = [R => 8.31, M => 28, A => 5, L => 10, dV => A*L/N, T => 300]
# Operating conditions
p0 = 1.2*p_a*ones(N)
x0 = [m .=> p0*(A*L/N)/(R*T)*M] # Specifying initial state(s)
u_p(t) = p0_const(t)   # Specifying input function
u_v(t) = v_const(t)    # Specifying input function
# Time span
tspan = (0.0,150)   # Specifying tuple of initial and final time
``````

… no error messages.

But the following:

``````prob = ODEProblem(pipe_simp, x0, tspan, par)
``````

fails miserably.

• Question: is it still a problem with array variable vs. array of variables?

Would you mind posting the full code in a copy-pasteable format to reproduce the error?

OK – here is the complete code:

``````using ModelingToolkit
using DifferentialEquations

# Number of volumes
N = 10
#
@variables t m(t)[1:N] n(t)[1:N] md_i(t)[1:N] md_e(t)[1:N] nd_i(t)[1:N] nd_e(t)[1:N] p(t)[1:N] p_i(t)[1:N] Vd(t) pi0(t) v(t)
#
# Converting symbolic array variables into arrays of variables
m, n, md_i, md_e, nd_i, nd_e, p, p_i = collect.( (m, n, md_i, md_e, nd_i, nd_e, p, p_i) )
@parameters R M A L T dV
Dt = Differential(t)

@register_symbolic u_p(t)
@register_symbolic u_v(t)

eqs = [ Dt.(m) .~ md_i - md_e;
m .~ n*M;
p*dV .~ n*R*T;
p_i ~ pi0;
p_i[2:N] .~ p[1:N-1];
Vd ~ v*A;
nd_i .~ Vd*p_i/(R*T);
nd_e .~ Vd*p/(R*T);
md_i .~ nd_i*M;
md_e .~ nd_e*M;
pi0 ~ u_p(t);
v ~ u_v(t) ]

@named pipe = ODESystem(eqs)

pipe_simp = structural_simplify(pipe)

p_a = 1.03e5 # Pa

p0_const(t) = 10*p_a
p0_step(t) = t<25 ? 10*p_a : 12*p_a

v_const(t) = 0.2 # m/s
v_step(t) = t<25 ? 0.2 : 0.3

# Parameters
par = [R => 8.31, M => 28, A => 5, L => 10, dV => A*L/N, T => 300]
# Operating conditions
p0 = 1.2*p_a*ones(N)
x0 = [m .=> p0*(A*L/N)/(R*T)*M] # Specifying initial state(s)
u_p(t) = p0_const(t)   # Specifying input function
u_v(t) = v_const(t)    # Specifying input function
# Time span
tspan = (0.0,150)   # Specifying tuple of initial and final time

prob = ODEProblem(pipe_simp, x0, tspan, par)
``````

Everything works fine until the last line. The error message comes from the last line.

This

``````x0 = [m .=> p0*(A*L/N)/(R*T)*M]
``````

creates a vector of vectors, simply

``````x0 = m .=> p0*(A*L/N)/(R*T)*M
``````

is what you want

1 Like

``````prob = ODEProblem(pipe_simp, x0, tspan, par)
sol = solve(prob; reltol=1e-4)

using Plots, LaTeXStrings
LW1 = 2.5
LW2 = 1.5
LS1 = :solid
LS2 = :dash
LC1 = :red
LC2 = :blue

plot(sol; idxs=p[1:N-1]/p_a, llw=LW2,ls=LS2, label=permutedims([L"\$p_{%\$i}(t)\$" for i in 1:N-1]))
plot!(sol; idxs=[pi0/p_a], lc=LC1, lw=LW1, label=L"p_{x=0}(t)")
plot!(sol; idxs=[p[N]/p_a], lc=LC2, lw=LW1, label=L"p_{x=L}(t)")
plot!(title="Gas Pipe, pressure distribution", ylabel=L"\$p\$ [atm]", xlabel=L"\$t\$ [s]")
``````

Result:

Quite nice! Thanks @baggepinnen

1 Like