Error with jac_prototype in Mass Matrix Form of ODEFunction (svd! on Sparse Matrix)

Specific Domains Numerics
1

Hi all,

I’m solving a DAE system using DifferentialEquations.jl in the mass matrix form. The issue arises when I include jac_prototype in the ODEFunction. If I define ODEFunction without jac_prototype, everything works fine. However, when I include jac_prototype, which is a sparse matrix, I get the following error during solve:

ERROR: MethodError: no method matching svd!(::SparseMatrixCSC{Float64, Int64}; full::Bool, alg::LinearAlgebra.DivideAndConquer)
The function `svd!` exists, but no method is defined for this combination of argument types.

This is my MWE

using LinearAlgebra
using SparseArrays
using ADTypes, SparseConnectivityTracer
# import DifferentialEquations as DE
using DifferentialEquations

function laplacian(n)
    main_diag = -2.0 * ones(n)
    off_diag = ones(n - 1)
    return spdiagm(0 => main_diag, -1 => off_diag, 1 => off_diag)
end

function landau_khalatnikov_debug!(du, u, p)
    poisson_residue_debug!(du, u, p)
    return nothing
end

function poisson_residue_debug!(du, u, p)
    l1, l2 = laplacian.(p)
    lap = kron(l2, I(p[1])) + kron(I(p[2]), l1)
    du, u = view(du, :), view(u, :)
    mul!(du, lap, u)
    return nothing
end

function system_residue_debug!(du, u, p, t)
    Ny, Nx, Nyfed = p
    num_nodes = Ny * Nx
    num_edges = Nyfed * Nx
    fill!(du, zero(eltype(u)))
    Py = view(u, 1:num_edges)
    dPy = view(du, 1:num_edges)
    ϕ = view(u, 1+num_edges:num_edges+num_nodes)
    dϕ = view(du, 1+num_edges:num_edges+num_nodes)

    Py = reshape(Py, Nyfed, Nx)
    dPy = reshape(dPy, Nyfed, Nx)
    ϕ = reshape(ϕ, Ny, Nx)
    dϕ = reshape(dϕ, Ny, Nx)

    landau_khalatnikov_debug!(dPy, Py, (Nyfed, Nx))
    poisson_residue_debug!(dϕ, ϕ, (Ny, Nx))
    return nothing
end

Nx, Ny, Nyfe = 32, 32, 13
params = Ny, Nx, Nyfe + 1

ϕ = zeros(Ny, Nx)
Py = rand(Nyfe + 1, Nx)
u0 = rand(length(Py) + length(ϕ))
du = similar(u0)

# Jacobian sparsity
system_residue_debug!(du, u0, params, 0.0)
jac_sparsity_system = ADTypes.jacobian_sparsity((du, u) -> system_residue_debug!(du, u, params, 0.0), similar(u0), copy(u0), TracerSparsityDetector())
jac_prototype = float.(jac_sparsity_system)

mass_matrix = Diagonal(vcat(ones(length(Py)), zeros(length(ϕ))))
ode_func = ODEFunction(system_residue_debug!, jac_prototype=jac_prototype, mass_matrix=mass_matrix)
ode_prob = ODEProblem(ode_func, u0, (0.0, 1e-19), params)
sol = solve(ode_prob, Rodas5(), reltol=1e-8, abstol=1e-8)

ode_func = ODEFunction(system_residue_debug!, mass_matrix=mass_matrix)
ode_prob = ODEProblem(ode_func, u0, (0.0, 1e-19), params)
sol = solve(ode_prob, Rodas5(), reltol=1e-8, abstol=1e-8)

The documentation suggests that this should be possible with Rodas5() algorithm.

Any help is appreciated.

1 Like
2

Can you open an issue for this one? It ended up being weirder than expected. The solution for now is to force a different nonlinear solver during initialization, but that shouldn’t be required. I believe @Oscar_Smith is looking into this but because it’s not in the bug tracker it’s not being tracked and will get lost.

1 Like
3

So the specific problem is an issue with solving a nonlinear problem with Broyden(; init_jacobian = Val(:true_jacobian), autodiff=AutoFiniteDiff()).

Specifically, the MWE is

using LinearAlgebra, SparseArrays, NonlinearSolveQuasiNewton, ForwardDiff
function f(du, u,p)
    @. du = u^2 - 2
    du[1] += u[4]
end

prob = NonlinearProblem(NonlinearFunction(f; jac_prototype=sparse(ones(4, 4))), rand(4))
solve(prob, Broyden(; init_jacobian = Val(:true_jacobian)))
4

I have opened the issue here.

Thanks, I will try to set the nonlinear solvers as NewtonRaphson and see if I can make it work.

I also read up on this, and it seems that solving mass matrix ODE is same as normal ODE, it is just that effectively during the stage computation of an implicit method, coefficent of y_next is modified from (I - γhJ) to (M - γhJ). So just for educational purposes I am going to implement it.