Slow access to solution via symbols

[spinnau]

Access to the solution object via symbols can be quite slow for large solutions with many time steps. For demonstration, I have modified the example RC circuit from the tutorials with a sinusoidal input voltage and simulated it for a long period of time to get a larger solution object (please see the example below).

Access to the solution for one variable takes approx. 1.7 seconds for this example. Maybe this issue is related to https://github.com/SciML/ModelingToolkit.jl/issues/1102. But there's not much difference here, if the accessed variable is a state or an observed variable:

julia> @benchmark sol[capacitor.v] # state variable

BenchmarkTools.Trial: 3 samples with 1 evaluation.
 Range (min … max):  1.665 s …    1.855 s  ┊ GC (min … max): 0.00% … 10.82%
 Time  (median):     1.698 s               ┊ GC (median):    0.00%
 Time  (mean ± σ):   1.740 s ± 101.193 ms  ┊ GC (mean ± σ):  3.85% ±  6.25%

  █         █                                              █  
  █▁▁▁▁▁▁▁▁▁█▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁█ ▁
  1.67 s         Histogram: frequency by time         1.85 s <

 Memory estimate: 760.17 MiB, allocs estimate: 21542941.

julia> @benchmark sol[resistor.v] # observed variable

BenchmarkTools.Trial: 3 samples with 1 evaluation.
 Range (min … max):  1.767 s …   1.902 s  ┊ GC (min … max): 0.00% … 11.14%
 Time  (median):     1.796 s              ┊ GC (median):    0.00%
 Time  (mean ± σ):   1.822 s ± 71.156 ms  ┊ GC (mean ± σ):  3.88% ±  6.43%

  █           █                                           █  
  █▁▁▁▁▁▁▁▁▁▁▁█▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁█ ▁
  1.77 s         Histogram: frequency by time         1.9 s <

 Memory estimate: 760.17 MiB, allocs estimate: 21542941.

It is much faster if the solution array will be directly accessed:

julia> @benchmark reduce(hcat, sol.u)[1, :]
BenchmarkTools.Trial: 34 samples with 1 evaluation.
 Range (min … max):  135.073 ms … 152.920 ms  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     148.863 ms               ┊ GC (median):    0.00%
 Time  (mean ± σ):   148.531 ms ±   3.524 ms  ┊ GC (mean ± σ):  0.20% ± 1.14%

                                              ▂ █▄    ▂          
  ▄▁▄▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁████▄█▄▄█▄▁▁▁▁▁▁▆ ▁
  135 ms           Histogram: frequency by time          153 ms <

 Memory estimate: 41.09 MiB, allocs estimate: 6.

However, to analyze the simulation results I also need access to the observed variables, that are not contained in the solution object. Thus, this faster access method can't be used for all cases.

Simulation of the example below takes approx. 8 seconds on my computer. Accessing and exporting all 19 observed variables would take longer than the simulation in this case.

Are there any other faster options to get the full solution with all variables and export it to a file?

Example:

using ModelingToolkit, DifferentialEquations, Plots

@variables t
@connector function Pin(;name)
    sts = @variables v(t)=1.0 i(t)=1.0 [connect = Flow]
    ODESystem(Equation[], t, sts, []; name=name)
end

function Ground(;name)
    @named g = Pin()
    eqs = [g.v ~ 0]
    compose(ODESystem(eqs, t, [], []; name=name), g)
end

function OnePort(;name)
    @named p = Pin()
    @named n = Pin()
    sts = @variables v(t)=1.0 i(t)=1.0
    eqs = [
           v ~ p.v - n.v
           0 ~ p.i + n.i
           i ~ p.i
          ]
    compose(ODESystem(eqs, t, sts, []; name=name), p, n)
end

function Resistor(;name, R = 1.0)
    @named oneport = OnePort()
    @unpack v, i = oneport
    ps = @parameters R=R
    eqs = [
           v ~ i * R
          ]
    extend(ODESystem(eqs, t, [], ps; name=name), oneport)
end

function Capacitor(;name, C=1.0)
    @named oneport = OnePort()
    @unpack v, i = oneport
    ps = @parameters C=C
    D = Differential(t)
    eqs = [
           D(v) ~ i / C
          ]
    extend(ODESystem(eqs, t, [], ps; name=name), oneport)
end

function Voltage(;name, V=1.0)
    @named oneport = OnePort()
    @unpack v = oneport
    ps = @parameters V=V
    eqs = [
           v ~ V * sin(t)
          ]
    extend(ODESystem(eqs, t, [], ps; name=name), oneport)
end

R = 1.0
C = 1.0
V = 1.0
@named resistor = Resistor(R=R)
@named capacitor = Capacitor(C=C)
@named source = Voltage(V=V)
@named ground = Ground()

rc_eqs = [
          connect(source.p, resistor.p)
          connect(resistor.n, capacitor.p)
          connect(capacitor.n, source.n)
          connect(capacitor.n, ground.g)
         ]

@named _rc_model = ODESystem(rc_eqs, t)
@named rc_model = compose(_rc_model,
                          [resistor, capacitor, source, ground])
sys = structural_simplify(rc_model)
u0 = [
      capacitor.v => 0.0
     ]
prob = ODAEProblem(sys, u0, (0, 900000.0))
sol = solve(prob, Rodas42(), maxiters=1e7)

[YingboMa]

I think sol[capacitor.v] is the wrong interface to get the solution in a performant way. It's inherently extremely dynamic. Maybe we should just have something like:

getcv = getfunction(sol, capacitor.v) # This can be a bit slow, but that's unavoidable
getcv.(sol.u) # This will be fast

[isaacsas]

Is it possible to just cache that function internally after the first call to sol[capacitor.v], so subsequent calls are fast?

For actual variables, not observables, why is it slow? Couldn't there just be a Dict mapping from the symbol to the index in sol, or is that too simplified?

(@YingboMa sorry for my naïveté, I know you've thought about this a bunch so there are probably good reasons the latter can't be done or is slow -- I'm just curious!)

[YingboMa]

Is it possible to just cache that function internally after the first call to sol[capacitor.v], so subsequent calls are fast?

That's exactly what it does, and it's still slow due to the dynamism. Currently, it's doing

map(sol.u) do u
    getcv = getfunction(sol, capacitor.v)
    getcv(u)
end

And we should do

map(getfunction(sol, capacitor.v), sol.u)

[ChrisRackauckas]

Fixed with the new SymbolicIndexingInterface functions.