`Rodas5P` Goes unstable with an increase in `abstol`

[bradcarman]

Take...

using ModelingToolkit
using ModelingToolkit: t_nounits as t, D_nounits as D
using OrdinaryDiffEqRosenbrock

function transition(x1, x2, y1, y2, x)
    u = (x - x1) / (x2 - x1)
    blend = u^2 * (3 - 2 * u)
    return (1 - blend) * y1 + blend * y2
end

liquid_density(p, rho_0, beta) = rho_0*(1 + p/beta)

full_density(p, rho_0, rho_gas, beta, p_0, p_gas) = ifelse(p>p_0, liquid_density(p, rho_0, beta), transition(p_0, p_gas, rho_0, rho_gas, p))

function System(;name)
    pars = @parameters begin
        rho_0 = 1000
        
        p_0 = 0
        p_gas = -1000

        rho_gas = 1
        beta = 2e9
        A = 0.1
        m = 100
        L = 0.01

        C = 1.35
        c = 1000

        force = 3e5
    end   
    vars = @variables begin
        p_1(t)=10e5
        p_2(t)=10e5
        x(t)=0
        dx(t)=0
        ddx(t), [guess=0]
        rho_1(t), [guess=rho_0]
        rho_2(t), [guess=rho_0]
        m_1(t), [guess=0]
        m_2(t), [guess=0]
    end
    
    eqs = [

        D(x) ~ dx
        D(dx) ~ ddx

        m_1 ~ rho_1*(L+x)*A
        m_2 ~ rho_2*(L-x)*A

        rho_1 ~ full_density(p_1, rho_0, rho_gas, beta, p_0, p_gas)
        rho_2 ~ full_density(p_2, rho_0, rho_gas, beta, p_0, p_gas)

        m*ddx ~ (p_1 - p_2)*A - c*dx + force*sin(2*pi*t/0.01)

        D(m_1) ~ 0
        D(m_2) ~ 0

    ]

    

    return ODESystem(eqs, t, vars, pars; name)
end

@mtkbuild sys = System() 
prob = ODEProblem(sys, [], (0, 0.01))

sol = solve(prob, Rodas5P(); abstol=1e-6, reltol=1e-9);  # Success
sol = solve(prob, Rodas5P(); abstol=1e-5, reltol=1e-9);  # Success
sol = solve(prob, Rodas5P(); abstol=1e-4, reltol=1e-9);  # Unstable

Note that increasing the abstol causes the solution to go unstable.

[gstein3m]

Equations

        rho_1 ~ full_density(p_1, rho_0, rho_gas, beta, p_0, p_gas)
        rho_2 ~ full_density(p_2, rho_0, rho_gas, beta, p_0, p_gas)

must be used to calculate p_1 and p_2 in the DAE system. If the derivative of rho with respect to p is zero, the index-1 condition is violated. This is the case with the transition function in the vicinity of p=p_0. Also the gradient in the liquid function is very small.
It is therefore more advantageous to formulate a function p(rho) instead of rho(p):

        p_1 ~ fp(rho_1,rho_0, rho_gas, beta, p_0, p_gas)
        p_2 ~ fp(rho_2,rho_0, rho_gas, beta, p_0, p_gas)

with
fp(rho, rho_0, rho_gas, beta, p_0, p_gas) = ifelse(rho > rho_0, beta*(rho/rho_0 - 1), p_0 + (p_gas - p_0)*(rho - rho_0)/(rho_gas - rho_0))

[bradcarman]

@gstein3m thank you so much!!! This is the magic that was needed that I've been struggling with for years. Here is the reformulation...

using ModelingToolkit
using ModelingToolkit: t_nounits as t, D_nounits as D
using OrdinaryDiffEqRosenbrock
using ToleranceTools
using Plots

beta(rho, rho_0, beta_liquid, beta_gas) = ifelse(rho > rho_0, beta_liquid, beta_gas)
full_pressure(rho, p_0, rho_0, beta_liquid, beta_gas) = (rho/rho_0 - 1)*beta(rho, rho_0, beta_liquid, beta_gas) + p_0


function System(;name)
    pars = @parameters begin
        rho_0 = 1000
        p_0 = 0
        beta_liquid = 2e9
        beta_gas = 1e3
        A = 0.1
        m = 100
        L = 0.01
        c = 1000
        force = 3e5
    end   
    vars = @variables begin
        p_1(t), [guess=10e5]
        p_2(t), [guess=10e5]
        x(t)=0
        dx(t)=0
        ddx(t), [guess=0]
        rho_1(t)=rho_0*(1 + 10e5/beta_liquid)
        rho_2(t)=rho_0*(1 + 10e5/beta_liquid)
        m_1(t), [guess=0]
        m_2(t), [guess=0]
    end
    eqs = [
        D(x) ~ dx
        D(dx) ~ ddx

        m_1 ~ rho_1*(L+x)*A
        m_2 ~ rho_2*(L-x)*A

        p_1 ~ full_pressure(rho_1, p_0, rho_0, beta_liquid, beta_gas)
        p_2 ~ full_pressure(rho_2, p_0, rho_0, beta_liquid, beta_gas)

        m*ddx ~ (p_1 - p_2)*A - c*dx + force*sin(2*pi*t/0.01)

        D(m_1) ~ 0
        D(m_2) ~ 0
    ]
    return ODESystem(eqs, t, vars, pars; name)
end

@mtkbuild sys = System() 
prob = ODEProblem(sys, [], (0, 0.03))
sol = solve(prob, Rodas5P(); abstol=1e-9, reltol=1e-9) #Success

With this new formulation I can use the tolerance settings and get reasonable results as I adjust abstol and reltol. My package ToleranceTools can be used to check the residual across the full grid from 1e-3 to 1e-9...

res = ToleranceTools.analysis(prob, Rodas5P())
savefig(res, "rodas.png")