From a2b392a91ea807aa4bb4e7b17a13429775637d83 Mon Sep 17 00:00:00 2001
From: Chris Elrod <elrodc@gmail.com>
Date: Thu, 6 Jun 2024 22:39:20 -0400
Subject: [PATCH 01/19] Add RC Circuit example to MTK benchmarks.

---
 benchmarks/ModelingToolkit/RCCircuit.jmd | 113 +++++++++++++++++++++++
 1 file changed, 113 insertions(+)
 create mode 100644 benchmarks/ModelingToolkit/RCCircuit.jmd

diff --git a/benchmarks/ModelingToolkit/RCCircuit.jmd b/benchmarks/ModelingToolkit/RCCircuit.jmd
new file mode 100644
index 000000000..452d51bc8
--- /dev/null
+++ b/benchmarks/ModelingToolkit/RCCircuit.jmd
@@ -0,0 +1,113 @@
+---
+title: RC Circuit
+author: Yingbo Ma, Chris Elrod
+---
+
+When a model is defined using repeated components, JuliaSimCompiler is able to take advantage of this to scale efficiently by rerolling equations into loops. This option can be disabled by setting `loop=false`. Here, we build an RC circuit model with variable numbers of components to show scaling of compile and runtimes of MTK vs JuliaSimCompiler's three backends with and without loop rerolling.
+```julia
+using JuliaSimCompiler, ModelingToolkit, BenchmarkTools, ModelingToolkitStandardLibrary, CairoMakie
+using ModelingToolkitStandardLibrary.Blocks
+using ModelingToolkitStandardLibrary.Electrical
+
+const t = Blocks.t
+
+function build_system(n)
+  systems = @named begin
+    sine = Sine(frequency = 10)
+    source = Voltage()
+    resistors[1:n] = Resistor()
+    capacitors[1:n] = Capacitor()
+    ground = Ground()
+  end
+  systems = reduce(vcat, systems)
+  eqs = [connect(sine.output, source.V)
+       connect(source.p, resistors[1].p)
+       [connect(resistors[i].n, resistors[i + 1].p, capacitors[i].p)
+        for i in 1:(n - 1)]
+       connect(resistors[end].n, capacitors[end].p)
+       [connect(capacitors[i].n, source.n) for i in 1:n]
+       connect(source.n, ground.g)]
+  @named sys = ODESystem(eqs, t; systems)
+  u0 = [capacitors[i].v => float(i) for i in 1:n];
+  ps = [[resistors[i].R => 1 / i for i in 1:n];
+        [capacitors[i].C => 1 / i^2 for i in 1:n]]
+  return sys, u0, ps
+end
+
+function compile_run_problem(sys, u0, ps; target=JuliaSimCompiler.JuliaTarget())
+  tspan = (0.0, 10.0)
+  t0 = time()
+  if target === JuliaSimCompiler.JuliaTarget()
+    prob = ODEProblem(sys, u0, tspan, ps)
+  else
+    prob = ODEProblem(sys, target, u0, tspan, ps)
+  end
+  (; f, u0, p) = prob
+  ff = f.f
+  u0 .= 12.0
+  du = similar(u0)
+  ff(du, u0, p, 0.0)
+  t_fode = time() - t0
+  t_run = @belapsed $ff($du, $u0, $p, 0.0)
+  t_fode, t_run
+end
+
+const C = JuliaSimCompiler.CTarget();
+const LLVM = JuliaSimCompiler.llvm.LLVMTarget();
+
+ss_times = Matrix{Float64}(undef, length(N), 3);
+times = Matrix{NTuple{2,Float64}}(undef, length(N), 7);
+
+function run_and_time!(ss_times, times, i, n)
+  sys, u0, ps = build_system(n)
+  ss_times[i, 1] = @elapsed sys_mtk = structural_simplify(sys)
+  t_cce = @elapsed ir_state = JuliaSimCompiler.compressed_connection_expansion(sys)
+  ss_times[i, 2] = t_cce + @elapsed sys_jsir_scalar = structural_simplify(ir_state, loop=false)
+  ss_times[i, 3] = t_cce + @elapsed sys_jsir_loop = structural_simplify(ir_state)
+  times[i, 1] = compile_run_problem(sys_mtk, u0, ps)
+  times[i, 2] = compile_run_problem(sys_jsir_scalar, u0, ps)
+  times[i, 3] = compile_run_problem(sys_jsir_scalar, u0, ps, target=C)
+  times[i, 4] = compile_run_problem(sys_jsir_scalar, u0, ps, target=LLVM)
+  times[i, 5] = compile_run_problem(sys_jsir_loop, u0, ps)
+  times[i, 6] = compile_run_problem(sys_jsir_loop, u0, ps, target=C)
+  times[i, 7] = compile_run_problem(sys_jsir_loop, u0, ps, target=LLVM)
+  @show n, ss_times[i, :], times[i, :]
+end
+
+@time run_and_time!(ss_times, times, 1, 4); # precompile
+
+N = [5, 10, 20, 40, 60, 80, 160, 320, 480];
+
+for (i, n) in enumerate(N)
+  run_and_time!(ss_times, times, i, n)
+end
+
+f = Figure();
+let ax = Axis(f[1, 1]; title="Structural Simplify Time")
+  for (j, label) in enumerate(("MTK", "JSIR-Scalar", "JSIR-Loop"))
+    ts = @view(ss_times[:, j])
+    lines!(N, ts; label)
+  end
+  axislegend(ax)
+end
+for (i, timecat) in enumerate(("ODEProblem + f!", "Run"))
+  title = timecat * " Time"
+  ax = Axis(f[i+1, 1]; title)
+  for (j, label) in enumerate(("MTK", "JSIR - Scalar - Julia", "JSIR - Scalar - C", "JSIR - Scalar - LLVM", "JSIR - Loop - Julia", "JSIR - Loop - C", "JSIR - Loop - LLVM"))
+    ts = getindex.(@view(times[:, j]), i)
+    lines!(N, ts; label)
+  end
+  axislegend(ax)
+end
+f
+```
+All three backends compiled more quickly with loops, but the C and LLVM backends are so much quicker to compile than the Julia backend that this made much less difference for them.
+The impact on runtime was more varied.
+
+## Appendix
+
+```julia, echo = false
+using SciMLBenchmarks
+SciMLBenchmarks.bench_footer(WEAVE_ARGS[:folder],WEAVE_ARGS[:file])
+```
+

From b3a62c3f05980d608006aa1966fcad81f06f56d3 Mon Sep 17 00:00:00 2001
From: Chris Elrod <elrodc@gmail.com>
Date: Thu, 6 Jun 2024 22:45:53 -0400
Subject: [PATCH 02/19] move `axislegend`s to topleft.

---
 benchmarks/ModelingToolkit/ThermalFluid.jmd | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/benchmarks/ModelingToolkit/ThermalFluid.jmd b/benchmarks/ModelingToolkit/ThermalFluid.jmd
index 6ef781fc9..55983e028 100644
--- a/benchmarks/ModelingToolkit/ThermalFluid.jmd
+++ b/benchmarks/ModelingToolkit/ThermalFluid.jmd
@@ -309,7 +309,7 @@ let ax = Axis(f[1, 1]; title="Structural Simplify Time")
     ts = @view(ss_times[:, j])
     lines!(N_states, ts; label)
   end
