IDASetUserData sets user_data as Ptr{Nothing}

[bertapedret]

Hi!

I want to use Sundials IDAS for Forward Sensitivity Analysis on a simple DAE problem. I have built 1st the DAE solver alone by following the steps on the documentation and it works fine without the parameters. To add Sensitivity Analysis I require the user_data pointer that is passed to all other functions.

I have tried passing this data with the function IDASetUserData(ida_mem, user_data) for many types of user_data array but all of them appear as Ptr{Nothing} inside the residual function (I'm currently printing user_data there as I can't find any function to retrieve the value of user_data).

The print can be seen either at IDACalcIC or IDASolve. This is the code I'm using:

# Goal - Complete the Sensitivity Analysis on the system of DAEs: isomerization kinetics A -> B <-> C
## ************** DAE model ******************
# dCA/dt = - k1*CA
# dCB/dt = k1*CA - k2*CB + k3*CC
# CA(t) + CB(t) + CC(t) = CA(0) + CB(0) + CC(0)

using Sundials

CA_0, CB_0, CC_0 = 1.0, 0.0, 0.0; # concentrations [mol/L]
k1, k2, k3 = 1000.0, 1.0, 1.0;    # rate constants [1/s]
 
# The function that computes the residual F in the DAE. Has the form res(t, yy, yp, resval, user_data)
function F_idas(t, yy_nv, yp_nv, rr_nv, user_data)
    println(user_data)
    yy = convert(Vector, yy_nv)
    yp = convert(Vector, yp_nv)
    rr = convert(Vector, rr_nv)
    CA_0, k1, k2, k3 = θ; #convert(Vector, user_data); # HOW TO TRANSFER/LOAD USER DATA???
    CA, CB, CC = yy
    dCA, dCB, dCC = yp
    rr[1] = -k1*CA - dCA
    rr[2] = k1*CA - k2*CB + k3*CC - dCB
    rr[3] = CA + CB + CC - (CA_0 + CB_0 + CC_0)
    return Sundials.IDA_SUCCESS
end

# 3. Set vectors of initial values:
t0 = 0.0
tout1 = 0.02
yy0 = [CA_0, CB_0, CC_0];
yp0 = [-10.0, 2.0, 2.0];
rtol = 1e-4
avtol = [1e-8, 1e-8, 1e-6]

θ = [CA_0; k1; k2; k3]; # Set of parameters and values - to do the SA:

# 4. Create idas object
ida_mem = Sundials.IDACreate();

# HERE I TRY TO CREATE THE USER_DATA POINTER TO THE ACTUAL PARAMETERS THAT I WANT TO PASS ************
user_data = convert(Sundials.N_Vector, θ) # create NVector 
Sundials.@checkflag Sundials.IDASetUserData(ida_mem, user_data) # Set up IDAS user data

Sundials.@checkflag Sundials.IDAInit(ida_mem, F_idas, t0, yy0, yp0) # 5. Initialize idas solver 
Sundials.@checkflag Sundials.IDASVtolerances(ida_mem, rtol, avtol) # 6. Specify integration tolerances
M = Sundials.SUNDenseMatrix(length(yy0), length(yy0)) # 7. Create matrix object
LS = Sundials.SUNLinSol_Dense(yy0, M) # 8. Create linear solver
Sundials.@checkflag Sundials.IDADlsSetLinearSolver(ida_mem, LS, M) # 10. Attach linear solver module

Sundials.@checkflag Sundials.IDACalcIC(ida_mem, Sundials.IDA_Y_INIT, tout1) # 15. Correct initial values

# 17. Advance solution in time
tout = tout1 # the next time at which a computed solution is desired.
tret = [1.0] # the time reached by the solver (output).
yy = similar(yy0) #the computed solution vector y
yp = similar(yp0) #the computed solution vector y'
retval = Sundials.IDASolve(ida_mem, tout, tret, yy, yp, Sundials.IDA_NORMAL)

# ** DEALLOCATE MEMORY:
# 19. Deallocate memory for solution vectors
Sundials.N_VDestroy(yy)
Sundials.N_VDestroy(yp)

# 20. Free solver memory
Sundials.IDAFree(ida_mem)

# 22. Free linear solver and matrix memory
Sundials.SUNLinSolFree_Dense(LS)
Sundials.SUNMatDestroy_Dense(M)

[ChrisRackauckas]

https://github.com/SciML/Sundials.jl/blob/master/test/cvodes_dns.jl is the closest we've gotten to wrapping the sensitivity parts.