Cannot beat OpenModelica score with ModelingToolkit.jl and OrdinaryDiffEq.jl

[LiborKudela]

Hi,

I am playing with MTK and OpenModelica (OMC). I have written a model of pipe (1D advection-diffusion-source PDE) that uses second order upwind scheme. The pipe is surrounded with two layers of solid matter (heat equation) and is connected to the pipe via heatport (connector) that has vector variables in it. The pipe and the solid are discretized into N pieces. I have written the exact same model in Modelica.
The model agrees with experimental data so it should be correct.

I have compared the performance and I cannot beat the OMC performance with MTK + OrdinaryDiffEq for higher N (x-N, y-cputime).
(the number of states is 4*N btw)

I have tried quite a bunch of options but Tsit5 and QNDF look most promising.
The Tsit5() scales the same as OMC but is approximately. twice slower - > I am happy.
The QNDF(autodiff=false) is quite faster for small N but scales quite bad.
Why does QNDF scale this bad with N?
Is there some other solvers I can try, or some other solver options etc.?
(I also tried IDA() from Sundials.jl, but it was no good either.)
Should I even expect to beat OMC with Julia like this?

Thank you for any response

The MTK model:

using ModelingToolkit, OrdinaryDiffEq, Symbolics, IfElse, XSteam, 
      Polynomials, Plots

#          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

@variables t
D = Differential(t)
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);
kin_visc_T = fit(T_vec, my_pT.(p_l, T_vec)./rho_pT.(p_l, T_vec), 5);
lambda_T = fit(T_vec, tc_pT.(p_l, T_vec), 3);
Pr_T = fit(T_vec, 1e3*Cp_pT.(p_l, T_vec).*my_pT.(p_l, T_vec)./tc_pT.(p_l, T_vec), 5);
rho_T = fit(T_vec, rho_pT.(p_l, T_vec), 4);
rhocp_T = fit(T_vec, 1000*rho_pT.(p_l, T_vec).*Cp_pT.(p_l, T_vec), 5);

@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)
    #TODO: Vector of flow vars warning but eqs are correct!!
    sts = @variables T[1:N](t)=T0 Q[1:N](t)=Q0 [connect=Flow]
    ODESystem(Equation[], t, [T; Q], []; name=name)
end

#Taylor-aris dispersion model
function Dxx_coeff(u, d, T)
    Re = abs(u)*d/kin_visc_T(T) + 0.1
    IfElse.ifelse(Re < 1000, (d^2/4)*u^2/48/0.14e-6, d*u*(1.17e9*Re^(-2.5) + 0.41))
end

#Nusselt number model
function Nusselt(Re, Pr, f)
    IfElse.ifelse(Re <= 2300, 3.66, 
    IfElse.ifelse(Re <= 3100, 3.5239*(Re/1000)^4-45.158*(Re/1000)^3+212.13*(Re/1000)^2-427.45*(Re/1000)+316.08,
                  f/8*((Re-1000)*Pr)/(1+12.7*(f/8)^(1/2)*(Pr^(2/3)-1))))
end

