Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Many events does not work with the saving callback enabled #34

Open
JKRT opened this issue Aug 5, 2021 · 0 comments
Open

Many events does not work with the saving callback enabled #34

JKRT opened this issue Aug 5, 2021 · 0 comments

Comments

@JKRT
Copy link
Owner

JKRT commented Aug 5, 2021

Attempting to handly many events using force save to convert it into an OMSolution fails

using DiffEqBase
using DifferentialEquations
using Plots
using Sundials
import OMBackend
Mode = [0]
function ManyEvents5StartConditions(aux, t)
    local x = zeros(5)
    local dx = zeros(5)
    local p = aux[1]
    local reals = aux[2]
    begin
        begin
            reals[6] = 0.0
            reals[7] = 0.0
            reals[8] = 0.0
            reals[9] = 0.0
            reals[10] = 0.0
        end
        begin
            reals[1] = 0.0
            reals[2] = 0.0
            reals[3] = 0.0
            reals[4] = 0.0
            reals[5] = 0.0
        end
        dx[1] = float(p[2]) / float(p[1])
        dx[2] = float(p[2]) / float(p[1] - 1)
        dx[3] = float(p[2]) / float(p[1] - 2)
        dx[4] = float(p[2]) / float(p[1] - 3)
        dx[5] = float(p[2]) / float(p[1] - 4)
    end
    x[1] = reals[1]
    x[2] = reals[2]
    x[3] = reals[3]
    x[4] = reals[4]
    x[5] = reals[5]
    return (x, dx)
end
function ManyEvents5AuxVarsHandler(res, dx, x, aux, t)
    local p = aux[1]
    local reals = aux[2]
end
function ManyEvents5DifferentialVars()
    return Bool[1, 1, 1, 1, 1]
end
function ManyEvents5DAE_equations(res, dx, x, aux, t)
    ManyEvents5AuxVarsHandler(res, dx, x, aux, t)
    local p = aux[1]
    local reals = aux[2]
    @info "Value of reals: $(reals)"
    res[1] = dx[1] - float(p[2]) / float((p[1] + 1) - 1)
    res[2] = dx[2] - float(p[2]) / float((p[1] + 1) - 2)
    res[3] = dx[3] - float(p[2]) / float((p[1] + 1) - 3)
    res[4] = dx[4] - float(p[2]) / float((p[1] + 1) - 4)
    res[5] = dx[5] - float(p[2]) / float((p[1] + 1) - 5)
    reals[1] = x[1]
    reals[2] = x[2]
    reals[3] = x[3]
    reals[4] = x[4]
    reals[5] = x[5]
end
function ManyEvents5ParameterVars()
    local aux = Array{Array{Float64}}(undef, 2)
    local p = Array{Float64}(undef, 2)
    local reals = Array{Float64}(undef, 10)
    aux[1] = p
    aux[2] = reals
    p[1] = 5
    p[2] = 5
    return aux
end
begin
    saved_values_ManyEvents5 = SavedValues(Float64, Tuple{Float64,Array})
    function ManyEvents5CallbackSet(aux)
        local p = aux[1]
        local reals = aux[2]
        begin
            function condition1(x, t, integrator)
                x[1] - 1.0
            end
            function affect1!(integrator)
                @info "Affect 1"
                reals[6] = 1.0
            end
            cb1 = ContinuousCallback(
                condition1,
                affect1!,
                rootfind = true,
                save_positions = (false, false),
                affect_neg! = affect1!,
            )
        end
        begin
            function condition2(x, t, integrator)
                x[2] - 1.0
            end
            function affect2!(integrator)
                 @info "Affect 2"
                reals[7] = 1.0
            end
            cb2 = ContinuousCallback(
                condition2,
                affect2!,
                rootfind = true,
                save_positions = (false, false),
                affect_neg! = affect2!,
            )
        end
        begin
            function condition3(x, t, integrator)
                x[3] - 1.0
            end
            function affect3!(integrator)
                @info "Affect 3"
                reals[8] = 1.0
            end
            cb3 = ContinuousCallback(
                condition3,
                affect3!,
                rootfind = true,
                save_positions = (false, false),
                affect_neg! = affect3!,
            )
        end
        begin
            function condition4(x, t, integrator)
                x[4] - 1.0
            end
            function affect4!(integrator)
                @info "Affect 4"
                reals[9] = 1.0
            end
            cb4 = ContinuousCallback(
                condition4,
                affect4!,
                rootfind = true,
                save_positions = (false, false),
                affect_neg! = affect4!,
            )
        end
        begin
            function condition5(x, t, integrator)
                x[5] - 1.0
            end
            function affect5!(integrator)
                @info "Affect 5"
                reals[10] = 1.0
            end
            cb5 = ContinuousCallback(
                condition5,
                affect5!,
                rootfind = true,
                save_positions = (false, false),
                affect_neg! = affect5!,
            )
        end            
        return CallbackSet(cb1, cb2, cb3, cb4, cb5)
    end
