Wrong dense output in the event-truncated step for mass-matrix systems with SymbolicContinuousCallback

[baggepinnen]

Summary

For a mass-matrix system (one algebraic variable) with a SymbolicContinuousCallback, sol(t) queried inside the step that was truncated by the event returns wrong values — even when the event's ImperativeAffect modifies nothing. Values at the step endpoints and at the event time itself are exact, as is everything after the event; only the dense output of the single event-truncated step is affected. On the MWE below (exact solution w = 2 - t, so any interpolant should be exact) the mid-step error is 0.32.

I'm building a Coulomb-friction component whose stick/slip mode switching uses continuous events, and validating trajectories against reference data through sol(t) interpolation, which is how I ran into this.

MWE

using ModelingToolkit
using ModelingToolkit: t_nounits as t, D_nounits as D
using ModelingToolkit: SymbolicContinuousCallback, ImperativeAffect, @discretes
using OrdinaryDiffEqDefault
import ModelingToolkit.SciMLBase as SciMLBase

mode = only(@discretes mode(t) = 1.0)
@variables w(t) phi(t) sa(t) tau(t) a(t)
eqs = [
    D(phi) ~ w,
    a ~ D(w),
    a ~ ifelse(mode > 0.5, sa - 1.0, ifelse(mode < -0.5, sa + 1.0, 0.0)),
    tau ~ ifelse(mode > 0.5, 1.0, ifelse(mode < -0.5, -1.0, sa)),
    1.0 * a ~ -tau,
]
# affect intentionally changes nothing
aff = ImperativeAffect(modified = (; mode)) do m, o, c, integ
    (; mode = m.mode)
end
cb = SymbolicContinuousCallback([0 ~ w], aff; affect_neg = nothing, reinitializealg = SciMLBase.NoInit())
@named m = ModelingToolkit.System(eqs, t, [w, phi, sa, tau, a], [mode]; continuous_events = [cb])
sys = mtkcompile(m)
prob = ODEProblem(sys, [w => 2.0, phi => 0.0], (0.0, 4.0))
sol = solve(prob)
println("steps: ", sol.t)
errs = [abs(sol(tt, idxs = w) - (2 - tt)) for tt in [0.5, 1.0, 1.5]]
println("max interpolation error in event-truncated step: ", maximum(errs))

Output:

steps: [0.0, 1.0e-6, 4.4076983415698676e-5, 0.00047484681757268546, 0.0047825451591425525, 0.04785952857484122, 0.43631724946651396, 1.9999999999999643, 1.9999999999999643, 4.0]
max interpolation error in event-truncated step: 0.3184817774959323

The event (root of w) is found at exactly t = 2.0 and the solution values at all solver steps are exact (w = 2 - t). But sol(1.0, idxs = w) returns 1.318 instead of 1.0 — the query lies inside the step [0.436, 2.0] that the event truncated. The same numbers appear with solve(prob, Rodas5P(autodiff = AutoFiniteDiff())), so it is not specific to the default-algorithm composite.

What I tried

• The affect above changes nothing; adding SciMLBase.u_modified!(integ, false), or actually flipping mode, or calling SciMLBase.addsteps!(integ) in the affect — all give the identical wrong values.
• Removing the algebraic variable (pure ODE D(w) ~ ifelse(mode > 0.5, -1.0, 0.0), same callback): interpolation exact. The mass matrix seems necessary.
• An equivalent raw SciMLBase.VectorContinuousCallback on a hand-written mass-matrix ODEFunction (same three unknowns, same branch structure, affect flips the parameter), solved with the same Rodas5P(autodiff = AutoFiniteDiff()): interpolation of the truncated step is exact. So the raw solver/callback path appears fine and the difference enters with the MTK-compiled callback.
• saveat makes things worse (saved values are computed from the same interpolant); shrinking dtmax shrinks the error along with the step size.

Versions

• ModelingToolkit v11.38.1 (latest at time of writing; ModelingToolkitBase v1.60.1)
• OrdinaryDiffEqDefault v2.4.3, OrdinaryDiffEqCore v4.13.0, OrdinaryDiffEqRosenbrock v2.6.3
• DiffEqBase v7.12.0, SciMLBase v3.43.0
• Julia 1.12

Is dense output within the event-truncated step expected to be valid here? If there is a recommended way to obtain accurate sol(t) near events in the meantime, I'd appreciate a pointer — and happy to test a fix.