#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 coeficients
    A = pi*dn^2/4
    
    p = @parameters C_shift = 0.0 Rw = 0.0 # stuff for latter
    @variables begin
        T[1:N](t) = fill(T0, N)
        Twall[1:N](t) = fill(T0, N)
        S[1:N](t) = fill(T0, N)
        C[1:N](t) = 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 = [T..., Twall..., S..., C..., u, Re, Dxx, Pr, alpha, f]
    
    
    eqs = [
        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

function Cn_circular_wall_inner(d, D, cp, ρ)
    C = pi/4*(D^2-d^2)*cp*ρ
    return C/2
end

function Cn_circular_wall_outer(d, D, cp, ρ)
    C = pi/4*(D^2-d^2)*cp*ρ
    return C/2
end

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[1:N, 1:Nn](t) = fill(T0, N, Nn)
        Qe[1:N, 1:Ne](t) = fill(T0, N, Ne)
    end
    sts = [vec(Tn); vec(Qe)]
    eqs = [
        [inner_heatport.T[i] ~ Tn[i,1] for i in 1:N]
        [outer_heatport.T[i] ~ Tn[i,Nn] for i in 1:N]
        [Qe[i,j] ~ Ke[j]*(-Tn[i,j+1]+Tn[i,j]) for i in 1:N for j in 1:Ne]...
        [D(Tn[i,1])*(Cn[1]+C_shift) ~ inner_heatport.Q[i]-Qe[i,1] for i in 1:N]
        [D(Tn[i,j])*Cn[j] ~ Qe[i,j-1]-Qe[i,j] for i in 1:N for j in 2:Nn-1]...
        [D(Tn[i,Nn])*Cn[Nn] ~ Qe[i,Ne]+outer_heatport.Q[i] for i in 1:N]
    ]
    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 test_speed(N; ODAE=true, solver=Tsit5())
    tspan = (0.0, 19*3600)
    @named testbench = TestBenchPreinsulated(L=470, N=N, dn=0.3127, t_layer=[0.0056, 0.058])
    sys = structural_simplify(testbench)
    if ODAE
        prob = ODAEProblem(sys, [], tspan, jac=true, sparse=true)
    else
        prob = ODEProblem(sys, [], tspan, jac=true, sparse=true)
    end
    prob.u0[:] .= 12.0
    #TODO: unrecognized keywords jac and sparse in solve??
    solve(prob, solver, reltol=1e-6, abstol=1e-6, saveat=100);
    return @elapsed solve(prob, solver, reltol=1e-6, abstol=1e-6, saveat=100);
end

# scaling tests
N_x = [5, 10, 20, 40, 80, 160, 320, 640, 1280]
N_states = 4 .* N_x
OMC_IDA_470 = [0.0125157, 0.0106602, 0.0172244, 0.0255715, 
               0.0567469, 0.126823, 0.247737, 0.622899, 1.32534]

JULIA_ODAE_Tsit5 = zeros(length(N_x))
for i in 1:length(N_x)
    JULIA_ODAE_Tsit5[i] = test_speed(N_x[i], solver=Tsit5())
end

JULIA_ODAE_QNDF = zeros(length(N_x))
for i in 1:length(N_x)
    JULIA_ODAE_QNDF[i] = test_speed(N_x[i], solver=QNDF(autodiff=false))
end

fig = plot(N_x, OMC_IDA, label="OMC IDA", 
           legend=:topleft, yscale=:log10, xscale=:log10,
           xlim=(1,10000), ylim=(1e-3, 1e2), marker=:circle)
plot!(N_x, JULIA_ODAE_Tsit5, label="Julia ODAE Tsit5", marker=:square)
plot!(N_x, JULIA_ODAE_QNDF, label="Julia ODAE QNDF", marker=:star)
savefig(fig,"Scaling experiment.png")

Modelica Version:

package DhnControl
  package Models
    connector FluidPort
      Real p;
      flow Real m;
      stream Real T;
      annotation(
        Icon(graphics = {Ellipse(fillColor = {1, 132, 255}, fillPattern = FillPattern.Solid, extent = {{-100, 100}, {100, -100}})}));
    end FluidPort;

    connector HeatPort
      Real T;
      flow Real Q;
      annotation(
        Icon(graphics = {Rectangle(fillColor = {255, 0, 0}, fillPattern = FillPattern.Solid, extent = {{-100, 100}, {100, -100}})}));
    end HeatPort;

    model FluidRegion
      parameter Real T0 = 0.0;
      parameter Real lumped_T = 20;
      parameter Real L = 100;
      parameter Integer N = 100;
      parameter Real dn = 0.05;
      parameter Real eps = 1e-4;
      parameter Real C_shift = 0.0;
      parameter Real Rw = 0.0;
      parameter Real dx = L / N;
      parameter Real A = pi * dn ^ 2 / 4;
      //variables
      Real T[N](each start = T0);
      Real Twall[N](each start = T0);
      Real S[N];
      Real C[N];
      Real u;
      Real Re;
      Real Dxx;
      Real Pr;
      Real alpha;
      Real f;
      Real Nu;
      constant Real c[3] = {-1 / 8, -3 / 8, -3 / 8};
      constant Real pi = Modelica.Constants.pi;
      DhnControl.Models.FluidPort inlet annotation(
        Placement(visible = true, transformation(origin = {-90, 0}, extent = {{-10, -10}, {10, 10}}, rotation = 0), iconTransformation(origin = {-90, 0}, extent = {{-10, -10}, {10, 10}}, rotation = 0)));
      DhnControl.Models.FluidPort outlet annotation(
        Placement(visible = true, transformation(origin = {90, 0}, extent = {{-10, -10}, {10, 10}}, rotation = 0), iconTransformation(origin = {90, 0}, extent = {{-10, -10}, {10, 10}}, rotation = 0)));
      DhnControl.Models.HeatPort heatport[N] annotation(
        Placement(visible = true, transformation(origin = {0, 30}, extent = {{-10, -10}, {10, 10}}, rotation = 0), iconTransformation(origin = {-2, 30}, extent = {{-10, -10}, {10, 10}}, rotation = 0)));
    equation
      Re = 0.1 + dn * abs(u) / (1.7631408357540536e-6 - 5.37008081108929e-8 * lumped_T + 9.71612371740893e-10 * lumped_T ^ 2 - 1.0026133706457947e-11 * lumped_T ^ 3 + 5.3314727276417887e-14 * lumped_T ^ 4 - 1.1234839851121354e-16 * lumped_T ^ 5);
      Pr = 13.167521706890668 - 0.4441712774475796 * lumped_T + 0.008404050163010567 * lumped_T ^ 2 - 8.863470332920757e-5 * lumped_T ^ 3 + 4.768646349109186e-7 * lumped_T ^ 4 - 1.0118106644874892e-9 * lumped_T ^ 5;
      if Re <= 2300 then
        Nu = 3.66;
      elseif Re > 2300 and Re <= 3100 then
        Nu = 3.5239 * (Re / 1000) ^ 4 - 45.158 * (Re / 1000) ^ 3 + 212.13 * (Re / 1000) ^ 2 - 427.45 * (Re / 1000) + 316.08;
      else
        Nu = f / 8 * ((Re - 1000) * Pr) / (1 + 12.7 * (f / 8) ^ (1 / 2) * (Pr ^ (2 / 3) - 1));
      end if;
      if Re < 1000 then
        Dxx = dn ^ 2 / 4 * u ^ 2 / 48 / 0.14e-6;
      else
        Dxx = dn * u * (1.17e9 * Re ^ (-2.5) + 0.41);
      end if;
      f = 8 * ((8 / Re) ^ 12 + 1 / ((-2.457 * log((7 / Re) ^ 0.9 + 0.27 * (eps / dn))) ^ 16 + (37530 / Re) ^ 16) ^ 1.5) ^ (1 / 12);
//churchill
      alpha = Nu * (0.5629006705766969 + 0.002027906870916878 * lumped_T - 1.0062563416770859e-5 * lumped_T ^ 2 + 1.2897392253800086e-8 * lumped_T ^ 3) / dn;
      inlet.m = -outlet.m;
      inlet.p = outlet.p;
      inlet.T = inStream(inlet.T);
      outlet.T = T[N];
      u = inlet.m / (1000.2325951116342 + 0.02341353916404164 * inlet.T - 0.006409744347791998 * inlet.T ^ 2 + 2.6080324835019334e-5 * inlet.T ^ 3 - 6.044425395404277e-8 * inlet.T ^ 4) / A;
      C[:] = {dx * A * (4.2146276665115245e6 - 2231.4603937980637 * T[i] + 24.383955576707972 * T[i] ^ 2 - 0.3896654168053555 * T[i] ^ 3 + 0.002537945351479517 * T[i] ^ 4 - 5.879294024916786e-6 * T[i] ^ 5) for i in 1:N};
      S[:] = {heatport[i].Q for i in 1:N};
      Twall[:] = {heatport[i].T for i in 1:N};
      S[:] = {1 / (1 / (alpha * dn * pi * dx) + abs(Rw / 1000)) * (Twall[i] - T[i]) for i in 1:N};
      der(T[1]) = u / dx * (inlet.T - T[1]) + Dxx * (T[2] - T[1]) / dx ^ 2 + S[1] / (C[1] - C_shift);
      der(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);
      der(T[3:N - 1]) = {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};
      der(T[N]) = u / dx * (T[N - 1] - T[N]) + Dxx * (T[N - 1] - T[N]) / dx ^ 2 + S[N] / (C[N] - C_shift);
      annotation(
        Icon(graphics = {Rectangle(fillColor = {50, 151, 213}, fillPattern = FillPattern.Solid, extent = {{-80, 20}, {80, -20}})}));
    end FluidRegion;

    model Source
      parameter Real t_start = 12 * 3600;
      parameter Real p_feed = 100000;
      parameter Real m_flow = 1.0;
      parameter Real T_out = 1.0;
      FluidPort outlet annotation(
        Placement(visible = true, transformation(origin = {90, 0}, extent = {{-10, -10}, {10, 10}}, rotation = 0), iconTransformation(origin = {90, 0}, extent = {{-10, -10}, {10, 10}}, rotation = 0)));
    equation
      if time > t_start then
        outlet.T = 56.0 + 12.0;
      else
        outlet.T = 12.0;
      end if;
      outlet.m = -m_flow;
      outlet.p = p_feed;
      annotation(
        Icon(graphics = {Rectangle(origin = {-10, 0}, fillColor = {204, 204, 204}, fillPattern = FillPattern.Solid, extent = {{-90, 100}, {90, -100}}), Text(origin = {-12, 52}, extent = {{-70, 24}, {70, -24}}, textString = "%name")}),
        experiment(StartTime = 0, StopTime = 1, Tolerance = 1e-6, Interval = 0.002));
    end Source;

    model Sink
      FluidPort inlet annotation(
        Placement(visible = true, transformation(origin = {-90, 0}, extent = {{-10, -10}, {10, 10}}, rotation = 0), iconTransformation(origin = {-90, 0}, extent = {{-10, -10}, {10, 10}}, rotation = 0)));
    equation
      inlet.T = inStream(inlet.T);
      annotation(
        Icon(graphics = {Rectangle(origin = {10, 0}, fillColor = {212, 212, 212}, fillPattern = FillPattern.Solid, extent = {{-90, 100}, {90, -100}}), Text(origin = {3, 53}, extent = {{-63, 25}, {63, -25}}, textString = "%name")}));
    end Sink;

    model CylindricalWallFem
      parameter Real T0 = 0.0;
      parameter Real d = 0.3;
      parameter Real L = 100;
      parameter Integer N = 100;
      parameter Real dx = L / N;
      parameter Real t_layer[2] = {0.00391, 0.013};
      parameter Integer Ne = size(t_layer,1);
      parameter Real cp[Ne] = {500, 1200};
      parameter Real rho[Ne] = {7850, 40};
      parameter Real lambda[Ne] = {50, 0.04};
      parameter Integer Nn = Ne + 1;
      parameter Real dn[Nn] = cat(1, {d}, d .+ 2.0.*{sum(t_layer[1:i]) for i in 1:Ne});
      parameter Real Cn[Nn] = cat(1, {pi/8*(dn[i+1]^2-dn[i]^2)*cp[i]*rho[i] for i in 1:Ne}, {0})+ cat(1, {0}, {pi/8*(dn[i+1]^2-dn[i]^2)*cp[i]*rho[i] for i in 1:Ne});
      parameter Real Ke[Ne] = {2*pi*lambda[i]/log(dn[i+1]/dn[i]) for i in 1:Ne};
      Real Tn[N,Nn](each start = T0);
      Real Qe[N,Ne];
      DhnControl.Models.HeatPort inner_heatport[N] annotation(
          Placement(visible = true, transformation(origin = {0, -90}, extent = {{-10, -10}, {10, 10}}, rotation = 0), iconTransformation(origin = {0, -90}, extent = {{-10, -10}, {10, 10}}, rotation = 0)));
      DhnControl.Models.HeatPort outer_heatport[N] annotation(
          Placement(visible = true, transformation(origin = {0, 92}, extent = {{-10, -10}, {10, 10}}, rotation = 0), iconTransformation(origin = {0, 90}, extent = {{-10, -10}, {10, 10}}, rotation = 0)));
      constant Real pi = Modelica.Constants.pi;
    equation
      for i in 1:N loop
        inner_heatport[i].T = Tn[i,1];
        outer_heatport[i].T = Tn[i,Nn];
        Qe[i,:] = {Ke[j]*(-Tn[i,j+1]+Tn[i,j]) for j in 1:Ne};
      end for;
      
      der(Tn[:,1])*(Cn[1]) = {inner_heatport[i].Q - Qe[i,1] for i in 1:N};
      for i in 1:N loop
        der(Tn[i,2:Ne]).*Cn[2:Ne] = {Qe[i,j-1]-Qe[i,j] for j in 2:Ne};
      end for;
      der(Tn[:,Nn])*Cn[Nn] = {Qe[i, Ne] + outer_heatport[i].Q for i in 1:N};
    
      annotation(
        Icon(graphics = {Rectangle(fillColor = {214, 197, 109}, fillPattern = FillPattern.Solid, extent = {{-100, 80}, {100, -80}}), Text(origin = {-4, 9}, extent = {{-72, 27}, {72, -27}}, textString = "%name")}),
        experiment(StartTime = 0, StopTime = 1, Tolerance = 1e-06, Interval = 0.002));
    end CylindricalWallFem;

    model PreinsulatedPipe
      parameter Real T0 = 0.0;
      parameter Real lumped_T = 20.0;
      parameter Real eps = 1e-4;
      parameter Real L =470;
      parameter Integer N = 100;
      parameter Real dn = 0.3;
      parameter Real t_layer[:] = {0.0056, 0.058};
      parameter Real cp[:] = {500, 1200};
      parameter Real rho[:] = {7850, 40};
      parameter Real lambda[:] = {50, 0.04};
      parameter Real alpha_out = 4.0;
      parameter Real Tenv = 18.0;
      DhnControl.Models.FluidPort inlet annotation(
        Placement(visible = true, transformation(origin = {-90, 0}, extent = {{-10, -10}, {10, 10}}, rotation = 0), iconTransformation(origin = {-90, 0}, extent = {{-10, -10}, {10, 10}}, rotation = 0)));
      DhnControl.Models.FluidPort outlet annotation(
        Placement(visible = true, transformation(origin = {90, 0}, extent = {{-10, -10}, {10, 10}}, rotation = 0), iconTransformation(origin = {90, 0}, extent = {{-10, -10}, {10, 10}}, rotation = 0)));
      DhnControl.Models.FluidRegion fluidRegion(L=L, N=N, dn=dn, T0=T0, lumped_T=lumped_T,eps=eps) annotation(
        Placement(visible = true, transformation(origin = {0, 0}, extent = {{-10, -10}, {10, 10}}, rotation = 0)));
      DhnControl.Models.CylindricalWallFem cylindricalWallFem(t_layer=t_layer, L=L, N=N, d=dn, rho=rho, lambda=lambda, cp=cp, T0=T0) annotation(
        Placement(visible = true, transformation(origin = {0, 24}, extent = {{-10, -10}, {10, 10}}, rotation = 0)));
  DhnControl.Models.CylindricalSurfaceConvection cylindricalSurfaceConvection(L=L, N=N, d=dn+2*sum(t_layer), alpha=alpha_out,Tenv=Tenv) annotation(
        Placement(visible = true, transformation(origin = {0, 50}, extent = {{-10, -10}, {10, 10}}, rotation = 0)));
    equation
      connect(fluidRegion.outlet, outlet) annotation(
        Line(points = {{9, 0}, {90, 0}}));
      connect(fluidRegion.inlet, inlet) annotation(
        Line(points = {{-9, 0}, {-90, 0}}));
      connect(fluidRegion.heatport, cylindricalWallFem.inner_heatport) annotation(
        Line(points = {{0, 4}, {0, 16}}, thickness = 0.5));
  connect(cylindricalWallFem.outer_heatport, cylindricalSurfaceConvection.heatport) annotation(
        Line(points = {{0, 34}, {0, 42}}, thickness = 0.5));
    annotation(
        Icon(graphics = {Rectangle( lineColor = {254, 247, 33}, fillColor = {108, 115, 77}, fillPattern = FillPattern.HorizontalCylinder, extent = {{-80, 34}, {80, -34}}),Rectangle(fillColor = {108, 108, 108}, fillPattern = FillPattern.HorizontalCylinder, extent = {{-80, 24}, {80, -24}}), Rectangle(lineColor = {206, 206, 206},fillColor = {26, 139, 238}, fillPattern = FillPattern.HorizontalCylinder, extent = {{-80, 20}, {80, -20}})}));
        end PreinsulatedPipe;

    model CylindricalSurfaceConvection
    parameter Real d = 0.3;
    parameter Real L = 470;
    parameter Integer N = 100;
    parameter Real Tenv = 18;
    parameter Real dx= L/N;
    parameter Real alpha = 4.0;
    constant Real pi = Modelica.Constants.pi;
    DhnControl.Models.HeatPort heatport[N] annotation(
        Placement(visible = true, transformation(origin = {0, -90}, extent = {{-10, -10}, {10, 10}}, rotation = 0), iconTransformation(origin = {0, -88}, extent = {{-10, -10}, {10, 10}}, rotation = 0)));
    equation
for i in 1:N loop
    heatport[i].Q = alpha*pi*d*dx*(heatport[i].T-Tenv);
    end for;
    annotation(
        Icon(graphics = {Line(origin = {-1.45, -60.62}, points = {{-100, 0}, {100, 0}}, thickness = 1), Line(origin = {-70.06, 0.679544}, points = {{9.76713, -58.9684}, {-10.2329, -20.9684}, {9.76713, 11.0316}, {-8.23287, 59.0316}}, color = {255, 0, 0}, thickness = 0.75, arrow = {Arrow.None, Arrow.Filled}, smooth = Smooth.Bezier), Line(origin = {-30.7239, 0.513569}, points = {{9.76713, -58.9684}, {-10.2329, -20.9684}, {9.76713, 11.0316}, {-8.23287, 59.0316}}, color = {255, 0, 0}, thickness = 0.75, arrow = {Arrow.None, Arrow.Filled}, smooth = Smooth.Bezier), Line(origin = {10.9359, 0.63805}, points = {{9.76713, -58.9684}, {-10.2329, -20.9684}, {9.76713, 11.0316}, {-8.23287, 59.0316}}, color = {255, 0, 0}, thickness = 0.75, arrow = {Arrow.None, Arrow.Filled}, smooth = Smooth.Bezier), Line(origin = {52.5957, 0.762531}, points = {{9.76713, -58.9684}, {-10.2329, -20.9684}, {9.76713, 11.0316}, {-8.23287, 59.0316}}, color = {255, 0, 0}, thickness = 0.75, arrow = {Arrow.None, Arrow.Filled}, smooth = Smooth.Bezier)}));
        end CylindricalSurfaceConvection;
  end Models;

  package Test
    model test_adiabatic_470
      DhnControl.Models.Source source(m_flow = 2.75) annotation(
        Placement(visible = true, transformation(origin = {-50, 10}, extent = {{-10, -10}, {10, 10}}, rotation = 0)));
      DhnControl.Models.FluidRegion fluidRegion(L = 470, N = 470, T0 = 12, dn = 0.3127) annotation(
        Placement(visible = true, transformation(origin = {-12, 10}, extent = {{-10, -10}, {10, 10}}, rotation = 0)));
      Models.Sink sink annotation(
        Placement(visible = true, transformation(origin = {28, 10}, extent = {{-10, -10}, {10, 10}}, rotation = 0)));
    equation
      connect(source.outlet, fluidRegion.inlet) annotation(
        Line(points = {{-40, 10}, {-21, 10}}));
      connect(fluidRegion.outlet, sink.inlet) annotation(
        Line(points = {{-3, 10}, {20, 10}}));
      annotation(
        experiment(StartTime = 0, StopTime = 68400, Tolerance = 1e-06, Interval = 100));
    end test_adiabatic_470;
    
    model test_preinsulated_470
      DhnControl.Models.Source source(m_flow = 2.75) annotation(
        Placement(visible = true, transformation(origin = {-50, 10}, extent = {{-10, -10}, {10, 10}}, rotation = 0)));
      Models.Sink sink annotation(
        Placement(visible = true, transformation(origin = {28, 10}, extent = {{-10, -10}, {10, 10}}, rotation = 0)));
    Models.PreinsulatedPipe preinsulatedPipe(L = 470, N = 1280, T0 = 12, dn = 0.3127)  annotation(
        Placement(visible = true, transformation(origin = {-12, 10}, extent = {{-10, -10}, {10, 10}}, rotation = 0)));
    equation
      connect(source.outlet, preinsulatedPipe.inlet) annotation(
        Line(points = {{-40, 10}, {-20, 10}}));
    connect(preinsulatedPipe.outlet, sink.inlet) annotation(
        Line(points = {{-2, 10}, {20, 10}}));
      annotation(
        experiment(StartTime = 0, StopTime = 68400, Tolerance = 1e-06, Interval = 100),
        __OpenModelica_simulationFlags(lv = "LOG_STATS", s = "ida"));
    end test_preinsulated_470;
    
  end Test;
end DhnControl;

]st