end
function ManyEvents5Simulate(tspan = (0.0, 1.0))
    local aux = ManyEvents5ParameterVars()
    (x0, dx0) = ManyEvents5StartConditions(aux, tspan[1])
    local differential_vars = ManyEvents5DifferentialVars()
    local problem = DAEProblem(
        ManyEvents5DAE_equations,
        dx0,
        x0,
        tspan,
        aux,
        differential_vars = differential_vars,
        callback = ManyEvents5CallbackSet(aux),
    )
    local solution = solve(problem, IDA())
    local savedSol = map(collect, saved_values_ManyEvents5.saveval)
    local t = [(savedSol[i])[1] for i = 1:length(savedSol)]
    local vars = [(savedSol[i])[2] for i = 1:length(savedSol)]
    local T = eltype(eltype(vars))
    local N = length(aux[2])
    local nsolution = DAESolution{
        Float64,
        N,
        typeof(vars),
        Nothing,
        Nothing,
        Nothing,
        typeof(t),
        typeof(problem),
        typeof(solution.alg),
        typeof(solution.interp),
        typeof(solution.destats),
    }(
        vars,
        nothing,
        nothing,
        nothing,
        t,
        problem,
        solution.alg,
        solution.interp,
        solution.dense,
        0,
        solution.destats,
        solution.retcode,
    )
    ht = Dict{Any,Any}(
        5 => "x₅",
        4 => "x₄",
        6 => "e₁",
        7 => "e₂",
        2 => "x₂",
        10 => "e₅",
        9 => "e₄",
        8 => "e₃",
        3 => "x₃",
        1 => "x₁",
    )
    omSolution = OMBackend.Runtime.OMSolution(nsolution, ht)
    return omSolution
end

Removing the saving callback seems to result in the correct behavior

using DiffEqBase
using DifferentialEquations
using Plots
using Sundials
import OMBackend
Mode = [0]
function ManyEvents5StartConditions(aux, t)
    local x = zeros(5)
    local dx = zeros(5)
    local p = aux[1]
    local reals = aux[2]
    begin
        begin
            reals[6] = 0.0
            reals[7] = 0.0
            reals[8] = 0.0
            reals[9] = 0.0
            reals[10] = 0.0
        end
        begin
            reals[1] = 0.0
            reals[2] = 0.0
            reals[3] = 0.0
            reals[4] = 0.0
            reals[5] = 0.0
        end
        dx[1] = float(p[2]) / float(p[1])
        dx[2] = float(p[2]) / float(p[1] - 1)
        dx[3] = float(p[2]) / float(p[1] - 2)
        dx[4] = float(p[2]) / float(p[1] - 3)
        dx[5] = float(p[2]) / float(p[1] - 4)
    end
    x[1] = reals[1]
    x[2] = reals[2]
    x[3] = reals[3]
    x[4] = reals[4]
    x[5] = reals[5]
    return (x, dx)
end
function ManyEvents5AuxVarsHandler(res, dx, x, aux, t)
    local p = aux[1]
    local reals = aux[2]
