Model Regression from v6.74.1 to v6.83.1

[bradcarman]

The following model worked previously but is now failing to solve with message: At t=0.0, dt was forced below floating point epsilon 5.0e-324, and step error estimate = 1.5295211089754415e8. Aborting.

Using non-adaptive ImplicitEuler gives the correct solution.

using ModelingToolkit
using DifferentialEquations
using Symbolics
using ModelingToolkit: t_nounits as t, D_nounits as D

# parameters -------
pars = @parameters begin
    r₀ = 876 #kg/m^3
    β = 1.2e9 #Pa
    A = 0.01 #m²
    x₀ = 1.0 #m
    M = 10_000 #kg
    g = 9.807 #m/s²
    amp = 5e-2 #m
    f = 15 #Hz    
end

dt = 1e-4 #s
t_end = 0.2 #s
time = 0:dt:t_end

x_fun(t,amp,f) = amp*sin(2π*t*f) + x₀

ẋ_fun = build_function(expand_derivatives(D(x_fun(t,amp,f))), t, amp, f; expression=false)

vars = @variables begin
    x(t) = x₀
    ẋ(t)
    ẍ(t)
    p(t) = M*g/A #Pa
    ṁ(t)
    r(t)
    ṙ(t)
end 

function get_base_equations(density_type) 

    eqs = [
        D(x) ~ ẋ 
        D(ẋ) ~ ẍ
        D(r) ~ ṙ

        r ~ r₀*(1 + p/β)

        ṁ ~ ṙ*x*A + (density_type)*ẋ*A
        M*ẍ ~ p*A - M*g
    ]

    return eqs
end
nothing # hide

eqs_x = [
    get_base_equations(r)...
    x ~ x_fun(t,amp,f) # (4) Input - target x 
]

@mtkbuild odesys_x = ODESystem(eqs_x, t, vars, pars)

prob_x = ODEProblem(odesys_x, [], (0, t_end))
sol_x = solve(prob_x; saveat=time)

eqs_ṁ1 = [
    get_base_equations(r₀)...
    ṁ ~ ẋ_fun(t,amp,f)*A*r # (4) Input - mass flow guess
]

@mtkbuild odesys_ṁ1 = ODESystem(eqs_ṁ1, t, vars, pars)
u0 = [sol_x[s][1] for s in unknowns(odesys_ṁ1)]
prob_ṁ1 = ODEProblem(odesys_ṁ1, u0, (0, t_end))

# BUG: At t=0.0, dt was forced below floating point epsilon 5.0e-324, and step error estimate = 1.5295211089754415e8. Aborting.
# Fails with OrdinaryDiffEq v6.83.1, works with v6.74.1
@time sol_ṁ1 = solve(prob_ṁ1; initializealg=NoInit()); # BUG
@time sol_ṁ1 = solve(prob_ṁ1, ImplicitEuler(); initializealg=NoInit(), adaptive=false, dt); #OK

[oscardssmith]

It appears that you have initialized your system incorrectly. If you tell it to initialize with BrownFullBasicInit() it works. Specifically, the problem is that your u0 has u0[4]==-30556.685668435468 but for your algebraic equations to be satisfied at t=0, you need it to be 33.736511073364454

julia> sol_ṁ1 = solve(prob_ṁ1, FBDF(), initializealg=BrownFullBasicInit())
retcode: Success
Interpolation: 3rd order Hermite
t: 160-element Vector{Float64}:

[ChrisRackauckas]

Consider any case of initializealg=NoInit() a user error. If you do not initialize the DAE properly, it will not necessarily be solvable!