Sundials KLU issue with Jacobian prototype

Specific Domains Modelling & Simulations
1

Hi everyone,

I am currently encountering issues when solving large (2814 states) stiff DAE problems (energy system time domain simulations) in Sundials IDA with the KLU linear solver (to exploit the inherent sparsity of my problem).

I constantly get the error message

exception = Sparsity Pattern in receiving SUNMatrix doesn't match sending SparseMatrix

Alternating the passed Jacobian prototype (as suggested in the thread linked below) did not resolve the problem.
However, after some primitive debugging, I have found that the exception arises from the copyto!() function in (Sundials.jl/src/types_and_consts_additions.jl, L18)
because the size of rowvals of the passed jacobian is (much) greater than the size of indexvals.
My jacobian should have 13480 non-zero entries but it has over 700k in the copyto!()function. To be exact, the passed prototype has these 13480 non-zeros until the idajac() function (Sundials.jl/src/common_interface/function_types.jl, L118) is called as seen in below.

function idajac(t::realtype,
    cj::realtype,
    u::N_Vector,
    du::N_Vector,
    res::N_Vector,
    _J::SUNMatrix,
    funjac::AbstractFunJac{<:SparseArrays.SparseMatrixCSC},
    tmp1::N_Vector,
    tmp2::N_Vector,
    tmp3::N_Vector)
    jac_prototype = funjac.jac_prototype
    N = ndims(funjac.u)
    funjac.u = unsafe_wrap(Array{Float64, N}, N_VGetArrayPointer_Serial(u), size(funjac.u))
    _u = funjac.u
    funjac.du = unsafe_wrap(Array{Float64, N}, N_VGetArrayPointer_Serial(du),
        size(funjac.du))
    _du = funjac.du

    println("DEBUG 1")
    println(size(jac_prototype.rowval))  # 13480
    println(size(jac_prototype.colptr))  # 2815

    funjac.jac(jac_prototype, _du, _u, funjac.p, cj, t) # issue comes from this line

    println("DEBUG 2")
    println(size(jac_prototype.rowval))  # 769927
    println(size(jac_prototype.colptr))  # 2815

    copyto!(_J, jac_prototype)

    return IDA_SUCCESS
end

My model is currently not solvable with a dense linear solver (somewhat ill-conditioned, the dense linear solver doesn’t converge) and therefore I opted for KLU to exploit the sparsity.

Is this exact problem a bug or issue with Sundials.jl or does it actually come from my model?

Thanks very much in advance!

(Also, sorry, posting direct links was not possible for some reason…)

2

This would come from your model. If your Jacobian function is creating new non-zeros in your sparsity pattern then this is the error that you would get. You need to make sure that your jac_prototype holds for all values of the integration.

3

Thank you very much for your help!

Then my sparsity detection is probably not working properly.