Add stiff heating system benchmark

benchmarks/StiffODE/HeatingSystem_wpd.jmd

@@ -0,0 +1,364 @@
+---
+title: Heating System Work-Precision Diagrams
+author: Chris Rackauckas
+---
+
+# Heating System
+
+This benchmark tests a heating system model with temperature control and hysteresis behavior. The system models a distribution circuit with multiple heated units, each having its own temperature controller with hysteretic behavior. The system exhibits stiffness due to the different time scales of the heat transfer processes and the rapid switching behavior of the controllers.
+
+The model consists of:
+- N heated units with individual temperature controllers
+- A distribution circuit that supplies heat to all units  
+- External temperature variations (daily cycle)
+- Hysteretic control logic for each heating unit
+- Saturation functions for control outputs
+
+This represents a realistic engineering system with multiple time scales, making it a good test for stiff ODE solvers.
+
+```julia
+using OrdinaryDiffEq, DiffEqDevTools, Sundials, Plots, ODEInterfaceDiffEq, LSODA, LinearSolve
+using PreallocationTools
+using ProfileSVG, BenchmarkTools, Profile
+gr() # gr(fmt=:png)
+using LinearAlgebra, StaticArrays, RecursiveFactorization
+
+# System parameters
+const N = 40  # Number of heated units
+
+# Physical parameters
+const Cu = N == 1 ? [2e7] : (ones(N) .+ range(0,1.348,length=N))*1e7  # Heat capacity of heated units
+const Cd = 2e6*N  # Heat capacity of distribution circuit
+const Gh = 200    # Thermal conductance of heating elements
+const Gu = 150    # Thermal conductance of heated units to the atmosphere
+const Qmax = N*3000  # Maximum power output of heat generation unit
+const Teps = 0.5     # Threshold of heated unit temperature controllers
+const Td0 = 343.15   # Set point of distribution circuit temperature
+const Tu0 = 293.15   # Heated unit temperature set point
+const Kp = Qmax/4    # Proportional gain of heat generation unit temperature controller
+const a = 50         # Gain of the hysteresis function
+const b = 15         # Slope of the saturation function at the origin
+
+# Helper functions for parameter and variable organization
+function set_up_params(p)
+    Cu   = @view p[1:N]
+    Cd   = p[N+1]
+    Gh   = p[N+2]
+    Gu   = p[N+3]
+    Qmax = p[N+4]
+    Teps = p[N+5]
+    Td0  = p[N+6]
+    Tu0  = p[N+7]
+    Kp   = p[N+8]
+    a    = p[N+9]
+    b    = p[N+10]
+
+    Cu, Cd, Gh, Gu, Qmax, Teps, Td0, Tu0, Kp, a, b
+end
+
+function set_up_vars(u)
+    Td = u[1]
+    Tu = @view u[2:N+1]
+    x  = @view u[N+2:2N+1]
+
+    Td, Tu, x
+end
+
+function combine_vars!(du, dTd, dTu, dx)
+    du[1] = dTd
+    du[2:N+1] .= dTu[1:N]
+    du[N+2:2N+1] .= dx[1:N]
+
+    nothing
+end
+
+# Hysteresis function for temperature controllers
+hist(x, p, e = 1) = -(x + 0.5)*(x - 0.5) * x * 1/(0.0474)*e + p
+
+# Saturation function for control outputs
+sat(x, xmin, xmax) = tanh(2*(x-xmin)/(xmax-xmin)-1) * (xmax-xmin)/2 + (xmax+xmin)/2
+
+# Main ODE function for the heating system
+function heating(du, u, (p, _Qh, _Que), t)
+    (Td, Tu, x) = set_up_vars(u)
+    (dTd, dTu, dx) = set_up_vars(du)
+    (Cu, Cd, Gh, Gu, Qmax, Teps, Td0, Tu0, Kp, a, b) = set_up_params(p)
+    Qh = PreallocationTools.get_tmp(_Qh, du)
+    Que = PreallocationTools.get_tmp(_Que, du)
+
+    # External temperature with daily variation
+    Text = 278.15 + 8*sin(2π*t/86400)
+
+    # Heat transfer from distribution to units (with control)
+    @inbounds for i in 1:length(Qh)
+        Qh[i] = Gh * (Td - Tu[i]) * (sat(b * x[i], -0.5, 0.5) + 0.5)
+    end
+
+    # Heat loss from units to environment
+    @inbounds for i in 1:length(Que)
+        Que[i] = Gu * (Tu[i] - Text)
+    end
+
+    # Heat input to distribution circuit (with saturation)
+    Qd = sat(Kp*(Td0 - Td), 0, Qmax)
+
+    # Temperature dynamics
+    dTd = (Qd - sum(Qh))/Cd
+    @inbounds for i in 1:length(dTu)
+        dTu[i] = (Qh[i] - Que[i]) / Cu[i]
+    end
+
+    # Controller states (hysteretic behavior)
+    @inbounds for i in 1:length(dx)
+        dx[i] = a * hist(x[i], Tu0 - Tu[i], Teps)
+    end
+
+    combine_vars!(du, dTd, dTu, dx)
+end
+
+# Initial conditions and parameters
+u0 = [Td0, fill(Tu0, N)..., fill(-0.5, N)...]
+Qh0 = DiffCache(zeros(N))
+Que0 = DiffCache(zeros(N))
+p = [Cu..., Cd, Gh, Gu, Qmax, Teps, Td0, Tu0, Kp, a, b]
+tspan = (0., 50_000.)
+
+# Problem setup
+prob = ODEProblem(heating, u0, tspan, (p, Qh0, Que0))
+
+# Reference solution using a robust stiff solver
+sol = solve(prob, Rodas5P(), abstol=1e-12, reltol=1e-12)
+test_sol = TestSolution(sol)
+abstols = 1.0 ./ 10.0 .^ (4:11)
+reltols = 1.0 ./ 10.0 .^ (1:8);
+```
+
+```julia
+plot(sol, vars=1, title="Distribution Circuit Temperature")
+```
+
+```julia
+plot(sol, tspan=(0.0, 5000.0))
+```
+
+## High Tolerances
+
+This is the speed when you just want the answer.
+
+```julia
+abstols = 1.0 ./ 10.0 .^ (5:8)
+reltols = 1.0 ./ 10.0 .^ (1:4);
+setups = [Dict(:alg=>Rosenbrock23()),
+          Dict(:alg=>FBDF()),
+          Dict(:alg=>QNDF()),
+          Dict(:alg=>TRBDF2()),
+          Dict(:alg=>CVODE_BDF()),
+          Dict(:alg=>rodas()),
+          Dict(:alg=>radau()),
+          Dict(:alg=>lsoda()),
+          Dict(:alg=>RadauIIA5()),
+          ]
+wp = WorkPrecisionSet(prob,abstols,reltols,setups;verbose=false,
+                      save_everystep=false,appxsol=test_sol,maxiters=Int(1e5),numruns=10)
+plot(wp)
+```
+
+```julia
+wp = WorkPrecisionSet(prob,abstols,reltols,setups;dense = false,verbose = false,
+                      appxsol=test_sol,maxiters=Int(1e5),error_estimate=:l2,numruns=10)
+plot(wp)
+```
+
+```julia
+wp = WorkPrecisionSet(prob,abstols,reltols,setups;verbose=false,
+                      appxsol=test_sol,maxiters=Int(1e5),error_estimate=:L2,numruns=10)
+plot(wp)
+```
+
+```julia
+setups = [Dict(:alg=>Rosenbrock23()),
+          Dict(:alg=>Kvaerno3()),
+          Dict(:alg=>CVODE_BDF()),
+          Dict(:alg=>KenCarp4()),
+          Dict(:alg=>TRBDF2()),
+          Dict(:alg=>KenCarp3()),
+          Dict(:alg=>Rodas4()),
+          Dict(:alg=>lsoda()),
+          Dict(:alg=>radau())]
+wp = WorkPrecisionSet(prob,abstols,reltols,setups; verbose=false,
+                      save_everystep=false,appxsol=test_sol,maxiters=Int(1e5),numruns=10)
+plot(wp)
+```
+
+```julia
+wp = WorkPrecisionSet(prob,abstols,reltols,setups;dense = false,verbose = false,
+                      appxsol=test_sol,maxiters=Int(1e5),error_estimate=:l2,numruns=10)
+plot(wp)
+```
+
+```julia
+wp = WorkPrecisionSet(prob,abstols,reltols,setups; verbose=false,
+                      appxsol=test_sol,maxiters=Int(1e5),error_estimate=:L2,numruns=10)
+plot(wp)
+```
+
+```julia
+setups = [Dict(:alg=>Rosenbrock23()),
+          Dict(:alg=>KenCarp5()),
+          Dict(:alg=>KenCarp4()),
+          Dict(:alg=>KenCarp3()),
+          Dict(:alg=>ARKODE(order=5)),
+          Dict(:alg=>ARKODE()),
+          Dict(:alg=>ARKODE(order=3))]
+names = ["Rosenbrock23" "KenCarp5" "KenCarp4" "KenCarp3" "ARKODE5" "ARKODE4" "ARKODE3"]
+wp = WorkPrecisionSet(prob,abstols,reltols,setups; verbose=false,
+                      names=names,save_everystep=false,appxsol=test_sol,maxiters=Int(1e5),numruns=10)
+plot(wp)
+```
+
+```julia
+setups = [Dict(:alg=>Rosenbrock23()),
+          Dict(:alg=>TRBDF2()),
+          Dict(:alg=>ImplicitEulerExtrapolation()),
+          Dict(:alg=>ImplicitEulerBarycentricExtrapolation()),
+          Dict(:alg=>ImplicitHairerWannerExtrapolation()),
+          Dict(:alg=>ABDF2()),
+          Dict(:alg=>FBDF()),
+          Dict(:alg=>QNDF()),
+]
+wp = WorkPrecisionSet(prob,abstols,reltols,setups; verbose=false,
+                      save_everystep=false,appxsol=test_sol,maxiters=Int(1e5))
+plot(wp)
+```
+
+```julia
+setups = [Dict(:alg=>Rosenbrock23()),
+          Dict(:alg=>TRBDF2()),
+          Dict(:alg=>ImplicitEulerExtrapolation(linsolve = RFLUFactorization())),
+          Dict(:alg=>ImplicitEulerBarycentricExtrapolation(linsolve = RFLUFactorization())),
+          Dict(:alg=>ImplicitHairerWannerExtrapolation(linsolve = RFLUFactorization())),
+          Dict(:alg=>ABDF2()),
+          Dict(:alg=>FBDF()),
+          Dict(:alg=>QNDF()),
+]
+wp = WorkPrecisionSet(prob,abstols,reltols,setups; verbose=false,
+                      save_everystep=false,appxsol=test_sol,maxiters=Int(1e5))
+plot(wp)
+```
+
+### Low Tolerances
+
+This is the speed at lower tolerances, measuring what's good when accuracy is needed.
+
+```julia
+abstols = 1.0 ./ 10.0 .^ (7:13)
+reltols = 1.0 ./ 10.0 .^ (4:10)
+
+setups = [
+          Dict(:alg=>FBDF()),
+          Dict(:alg=>QNDF()),
+          Dict(:alg=>Rodas4P()),
+          Dict(:alg=>CVODE_BDF()),
+          Dict(:alg=>ddebdf()),
+          Dict(:alg=>Rodas4()),
+          Dict(:alg=>Rodas5P()),
+          Dict(:alg=>rodas()),
+          Dict(:alg=>radau()),
+          Dict(:alg=>lsoda())
+          ]
+wp = WorkPrecisionSet(prob,abstols,reltols,setups;verbose=false,
+                      save_everystep=false,appxsol=test_sol,maxiters=Int(1e5),numruns=10)
+plot(wp)
+```
+
+```julia
+wp = WorkPrecisionSet(prob,abstols,reltols,setups;verbose=false,
+                      dense=false,appxsol=test_sol,maxiters=Int(1e5),error_estimate=:l2,numruns=10)
+plot(wp)
+```
+
+```julia
+wp = WorkPrecisionSet(prob,abstols,reltols,setups;verbose=false,
+                      appxsol=test_sol,maxiters=Int(1e5),error_estimate=:L2,numruns=10)
+plot(wp)
+```
+
+```julia
+setups = [Dict(:alg=>GRK4A()),
+          Dict(:alg=>Rodas5()),
+          Dict(:alg=>Kvaerno4()),
+          Dict(:alg=>Kvaerno5()),
+          Dict(:alg=>CVODE_BDF()),
+          Dict(:alg=>KenCarp4()),
+          Dict(:alg=>KenCarp5()),
+          Dict(:alg=>Rodas4()),
+          Dict(:alg=>Rodas5P()),
+          Dict(:alg=>radau()),
+          Dict(:alg=>ImplicitEulerExtrapolation(min_order = 3)),
+          Dict(:alg=>ImplicitEulerBarycentricExtrapolation()),
+          Dict(:alg=>ImplicitHairerWannerExtrapolation()),
+          ]
+wp = WorkPrecisionSet(prob,abstols,reltols,setups; verbose=false,
+                      save_everystep=false,appxsol=test_sol,maxiters=Int(1e5),numruns=10)
+plot(wp)
+```
+
+```julia
+wp = WorkPrecisionSet(prob,abstols,reltols,setups;verbose=false,
+                      dense=false,appxsol=test_sol,maxiters=Int(1e5),error_estimate=:l2,numruns=10)
+plot(wp)
+```
+
+```julia
+wp = WorkPrecisionSet(prob,abstols,reltols,setups; verbose=false,
+                      appxsol=test_sol,maxiters=Int(1e5),error_estimate=:L2,numruns=10)
+plot(wp)
+```
+
+Multithreading benchmarks with Parallel Extrapolation Methods
+
+```julia
+#Setting BLAS to one thread to measure gains
+LinearAlgebra.BLAS.set_num_threads(1)
+
+abstols = 1.0 ./ 10.0 .^ (11:13)
+reltols = 1.0 ./ 10.0 .^ (8:10)
+
+setups = [
+            Dict(:alg=>CVODE_BDF()),
+            Dict(:alg=>KenCarp4()),
+            Dict(:alg=>Rodas4()),
+            Dict(:alg=>Rodas5()),
+            Dict(:alg=>Rodas5P()),
+            Dict(:alg=>QNDF()),
+            Dict(:alg=>lsoda()),
+            Dict(:alg=>radau()),
+            Dict(:alg=>seulex()),
+            Dict(:alg=>ImplicitEulerExtrapolation(min_order = 5, init_order = 3,threading = OrdinaryDiffEq.PolyesterThreads())),
+            Dict(:alg=>ImplicitEulerExtrapolation(min_order = 5, init_order = 3,threading = false)),
+            Dict(:alg=>ImplicitEulerBarycentricExtrapolation(min_order = 5, threading = OrdinaryDiffEq.PolyesterThreads())),
+            Dict(:alg=>ImplicitEulerBarycentricExtrapolation(min_order = 5, threading = false)),
+            Dict(:alg=>ImplicitHairerWannerExtrapolation(threading = OrdinaryDiffEq.PolyesterThreads())),
+            Dict(:alg=>ImplicitHairerWannerExtrapolation(threading = false)),
+            ]
+
+solnames = ["CVODE_BDF","KenCarp4","Rodas4","Rodas5","Rodas5P","QNDF","lsoda","radau","seulex","ImplEulerExtpl (threaded)", "ImplEulerExtpl (non-threaded)",
+            "ImplEulerBaryExtpl (threaded)","ImplEulerBaryExtpl (non-threaded)","ImplHWExtpl (threaded)","ImplHWExtpl (non-threaded)"]
+
+wp = WorkPrecisionSet(prob,abstols,reltols,setups; verbose=false,
+                    names = solnames,save_everystep=false,appxsol=test_sol,maxiters=Int(1e5),numruns=10)
+
+plot(wp, title = "Implicit Methods: Heating System",legend=:outertopleft,size = (1000,500),
+     xticks = 10.0 .^ (-15:1:1),
+     yticks = 10.0 .^ (-6:0.3:5),
+     bottom_margin= 5Plots.mm)
+```
+
+### Conclusion
+
+This heating system benchmark tests stiff ODE solvers on a realistic engineering system with hysteretic control behavior, multiple time scales, and complex thermal dynamics. The switching behavior of the controllers creates challenges for adaptive time stepping algorithms.
+
+```julia, echo = false
+using SciMLBenchmarks
+SciMLBenchmarks.bench_footer(WEAVE_ARGS[:folder],WEAVE_ARGS[:file])
+```

