Explicit Jacobian on DAEs?

Specific Domains Numerics
2

Hey,
So after getting rid of the dispatch and changing back to the default linear solver it works:

using Sundials

p=[2,1.2,3]


function f(res,du,u,p,t)
  res[1] = du[1] - p[1] * u[1] + p[2] * u[1]*u[2]
  res[2] = du[2] -p[3] * u[2] - u[1]*u[2]
  res[3]=  u[2] + u[1] - u[3]
end

function g(J,du,u,p,gamma,t)
  J[1,1] = gamma - p[1] + p[2] * u[2]
  J[1,2] = p[2]* u[1]
  J[1,3] = 0
  J[2,1] = - 1 * u[2]
  J[2,2] = gamma - p[3] - u[1]
  J[2,3] = 0
  J[3,1] = 1
  J[3,2] = 1
  J[3,3] = -1
  nothing
end

fg = DAEFunction(f;jac=g)

dIC=[-0.8,-4]
IC =[1.0,1.0,2.0]
tspan=(0.0,10.0)
diff_vars = [true,true, false]

testProb2 = DAEProblem(fg,dIC,IC,tspan,p,differential_vars=diff_vars)
sol = solve(testProb2,IDA())

So that’s using your Jacobian. You can check by putting a print in there. But now you want to scale to a large problem. Instead of going straight to GMRES, if you’re willing to define your Jacobian then it’s usually more efficient to define a sparse Jacobian. You do that by defining jac_prototype and pass in a matrix with the right sparsity pattern.
Then you can directly write the Jacobian pieces in an call a sparse factorization method (:KLU). You can then directly loop through J.data to speed up the Jacobian updates instead of using Cartesian indexing.

If you are willing to do this work though, sparse factorization with a 5th order Rosenbrock method, writing the DAE in mass matrix form, can be more efficient (it’ll get more miles per factorization!), but of course than can be a pretty substantial rewrite.