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 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?)
@nonamespacemacro, 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?