Chaostools tangent_integrator stiff solvers

Specific Domains Modelling & Simulations
1

Hi,
I want to use the tangent_integrator function from the ChaosTools package with solvers for stiff problems:

integ = tangent_integrator(cds, w0; alg=Rodas5(), abstol=1e-9, reltol=1e-9)
step!(integ, 1.0)

but this gives me:

MethodError: no method matching Float64(::ForwardDiff.Dual{ForwardDiff.Tag{DiffEqDiffTools.UJacobianWrapper{ODEFunction{true,getfield(DynamicalSystemsBase, Symbol(“##21#22”)){6,typeof(f),typeof(jac)},LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Float64,Array{Float64,1}},Float64},Float64,11})
Closest candidates are:
Float64(::Real, !Matched::RoundingMode) where T<:AbstractFloat at rounding.jl:194
Float64(::T<:Number) where T<:Number at boot.jl:741
Float64(!Matched::Int8) at float.jl:60
…

I had this problem before when calling radau-fortran solver, see https://discourse.julialang.org/t/lyapunov-exponents-using-chaostools-and-fortran-radau-solver/26007.

I also tried, as suggested in the documentation, to chose a solver by alg_hint keyword, but I had a look into the implementation of tangent_integrator and this is not making use of the hint-keyword. How can I use e.g. Rodas5() with the tangent:integratorfunction ?

2

Hey there.

  1. The alg_hint is not used by tangent_integrator, it is DifferentialEquations.jl exclusive, sorry.
  2. To be able to reproduce your problem, it is useful to provide a MinimalWorkingExample, some runnable code that leads to the same error.
  3. Depending on your dynamical system, it may just be the case that Rodas5 does not support the state type that the tangent integrator uses, since this can vary from matrix to vector-of-vectors.
3

@Datseris I will try to make up a minimal example quickly.

What is the “state type”?

Do you mean by state type: in-place or out-of-place versions of the functions?

4

The “state type” is simple typeof(integ.u) , i.e. what is the
Julia structure that represents the state of the integrator. This
is of course affected by whether the integrator is in-place or
not, but there are other factors as well.

5

So I made a working example with the Lorenz system:

function lorenz_iip(du,u,p,t)
    sigma = p[1]
    rho = p[2]
    beta = p[3]

    # f
    du[1] = sigma*(u[2]-u[1])
    du[2] = u[1]*(rho-u[3])-u[2]
    du[3] = u[1]*u[2]-beta*u[3]

end

function lorenz_jac_iip(J,u,p,t)

    sigma = p[1]
    rho = p[2]
    beta = p[3]

    J[1,1] = -sigma
    J[1,2] = sigma
    J[1,3] = 0
    J[2,1] = rho-u[3]
    J[2,2] = -1
    J[2,3] = -u[1]
    J[3,1] = u[2]
    J[3,2] = u[1]
    J[3,3] = -beta
end

Random.seed!(42)
N=10000
tol = 1e-12
w0 = Matrix(LinearAlgebra.qr(Random.rand(3,3)).Q)
param = [16.,45.92,4.]
x_start = [19.,20.,50.]

cds = ContinuousDynamicalSystem(lorenz_iip,x_start,param,lorenz_jac_iip)
integ_rodas = tangent_integrator(cds; alg=Rodas5(), abstol=tol, reltol=tol)
integ = tangent_integrator(cds; abstol=tol, reltol=tol)

and when I call step!(integ) this works but calling step!(integ_rodas) gives me the above error message:

MethodError: no method matching Float64(::ForwardDiff.Dual{ForwardDiff.Tag{DiffEqDiffTools.UJacobianWrapper{ODEFunction{true,getfield(DynamicalSystemsBase, Symbol(“##21#22”)){3,typeof(lorenz_iip),typeof(lorenz_jac_iip)},LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Float64,Array{Float64,1}},Float64},Float64,12})
Closest candidates are:
Float64(::Real, !Matched::RoundingMode) where T<:AbstractFloat at rounding.jl:194
Float64(::T<:Number) where T<:Number at boot.jl:741
Float64(!Matched::Int8) at float.jl:60
…

1 Like
6

I checked the states:
typeof(integ.u) and typeof(integ_rodas.u) both give

Array{Float64,2}

7

Okay, a first question to ask before trying anything else is: @ChrisRackauckas , does Rodas5 work with Matrix states?

8

Yes

9

Are you creating any cache variables? If you are, are they dual cached?

10

I tested some more algorithms, and so far none of Kvaerno5, `KenCarp4, Rodas4P worked.

cds = ContinuousDynamicalSystem(lorenz_iip,x_start,param)
  • leaving out the Jacobian (I assume its auto diffed then), does not work either.
11

I am just running the code above in a Jupyter notebook. What do you mean by cached/dual cached in this context?

12

To the best of my understanding, tangent_integrator propagates everything to __init and ODEProblem. See here: https://github.com/JuliaDynamics/DynamicalSystemsBase.jl/blob/master/src/continuous.jl#L50-L68 :

function tangent_integrator(ds::CDS{IIP}, Q0::AbstractMatrix;
    u0 = ds.u0, diffeq...) where {IIP}

    t0 = ds.t0
    Q = safe_matrix_type(Val{IIP}(), Q0)
    u = safe_state_type(Val{IIP}(), u0)

    k = size(Q)[2]
    k > dimension(ds) && throw(ArgumentError(
    "It is not possible to evolve more tangent vectors than the system's dimension!"
    ))

    tangentf = create_tangent(ds.f, ds.jacobian, ds.J, Val{IIP}(), Val{k}())
    tanprob = ODEProblem{IIP}(tangentf, hcat(u, Q), (t0, typeof(t0)(Inf)), ds.p)

    solver = _get_solver(diffeq)
    return __init(tanprob, solver; DEFAULT_DIFFEQ_KWARGS..., internalnorm = _tannorm,
                  save_everystep = false, diffeq...)
end

(also, since the problem comes from using different integration algorithms, it cannot stem from DynamicalSystems.jl since the step! functionality is beyond our control; we only initialize the states and functions and give everything to init)

13

Can you point to the implementation of create_tangent? How do you tie together the “f” and “jacobian”? Could this be something, the stiff solvers do not like?

14

The tangent space function is created generally in this part of the code: https://github.com/JuliaDynamics/DynamicalSystemsBase.jl/blob/master/src/dynamicalsystem.jl#L342

e.g. for the case of in-place system with provided jacobian function one has (https://github.com/JuliaDynamics/DynamicalSystemsBase.jl/blob/master/src/dynamicalsystem.jl#L343-L355 )

function create_tangent(@nospecialize(f::F), @nospecialize(jacobian::JAC), J::JM,
    ::Val{true}, ::Val{k}) where {F, JAC, JM, k}
    J = deepcopy(J)
    tangentf = (du, u, p, t) -> begin
        uv = @view u[:, 1]
        f(view(du, :, 1), uv, p, t)
        jacobian(J, uv, p, t)
        mul!((@view du[:, 2:(k+1)]), J, (@view u[:, 2:(k+1)]))
        nothing
    end
    return tangentf
end
15

I can not post the whole stack trace here, as it is too long (over 32000 characters). I can post it somewhere else, if someone is interested :slight_smile:

16

@Datseris So if i switch to Rodas5(autodiff=false) this works. So it seems that differentiating tangentf fails for some reason.

1 Like
17

Hm, nice that you found the root of the problem. I’ve cited above how the tangent function is created, but I don’t really have a clue of what Rodas does (or what it is in general), so I can’t really help on why this fails.

Maybe Chris has some ideas?

18

Hmm, I do know way too little to track down this error further. Are the non stiff algorithms not using ForwardDiff ? Because they work flawlessly here. Maybe @ChrisRackauckas can help?

19

Yup, you’re creating cache variables which aren’t dual cached which violates the assumptions of autodiff.

http://docs.juliadiffeq.org/latest/basics/faq.html#I-get-Dual-number-errors-when-I-solve-my-ODE-with-Rosenbrock-or-SDIRK-methods…?-1

1 Like
20

Ah cool, this seems like something that can be fixed!

@jamblejoe if you want to open a PR for this on DynamicalSystemsBase you are welcomed! (I won’t have time to fix it myself)