Merge pull request #949 from SciML/ChrisRackauckas-patch-1

Add Modelica code to ThermalFluid Benchmark

.buildkite/build_benchmark.sh

@@ -3,6 +3,7 @@
 set -euo pipefail
 
 JULIAHUBREGISTRY_BENCHMARK_TARGETS=(benchmarks/ModelingToolkit/)
+OPENMODELICA_BENCHMARK_TARGETS=(benchmarks/ModelingToolkit/)
 
 if [[ "${JULIAHUBREGISTRY_BENCHMARK_TARGETS[*]}" =~ "${1}" ]]; then
 	echo "--- :julia: Adding JuliaHubRegistry"

.buildkite/launch_benchmarks.yml.signature

@@ -1,2 +1 @@
-Salted__�1�����6�H�	�z)��-h�}[$>�W���^������f�ώfvc
-L��x�9�n�S�m��Ý(e^��Fg��u�����qD��C
+Salted__���O�������*��Ԏ�J�>X:p��L`�8��x��*)ޡ�Y� �r�wj�O�C��U�~^�¥M�VɃ�S�v�E%�

.buildkite/run_benchmark.yml

@@ -19,8 +19,10 @@ steps:
       - JuliaCI/julia#v1:
           version: "1.10"
       - staticfloat/sandbox:
-          rootfs_url: "https://github.com/ven-k/Placeholder/releases/download/v0.23.0/aws_uploader.x86_64.tar.gz"
-          rootfs_treehash: "d46b35aa927024de8729d63fde18442a0a590e62"
+          rootfs_url: "https://github.com/thazhemadam/openmodelica-rootfs-image/releases/download/v1.23.0/rootfs-openmodelica-v1.23.0.amd64.tar.gz"
+          rootfs_treehash: "82970243dc4f188e599a976abc20951f4aba2912"
+          uid: 1000
+          gid: 1000
           workspaces:
             # Include the julia we just downloaded
             - "/cache/julia-buildkite-plugin:/cache/julia-buildkite-plugin"

benchmarks/ModelingToolkit/DymolaTiming.mos

@@ -0,0 +1,53 @@
+import Modelica.Utilities.Streams.print;
+import DhnControl.TimingFunctions.tic;
+dir = "ThermalFluid_Benchmark"
+// Load the test model file
+openModel(dir + "DhnControl.mo");
+testFolder = "DhnControl.Test.";
+testModels = {
+    "test_preinsulated_470_5",
+    "test_preinsulated_470_10",
+    "test_preinsulated_470_20",
+    "test_preinsulated_470_40",
+    "test_preinsulated_470_60",
+    "test_preinsulated_470_80",
+    "test_preinsulated_470_160",
+    "test_preinsulated_470_320",
+    "test_preinsulated_470_480",
+    "test_preinsulated_470_640",
+    "test_preinsulated_470_800",
+    "test_preinsulated_470_960",
+    "test_preinsulated_470_1280"
+};
+
+outputFile = dir + "benchmark_results.txt"
+
+// Clear the output file at the beginning
+print("", fileName=outputFile);
+// Loop over each test model
+for i in 1:size(testModels, 1) loop
+  testName = testFolder + testModels[i];
+
+  // Start the timer
+  startTime := tic();
+
+  // Translate the model
+  translateModel(testName);
+
+  translateEndTime := tic();
+  // Simulate the model
+  simulateModel(testName);
+
+  // Stop the timer
+  endTime := tic();
+
+  // Calculate the elapsed time
+  translationTime := translateEndTime - startTime;
+  simulationTime := endTime - startTime;
+
+  // Prepare the output string
+  outputString = "Test: " + testModels[i] + ", " + String(translationTime) + ", " + String(simulationTime);
+
+  // Print the test name and simulation time to the output file
+  print(outputString, fileName=outputFile);
+end for;

benchmarks/ModelingToolkit/Project.toml

@@ -5,18 +5,24 @@ JuliaSimCompiler = "0.1.9"
 JuliaSimCompilerRuntime = "1.0.0"
 LinearSolve = "2.30.0"
 ModelingToolkit = "9.15.0"
+OMJulia = "0.3.1"
+OrdinaryDiffEq = "6.80.0"
 Polynomials = "4.0.8"
 PreferenceTools = "0.1.2"
 SciMLBenchmarks = "0.1.3"
 Symbolics = "5.28.0"
 
 [deps]
 BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
+CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
 CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
+DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
 JuliaSimCompiler = "8391cb6b-4921-5777-4e45-fd9aab8cb88d"
 JuliaSimCompilerRuntime = "9cbdfd5a-25c0-4dde-8b55-206f91b28bd9"
 LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
 ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
+OMJulia = "0f4fe800-344e-11e9-2949-fb537ad918e1"
+OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 Polynomials = "f27b6e38-b328-58d1-80ce-0feddd5e7a45"
 PreferenceTools = "ba661fbb-e901-4445-b070-854aec6bfbc5"
 SciMLBenchmarks = "31c91b34-3c75-11e9-0341-95557aab0344"

benchmarks/ModelingToolkit/ThermalFluid.jmd

@@ -1,16 +1,22 @@
 ---
 title: Thermal Fluid ODE Compilation and Perf
-author: Libor Kudela, Yingbo Ma, Chris Elrod
+author: Libor Kudela, Yingbo Ma, Chris Elrod, Chris Rackauckas
 ---
 
-This is a 1D advection-diffusion-source PDE that uses a second order upwind scheme. 
+This is a 1D advection-diffusion-source PDE that uses a second order upwind scheme.
+
 ```julia
 using Pkg
 # Rev fixes precompilation https://github.com/hzgzh/XSteam.jl/pull/2
 Pkg.add(Pkg.PackageSpec(;name="XSteam", rev="f2a1c589054cfd6bba307985a3a534b6f5a1863b"))
 
 using ModelingToolkit, JuliaSimCompiler, Symbolics, XSteam, Polynomials, BenchmarkTools, CairoMakie
+using OMJulia, OrdinaryDiffEq
+```
+
+## Setup Julia Code
 
+```julia
 #          o  o  o  o  o  o  o < heat capacitors
 #          |  |  |  |  |  |  | < heat conductors
 #          o  o  o  o  o  o  o
@@ -53,7 +59,7 @@ end
 #Taylor-aris dispersion model
 function Dxx_coeff(u, d, T)
   Re = abs(u) * d / kin_visc_T(T) + 0.1
-  if Re < 1000
+  if Re < 1000.0
     (d^2 / 4) * u^2 / 48 / 0.14e-6
   else
     d * u * (1.17e9 * Re^(-2.5) + 0.41)
@@ -63,9 +69,9 @@ end
 @register_symbolic Nusselt(Re, Pr, f)
 #Nusselt number model
 function Nusselt(Re, Pr, f)
-  if Re <= 2300
+  if Re <= 2300.0
     3.66
-  elseif Re <= 3100
+  elseif Re <= 3100.0
     3.5239 * (Re / 1000)^4 - 45.158 * (Re / 1000)^3 + 212.13 * (Re / 1000)^2 - 427.45 * (Re / 1000) + 316.08
   else
     f / 8 * ((Re - 1000) * Pr) / (1 + 12.7 * (f / 8)^(1 / 2) * (Pr^(2 / 3) - 1))
@@ -242,8 +248,6 @@ function TestBenchPreinsulated(; name, L=1.0, dn=0.05, t_layer=[0.0056, 0.013],
   compose(ODESystem(eqs, t, [], []; name=name), subs)
 end
 
-
-
 function build_system(fsys, N)
   N >= 4 || throw("Problem sizes smaller than 4 not supported; received $N.")
   @named testbench = TestBenchPreinsulated(; L=470, N, dn=0.3127, t_layer=[0.0056, 0.058])
@@ -278,6 +282,33 @@ function build_run_problem(fsys, N; target=JuliaSimCompiler.JuliaTarget())
   t_ss, t_fode, t_run
 end
 
+function test_speed(fsys, N; solver=FBDF(), target=JuliaSimCompiler.JuliaTarget())
+    tspan = (0.0, 19*3600)
+    t_total = @elapsed begin
+        @named testbench = TestBenchPreinsulated(L=470, N=N, dn=0.3127, t_layer=[0.0056, 0.058])
+        sys = structural_simplify(fsys(testbench))
+        if target === JuliaSimCompiler.JuliaTarget()
+            prob = ODEProblem(sys, [], tspan, sparse=true)
+        else
+            prob = ODEProblem(sys, target, [], tspan, sparse=true)
+        end
+        prob.u0 .= 12.0
+        solve(prob, solver, reltol=1e-6, abstol=1e-6, saveat=100);
+    end
+end
+```
+
+```julia
+N_x = [5, 10, 20, 40, 60, 80, 160, 320, 480, 640, 800, 960, 1280];
+N_states = 4 .* N_x; # x-axis for plots
+ss_times = Matrix{Float64}(undef, length(N_x), 2);
+times = Matrix{NTuple{2,Float64}}(undef, length(N_x), 4);
+total_times = Matrix{Float64}(undef, length(N_x), 6);
+```
+
+## Time Julia
+
+```julia
 using JuliaSimCompiler: IRSystem
 const MTKSystem = identity
 const CBackend = JuliaSimCompiler.CTarget();
@@ -288,41 +319,117 @@ const LLVMBackend = JuliaSimCompiler.llvm.LLVMTarget();
 @time build_run_problem(IRSystem, 4; target=CBackend)
 @time build_run_problem(IRSystem, 4; target=LLVMBackend)
 
-N_x = [5, 10, 20, 40, 60, 80, 160, 320, 480, 640, 800, 960];
-N_states = 4 .* N_x; # x-axis for plots
+@show "Start Julia Timings"
+
 
-ss_times = Matrix{Float64}(undef, length(N_x), 2);
-times = Matrix{NTuple{2,Float64}}(undef, length(N_x), 4);
 for (i, N_x_i) in enumerate(N_x)
+  @show i
   ss_times[i, 1], sys_mtk = build_system(MTKSystem, N_x_i)
   ss_times[i, 2], sys_jsir = build_system(IRSystem, N_x_i)
   times[i, 1] = compile_run_problem(sys_mtk)
   times[i, 2] = compile_run_problem(sys_jsir)
   times[i, 3] = compile_run_problem(sys_jsir, target=CBackend)
   times[i, 4] = compile_run_problem(sys_jsir, target=LLVMBackend)
-  @show N_x_i, ss_times[i, :], times[i, :]
+
+  if N_x_i >= 480
+    total_times[i, 1] = NaN
+  else
+    total_times[i, 1] = test_speed(MTKSystem, N_x_i)
+  end
+  total_times[i, 2] = test_speed(IRSystem, N_x_i)
+  total_times[i, 3] = test_speed(IRSystem, N_x_i, target=CBackend)
+  total_times[i, 4] = test_speed(IRSystem, N_x_i, target=LLVMBackend)
+
+  @show N_x_i, ss_times[i, :], times[i, :], total_times[i, :]
 end
+```
+
+## Time OpenModelica
+
+```julia
+using OMJulia, CSV, DataFrames
+mod = OMJulia.OMCSession();
+OMJulia.sendExpression(mod, "getVersion()")
+OMJulia.sendExpression(mod, "installPackage(Modelica)")
+modelicafile = "../../benchmarks/ModelingToolkit/DhnControl.mo"
+resultfile = "modelica_res.csv"
+
+@show "Start OpenModelica Timings"
+
+for i in 1:length(N_x)
+    N = N_x[i]
+    @show N
+    totaltime = @elapsed res = begin
+        @sync ModelicaSystem(mod, modelicafile, "DhnControl.Test.test_preinsulated_470_$N")
+        sendExpression(mod, "simulate(DhnControl.Test.test_preinsulated_470_$N)")
+    end
+    #runtime = res["timeTotal"]
+    @assert res["messages"][1:11] == "LOG_SUCCESS"
+    total_times[i, 5] = totaltime
+end
+
+OMJulia.quit(mod)
+total_times[:, 5]
+```
+
+## Time Dymola
 
-f = Figure();
-let ax = Axis(f[1, 1]; title="Structural Simplify Time")
-  for (j, label) in enumerate(("MTK", "JSIR"))
-    ts = @view(ss_times[:, j])
-    lines!(N_states, ts; label)
+Dymola requires a license server and thus cannot be hosted. This was run locally for the
+following times:
+
+```julia
+translation_and_total_times = [1.802 1.921
+                               1.78 1.846
+                               1.84 1.877
+                               2.028 2.075
+                               2.221 2.283
+                               2.409 2.496
+                               3.189 3.427
+                               4.758 5.577
+                               6.39 8.128
+                               8.052 11.026
+                               9.707 14.393
+                               11.411 17.752
+                               15.094 27.268]
+total_times[:, 6] = translation_and_total_times[1:length(N_x),2]
+```
+
+## Generate Final Plots
+
+```julia
+f = Figure(size = (800, 1200));
+let ax = Axis(f[1, 1]; yscale = log10, xscale = log10, title="Structural Simplify Time")
+  names = ["MTK", "JSIR"]
+  _lines = map(eachcol(ss_times)) do ts
+    lines!(N_states, ts)
   end
-  axislegend(ax, position = :lt)
+  Legend(f[1,2], _lines, names)
 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 - Julia", "JSIR - C", "JSIR - LLVM"))
-    ts = getindex.(@view(times[:, j]), i)
-    lines!(N_states, ts; label)
+  ax = Axis(f[i+1, 1]; yscale = log10, xscale = log10, title)
+  names = ["MTK", "JSIR - Julia", "JSIR - C", "JSIR - LLVM"]
+  _lines = map(eachcol(times)) do ts
+    lines!(N_states, getindex.(ts, i))
   end
-  axislegend(ax, position = :lt)
+  Legend(f[i+1,2], _lines, names)
 end
 f
 ```
 
+```julia
+f2 = Figure(size = (800, 400));
+title = "Total Time: Thermal Fluid Benchmark"
+ax = Axis(f2[1, 1]; yscale = log10, xscale = log10, title)
+names = ["MTK", "JSIR - Julia", "JSIR - C", "JSIR - LLVM", "OpenModelica", "Dymola"]
+_lines = map(enumerate(names)) do (j, label)
+    ts = @view(total_times[:, j])
+    lines!(N_states, ts)
+end
+Legend(f2[1,2], _lines, names)
+f2
+```
+
 ## Appendix
 
 ```julia, echo = false