[6e4b80f9] BenchmarkTools v1.3.1
[336ed68f] CSV v0.10.4
[052768ef] CUDA v3.11.0
[a93c6f00] DataFrames v1.3.4
[82cc6244] DataInterpolations v3.9.2
[39dd38d3] Dierckx v0.5.2
[aae7a2af] DiffEqFlux v1.51.2
[41bf760c] DiffEqSensitivity v6.79.0
[0c46a032] DifferentialEquations v7.2.0
[f6369f11] ForwardDiff v0.10.30
[615f187c] IfElse v0.1.1
[a98d9a8b] Interpolations v0.14.0
[ef3ab10e] KLU v0.3.0
[7f56f5a3] LSODA v0.7.0
[b964fa9f] LaTeXStrings v1.3.0
[7ed4a6bd] LinearSolve v1.20.0
[eff96d63] Measurements v2.7.2
[961ee093] ModelingToolkit v8.15.1
[7f7a1694] Optimization v3.7.1
[1dea7af3] OrdinaryDiffEq v6.18.1
[f0f68f2c] PlotlyJS v0.18.8
[91a5bcdd] Plots v1.31.2
[f27b6e38] Polynomials v3.1.4
[1fd47b50] QuadGK v2.4.2
[731186ca] RecursiveArrayTools v2.31.1
[276daf66] SpecialFunctions v2.1.7
[43dc94dd] SteamTables v1.3.0
[c3572dad] Sundials v4.9.4
[0c5d862f] Symbolics v4.9.0
[a759f4b9] TimerOutputs v0.5.20
[95ff35a0] XSteam v0.3.0

