Apparent failure to catch a zero cross of a condition

[klaff]

The following Pluto notebook shows a missed zero cross of a condition.

### A Pluto.jl notebook ###
# v0.12.18

using Markdown
using InteractiveUtils

# ╔═╡ 8d429f70-5678-11eb-30d3-75441b2f27c6
begin # define independent package environment
	import Pkg
	Pkg.activate(mktempdir())
	Pkg.add("Plots"); using Plots
	Pkg.add("DifferentialEquations"); using DifferentialEquations
	Pkg.add("Parameters"); using Parameters
	Pkg.add("StaticArrays"); using StaticArrays
	Pkg.add("PlutoUI"); using PlutoUI
end

# ╔═╡ 738f95f2-567e-11eb-072f-491891d11092
md"""
###### MWE for possible missed callback by VectorContinuousCallback
"""

# ╔═╡ 16576570-5679-11eb-0ef8-e14be3d37eaa
LEFT_STOP, SLIDING_RIGHT, STOPPED, SLIDING_LEFT, RIGHT_STOP = 1:5;

# ╔═╡ bc3a1740-5678-11eb-00da-5d305b9dfb11
@with_kw mutable struct Params # do the ODE in SI units
	M::Float64 = 0.0003361185937456267
	kf::Float64 = 1488.578099595049
	Fsp::Float64 = 8.54503372291542
	xls::Float64 = 0.0
	xrs::Float64 = 0.001016
	Fsf::Float64 = 1.0 # static friction force, N
	Fkf::Float64 = 1.0 # dynamic friction force, N
	cor::Float64 = 0.3 # coefficient of restitution
	state::Int = LEFT_STOP
	A::Float64 = 7.373657906574696e-5
	minbouncespeed::Float64 = 0.001
end;

# ╔═╡ 0537d9f0-5679-11eb-1c65-1b4f21e998fd
P(t) = t>0.0 ? 6895*(18+4*sin(2π*250*t)) : 0.0;

# ╔═╡ ff171860-5678-11eb-026e-4750865dc472
function ode(u,p,t)
    v, x = u
	Params
	@unpack_Params p
    Fspring = kf*x+Fsp
	if state==LEFT_STOP 
		dv = 0.0 # not moving
		dx = 0.0 # not moving
	elseif state==SLIDING_RIGHT
		dv = 1/M*(P(t)*A-Fkf-Fspring)
		dx = v
	elseif state==STOPPED
		dv = 0.0
		dx = 0.0
	elseif state==SLIDING_LEFT
		dv = 1/M*(P(t)*A+Fkf-Fspring)
		dx = v
	else # state==RIGHT_STOP
		dv = 0.0
		dx = 0.0
	end
	@SVector [dv, dx]
end;

# ╔═╡ f7f78470-5678-11eb-2975-15d58ef4dbf0
function FricSMCondx(out, u, t, integ)
	println("[condx ] ",t," ", u)
	Params
	@unpack_Params integ.p
    v, x = u
	Fspring = kf*x+Fsp
	out[1] = (state==LEFT_STOP)*(P(t)*A-Fsf-Fspring)
	out[2] = (state==SLIDING_RIGHT)*(-v)
	out[3] = (state==SLIDING_RIGHT)*(x-xrs)
	out[4] = (state==STOPPED)*(P(t)*A-Fsf-Fspring)
	out[5] = (state==STOPPED)*(-(P(t)*A+Fsf-Fspring))
	out[6] = (state==SLIDING_LEFT)*v
	out[7] = (state==SLIDING_LEFT)*(-(x-xls))
	out[8] = (state==RIGHT_STOP)*(-(P(t)*A-Fspring+Fsf))
end;

# ╔═╡ f0b30362-5678-11eb-1ae9-4bbaf2175e41
function FricSMaffect!(integ, idx)
	print("[affect] ",integ.t,", ", integ.p.state,", ")
	Params
	@unpack_Params integ.p
	v, x = integ.u
	t = integ.t
	Fspring = Fsp+x*kf
	if     idx==1 # (state==LEFT_STOP)*(P(t)*A-Fsf-Fspring)
		integ.p.state = SLIDING_RIGHT
	elseif idx==2 # (state==SLIDING_RIGHT)*(-v)
		if P(t)*A-Fspring < -Fkf # force great enough to left that we don't stop
			integ.p.state = SLIDING_LEFT
		else
			integ.p.state = STOPPED
		end
	elseif idx==3 # (state==SLIDING_RIGHT)*(x-xrs)
		if cor==0.0 || cor*v<minbouncespeed
			integ.p.state = RIGHT_STOP
			integ.u[1] = 0.0 # set velocity to zero
		else
			integ.p.state = SLIDING_LEFT
			integ.u[1] = -cor*v 
		end
	elseif idx==4 # (state==STOPPED)*(P(t)*A-Fsf-Fspring)
		integ.p.state = SLIDING_RIGHT
	elseif idx==5 # (state==STOPPED)*(-(P(t)*A+Fsf-Fspring))
		integ.p.state = SLIDING_LEFT
	elseif idx==6 # (state==SLIDING_LEFT)*v
		if P(t)*A-Fspring > Fkf
			integ.p.state = SLIDING_RIGHT
		else
			integ.p.state = STOPPED
			integ.u[1] = 0.0
		end
	elseif idx==7 # (state==SLIDING_LEFT)*(-(x-xls))
		if cor==0.0 || -cor*v<minbouncespeed
			integ.p.state = LEFT_STOP
			integ.u[1] = 0.0
		else
			integ.p.state = SLIDING_RIGHT
			integ.u[1] = -cor*v
		end
	elseif idx==8 # (state==RIGHT_STOP)*(-(P(t)*A+Fsf-Fspring))
		integ.p.state = SLIDING_LEFT
	end
	println(join([idx, integ.p.state], ", "))
