# Preconditioner for Sundials IDA solver in DifferentialEquations.jl

**URL:** <https://discourse.julialang.org/t/preconditioner-for-sundials-ida-solver-in-differentialequations-jl/105124>\
**Category:** Numerics\
**Tags:** package, differentialequation\
**Created:** [October 18, 2023, 12:53pm UTC](https://discourse.julialang.org/t/preconditioner-for-sundials-ida-solver-in-differentialequations-jl/105124 "2023-10-18T12:53:17Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![jespfra](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jespfra/32/202556_2.png) [@jespfra](https://discourse.julialang.org/u/jespfra)\
**Post date:** [October 18, 2023, 12:53pm UTC](https://discourse.julialang.org/t/preconditioner-for-sundials-ida-solver-in-differentialequations-jl/105124/1 "2023-10-18T12:53:17Z")

</div>

Hello Everyone,

In DifferentialEquations.jl, I have not found an example of how to should set a preconditioner for the Sundials IDA solver. I have found examples of how to set it up for ODE systems for the CVODE\_BDF solver ( [Solving Large Stiff Equations · DifferentialEquations.jl (sciml.ai)](https://docs.sciml.ai/DiffEqDocs/latest/tutorials/advanced_ode_example/#stiff)) but I am also aware that the arguments for the preconditioner and the psetup are different for the IDA solver ( [DAE Problems · DifferentialEquations.jl (sciml.ai)](https://docs.sciml.ai/DiffEqDocs/latest/types/dae_types/)).  
Therefore, if any of you guys have a (simple) test case that demonstrate how to setup a precilu and psetupilu for the IDA solver, I would be grateful!

Thanks in advance.

Best,  
Jesper

---

<div class="post-metadata">

**Author:** ![ufechner7](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ufechner7/32/51363_2.png) [@ufechner7](https://discourse.julialang.org/u/ufechner7)\
**Post date:** [October 18, 2023, 1:26pm UTC](https://discourse.julialang.org/t/preconditioner-for-sundials-ida-solver-in-differentialequations-jl/105124/2 "2023-10-18T13:26:42Z")

</div>

Good question! Would also be interested to understand how to use the preconditioner…

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [October 18, 2023, 2:01pm UTC](https://discourse.julialang.org/t/preconditioner-for-sundials-ida-solver-in-differentialequations-jl/105124/3 "2023-10-18T14:01:05Z")

</div>

In theory it’s all the same as with CVODE, so that tutorial should be all that’s needed. Looking at the source code, it is all setup and wrapped. That said, all of the tests are on CVODE, not with IDA, and so something could be missing. Can you setup an MWE? I’d be happy to dig into it.

---

<div class="post-metadata">

**Author:** ![jespfra](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jespfra/32/202556_2.png) [@jespfra](https://discourse.julialang.org/u/jespfra)\
**Post date:** [October 19, 2023, 9:43am UTC](https://discourse.julialang.org/t/preconditioner-for-sundials-ida-solver-in-differentialequations-jl/105124/4 "2023-10-19T09:43:32Z")

</div>

Hi Again,

Great to hear that it has interest for others too!

Below, I have set up a case for the Rober equation where I first solve it using default ODE solver, then solve it using the CVODE preconditioning, then solve it using IDA without preconditioning and finally a failed attempt to solve the equations using IDA preconditioning according to the documentation ( [Sundials.jl · DifferentialEquations.jl (sciml.ai)](https://docs.sciml.ai/DiffEqDocs/stable/api/sundials/)). The solver fails because it reaches the maximum number of iterations. Thank you very much in advance.

```julia
# The rober example
#Load packages
using DifferentialEquations
using Sundials
using IncompleteLU
using SparseArrays
using Plots
using LinearAlgebra

################## Solving ODE ##################

#ODE repressentation
function rober!(du, u, p, t)
    y₁, y₂, y₃ = u
    k₁, k₂, k₃ = p
    du[1] = -k₁ * y₁ + k₃ * y₂ * y₃
    du[2] = k₁ * y₁ - k₂ * y₂^2 - k₃ * y₂ * y₃
    du[3] = k₂ * y₂^2
    nothing
end

u0 = [1.0; 0.0; 0.0]
tspan = (0.0, 1e4)
p = [0.04, 3e7, 1e4]
prob = ODEProblem(rober!, u0, tspan,p)
sol = solve(prob)
p1 = plot();
plot!(p1,sol);

################## Solving ODE with CVODE preconditioning ##################
#Jacobian
function rober_jac!(J, u, p, t)
    y₁, y₂, y₃ = u
    k₁, k₂, k₃ = p
    J[1, 1] = k₁ * -1
    J[2, 1] = k₁
    J[3, 1] = 0
    J[1, 2] = y₃ * k₃
    J[2, 2] = y₂ * k₂ * -2 + y₃ * k₃ * -1
    J[3, 2] = y₂ * 2 * k₂
    J[1, 3] = k₃ * y₂
    J[2, 3] = k₃ * y₂ * -1
    J[3, 3] = 0
    return J
end

#psetup
function psetupilu(p, t, u, du, jok, jcurPtr, gamma)
    if !jok
        rober_jac!(jaccache, u, p, t) #Sundials handling
        jcurPtr[] = true

        # W = I - gamma*J
        @. W = -gamma * jaccache
        idxs = diagind(W)
        @. @view(W[idxs]) = @view(W[idxs]) + 1

        # Build preconditioner on W
        preccache[] = ilu(W, τ = 5.0)
    end
end

# preconditioner 
function precilu(z,r,p,t,y,fy,gamma,delta,lr)
    ldiv!(z,preccache[],r)
end

jacc = sparse(zeros(length(u0),length(u0)))
#Define constants
const jaccache = rober_jac!(jacc, u0, p, 0.0)
const W = I - 1.0*jaccache
prectmp = ilu(W, τ = 50.0)
const preccache = Ref(prectmp)

fun = ODEFunction(rober!) #
prob = ODEProblem(fun, u0, tspan,p) #
sol1 = solve(prob,CVODE_BDF(linear_solver = :GMRES, prec = precilu, psetup = psetupilu, prec_side = 1))
plot!(p1,sol1)

################## DAE repressentation ##################
function rober_DAE!(out,RHS, u, p, t)
    k1, k2, k3 = p

    out[1] = -k1 * u[1] + k3 * u[2]*u[3] - RHS[1]
    out[2] = k1 * u[1] - k2 * u[2]^2 -k3 * u[2]*u[3] - RHS[2]
    out[3] = k2 * u[2]^2 - RHS[3]
    nothing
end

differential_vars = [true,true,true]

du0 = zeros(Float64,3); RHS = zeros(Float64,3)
rober_DAE!(du0,RHS, u0, p, 0.0)

prob = DAEProblem(rober_DAE!, du0, u0, tspan,p, differential_vars = differential_vars)
sol2 = solve(prob,IDA())
plot!(p1,sol2)

################## Preconditioning with IDA solver ##################
#prec
function preciluDAE(z,r,p,t,y,fy,gamma,delta,lr)
    ldiv!(z,preccacheDAE[],r)
end

#psetup
function psetupiluDAE(p,t,resid,u,du, gamma)
    
    rober_jac!(jaccacheDAE, u, p, t) #
    # jcurPtr[] = true

    # W = I - gamma*J
    @. W = -gamma * jaccacheDAE
    idxs = diagind(W)
    @. @view(W[idxs]) = @view(W[idxs]) + 1

    # Build preconditioner on W
    preccacheDAE[] = ilu(W, τ = 5.0)
end

#Defining constants
jacc = sparse(zeros(length(u0),length(u0)))
const jaccacheDAE = rober_jac!(jacc, u0, p, 0.0)
const W = I - 1.0*jaccacheDAE
prectmp = ilu(W, τ = 50.0)
const preccacheDAE = Ref(prectmp)

prob = DAEProblem(rober_DAE!, du0, u0, tspan,p, differential_vars = differential_vars)
sol3 = solve(prob,IDA(linear_solver = :GMRES,prec = preciluDAE, psetup = psetupiluDAE)) #Fails with retcode: MaxIters
```

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [April 6, 2024, 10:11pm UTC](https://discourse.julialang.org/t/preconditioner-for-sundials-ida-solver-in-differentialequations-jl/105124/5 "2024-04-06T22:11:24Z")

</div>

@jespfra sorry it took awhile to get back to this. The issue was with your Jacobain definition. The follow code works:

```julia
using Sundials
using IncompleteLU
using SparseArrays
using Plots
using LinearAlgebra

################## Solving ODE ##################

#ODE repressentation
function rober!(du, u, p, t)
    y₁, y₂, y₃ = u
    k₁, k₂, k₃ = p
    du[1] = -k₁ * y₁ + k₃ * y₂ * y₃
    du[2] = k₁ * y₁ - k₂ * y₂^2 - k₃ * y₂ * y₃
    du[3] = k₂ * y₂^2
    nothing
end

u0 = [1.0; 0.0; 0.0]
tspan = (0.0, 1e4)
p = [0.04, 3e7, 1e4]
prob = ODEProblem(rober!, u0, tspan,p)
sol = solve(prob)
p1 = plot();
plot!(p1,sol);

################## Solving ODE with CVODE preconditioning ##################
#Jacobian
function rober_jac!(J, u, p, t)
    y₁, y₂, y₃ = u
    k₁, k₂, k₃ = p
    J[1, 1] = k₁ * -1
    J[2, 1] = k₁
    J[3, 1] = 0
    J[1, 2] = y₃ * k₃
    J[2, 2] = y₂ * k₂ * -2 + y₃ * k₃ * -1
    J[3, 2] = y₂ * 2 * k₂
    J[1, 3] = k₃ * y₂
    J[2, 3] = k₃ * y₂ * -1
    J[3, 3] = 0
    return J
end

#psetup
function psetupilu(p, t, u, du, jok, jcurPtr, gamma)
    if !jok
        rober_jac!(jaccache, u, p, t) #Sundials handling
        jcurPtr[] = true

        # W = I - gamma*J
        @. W = -gamma * jaccache
        idxs = diagind(W)
        @. @view(W[idxs]) = @view(W[idxs]) + 1

        # Build preconditioner on W
        preccache[] = ilu(W, τ = 5.0)
    end
end

# preconditioner 
function precilu(z,r,p,t,y,fy,gamma,delta,lr)
    ldiv!(z,preccache[],r)
end

jacc = sparse(zeros(length(u0),length(u0)))
#Define constants
const jaccache = rober_jac!(jacc, u0, p, 0.0)
const W = I - 1.0*jaccache
prectmp = ilu(W, τ = 50.0)
const preccache = Ref(prectmp)

fun = ODEFunction(rober!) #
prob = ODEProblem(fun, u0, tspan,p) #
sol1 = solve(prob,CVODE_BDF(linear_solver = :GMRES, prec = precilu, psetup = psetupilu, prec_side = 1))
plot!(p1,sol1)

################## DAE repressentation ##################
function rober_DAE!(out,RHS, u, p, t)
    k1, k2, k3 = p

    out[1] = -k1 * u[1] + k3 * u[2]*u[3] - RHS[1]
    out[2] = k1 * u[1] - k2 * u[2]^2 -k3 * u[2]*u[3] - RHS[2]
    out[3] = k2 * u[2]^2 - RHS[3]
    nothing
end

differential_vars = [true,true,true]

du0 = zeros(Float64,3); RHS = zeros(Float64,3)
rober_DAE!(du0,RHS, u0, p, 0.0)

prob = DAEProblem(rober_DAE!, du0, u0, tspan,p, differential_vars = differential_vars)
sol2 = solve(prob,IDA())
plot!(p1,sol2)

################## Preconditioning with IDA solver ##################
#prec

function psetupilu(p, t, u, du, jok, jcurPtr, gamma)
  if !jok
      rober_jac!(jaccache, u, p, t) #Sundials handling
      jcurPtr[] = true

      # W = I - gamma*J
      @. W = -gamma * jaccache
      idxs = diagind(W)
      @. @view(W[idxs]) = @view(W[idxs]) + 1

      # Build preconditioner on W
      preccache[] = ilu(W, τ = 5.0)
  end
end

function preciluDAE(z,r,p,t,y,fy,gamma,delta,lr)
  ldiv!(z,preccacheDAE[],r)
end

#psetup

function psetupilu(p, t, u, du, jok, jcurPtr, gamma)
  if !jok
      rober_jac!(jaccache, u, p, t) #Sundials handling
      jcurPtr[] = true

      # W = I - gamma*J
      @. W = -gamma * jaccache
      idxs = diagind(W)
      @. @view(W[idxs]) = @view(W[idxs]) + 1

      # Build preconditioner on W
      preccache[] = ilu(W, τ = 5.0)
  end
end

function psetupiluDAE(p,t,resid,u,du, gamma)

    rober_jac!(jaccacheDAE, u, p, t) #
    # jcurPtr[] = true

    # W = df/du + gamma * df/d(du)
    @. W = jaccacheDAE
    idxs = diagind(W)
    @. @view(W[idxs]) = -gamma

    # Build preconditioner on W
    preccacheDAE[] = ilu(W, τ = 5.0)
end

#Defining constants
jacc = sparse(zeros(length(u0),length(u0)))
const jaccacheDAE = rober_jac!(jacc, u0, p, 0.0)
const W = I - 1.0*jaccacheDAE
prectmp = ilu(W, τ = 50.0)
const preccacheDAE = Ref(prectmp)

prob = DAEProblem(rober_DAE!, du0, u0, tspan,p, differential_vars = differential_vars)
sol3 = solve(prob,IDA(linear_solver = :GMRES,prec = preciluDAE, psetup = psetupiluDAE))

```

IDA defines the W matrix differently than CVODE. In particular, it’s of the form `W = df/du + gamma * df/d(du)`. You were building `W = df/d(du) - gamma*df/du`, which is what matches CVODE and is thus why it was diverging.