Zero crossing termination not correct with zero crossing functions crossing closely together

[MartinOtter]

I have a simulation where zero crossing gives a wrong result. Its hard to construct a simple example with this effect, but the reason for this issue can be clearly described. Hope this is sufficient to find the bug:

The red, green and blue cylinders are falling down and hit the grey cylinder at the same time instant (analytically they would exactly hit at the same time instant, but in the simulation, the time instants are slightly different):

The issue occurs with integrators CVODE_BDF, QBDF, Tsit5 using option rootfind=DifferentialEquations.SciMLBase.RightRootFind, so the problem is not related to the integrator but to the root finding strategy. Furthermore, DifferentialEquations.VectorContinuousCallback(..) is used, so all zero-crossing functions are packed into one vector.

Log of the zero crossing:

State event (zero-crossing) at time = 0.18240210174165705 s
    distance(cylinder1,redCylinder) = -4.884580354174032e-17 became <= 0

State event (zero-crossing) at time = 0.18240210174169105 s
    distance(cylinder1,blueCylinder) = -2.837887238426206e-17 became <= 0

State event (zero-crossing) at time = 0.18244699179731702 s                   #  !!!! wrong
    distance(cylinder1,greenCylinder) = -0.00014384381616780407 became <= 0   #  !!!! wrong
    distance(cylinder1,blueCylinder) = 1.7002864505158043e-16 became > 0

redCylinder and blueCylinder have a zero crossing from positive to negative (negative means they penetrate the grey cylinder1). The values at the RightRootFind are in the order of -1e-17. This is fine.

At some further time instant, there is a zero crossing from the greenCylinder from positive to negative and from the blueCylinder from negative to positive. DifferentialEquations returns with the time instant where the blueCylinder has the zero crossing. But this is wrong, since the greenCylinder has also a zero crossing at this time instant, but with the huge value at the RightRootFind in the order of -1e.4. The root finding iteration should continue, until the greenCylinder root function is in the order of -1e-17, resulting in an earlier time instant of the zero crossing (and blueCylinder has no zero crossing here).

Since an elastic response calculation is used that depends on the penetration depth, completely nonphysical behavior occurs.

[ChrisRackauckas]

Yeah, this case just isn't handled. Not as much of a bug as a not handled case. If you have a VectorContinuousCallback and you have multiple events hit simultaneously (or close to simultaneously), it will still return one index. You could manually check the conditions, but that's not great. We need to actually expand the interface here, have affect!(integrator,idx::Int) and require a separate dispatch affect!(integrator,idxs::Vector). But we want to set this up so it's not allocating or anything, so it's a bit tricky to do well. This is planned for the "not too far future" though, but for now yes it's a known unhandled case.

[MartinOtter]

Any progress here?

Not as much of a bug as a not handled case. If you have a VectorContinuousCallback and you have multiple events hit simultaneously (or close to simultaneously), it will still return one index. You could manually check the conditions, but that's not great.

Sorry, I do not understand this. The issue is that DifferentialEquations does not call VectorContinuousCallback at the right time instant (it is not called when the zero-crossing function of the greenCylinder crosses zero), so I see not any workaround to cope with it. The index vector provided by VectorContinuousCallback is not used by Modia (and would not be used by Modelica-related tools), because anyway the whole model needs to be evaluated to figure out the right branch of every if-clause. The index-vector of VectorContinuousCallback is helpful only in very simple cases.

We have now a more realistic application where this issue occurs and I do not know what to do.

[ChrisRackauckas]

Can you share an example which has the issue? I still haven't seen an example to work on.

[ThummeTo]

Dear @MartinOtter and @ChrisRackauckas,

unfortunately this is also a issue when simulationg FMUs, because we are forced to use the right root for event handling and many FMUs fire multiple events for the very same event condition. They are often implemented to e.g. change the discrete state by triggering a first event and then the continuous state by triggering another one with the very same condition. This is actually not a good way of implementing it, but it is like it is ...

