I’m building fairly large models with ModelingToolkit.jl. While I think the library is very cool, I’m stumbling through how to set up large models effectively, both from a numerical perspective and from a developer perspective (making it easy to use the models I build). Looking for some guidance here. I keep reading articles about MTK being used for 10,000 variables and equations and I’m wondering “How?”–so clearly I’m missing something.
Requirements for my models:
@extendonly works for one layer, so I resorted to model composition to build nested components, i.e., define the physics in one model and build specializations on top of that model;- I only want to define core physics in one place, same goes for geometry definitions. I want to be able to take those core definitions and then use them in more specialized places.
- I need adjoints with respect to the parameters in the geometry (
Ain this case).
I have tried two approaches.
Option 1: Define geometry as a model
A MWE with two nested models:
using ModelingToolkit
using ModelingToolkit: t_nounits as t, D_nounits as D
using OrdinaryDiffEq
@component function ConduitGeometry(;
s_nodes = [0.0, 1.0], A_nodes = [1.0, 1.0], name = :ConduitGeometry)
N = length(s_nodes)
function trapz(sv, Av)
sum(0.5 * (Av[i] + Av[i + 1]) * (sv[i + 1] - sv[i]) for i in 1:(length(sv) - 1))
end
params = @parameters begin
s[1:N] = s_nodes
A[1:N] = A_nodes
A_start = A[1]
A_end = A[end]
H = s[end] - s[1]
V_total = trapz(s_nodes, A_nodes)
end
System(Equation[], t, Num[], params; name)
end
@mtkmodel HydroCore begin
@components begin
geometry = ConduitGeometry(; s_nodes, A_nodes) # base owns the one geometry
end
@parameters begin
k_out = 0.2
end
@variables begin
h(t)
q_in(t)
q_out(t)
end
@equations begin
q_out ~ k_out * h
D(h) ~ (q_in - q_out) / geometry.A_start
end
end
@mtkmodel TankUnit begin
@components begin
core = HydroCore(; geometry__s_nodes, geometry__A_nodes)
end
@parameters begin
Qin = 0.01
end
@equations begin
core.q_in ~ Qin
end
end
@mtkmodel PlantSystem begin
@components begin
unit = TankUnit(; Qin, core__geometry__s_nodes, core__geometry__A_nodes)
end
end
@mtkcompile sys = PlantSystem(;
unit.core.geometry.s_nodes = [0.0, 7.0],
unit.core.geometry.A_nodes = [0.1, 0.1],
unit.Qin = 8e-3
)
u0 = [sys.unit.core.h => 0.0]
tspan = (0.0, 50.0)
sol = solve(ODEProblem(sys, u0, tspan), Rodas4(); reltol = 1e-6, abstol = 1e-9)
Option 2: Represent geometry as a non-numeric parameter and pass through
Another alternative I thought about was making the geometry a non-numeric parameter in each model and pass them through (how you would do with a regular scalar), but I can’t get that to work at all. Note that I am using @parameters here because I want adjoints.
using ModelingToolkit
using ModelingToolkit: t_nounits as t, D_nounits as D
using OrdinaryDiffEq
function trapz(sv, Av)
sum(0.5 * (Av[i] + Av[i + 1]) * (sv[i + 1] - sv[i]) for i in 1:(length(sv) - 1))
end
struct Geometry
s_m::Vector{Float64}
A_m2::Vector{Float64}
height_m::Float64
volume_m3::Float64
function Geometry(s_m, A_m2)
height_m = s_m[end] - s_m[1]
volume_m3 = trapz(s_m, A_m2)
new(s_m, A_m2, height_m, volume_m3)
end
end
@mtkmodel BaseStage begin
@parameters begin
geometry::Geometry = Geometry([0.0, 7.0], [0.1, 0.1])
A0_m2 = 0.1
H_m = 7.0
kout = 0.2
end
@variables begin
h(t)
qin(t)
qout(t)
end
@equations begin
qout ~ kout * h
D(h) ~ (qin - qout) / A0_m2
end
end
@mtkmodel MidStage begin
@parameters begin
geometry::Geometry = Geometry([0.0, 7.0], [0.1, 0.1])
Qin = 0.01
end
@components begin
base = BaseStage(;
geometry = geometry,
)
end
@equations base.qin ~ Qin
end
@mtkmodel TopSystem begin
@parameters begin
geometry::Geometry = Geometry([0.0, 7.0], [0.1, 0.1])
end
@components begin
mid = MidStage(;
geometry,
)
end
end
@mtkcompile sys = TopSystem(
geometry = Geometry([0.0, 9.0], [0.1, 0.1]),
mid.Qin = 8e-3
)
u0 = [sys.mid.base.h => 0.0]
sol = solve(ODEProblem(sys, u0, (0.0, 50.0)), Rodas4(); reltol = 1e-6, abstol = 1e-9)
The error I get is
ERROR: LoadError: MethodError: Cannot `convert` an object of type Symbol to an object of type Expr
The function `convert` exists, but no method is defined for this combination of argument types.
Closest candidates are:
Expr(::Any...)
@ Core boot.jl:271
convert(::Type{T}, ::T) where T
@ Base Base.jl:126
convert(::Type{Expr}, ::Revise.RelocatableExpr)
@ Revise ~/.julia/packages/Revise/7QZSP/src/relocatable_exprs.jl:22
Stacktrace:
[1] push!(a::Vector{Expr}, item::Symbol)
@ Base ./array.jl:1260
[2] parse_equations!(exprs::Vector{Any}, eqs::Vector{Expr}, dict::Dict{Symbol, Any}, body::Expr)
@ ModelingToolkit ~/.julia/packages/ModelingToolkit/vyvd9/src/systems/model_parsing.jl:1107
[3] parse_model!(exprs::Vector{…}, comps::Vector{…}, ext::Vector{…}, eqs::Vector{…}, icon::Base.RefValue{…}, vs::Vector{…}, ps::Vector{…}, sps::Vector{…}, c_evts::Vector{…}, d_evts::Vector{…}, cons::Vector{…}, costs::Vector{…}, dict::Dict{…}, mod::Module, arg::Expr, kwargs::OrderedCollections.OrderedSet{…}, where_types::Vector{…})
@ ModelingToolkit ~/.julia/packages/ModelingToolkit/vyvd9/src/systems/model_parsing.jl:644
[4] _model_macro(mod::Module, fullname::Symbol, expr::Expr, isconnector::Bool)
@ ModelingToolkit ~/.julia/packages/ModelingToolkit/vyvd9/src/systems/model_parsing.jl:82
[5] var"@mtkmodel"(__source__::LineNumberNode, __module__::Module, fullname::Union{Expr, Symbol}, body::Any)
@ ModelingToolkit ~/.julia/packages/ModelingToolkit/vyvd9/src/systems/model_parsing.jl:33
[6] include(fname::String)
@ Main ./sysimg.jl:38
[7] top-level scope
@ REPL[1]:1
My questions specifically are:
- Option 1 “works” (compiles and runs), but having to pass those keywords up through each subsequent model feels hacky. Is that really the best way of doing this (I’m just imaging this with a 10-level nested model for example)?
- Option 2 does not work. Is this a better option if we could get it to work? I like the syntax and set up of it more. Also not sure how that would work with adjoints.