end
function ManyEvents5DifferentialVars()
    return Bool[1, 1, 1, 1, 1]
end
function ManyEvents5DAE_equations(res, dx, x, aux, t)
    ManyEvents5AuxVarsHandler(res, dx, x, aux, t)
    local p = aux[1]
    local reals = aux[2]
    @info "Value of reals: $(reals)"
    res[1] = dx[1] - float(p[2]) / float((p[1] + 1) - 1)
    res[2] = dx[2] - float(p[2]) / float((p[1] + 1) - 2)
    res[3] = dx[3] - float(p[2]) / float((p[1] + 1) - 3)
    res[4] = dx[4] - float(p[2]) / float((p[1] + 1) - 4)
    res[5] = dx[5] - float(p[2]) / float((p[1] + 1) - 5)
    reals[1] = x[1]
    reals[2] = x[2]
    reals[3] = x[3]
    reals[4] = x[4]
    reals[5] = x[5]
end
function ManyEvents5ParameterVars()
    local aux = Array{Array{Float64}}(undef, 2)
    local p = Array{Float64}(undef, 2)
    local reals = Array{Float64}(undef, 10)
    aux[1] = p
    aux[2] = reals
    p[1] = 5
    p[2] = 5
    return aux
end
begin
    saved_values_ManyEvents5 = SavedValues(Float64, Tuple{Float64,Array})
    function ManyEvents5CallbackSet(aux)
        local p = aux[1]
        local reals = aux[2]
        begin
            function condition1(x, t, integrator)
                x[1] - 1.0
            end
            function affect1!(integrator)
                @info "Affect 1"
                reals[6] = 1.0
            end
            cb1 = ContinuousCallback(
                condition1,
                affect1!,
                rootfind = true,
                save_positions = (false, false),
                affect_neg! = affect1!,
            )
        end
        begin
            function condition2(x, t, integrator)
                x[2] - 1.0
            end
            function affect2!(integrator)
                 @info "Affect 2"
                reals[7] = 1.0
            end
            cb2 = ContinuousCallback(
                condition2,
                affect2!,
                rootfind = true,
                save_positions = (false, false),
                affect_neg! = affect2!,
            )
        end
        begin
            function condition3(x, t, integrator)
                x[3] - 1.0
            end
            function affect3!(integrator)
                @info "Affect 3"
                reals[8] = 1.0
            end
            cb3 = ContinuousCallback(
                condition3,
                affect3!,
                rootfind = true,
                save_positions = (false, false),
                affect_neg! = affect3!,
            )
        end
        begin
            function condition4(x, t, integrator)
                x[4] - 1.0
            end
            function affect4!(integrator)
                @info "Affect 4"
                reals[9] = 1.0
            end
            cb4 = ContinuousCallback(
                condition4,
                affect4!,
                rootfind = true,
                save_positions = (false, false),
                affect_neg! = affect4!,
            )
        end
        begin
            function condition5(x, t, integrator)
                x[5] - 1.0
            end
            function affect5!(integrator)
                @info "Affect 5"
                reals[10] = 1.0
            end
            cb5 = ContinuousCallback(
                condition5,
                affect5!,
                rootfind = true,
                save_positions = (false, false),
                affect_neg! = affect5!,
            )
        end            
        return CallbackSet(cb1, cb2, cb3, cb4, cb5)
    end
end
function ManyEvents5Simulate(tspan = (0.0, 1.0))
    local aux = ManyEvents5ParameterVars()
    (x0, dx0) = ManyEvents5StartConditions(aux, tspan[1])
    local differential_vars = ManyEvents5DifferentialVars()
    local problem = DAEProblem(
        ManyEvents5DAE_equations,
        dx0,
        x0,
        tspan,
        aux,
        differential_vars = differential_vars,
        callback = ManyEvents5CallbackSet(aux),
    )
    solution = solve(problem, IDA())
    return solution
end

image

Probably related to #33

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

1 participant