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

**URL:** <https://discourse.julialang.org/t/mtk-vector-equations-get-error-message-on-axes-of-broadcast-not-known/94343>\
**Category:** Modelling & Simulations\
**Created:** [February 9, 2023, 12:50pm UTC](https://discourse.julialang.org/t/mtk-vector-equations-get-error-message-on-axes-of-broadcast-not-known/94343 "2023-02-09T12:50:16Z")\
**Posts on this page:** 10\
**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:** [February 9, 2023, 12:50pm UTC](https://discourse.julialang.org/t/mtk-vector-equations-get-error-message-on-axes-of-broadcast-not-known/94343/1 "2023-02-09T12:50:16Z")

</div>

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:

```julia
# 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[1] ~ 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:

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

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

```

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

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [February 9, 2023, 1:40pm UTC](https://discourse.julialang.org/t/mtk-vector-equations-get-error-message-on-axes-of-broadcast-not-known/94343/2 "2023-02-09T13:40:48Z")

</div>

> [@BLI](#):
>
> `Dt(m) ~ md_i - md_e`

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

```julia
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

```julia
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[1] ~ 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)] 

```

---

<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:** [February 9, 2023, 2:29pm UTC](https://discourse.julialang.org/t/mtk-vector-equations-get-error-message-on-axes-of-broadcast-not-known/94343/3 "2023-02-09T14:29:11Z")

</div>

I followed your advice on broadcasting.

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?

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [February 9, 2023, 2:37pm UTC](https://discourse.julialang.org/t/mtk-vector-equations-get-error-message-on-axes-of-broadcast-not-known/94343/4 "2023-02-09T14:37:02Z")

</div>

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.

---

<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:** [February 9, 2023, 2:43pm UTC](https://discourse.julialang.org/t/mtk-vector-equations-get-error-message-on-axes-of-broadcast-not-known/94343/5 "2023-02-09T14:43:33Z")

</div>

AHA! I see 🙂

---

<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:** [February 9, 2023, 3:01pm UTC](https://discourse.julialang.org/t/mtk-vector-equations-get-error-message-on-axes-of-broadcast-not-known/94343/6 "2023-02-09T15:01:46Z")

</div>

Hm…

```julia
pipe_simp = structural_simplify(pipe)

```

works fine.

```julia
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).

```julia
# 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:

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

```

fails miserably.

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

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [February 9, 2023, 3:32pm UTC](https://discourse.julialang.org/t/mtk-vector-equations-get-error-message-on-axes-of-broadcast-not-known/94343/7 "2023-02-09T15:32:21Z")

</div>

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

---

<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:** [February 9, 2023, 3:44pm UTC](https://discourse.julialang.org/t/mtk-vector-equations-get-error-message-on-axes-of-broadcast-not-known/94343/8 "2023-02-09T15:44:03Z")

</div>

OK – here is the complete code:

```julia
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[1] ~ 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.

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [February 9, 2023, 5:30pm UTC](https://discourse.julialang.org/t/mtk-vector-equations-get-error-message-on-axes-of-broadcast-not-known/94343/9 "2023-02-09T17:30:20Z")

</div>

This

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

```

creates a vector of vectors, simply

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

```

is what you want

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

---

<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:** [February 9, 2023, 9:15pm UTC](https://discourse.julialang.org/t/mtk-vector-equations-get-error-message-on-axes-of-broadcast-not-known/94343/10 "2023-02-09T21:15:23Z")

</div>

With your fix, I did:

```julia
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:

 ![image](https://global.discourse-cdn.com/julialang/original/3X/c/d/cd2dce736013e9824ee657d4555a625991b10ed4.png)

Quite nice! Thanks @baggepinnen