Preconditioner for Sundials IDA solver in DifferentialEquations.jl

Specific Domains Numerics
,
1

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)) 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)).
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

2 Likes
2

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

3

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.

4

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)). The solver fails because it reaches the maximum number of iterations. Thank you very much in advance.

# 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
5

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

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.

4 Likes