end;

# ╔═╡ 104a7a00-5679-11eb-2a5a-fbe375f59348
vcb = VectorContinuousCallback(FricSMCondx, FricSMaffect!, nothing, 
	save_positions = (true,true), 8);

# ╔═╡ 0a7910f0-5679-11eb-1deb-417a8d657bad
function sim()
	u0 = [0.0, 0.0]
	tspan = (-10e-6, 0.0025)
	p = Params()
	prob = ODEProblem(ode, u0, tspan, p)
	solve(prob, Rodas5(), callback=vcb, reltol=1e-6, abstol=1e-6, saveeverystep=true)
end;

# ╔═╡ 625b8d80-567d-11eb-34a8-5d3c85e96fa9
md"""
table of callback function calls by print statements in the condition and affect functions:

[condx ] time [velocity, position]
[affect] time prev_state, affect_idx, next_stat

"""

# ╔═╡ 218380f0-5679-11eb-2b30-7fafdb9e1dc9
with_terminal() do
	global sol = sim()
end

# ╔═╡ 024dc770-567b-11eb-3c68-693011d452c7
md"""
In the plot, we see a sign change in velocity at 0.00124.
Our state was 4 (`SLIDING_LEFT`) and was changed to 2 (`SLIDING_RIGHT`).
The change in state changes the direction of friction, and at this point in time
the driving force is ≈ 3% greater than the friction, which is why the deceleration
slope is ≈60x steeper than the acceleration slope after the cross.
So far so good.

At 0.001291 the velocity again crosses through zero. Expected result here is
that the condition `out[2] = (state==SLIDING_RIGHT)*(-v)` should trigger at that
point in time.

The interval between these times is $(round(.001291-.00124, sigdigits=3)) which is
much greater than `100eps(t)` ≈ $(round(100eps(.00124), sigdigits=3)) so I expected
this cross to be detected.

Looking at the table of condx/affect calls above, we see that immediately after the
`affect` call at 0.001239975..., the `condx` is called at the same point, then
at a slightly (maybe `100eps(t)` later point (where the velocity is positive)
then at 0.00135 where the velocity is negative. Subsequently there are 8 or 9
calls to `condx` between those points (as if it's root finding?) but no resulting
affect call.
"""

# ╔═╡ 323b7100-5679-11eb-2f4a-edcb9d254563
begin
	function custplot(sol, tspan)
		p1 = plot(sol, vars=[1], ylabel = "velocity, m/s", tspan=tspan)
		p2 = plot(sol, vars=[2], ylabel = "position, m", tspan=tspan)
		plot(p1,p2,layout=(2,1),leg=false, link=:x)
	end
	plotcutter(sol) = plotcutter(sol, sol.prob.tspan)
end

# ╔═╡ 351b8680-5679-11eb-3168-714d16feba2f
custplot(sol, (.001235, .00135))

# ╔═╡ 97f9b530-567c-11eb-0cfd-f1d95d8fa486
with_terminal() do
	Pkg.status(Pkg.PKGMODE_MANIFEST)
end

# ╔═╡ Cell order:
# ╟─738f95f2-567e-11eb-072f-491891d11092
# ╠═16576570-5679-11eb-0ef8-e14be3d37eaa
# ╠═bc3a1740-5678-11eb-00da-5d305b9dfb11
# ╠═ff171860-5678-11eb-026e-4750865dc472
# ╠═f7f78470-5678-11eb-2975-15d58ef4dbf0
# ╠═f0b30362-5678-11eb-1ae9-4bbaf2175e41
# ╠═104a7a00-5679-11eb-2a5a-fbe375f59348
# ╠═0537d9f0-5679-11eb-1c65-1b4f21e998fd
# ╠═0a7910f0-5679-11eb-1deb-417a8d657bad
# ╟─625b8d80-567d-11eb-34a8-5d3c85e96fa9
# ╠═218380f0-5679-11eb-2b30-7fafdb9e1dc9
# ╟─024dc770-567b-11eb-3c68-693011d452c7
# ╠═351b8680-5679-11eb-3168-714d16feba2f
# ╠═323b7100-5679-11eb-2f4a-edcb9d254563
# ╠═8d429f70-5678-11eb-30d3-75441b2f27c6
# ╠═97f9b530-567c-11eb-0cfd-f1d95d8fa486

[klaff]

This may be related to #632 (comment).