OMC
OpenModelica v1.19.2

[baggepinnen]

Hi @LiborKudela !

Are you solving to the same tolerances with both MTK and modelica?

[ChrisRackauckas]

1. is sparse truly better here?
2. autodiff=false shouldn't matter if jac=true?
3. Did you try FBDF? That's close to CVODE/IDA
4. Which linear solver? Did you try KLU?
5. How does CVODE do? Or IDA from Julia?

[LiborKudela]

Sorry for the delay.. A lot of compilation time with some of these variants..
I had to terminate many I think they would not come to an end

Are you solving to the same tolerances with both MTK and modelica?

Yes, both are set to 1e-6. As well as the saveat=100 equivalent in OMC.

2. autodiff=false shouldn't matter if jac=true?

when autodiff=true first call to solve takes much longer time (for higher N) otherwise no effect -> that might be a bug

1. is sparse truly better here?

Well, the derivations of the states need only few other states around it so it should be better, but I will check.

3. Did you try FBDF? That's close to CVODE/IDA

sparse= false -> no effect, which I actually did not expect

5. How does CVODE do? Or IDA from Julia?

sparse=true, I could not make any DAE form solve faster inluding IDA. I have got a lot of crashes with dae form.
I might not be using the IDA correctly.
the CVODE_BDF works and also scales bad.