-  axislegend(ax)
+  axislegend(ax, position = :lt)
 end
 for (i, timecat) in enumerate(("ODEProblem + f!", "Run"))
   title = timecat * " Time"
@@ -318,7 +318,7 @@ for (i, timecat) in enumerate(("ODEProblem + f!", "Run"))
     ts = getindex.(@view(times[:, j]), i)
     lines!(N_states, ts; label)
   end
-  axislegend(ax)
+  axislegend(ax, position = :lt)
 end
 f
 ```

From 07ef5e4c037891443b0b498540a271f5237809b6 Mon Sep 17 00:00:00 2001
From: Chris Elrod <elrodc@gmail.com>
Date: Fri, 7 Jun 2024 11:52:31 -0400
Subject: [PATCH 03/19] fix order

---
 benchmarks/ModelingToolkit/RCCircuit.jmd | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/benchmarks/ModelingToolkit/RCCircuit.jmd b/benchmarks/ModelingToolkit/RCCircuit.jmd
index 452d51bc8..09511e3a5 100644
--- a/benchmarks/ModelingToolkit/RCCircuit.jmd
+++ b/benchmarks/ModelingToolkit/RCCircuit.jmd
@@ -55,9 +55,6 @@ end
 const C = JuliaSimCompiler.CTarget();
 const LLVM = JuliaSimCompiler.llvm.LLVMTarget();
 
-ss_times = Matrix{Float64}(undef, length(N), 3);
-times = Matrix{NTuple{2,Float64}}(undef, length(N), 7);
-
 function run_and_time!(ss_times, times, i, n)
   sys, u0, ps = build_system(n)
   ss_times[i, 1] = @elapsed sys_mtk = structural_simplify(sys)
@@ -74,9 +71,12 @@ function run_and_time!(ss_times, times, i, n)
   @show n, ss_times[i, :], times[i, :]
 end
 
+N = [5, 10, 20, 40, 60, 80, 160, 320, 480];
+ss_times = Matrix{Float64}(undef, length(N), 3);
+times = Matrix{NTuple{2,Float64}}(undef, length(N), 7);
+
 @time run_and_time!(ss_times, times, 1, 4); # precompile
 
-N = [5, 10, 20, 40, 60, 80, 160, 320, 480];
 
 for (i, n) in enumerate(N)
   run_and_time!(ss_times, times, i, n)

From c792abd510fba1b0ef0d338d27abf446daf0bc25 Mon Sep 17 00:00:00 2001
From: Chris Elrod <elrodc@gmail.com>
Date: Fri, 7 Jun 2024 12:18:26 -0400
Subject: [PATCH 04/19] Add mechanism for setting size limit per method.
 Currently, disable MTK > 1000.

---
 benchmarks/ModelingToolkit/RCCircuit.jmd | 34 ++++++++++++++----------
 1 file changed, 20 insertions(+), 14 deletions(-)

diff --git a/benchmarks/ModelingToolkit/RCCircuit.jmd b/benchmarks/ModelingToolkit/RCCircuit.jmd
index 09511e3a5..8573dd280 100644
--- a/benchmarks/ModelingToolkit/RCCircuit.jmd
+++ b/benchmarks/ModelingToolkit/RCCircuit.jmd
@@ -55,31 +55,37 @@ end
 const C = JuliaSimCompiler.CTarget();
 const LLVM = JuliaSimCompiler.llvm.LLVMTarget();
 
-function run_and_time!(ss_times, times, i, n)
+function run_and_time!(ss_times, times, max_sizes, i, n)
   sys, u0, ps = build_system(n)
-  ss_times[i, 1] = @elapsed sys_mtk = structural_simplify(sys)
+  if n <= max_sizes[1]
+    ss_times[i, 1] = @elapsed sys_mtk = structural_simplify(sys)
+    times[i, 1] = compile_run_problem(sys_mtk, u0, ps)
+  end
   t_cce = @elapsed ir_state = JuliaSimCompiler.compressed_connection_expansion(sys)
   ss_times[i, 2] = t_cce + @elapsed sys_jsir_scalar = structural_simplify(ir_state, loop=false)
   ss_times[i, 3] = t_cce + @elapsed sys_jsir_loop = structural_simplify(ir_state)
-  times[i, 1] = compile_run_problem(sys_mtk, u0, ps)
-  times[i, 2] = compile_run_problem(sys_jsir_scalar, u0, ps)
-  times[i, 3] = compile_run_problem(sys_jsir_scalar, u0, ps, target=C)
-  times[i, 4] = compile_run_problem(sys_jsir_scalar, u0, ps, target=LLVM)
-  times[i, 5] = compile_run_problem(sys_jsir_loop, u0, ps)
-  times[i, 6] = compile_run_problem(sys_jsir_loop, u0, ps, target=C)
-  times[i, 7] = compile_run_problem(sys_jsir_loop, u0, ps, target=LLVM)
+  n <= max_sizes[2] && (times[i, 2] = compile_run_problem(sys_jsir_scalar, u0, ps))
+  n <= max_sizes[3] && (times[i, 3] = compile_run_problem(sys_jsir_scalar, u0, ps, target=C))
+  n <= max_sizes[4] && (times[i, 4] = compile_run_problem(sys_jsir_scalar, u0, ps, target=LLVM))
+  n <= max_sizes[5] && (times[i, 5] = compile_run_problem(sys_jsir_loop, u0, ps))
+  n <= max_sizes[6] && (times[i, 6] = compile_run_problem(sys_jsir_loop, u0, ps, target=C))
+  n <= max_sizes[7] && (times[i, 7] = compile_run_problem(sys_jsir_loop, u0, ps, target=LLVM))
   @show n, ss_times[i, :], times[i, :]
 end
 
-N = [5, 10, 20, 40, 60, 80, 160, 320, 480];
-ss_times = Matrix{Float64}(undef, length(N), 3);
-times = Matrix{NTuple{2,Float64}}(undef, length(N), 7);
+N = vcat([5, 10, 20, 40, 60, 80, 160, 320, 480, 640, 800], 1_000:1_000:10_000);
+
+# max size we test per method
+max_sizes = [1000, 10_000, 10_000, 10_000, 10_000, 10_000, 10_000];
 
-@time run_and_time!(ss_times, times, 1, 4); # precompile
+# NaN-initialize so Makie will ignore incomplete
+ss_times = fill(NaN, length(N), 3);
+times = fill((NaN,NaN), length(N), length(max_sizes));
 
+@time run_and_time!(ss_times, times, max_sizes, 1, 4); # precompile
 
 for (i, n) in enumerate(N)
-  run_and_time!(ss_times, times, i, n)
+  run_and_time!(ss_times, times, max_sizes, i, n)
 end
 
 f = Figure();

From 4af003a94a5cd8db72dcde62714101492aad5b56 Mon Sep 17 00:00:00 2001
From: Chris Elrod <elrodc@gmail.com>
Date: Fri, 7 Jun 2024 13:14:52 -0400
Subject: [PATCH 05/19] set better limits, remove `.= 12` initialization

---
 benchmarks/ModelingToolkit/RCCircuit.jmd | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/benchmarks/ModelingToolkit/RCCircuit.jmd b/benchmarks/ModelingToolkit/RCCircuit.jmd
index 8573dd280..4ac5e88f0 100644
--- a/benchmarks/ModelingToolkit/RCCircuit.jmd
+++ b/benchmarks/ModelingToolkit/RCCircuit.jmd
@@ -44,7 +44,6 @@ function compile_run_problem(sys, u0, ps; target=JuliaSimCompiler.JuliaTarget())
   end
   (; f, u0, p) = prob
   ff = f.f
-  u0 .= 12.0
   du = similar(u0)
   ff(du, u0, p, 0.0)
   t_fode = time() - t0
@@ -76,7 +75,7 @@ end
 N = vcat([5, 10, 20, 40, 60, 80, 160, 320, 480, 640, 800], 1_000:1_000:10_000);
 
 # max size we test per method
-max_sizes = [1000, 10_000, 10_000, 10_000, 10_000, 10_000, 10_000];
+max_sizes = [1000, 2_000, 4_000, 4_000, 10_000, 10_000, 10_000];
 
 # NaN-initialize so Makie will ignore incomplete
 ss_times = fill(NaN, length(N), 3);

From 1399ab03e1bfa554ed2cad90665a11c8b31e4e6f Mon Sep 17 00:00:00 2001
From: Chris Elrod <elrodc@gmail.com>
Date: Fri, 7 Jun 2024 14:44:38 -0400
Subject: [PATCH 06/19] axislegend position

---
 benchmarks/ModelingToolkit/RCCircuit.jmd | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/benchmarks/ModelingToolkit/RCCircuit.jmd b/benchmarks/ModelingToolkit/RCCircuit.jmd
index 4ac5e88f0..58ff63151 100644
--- a/benchmarks/ModelingToolkit/RCCircuit.jmd
+++ b/benchmarks/ModelingToolkit/RCCircuit.jmd
@@ -93,7 +93,7 @@ let ax = Axis(f[1, 1]; title="Structural Simplify Time")
     ts = @view(ss_times[:, j])
     lines!(N, ts; label)
   end
-  axislegend(ax)
+  axislegend(ax, position = :lt)
 end
 for (i, timecat) in enumerate(("ODEProblem + f!", "Run"))
   title = timecat * " Time"
@@ -102,7 +102,7 @@ for (i, timecat) in enumerate(("ODEProblem + f!", "Run"))
     ts = getindex.(@view(times[:, j]), i)
     lines!(N, ts; label)
   end
-  axislegend(ax)
+  axislegend(ax, position = :lt)
 end
 f
 ```

