`jac` and `sparse` are not implemented for DAEProblem

[ChrisRackauckas]

using OrdinaryDiffEq, ModelingToolkit
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)
prob3 = DAEProblem(sys, du0, u0, tspan, parameters, jac=true, sparse=true)
prob3.f.jac_prototype == nothing
prob3.f.jac == nothing

We might want to make this throw a warning until it's handled.

https://discourse.julialang.org/t/modelingtoolkit-building-a-daeproblem-and-solving-it-with-ida-with-klu-jacobian-usage/77814

[ChrisRackauckas]

The DAE definition is https://diffeq.sciml.ai/stable/features/performance_overloads/#Declaring-Explicit-Jacobians-for-DAEs
so you'd need to compute the Jacobian w.r.t. the du's, and the u's, then build the function from those two.

Just generate the sparse Jacobian for both of the two, and then add them
so it's very similar to the ODEProblem case, just where you compute two Jacobians and symbolically add them before you build the function from it.

See this
https://github.com/SciML/ModelingToolkit.jl/blob/master/src/systems/diffeqs/abstractodesystem.jl#L348-L354
it creates the Jacobian function from the symbolic Jacobian
https://github.com/SciML/ModelingToolkit.jl/blob/master/src/systems/diffeqs/abstractodesystem.jl#L69-L73
https://github.com/SciML/ModelingToolkit.jl/blob/master/src/systems/diffeqs/abstractodesystem.jl#L20-L39
calculate_jacobian is how df/du is symbolically calculated
because dvs = states(sys) is the states
rhs is the symbolic f
I would do it like this.

1. Extend calculate_jacobian to include a default arguement dvs = states(sys).
2. Create generate_dae_jacobian by using calculate_jacobian with dvs = states(sys) and dvs = Differential(independent_variable(sys)).(states(sys)) to create the two parts, and some gamma variable.
3. Add Jacobian stuff to https://github.com/SciML/ModelingToolkit.jl/blob/master/src/systems/diffeqs/abstractodesystem.jl#L446-L450

[ChrisRackauckas]

Not quite resolved, but much closer. #1510 (comment)