Adaptivity fails for DAEs with non-continuous solution

[termi-official]

Describe the bug 🐞

I think this is technically not a bug, but the title should explain the issue.

Expected behavior

A clear and concise description of what you expected to happen.

Minimal Reproducible Example 👇

It follows a simple stokes problem in 1D with different boundary conditions.
A) ramping up the inflow velocity linearly from t=0 to t=2 -> breaks adaptivity at t=2, but problem when adaptivity is turned off
B) ramping up the inflow velocity smoothly from t=0 to t=2 -> works adaptively and without adaptivity

Just comment in and out the first few lines in the combinations described above to reproduce the different behaviors.

using OrdinaryDiffEq, UnPack

# Boundary condition - Issue 1
vᵢₙ(t) = min(t*1.5/2.0, 1.5)
adaptive = true  # fails at t=2.0
# adaptive = false # works

# vᵢₙ(t) = t < 2.0 ? 1.5*(sin(-π/2 + t*π/2)+1)/2 : 1.5
# adaptive = false # works
# adaptive = true # works

T = 6.0
Δt₀ = 0.1
Δt_save = 0.1

# Mass matrix and stiffness matrix for stokes problem
M = [0.07272727272727277 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.08888888888888902 0.022222222222222303 0.0 0.0 -0.01111111111111112 0.022222222222222185 0.0 0.0 0.0 0.0; 0.0 0.022222222222222303 0.17777777777777778 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 -0.01111111111111112 0.0 0.0 0.0 0.08888888888888902 0.022222222222222303 0.0 -0.011111111111111124 0.02222222222222219 0.0; 0.0 0.022222222222222185 0.0 0.0 0.0 0.022222222222222303 0.17777777777777778 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 -0.011111111111111124 0.0 0.0 0.04444444444444456 0.0222222222222223 0.0; 0.0 0.0 0.0 0.0 0.0 0.02222222222222219 0.0 0.0 0.0222222222222223 0.17777777777777778 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0]
K = [-0.007000000000000001 -0.001000000000000003 0.008000000000000004 -0.8333333333333333 -0.16666666666666624 0.0 0.0 0.0 0.0 0.0 0.0; -0.001000000000000003 -0.014000000000000016 0.008000000000000018 0.16666666666666682 1.7763568394002505e-15 -0.001000000000000003 0.008000000000000004 -0.16666666666666624 0.0 0.0 0.0; 0.008000000000000004 0.008000000000000018 -0.01600000000000002 0.6666666666666663 -0.6666666666666685 0.0 0.0 0.0 0.0 0.0 0.0; -0.8333333333333333 0.16666666666666682 0.6666666666666663 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; -0.16666666666666624 1.7763568394002505e-15 -0.6666666666666685 0.0 0.0 0.16666666666666682 0.6666666666666663 0.0 0.0 0.0 0.0; 0.0 -0.001000000000000003 0.0 0.0 0.16666666666666682 -0.014000000000000014 0.008000000000000018 1.6653345369377348e-15 -0.0010000000000000033 0.008000000000000002 -0.16666666666666624; 0.0 0.008000000000000004 0.0 0.0 0.6666666666666663 0.008000000000000018 -0.01600000000000002 -0.6666666666666685 0.0 0.0 0.0; 0.0 -0.16666666666666624 0.0 0.0 0.0 1.6653345369377348e-15 -0.6666666666666685 0.0 0.16666666666666682 0.6666666666666664 0.0; 0.0 0.0 0.0 0.0 0.0 -0.0010000000000000035 0.0 0.16666666666666682 -0.007000000000000015 0.008000000000000018 0.833333333333335; 0.0 0.0 0.0 0.0 0.0 0.008000000000000002 0.0 0.6666666666666664 0.008000000000000018 -0.01600000000000002 -0.6666666666666687; 0.0 0.0 0.0 0.0 0.0 -0.16666666666666624 0.0 0.0 0.833333333333335 -0.6666666666666687 0.0]
jac_sparsity = copy(K);

# Consistent initial condition
u₀ = zeros(11)

# RHS helper
struct RHSparams
   K::Matrix
   u::Vector
end
p = RHSparams(K, copy(u₀))

function ferrite_limiter!(u, _, p, t)
    u[1] = vᵢₙ(t)
end

function stokes!(du,u_uc,p::RHSparams,t)
    @unpack K,u = p

    u .= u_uc
    u[1] = vᵢₙ(t) # Boundary condition

    du .= K * u

    # du[1] = 0.0 # Prescribed, so what to do here?
end;


function stokes_jac!(J,u_uc,p,t)
   @unpack K, u = p

   u .= u_uc
   u[1] = vᵢₙ(t) # Boundary condition

   J .= K
   # Eliminate dof
   J[1,:] .= 0.0
   J[:,1] .= 0.0
   J[1,1] = 1.0
end;

rhs = ODEFunction(stokes!, mass_matrix=M; jac=stokes_jac!, jac_prototype=jac_sparsity)
problem = ODEProblem(rhs, u₀, (0.0,T), p);

struct FreeDofErrorNorm
end
(::FreeDofErrorNorm)(u::AbstractArray, t) = DiffEqBase.ODE_DEFAULT_NORM(u[2:end], t)

timestepper = Rodas5P(autodiff=false, step_limiter! = ferrite_limiter!);

using Logging
debuglogger = ConsoleLogger(stderr, Logging.Debug)
# with_logger(debuglogger) do
    integrator = init(
        problem, timestepper, initializealg=NoInit(), dt=Δt₀,
        adaptive=adaptive, abstol=1e-4, reltol=1e-5,
        progress=true, progress_steps=1,
        verbose=true
    );

    for (u,t) in tuples(integrator)
        @info t
        velocity_dof_indices = [1,2,3,6,7,9,10]
        @assert all(u[velocity_dof_indices] .≤ 2.0) "$(u[velocity_dof_indices])"
    end
# end

Error & Stacktrace ⚠️

...
[ Info: 1.9999991627202856
[ Info: 1.9999996552983934
[ Info: 1.9999998940463428
[ Info: 1.9999999317284627
[ Info: 1.9999999735509675
┌ Warning: dt(2.220446049250313e-16) <= eps(t)(1.9999999735509675) , and step error estimate = 15299.784213961508. Aborting. There is either an error in your model specification or the true solution is unstable (or the true solution can not be represented in the precision of Float64).
└ @ SciMLBase ~/.julia/packages/SciMLBase/JUp1I/src/integrator_interface.jl:632
[ Info: 1.9999999735509675

Environment (please complete the following information):

• Output of using Pkg; Pkg.status()

  [1dea7af3] OrdinaryDiffEq v6.80.1

[oscardssmith]

This is expected behavior. Setting d_discontinuities=[2.0] as a kwarg to init to tell the solver where the discontinuity occurs (or using a callback) solves it.