diff --git a/benchmarks/Symbolics/BCR.jmd b/benchmarks/Symbolics/BCR.jmd new file mode 100644 index 000000000..ad229da89 --- /dev/null +++ b/benchmarks/Symbolics/BCR.jmd @@ -0,0 +1,244 @@ +--- +title: BCR Symbolic Jacobian +author: Aayush Sabharwal, Bowen Zhu, Chris Rackauckas +--- + +The following benchmark is of 1122 ODEs with 24388 terms that describe a stiff +chemical reaction network modeling the BCR signaling network from [Barua et +al.](https://doi.org/10.4049/jimmunol.1102003). We use +[`ReactionNetworkImporters`](https://github.com/isaacsas/ReactionNetworkImporters.jl) +to load the BioNetGen model files as a +[Catalyst](https://github.com/SciML/Catalyst.jl) model, and then use +[ModelingToolkit](https://github.com/SciML/ModelingToolkit.jl) to convert the +Catalyst network model to ODEs. + +The resultant large model is used to benchmark the time taken to compute a symbolic +jacobian, generate a function to calculate it and call the function. + + +```julia +using Catalyst, ReactionNetworkImporters, + TimerOutputs, LinearAlgebra, ModelingToolkit, Chairmarks, + LinearSolve, Symbolics, SymbolicUtils, SymbolicUtils.Code, SparseArrays, CairoMakie, + PrettyTables + +using Pkg +Pkg.add(Pkg.PackageSpec(; name = "SymbolicUtils", rev="9d52810c547276a60c0aae502ffee7f93ff99143")) + +datadir = joinpath(dirname(pathof(ReactionNetworkImporters)),"../data/bcr") +const to = TimerOutput() +tf = 100000.0 + +# generate ModelingToolkit ODEs +prnbng = loadrxnetwork(BNGNetwork(), joinpath(datadir, "bcr.net")) +show(to) +rn = complete(prnbng.rn; split = false) +obs = [eq.lhs for eq in observed(rn)] +osys = convert(ODESystem, rn) + +rhs = [eq.rhs for eq in full_equations(osys)] +vars = unknowns(osys) +pars = parameters(osys) +``` + +The `sparsejacobian` function in Symbolics.jl is optimized for hashconsing and caching, and as such +performs very poorly without either of those features. We use the old implementation, optimized without +hashconsing, to benchmark performance without hashconsing and without caching to avoid biasing the results. + +```julia +include("old_sparse_jacobian.jl") +``` + +```julia +SymbolicUtils.ENABLE_HASHCONSING[] = false +@timeit to "Calculate jacobian - without hashconsing" jac_nohc = old_sparsejacobian(rhs, vars); +SymbolicUtils.ENABLE_HASHCONSING[] = true +Symbolics.toggle_derivative_caching!(false) +Symbolics.clear_derivative_caches!() +@timeit to "Calculate jacobian - hashconsing, without caching" jac_hc_nocache = old_sparsejacobian(rhs, vars); +Symbolics.toggle_derivative_caching!(true) +for fn in Symbolics.cached_derivative_functions() + stats = SymbolicUtils.get_stats(fn) + @assert stats.hits == stats.misses == 0 +end +Symbolics.clear_derivative_caches!() +@timeit to "Calculate jacobian - hashconsing and caching" jac_hc_cache = Symbolics.sparsejacobian(rhs, vars); + +@assert isequal(jac_nohc, jac_hc_nocache) +@assert isequal(jac_hc_nocache, jac_hc_cache) + +jac = jac_hc_cache +args = (vars, pars, ModelingToolkit.get_iv(osys)) +# out of place versions run into an error saying the expression is too large +# due to the `SymbolicUtils.Code.create_array` call. `iip_config` prevents it +# from trying to build the function. +kwargs = (; iip_config = (false, true), expression = Val{true}) +@timeit to "Build jacobian - no CSE" _, jac_nocse_iip = build_function(jac, args...; cse = false, kwargs...); +@timeit to "Build jacobian - CSE" _, jac_cse_iip = build_function(jac, args...; cse = true, kwargs...); + +jac_nocse_iip = eval(jac_nocse_iip) +jac_cse_iip = eval(jac_cse_iip) + +defs = defaults(osys) +u = Float64[Symbolics.fixpoint_sub(var, defs) for var in vars] +buffer_cse = similar(jac, Float64) +buffer_nocse = similar(jac, Float64) +p = Float64[Symbolics.fixpoint_sub(par, defs) for par in pars] +tt = 0.0 + +@timeit to "Compile jacobian - CSE" jac_cse_iip(buffer_cse, u, p, tt) +@timeit to "Compute jacobian - CSE" jac_cse_iip(buffer_cse, u, p, tt) + +@timeit to "Compile jacobian - no CSE" jac_nocse_iip(buffer_nocse, u, p, tt) +@timeit to "Compute jacobian - no CSE" jac_nocse_iip(buffer_nocse, u, p, tt) + +@assert isapprox(buffer_cse, buffer_nocse, rtol = 1e-10) + +show(to) +``` + +We'll also measure scaling. + + +```julia +function run_and_time_construct!(rhs, vars, pars, iv, N, i, jac_times, jac_allocs, build_times, functions) + outputs = rhs[1:N] + SymbolicUtils.ENABLE_HASHCONSING[] = false + jac_result = @be old_sparsejacobian(outputs, vars) + jac_times[1][i] = minimum(x -> x.time, jac_result.samples) + jac_allocs[1][i] = minimum(x -> x.bytes, jac_result.samples) + + SymbolicUtils.ENABLE_HASHCONSING[] = true + jac_result = @be (Symbolics.clear_derivative_caches!(); Symbolics.sparsejacobian(outputs, vars)) + jac_times[2][i] = minimum(x -> x.time, jac_result.samples) + jac_allocs[2][i] = minimum(x -> x.bytes, jac_result.samples) + + jac = Symbolics.sparsejacobian(outputs, vars) + args = (vars, pars, iv) + kwargs = (; iip_config = (false, true), expression = Val{true}) + + build_result = @be build_function(jac, args...; cse = false, kwargs...); + build_times[1][i] = minimum(x -> x.time, build_result.samples) + jacfn_nocse = eval(build_function(jac, args...; cse = false, kwargs...)[2]) + + build_result = @be build_function(jac, args...; cse = true, kwargs...); + build_times[2][i] = minimum(x -> x.time, build_result.samples) + jacfn_cse = eval(build_function(jac, args...; cse = true, kwargs...)[2]) + + functions[1][i] = let buffer = similar(jac, Float64), fn = jacfn_nocse + function nocse(u, p, t) + fn(buffer, u, p, t) + buffer + end + end + functions[2][i] = let buffer = similar(jac, Float64), fn = jacfn_cse + function cse(u, p, t) + fn(buffer, u, p, t) + buffer + end + end + + return nothing +end + +function run_and_time_call!(i, u, p, tt, functions, first_call_times, second_call_times) + jacfn_nocse = functions[1][i] + jacfn_cse = functions[2][i] + + call_result = @timed jacfn_nocse(u, p, tt) + first_call_times[1][i] = call_result.time + call_result = @timed jacfn_cse(u, p, tt) + first_call_times[2][i] = call_result.time + + call_result = @be jacfn_nocse(u, p, tt) + second_call_times[1][i] = minimum(x -> x.time, call_result.samples) + call_result = @be jacfn_cse(u, p, tt) + second_call_times[2][i] = minimum(x -> x.time, call_result.samples) +end +``` + +# Run benchmark + +```julia +Chairmarks.DEFAULTS.seconds = 15.0 +N = [10, 20, 40, 80, 160, 320] +jacobian_times = [zeros(Float64, length(N)), zeros(Float64, length(N))] +jacobian_allocs = copy.(jacobian_times) +functions = [Vector{Any}(undef, length(N)), Vector{Any}(undef, length(N))] +# [without_cse_times, with_cse_times] +build_times = copy.(jacobian_times) +first_call_times = copy.(jacobian_times) +second_call_times = copy.(jacobian_times) + +iv = ModelingToolkit.get_iv(osys) +run_and_time_construct!(rhs, vars, pars, iv, 10, 1, jacobian_times, jacobian_allocs, build_times, functions) +run_and_time_call!(1, u, p, tt, functions, first_call_times, second_call_times) +for (i, n) in enumerate(N) + @info i n + run_and_time_construct!(rhs, vars, pars, iv, n, i, jacobian_times, jacobian_allocs, build_times, functions) +end +for (i, n) in enumerate(N) + @info i n + run_and_time_call!(i, u, p, tt, functions, first_call_times, second_call_times) +end +``` + +# Plot figures + +```julia +tabledata = hcat(N, jacobian_times..., jacobian_allocs..., build_times..., first_call_times..., second_call_times...) +header = ["N", "Jacobian time (no hashconsing)", "Jacobian time (hashconsing)", "Jacobian allocated memory (no hashconsing) (B)", "Jacobian allocated memory (hashconsing) (B)", "`build_function` time (no CSE)", "`build_function` time (CSE)", "First call time (no CSE)", "First call time (CSE)", "Second call time (no CSE)", "Second call time (CSE)"] +pretty_table(tabledata; header, backend = Val(:html)) +``` + +```julia +f = Figure(size = (750, 400)) +titles = [ + "Jacobian symbolic computation", "Jacobian symbolic computation", "Code generation", + "Numerical function compilation", "Numerical function evaluation"] +labels = ["Time (seconds)", "Allocated memory (bytes)", + "Time (seconds)", "Time (seconds)", "Time (seconds)"] +times = [jacobian_times, jacobian_allocs, build_times, first_call_times, second_call_times] +axes = Axis[] +for i in 1:2 + label = labels[i] + data = times[i] + ax = Axis(f[1, i], xscale = log10, yscale = log10, xlabel = "model size", + xlabelsize = 10, ylabel = label, ylabelsize = 10, xticks = N, + title = titles[i], titlesize = 12, xticklabelsize = 10, yticklabelsize = 10) + push!(axes, ax) + l1 = scatterlines!(ax, N, data[1], label = "without hashconsing") + l2 = scatterlines!(ax, N, data[2], label = "with hashconsing") +end +Legend(f[1, 3], axes[1], "Methods", tellwidth = false, labelsize = 12, titlesize = 15) +axes2 = Axis[] +# make equal y-axis unit length +mn3, mx3 = extrema(reduce(vcat, times[3])) +xn3 = log10(mx3 / mn3) +mn4, mx4 = extrema(reduce(vcat, times[4])) +xn4 = log10(mx4 / mn4) +mn5, mx5 = extrema(reduce(vcat, times[5])) +xn5 = log10(mx5 / mn5) +xn = max(xn3, xn4, xn5) +xn += 0.2 +hxn = xn / 2 +hxn3 = (log10(mx3) + log10(mn3)) / 2 +hxn4 = (log10(mx4) + log10(mn4)) / 2 +hxn5 = (log10(mx5) + log10(mn5)) / 2 +ylims = [(exp10(hxn3 - hxn), exp10(hxn3 + hxn)), (exp10(hxn4 - hxn), exp10(hxn4 + hxn)), + (exp10(hxn5 - hxn), exp10(hxn5 + hxn))] +for i in 1:3 + ir = i + 2 + label = labels[ir] + data = times[ir] + ax = Axis(f[2, i], xscale = log10, yscale = log10, xlabel = "model size", + xlabelsize = 10, ylabel = label, ylabelsize = 10, xticks = N, + title = titles[ir], titlesize = 12, xticklabelsize = 10, yticklabelsize = 10) + ylims!(ax, ylims[i]...) + push!(axes2, ax) + l1 = scatterlines!(ax, N, data[1], label = "without hashconsing") + l2 = scatterlines!(ax, N, data[2], label = "with hashconsing") +end +save("bcr.pdf", f) +f +``` diff --git a/benchmarks/Symbolics/Project.toml b/benchmarks/Symbolics/Project.toml new file mode 100644 index 000000000..f434b6b74 --- /dev/null +++ b/benchmarks/Symbolics/Project.toml @@ -0,0 +1,25 @@ +[deps] +BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" +CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0" +Catalyst = "479239e8-5488-4da2-87a7-35f2df7eef83" +Chairmarks = "0ca39b1e-fe0b-4e98-acfc-b1656634c4de" +DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" +LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae" +ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78" +OrderedCollections = "bac558e1-5e72-5ebc-8fee-abe8a469f55d" +OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" +Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" +Polynomials = "f27b6e38-b328-58d1-80ce-0feddd5e7a45" +PrettyTables = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d" +ReactionNetworkImporters = "b4db0fb7-de2a-5028-82bf-5021f5cfa881" +SciMLBenchmarks = "31c91b34-3c75-11e9-0341-95557aab0344" +SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" +SymbolicUtils = "d1185830-fcd6-423d-90d6-eec64667417b" +Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7" +TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f" +XSteam = "95ff35a0-be81-11e9-2ca3-5b4e338e8476" + +[compat] +ModelingToolkit = "9" +SymbolicUtils = "3.26" +Symbolics = "6.34" diff --git a/benchmarks/Symbolics/ThermalFluid.jmd b/benchmarks/Symbolics/ThermalFluid.jmd new file mode 100644 index 000000000..22066a25a --- /dev/null +++ b/benchmarks/Symbolics/ThermalFluid.jmd @@ -0,0 +1,443 @@ +--- +title: Thermal Fluid Symbolic Jacobian Scaling +author: Aayush Sabharwal, Bowen Zhu, Chris Rackauckas +--- + +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")) +Pkg.add(Pkg.PackageSpec(; name = "SymbolicUtils", rev="9d52810c547276a60c0aae502ffee7f93ff99143")) + +using ModelingToolkit, Symbolics, SymbolicUtils, XSteam, Polynomials, CairoMakie, PrettyTables +using SparseArrays, Chairmarks, Statistics +using ModelingToolkit: t_nounits as t, D_nounits as D +``` + +## Setup Julia Code + +```julia +# o o o o o o o < heat capacitors +# | | | | | | | < heat conductors +# o o o o o o o +# | | | | | | | +#Source -> o--o--o--o--o--o--o -> Sink +# advection diff source PDE + +m_flow_source(t) = 2.75 +T_source(t) = (t > 12 * 3600) * 56.0 + 12.0 +# @register_symbolic m_flow_source(t) +# @register_symbolic T_source(t) + +#build polynomial liquid-water property only dependent on Temperature +p_l = 5 #bar +T_vec = collect(1:1:150); +@generated kin_visc_T(t) = :(Base.evalpoly(t, $(fit(T_vec, my_pT.(p_l, T_vec) ./ rho_pT.(p_l, T_vec), 5).coeffs...,))) +@generated lambda_T(t) = :(Base.evalpoly(t, $(fit(T_vec, tc_pT.(p_l, T_vec), 3).coeffs...,))) +@generated Pr_T(t) = :(Base.evalpoly(t, $(fit(T_vec, 1e3 * Cp_pT.(p_l, T_vec) .* my_pT.(p_l, T_vec) ./ tc_pT.(p_l, T_vec), 5).coeffs...,))) +@generated rho_T(t) = :(Base.evalpoly(t, $(fit(T_vec, rho_pT.(p_l, T_vec), 4).coeffs...,))) +@generated rhocp_T(t) = :(Base.evalpoly(t, $(fit(T_vec, 1000 * rho_pT.(p_l, T_vec) .* Cp_pT.(p_l, T_vec), 5).coeffs...,))) +# @register_symbolic kin_visc_T(t) +# @register_symbolic lambda_T(t) +# @register_symbolic Pr_T(t) +# @register_symbolic rho_T(t) +# @register_symbolic rhocp_T(t) + +@connector function FluidPort(; name, p=101325.0, m=0.0, T=0.0) + sts = @variables p(t) = p m(t) = m [connect = Flow] T(t) = T [connect = Stream] + ODESystem(Equation[], t, sts, []; name=name) +end + +@connector function VectorHeatPort(; name, N=100, T0=0.0, Q0=0.0) + sts = @variables (T(t))[1:N] = T0 (Q(t))[1:N] = Q0 [connect = Flow] + ODESystem(Equation[], t, [T; Q], []; name=name) +end + +@register_symbolic Dxx_coeff(u, d, T) +#Taylor-aris dispersion model +function Dxx_coeff(u, d, T) + Re = abs(u) * d / kin_visc_T(T) + 0.1 + if Re < 1000.0 + (d^2 / 4) * u^2 / 48 / 0.14e-6 + else + d * u * (1.17e9 * Re^(-2.5) + 0.41) + end +end + +@register_symbolic Nusselt(Re, Pr, f) +#Nusselt number model +function Nusselt(Re, Pr, f) + if Re <= 2300.0 + 3.66 + 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)) + end +end + +# @register_symbolic Churchill_f(Re, epsilon, d) +#Darcy weisbach friction factor +function Churchill_f(Re, epsilon, d) + theta_1 = (-2.457 * log(((7 / Re)^0.9) + (0.27 * (epsilon / d))))^16 + theta_2 = (37530 / Re)^16 + 8 * ((((8 / Re)^12) + (1 / ((theta_1 + theta_2)^1.5)))^(1 / 12)) +end + +function FluidRegion(; name, L=1.0, dn=0.05, N=100, T0=0.0, + lumped_T=50, diffusion=true, e=1e-4) + @named inlet = FluidPort() + @named outlet = FluidPort() + @named heatport = VectorHeatPort(N=N) + + dx = L / N + c = [-1 / 8, -3 / 8, -3 / 8] # advection stencil coefficients + A = pi * dn^2 / 4 + + p = @parameters C_shift = 0.0 Rw = 0.0 # stuff for latter + @variables begin + (T(t))[1:N] = fill(T0, N) + Twall(t)[1:N] = fill(T0, N) + (S(t))[1:N] = fill(T0, N) + (C(t))[1:N] = fill(1.0, N) + u(t) = 1e-6 + Re(t) = 1000.0 + Dxx(t) = 0.0 + Pr(t) = 1.0 + alpha(t) = 1.0 + f(t) = 1.0 + end + + sts = vcat(T, Twall, S, C, Num[u], Num[Re], Num[Dxx], Num[Pr], Num[alpha], Num[f]) + + eqs = Equation[ + Re ~ 0.1 + dn * abs(u) / kin_visc_T(lumped_T) + Pr ~ Pr_T(lumped_T) + f ~ Churchill_f(Re, e, dn) #Darcy-weisbach + alpha ~ Nusselt(Re, Pr, f) * lambda_T(lumped_T) / dn + Dxx ~ diffusion * Dxx_coeff(u, dn, lumped_T) + inlet.m ~ -outlet.m + inlet.p ~ outlet.p + inlet.T ~ instream(inlet.T) + outlet.T ~ T[N] + u ~ inlet.m / rho_T(inlet.T) / A + [C[i] ~ dx * A * rhocp_T(T[i]) for i in 1:N] + [S[i] ~ heatport.Q[i] for i in 1:N] + [Twall[i] ~ heatport.T[i] for i in 1:N] + + #source term + [S[i] ~ (1 / (1 / (alpha * dn * pi * dx) + abs(Rw / 1000))) * (Twall[i] - T[i]) for i in 1:N] + + #second order upwind + diffusion + source + D(T[1]) ~ u / dx * (inlet.T - T[1]) + Dxx * (T[2] - T[1]) / dx^2 + S[1] / (C[1] - C_shift) + D(T[2]) ~ u / dx * (c[1] * inlet.T - sum(c) * T[1] + c[2] * T[2] + c[3] * T[3]) + Dxx * (T[1] - 2 * T[2] + T[3]) / dx^2 + S[2] / (C[2] - C_shift) + [D(T[i]) ~ u / dx * (c[1] * T[i-2] - sum(c) * T[i-1] + c[2] * T[i] + c[3] * T[i+1]) + Dxx * (T[i-1] - 2 * T[i] + T[i+1]) / dx^2 + S[i] / (C[i] - C_shift) for i in 3:N-1] + D(T[N]) ~ u / dx * (T[N-1] - T[N]) + Dxx * (T[N-1] - T[N]) / dx^2 + S[N] / (C[N] - C_shift) + ] + + ODESystem(eqs, t, sts, p; systems=[inlet, outlet, heatport], name=name) +end + +# @register_symbolic Cn_circular_wall_inner(d, D, cp, ρ) +function Cn_circular_wall_inner(d, D, cp, ρ) + C = pi / 4 * (D^2 - d^2) * cp * ρ + return C / 2 +end + +# @register_symbolic Cn_circular_wall_outer(d, D, cp, ρ) +function Cn_circular_wall_outer(d, D, cp, ρ) + C = pi / 4 * (D^2 - d^2) * cp * ρ + return C / 2 +end + +# @register_symbolic Ke_circular_wall(d, D, λ) +function Ke_circular_wall(d, D, λ) + 2 * pi * λ / log(D / d) +end + +function CircularWallFEM(; name, L=100, N=10, d=0.05, t_layer=[0.002], + λ=[50], cp=[500], ρ=[7850], T0=0.0) + @named inner_heatport = VectorHeatPort(N=N) + @named outer_heatport = VectorHeatPort(N=N) + dx = L / N + Ne = length(t_layer) + Nn = Ne + 1 + dn = vcat(d, d .+ 2.0 .* cumsum(t_layer)) + Cn = zeros(Nn) + Cn[1:Ne] += Cn_circular_wall_inner.(dn[1:Ne], dn[2:Nn], cp, ρ) .* dx + Cn[2:Nn] += Cn_circular_wall_outer.(dn[1:Ne], dn[2:Nn], cp, ρ) .* dx + p = @parameters C_shift = 0.0 + Ke = Ke_circular_wall.(dn[1:Ne], dn[2:Nn], λ) .* dx + @variables begin + (Tn(t))[1:N, 1:Nn] = fill(T0, N, Nn) + (Qe(t))[1:N, 1:Ne] = fill(T0, N, Ne) + end + sts = [vec(Tn); vec(Qe)] + e0 = Equation[inner_heatport.T[i] ~ Tn[i, 1] for i in 1:N] + e1 = Equation[outer_heatport.T[i] ~ Tn[i, Nn] for i in 1:N] + e2 = Equation[Qe[i, j] ~ Ke[j] * (-Tn[i, j+1] + Tn[i, j]) for i in 1:N for j in 1:Ne] + e3 = Equation[D(Tn[i, 1]) * (Cn[1] + C_shift) ~ inner_heatport.Q[i] - Qe[i, 1] for i in 1:N] + e4 = Equation[D(Tn[i, j]) * Cn[j] ~ Qe[i, j-1] - Qe[i, j] for i in 1:N for j in 2:Nn-1] + e5 = Equation[D(Tn[i, Nn]) * Cn[Nn] ~ Qe[i, Ne] + outer_heatport.Q[i] for i in 1:N] + eqs = vcat(e0, e1, e2, e3, e4, e5) + ODESystem(eqs, t, sts, p; systems=[inner_heatport, outer_heatport], name=name) +end + +function CylindricalSurfaceConvection(; name, L=100, N=100, d=1.0, α=5.0) + dx = L / N + S = pi * d * dx + @named heatport = VectorHeatPort(N=N) + sts = @variables Tenv(t) = 0.0 + eqs = [ + Tenv ~ 18.0 + [heatport.Q[i] ~ α * S * (heatport.T[i] - Tenv) for i in 1:N] + ] + + ODESystem(eqs, t, sts, []; systems=[heatport], name=name) +end + +function PreinsulatedPipe(; name, L=100.0, N=100.0, dn=0.05, T0=0.0, t_layer=[0.004, 0.013], + λ=[50, 0.04], cp=[500, 1200], ρ=[7800, 40], α=5.0, + e=1e-4, lumped_T=50, diffusion=true) + @named inlet = FluidPort() + @named outlet = FluidPort() + @named fluid_region = FluidRegion(L=L, N=N, dn=dn, e=e, lumped_T=lumped_T, diffusion=diffusion) + @named shell = CircularWallFEM(L=L, N=N, d=dn, t_layer=t_layer, λ=λ, cp=cp, ρ=ρ) + @named surfconv = CylindricalSurfaceConvection(L=L, N=N, d=dn + 2.0 * sum(t_layer), α=α) + systems = [fluid_region, shell, inlet, outlet, surfconv] + eqs = [ + connect(fluid_region.inlet, inlet) + connect(fluid_region.outlet, outlet) + connect(fluid_region.heatport, shell.inner_heatport) + connect(shell.outer_heatport, surfconv.heatport) + ] + ODESystem(eqs, t, [], []; systems=systems, name=name) +end + +function Source(; name, p_feed=100000) + @named outlet = FluidPort() + sts = @variables m_flow(t) = 1e-6 + eqs = [ + m_flow ~ m_flow_source(t) + outlet.m ~ -m_flow + outlet.p ~ p_feed + outlet.T ~ T_source(t) + ] + compose(ODESystem(eqs, t, sts, []; name=name), [outlet]) +end + +function Sink(; name) + @named inlet = FluidPort() + eqs = [ + inlet.T ~ instream(inlet.T) + ] + compose(ODESystem(eqs, t, [], []; name=name), [inlet]) +end + +function TestBenchPreinsulated(; name, L=1.0, dn=0.05, t_layer=[0.0056, 0.013], N=100, diffusion=true, lumped_T=20) + @named pipe = PreinsulatedPipe(L=L, dn=dn, N=N, diffusion=diffusion, t_layer=t_layer, lumped_T=lumped_T) + @named source = Source() + @named sink = Sink() + subs = [source, pipe, sink] + eqs = [ + connect(source.outlet, pipe.inlet) + connect(pipe.outlet, sink.inlet) + ] + compose(ODESystem(eqs, t, [], []; name=name), subs) +end + +function call(fn, args...) + fn(args...) +end +``` + +The `sparsejacobian` function in Symbolics.jl is optimized for hashconsing and caching, and as such +performs very poorly without either of those features. We use the old implementation, optimized without +hashconsing, to benchmark performance without hashconsing and without caching to avoid biasing the results. + +```julia +include("old_sparse_jacobian.jl") + +function run_and_time_construction!(jacobian_times, jacobian_gctimes, jacobian_allocs, build_times, functions, i, N) + @mtkbuild sys = TestBenchPreinsulated(L=470, N=N, dn=0.3127, t_layer=[0.0056, 0.058]) + rhs = [eq.rhs for eq in full_equations(sys)] + dvs = unknowns(sys) + + @info "Built system" + SymbolicUtils.ENABLE_HASHCONSING[] = false + jac_result = @be old_sparsejacobian(rhs, dvs) + @info "No hashconsing benchmark" + jac_nocse = old_sparsejacobian(rhs, dvs) + @info "No hashconsing result" + jacobian_times[1][i] = mean(x -> x.time, jac_result.samples) + jacobian_gctimes[1][i] = mean(x -> x.time * x.gc_fraction, jac_result.samples) + jacobian_allocs[1][i] = mean(x -> x.bytes, jac_result.samples) + @info "times" jacobian_times[1][i] jacobian_gctimes[1][i] jacobian_allocs[1][i] + SymbolicUtils.ENABLE_HASHCONSING[] = true + jac_result = @be (Symbolics.clear_derivative_caches!(); Symbolics.sparsejacobian(rhs, dvs)) + @info "Hashconsing benchmark" + jac_cse = Symbolics.sparsejacobian(rhs, dvs) + @info "Hashconsing result" + jacobian_times[2][i] = mean(x -> x.time, jac_result.samples) + jacobian_gctimes[2][i] = mean(x -> x.time * x.gc_fraction, jac_result.samples) + jacobian_allocs[2][i] = mean(x -> x.bytes, jac_result.samples) + @info "times" jacobian_times[2][i] jacobian_gctimes[2][i] jacobian_allocs[2][i] + @assert isequal(jac_nocse, jac_cse) + jac = jac_cse + + ps = parameters(sys) + defs = defaults(sys) + u0 = Float64[Symbolics.fixpoint_sub(v, defs) for v in dvs] + p = Float64[Symbolics.fixpoint_sub(v, defs) for v in ps] + t0 = 0.0 + buffer_nocse = similar(jac, Float64) + buffer_nocse.nzval .= 0.0 + buffer_cse = similar(jac, Float64) + buffer_cse.nzval .= 0.0 + + f_jac_nocse = eval(build_function(jac, dvs, ps, t; iip_config = (false, true), expression = Val{true}, cse = false)[2]) + functions[1][i] = let buffer_nocse = buffer_nocse, u0 = u0, p = p, t0 = t0, f_jac_nocse = f_jac_nocse + function nocse() + f_jac_nocse(buffer_nocse, u0, p, t0) + buffer_nocse + end + end + @info "No CSE build_function result" + build_result_nocse = @be build_function(jac, dvs, ps, t; iip_config = (false, true), expression = Val{true}, cse = false) + @info "No CSE build_function benchmark" + build_times[1][i] = mean(x -> x.time, build_result_nocse.samples) + @info "build_function time" build_times[1][i] + + f_jac_cse = eval(build_function(jac, dvs, ps, t; iip_config = (false, true), expression = Val{true}, cse = true)[2]) + functions[2][i] = let buffer_cse = buffer_cse, u0 = u0, p = p, t0 = t0, f_jac_cse = f_jac_cse + function nocse() + f_jac_cse(buffer_cse, u0, p, t0) + buffer_cse + end + end + @info "CSE build_function result" + build_result_cse = @be build_function(jac, dvs, ps, t; iip_config = (false, true), expression = Val{true}, cse = true) + @info "CSE build_function benchmark" + build_times[2][i] = mean(x -> x.time, build_result_cse.samples) + @info "build_function time" build_times[2][i] + + return nothing +end + +function run_and_time_call!(functions, first_call_times, second_call_times, i) + fnocse = functions[1][i] + fcse = functions[2][i] + first_call_result_nocse = @timed fnocse() + first_call_times[1][i] = first_call_result_nocse.time + @info "First call time" first_call_times[1][i] + second_call_result_nocse = @be fnocse() + second_call_times[1][i] = mean(x -> x.time, second_call_result_nocse.samples) + @info "Runtime" second_call_times[1][i] + + first_call_result_cse = @timed fcse() + first_call_times[2][i] = first_call_result_cse.time + @info "First call time" first_call_times[2][i] + second_call_result_cse = @be fcse() + second_call_times[2][i] = mean(x -> x.time, second_call_result_cse.samples) + @info "Runtime" second_call_times[2][i] +end +``` + +```julia +N = [5, 10, 20, 40, 80, 160, 320]; +jacobian_times = [zeros(Float64, length(N)), zeros(Float64, length(N))] +functions = [Vector{Any}(undef, length(N)), Vector{Any}(undef, length(N))] +jacobian_gctimes = copy.(jacobian_times) +jacobian_allocs = copy.(jacobian_times) +# [without_cse_times, with_cse_times] +build_times = copy.(jacobian_times) +first_call_times = copy.(jacobian_times) +second_call_times = copy.(jacobian_times) +``` + +## Timings + +```julia +Chairmarks.DEFAULTS.seconds = 15.0 +# compile +run_and_time_construction!(jacobian_times, jacobian_gctimes, jacobian_allocs, build_times, functions, 1, 5) +run_and_time_call!(functions, first_call_times, second_call_times, 1) +for (i, n) in enumerate(N) + @info i n + @time run_and_time_construction!(jacobian_times, jacobian_gctimes, jacobian_allocs, build_times, functions, i, n) +end + +for (i, n) in enumerate(N) + @info i n + run_and_time_call!(functions, first_call_times, second_call_times, i) +end +``` + +## Results + +```julia +tabledata = hcat(N, jacobian_times..., jacobian_gctimes..., jacobian_allocs..., build_times..., first_call_times..., second_call_times...) +header = ["N", "Jacobian time (no hashconsing)", "Jacobian time (hashconsing)", "Jacobian GC time (no hashconsing)", "Jacobian GC time (hashconsing)", "Jacobian allocated memory (no hashconsing) (B)", "Jacobian allocated memory (hashconsing) (B)", "`build_function` time (no CSE)", "`build_function` time (CSE)", "First call time (no CSE)", "First call time (CSE)", "Second call time (no CSE)", "Second call time (CSE)"] +pretty_table(tabledata; header, backend = Val(:html)) +``` + +```julia +f = Figure(size = (750, 400)) +titles = [ + "Jacobian symbolic computation", "Jacobian symbolic computation", "Code generation", + "Numerical function compilation", "Numerical function evaluation"] +labels = ["Time (seconds)", "Allocated memory (bytes)", + "Time (seconds)", "Time (seconds)", "Time (seconds)"] +times = [jacobian_times, jacobian_allocs, build_times, first_call_times, second_call_times] +axes = Axis[] +for i in 1:2 + label = labels[i] + data = times[i] + ax = Axis(f[1, i], xscale = log10, yscale = log10, xlabel = "model size", + xlabelsize = 10, ylabel = label, ylabelsize = 10, xticks = N, + title = titles[i], titlesize = 12, xticklabelsize = 10, yticklabelsize = 10) + push!(axes, ax) + l1 = scatterlines!(ax, N, data[1], label = "without hashconsing") + l2 = scatterlines!(ax, N, data[2], label = "with hashconsing") +end +Legend(f[1, 3], axes[1], "Methods", tellwidth = false, labelsize = 12, titlesize = 15) +axes2 = Axis[] +# make equal y-axis unit length +mn3, mx3 = extrema(reduce(vcat, times[3])) +xn3 = log10(mx3 / mn3) +mn4, mx4 = extrema(reduce(vcat, times[4])) +xn4 = log10(mx4 / mn4) +mn5, mx5 = extrema(reduce(vcat, times[5])) +xn5 = log10(mx5 / mn5) +xn = max(xn3, xn4, xn5) +xn += 0.2 +hxn = xn / 2 +hxn3 = (log10(mx3) + log10(mn3)) / 2 +hxn4 = (log10(mx4) + log10(mn4)) / 2 +hxn5 = (log10(mx5) + log10(mn5)) / 2 +ylims = [(exp10(hxn3 - hxn), exp10(hxn3 + hxn)), (exp10(hxn4 - hxn), exp10(hxn4 + hxn)), + (exp10(hxn5 - hxn), exp10(hxn5 + hxn))] +for i in 1:3 + ir = i + 2 + label = labels[ir] + data = times[ir] + ax = Axis(f[2, i], xscale = log10, yscale = log10, xlabel = "model size", + xlabelsize = 10, ylabel = label, ylabelsize = 10, xticks = N, + title = titles[ir], titlesize = 12, xticklabelsize = 10, yticklabelsize = 10) + ylims!(ax, ylims[i]...) + push!(axes2, ax) + l1 = scatterlines!(ax, N, data[1], label = "without hashconsing") + l2 = scatterlines!(ax, N, data[2], label = "with hashconsing") +end +save("thermal_fluid.pdf", f) +f +``` + +## Appendix + +```julia, echo = false +using SciMLBenchmarks +SciMLBenchmarks.bench_footer(WEAVE_ARGS[:folder],WEAVE_ARGS[:file]) +``` diff --git a/benchmarks/Symbolics/old_sparse_jacobian.jl b/benchmarks/Symbolics/old_sparse_jacobian.jl new file mode 100644 index 000000000..9a009473c --- /dev/null +++ b/benchmarks/Symbolics/old_sparse_jacobian.jl @@ -0,0 +1,227 @@ +function old_sparsejacobian(ops::AbstractVector, vars::AbstractVector) + sp = Symbolics.jacobian_sparsity(ops, vars) + I,J,_ = findnz(sp) + + exprs = old_sparsejacobian_vals(ops, vars, I, J) + + sparse(I, J, exprs, length(ops), length(vars)) +end + +function old_sparsejacobian_vals(ops::AbstractVector, vars::AbstractVector, I::AbstractVector, J::AbstractVector; simplify::Bool=false, kwargs...) + exprs = Num[] + sizehint!(exprs, length(I)) + + for (i,j) in zip(I, J) + push!(exprs, Num(old_expand_derivatives(Differential(vars[j])(ops[i]), simplify; kwargs...))) + end + exprs +end + + +function old_expand_derivatives(O::SymbolicUtils.Symbolic, simplify=false; throw_no_derivative=false) + if iscall(O) && isa(operation(O), Differential) + arg = only(arguments(O)) + arg = old_expand_derivatives(arg, false; throw_no_derivative) + return old_executediff(operation(O), arg, simplify; throw_no_derivative) + elseif iscall(O) && isa(operation(O), Integral) + return operation(O)(old_expand_derivatives(arguments(O)[1]; throw_no_derivative)) + elseif !Symbolics.hasderiv(O) + return O + else + args = map(a->old_expand_derivatives(a, false; throw_no_derivative), arguments(O)) + O1 = operation(O)(args...) + return simplify ? SymbolicUtils.simplify(O1) : O1 + end +end +function old_expand_derivatives(n::Num, simplify=false; kwargs...) + Symbolics.wrap(old_expand_derivatives(Symbolics.value(n), simplify; kwargs...)) +end + +function old_occursin_info(x, expr, fail = true) + if SymbolicUtils.symtype(expr) <: AbstractArray + if fail + error("Differentiation with array expressions is not yet supported") + else + return occursin(x, expr) + end + end + + # Allow scalarized expressions + function is_scalar_indexed(ex) + (iscall(ex) && operation(ex) == getindex && !(SymbolicUtils.symtype(ex) <: AbstractArray)) || + (iscall(ex) && (SymbolicUtils.issym(operation(ex)) || iscall(operation(ex))) && + is_scalar_indexed(operation(ex))) + end + + # x[1] == x[1] but not x[2] + if is_scalar_indexed(x) && is_scalar_indexed(expr) && + isequal(first(arguments(x)), first(arguments(expr))) + return isequal(operation(x), operation(expr)) && + isequal(arguments(x), arguments(expr)) + end + + if is_scalar_indexed(x) && is_scalar_indexed(expr) && + !occursin(first(arguments(x)), first(arguments(expr))) + return false + end + + if is_scalar_indexed(expr) && !is_scalar_indexed(x) && !occursin(x, expr) + return false + end + + !iscall(expr) && return isequal(x, expr) + if isequal(x, expr) + true + else + args = map(a->old_occursin_info(x, a, operation(expr) !== getindex), arguments(expr)) + if all(_isfalse, args) + return false + end + Term{Real}(true, args) + end +end + +function old_occursin_info(x, expr::Sym, fail) + if SymbolicUtils.symtype(expr) <: AbstractArray && fail + error("Differentiation of expressions involving arrays and array variables is not yet supported.") + end + isequal(x, expr) +end + +_isfalse(occ::Bool) = occ === false +_isfalse(occ::SymbolicUtils.Symbolic) = iscall(occ) && _isfalse(operation(occ)) + +_iszero(x) = false +_isone(x) = false +_iszero(x::Number) = iszero(x) +_isone(x::Number) = isone(x) +_iszero(::SymbolicUtils.Symbolic) = false +_isone(::SymbolicUtils.Symbolic) = false +_iszero(x::Num) = _iszero(value(x))::Bool +_isone(x::Num) = _isone(value(x))::Bool + + +function old_executediff(D, arg, simplify=false; occurrences=nothing, throw_no_derivative=false) + if occurrences == nothing + occurrences = old_occursin_info(D.x, arg) + end + + _isfalse(occurrences) && return 0 + occurrences isa Bool && return 1 # means it's a `true` + + if !iscall(arg) + return D(arg) # Cannot expand + elseif (op = operation(arg); SymbolicUtils.issym(op)) + inner_args = arguments(arg) + if any(isequal(D.x), inner_args) + return D(arg) # base case if any argument is directly equal to the i.v. + else + return sum(inner_args, init=0) do a + return old_executediff(Differential(a), arg; throw_no_derivative) * + old_executediff(D, a; throw_no_derivative) + end + end + elseif op === getindex + inner_args = arguments(arguments(arg)[1]) + c = 0 + for a in inner_args + if isequal(a, D.x) + return D(arg) + else + c += Differential(a)(arg) * D(a) + end + end + return old_expand_derivatives(c) + elseif op === ifelse + args = arguments(arg) + O = op(args[1], + old_executediff(D, args[2], simplify; occurrences=arguments(occurrences)[2], throw_no_derivative), + old_executediff(D, args[3], simplify; occurrences=arguments(occurrences)[3], throw_no_derivative)) + return O + elseif isa(op, Differential) + # The recursive expand_derivatives was not able to remove + # a nested Differential. We can attempt to differentiate the + # inner expression wrt to the outer iv. And leave the + # unexpandable Differential outside. + if isequal(op.x, D.x) + return D(arg) + else + inner = old_executediff(D, arguments(arg)[1], false; throw_no_derivative) + # if the inner expression is not expandable either, return + if iscall(inner) && operation(inner) isa Differential + return D(arg) + else + # otherwise give the nested Differential another try + return old_executediff(op, inner, simplify; throw_no_derivative) + end + end + elseif isa(op, Integral) + if isa(op.domain.domain, Symbolics.AbstractInterval) + domain = op.domain.domain + a, b = Symbolics.DomainSets.endpoints(domain) + c = 0 + inner_function = arguments(arg)[1] + if iscall(value(a)) + t1 = SymbolicUtils.substitute(inner_function, Dict(op.domain.variables => value(a))) + t2 = D(a) + c -= t1*t2 + end + if iscall(value(b)) + t1 = SymbolicUtils.substitute(inner_function, Dict(op.domain.variables => value(b))) + t2 = D(b) + c += t1*t2 + end + inner = old_executediff(D, arguments(arg)[1]; throw_no_derivative) + c += op(inner) + return Symbolics.value(c) + end + end + + inner_args = arguments(arg) + l = length(inner_args) + exprs = [] + c = 0 + + for i in 1:l + t2 = old_executediff(D, inner_args[i],false; occurrences=arguments(occurrences)[i], throw_no_derivative) + + x = if _iszero(t2) + t2 + elseif _isone(t2) + d = Symbolics.derivative_idx(arg, i) + if d isa Symbolics.NoDeriv + throw_no_derivative && error((arg, i)) + D(arg) + else + d + end + else + t1 = Symbolics.derivative_idx(arg, i) + t1 = if t1 isa Symbolics.NoDeriv + throw_no_derivative && error((arg, i)) + D(arg) + else + t1 + end + t1 * t2 + end + + if _iszero(x) + continue + elseif x isa SymbolicUtils.Symbolic + push!(exprs, x) + else + c += x + end + end + + if isempty(exprs) + return c + elseif length(exprs) == 1 + term = (simplify ? SymbolicUtils.simplify(exprs[1]) : exprs[1]) + return _iszero(c) ? term : c + term + else + x = +((!_iszero(c) ? vcat(c, exprs) : exprs)...) + return simplify ? SymbolicUtils.simplify(x) : x + end +end