Convergence issues with Theseus DIRK and IMEX solvers

[ranocha]

I enabled the option to pass more arguments to the Newton solver in #79. Trying to reduce the relative tolerance to fix the erratic error behavior described in #78, I encountered convergence issues like

julia> using Theseus

julia> ode = ODEProblem([1.0, 0.0, 1.0], (0.0, 1.0)) do du, u, p, t
           du[1] = -u[2]
           du[2] = u[1]
           du[3] = -u[3]
           return nothing
       end

julia> solve(ode, Theseus.DIRK43(); 
             adaptive = false, dt = 0.25, 
             newton_tol_abs = 1.0e-10, newton_tol_rel = 0.0, 
             newton_max_niter = 10)
┌ Warning: Newton did not converge
│   stats = (solved = false, stats = Ariadne.Stats(11, 8, 4.352315131535924e-10), t = 1.3791e-5)
│   integrator.t = 0.0
└ @ Theseus ~/.julia/dev/Ariadne/libs/Theseus/src/dirk/dirk.jl:273
ERROR: Newton did not converge

julia> solve(ode, Theseus.DIRK43(); 
             adaptive = false, dt = 0.25, 
             newton_tol_abs = 1.0e-10, newton_tol_rel = 0.0, 
             newton_max_niter = 100)
┌ Warning: Newton did not converge
│   stats = (solved = false, stats = Ariadne.Stats(101, 8, 4.352315131535924e-10), t = 0.000106334)
│   integrator.t = 0.0
└ @ Theseus ~/.julia/dev/Ariadne/libs/Theseus/src/dirk/dirk.jl:273
ERROR: Newton did not converge

julia> solve(ode, Theseus.DIRK43(); 
             adaptive = false, dt = 0.25, 
             newton_tol_abs = 1.0e-10, newton_tol_rel = 0.0, 
             newton_max_niter = 1000)
┌ Warning: Newton did not converge
│   stats = (solved = false, stats = Ariadne.Stats(1001, 8, 4.352315131535924e-10), t = 0.000598125)
│   integrator.t = 0.0
└ @ Theseus ~/.julia/dev/Ariadne/libs/Theseus/src/dirk/dirk.jl:273
ERROR: Newton did not converge

Note that the residual norm reported as the last field of the Stats does not change when allowing more Newton iterations. In this case, we just have to solve a linear system with three components.

[ranocha]

In this case, this is fixed by passing stricter tolerances to the Krylov solver:

julia> solve(ode, Theseus.DIRK43(); 
             adaptive = false, dt = 0.25, 
             newton_tol_abs = 1.0e-10, newton_tol_rel = 0.0, 
             krylov_kwargs = (; atol = 1.0e-11, rtol = 1.0e-11))

We should reconsider the default tolerances passed to the Krylov solver depending on the outer Newton tolerances.