SavingCallback for DAEs: Not clear how to get du

[MartinOtter]

The callback save_func(u, t, integrator) can be provided to compute additional outputs y based on the current time t and the current state u (and the model parameters integrator.p): y=f(u,t,p).

To compute additional outputs y for DAEs, also the derivatives of the states du are in general needed (y = f(du, u, t, p)), but du is not in the interface of save_func. Is it possible to get du somehow from the integrator (e.g. integrator.du), or is a function call needed?

[ChrisRackauckas]

integrator(t,Val{1}) probably does it (Val{1} means first derivative from the interpolant).

[MartinOtter]

Thanks. In principal this works. However, there are the following issues:

Fails for first time point

At the first time point

du = integrator(t,Val{1})

triggers the following Sundials error:

[IDAS ERROR]  IDAGetDky
  Illegal value for k.

┌ Warning: IDAGetDky failed with error code =
│   flag = -25
└ @ Sundials D:\otter\home\.julia\dev\Sundials\src\simple.jl:20

I guess, this can be seen as a bug in IDA, because IDA has du at the initial point internally stored and could just return it.

Temporarily, I use the code:

function save_func(u, t, integrator)
    p = integrator.p
    if t==0
        du = integrator.du
    else
        du = integrator(t, Val{1})   # allocates memory for der_x
    end
    ...
end

Would be better if this bug-fix would be included in function integrator(...).

I tried whether this problem is present with DASKR, but unfortunately DASKR does not support callbacks.

Allocates memory at every call

du = integrator(t,Val{1})

always allocates (unnecessary) memory for vector du. It would be better to use the call:

integrator(du, t, Val{1})

However, then "somehow" a cache-vector du needs to be provided that is specific for the used integrator (e.g. IDA needs an N_Vector array). Temporarily, I use the code:

function save_func(u, t, integrator)
    p = integrator.p
    if t==0
        p.du = copy(integrator.du)
    else
        integrator(p.du, t, Val{1})
    end
    du = p.du
    ...
end

This is also not nice (e.g. requires to define an Any vector du in p)

It would be better, if a second tmp vector of length(u) could be stored in integrator and this tmp vector could be utilized at this place (in a similar way as for u)

[ChrisRackauckas]

I guess, this can be seen as a bug in IDA, because IDA has du at the initial point internally stored and could just return it.

Yeah 🤷‍♂

However, then "somehow" a cache-vector du needs to be provided that is specific for the used integrator (e.g. IDA needs an N_Vector array). Temporarily, I use the code:

https://docs.juliadiffeq.org/latest/basics/integrator/#Caches-1

first(get_tmp_cache(integrator)) is likely what you want.

[MartinOtter]

Unfortunately, this does not work for IDA:

When printing the size of the cache vector:

function save_func(u, t, integrator)
    tmps = get_tmp_cache(integrator)
    println("typeof(tmps) = ", typeof(tmps), ", length(tmps) = ", length(tmps))
    ...
end
...
typeof(tmps) = Tuple{Array{Float64,1}}, length(tmps) = 1

So, there is only one cache vector and this cache vector is already used in function (affect!::SavingAffect) in file DiffEqCallbacks\src\saving.jl:

            # it should really be a mutability trait
            if typeof(integrator.u) <: Union{Number,SArray, ArrayPartition{T, Tuple{SArray{R, T, N, M}, SArray{R, T, N, M}}} where {T, R, N, M}}
                curu = integrator(curt)
            else
                curu = first(get_tmp_cache(integrator))
                integrator(curu,curt) # inplace since save_func allocates
            end
            copyat_or_push!(affect!.saved_values.t, affect!.saveiter, curt)
            copyat_or_push!(affect!.saved_values.saveval, affect!.saveiter,
                            affect!.save_func(curu, curt, integrator),Val{false})

curu = first(get_tmp_cache(integrator)) gets the cache vector and uses it in affect!.safe_func(curu, ...).

Could you add more cache vectors to IDA.

[ChrisRackauckas]

I think we need to dig into Sundials and see if we can grab more cache vectors that it's already made. That's probably better than adding another register, since there's gotta be one in there.

[MartinOtter]

function save_func(u, t, integrator)
   p = integrator.p
   if t==0
       p.du = copy(integrator.du)
   else
       integrator(p.du, t, Val{1})
   end
   du = p.du
   ...
end

This is also not nice (e.g. requires to define an Any vector du in p)

I tried to use a Float64 vector instead of an Any vector:

mutable struct Model
   ...
   du::Vector{Float64}
end
p = Model(...)

Surprisingly this worked. Since IDA is using an N_Vector array, this could mean that at every model evaluation there is some overhead associated to copy the N_Vector to a Julia vector.

[ChrisRackauckas]

It should just be doing an unsafe_wrap, but FFI is not my specialty so it would be good if someone double check this:

https://github.com/JuliaDiffEq/Sundials.jl/blob/master/src/common_interface/function_types.jl#L22-L23