benchmarks/StiffODE/Manifest.toml

@@ -1,8 +1,8 @@
 # This file is machine-generated - editing it directly is not advised
 
-julia_version = "1.10.9"
+julia_version = "1.10.10"
 manifest_format = "2.0"
-project_hash = "3916e167d6740e94237c780dce42c34ea9dc912e"
+project_hash = "051150c3c21cedc0c14e999dae64db139a481a27"
 
 [[deps.ADTypes]]
 git-tree-sha1 = "be7ae030256b8ef14a441726c4c37766b90b93a3"
@@ -1803,7 +1803,7 @@ version = "0.3.23+4"
 [[deps.OpenLibm_jll]]
 deps = ["Artifacts", "Libdl"]
 uuid = "05823500-19ac-5b8b-9628-191a04bc5112"
-version = "0.8.1+4"
+version = "0.8.5+0"
 
 [[deps.OpenMPI_jll]]
 deps = ["Artifacts", "CompilerSupportLibraries_jll", "Hwloc_jll", "JLLWrappers", "LazyArtifacts", "Libdl", "MPIPreferences", "TOML", "Zlib_jll"]

benchmarks/StiffODE/Project.toml

@@ -14,6 +14,7 @@ ODEInterfaceDiffEq = "09606e27-ecf5-54fc-bb29-004bd9f985bf"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 ParameterizedFunctions = "65888b18-ceab-5e60-b2b9-181511a3b968"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+PreallocationTools = "d236fae5-4411-538c-8e31-a6e3d9e00b46"
 ProfileSVG = "132c30aa-f267-4189-9183-c8a63c7e05e6"
 RecursiveFactorization = "f2c3362d-daeb-58d1-803e-2bc74f2840b4"
 SciMLBenchmarks = "31c91b34-3c75-11e9-0341-95557aab0344"