From dbe30910ae5981749dcccf298f87ee8ead1ce0eb Mon Sep 17 00:00:00 2001
From: Chris Elrod <elrodc@gmail.com>
Date: Mon, 10 Jun 2024 01:49:53 -0400
Subject: [PATCH 07/19] add `solve` times

---
 benchmarks/ModelingToolkit/RCCircuit.jmd | 35 +++++++++++++-----------
 1 file changed, 19 insertions(+), 16 deletions(-)

diff --git a/benchmarks/ModelingToolkit/RCCircuit.jmd b/benchmarks/ModelingToolkit/RCCircuit.jmd
index 58ff63151..22213969a 100644
--- a/benchmarks/ModelingToolkit/RCCircuit.jmd
+++ b/benchmarks/ModelingToolkit/RCCircuit.jmd
@@ -5,7 +5,7 @@ author: Yingbo Ma, Chris Elrod
 
 When a model is defined using repeated components, JuliaSimCompiler is able to take advantage of this to scale efficiently by rerolling equations into loops. This option can be disabled by setting `loop=false`. Here, we build an RC circuit model with variable numbers of components to show scaling of compile and runtimes of MTK vs JuliaSimCompiler's three backends with and without loop rerolling.
 ```julia
-using JuliaSimCompiler, ModelingToolkit, BenchmarkTools, ModelingToolkitStandardLibrary, CairoMakie
+using JuliaSimCompiler, ModelingToolkit, OrdinaryDiffEq, BenchmarkTools, ModelingToolkitStandardLibrary, CairoMakie
 using ModelingToolkitStandardLibrary.Blocks
 using ModelingToolkitStandardLibrary.Electrical
 
@@ -34,41 +34,44 @@ function build_system(n)
   return sys, u0, ps
 end
 
-function compile_run_problem(sys, u0, ps; target=JuliaSimCompiler.JuliaTarget())
+function compile_run_problem(sys, u0, ps; target=JuliaSimCompiler.JuliaTarget(), duref=nothing)
   tspan = (0.0, 10.0)
   t0 = time()
-  if target === JuliaSimCompiler.JuliaTarget()
-    prob = ODEProblem(sys, u0, tspan, ps)
+  prob = if target === JuliaSimCompiler.JuliaTarget()
+    ODEProblem(sys, u0, tspan, ps)
   else
-    prob = ODEProblem(sys, target, u0, tspan, ps)
+    ODEProblem(sys, target, u0, tspan, ps)
   end
   (; f, u0, p) = prob
   ff = f.f
   du = similar(u0)
   ff(du, u0, p, 0.0)
   t_fode = time() - t0
+  duref === nothing || @assert duref ≈ du
   t_run = @belapsed $ff($du, $u0, $p, 0.0)
-  t_fode, t_run
+  t_solve = @elapsed solve(prob, Rodas5(autodiff = false))
+  (t_fode, t_run, t_solve), du
 end
 
 const C = JuliaSimCompiler.CTarget();
 const LLVM = JuliaSimCompiler.llvm.LLVMTarget();
 
 function run_and_time!(ss_times, times, max_sizes, i, n)
-  sys, u0, ps = build_system(n)
+  sys, u0, ps = build_system(n);
   if n <= max_sizes[1]
     ss_times[i, 1] = @elapsed sys_mtk = structural_simplify(sys)
-    times[i, 1] = compile_run_problem(sys_mtk, u0, ps)
+    times[i, 1], _ = compile_run_problem(sys_mtk, u0, ps)
   end
   t_cce = @elapsed ir_state = JuliaSimCompiler.compressed_connection_expansion(sys)
   ss_times[i, 2] = t_cce + @elapsed sys_jsir_scalar = structural_simplify(ir_state, loop=false)
   ss_times[i, 3] = t_cce + @elapsed sys_jsir_loop = structural_simplify(ir_state)
-  n <= max_sizes[2] && (times[i, 2] = compile_run_problem(sys_jsir_scalar, u0, ps))
-  n <= max_sizes[3] && (times[i, 3] = compile_run_problem(sys_jsir_scalar, u0, ps, target=C))
-  n <= max_sizes[4] && (times[i, 4] = compile_run_problem(sys_jsir_scalar, u0, ps, target=LLVM))
-  n <= max_sizes[5] && (times[i, 5] = compile_run_problem(sys_jsir_loop, u0, ps))
-  n <= max_sizes[6] && (times[i, 6] = compile_run_problem(sys_jsir_loop, u0, ps, target=C))
-  n <= max_sizes[7] && (times[i, 7] = compile_run_problem(sys_jsir_loop, u0, ps, target=LLVM))
+  duref = nothing
+  n <= max_sizes[2] && ((times[i, 2], duref) = compile_run_problem(sys_jsir_scalar, u0, ps))
+  n <= max_sizes[3] && ((times[i, 3], duref) = compile_run_problem(sys_jsir_scalar, u0, ps; target=C, duref))
+  n <= max_sizes[4] && ((times[i, 4], duref) = compile_run_problem(sys_jsir_scalar, u0, ps; target=LLVM, duref))
+  n <= max_sizes[5] && ((times[i, 5], duref) = compile_run_problem(sys_jsir_loop, u0, ps; duref))
+  n <= max_sizes[6] && ((times[i, 6], duref) = compile_run_problem(sys_jsir_loop, u0, ps; target=C, duref))
+  n <= max_sizes[7] && ((times[i, 7], duref) = compile_run_problem(sys_jsir_loop, u0, ps; target=LLVM, duref))
   @show n, ss_times[i, :], times[i, :]
 end
 
@@ -79,7 +82,7 @@ max_sizes = [1000, 2_000, 4_000, 4_000, 10_000, 10_000, 10_000];
 
 # NaN-initialize so Makie will ignore incomplete
 ss_times = fill(NaN, length(N), 3);
-times = fill((NaN,NaN), length(N), length(max_sizes));
+times = fill((NaN,NaN,NaN), length(N), length(max_sizes));
 
 @time run_and_time!(ss_times, times, max_sizes, 1, 4); # precompile
 
@@ -95,7 +98,7 @@ let ax = Axis(f[1, 1]; title="Structural Simplify Time")
   end
   axislegend(ax, position = :lt)
 end
-for (i, timecat) in enumerate(("ODEProblem + f!", "Run"))
+for (i, timecat) in enumerate(("ODEProblem + f!", "Run", "Solve"))
   title = timecat * " Time"
   ax = Axis(f[i+1, 1]; title)
   for (j, label) in enumerate(("MTK", "JSIR - Scalar - Julia", "JSIR - Scalar - C", "JSIR - Scalar - LLVM", "JSIR - Loop - Julia", "JSIR - Loop - C", "JSIR - Loop - LLVM"))

From dce51ccccfb27568e62e438e2c9395c1c6c1fef6 Mon Sep 17 00:00:00 2001
From: Chris Elrod <elrodc@gmail.com>
Date: Mon, 10 Jun 2024 02:26:16 -0400
Subject: [PATCH 08/19] sparse=true

---
 benchmarks/ModelingToolkit/RCCircuit.jmd | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/benchmarks/ModelingToolkit/RCCircuit.jmd b/benchmarks/ModelingToolkit/RCCircuit.jmd
index 22213969a..7bc5646fd 100644
--- a/benchmarks/ModelingToolkit/RCCircuit.jmd
+++ b/benchmarks/ModelingToolkit/RCCircuit.jmd
@@ -38,9 +38,9 @@ function compile_run_problem(sys, u0, ps; target=JuliaSimCompiler.JuliaTarget(),
   tspan = (0.0, 10.0)
   t0 = time()
   prob = if target === JuliaSimCompiler.JuliaTarget()
-    ODEProblem(sys, u0, tspan, ps)
+    ODEProblem(sys, u0, tspan, ps; sparse = true)
   else
-    ODEProblem(sys, target, u0, tspan, ps)
+    ODEProblem(sys, target, u0, tspan, ps; sparse = true)
   end
   (; f, u0, p) = prob
   ff = f.f

From 9de8afc1a4fd0f7cdc1a3f1c54831528d95d7a22 Mon Sep 17 00:00:00 2001
From: avinashresearch1 <avinash.subramanian.research@gmail.com>
Date: Tue, 11 Jun 2024 03:12:00 +0200
Subject: [PATCH 09/19] feat: add openmodelica comparison times

---
 Project.toml                             |   3 +
 benchmarks/ModelingToolkit/RCCircuit.jmd |  49 +++++++++
 benchmarks/ModelingToolkit/RC_Circuit.mo | 128 +++++++++++++++++++++++
 3 files changed, 180 insertions(+)
 create mode 100644 benchmarks/ModelingToolkit/RC_Circuit.mo

diff --git a/Project.toml b/Project.toml
index 8c359c281..6644e4029 100644
--- a/Project.toml
+++ b/Project.toml
@@ -4,10 +4,13 @@ authors = ["Chris Rackauckas, Yingbo Ma, Vaibhav Dixit"]
 version = "0.1.3"
 
 [deps]
+CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
+DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
 Git = "d7ba0133-e1db-5d97-8f8c-041e4b3a1eb2"
 IJulia = "7073ff75-c697-5162-941a-fcdaad2a7d2a"
 InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
 Markdown = "d6f4376e-aef5-505a-96c1-9c027394607a"
+OMJulia = "0f4fe800-344e-11e9-2949-fb537ad918e1"
 Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 Weave = "44d3d7a6-8a23-5bf8-98c5-b353f8df5ec9"
 
diff --git a/benchmarks/ModelingToolkit/RCCircuit.jmd b/benchmarks/ModelingToolkit/RCCircuit.jmd
index 7bc5646fd..90f6bc6a5 100644
--- a/benchmarks/ModelingToolkit/RCCircuit.jmd
+++ b/benchmarks/ModelingToolkit/RCCircuit.jmd
@@ -119,3 +119,52 @@ using SciMLBenchmarks
 SciMLBenchmarks.bench_footer(WEAVE_ARGS[:folder],WEAVE_ARGS[:file])
 ```
 
