SciML / Sundials.jl

Julia interface to Sundials, including a nonlinear solver (KINSOL), ODE's (CVODE and ARKODE), and DAE's (IDA) in a SciML scientific machine learning enabled manner
https://diffeq.sciml.ai
BSD 2-Clause "Simplified" License
208 stars 78 forks source link

IDA-key error when saving positions #290

Open JKRT opened 3 years ago

JKRT commented 3 years ago

Hi!

I have the following system (see below) describing a bouncing ball. However, when running this with julia> plot(BouncingBallRealsSimulate((0.0, 3.0)), vars=(2))

I get the the following error:

┌ Warning: IDAGetDky failed with error code =
│   flag = -25
└ @ Sundials C:\Users\John\.julia\packages\Sundials\xAoTr\src\simple.jl:20

The plot however, seems to be as I would expect

image

I found one post by Martin Otter, in the Sundials package. However, I am not sure..

using DiffEqBase
using DifferentialEquations
using Plots
using Sundials
using Revise

function BouncingBallRealsStartConditions(p, t0)
  local x0 = Array{Float64}(undef, 2)
  local dx0 = Array{Float64}(undef, 2)
  x0[2] = (1.0) #= h =#
  return x0, dx0
end

function BouncingBallRealsDifferentialVars()
  return Bool[1, 1]
end

function BouncingBallRealsDAE_equations(res, dx, x, p, t) #=time=#
  res[1] = ((dx[1]) - ((-(p[1])))) #= g =#
  res[2] = ((dx[2]) - (x[1])) #= v =#
end

function BouncingBallRealsParameterVars()
  p = Array{Float64}(undef, 2)
  p[1] = (9.81) #= g =#
  p[2] = (0.7) #= e =#
  return p
end

function BouncingBallRealsCallbackSet(p)#=parameters=#
  condition1(x, t, integrator) = begin
    test = x[2] - 0
    return test
  end
  affect1!(integrator) =
    begin
      integrator.u[1] = -((p[2]) * ((integrator.u[1]))) #= e =#
    end
  cb1 = ContinuousCallback(condition1,
                           affect1!,
                           rootfind=true,
                           save_positions=(true, true),
                           affect_neg! = affect1!)
  return CallbackSet(cb1)
end

function BouncingBallRealsSimulate(tspan = (0.0, 5.0))
  # Define problem
  p = BouncingBallRealsParameterVars()
  (x0, dx0) = BouncingBallRealsStartConditions(p, tspan[1])
  differential_vars = BouncingBallRealsDifferentialVars()
  #= Pass the residual equations =#
  problem = DAEProblem(
    BouncingBallRealsDAE_equations,
    dx0,
    x0,
    tspan,
    p,
    differential_vars = differential_vars,
    callback = BouncingBallRealsCallbackSet(p),
  )
  solution = solve(problem, IDA())
  return solution
end

Cheers, John

ChrisRackauckas commented 3 years ago

It means that the interpolation did something funny, and we can investigate that, but the integration seems fine. Transferring to Sundials.jl

JKRT commented 3 years ago

Forgot to update:

using DiffEqBase
using DifferentialEquations
using Plots
using Sundials
using Revise

function BouncingBallRealsCallbackSet(p)#=parameters=#
  condition1(x, t, integrator) = begin
    test = x[2] - 0
    return test
  end
  affect1!(integrator) =
    begin
      integrator.u[1] = -((p[2]) * ((integrator.u[1]))) #= e =#
    end
  cb1 = ContinuousCallback(condition1,
                           affect1!,
                           rootfind=true,
                           save_positions=(false, false),
                           affect_neg! = affect1!)
  return CallbackSet(cb1)
end

Setting save_positions to (false, false) does not give the warning, plot still does look alright.

using DiffEqBase
using DifferentialEquations
using Plots
using Sundials
using Revise

function BouncingBallRealsStartConditions(p, t0)
  local x0 = Array{Float64}(undef, 2)
  local dx0 = Array{Float64}(undef, 2)
  x0[2] = (1.0) #= h =#
  return x0, dx0
end

function BouncingBallRealsDifferentialVars()
  return Bool[1, 1]
end

function BouncingBallRealsDAE_equations(res, dx, x, p, t) #=time=#
  res[1] = ((dx[1]) - ((-(p[1])))) #= g =#
  res[2] = ((dx[2]) - (x[1])) #= v =#
end

function BouncingBallRealsParameterVars()
  p = Array{Float64}(undef, 2)
  p[1] = (9.81) #= g =#
  p[2] = (0.7) #= e =#
  return p
end

function BouncingBallRealsCallbackSet(p)#=parameters=#
  condition1(x, t, integrator) = begin
    test = x[2] - 0
    return test
  end
  affect1!(integrator) =
    begin
      integrator.u[1] = -((p[2]) * ((integrator.u[1]))) #= e =#
    end
  cb1 = ContinuousCallback(condition1,
                           affect1!,
                           rootfind=true,
                           save_positions=(true, true),
                           affect_neg! = affect1!)
  return CallbackSet(cb1)
end

function BouncingBallRealsSimulate(tspan = (0.0, 5.0))
  # Define problem
  p = BouncingBallRealsParameterVars()
  (x0, dx0) = BouncingBallRealsStartConditions(p, tspan[1])
  differential_vars = BouncingBallRealsDifferentialVars()
  #= Pass the residual equations =#
  problem = DAEProblem(
    BouncingBallRealsDAE_equations,
    dx0,
    x0,
    tspan,
    p,
    differential_vars = differential_vars,
    callback = BouncingBallRealsCallbackSet(p),
  )
  solution = solve(problem, IDA())
  return solution
end

Looks alright image

However, if I add a saving callback:

    begin
        saved_values_BouncingBallReals = SavedValues(Float64, Tuple{Float64,Array})
        function BouncingBallRealsCallbackSet(aux)
            local p = aux[1]
            begin
                function condition1(x, t, integrator)
                    x[1] - 0.0
                end
                function affect1!(integrator)
                    integrator.u[2] = -(p[1] * integrator.u[2])
                end
                cb1 = ContinuousCallback(
                    condition1,
                    affect1!,
                    rootfind = true,
                    save_positions = (false, false),
                    affect_neg! = affect1!,
                )
            end
            begin
                savingFunction(u, t, integrator) =
                    let
                        (t, deepcopy(integrator.p[2]))
                    end
                cb2 = SavingCallback(savingFunction, saved_values_BouncingBallReals)
            end
            return CallbackSet(cb1, cb2)
        end
    end

The error occurs again. I think it is related to this problem: https://github.com/SciML/Sundials.jl/issues/292 but also https://github.com/SciML/DifferentialEquations.jl/issues/547 it seems Martin Otter got a similar issue with saving callbacks and Sundials

ChrisRackauckas commented 3 years ago

Yes, SavingCallback will not be able to solve this. We need something different in saving if we want to merge that in.

JKRT commented 3 years ago

Ah sorry my mistake, I was imprecise as usual 💯

I did not indent for the saving callback to solve the issue, rather I wanted to highlight that the issue seems to be with saving in general since the error only occurs for me when I:

A) Either run save positions (Which I should) B) Use a saving callback