# Setting up implicit solvers to beat the performance of explicit solvers

**URL:** <https://discourse.julialang.org/t/setting-up-implicit-solvers-to-beat-the-performance-of-explicit-solvers/102392>\
**Category:** Numerics\
**Tags:** diffeq, forwarddiff\
**Created:** [August 2, 2023, 11:57am UTC](https://discourse.julialang.org/t/setting-up-implicit-solvers-to-beat-the-performance-of-explicit-solvers/102392 "2023-08-02T11:57:16Z")\
**Posts on this page:** 11\
**Page:** 1

<div class="post-metadata">

**Author:** ![JordiBolibar](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jordibolibar/32/24307_2.png) [@JordiBolibar](https://discourse.julialang.org/u/JordiBolibar)\
**Post date:** [August 2, 2023, 11:57am UTC](https://discourse.julialang.org/t/setting-up-implicit-solvers-to-beat-the-performance-of-explicit-solvers/102392/1 "2023-08-02T11:57:16Z")

</div>

We are trying to benchmark implicit solvers for an [ice flow diffusivity PDE](https://github.com/ODINN-SciML/iceflow_sandbox/blob/main/scripts/1D_SIA.jl). So far we’ve stuck to explicit ones, but now we would like to cover all of them. Since many implicit solvers from DifferentialEquations.jl require AD, we’re currently having the following issue:

```julia
First call to automatic differentiation for the Jacobian
failed. This means that the user `f` function is not compatible
with automatic differentiation. Methods to fix this include:

1. Turn off automatic differentiation (e.g. Rosenbrock23() becomes
   Rosenbrock23(autodiff=false)). More details can befound at
   https://docs.sciml.ai/DiffEqDocs/stable/features/performance_overloads/
2. Improving the compatibility of `f` with ForwardDiff.jl automatic 
   differentiation (using tools like PreallocationTools.jl). More details
   can be found at https://docs.sciml.ai/DiffEqDocs/stable/basics/faq/#Autodifferentiation-and-Dual-Numbers
3. Defining analytical Jacobians. More details can be
   found at https://docs.sciml.ai/DiffEqDocs/stable/types/ode_types/#SciMLBase.ODEFunction

```

Since our `f` function is really optimized to avoid almost all memory allocations, I guess the problem comes from the usual fact that AD doesn’t like mutation. We have tried deactivating AD, but it’s terribly slow.

So my question is: what should one do in this case? The documentation is really scarce, and we haven’t found any clear examples on how to work around this issue. The code can be found [here](https://github.com/ODINN-SciML/iceflow_sandbox/blob/main/scripts/1D_SIA.jl), and [here](https://github.com/ODINN-SciML/iceflow_sandbox/tree/main)’s the repository with all the benchmarks (with explicit solvers so far).

Thanks in advance!

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [August 2, 2023, 12:05pm UTC](https://discourse.julialang.org/t/setting-up-implicit-solvers-to-beat-the-performance-of-explicit-solvers/102392/2 "2023-08-02T12:05:39Z")

</div>

the easy answer is to pass autodiff=false to the solver. it will then use finite different

---

<div class="post-metadata">

**Author:** ![JordiBolibar](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jordibolibar/32/24307_2.png) [@JordiBolibar](https://discourse.julialang.org/u/JordiBolibar)\
**Post date:** [August 2, 2023, 12:12pm UTC](https://discourse.julialang.org/t/setting-up-implicit-solvers-to-beat-the-performance-of-explicit-solvers/102392/3 "2023-08-02T12:12:58Z")

</div>

Yes, as I mentioned, we tried this and it’s horribly slow. We’re looking into implicit solvers to beat the performance of our best explicit solver.

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [August 2, 2023, 12:28pm UTC](https://discourse.julialang.org/t/setting-up-implicit-solvers-to-beat-the-performance-of-explicit-solvers/102392/4 "2023-08-02T12:28:23Z")

</div>

It looks like the error message suggests improving compatiblity with ForwardDiff.jl. The problem then is likely that you allocate float arrays instead of generically typed arrays. How do your allocations look? Can you manually compute the Jacobian using ForwardDiff?

> <https://github.com/ODINN-SciML/iceflow_sandbox/blob/ad205e17e6141520c1d66c27f8228566f2bca182/scripts/1D_SIA.jl#L106>

It does indeed look like you hard code types everywhere, don’t do that 😉

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [August 2, 2023, 12:38pm UTC](https://discourse.julialang.org/t/setting-up-implicit-solvers-to-beat-the-performance-of-explicit-solvers/102392/5 "2023-08-02T12:38:57Z")

</div>

You can use the implicit solvers without AD. For example pass `FBDF(autodiff=false)` as the solver.

---

<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:** [August 2, 2023, 1:20pm UTC](https://discourse.julialang.org/t/setting-up-implicit-solvers-to-beat-the-performance-of-explicit-solvers/102392/6 "2023-08-02T13:20:21Z")

</div>

> [@JordiBolibar](#):
>
> Yes, as I mentioned, we tried this and it’s horribly slow. We’re looking into implicit solvers to beat the performance of our best explicit solver.

Note that finite difference of the solver does not make a huge impact on performance unless it’s a Rodas method. It’s like a 2x-4x performance thing because it’s forward mode AD and has nothing to do with the adjoints of the back solve which is what really matters.

That said, for PDEs the issue is that Rosenbrock23 is a bad idea. As the docs mention, it’s not a method that scales to larger systems well. Did you try FBDF or KenCarp47? Those are more sensible algorithms for large equations. And then when optimizing that, you should look into the tutorial on handling large systems:

> **[Solving Large Stiff Equations · DifferentialEquations.jl](https://docs.sciml.ai/DiffEqDocs/stable/tutorials/advanced_ode_example/)**
>
> Documentation for DifferentialEquations.jl.

Setting up sparse Jacobians and iLU or multigrid preconditioning is such a huge boost that you should always do it for PDEs.

---

<div class="post-metadata">

**Author:** ![LucilleGimenes](https://avatars.discourse-cdn.com/v4/letter/l/b19c9b/32.png) [@LucilleGimenes](https://discourse.julialang.org/u/LucilleGimenes)\
**Post date:** [August 17, 2023, 9:15am UTC](https://discourse.julialang.org/t/setting-up-implicit-solvers-to-beat-the-performance-of-explicit-solvers/102392/7 "2023-08-17T09:15:26Z")

</div>

I am working with @JordiBolibar on this ice flow diffusivity PDE issue. Actually, the latest version of the `f` function that we use can be found [here](https://github.com/lucillegimenes/oggm/blob/master/oggm/core/SIA1D_utils.jl) and the latest benchmark is available [here](https://github.com/lucillegimenes/iceflow_sandbox).

We tried using FBDF and KenCarp47, and it is still at least 10 times slower than when using explicit solvers.

Also, when trying to set up sparse Jacobians with  
`u0 = iceflow_model.H`  
`du0 = zeros(Float64,iceflow_model.nx)`  
`jac_sparsity = Symbolics.jacobian_sparsity((du, u) -> SIA1D!(du, u, iceflow_model, 0.0), du0,u0)`  
we get the following error :

```julia
 MethodError: no method matching SIA1D!(::Vector{Num}, ::Vector{Num}, ::SIA1Dmodel{Float64, Int64}, ::Float64)
Closest candidates are:
  SIA1D!(!Matched::Vector{Float64}, !Matched::Vector{Float64}, ::SIA1Dmodel{Float64, Int64}, ::Float64) at ~/oggm/oggm/core/SIA1D_utils.jl:14
Stacktrace:
 [1] (::var"#10#13"{SIA1Dmodel{Float64, Int64}})(du::Vector{Num}, u::Vector{Num})
   @ Main ~/oggm/oggm/core/SIA1D_utils.jl:174
 [2] jacobian_sparsity(::var"#10#13"{SIA1Dmodel{Float64, Int64}}, ::Vector{Float64}, ::Vector{Float64}; kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
   @ Symbolics ~/.julia/packages/Symbolics/BQlmn/src/diff.jl:584
 [3] jacobian_sparsity(::Function, ::Vector{Float64}, ::Vector{Float64})
   @ Symbolics ~/.julia/packages/Symbolics/BQlmn/src/diff.jl:579

```

---

<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:** [August 17, 2023, 9:27am UTC](https://discourse.julialang.org/t/setting-up-implicit-solvers-to-beat-the-performance-of-explicit-solvers/102392/8 "2023-08-17T09:27:44Z")

</div>

Your caches force `Float64`, so that’s of course going to fail on AD and sparsity detection. That’s what PreallocationTools.jl is for:

> **[GitHub - SciML/PreallocationTools.jl: Tools for building non-allocating...](https://github.com/SciML/PreallocationTools.jl)**
>
> Tools for building non-allocating pre-cached functions in Julia, allowing for GC-free usage of automatic differentiation in complex codes - GitHub - SciML/PreallocationTools.jl: Tools for building ...

To see if this is a direction you should go, did you try using GMRES without a preconditioner and see how that does?

---

<div class="post-metadata">

**Author:** ![LucilleGimenes](https://avatars.discourse-cdn.com/v4/letter/l/b19c9b/32.png) [@LucilleGimenes](https://discourse.julialang.org/u/LucilleGimenes)\
**Post date:** [August 24, 2023, 2:54pm UTC](https://discourse.julialang.org/t/setting-up-implicit-solvers-to-beat-the-performance-of-explicit-solvers/102392/9 "2023-08-24T14:54:37Z")

</div>

Thanks for your message; by using PreallocationTools.jl we were able to use the AD of implicit solvers such as KenCarp47 or FBDF.

However, they still perform at least 3 times slower that the explicit solver we were using previously (RDPK3Sp35), even when setting up sparse jacobians and using GMRES (i.e adding `linsolve = KrylovJL_GMRES()` as a solver agument if that’s what you meant).

---

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [August 24, 2023, 3:52pm UTC](https://discourse.julialang.org/t/setting-up-implicit-solvers-to-beat-the-performance-of-explicit-solvers/102392/10 "2023-08-24T15:52:54Z")

</div>

You need to tune `KrylovJL_GMRES()`, pass in verbose mode and see how fast it converges

---

<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:** [August 24, 2023, 4:16pm UTC](https://discourse.julialang.org/t/setting-up-implicit-solvers-to-beat-the-performance-of-explicit-solvers/102392/11 "2023-08-24T16:16:45Z")

</div>

Are you sure the equation is stiff? What makes you think so?

If it’s on the edge, did you try `ROCK2()` or `ROCK4()`? There are some PDE cases which have “not tiny but not large real valued eigenvalues” (i.e. a Laplacian) where this is the most efficient solver.

> [@LucilleGimenes](#):
>
> (i.e adding `linsolve = KrylovJL_GMRES()` as a solver agument if that’s what you meant).

What preconditioner? Did you ilu and tune the cutoff?

> **[Solving Large Stiff Equations · DifferentialEquations.jl](https://docs.sciml.ai/DiffEqDocs/stable/tutorials/advanced_ode_example/)**
>
> Documentation for DifferentialEquations.jl.