+## Time OpenModelica
+
+```julia
+using OMJulia, CSV, DataFrames
+mod = OMJulia.OMCSession();
+
+OMJulia.sendExpression(mod, "getVersion()")
+
+modelicafile = joinpath(@__DIR__, "RC_Circuit.mo")
+
+resultfile = joinpath(@__DIR__, "RC_Circuit_modelica_res.csv")
+
+@show "Start OpenModelica Timings"
+
+# TODO: Sync these test cases with the N used in MTK: 
+N_x = [5, 10, 20, 40, 60, 80, 160, 320, 480, 640, 800, 1000, 2000, 3000, 4000, 5000, 6000, 7000, 8000, 9000, 10000, 20000]
+
+for N in N_x
+    @show N
+    totaltime = @elapsed res = begin
+        ModelicaSystem(mod, modelicafile, "RC_Circuit.Test.RC_Circuit_MTK_test_$N")
+        sendExpression(mod, "simulate(RC_Circuit.Test.RC_Circuit_MTK_test_$N)")
+    end
+    #runtime = res["timeTotal"]
+    @assert res["messages"][1:11] == "LOG_SUCCESS"
+    @show totaltime
+    # total_times[i, 5] = totaltime TODO: Add to comparison table
+end
+
+OMJulia.quit(mod) 
+```
+
+Dymola requires a license server and thus cannot be hosted. This was run locally for the
+following times:
+
+```julia
+translation_and_total_times = [
+7.027, 7.237
+11.295, 11.798
+16.681, 17.646
+22.125, 23.839
+27.529, 29.82
+33.282, 36.622
+39.007, 43.088
+44.825, 51.601
+50.281, 56.676
+] # TODO: I will add other times once the Dymola license server is back up.
+#total_times[:, 6] = translation_and_total_times[1:length(N_x),2]
+```
\ No newline at end of file
diff --git a/benchmarks/ModelingToolkit/RC_Circuit.mo b/benchmarks/ModelingToolkit/RC_Circuit.mo
new file mode 100644
index 000000000..b76068095
--- /dev/null
+++ b/benchmarks/ModelingToolkit/RC_Circuit.mo
@@ -0,0 +1,128 @@
+within ;
+package RC_Circuit
+
+  model RC_Circuit_Base
+    parameter Integer N=1000;
+    Modelica.Electrical.Analog.Sources.SineVoltage source(f=0.05, V=1);
+    Modelica.Electrical.Analog.Basic.Ground g;
+    Modelica.Electrical.Analog.Basic.Resistor R[N](each R=1.0);
+    Modelica.Electrical.Analog.Basic.Capacitor C[N](each C=1.0, each v(fixed=true));
+
+  equation
+    connect(g.p, source.n);
+    connect(R[1].p, source.p);
+    connect(C[1].p, R[1].n);
+    connect(C[1].n, g.p);
+    for i in 2:N loop
+      connect(R[i].n, C[i].p);
+      connect(R[i].p, R[i-1].n);
+      connect(C[i].n, g.p);
+    end for;
+            annotation (Icon(coordinateSystem(preserveAspectRatio=false)), Diagram(
+          coordinateSystem(preserveAspectRatio=false)),
+      experiment(StopTime=10, __Dymola_Algorithm="Dassl"));
+  end RC_Circuit_Base;
+
+  package Test
+    model RC_Circuit_MTK_test_5
+      extends RC_Circuit_Base(N=5);
+    end RC_Circuit_MTK_test_5;
+
+    model RC_Circuit_MTK_test_10
+      extends RC_Circuit_Base(N=10);
+    end RC_Circuit_MTK_test_10;
+
+    model RC_Circuit_MTK_test_20
+      extends RC_Circuit_Base(N=20);
+    end RC_Circuit_MTK_test_20;
+
+    model RC_Circuit_MTK_test_40
+      extends RC_Circuit_Base(N=40);
+    end RC_Circuit_MTK_test_40;
+
+    model RC_Circuit_MTK_test_60
+      extends RC_Circuit_Base(N=60);
+    end RC_Circuit_MTK_test_60;
+
+    model RC_Circuit_MTK_test_80
+      extends RC_Circuit_Base(N=80);
+    end RC_Circuit_MTK_test_80;
+
+    model RC_Circuit_MTK_test_160
+      extends RC_Circuit_Base(N=160);
+    end RC_Circuit_MTK_test_160;
+
+    model RC_Circuit_MTK_test_320
+      extends RC_Circuit_Base(N=320);
+    end RC_Circuit_MTK_test_320;
+
+    model RC_Circuit_MTK_test_480
+      extends RC_Circuit_Base(N=480);
+    end RC_Circuit_MTK_test_480;
+
+    model RC_Circuit_MTK_test_640
+      extends RC_Circuit_Base(N=640);
+    end RC_Circuit_MTK_test_640;
+
+    model RC_Circuit_MTK_test_800
+      extends RC_Circuit_Base(N=800);
+    end RC_Circuit_MTK_test_800;
+
+    model RC_Circuit_MTK_test_1000
+      extends RC_Circuit_Base(N=1000);
+    end RC_Circuit_MTK_test_1000;
+
+    model RC_Circuit_MTK_test_2000
+      extends RC_Circuit_Base(N=2000);
+    end RC_Circuit_MTK_test_2000;
+
+    model RC_Circuit_MTK_test_3000
+      extends RC_Circuit_Base(N=3000);
+    end RC_Circuit_MTK_test_3000;
+
+    model RC_Circuit_MTK_test_4000
+      extends RC_Circuit_Base(N=4000);
+    end RC_Circuit_MTK_test_4000;
+
+    model RC_Circuit_MTK_test_5000
+      extends RC_Circuit_Base(N=5000);
+    end RC_Circuit_MTK_test_5000;
+
+    model RC_Circuit_MTK_test_6000
+      extends RC_Circuit_Base(N=6000);
+    end RC_Circuit_MTK_test_6000;
+
+    model RC_Circuit_MTK_test_7000
+      extends RC_Circuit_Base(N=7000);
+    end RC_Circuit_MTK_test_7000;
+
+    model RC_Circuit_MTK_test_8000
+      extends RC_Circuit_Base(N=8000);
+    end RC_Circuit_MTK_test_8000;
+
+    model RC_Circuit_MTK_test_9000
+      extends RC_Circuit_Base(N=9000);
+    end RC_Circuit_MTK_test_9000;
+
+    model RC_Circuit_MTK_test_10000
+      extends RC_Circuit_Base(N=10000);
+    end RC_Circuit_MTK_test_10000;
+
+    model RC_Circuit_MTK_test_20000
+      extends RC_Circuit_Base(N=20000);
+    end RC_Circuit_MTK_test_20000;
+  end Test;
+
+  package TimingFunctions
+
+    function tic
+      "Function to record the internal time in [s]"
+      output Real tic;
+    algorithm
+      (ms_tic, sec_tic, min_tic, hour_tic, day_tic, mon_tic, year_tic) := Modelica.Utilities.System.getTime();
+      tic := (ms_tic*0.001) + sec_tic + (min_tic*60) + (hour_tic*3600) + (day_tic*86400);
+    end tic;
+
+  end TimingFunctions;
+    
+end RC_Circuit;
\ No newline at end of file

