# ModelingToolkit: Some questions on usage & best practices

Based on an earlier discussion here on discourse formulated a physical system in Modeling toolkit. The system consists of a discretized 1D PDE, ODEs and an algebraic equation. The physical nature of the system is further described in my earlier posting, if anyone is interested.

I am still in my very early days of using Julia and ModelingToolkit, so I have a few hopefully simple questions about usage of MTK.

My code is as follows:

``````using OrdinaryDiffEq, ModelingToolkit
using DASSL
using DifferentialEquations
using BenchmarkTools
using Plots
using Sundials

function build_system(n)

@parameters t # or would @variables t be better!?
@variables Tgas[1:n](t) Tsol[1:n](t) cR[1:n](t) mflux(t)
@parameters params[1:3]
D = Differential(t)

dp, Tin, h = params

dx = h / n

############################################################
# Build equations ##########################################
############################################################

eqns = Array{Equation}(undef, 3 * n + 1)
count = 1

# specification as function or intermediate variabel does not matter. I could also make this an
# actual @variable, then I could retrieve it later.
rR = 5.5e6 .* exp.(.-5680 ./(273 .+ Tsol)) .* (1 .- exp.(.-cR./50))
q = 1e4 * (Tgas - Tsol)

# gas energy balance

eqns[count] = 0 ~  mflux * 1000 * (Tin - Tgas) - q * dx
count += 1
for i in 2:n
eqns[count] = 0 ~ mflux * 1000 * (Tgas[i-1] - Tgas[i]) - q[i] * dx
count += 1
end

# solids energy balance
for i in 1:n
eqns[count] = D(Tsol[i]) ~ q[i] / (1000 * 2000)
count += 1
end

# gas momentum balance
eqns[count] = 0 ~ - mflux + 5e-3 * sqrt(dp)
count += 1

# reactant concentration equations
for i in 1:n
eqns[count] = D(cR[i]) ~ -rR[i]
count += 1
end

@named sys = ODESystem(eqns, t, vcat(Tgas, Tsol, cR, [mflux]), params)

return sys

end

function build_problem(sys; tspan = (0.0, 1.5), param_values = [100e2, 200, 0.5])

n = Int((length(sys.states) - 1)/3)
u0 = ones(length(sys.states))
u0[1:n] .= param_values
u0[n+1 : 2*n] .= 20
u0[2*n+1 : 3*n] .= 3000
u0[3*n+1] = 0.5

u0map = [
sys.Tgas => 100.0,
sys.Tsol => 20.0,
sys.cR => 3000.0,
sys.mflux => 0.5
]

prob = ODEProblem(sys, u0, tspan, param_values, jac=true, sparse=true)
return prob

end

dp = 100e2
Tin = 200
h = 0.5
parameters = [dp, Tin, h]

N = 100
tspan = (0.0, 1500.0)

sys = build_system(N)
prob = build_problem(sys, tspan=tspan, param_values = parameters)
sol = solve(prob, QBDF())

``````

As you see, my code is structured such that all the `@variables` are defined in the local scope of function `build_system`. Since that is the case, how would I later, in another function or in global scope access parts of the solutions, like for example in the plotting recipe?

• Tutorials typically show something like `plot(sol, vars=(Tgas))`, but that won’t work in my case because of the scoping.

• Doing something like `plot(sol, vars=[sys.Tgas])` does not work either, it gives the error: `ArgumentError: sys₊Tgas(t) is either an observed nor a state variable.` (BTW, there seems to be a typo in the error message, I guess it should say “neither” instead of “either”)

• What works is: `plot(sol, vars=[states(sys)])`, but that does not allow me to specify the variable name and is not best practice, I guess

• in some thread I also found a hint on using the (undocumented?) `@nonamespace` macro, which works: `plot(sol, vars=[@nonamespace sys.Tgas])`, but is that really the intended way?

Is there a good way to do it? A question very related to this would be how to define `u0map` in the function `build_problem`. Using the definition in my code listing above, building the `ODEProblem` throws an error, which I guess is again due to namespacing of the MTK variables.

So, what are the best practices here? Is it a non-optimal idea to define the variables, build the ODE problem and plotting the solution in different functions / scopes in the first place?

`plot(sol, vars=(sys.Tgas))` should work on the latest versions. What version are you using?

Thanks for the hint. I updated MTK to 8.5.1 (sorry, forgot to check which version I had before updating).
But still, I get the following error:

``````plot(sol, vars=(sys.Tgas))
# ArgumentError: sys₊Tgas(t) is either an observed nor a state variable.
# (should it read "_n_either... nor"?)
``````

Further info:

``````states(sys)
#301-element Vector{Term{Real, Base.ImmutableDict{DataType, Any}}}:
# Tgas(t)
# Tgas(t)
# Tgas(t)
# ⋮
# cR(t)
# mflux(t)
``````
``````sys.Tgas
# sys₊Tgas(t)
``````
``````Pkg.status()
#  [e993076c] DASSL v2.6.0
#  [0c46a032] DifferentialEquations v7.1.0
#  [28b8d3ca] GR v0.64.0
#  [33e6dc65] MKL v0.5.0
#  [961ee093] ModelingToolkit v8.5.1
#  [1dea7af3] OrdinaryDiffEq v6.7.0
#  [91a5bcdd] Plots v1.26.0
#  [d236fae5] PreallocationTools v0.2.4