4. Which linear solver? Did you try KLU?

I did try KLU but it gave me an Error so I moved on. Now when I se the Error again I have questions.

solve(prob, QNDF(autodiff=false, linsolve = KLUFactorization()), reltol=1e-6, abstol=1e-6, saveat=100)

ERROR: MethodError: no method matching KLU.KLUFactorization(::Matrix{Float64})

Which is weird because i have set sparse=true in ODAEProblem so why it passes in non sparse matrix?

Btw I also get this warnings comming from solve

Unrecognized keyword arguments: [:jac, :sparse]

Is it posible that the sparsity matrix and jacobian are not comming through from prob or something?
How can we check that?

[ChrisRackauckas]

Share ]st.

[LiborKudela]

[6e4b80f9] BenchmarkTools v1.3.1
[336ed68f] CSV v0.10.4
[052768ef] CUDA v3.11.0
[a93c6f00] DataFrames v1.3.4
[82cc6244] DataInterpolations v3.9.2
[39dd38d3] Dierckx v0.5.2
[aae7a2af] DiffEqFlux v1.51.2
[41bf760c] DiffEqSensitivity v6.79.0
[0c46a032] DifferentialEquations v7.2.0
[f6369f11] ForwardDiff v0.10.30
[615f187c] IfElse v0.1.1
[a98d9a8b] Interpolations v0.14.0
[ef3ab10e] KLU v0.3.0
[7f56f5a3] LSODA v0.7.0
[b964fa9f] LaTeXStrings v1.3.0
[7ed4a6bd] LinearSolve v1.20.0
[eff96d63] Measurements v2.7.2
[961ee093] ModelingToolkit v8.15.1
[7f7a1694] Optimization v3.7.1
[1dea7af3] OrdinaryDiffEq v6.18.1
[f0f68f2c] PlotlyJS v0.18.8
[91a5bcdd] Plots v1.31.2
[f27b6e38] Polynomials v3.1.4
[1fd47b50] QuadGK v2.4.2
[731186ca] RecursiveArrayTools v2.31.1
[276daf66] SpecialFunctions v2.1.7
[43dc94dd] SteamTables v1.3.0
[c3572dad] Sundials v4.9.4
[0c5d862f] Symbolics v4.9.0
[a759f4b9] TimerOutputs v0.5.20
[95ff35a0] XSteam v0.3.0

