Merge pull request #954 from chriselrod/rccircuit

RC Circuit MTK and JSIR example.

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"
 

benchmarks/ModelingToolkit/RCCircuit.jmd

@@ -0,0 +1,174 @@
+---
+title: RC Circuit
+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
+using ModelingToolkitStandardLibrary.Electrical
+
+# ModelingToolkit and JuliaSimCompiler
+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(), duref=nothing)
+  tspan = (0.0, 10.0)
+  t0 = time()
+  prob = if target === JuliaSimCompiler.JuliaTarget()
+    ODEProblem(sys, u0, tspan, ps; sparse = true)
+  else
+    ODEProblem(sys, target, u0, tspan, ps; sparse = true)
+  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_solve = @elapsed solve(prob, Rodas5(autodiff = false))
+  (t_fode, t_run, t_solve), du
+end
+
+const C = JuliaSimCompiler.CTarget();
+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)
+  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)
+    times[i, 1], _ = compile_run_problem(sys_mtk, u0, ps)
+  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))
+  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))
+  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
+
+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 = [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) - 1);
+total_times = fill(NaN, length(N), length(max_sizes));
+
+@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)
+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
+  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)
+  _lines = map(eachcol(times)) do ts
+    lines!(N, getindex.(ts, i))
+  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.
+
+## Appendix
+
+```julia, echo = false
+using SciMLBenchmarks
+SciMLBenchmarks.bench_footer(WEAVE_ARGS[:folder],WEAVE_ARGS[:file])
+```
+
+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]
+```

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;

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
 ```