# ModelingToolkit: Some questions on usage & best practices
**URL:** https://discourse.julialang.org/t/modelingtoolkit-some-questions-on-usage-best-practices/77426
**Category:** Modelling & Simulations
**Tags:** modelingtoolkit
**Created:** [March 4, 2022, 9:30pm UTC](https://discourse.julialang.org/t/modelingtoolkit-some-questions-on-usage-best-practices/77426 "2022-03-04T21:30:37Z")
**Posts on this page:** 6
**Page:** 1
### Author:  [@RobertS](https://discourse.julialang.org/u/RobertS)
#### Post date: [March 4, 2022, 9:30pm UTC](https://discourse.julialang.org/t/modelingtoolkit-some-questions-on-usage-best-practices/77426/1 "2022-03-04T21:30:37Z")
Based on an earlier [discussion](https://discourse.julialang.org/t/rodas5-and-modelingtoolkit-much-slower-than-sundials-ida-on-large-systems/60712) 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](https://discourse.julialang.org/t/help-me-outperform-matlab-in-numerical-solution-of-semidiscretized-pde-system-as-much-as-possible/76777/20), 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:
```julia
using OrdinaryDiffEq, ModelingToolkit
using DASSL
using DASKR
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[1]) - q[1] * 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[1] + 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[2]
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[1]))`, but that won’t work in my case because of the scoping.
- Doing something like `plot(sol, vars=[sys.Tgas[1]])` does not work either, it gives the error: `ArgumentError: sys₊Tgas[1](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)[5]])`, 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[1]])`, 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?
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [March 6, 2022, 2:11am UTC](https://discourse.julialang.org/t/modelingtoolkit-some-questions-on-usage-best-practices/77426/2 "2022-03-06T02:11:48Z")
`plot(sol, vars=(sys.Tgas[1]))` should work on the latest versions. What version are you using?
---
### Author:  [@RobertS](https://discourse.julialang.org/u/RobertS)
#### Post date: [March 6, 2022, 6:54pm UTC](https://discourse.julialang.org/t/modelingtoolkit-some-questions-on-usage-best-practices/77426/3 "2022-03-06T18:54:26Z")
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:
```julia
plot(sol, vars=(sys.Tgas[1]))
# ArgumentError: sys₊Tgas[1](t) is either an observed nor a state variable.
# (should it read "_n_either... nor"?)
```
Further info:
```julia
states(sys)
#301-element Vector{Term{Real, Base.ImmutableDict{DataType, Any}}}:
# Tgas[1](t)
# Tgas[2](t)
# Tgas[3](t)
# ⋮
# cR[100](t)
# mflux(t)
```
```julia
sys.Tgas[1]
# sys₊Tgas[1](t)
```
```julia
Pkg.status()
# [165a45c3] DASKR v2.8.0
# [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
# [c3572dad] Sundials v4.9.2
# [0c5d862f] Symbolics v4.3.0
```
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [March 6, 2022, 9:15pm UTC](https://discourse.julialang.org/t/modelingtoolkit-some-questions-on-usage-best-practices/77426/4 "2022-03-06T21:15:02Z")
@YingboMa it should be released right?
---
### Author:  [@RobertS](https://discourse.julialang.org/u/RobertS)
#### Post date: [March 10, 2022, 10:02am UTC](https://discourse.julialang.org/t/modelingtoolkit-some-questions-on-usage-best-practices/77426/5 "2022-03-10T10:02:07Z")
I checked the release notes of the past view MTK releases but did not find anything that seems to relate to the namespacing of MTK variables.
---
### Author:  [@RobertS](https://discourse.julialang.org/u/RobertS)
#### Post date: [March 20, 2022, 8:49am UTC](https://discourse.julialang.org/t/modelingtoolkit-some-questions-on-usage-best-practices/77426/6 "2022-03-20T08:49:55Z")
Created an [issue](https://github.com/SciML/ModelingToolkit.jl/issues/1488).