I’m seeing a whole lot of commits and minor version increases the past week or two. It seems like a lot of it is internal reworking / refactoring. Just curious to know what it’s all about
Lots of work has gone into improving compilation and runtime of mtkcompile and ODEproblem.
Great! Looking forward to trying it out
I am seeing these improvements, which are really great! Now, for my MTK model, the JIT-compilation of the RHS function has become the bottleneck (takes around 15 minutes while the structural simplification takes just a minute). Is this something that will be looked at in the future?
Everything between “model change” and “plot in hand” is being worked at 
15 minutes is very high, do you have an example that can be scaled down in size that can be used as benchmark target?
Here is a minimal, scalable reproducer for the slow model-build of a large MTK system:
Builds a chain of N point masses. Each particle carries a chain of M per-particle observed intermediates (a shared body-frame rotation R + smooth_norm, each referencing the previous one, like a real computation pipeline), plus drag, and each of the N-1 segments a nonlinear spring. M sets the observed/state density. M≈58 gives ~31 observed per state, matching my
real model.
All times are full first-run wall-clock (seconds) per phase:
simplify : mtkcompile
odeprob : ODEProblem, build_initializeprob=false (build_function + RHS spec.)
rhs 1st : first prob.f call (RHS is largely specialized during odeprob)
getter : build + first call of a getter over all observed
The rhs/getter costs are dominated by single-threaded first-call JIT (inference+LLVM).
Sidenote: a parameter-heavy version crashes mtkcompile (MethodError: isless(::Type{...}) in the tearing sort), a separate MTK bug. Literals sidestep it.
using ModelingToolkit
using ModelingToolkit: t_nounits as t, D_nounits as D
using LinearAlgebra
smooth_norm(v) = sqrt(v[1]^2 + v[2]^2 + v[3]^2 + 1e-6)
calc_rho(height) = 1.225 * exp(-height / 8400.0)
wind_factor(height) = (max(height, 1.0) / 6.0)^0.14
"""
build_chain(N; M)
Return an `mtkcompile`d `System`: `N` dynamic point masses, each with `M` chained
intermediate observed 3-vectors reusing a shared rotation `R`. Particle 1 pinned
by a stiff restoring force.
"""
function build_chain(N; M)
@variables (pos(t))[1:3, 1:N] (vel(t))[1:3, 1:N] (acc(t))[1:3, 1:N]
@variables (drag(t))[1:3, 1:N] (va(t))[1:3, 1:N]
@variables (aux(t))[1:3, 1:(M * N)] (spring(t))[1:3, 1:(N - 1)] (ang(t))[1:3]
auxidx(k, i) = (i - 1) * M + k
ca, cb, cc = cos(ang[1]), cos(ang[2]), cos(ang[3])
sa, sb, sc = sin(ang[1]), sin(ang[2]), sin(ang[3])
Rt = permutedims([ cb*cc (sa*sb*cc - ca*sc) (ca*sb*cc + sa*sc);
cb*sc (sa*sb*sc + ca*cc) (ca*sb*sc - sa*cc);
-sb sa*cb ca*cb ])
eqs = Equation[]
for j in 1:(N - 1)
d = pos[:, j + 1] - pos[:, j]
len = smooth_norm(d)
push!(eqs, spring[:, j] ~ 5.0e3 * (len - 1.0) * (d / len))
end
for i in 1:N
push!(eqs, va[:, i] ~ [12.0, 0.0, 0.0] - vel[:, i])
height = max(0.0, pos[3, i])
push!(eqs, drag[:, i] ~ 0.5 * calc_rho(height) * 1.2 *
wind_factor(height) * smooth_norm(va[:, i]) * 0.05 * va[:, i])
# Each aux is a DIRECT function of the states (no chaining), so mtkcompile
# eliminates it to `observed` instead of tearing it into the state vector.
for k in 1:M
push!(eqs, aux[:, auxidx(k, i)] ~
Rt * vel[:, i] + (0.001 * k) * smooth_norm(pos[:, i]) * pos[:, i])
end
F = drag[:, i] - 5.0 * vel[:, i] + [0.0, 0.0, -9.81]
for k in 1:M
F = F + 0.0001 * aux[:, auxidx(k, i)]
end
i > 1 && (F = F - spring[:, i - 1])
i < N && (F = F + spring[:, i])
i == 1 && (F = F - 1.0e5 * pos[:, i])
push!(eqs, acc[:, i] ~ F)
push!(eqs, D(pos[:, i]) ~ vel[:, i])
push!(eqs, D(vel[:, i]) ~ acc[:, i])
end
push!(eqs, D(ang) ~ 0.05 * ang)
@named sys = System(eqs, t)
return mtkcompile(sys)
end
function bench(N; M)
GC.gc()
t_simplify = @elapsed sys = build_chain(N; M)
sts = ModelingToolkit.unknowns(sys)
u0map = [s => 0.5 + 0.01 * i for (i, s) in enumerate(sts)]
t_odeprob = @elapsed prob =
ODEProblem(sys, u0map, (0.0, 1.0); build_initializeprob=false)
du = similar(prob.u0)
t_rhs = @elapsed prob.f(du, prob.u0, prob.p, 0.0)
all_obs = ModelingToolkit.observed(sys)
obs_syms = [eq.lhs for eq in all_obs]
t_getter = @elapsed begin
getter = ModelingToolkit.getu(sys, obs_syms)
getter(prob)
end
n_obs = length(all_obs)
n_states = length(ModelingToolkit.unknowns(sys))
total = t_simplify + t_odeprob + t_rhs + t_getter
return (; N, M, n_states, n_obs, ratio=round(n_obs / n_states; digits=1),
t_simplify, t_odeprob, t_rhs, t_getter, total)
end
# Sweep intermediate density M at fixed N: isolates the effect of observed count
# (does reducing intermediate @variables speed up the build?).
N_FIXED = 40
MS = [2, 10, 20, 40]
println("warming up…"); flush(stdout)
bench(8; M=2)
rows = NamedTuple[]
for M in MS
println("--- N=$N_FIXED, M=$M ---"); flush(stdout)
r = bench(N_FIXED; M)
push!(rows, r)
println(" states=$(r.n_states) observed=$(r.n_obs) ratio=$(r.ratio) ",
"total=$(round(r.total; digits=1)) s")
flush(stdout)
end
println("\n\n===== FIRST-RUN WALL TIME vs observed density (N=$N_FIXED) =====")
println(rpad("M", 5), rpad("observed", 10), rpad("ratio", 7),
rpad("simplify", 10), rpad("odeprob", 9), rpad("rhs1st", 8),
rpad("getter", 8), "TOTAL")
for r in rows
println(rpad(r.M, 5), rpad(r.n_obs, 10), rpad(r.ratio, 7),
rpad(round(r.t_simplify; digits=1), 10),
rpad(round(r.t_odeprob; digits=1), 9),
rpad(round(r.t_rhs; digits=1), 8),
rpad(round(r.t_getter; digits=1), 8),
round(r.total; digits=1))
end
println("================================================================")
flush(stdout)
Output:
warming up…
--- N=40, M=2 ---
states=243 observed=722 ratio=3.0 total=4.5 s
--- N=40, M=10 ---
states=243 observed=1682 ratio=6.9 total=18.1 s
--- N=40, M=20 ---
states=243 observed=2882 ratio=11.9 total=48.9 s
--- N=40, M=40 ---
states=243 observed=5282 ratio=21.7 total=165.7 s
===== FIRST-RUN WALL TIME vs observed density (N=40) =====
M observed ratio simplify odeprob rhs1st getter TOTAL
2 722 3.0 0.3 0.7 0.8 2.7 4.5
10 1682 6.9 0.9 6.5 1.5 9.2 18.1
20 2882 11.9 1.5 23.5 3.1 20.8 48.9
40 5282 21.7 3.3 95.7 5.7 61.0 165.7
================================================================
@dhairyagandhi96 any chance that your DCP passes identifies the repeated structure in this model?
Very likely yes! The difference is that I need a better system to identify matmul-expansions without needing to walk sub systems. This model doesn’t have algebraic equations, so our first “seeding” of matrix discovery would need a slightly different approach.
The affine transforms already should work, and we get some nice clusters emerging as well!
julia> SC.get_candidate_expr(ir_ref, clusters_sorted[2]) # the smooth_norm function
sqrt(1.0e-6 + (-(pos(t))[1, 2] + (pos(t))[1, 3])^2 + (-(pos(t))[2, 2] + (pos(t))[2, 3])^2 + (-(pos(t))[3, 2] + (pos(t))[3, 3])^2)
julia> SC.get_candidate_expr(ir_ref, clusters_sorted[10]) # the expanded wind_factor
(max(max(0.0, (pos(t))[3, 3]), 1.0) / 6.0)^0.14
julia> SC.get_candidate_expr(ir_ref, clusters_sorted[1]) # aux expanded
0.002(pos(t))[2, 3]*sqrt(1.0e-6 + (pos(t))[1, 3]^2 + (pos(t))[2, 3]^2 + (pos(t))[3, 3]^2) + cos((ang(t))[2])*sin((ang(t))[1])*(vel(t))[3, 3] + (-cos((ang(t))[1])*sin((ang(t))[3]) + cos((ang(t))[3])*sin((ang(t))[2])*sin((ang(t))[1]))*(vel(t))[1, 3] + (cos((ang(t))[1])*cos((ang(t))[3]) + sin((ang(t))[2])*sin((ang(t))[1])*sin((ang(t))[3]))*(vel(t))[2, 3]
@baggepinnen So the idea here would be to find the most common operations like smooth_norm and extract it into shared methods, such that the RHS function will become smaller?
That’s part of it. We have a full structural recovery pass which handles recurrences, loops etc and allows scaling. The clusters were an example of some artifacts we produce along the way.
Although perhaps a distraction to your primary Julia-related question, may I ask you to describe the physical/engineering context here? Purely for my Julia-unrelated curiosity.
Take a look at SymbolicAWEModels.jl which allows you to generate models of kite power systems (airborne wind energy). An example of an implemented kite model of an open source kite: V3Kite.jl. This kite model is still relatively simple, it doesn’t have a lot of bridle points, and has around 100-200 differential states. But for a private company I made a model of a kite with a lot more bridle points, which is why I am hitting these performance issues now.
If I understand it correctly, the large JIT compile times are because everything is inlined. Could we have an option to do less inlining and create shared methods to reduce JIT compile time?
This was a major problem for us with Neuroblox. For us though, we were interested in systems which are so big that the solutions that @dhairyagandhi96 is talking about aren’t really applicable, because even building the symbolic representation in the first place was too costly, so it would be too late to try and ‘rediscover’ repeated structures and factor them out.
Rather, what we did was create something where the loop/graph structure is in the problem directly, and doesn’t need to be “rediscovered” post-hoc. We made a graph-compiler for this in GraphDynamics.jl. I gave a Juliacon talk about it last year here:
We’re still actively developing it, but for now our current versions are closed source, and the open-source release is a bit outdated (uses old versions of SciML and OrderedCollections).
Regardless though, I had Claude turn your MTK problem into an equivalent GraphDynamics problem and then solve it (accuracy of the final results agree to within 1e-5):
❯ julia --project=. chain_gd_demo.jl
warming up…
--- N=40, M=2 ---
states=240 observed=80 ratio=0.3 total=0.009 s
--- N=40, M=10 ---
states=240 observed=400 ratio=1.7 total=1.055 s
--- N=40, M=20 ---
states=240 observed=800 ratio=3.3 total=1.222 s
--- N=40, M=40 ---
states=240 observed=1600 ratio=6.7 total=1.1 s
--- N=40, M=100 ---
states=240 observed=4000 ratio=16.7 total=1.172 s
--- N=40, M=1000 ---
states=240 observed=40000 ratio=166.7 total=13.975 s
===== GraphDynamics FIRST-RUN WALL TIME vs observed density (N=40) =====
M observed ratio build odeprob rhs1st getter TOTAL
2 80 0.3 0.008 0.001 0.0 0.0 0.009
10 400 1.7 0.243 0.489 0.088 0.236 1.055
20 800 3.3 0.254 0.589 0.096 0.282 1.222
40 1600 6.7 0.278 0.473 0.088 0.26 1.1
100 4000 16.7 0.244 0.473 0.086 0.369 1.172
1000 40000 166.7 0.24 0.496 0.091 13.148 13.975
================================================================
For fun, I added M = 100 and M = 1000 to the benchmark. Main takeaway is that GraphDynamics only compiles one method per set of node-connection-node combinations, and the compilation is independant of the size of your problem. That’s why the build time doesn’t scale with your problem size, and we’re able to go up to M = 1000 without any issue.
The GraphDynamics version with M=1000 is over 10x faster than the MTK M = 40 example.
Here’s the benchmark script (when you install GraphDynamics make sure you get v0.7 to make this work, this might require downgrading some other packages):
using GraphDynamics
using OrdinaryDiffEqTsit5
using SymbolicIndexingInterface: getu
# ------------------------------------------------------------------
# Same physics as the MTK reproducer:
# - N point masses in a chain, nonlinear springs between neighbours
# - each particle carries a drag force
# - each particle carries M chained "observed" auxiliary 3-vectors,
# reusing a shared body-frame rotation R, which also (weakly) feed
# back into the particle's own force balance
# - particle 1 is pinned by a stiff restoring force
#
# The key structural difference vs. the MTK version: M is a *runtime*
# loop bound baked into ONE compiled method per M value, not M separate
# symbolic equations that have to be tearing-sorted and code-generated
# individually.
# ------------------------------------------------------------------
function rotation_matrix(ang)
a, b, c = ang
ca, cb, cc = cos(a), cos(b), cos(c)
sa, sb, sc = sin(a), sin(b), sin(c)
# rows of R (already the transpose, matching `Rt` in the MTK version)
(
(cb*cc, sa*sb*cc - ca*sc, ca*sb*cc + sa*sc),
(cb*sc, sa*sb*sc + ca*cc, ca*sb*sc - sa*cc),
(-sb, sa*cb, ca*cb),
)
end
matvec3((r1, r2, r3), (x, y, z)) =
(r1[1]*x + r1[2]*y + r1[3]*z,
r2[1]*x + r2[2]*y + r2[3]*z,
r3[1]*x + r3[2]*y + r3[3]*z)
smooth_norm3(x, y, z) = sqrt(x^2 + y^2 + z^2 + 1e-6)
calc_rho(height) = 1.225 * exp(-height / 8400.0)
wind_factor(height) = (max(height, 1.0) / 6.0)^0.14
# k-th chained auxiliary 3-vector for a particle at (px,py,pz),(vx,vy,vz)
@inline function aux_term(k, px, py, pz, vx, vy, vz)
RT = rotation_matrix((0.1, 0.2, 0.3))
rvx, rvy, rvz = matvec3(RT, (vx, vy, vz))
n = smooth_norm3(px, py, pz)
c = 0.001 * k * n
(rvx + c*px, rvy + c*py, rvz + c*pz)
end
# ------------------------------------------------------------------
# Particle subsystem. `M` (number of chained aux observables) is a type
# parameter so that `computed_properties` can expose `aux1 .. auxM` as
# genuine named "observed" quantities queryable via the standard
# SymbolicIndexingInterface machinery (`sol[:p3₊aux17]`, `getu`, etc.),
# exactly mirroring what `observed(sys)` gives you in the MTK version.
# ------------------------------------------------------------------
struct Particle{M} end
@kwdef struct ParticleSpec{M}
name::Symbol
m::Float64 = 1.0
k_anchor::Float64 = 0.0
pos_init::NTuple{3,Float64}
vel_init::NTuple{3,Float64} = (0.0, 0.0, 0.0)
end
function GraphDynamics.to_subsystem(p::ParticleSpec{M}) where {M}
(; name, m, k_anchor, pos_init, vel_init) = p
(px, py, pz) = pos_init
(vx, vy, vz) = vel_init
states = SubsystemStates{Particle{M}}(; px, py, pz, vx, vy, vz)
params = SubsystemParams{Particle{M}}(; m, k_anchor)
Subsystem(states, params)
end
GraphDynamics.initialize_input(::Subsystem{Particle{M}}) where {M} = (; Fx=0.0, Fy=0.0, Fz=0.0)
function GraphDynamics.subsystem_differential(sys::Subsystem{Particle{M}}, input, t) where {M}
(; px, py, pz, vx, vy, vz, m, k_anchor) = sys
(; Fx, Fy, Fz) = input
# drag, purely a function of this particle's own state
wx, wy, wz = 12.0 - vx, -vy, -vz
height = max(0.0, pz)
dragfac = 0.5 * calc_rho(height) * 1.2 * wind_factor(height) * smooth_norm3(wx, wy, wz) * 0.05
Fx += dragfac * wx
Fy += dragfac * wy
Fz += dragfac * wz
# M chained aux observables feeding weakly back into the force,
# same as `F = F + 0.0001 * aux[:, k]` in the MTK version.
ax, ay, az = 0.0, 0.0, 0.0
for k in 1:M
axk, ayk, azk = aux_term(k, px, py, pz, vx, vy, vz)
ax += axk; ay += ayk; az += azk
end
Fx += 0.0001 * ax
Fy += 0.0001 * ay
Fz += 0.0001 * az
# gravity + linear damping + optional pin-to-origin anchor
Fx += -5.0 * vx - k_anchor * px
Fy += -5.0 * vy - k_anchor * py
Fz += -5.0 * vz - 9.81 - k_anchor * pz
SubsystemStates{Particle{M}}(; px=vx, py=vy, pz=vz, vx=Fx/m, vy=Fy/m, vz=Fz/m)
end
function GraphDynamics.computed_properties(::Type{Particle{M}}) where {M}
names = ntuple(k -> Symbol(:aux, k), M)
funcs = ntuple(M) do k
sys -> aux_term(k, sys.px, sys.py, sys.pz, sys.vx, sys.vy, sys.vz)
end
NamedTuple{names}(funcs)
end
# ------------------------------------------------------------------
# Nonlinear spring connection between chain neighbours
# ------------------------------------------------------------------
struct Spring3D <: ConnectionRule
k::Float64
L0::Float64
end
Base.zero(::Type{Spring3D}) = Spring3D(0.0, 0.0)
function ((; k, L0)::Spring3D)(src::Subsystem, dst::Subsystem, t)
dx = src.px - dst.px
dy = src.py - dst.py
dz = src.pz - dst.pz
len = smooth_norm3(dx, dy, dz)
f = k * (len - L0) / len
(; Fx = f*dx, Fy = f*dy, Fz = f*dz)
end
# ------------------------------------------------------------------
# Build & bench, mirroring the MTK script's phases:
# build : constructing the GraphSystem (~ mtkcompile)
# odeprob : ODEProblem construction
# rhs 1st : first prob.f call
# getter : build + first call of a getter over all `auxK` observables
# ------------------------------------------------------------------
function build_chain(N; M)
g = GraphSystem()
specs = [ParticleSpec{M}(; name=Symbol(:p, i), m=1.0,
k_anchor = i == 1 ? 1.0e5 : 0.0,
pos_init=(Float64(i), 0.0, 0.0))
for i in 1:N]
for s in specs
add_node!(g, s)
end
for i in 1:N-1
add_connection!(g, specs[i] => specs[i+1]; conn=Spring3D(5.0e3, 1.0))
add_connection!(g, specs[i+1] => specs[i]; conn=Spring3D(5.0e3, 1.0))
end
g
end
function bench(N; M)
GC.gc()
t_build = @elapsed g = build_chain(N; M)
u0 = Pair[]
for i in 1:N
push!(u0, Symbol(:p, i, :₊, :px) => 0.5 + 0.01*i)
push!(u0, Symbol(:p, i, :₊, :py) => 0.5 + 0.01*i)
push!(u0, Symbol(:p, i, :₊, :pz) => 0.5 + 0.01*i)
end
t_odeprob = @elapsed prob = ODEProblem(g, u0, (0.0, 1.0))
du = similar(prob.u0)
t_rhs = @elapsed prob.f(du, prob.u0, prob.p, 0.0)
obs_syms = Symbol[Symbol(:p, i, :₊, :aux, k) for i in 1:N for k in 1:M]
t_getter = @elapsed begin
getter = getu(prob, obs_syms)
getter(prob)
end
n_states = length(prob.u0)
n_obs = length(obs_syms)
total = t_build + t_odeprob + t_rhs + t_getter
return (; N, M, n_states, n_obs, ratio=round(n_obs / n_states; digits=1),
t_build, t_odeprob, t_rhs, t_getter, total)
end
function (@main)(args)
N_FIXED = 40
MS = [2, 10, 20, 40, 100, 1000]
println("warming up…"); flush(stdout)
bench(8; M=2)
rows = NamedTuple[]
for M in MS
println("--- N=$N_FIXED, M=$M ---"); flush(stdout)
r = bench(N_FIXED; M)
push!(rows, r)
println(" states=$(r.n_states) observed=$(r.n_obs) ratio=$(r.ratio) ",
"total=$(round(r.total; digits=3)) s")
flush(stdout)
end
println("\n\n===== GraphDynamics FIRST-RUN WALL TIME vs observed density (N=$N_FIXED) =====")
println(rpad("M", 5), rpad("observed", 10), rpad("ratio", 7),
rpad("build", 10), rpad("odeprob", 9), rpad("rhs1st", 8),
rpad("getter", 8), "TOTAL")
for r in rows
println(rpad(r.M, 5), rpad(r.n_obs, 10), rpad(r.ratio, 7),
rpad(round(r.t_build; digits=3), 10),
rpad(round(r.t_odeprob; digits=3), 9),
rpad(round(r.t_rhs; digits=3), 8),
rpad(round(r.t_getter; digits=3), 8),
round(r.total; digits=3))
end
println("================================================================")
end
Sidenote, but heterogeneous systems are very important to us, so extra care was taken in the implementation of GraphDynamics to make sure you can still maintain optimal performance while having many different node types with different dynamics, and many different connection types with different rules. so the numbers you see here aren’t a special case of the fact that you only had one type of particle, and one type of spring.
That’s really interesting, thanks for sharing! And yes, i agree that in cases where the cost of model simplification/ symbolic code generation is high, recovery after the fact isn’t optimal. I’ve been working on methods of walking raw (unsimplified) models and constructing the graph directly on shared structural information from the sub-systems too!
It’s very important to note the scope difference though. GraphDynamics is for causally connected graph dynamical systems. Nodes in GraphDynamics can be more complex objects, i.e. you can use acausal modeling or something like Catalyst to construct a flat set of equations for the node, but the overlaid connector system of GraphDynamics is a very small subset of what MTK supports. Additionally, it doesn’t handle most of the sophisticated DAE machinery and the pieces required for numerically difficult problems, such as the initialization systems, homotopy codegen, index reduction, tearing, customized nonlinear solvers, etc., along with not supporting the synchronous clocking elements and the other things in the long tail of features industrial modelers tend to use. And of course, optimization problems, dynamic optimization problem to BVP codegen, etc. again the long tail of “not just an ODE things”.
Because of this, GraphDynamics won’t handle “most” DAEs effectively and should not be recommended for most of the workflows that people use MTK for.
But, that’s not to downplay that GraphDynamics is a very useful tool. If you have a set of ODEs which scale to large systems by repeating the same set in a causally connected graph of interactions, then GraphDynamics will scale much better. There are many applications which fall into this exact domain: neuroscience was what it was made for, but spatial epidemic models, ecological models, etc. are all other areas I have seen match similar constructions. For those applications, GraphDynamics will always be a much better front end for scaling by narrowing the design scope and effectively eliminating most of the compiler work that MTK and going straight to the optimal codegen for this form. But it should be careful because not everything falls into this domain just because it has a graph: power system dynamics is a good example where the DAEs need to have alternative codegen elements like limiters and homotopies which are not handled in the GraphDynamics system (nor should it probably expand in that way, because then it would need to start doing analyses that fundamentally impact what it is good at!)
So that’s just a word of caution that should come along with it: probably like 10-20% of what I see done in MTK is a good fit for GraphDynamics, and for that subset it will be a better tool for scaling. But that caveat of what it’s for should go along with the statement that it’s better at scaling.
A lot of caution should be taken with the interpretation of these results. One example of a GraphDynamics system is a partial differential equation system discretization: the semidiscretization is a graph with the ODEs of a given location and the graph given by the strength of the connection, which is just the operator discretization weights given by FDM/FVM/FEM/… So for an isolated part of a multibody system, it is a PDE semi-discretization and GraphDynamics will scale to making that large. GraphDynamics will do really well for these PDE semi-discretizations which do not have algebraic equations (i.e. so not Navier-Stokes) and have simple boundary conditions (linear boundary conditions) that can be restated explicitly rather than implicitly. So one PDE like a reaction-diffusion equation with 100,000 spatial points, GD is a nice way to represent it.
However, the complexity of multibody systems is generally not from the single components but from the interactions between all of the components. Those acausal components mean you effectively have lots of little PDEs, some with even like 4 or 5 spatial points. The complexity is that you have lots of them and worse, you have a state selection problem, i.e. it’s a DAE with many different possible choices for what variables should be the ones to solve for, and the choice will often be different between the different components. For this kind of problem, solving the linear PDEs inline (i.e. not through the solver but reducing it to a linear subsystem), and finding the SCC decomposition to solve connected pieces sequentially, can be required in order to get something that’s a stable solve. That’s not something GraphDynamics can do.
But also, nor should it handle that! Its fast compilation is fast simply because it doesn’t do state selection, index reduction, SCC restructuring, inline linear solves, etc. which are some of the difficult parts of the MTK pipeline that lead to requiring scalar analyses.
So point is, right tool for the right job. GraphDynamics is a really good tool for where it applies, but we should really use caution when talking about it say as a thing that can handle a multibody system: it really does a single body system 
Yes of course. The point is that both the problems that I tend to look at, and the one that @Bart_van_de_Lint presented don’t use any of that DAE machinery other than just (ab)using it to write very simple transmission of forces or voltages between nodes. For those cases, all of the (very impressive and useful) stuff that MTK does is superfluous.
I don’t mean to knock MTK here, it’s an incredible tool and I should have provided more caveats, but I thought it was clear from the context of Bart’s problem that he was also in a situation like the one Neuroblox was in where he was more using MTK as a way of composing many structurally identical building blocks together, rather than using it for it’s powerful simplification and transformation tools. For that, MTK is just not (currently) the right tool.
Would there maybe be a way to get the best of both worlds and give hints or directives to MTK that certain parts of the problem can be handled differently because they are just repeated structures, and maybe skip some work there?
Like the example of the transmission line consisting of capacitors and inductors (I saw it somewhere).
In my case it would be things like fluid pipes, where you repeat flow resistances and volumes connected in a chain, or more general square grids of volumes and fluid resistances.
I have had good success so far with MTK-- just the fact you have access to much, much better nonlinear and ODE solvers than are available with Python has made so that problems that previously would just not converge now do so with no problem.
The composability of the components also is very handy since now I don’t have to hand craft so much.
But of course I want to do more. I had seen GraphDynamics and thought maybe it would be good for those sort of many node repeated volume and flow resistance networks, but decided that both the flexibility and generality of MTK would be better for taking on the immediate problems.
Yep, that would be the dream.
Not sure how relevant this is, but NetworkDynamics.jl also exists, which is a backend for PowerDynamics.jl (power system dynamics over a graph)