From a71c5dc0c102daf549121f52fac3bf79510b84fb Mon Sep 17 00:00:00 2001
From: Chris Elrod <elrodc@gmail.com>
Date: Tue, 11 Jun 2024 00:52:54 -0400
Subject: [PATCH 10/19] no classes when not rerolling loops so we can avoid the
 mass matrix

---
 benchmarks/ModelingToolkit/RCCircuit.jmd | 5 ++---
 1 file changed, 2 insertions(+), 3 deletions(-)

diff --git a/benchmarks/ModelingToolkit/RCCircuit.jmd b/benchmarks/ModelingToolkit/RCCircuit.jmd
index 90f6bc6a5..8cb9ff432 100644
--- a/benchmarks/ModelingToolkit/RCCircuit.jmd
+++ b/benchmarks/ModelingToolkit/RCCircuit.jmd
@@ -62,9 +62,8 @@ function run_and_time!(ss_times, times, max_sizes, i, n)
     ss_times[i, 1] = @elapsed sys_mtk = structural_simplify(sys)
     times[i, 1], _ = compile_run_problem(sys_mtk, u0, ps)
   end
-  t_cce = @elapsed ir_state = JuliaSimCompiler.compressed_connection_expansion(sys)
-  ss_times[i, 2] = t_cce + @elapsed sys_jsir_scalar = structural_simplify(ir_state, loop=false)
-  ss_times[i, 3] = t_cce + @elapsed sys_jsir_loop = structural_simplify(ir_state)
+  ss_times[i, 2] = @elapsed sys_jsir_scalar = structural_simplify(IRSystem(sys), loop=false)
+  ss_times[i, 3] = @elapsed sys_jsir_loop = structural_simplify(JuliaSimCompiler.compressed_connection_expansion(sys))
   duref = nothing
   n <= max_sizes[2] && ((times[i, 2], duref) = compile_run_problem(sys_jsir_scalar, u0, ps))
   n <= max_sizes[3] && ((times[i, 3], duref) = compile_run_problem(sys_jsir_scalar, u0, ps; target=C, duref))

From 1547edc961ac89da5437035803c6cce2538572eb Mon Sep 17 00:00:00 2001
From: Chris Elrod <elrodc@gmail.com>
Date: Tue, 11 Jun 2024 01:10:10 -0400
Subject: [PATCH 11/19] a few fixes

---
 benchmarks/ModelingToolkit/RCCircuit.jmd | 16 ++++++++--------
 1 file changed, 8 insertions(+), 8 deletions(-)