I built a MWE for that, that handles four cases: Simulation of the well known bouncing ball (a) with left or right root find (this results also in slightly different affect functions of course) and (b) with (good) single conditions and a (bad) redundant double condition. In the double event case, only the SECOND event is used to change the system continuous state. All configurations work, except the double condition with right root find.

[EDIT] removed some unnecessary lines from the "M"WE ...

using DifferentialEquations
using SciMLSensitivity.SciMLBase: RightRootFind, LeftRootFind
using Plots

# Parameters for the well known bouncing ball
ENERGY_LOSS = 0.7 
RADIUS = 0.0
GRAVITY = 9.81
DBL_MIN = 2.2250738585072013830902327173324040642192159804623318306e-308
NUM_STATES = 2

t_start = 0.0
t_step = 0.05
t_stop = 2.0
tData = t_start:t_step:t_stop
u0 = [1.0, 0.0] # start state: ball position (1.0) and velocity (0.0)

solver = Tsit5()

# setup BouncingBallODE 
function fx(u, p, t)
    return [u[2], -GRAVITY] 
end

ff = ODEFunction{false}(fx) 
prob = ODEProblem{false}(ff, u0, (t_start, t_stop), [])

function condition_single(out, u, t, integrator)
    out[1] = u[1]-RADIUS
end

function condition_double(out, u, t, integrator)
    out[1] = u[1]-RADIUS
    out[2] = u[1]-RADIUS 
end

# affect for right hand event handling (affect AFTER event)
function affect_right!(integrator, idx, numIndicators)
    @info "affect_right! triggered by #$(idx)"

    if idx != numIndicators
        # event #1 is handeled as "dummy" (e.g. discrete state changes in many FMUs)
        return 
    end

    s_new = RADIUS + DBL_MIN
    v_new = -integrator.u[2]*ENERGY_LOSS
    u_new = [s_new, v_new]

    @info "New state at $(integrator.t) is $(u_new) triggered by #$(idx)"

    integrator.u .= u_new
end

# affect for left hand event handling (affect BEFORE event)
function affect_left!(integrator, idx, numIndicators)
    @info "affect_left! triggered by #$(idx)"

    if idx != numIndicators
        # event #1 is handeled as "dummy" (e.g. discrete state changes in many FMUs)
        return 
    end

    s_new = integrator.u[1]
    v_new = -integrator.u[2]*ENERGY_LOSS
    u_new = [s_new, v_new]

    @info "New state at $(integrator.t) is $(u_new)"

    integrator.u .= u_new
end

conditions = [:single, :double]
affects = [:left, :right]

for condition in conditions
    for affect in affects

        numEventIndicators = 0
        conditionFct = nothing
        rootfind = nothing
        affectFct = nothing

        if condition == :single 
            numEventIndicators = 1
            conditionFct = condition_single

        else # if condition == :double 
            numEventIndicators = 2
            conditionFct = condition_double

        end

        if affect == :left
            affectFct = affect_left!
            rootfind = LeftRootFind

        else # if affect == :right 
            affectFct = affect_right!
            rootfind = RightRootFind

        end

        cb = VectorContinuousCallback(conditionFct,
                                   (integrator, idx) -> affectFct(integrator, idx, numEventIndicators),
                                   numEventIndicators;
                                   rootfind=rootfind, save_positions=(false, false))

        solution = solve(prob; p=(), alg=solver, saveat=tData, callback=cb)
        fig = plot(solution; title="RootFind=$(affect) | Condition=$(condition)")
        display(fig)
        savefig(fig, joinpath(@__DIR__, "RootFind=$(affect)Condition=$(condition).png"))
    end
end

[ThummeTo]

Maybe I am thinking to easy, but couldn't this be solved by calling affect for all found zero crossing indices one after another? (this way it is non-allocating because no indices-vector is needed, if this is the allocation-issue you mentioned @ChrisRackauckas)

I am currently hotfixing that by making an additional condition evaluation at the beginning of affect, if there are multiple zero crossings (I just check all signs since the last condition evaluation before the event instance), I just call the affect function also for every other found zero crossing. But as you said, this is not the fine way.

[ThummeTo]

Output: