Erratic convergence behavior of DIRK and IMEX methods

[ranocha]

I get the following results:

julia> using Theseus, LinearAlgebra

julia> function compute_errors(prob, u_ana, alg, dts; kwargs...)
           errors = similar(dts)
           for (i, dt) in enumerate(dts)
               sol = solve(prob, alg;
                           dt, adaptive = false, kwargs...)
               u_num = sol.u[end]
               errors[i] = norm(u_num - u_ana)
           end
           return errors
       end

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> u_ana = [cos(ode.tspan[end]),
                sin(ode.tspan[end]),
                exp(-ode.tspan[end])]

julia> compute_errors(ode, u_ana, Theseus.DIRK43(), 2.0 .^ (-2:-1:-18);
                      krylov_tol_abs = 1.0e-14)
┌ Warning: TODO forward zero-set of memorycopy used memset rather than runtime type 
│ Caused by:
│ Stacktrace:
│   [1] copy
│     @ ./array.jl:350
│   [2] unaliascopy
│     @ ./abstractarray.jl:1516
│   [3] unalias
│     @ ./abstractarray.jl:1500
│   [4] broadcast_unalias
│     @ ./broadcast.jl:946
│   [5] preprocess
│     @ ./broadcast.jl:953
│   [6] preprocess_args
│     @ ./broadcast.jl:956
│   [7] preprocess_args
│     @ ./broadcast.jl:955
│   [8] preprocess
│     @ ./broadcast.jl:952
│   [9] preprocess_args
│     @ ./broadcast.jl:956
│  [10] preprocess_args
│     @ ./broadcast.jl:955
│  [11] preprocess
│     @ ./broadcast.jl:952
│  [12] override_bc_copyto!
│     @ ~/.julia/packages/Enzyme/OhLFb/src/compiler/interpreter.jl:818
│  [13] copyto!
│     @ ./broadcast.jl:925
│  [14] materialize!
│     @ ./broadcast.jl:883
│  [15] materialize!
│     @ ./broadcast.jl:880
│  [16] DIRK
│     @ ~/.julia/dev/Ariadne/libs/Theseus/src/dirk/dirk.jl:33
└ @ Enzyme.Compiler ~/.julia/packages/Enzyme/OhLFb/src/rules/llvmrules.jl:812
17-element Vector{Float64}:
 0.0009095839842415547
 6.902787265585153e-5
 4.536884404070524e-6
 2.312461892802833e-7
 2.1474333974267422e-8
 6.248482209851633e-8
 7.765078477201435e-9
 7.009412413653584e-9
 1.453399550413022e-7
 3.635034002432725e-8
 1.2854643939353188e-8
 3.213969651607027e-9
 8.035280980953472e-10
 2.008842509197221e-10
 5.020962834785772e-11
 1.2559764774477285e-11
 3.1235541846379986e-12

Even after ignoring the annoying warning (see #69), it's strange that the error remains around 1.0e-8 for quite a long time before decreasing again.

[ranocha]

This also happens for IMEX methods: After executing the code from the first post, I get

julia> function rhs_split!(du, u, p, t)
           du[1] = -u[2] / 2
           du[2] = u[1] / 2
           du[3] = -u[3] / 2
           return nothing
       end

julia> ode_split = SplitODEProblem(rhs_split!, rhs_split!,
                                   ode.u0, ode.tspan)

julia> compute_errors(ode_split, u_ana, Theseus.ARS443(), 2.0 .^ (-2:-1:-18);
                      krylov_tol_abs = 1.0e-14)
17-element Vector{Float64}:
 0.00019949612870946836
 2.600103439413053e-5
 3.2022976187164364e-6
 1.036933672166747e-6
 1.298766785151961e-7
 1.6250467769254345e-8
 1.3629445498641903e-6
 3.407871188935382e-7
 8.520315094817657e-8
 2.130158232832665e-8
 5.325496500480966e-9
 1.331388184274275e-9
 3.328534156348306e-10
 8.322105499780579e-11
 2.080304660347559e-11
 5.195664399238211e-12
 1.3634273316631619e-6

(together with warnings similar to the one above).