Facing problems with Sundials.jl (ARKODE) with discrete callback

humhongekamyaab · September 11, 2024, 10:59am

Hello everyone,

I want to use CashKarp RK explicit algorithm (Numerical Recipe 2nd Edition RK) from Sundials.jl for my particular ODE system and need a discrete callback to do some manipulations after each time step.
https://sundials.readthedocs.io/en/latest/arkode/Butcher_link.html#cash-karp-6-4-5

J.R. Cash and A.H. Karp. A variable order Runge-Kutta method for initial value problems with rapidly varying right-hand sides. ACM Transactions on Mathematical Software, 16(3):201–222, 1990. doi:10.1145/79505.79507.

I am facing a bug when I am using Sundials (ARKODE) with the discrete callback. Following is the minimal code to reproduce the bug (I have modified the example from the Documentation and added a discrete callback).

using DifferentialEquations
using Sundials
function lorenz(du,u,p,t)
 du[1] = 10.0(u[2]-u[1])
 du[2] = u[1]*(28.0-u[3]) - u[2]
 du[3] = u[1]*u[2] - (8/3)*u[3]
end
u0 = [1.0;0.0;0.0]
tspan = (0.0,100.0)
prob = ODEProblem(lorenz,u0,tspan)

condition(u, t, integrator) = true
function affect!(integrator)
    @info "------------------------------------------------------------------"
end
cb = DiscreteCallback(condition, affect!)

alg = ARKODE(Sundials.Explicit(),etable = Sundials.CASH_KARP_6_4_5)
sol = solve(prob,alg;callback=cb)
using Plots; plot(sol,vars=(1,2,3))

The above minimal code gives the following error:

ERROR: MethodError: no method matching ARKStepReInit(::Sundials.Handle{…}, ::Nothing, ::ODEFunction{…}, ::Float64, ::Vector{…})

Closest candidates are:
  ARKStepReInit(::Any, ::Ptr{Nothing}, ::Ptr{Nothing}, ::Any, ::Any)
   @ Sundials ~/.julia/packages/Sundials/jrHAF/lib/libsundials_api.jl:42
  ARKStepReInit(::Any, ::Ptr{Nothing}, ::Ptr{Nothing}, ::Float64, ::Union{Ptr{Sundials._generic_N_Vector}, Sundials.NVector})
   @ Sundials ~/.julia/packages/Sundials/jrHAF/lib/libsundials_api.jl:35

I didn’t find the direct implementation of this algorithm in DifferentialEquations.jl and that’s why I want to use the Sundials interface.

Also, it seems that there is no unit test for ARKODE with callback in Sundials.jl unit tests …

github.com

SciML/Sundials.jl/blob/master/test/common_interface/callbacks.jl

using Sundials, Test

callback_f = function (du, u, p, t)
    du[1] = u[2]
    du[2] = -9.81
end

condtion = function (u, t, integrator) # Event when event_f(u,t,k) == 0
    u[1]
end

affect! = nothing
affect_neg! = function (integrator)
    integrator.u[2] = -integrator.u[2]
end

callback = ContinuousCallback(condtion, affect!, affect_neg!)

u0 = [50.0, 0.0]
tspan = (0.0, 15.0)

This file has been truncated. show original

Could anyone help me to fix the above bug? I am not sure if the above interface to call Sundials ARKODE is correct, maybe I am missing something …

Thanks

ChrisRackauckas · September 12, 2024, 1:19am

Is there a specific reason behind this choice? We have not found a single problem for which ARKODE in non-IMEX mode for something that is not large PDE actually benchmarks well, so we generally never recommend actually using ARKODE over the native Julia methods except for academic benchmarks of ODE solvers. Specifically for Cash Karp, Tsit5 or BS5 are generally much more optimized for the same scenario so I’m surprised this one is coming up.

Note that the Sundials ARK formulation of the CashKarp method is not the adaptive order one of the paper, but instead just a tableau form of the 4th and 5th order methods from the paper. So in that context it’s simply the tableau and not the extra bits.

