heating_equation.jl

using OrdinaryDiffEq, LabelledArrays

N = 40  # Number of heated units

Cu = N == 1 ? [2e7] : (ones(N) .+ range(0,1.348,length=N))*1e7 # "Heat capacity of heated units";
Cd = 2e6*N # "Heat capacity of distribution circuit";
Gh = 200 # "Thermal conductance of heating elements";
Gu = 150 # "Thermal conductance of heated units to the atmosphere";
Qmax = N*3000 # "Maximum power output of heat generation unit";
Teps = 0.5 # "Threshold of heated unit temperature controllers";
Td0 = 343.15 # "Set point of distribution circuit temperature";
Tu0 = 293.15 # "Heated unit temperature set point";
Kp = Qmax/4 # "Proportional gain of heat generation unit temperature controller";
a = 50 # "Gain of the histeresis function";
b = 15 # "Slope of the saturation function at the origin";

hist(x, p, e = 1) = -(x + 0.5)*(x - 0.5) *
                        x * 1/(0.0474)*e + p

sat(x, xmin, xmax) = tanh(2*(x-xmin)/(xmax-xmin)-1)*
                          (xmax-xmin)/2 +
                          (xmax+xmin)/2

function heating(du, u, (Cu, Cd, Gh, Gu, Qmax, Teps, Td0, Tu0, Kp, a, b, N), t)

    @views Td, Tu, x = u[1], u[2:N+1], u[N+2:2N+1]
    @views dTu, dx = du[2:N+1], du[N+2:2N+1]

    text = 278.15 + 8*sin(2π*t/86400)
    Qh = Gh .* (Td .- Tu) .* (sat.(b .* x, -0.5, 0.5) .+ 0.5)
    Que = Gu .* (Tu .- text)

    Qd = sat(Kp*(Td0 - Td),0, Qmax)
    du[1] = (Qd - sum(Qh))/Cd # dTd
    dTu .= (Qh .- Que) ./ Cu
    dx  .=  a .* hist.(x, Tu0 .- Tu, Teps)
    return nothing
end
tspan = (0., 50_000.)

p = (Cu=Cu, Cd=Cd, Gh=Gh, Gu=Gu, Qmax=Qmax, Teps=Teps, Td0=Td0, Tu0=Tu0, Kp=Kp, a=a, b=b, N=N)
u0 = [Td0, fill(Tu0, N)..., fill(-0.5, N)...]

prob = ODEProblem(heating, u0, tspan, p);
solve(prob, KenCarp4(), rtol=1e-7);
@time sol1 = solve(prob, KenCarp4(), rtol=1e-7);
@time sol1 = solve(prob, KenCarp4(), rtol=1e-7);
@time sol1 = solve(prob, KenCarp4(), rtol=1e-7);
solve(prob, RadauIIA5(), rtol=1e-7);
@time sol2 = solve(prob, RadauIIA5(), rtol=1e-7);
@time sol2 = solve(prob, RadauIIA5(), rtol=1e-7);
@time sol2 = solve(prob, RadauIIA5(), rtol=1e-7);
solve(prob, Rosenbrock23(), rtol=1e-7);
@time sol3 = solve(prob, Rosenbrock23(), rtol=1e-7);
@time sol3 = solve(prob, Rosenbrock23(), rtol=1e-7);
@time sol3 = solve(prob, Rosenbrock23(), rtol=1e-7);

#=
julia> solve(prob, KenCarp4(), rtol=1e-7);

julia> @time sol1 = solve(prob, KenCarp4(), rtol=1e-7);
  4.648946 seconds (4.05 M allocations: 1.762 GiB, 3.83% gc time)

julia> @time sol1 = solve(prob, KenCarp4(), rtol=1e-7);
  4.703236 seconds (4.05 M allocations: 1.762 GiB, 3.71% gc time)

julia> @time sol1 = solve(prob, KenCarp4(), rtol=1e-7);
  4.552906 seconds (4.05 M allocations: 1.762 GiB, 3.56% gc time)

julia> solve(prob, RadauIIA5(), rtol=1e-7);

julia> @time sol2 = solve(prob, RadauIIA5(), rtol=1e-7);
  5.216250 seconds (2.69 M allocations: 878.738 MiB, 5.67% gc time)

julia> @time sol2 = solve(prob, RadauIIA5(), rtol=1e-7);
  4.816831 seconds (2.69 M allocations: 878.738 MiB, 1.91% gc time)

julia> @time sol2 = solve(prob, RadauIIA5(), rtol=1e-7);
  4.816324 seconds (2.69 M allocations: 878.738 MiB, 1.73% gc time)

julia> solve(prob, Rosenbrock23(), rtol=1e-7);

julia> @time sol3 = solve(prob, Rosenbrock23(), rtol=1e-7);
  1.896335 seconds (1.58 M allocations: 938.699 MiB, 5.20% gc time)

julia> @time sol3 = solve(prob, Rosenbrock23(), rtol=1e-7);
  1.895587 seconds (1.58 M allocations: 938.699 MiB, 4.92% gc time)

julia> @time sol3 = solve(prob, Rosenbrock23(), rtol=1e-7);
  1.889226 seconds (1.58 M allocations: 938.699 MiB, 5.12% gc time)
=#

using Sundials, DiffEqDevTools, Plots
sol = solve(prob, CVODE_BDF(), reltol=1e-12, abstol=1e-12)
test_sol = TestSolution(sol)
abstols = @. 1/10^(6:2:12)
reltols = @. 1/10^(2:2:8)
setups = [
          Dict(:alg=>Rosenbrock23()),
          Dict(:alg=>TRBDF2()),
          Dict(:alg=>KenCarp4()),
          Dict(:alg=>CVODE_BDF()),
         ]
wp1 = WorkPrecisionSet(prob,abstols,reltols,setups; save_everystep=false,appxsol=test_sol,seconds=1,numruns=1,print_names=true)
wp2 = WorkPrecisionSet(prob,abstols,reltols,setups; save_everystep=false,appxsol=test_sol,seconds=1,numruns=1,print_names=true,error_estimate=:l2)
plot(wp1)
savefig("$(homedir())/Desktop/heating_bench_finial.svg")
plot(wp2)
savefig("$(homedir())/Desktop/heating_bench_l2.svg")