[LiborKudela]

Great news!! (sort of)
The problem with sparsity and KLU not working is only happening with ODAEProblem.
I have tried pure ODEProblem so we have a non-singular mass matrix:

tspan = (0.0, 19*3600)
@named testbench = TestBenchPreinsulated(L=470, N=160, dn=0.3127, t_layer=[0.0056, 0.058])
sys = structural_simplify(testbench)
prob = ODEProblem(sys, [], tspan, jac=true, sparse=true)
prob.u0 .= 12
@time sol = solve(prob, QBDF(linsolve=KLUFactorization()), reltol=1e-6, abstol=1e-6, saveat=100);
@time sol = solve(prob, CVODE_BDF(linear_solver=:KLU), reltol=1e-6, abstol=1e-6, saveat=100);

Result: 😄 I knew that Julia should be faster.

But when I go from N=160 to N=320 it gets stuck at compilation and never actually runs.
My computer at home has only 16GB of RAM and it gets pushed to limit.
I will try it at school tomorrow on bigger machine, but this should be probably addressed?
OpenModelica compiles this sizes easily even on my home PC.

[ChrisRackauckas]

Modelica defaults to KLU. I just did a bunch of tests and it makes sense for us too, so I changed the default linear solver. UMFPACK is rarely better.

As for the compilation time, this is a known issue with MTK that is being worked on.

