ModelingToolkit: Building a DAEProblem and solving it with IDA with KLU & Jacobian usage

Specific Domains Modelling & Simulations
1

I have coded an MTK system that I build a DAEProblem from. I am now trying to find out how to solve that problem using Sundials IDA with KLU as linear solver. As far as I know, this requires that the Jacobian matrix is given. Here is what I tried:

using OrdinaryDiffEq, ModelingToolkit
using DASSL
using DASKR
using DifferentialEquations
using BenchmarkTools
using Plots
using Sundials

function build_system(n)

	@parameters t # or would @variables t be better!?
    @variables Tgas[1:n](t) Tsol[1:n](t) cR[1:n](t) mflux(t)
	@parameters params[1:3]
	D = Differential(t)
    

    dp, Tin, h = params 

    dx = h / n

	############################################################
	# Build equations ##########################################
	############################################################

	eqns = Array{Equation}(undef, 3 * n + 1)
	count = 1

    # specification as function or intermediate variabel does not matter. I could also make this an 
    # actual @variable, then I could retrieve it later.
    rR = 5.5e6 .* exp.(.-5680 ./(273 .+ Tsol)) .* (1 .- exp.(.-cR./50))
    q = 1e4 * (Tgas - Tsol)
    
    # gas energy balance

    eqns[count] = 0 ~  mflux * 1000 * (Tin - Tgas[1]) - q[1] * dx
    count += 1
	for i in 2:n
	        eqns[count] = 0 ~ mflux * 1000 * (Tgas[i-1] - Tgas[i]) - q[i] * dx
		count += 1
	end

	# solids energy balance
	for i in 1:n
        	eqns[count] = D(Tsol[i]) ~ q[i] / (1000 * 2000)
		count += 1
	end

	# gas momentum balance
	eqns[count] = 0 ~ - mflux[1] + 5e-3 * sqrt(dp)
	count += 1

	# reactant concentration equations
	for i in 1:n
        	eqns[count] = D(cR[i]) ~ -rR[i]
		count += 1
	end

    @named sys = ODESystem(eqns, t, vcat(Tgas, Tsol, cR, [mflux]), params)

    return sys

end


dp = 100e2
Tin = 200
h = 0.5
parameters = [dp, Tin, h]

N = 100
tspan = (0.0, 1500.0)


sys = build_system(N)

# quick and dirty definition of u0, should be Pairs instead
u0 = ones(length(sys.states))
u0[1:N] .= parameters[2]
u0[N+1 : 2*N] .= 20
u0[2*N+1 : 3*N] .= 3000
u0[3*N+1] = 0.5

du0 = 0 .* copy(u0)

prob2 = DAEProblem(sys, du0, u0, tspan, parameters)
prob3 = DAEProblem(sys, du0, u0, tspan, parameters, jac=true, sparse=true)

jac = eval(ModelingToolkit.generate_jacobian(sys)[2])
prob4 = DAEProblem(sys, du0, u0, tspan, parameters, jac=jac, sparse=true)

diff_vars = vcat(zeros(N),ones(N), ones(N), [0])
prob5 = DAEProblem(sys, du0, u0, tspan, parameters, differential_vars=diff_vars, jac=jac, sparse=true)

prob6 = DAEProblem(sys, du0, u0, tspan, parameters, differential_vars=diff_vars, jac=true, sparse=true)


# 124ms, 54285 allocs
@benchmark solve(prob2, IDA())

# 120ms, 54289 allocs
@benchmark solve(prob3, IDA())

# 123ms, 54289 allocs
@benchmark solve(prob4, IDA())

# 125ms, 54289 allocs
@benchmark solve(prob5, IDA()) 

# 124ms, 54289 allocs
@benchmark solve(prob6, IDA()) 

# Errors with: MethodError: no method matching nonzeros(::Nothing)
@benchmark solve(prob5, IDA(linear_solver = :KLU))

I have two problems with that code:

  • First thing that is suspicious is that all solves on prob2 through prob6 seem to do pretty much the same with practically exact same performance and allocations, although the problems are set up in quite different ways (with/without jac and sparsity, with/without explicit jac, etc). This finding is similar to what was discussed a year ago here.

  • Second thing is to solve the DAEProblem with IDA using :KLU as linear solver. What do I have to do to get this working? As you see at the end of my code, my trials are erroring with the same error message as reported earlier here.

Any help is greatly appreciated.

2

Everything beyond prob3 is trying too hard. This is just a current issue with the DAEProblem constructor.

https://github.com/SciML/ModelingToolkit.jl/issues/1482

3

I was able to get it working by specifying the Jacobian and sparsity in the undocumented(?) DAEFunction. My code is along the lines of

dae_f = DAEFunction(model!; jac = jac!, jac_prototype = build_sparse(n))
prob = DAEProblem(dae_f, dICs, ICs, tspan, params, differential_vars = differential_vars)
sol = solve(prob, IDA(linear_solver = :KLU))

In my implementation jac! numerically computes the Jacobian with FiniteDiff.jl, and build_sparse(n) specifies the sparsity of the Jacobian. You could probably do this automatically, but it was very slow for me when I tried it.

You will probably need to rewrite most of your code for it to be compatible with this setup, but I don’t know of any other way.

4

Thank you, I will take a deeper look on your proposed workaround.