diff --git a/benchmarks/ModelingToolkit/RCCircuit.jmd b/benchmarks/ModelingToolkit/RCCircuit.jmd
index 8cb9ff432..0216e8bb9 100644
--- a/benchmarks/ModelingToolkit/RCCircuit.jmd
+++ b/benchmarks/ModelingToolkit/RCCircuit.jmd
@@ -64,13 +64,13 @@ function run_and_time!(ss_times, times, max_sizes, i, n)
   end
   ss_times[i, 2] = @elapsed sys_jsir_scalar = structural_simplify(IRSystem(sys), loop=false)
   ss_times[i, 3] = @elapsed sys_jsir_loop = structural_simplify(JuliaSimCompiler.compressed_connection_expansion(sys))
-  duref = nothing
-  n <= max_sizes[2] && ((times[i, 2], duref) = compile_run_problem(sys_jsir_scalar, u0, ps))
-  n <= max_sizes[3] && ((times[i, 3], duref) = compile_run_problem(sys_jsir_scalar, u0, ps; target=C, duref))
-  n <= max_sizes[4] && ((times[i, 4], duref) = compile_run_problem(sys_jsir_scalar, u0, ps; target=LLVM, duref))
-  n <= max_sizes[5] && ((times[i, 5], duref) = compile_run_problem(sys_jsir_loop, u0, ps; duref))
-  n <= max_sizes[6] && ((times[i, 6], duref) = compile_run_problem(sys_jsir_loop, u0, ps; target=C, duref))
-  n <= max_sizes[7] && ((times[i, 7], duref) = compile_run_problem(sys_jsir_loop, u0, ps; target=LLVM, duref))
+  oderef = daeref = nothing
+  n <= max_sizes[2] && ((times[i, 2], oderef) = compile_run_problem(sys_jsir_scalar, u0, ps; duref = oderef))
+  n <= max_sizes[3] && ((times[i, 3], oderef) = compile_run_problem(sys_jsir_scalar, u0, ps; target=C, duref = oderef))
+  n <= max_sizes[4] && ((times[i, 4], oderef) = compile_run_problem(sys_jsir_scalar, u0, ps; target=LLVM, duref = oderef))
+  n <= max_sizes[5] && ((times[i, 5], deeref) = compile_run_problem(sys_jsir_loop, u0, ps; duref = daeref))
+  n <= max_sizes[6] && ((times[i, 6], daeref) = compile_run_problem(sys_jsir_loop, u0, ps; target=C, duref = daeref))
+  n <= max_sizes[7] && ((times[i, 7], daeref) = compile_run_problem(sys_jsir_loop, u0, ps; target=LLVM, duref = daeref))
   @show n, ss_times[i, :], times[i, :]
 end
 
@@ -166,4 +166,4 @@ translation_and_total_times = [
 50.281, 56.676
 ] # TODO: I will add other times once the Dymola license server is back up.
 #total_times[:, 6] = translation_and_total_times[1:length(N_x),2]
-```
\ No newline at end of file
+```

From a5cc0540a61ccec7b52afd3d9cb8284b47b05671 Mon Sep 17 00:00:00 2001
From: Chris Elrod <elrodc@gmail.com>
Date: Tue, 11 Jun 2024 01:29:24 -0400
Subject: [PATCH 12/19] increase max sizes

---
 benchmarks/ModelingToolkit/RCCircuit.jmd | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/benchmarks/ModelingToolkit/RCCircuit.jmd b/benchmarks/ModelingToolkit/RCCircuit.jmd
index 0216e8bb9..13d142747 100644
--- a/benchmarks/ModelingToolkit/RCCircuit.jmd
+++ b/benchmarks/ModelingToolkit/RCCircuit.jmd
@@ -77,7 +77,7 @@ end
 N = vcat([5, 10, 20, 40, 60, 80, 160, 320, 480, 640, 800], 1_000:1_000:10_000);
 
 # max size we test per method
-max_sizes = [1000, 2_000, 4_000, 4_000, 10_000, 10_000, 10_000];
+max_sizes = [1000, 4_000, 10_000, 10_000, 10_000, 10_000, 10_000];
 
 # NaN-initialize so Makie will ignore incomplete
 ss_times = fill(NaN, length(N), 3);

From 9db1429602305725d8b3635d4477642caf17774e Mon Sep 17 00:00:00 2001
From: Chris Elrod <elrodc@gmail.com>
Date: Tue, 11 Jun 2024 15:28:49 -0400
Subject: [PATCH 13/19] adjust axis scales, place legend to right

---
 benchmarks/ModelingToolkit/RCCircuit.jmd | 22 +++++++++++-----------
 1 file changed, 11 insertions(+), 11 deletions(-)

diff --git a/benchmarks/ModelingToolkit/RCCircuit.jmd b/benchmarks/ModelingToolkit/RCCircuit.jmd
index 13d142747..2786db94c 100644
--- a/benchmarks/ModelingToolkit/RCCircuit.jmd
+++ b/benchmarks/ModelingToolkit/RCCircuit.jmd
@@ -89,22 +89,22 @@ for (i, n) in enumerate(N)
   run_and_time!(ss_times, times, max_sizes, i, n)
 end
 
-f = Figure();
-let ax = Axis(f[1, 1]; title="Structural Simplify Time")
-  for (j, label) in enumerate(("MTK", "JSIR-Scalar", "JSIR-Loop"))
-    ts = @view(ss_times[:, j])
-    lines!(N, ts; label)
+f = Figure(size=(800,1200));
+let ax = Axis(f[1, 1]; yscale = log10, xscale = log10, title="Structural Simplify Time")
+  _lines = map(eachcol(ss_times)) do ts
+    lines!(N, ts)
   end