[YingboMa]

Hey, could you run the above benchmark again using the latest MTK? It should be able to handle larger systems now.

[LiborKudela]

Hi, I will run it later today. Thank you!

[LiborKudela]

@YingboMa, I have run the benchmark.

First I have updated everything.

Julia version:
1.8.2

]st :
[615f187c] IfElse v0.1.1
[7ed4a6bd] LinearSolve v1.27.0
[961ee093] ModelingToolkit v8.29.1
[1dea7af3] OrdinaryDiffEq v6.29.3
[91a5bcdd] Plots v1.35.4
[f27b6e38] Polynomials v3.2.0
[0c5d862f] Symbolics v4.13.0
[95ff35a0] XSteam v0.3.

I have also changed the symbolic arrays from A[1:N](t) to (A(t))[1:N] since the former is deprecated in Symbolics.
I have run only the ODEProblem version of the benchmark.
The test machine is a Linux machine with 16GB of RAM, the OS is Ubuntu 20.04 (freshly updated).

I am unable to run N=1280 with Tsit5() solver - it crashes the process (in VScode), because it eats all the RAM
I am unable to run N=320 with QBDF(linsolve=KLUFactorization()) solver with the same result as the above
I was able to run the OpenModelica version of the model with N=2560 on this machine relative fast (simulation time 10.5s)

Should I change something more in the script to take advantage of the new version of MTK?

The script I used for this test:

using ModelingToolkit, OrdinaryDiffEq, Symbolics, IfElse, XSteam, 
      Polynomials, Plots, LinearSolve

#          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

@variables t
D = Differential(t)
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);
kin_visc_T = fit(T_vec, my_pT.(p_l, T_vec)./rho_pT.(p_l, T_vec), 5);
lambda_T = fit(T_vec, tc_pT.(p_l, T_vec), 3);
Pr_T = fit(T_vec, 1e3*Cp_pT.(p_l, T_vec).*my_pT.(p_l, T_vec)./tc_pT.(p_l, T_vec), 5);
rho_T = fit(T_vec, rho_pT.(p_l, T_vec), 4);
rhocp_T = fit(T_vec, 1000*rho_pT.(p_l, T_vec).*Cp_pT.(p_l, T_vec), 5);

@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)
    #TODO: Vector of flow vars warning but eqs are correct!!
    sts = @variables (T(t))[1:N]=T0 (Q(t))[1:N]=Q0 [connect=Flow]
    ODESystem(Equation[], t, [T; Q], []; name=name)
end

#Taylor-aris dispersion model
function Dxx_coeff(u, d, T)
    Re = abs(u)*d/kin_visc_T(T) + 0.1
    IfElse.ifelse(Re < 1000, (d^2/4)*u^2/48/0.14e-6, d*u*(1.17e9*Re^(-2.5) + 0.41))
end

#Nusselt number model
function Nusselt(Re, Pr, f)
    IfElse.ifelse(Re <= 2300, 3.66, 
    IfElse.ifelse(Re <= 3100, 3.5239*(Re/1000)^4-45.158*(Re/1000)^3+212.13*(Re/1000)^2-427.45*(Re/1000)+316.08,
                  f/8*((Re-1000)*Pr)/(1+12.7*(f/8)^(1/2)*(Pr^(2/3)-1))))
end

#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 coeficients
    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 = [T..., Twall..., S..., C..., u, Re, Dxx, Pr, alpha, f]
    
    
    eqs = [
        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

function Cn_circular_wall_inner(d, D, cp, ρ)
    C = pi/4*(D^2-d^2)*cp*ρ
    return C/2
end

function Cn_circular_wall_outer(d, D, cp, ρ)
    C = pi/4*(D^2-d^2)*cp*ρ
    return C/2
end

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)]
    eqs = [
        [inner_heatport.T[i] ~ Tn[i,1] for i in 1:N]
        [outer_heatport.T[i] ~ Tn[i,Nn] for i in 1:N]
        [Qe[i,j] ~ Ke[j]*(-Tn[i,j+1]+Tn[i,j]) for i in 1:N for j in 1:Ne]...
        [D(Tn[i,1])*(Cn[1]+C_shift) ~ inner_heatport.Q[i]-Qe[i,1] for i in 1:N]
        [D(Tn[i,j])*Cn[j] ~ Qe[i,j-1]-Qe[i,j] for i in 1:N for j in 2:Nn-1]...
        [D(Tn[i,Nn])*Cn[Nn] ~ Qe[i,Ne]+outer_heatport.Q[i] for i in 1:N]
    ]
    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 test_speed(N; solver=Tsit5())
    tspan = (0.0, 19*3600)
    @named testbench = TestBenchPreinsulated(L=470, N=N, dn=0.3127, t_layer=[0.0056, 0.058])
    sys = structural_simplify(testbench)
    prob = ODEProblem(sys, [], tspan, jac=true, sparse=true)
    prob.u0[:] .= 12.0
    #TODO: unrecognized keywords jac and sparse in solve??
    solve(prob, solver, reltol=1e-6, abstol=1e-6, saveat=100);
    return @elapsed solve(prob, solver, reltol=1e-6, abstol=1e-6, saveat=100);
