# OrdinaryDiffEq throws error when using sparse jacobian without specifying linear solver
**URL:** https://discourse.julialang.org/t/ordinarydiffeq-throws-error-when-using-sparse-jacobian-without-specifying-linear-solver/89141
**Category:** Modelling & Simulations
**Tags:** ordinarydiffeq
**Created:** [October 23, 2022, 2:15pm UTC](https://discourse.julialang.org/t/ordinarydiffeq-throws-error-when-using-sparse-jacobian-without-specifying-linear-solver/89141 "2022-10-23T14:15:06Z")
**Posts on this page:** 4
**Page:** 1
### Author:  [@parb](https://discourse.julialang.org/u/parb)
#### Post date: [October 23, 2022, 2:15pm UTC](https://discourse.julialang.org/t/ordinarydiffeq-throws-error-when-using-sparse-jacobian-without-specifying-linear-solver/89141/1 "2022-10-23T14:15:06Z")
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](https://github.com/SciML/OrdinaryDiffEq.jl/blob/643f96145354c5a522f8ec6c091202eacfbeec4a/src/derivative_utils.jl#L550)):
```julia
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.
```julia
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!
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [October 23, 2022, 2:49pm UTC](https://discourse.julialang.org/t/ordinarydiffeq-throws-error-when-using-sparse-jacobian-without-specifying-linear-solver/89141/2 "2022-10-23T14:49:37Z")
> [@parb](#):
>
> `test_modelparamCSV`
I cannot run this code since that (and some other variables) are not defined.
---
### Author:  [@parb](https://discourse.julialang.org/u/parb)
#### Post date: [October 25, 2022, 7:00am UTC](https://discourse.julialang.org/t/ordinarydiffeq-throws-error-when-using-sparse-jacobian-without-specifying-linear-solver/89141/3 "2022-10-25T07:00:12Z")
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?
```julia
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);
```
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [December 24, 2022, 7:14am UTC](https://discourse.julialang.org/t/ordinarydiffeq-throws-error-when-using-sparse-jacobian-without-specifying-linear-solver/89141/4 "2022-12-24T07:14:49Z")
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](https://github.com/SciML/OrdinaryDiffEq.jl/pull/1820). 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.