OrdinaryDiffEq throws error when using sparse jacobian without specifying linear solver

Specific Domains Modelling & Simulations
1

Hi All,

I nearly successfully created a sparse ODEProblem by using Symbolics.jacobian_sparsity, but the solver throws an error unless I use TRBDF2(), (or one of the other linear solvers), coming from OrdinaryDiffEq (here):

ERROR: AssertionError: J.colptr == W.colptr

I constructed the ODEProblem with the sparse Jacobian from the simpler ODEProblem, in the following way, where the error occurs only on executing the last line.

parameters = loadparameters(
    test_modelparamCSV,
    test_computeparamCSV,
)
dudt = enclosethetimedifferential(parameters)
IC = ones(54) # IC = makebellcurveIC(parameters)
odeprob = ODEProblem(
    dudt,
    IC,
    (0, 2.1),
    parameters.prior,
);
du0 = copy(odeprob.u0);
jac_sparsity = Symbolics.jacobian_sparsity(
    (du,u)->dudt(du, u, parameters.prior, 0.0),
    du0,
    odeprob.u0,
);
f = ODEFunction(
    dudt;
    jac_prototype=float.(jac_sparsity),
);
sparseodeprob = ODEProblem(
    f,
    odeprob.u0,
    (0, 2.1),
    parameters.prior,
);

solve(odeprob);
solve(odeprob, TRBDF2());

solve(sparseodeprob, TRBDF2());
solve(sparseodeprob, KenCarp47(linsolve=KLUFactorization()));
solve(sparseodeprob, KenCarp47(linsolve=KrylovJL_GMRES()));
solve(sparseodeprob); # ERROR : AssertionError: J.colptr == W.colptr

Is this a problem with solve using a bad solver, or have I misspecified the Jacobian somehow?

Thanks!

2

I cannot run this code since that (and some other variables) are not defined.

3

My apologies, I’m a bit new to this, I thought I’d be expecting too much to think anybody would run the code rather than read it over. But the Julia community is great - and you’re a great example of that @ChrisRackauckas, and thanks for all the tutorials and inspiring work.

The code below is sufficient to reproduce the error on my Julia 1.8.0. I’ve just now found that if I let the boundary conditions stay zero, then there’s no issue. Perhaps changing the boundary conditions during the solving messes with the Jacobian?

using DifferentialEquations, DiffEqOperators, SparseArrays
using Parameters
using ComponentArrays
using Symbolics

function enclosethetimedifferential(parameters::NamedTuple)::Function
    @info "Enclosing the time differential"

    @unpack Δr, r_space, countorderapprox = parameters.compute
    countdiscretizationsteps = length(r_space)

    ∇ = CenteredDifference(1, countorderapprox, Δr, countdiscretizationsteps)
    Δ = CenteredDifference(2, countorderapprox, Δr, countdiscretizationsteps)
    # To concretize as arrays
    tmpbc = NeumannBC((3.14, 1.0), Δr)
    Δ, _ = Array(Δ * tmpbc)
    ∇, _ = Array(∇ * tmpbc)
    bc_x = zeros(Real, countdiscretizationsteps)
    bc_xx = zeros(Real, countdiscretizationsteps)

    function timedifferentialclosure!(du, u, p, t)
        @unpack (α, D, v, k_p, V_c, Q_l, Q_r, V_b,
        S, Lm, Dm, V_v,) = p

        c = u[1:end - 3]
        c_v = u[end - 2]
        c_c = u[end - 1]
        c_b = u[end]

        J_B0 = (Dm/Lm) * (α * c_v - c[1])
        J_BL = (Dm/Lm) * (c[end] - α * c_c)
        grad_0 = (v./D) .* c[1] .- J_B0 ./ D
        grad_L = (v./D) .* c[end] .- J_BL ./ D

        bc_x[1] = grad_0 / 2
        bc_x[end] = grad_L / 2
        grad_c = ∇ * c + bc_x

        bc_xx[1] = -grad_0 / Δr
        bc_xx[end] = grad_L / Δr
        Lap_c = Δ * c + bc_xx

        C = sum(Δr .* S * (k_p * (c .- c_b)))

        dc_dt = D * Lap_c - v * grad_c .- k_p * (c .- c_b)
        du[1:end-3] = dc_dt[1:end]

        dcv_dt = -S * J_B0 / V_v - (Q_l / V_v) * c_v
        du[end-2] = dcv_dt

        dcc_dt = S * α * J_BL / V_c + (Q_l / V_c) * c_v - (Q_l/V_c) * c_c
        du[end-1] = dcc_dt

        dcb_dt = (Q_l/V_b) * c_c + C / V_b
        du[end] = dcb_dt
    end

    return timedifferentialclosure!
end

prior = ComponentArray(;
    α = 0.2,
    D = 0.46,
    v = 0.0,
    k_p = 0.0,
    V_c = 18,
    Q_l = 20,
    Q_r = 3.6,
    V_b = 1490,
    S = 52,
    Lm = 0.05,
    Dm = 0.046,
    V_v = 18.0,
)

r_space = collect(range(start=0.0, stop=2.0, length=15))
computeparams = (
    Δr = r_space[2],
    r_space = r_space,
    countorderapprox = 2,
)
parameters = (
    prior=prior,
    compute=computeparams,
)

dudt = enclosethetimedifferential(parameters)
IC = ones(length(r_space) + 3)
odeprob = ODEProblem(
    dudt,
    IC,
    (0, 2.1),
    parameters.prior,
);
du0 = copy(odeprob.u0);
jac_sparsity = Symbolics.jacobian_sparsity(
    (du,u)->dudt(du, u, parameters.prior, 0.0),
    du0,
    odeprob.u0,
);
f = ODEFunction(
    dudt;
    jac_prototype=float.(jac_sparsity),
);
sparseodeprob = ODEProblem(
    f,
    odeprob.u0,
    (0, 2.1),
    parameters.prior,
);

#@btime
solve(odeprob);
solve(odeprob, TRBDF2());

solve(sparseodeprob, TRBDF2());
#solve(sparseodeprob, KenCarp47(linsolve=KLUFactorization()));
#solve(sparseodeprob, KenCarp47(linsolve=KrylovJL_GMRES()));
solve(sparseodeprob);
1 Like
4

Thanks, this is now fixed by Handle sparsity patterns which are specified as missing some diagonal by ChrisRackauckas · Pull Request #1820 · SciML/OrdinaryDiffEq.jl · GitHub. The issue is that, as documented, the sparsity pattern should be the sparsity of W = M - gamma*J for arbitrary gamma. But that’s a bit hard to always remember, so that PR adds any of the missing mass matrix parts (usually just I, so the diagonal) automatically, which fixes the issue.

1 Like