humhongekamyaab · September 12, 2024, 5:13am

Thanks for the quick reply.

I have a very large ode system of around 10000 dofs, and I have already tried DP and Tsi5 algorithms but they didnt work.

My adhoc RK CashKarp (numerical recipe 2nd) implementation magically works. I already knew that CashKarp alg work for these systems from my colleagues. They know it from their experience rather than any theoretical reasons.

I opted for Sundials as Julia doesnt have CashKarp implemented, otherwise I need to write my own CashKarp within DiffEq package framework.

ChrisRackauckas · September 12, 2024, 9:08am

Is that ad-hoc one with error control? My guess is that with the large system, it might be stiff and the adaptivity could be turning the dt down?

There is a CashKarp in the ExplicitRK form. ExplicitRK(tableau=constructCashKarp()) (requires using DiffEqDevTools).

github.com

SciML/DiffEqDevTools.jl/blob/v2.45.0/src/ode_tableaus.jl#L1227-L1248


      
          """
          constructCashKarp()
          
          Constructs the tableau object for the Cash-Karp Order 4/5 method.
          """
          function constructCashKarp(T::Type = Float64)
              A = [0 0 0 0 0 0
                   1//5 0 0 0 0 0
                   3//40 9//40 0 0 0 0
                   3//10 -9//10 6//5 0 0 0
                   -11//54 5//2 -70//27 35//27 0 0
                   1631//55296 175//512 575//13824 44275//110592 253//4096 0]
              c = [0; 1 // 5; 3 // 10; 3 // 5; 1; 7 // 8]
              α = [37 // 378; 0; 250 // 621; 125 // 594; 0; 512 // 1771]
              αEEst = [2825 // 27648; 0; 18575 // 48384; 13525 // 55296; 277 // 14336; 1 // 4]
              A = map(T, A)
              α = map(T, α)
              c = map(T, c)
              αEEst = map(T, αEEst)
              return (ExplicitRKTableau(A, c, α, 5, αEEst = αEEst, adaptiveorder = 4,

This file has been truncated. show original

For examples of ExplicitRK in use with chosen tableaus, you can see for example:

github.com

SciML/OrdinaryDiffEq.jl/blob/v6.89.0/test/interface/ode_initdt_tests.jl#L1-L10


      
          using OrdinaryDiffEq, DiffEqDevTools, Test
          import ODEProblemLibrary: prob_ode_linear, prob_ode_2Dlinear
          
          prob = prob_ode_linear
          sol = solve(prob, Rosenbrock32())
          dt₀ = sol.t[2]
          
          prob = prob_ode_2Dlinear
          sol = solve(prob, ExplicitRK(tableau = constructBogakiShampine3()))
          dt₀ = sol.t[2]

Though I can almost guarantee you that it’s not a tableau thing and is likely more to do with adaptivity trying to use error controls and your adhoc version not being error corrected. But you can give this a try and see.

humhongekamyaab · September 13, 2024, 9:11am

Hi @ChrisRackauckas,
Thanks for pointing it out. It was difficult to find in the documentation.

CashKarp ExplicitRK(tableau = constructCashKarp()) magically works for me. I don’t know the reason.

I also tried the following algorithms keeping the same error tolerance

ExplicitRK(tableau = constructDormandPrince())
DP5()
Tsit5()

but all of them fail at around the same time step. Failure occurs due to u becoming out of the domain (I have a \log function in du = f(u, p, t)) during intermediate stages of RK. I also tried specifying ifoutofdomain function in solver options (Common Solver Options (Solve Keyword Arguments) · DifferentialEquations.jl) to reject such time steps, but it only works for time step not intermediate RK stages. I don’t know if such an option exists, but it would be great to have such a function.

ChrisRackauckas · September 14, 2024, 6:24am

I’d love to get a hand on your case there because there’s nothing in the theory that would suggest why this would be the case. I’d want to really dig in and see.

As a try, did you also try VCABM?