Allow taking derivatives of algebraic variables

[bradcarman]

The following example fails the structural_transformation() function with the error:

... Internal error in Tearing.jl: vs[1] = 0.

This examples comes from connecting a static volume to a moving volume (i.e. a pipe connected to an actuator). The pressure and density of the 2 volumes are equal and therefore the equations are essentially repeated.

using ModelingToolkit
using DifferentialEquations

@variables t
vars = @variables p1(t) p2(t) rho1(t) drho1(t) rho2(t) drho2(t) V1(t) dV1(t) x(t) dx(t) ddx(t) dm(t)
params = @parameters F A l m rho_0 bulk V2
D = Differential(t)

defs = [F => 1, A => 1, l => 1, m => 1, rho_0 => 876, bulk => 1.2e9, l=>0.1, V1=>0.1, V2=>0.1, p1=>100e5, p2=>100e5, rho1=>876*(1 + 100e5/1.2e9), rho2=>876*(1 + 100e5/1.2e9), x=>0, dx=>0, ddx=>0, dV1=>0, dm=>0, drho1=>0, drho2=>0]

eqs = [
    p1 ~ p2
    m*ddx ~ F - p1*A
    -dm ~ V1*drho1 + dV1*rho1
    +dm ~ V2*drho2
    rho1 ~ rho_0*(1 + p1/bulk)
    rho2 ~ rho_0*(1 + p2/bulk)
    D(rho1) ~ drho1
    D(rho2) ~ drho2
    D(V1) ~ dV1
    dV1 ~ dx*A
    D(x) ~ dx
    D(dx) ~ ddx
]

original = ODESystem(eqs, t, vars, params, defaults = defs)
mod = alias_elimination(original)
prob = ODEProblem(mod, [], (0.0, 0.1))
sol = solve(prob, Rodas5(),adaptive=false, dt=4e-4) #verify that the problem solves and is valid

sys = structural_simplify(original) #Error

If I simplify the equations by hand with p = p1 = p2 and rho = rho1 = rho2 giving the following then structural_transofrmation() works.

eqs = [
    m*ddx ~ F - p*A
    -dm ~ V1*drho + dV1*rho
    +dm ~ V2*drho
    rho ~ rho_0*(1 + p/bulk)
    D(rho) ~ drho
    D(V1) ~ dV1
    dV1 ~ dx*A
    D(x) ~ dx
    D(dx) ~ ddx
]

[bradcarman]

I was able to find a even simpler example which produces the Tearing error:

eqs = [
    x ~ sin(t)
    y~ 1 + x
    D(y) ~ dy
]

A workaround which resolves the error is to include a derivative of the variable (x) which is in the algebraic expression (for y) ...

eqs = [
    x ~ sin(t)
    y~ 1 + x
    D(y) ~ dy
    D(x) ~ dx
]

Therefore, it seems that the tearing algorithm sees that the variable x is involved in a derivative and is looking for a direct derivative expression for it.

[ChrisRackauckas]

That doesn't make too much sense: if you've defined y ~ ... explicitly, how is there also a differential equation D(y)? This isn't a differential equation.

[bradcarman]

Please see my pervious example. My second equation (for mass flow (dm)) is equal to the differential of density (drho) and volume (dV). Density is a function of pressure (p). Therefore I’m defining an algebraic equation (for density: rho ~ f(p)) and then taking a differential of rho. My simple example demonstrates that taking a differential of an algebraic term causes a tearing error. I would expect that one should be able to define equations in this form.

[ChrisRackauckas]

Oh I see what you're trying to do now. Yeah it's not robust to this case yet.