end


test_speed(1280)

# scaling tests
N_x = [5, 10, 20, 40, 80, 160, 320, 640, 1280]
N_states = 4 .* N_x
OMC_IDA_470 = [0.0125157, 0.0106602, 0.0172244, 0.0255715, 
               0.0567469, 0.126823, 0.247737, 0.622899, 1.32534]

JULIA_ODAE_Tsit5 = zeros(length(N_x))
for i in 1:length(N_x)
    JULIA_ODAE_Tsit5[i] = test_speed(N_x[i], solver=Tsit5())
end
println(JULIA_ODAE_Tsit5)
#[0.071787957, 0.077838405, 0.095070696, 0.136035788, 0.227302308, 0.422744682, 0.90749644, 1.906403591, NaN]

JULIA_ODAE_QNDF = zeros(length(N_x))
for i in 1:length(N_x)
    JULIA_ODAE_QNDF[i] = test_speed(N_x[i], solver=QBDF(linsolve=KLUFactorization()))
end
println(JULIA_ODAE_QNDF)

fig = plot(N_x, OMC_IDA, label="OMC IDA", 
           legend=:topleft, yscale=:log10, xscale=:log10,
           xlim=(1,10000), ylim=(1e-3, 1e2), marker=:circle)
plot!(N_x, JULIA_ODAE_Tsit5, label="Julia ODAE Tsit5", marker=:square)
plot!(N_x, JULIA_ODAE_QNDF, label="Julia ODAE QNDF", marker=:star)
savefig(fig,"Scaling experiment.png")

[YingboMa]

Thanks for taking the time to run it again. For large and complex systems, symbolically computing the Jacobian won't be a good idea because the size of the expression will blowup, so I suggest you to use

prob = ODEProblem(sys, [], tspan, sparse=true)

[YingboMa]

Also, we are now defaulting to KLU for sparse matrices with moderate size, so you can just use QBDF(autodiff=false), QNDF(autodiff=false), or FBDF(autodiff=false).

[YingboMa]

BTW, the warning is fixed in #1895

[LiborKudela]

Results with jac=false.
Julia run from terminal/bash (no VScode involved).

N=1280 and above still causes crashes for both the Tsit5() and QNDF(autodiff=false).
The process runs out of RAM and gets killed.

[YingboMa]

I have modified the script a bit to make it friendlier to the compiler. Could you try this as well?

using ModelingToolkit, OrdinaryDiffEq, Symbolics, IfElse, XSteam,
      Polynomials, Plots, LinearSolve

#          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

@variables t
D = Differential(t)
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);
@noinline @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...,)))
@noinline @generated lambda_T(t)   = :(Base.evalpoly(t, $(fit(T_vec, tc_pT.(p_l, T_vec), 3).coeffs...,)))
@noinline @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...,)))
@noinline @generated rho_T(t)      = :(Base.evalpoly(t, $(fit(T_vec, rho_pT.(p_l, T_vec), 4).coeffs...,)))
@noinline @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
@noinline function Dxx_coeff(u, d, T)
    Re = abs(u)*d/kin_visc_T(T) + 0.1
    IfElse.ifelse(Re < 1000, (d^2/4)*u^2/48/0.14e-6, d*u*(1.17e9*Re^(-2.5) + 0.41))
end

@register_symbolic Nusselt(Re, Pr, f)
#Nusselt number model
@noinline function Nusselt(Re, Pr, f)
    if Re <= 2300
        3.66
    elseif Re <= 3100
        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
@noinline 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 coeficients
    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 = [T..., Twall..., S..., C..., u, Re, Dxx, Pr, alpha, f]
    
    
    eqs = [
        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)]
    eqs = [
        [inner_heatport.T[i] ~ Tn[i,1] for i in 1:N]
        [outer_heatport.T[i] ~ Tn[i,Nn] for i in 1:N]
        [Qe[i,j] ~ Ke[j]*(-Tn[i,j+1]+Tn[i,j]) for i in 1:N for j in 1:Ne]...
        [D(Tn[i,1])*(Cn[1]+C_shift) ~ inner_heatport.Q[i]-Qe[i,1] for i in 1:N]
        [D(Tn[i,j])*Cn[j] ~ Qe[i,j-1]-Qe[i,j] for i in 1:N for j in 2:Nn-1]...
        [D(Tn[i,Nn])*Cn[Nn] ~ Qe[i,Ne]+outer_heatport.Q[i] for i in 1:N]
    ]
    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