-  axislegend(ax, position = :lt)
+  names = ("MTK", "JSIR-Scalar", "JSIR-Loop")
+  Legend(f([1,2], _lines, names)
 end
 for (i, timecat) in enumerate(("ODEProblem + f!", "Run", "Solve"))
   title = timecat * " Time"
-  ax = Axis(f[i+1, 1]; title)
-  for (j, label) in enumerate(("MTK", "JSIR - Scalar - Julia", "JSIR - Scalar - C", "JSIR - Scalar - LLVM", "JSIR - Loop - Julia", "JSIR - Loop - C", "JSIR - Loop - LLVM"))
-    ts = getindex.(@view(times[:, j]), i)
-    lines!(N, ts; label)
+  ax = Axis(f[i+1, 1]; yscale = log10, xscale = log10, title)
+  names = ("MTK", "JSIR - Scalar - Julia", "JSIR - Scalar - C", "JSIR - Scalar - LLVM", "JSIR - Loop - Julia", "JSIR - Loop - C", "JSIR - Loop - LLVM")
+  _lines = map(eachcol(times)) do ts
+    lines!(N, getindex.(ts, i))
   end
-  axislegend(ax, position = :lt)
+  Legend(f[i+1, 2], _lines, names)
 end
 f
 ```

From 6667d4ee382580c3829a0fe18640cdadda345e0e Mon Sep 17 00:00:00 2001
From: Chris Elrod <elrodc@gmail.com>
Date: Tue, 11 Jun 2024 17:50:29 -0400
Subject: [PATCH 14/19] fix plot

---
 benchmarks/ModelingToolkit/RCCircuit.jmd | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/benchmarks/ModelingToolkit/RCCircuit.jmd b/benchmarks/ModelingToolkit/RCCircuit.jmd
index 2786db94c..62218def3 100644
--- a/benchmarks/ModelingToolkit/RCCircuit.jmd
+++ b/benchmarks/ModelingToolkit/RCCircuit.jmd
@@ -86,25 +86,25 @@ times = fill((NaN,NaN,NaN), length(N), length(max_sizes));
 @time run_and_time!(ss_times, times, max_sizes, 1, 4); # precompile
 
 for (i, n) in enumerate(N)
-  run_and_time!(ss_times, times, max_sizes, i, n)
+  @time run_and_time!(ss_times, times, max_sizes, i, n)
 end
 
 f = Figure(size=(800,1200));
+ss_names = ["MTK", "JSIR-Scalar", "JSIR-Loop"];
 let ax = Axis(f[1, 1]; yscale = log10, xscale = log10, title="Structural Simplify Time")
   _lines = map(eachcol(ss_times)) do ts
     lines!(N, ts)
   end
-  names = ("MTK", "JSIR-Scalar", "JSIR-Loop")
-  Legend(f([1,2], _lines, names)
+  Legend(f[1,2], _lines, ss_names)
 end
+method_names = ["MTK", "JSIR - Scalar - Julia", "JSIR - Scalar - C", "JSIR - Scalar - LLVM", "JSIR - Loop - Julia", "JSIR - Loop - C", "JSIR - Loop - LLVM"];
 for (i, timecat) in enumerate(("ODEProblem + f!", "Run", "Solve"))
   title = timecat * " Time"
   ax = Axis(f[i+1, 1]; yscale = log10, xscale = log10, title)
-  names = ("MTK", "JSIR - Scalar - Julia", "JSIR - Scalar - C", "JSIR - Scalar - LLVM", "JSIR - Loop - Julia", "JSIR - Loop - C", "JSIR - Loop - LLVM")
   _lines = map(eachcol(times)) do ts
     lines!(N, getindex.(ts, i))
   end
-  Legend(f[i+1, 2], _lines, names)
+  Legend(f[i+1, 2], _lines, method_names)
 end
 f
 ```

From a3186b135867f6d937572c9b2a70c6e748600fd4 Mon Sep 17 00:00:00 2001
From: Chris Elrod <elrodc@gmail.com>
Date: Wed, 12 Jun 2024 12:26:56 -0400
Subject: [PATCH 15/19] add author

---
 benchmarks/ModelingToolkit/RCCircuit.jmd | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/benchmarks/ModelingToolkit/RCCircuit.jmd b/benchmarks/ModelingToolkit/RCCircuit.jmd
index 62218def3..f4f6b05fa 100644
--- a/benchmarks/ModelingToolkit/RCCircuit.jmd
+++ b/benchmarks/ModelingToolkit/RCCircuit.jmd
@@ -1,6 +1,6 @@
 ---
 title: RC Circuit
-author: Yingbo Ma, Chris Elrod
+author: Avinash Subramanian, Yingbo Ma, Chris Elrod
 ---
 
 When a model is defined using repeated components, JuliaSimCompiler is able to take advantage of this to scale efficiently by rerolling equations into loops. This option can be disabled by setting `loop=false`. Here, we build an RC circuit model with variable numbers of components to show scaling of compile and runtimes of MTK vs JuliaSimCompiler's three backends with and without loop rerolling.

From c9f3407a168338d7ef79640fbb886163f6453e70 Mon Sep 17 00:00:00 2001
From: Chris Elrod <elrodc@gmail.com>
Date: Thu, 13 Jun 2024 13:17:49 -0400
Subject: [PATCH 16/19] run OpenModelica benchmarks; currently appears to be
 broken

---
 benchmarks/ModelingToolkit/RCCircuit.jmd | 80 ++++++++++++------------
 1 file changed, 41 insertions(+), 39 deletions(-)

diff --git a/benchmarks/ModelingToolkit/RCCircuit.jmd b/benchmarks/ModelingToolkit/RCCircuit.jmd
index f4f6b05fa..10efd7a17 100644
--- a/benchmarks/ModelingToolkit/RCCircuit.jmd
+++ b/benchmarks/ModelingToolkit/RCCircuit.jmd
@@ -5,10 +5,11 @@ author: Avinash Subramanian, Yingbo Ma, Chris Elrod
 
 When a model is defined using repeated components, JuliaSimCompiler is able to take advantage of this to scale efficiently by rerolling equations into loops. This option can be disabled by setting `loop=false`. Here, we build an RC circuit model with variable numbers of components to show scaling of compile and runtimes of MTK vs JuliaSimCompiler's three backends with and without loop rerolling.
 ```julia
-using JuliaSimCompiler, ModelingToolkit, OrdinaryDiffEq, BenchmarkTools, ModelingToolkitStandardLibrary, CairoMakie
+using JuliaSimCompiler, ModelingToolkit, OrdinaryDiffEq, BenchmarkTools, ModelingToolkitStandardLibrary, OMJulia, CairoMakie
 using ModelingToolkitStandardLibrary.Blocks
 using ModelingToolkitStandardLibrary.Electrical
 
+# ModelingToolkit and JuliaSimCompiler
 const t = Blocks.t
 
 function build_system(n)
@@ -56,7 +57,21 @@ end
 const C = JuliaSimCompiler.CTarget();
 const LLVM = JuliaSimCompiler.llvm.LLVMTarget();
 
-function run_and_time!(ss_times, times, max_sizes, i, n)
+# OMJ
+const omod = OMJulia.OMCSession();
+OMJulia.sendExpression(omod, "getVersion()")
+const modelicafile = joinpath(@__DIR__, "RC_Circuit.mo")
+
+function time_open_modelica(n::Int)
+  totaltime = @elapsed res = begin
+    ModelicaSystem(omod, modelicafile, "RC_Circuit.Test.RC_Circuit_MTK_test_$n")
+    sendExpression(omod, "simulate(RC_Circuit.Test.RC_Circuit_MTK_test_$n)")
+  end
+  @assert res["messages"][1:11] == "LOG_SUCCESS"
+  return totaltime
+end
+
+function run_and_time_julia!(ss_times, times, max_sizes, i, n)
   sys, u0, ps = build_system(n);
   if n <= max_sizes[1]
     ss_times[i, 1] = @elapsed sys_mtk = structural_simplify(sys)
@@ -71,19 +86,32 @@ function run_and_time!(ss_times, times, max_sizes, i, n)
   n <= max_sizes[5] && ((times[i, 5], deeref) = compile_run_problem(sys_jsir_loop, u0, ps; duref = daeref))
   n <= max_sizes[6] && ((times[i, 6], daeref) = compile_run_problem(sys_jsir_loop, u0, ps; target=C, duref = daeref))
   n <= max_sizes[7] && ((times[i, 7], daeref) = compile_run_problem(sys_jsir_loop, u0, ps; target=LLVM, duref = daeref))
-  @show n, ss_times[i, :], times[i, :]
+  for j = 1:7
+    ss_time = j == 1 ? ss_times[i,1] : ss_times[i, 2 + (j >= 5)]
+    t_fode, t_run, t_solve = times[i,j]
+    total_times[i, j] = ss_time + t_fode + t_solve
+  end
 end
 
-N = vcat([5, 10, 20, 40, 60, 80, 160, 320, 480, 640, 800], 1_000:1_000:10_000);
+function run_and_time!(ss_times, times, max_sizes, i, n)
+  run_and_time_julia!(ss_times, times, max_sizes, i, n)
+  if n <= max_sizes[8]
+    total_times[i, 8] = time_open_modelica(n)
+  end 
+  @views println("n = $(n)\nstructural_simplify_times = $(ss_times[i,:])\ncomponent times = $(times[i, :])\ntotal times = $(total_times[i, :])")
+end
+
+N = [5, 10, 20, 40, 60, 80, 160, 320, 480, 640, 800, 1000, 2000, 3000, 4000, 5000, 6000, 7000, 8000, 9000, 10000, 20000];
 
 # max size we test per method
-max_sizes = [1000, 4_000, 10_000, 10_000, 10_000, 10_000, 10_000];
+max_sizes = [4_000, 8_000, 20_000, 20_000, 20_000, 20_000, 20_000, 9000];
 
 # NaN-initialize so Makie will ignore incomplete
 ss_times = fill(NaN, length(N), 3);
-times = fill((NaN,NaN,NaN), length(N), length(max_sizes));
+times = fill((NaN,NaN,NaN), length(N), length(max_sizes) - 1);
+total_times = fill(NaN, length(N), length(max_sizes));
 
-@time run_and_time!(ss_times, times, max_sizes, 1, 4); # precompile
+@time run_and_time_julia!(ss_times, times, max_sizes, 1, 4); # precompile
 
 for (i, n) in enumerate(N)
   @time run_and_time!(ss_times, times, max_sizes, i, n)
@@ -106,7 +134,13 @@ for (i, timecat) in enumerate(("ODEProblem + f!", "Run", "Solve"))
   end
   Legend(f[i+1, 2], _lines, method_names)
 end
+let method_names_m = vcat(method_names, "OpenModelica");
+  ax = Axis(f[5, 1]; yscale = log10, xscale = log10, title = "Total Time")
+  _lines = map(Base.Fix1(lines!, N), eachcol(total_times))
+  Legend(f[5, 2], _lines, method_names_m)
+end
 f
+OMJulia.quit(omod) 
 ```
 All three backends compiled more quickly with loops, but the C and LLVM backends are so much quicker to compile than the Julia backend that this made much less difference for them.
 The impact on runtime was more varied.
@@ -118,38 +152,6 @@ using SciMLBenchmarks
 SciMLBenchmarks.bench_footer(WEAVE_ARGS[:folder],WEAVE_ARGS[:file])
 ```
 
-## Time OpenModelica
-
-```julia
-using OMJulia, CSV, DataFrames
-mod = OMJulia.OMCSession();
-
-OMJulia.sendExpression(mod, "getVersion()")
-
-modelicafile = joinpath(@__DIR__, "RC_Circuit.mo")
-
-resultfile = joinpath(@__DIR__, "RC_Circuit_modelica_res.csv")
-
-@show "Start OpenModelica Timings"
-
-# TODO: Sync these test cases with the N used in MTK: 
-N_x = [5, 10, 20, 40, 60, 80, 160, 320, 480, 640, 800, 1000, 2000, 3000, 4000, 5000, 6000, 7000, 8000, 9000, 10000, 20000]
-
-for N in N_x
-    @show N
-    totaltime = @elapsed res = begin
-        ModelicaSystem(mod, modelicafile, "RC_Circuit.Test.RC_Circuit_MTK_test_$N")
-        sendExpression(mod, "simulate(RC_Circuit.Test.RC_Circuit_MTK_test_$N)")
-    end
-    #runtime = res["timeTotal"]
-    @assert res["messages"][1:11] == "LOG_SUCCESS"
-    @show totaltime
-    # total_times[i, 5] = totaltime TODO: Add to comparison table
-end
-
-OMJulia.quit(mod) 
-```
-
 Dymola requires a license server and thus cannot be hosted. This was run locally for the
 following times:
 

From f3d65eba3aab324b322ff225f44fc0aa846bae82 Mon Sep 17 00:00:00 2001
From: Chris Elrod <elrodc@gmail.com>
Date: Thu, 13 Jun 2024 18:57:09 -0400
Subject: [PATCH 17/19] Install Modelica library

---
 benchmarks/ModelingToolkit/RCCircuit.jmd | 1 +
 1 file changed, 1 insertion(+)

diff --git a/benchmarks/ModelingToolkit/RCCircuit.jmd b/benchmarks/ModelingToolkit/RCCircuit.jmd
index 10efd7a17..704f59b27 100644
--- a/benchmarks/ModelingToolkit/RCCircuit.jmd
+++ b/benchmarks/ModelingToolkit/RCCircuit.jmd
@@ -60,6 +60,7 @@ const LLVM = JuliaSimCompiler.llvm.LLVMTarget();
 # OMJ
 const omod = OMJulia.OMCSession();
 OMJulia.sendExpression(omod, "getVersion()")
+OMJulia.sendExpression(omod, "installPackage(Modelica)")
 const modelicafile = joinpath(@__DIR__, "RC_Circuit.mo")
 
 function time_open_modelica(n::Int)

From 368251860ad1b486691ab4c1ee6787d9cab63dd5 Mon Sep 17 00:00:00 2001
From: Christopher Rackauckas <accounts@chrisrackauckas.com>
Date: Thu, 13 Jun 2024 20:35:24 -0400
Subject: [PATCH 18/19] Update RCCircuit.jmd

---
 benchmarks/ModelingToolkit/RCCircuit.jmd | 1 +
 1 file changed, 1 insertion(+)

diff --git a/benchmarks/ModelingToolkit/RCCircuit.jmd b/benchmarks/ModelingToolkit/RCCircuit.jmd
index 704f59b27..f1edaee20 100644
--- a/benchmarks/ModelingToolkit/RCCircuit.jmd
+++ b/benchmarks/ModelingToolkit/RCCircuit.jmd
@@ -4,6 +4,7 @@ author: Avinash Subramanian, Yingbo Ma, Chris Elrod
 ---
 
 When a model is defined using repeated components, JuliaSimCompiler is able to take advantage of this to scale efficiently by rerolling equations into loops. This option can be disabled by setting `loop=false`. Here, we build an RC circuit model with variable numbers of components to show scaling of compile and runtimes of MTK vs JuliaSimCompiler's three backends with and without loop rerolling.
+
 ```julia
 using JuliaSimCompiler, ModelingToolkit, OrdinaryDiffEq, BenchmarkTools, ModelingToolkitStandardLibrary, OMJulia, CairoMakie
 using ModelingToolkitStandardLibrary.Blocks

From 1cf905b1c6950093e7cf211a2a313bd9a6e3a029 Mon Sep 17 00:00:00 2001
From: Christopher Rackauckas <accounts@chrisrackauckas.com>
Date: Fri, 14 Jun 2024 06:02:07 -0400
Subject: [PATCH 19/19] Update RCCircuit.jmd

---
 benchmarks/ModelingToolkit/RCCircuit.jmd | 1 +
 1 file changed, 1 insertion(+)

diff --git a/benchmarks/ModelingToolkit/RCCircuit.jmd b/benchmarks/ModelingToolkit/RCCircuit.jmd
index f1edaee20..3de00a82d 100644
--- a/benchmarks/ModelingToolkit/RCCircuit.jmd
+++ b/benchmarks/ModelingToolkit/RCCircuit.jmd
@@ -144,6 +144,7 @@ end
 f
 OMJulia.quit(omod) 
 ```
+
 All three backends compiled more quickly with loops, but the C and LLVM backends are so much quicker to compile than the Julia backend that this made much less difference for them.
 The impact on runtime was more varied.