Non-numeric Dual{ForwardDiff.Tag...} output when solving DAE

General Usage
1

While solving a system of DAE using the integrator interface, I face a weird behavior with derivatives of the solution vector. At random steps, the resultant value (da here) is returned as non-numeric but a “Dual”. Each time that this happens the solution process freezes. BTW this behavior happens only when autodiff=true.

Anyone got an idea what’s going on?

The first few values of da (dx[2] in the DAE) printed from within Pᵥ whenever it is called are the following:

da=0.0
da=0.0
da=6.9529826154277e-310
da=3.84996053946e-313
da=Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}}(6.649515459556e-312,4.2e-322,6.64930681864e-312,6.64951549331e-312,4.45e-322,6.64930685026e-312,6.649515494497e-312,4.7e-322,6.649306858163e-312,6.64951549568e-312,4.94e-322)
da=0.0
da=0.0
da=0.0
da=0.0
da=0.0
da=0.0
da=Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}}(3.682625513e-315,1.0,3.682625513e-315,3.682625513e-315,3.682625513e-315,3.682625513e-315,3.682625513e-315,3.682625513e-315,3.682625513e-315,3.682625513e-315,3.682625513e-315)
da=-0.06745357092515131
da=-0.27408187162467584
da=-0.41785728598205224

As you can see, sometimes da is non-numeric.

The DAE system is:

function TwoPhaseFlow(dx, x, params, t)
    p, T, r = params
    Tᵢ = T[1]
    φ = porosity(x[2], x[3])
    rate = -A_g*x[7]*x[8]^(1.22)*exp(-E_g/(R_gas*x[6]))
    δ = x[2]/5
        

    dx[1] = (1.0/(ρₛ*x[2]*(1-cbrt(φ))))*(x[8] -p(t) -Pᵥ(x[2], c_radius(α(x[2], x[3]), Tᵢ, p(t)), dx[2], dx[5], x[3], α(x[2], x[3]), T, p(t), r)
            -sign(x[1])*Pᵧ(α(x[2], x[3]), T, p(t), r) -2*ρₛ*x[1]*((x[1] + dx[5]/ρₛ)*(1-cbrt(φ))) + 0.5*x[1]^2*(1-cbrt(φ)^4) -dx[5]^2*(1/ρₛ - 1/x[4]))
    dx[2] = x[1] + dx[5]/ρₛ
    dx[3] = (x[2]/x[3])^2 * (dx[2] - dx[5]/ρₛ)
    dx[4] = -3 * (x[4]/x[2]) * (dx[2] - dx[5]/x[4])
    dx[5] = (ρₛ * λₛ * Tᵢ * (R_gas*Tᵢ/(E_s)) * A_s * exp(-E_s/(R_gas*Tᵢ))/(q_s*(1-G_unreacted(Tᵢ)) - (Cp_s*(Tᵢ-Tg₀)-q_s)*log(1/G_unreacted(Tᵢ))))^0.5
    dx[6] = (1/(x[4]*Cp_mixg)) *  ( (3/x[2])*(λg*((Tᵢ - x[6])/δ) + dx[5]*Cp_gR*(Tᵢ-x[6])) - rate*(ΔhfRg - ΔhfP2 + (Cp_gR - Cp_gP2)*(x[6]-Tg₀)) + (1-η_abel*x[4])*dx[8]  )
    dx[7] = rate/x[4] + (3/(x[2]*x[4]))*dx[5]*(G_unreacted(Tᵢ) - x[7])
    dx[8] = x[4]*(R_gas_air)*x[6] / (1-η_abel*x[4]) - x[8]
    dx[9] = -rate/x[4] + (3/(x[2]*x[4]))*dx[5]*(1.0 - G_unreacted(Tᵢ) - x[9])
    dx[10] = -(3/(x[2]*x[4]))*dx[5]*x[10]
    
    nothing
end
2

That’s from automatic differentiation. For more information, see Forward-Mode Automatic Differentiation (AD) via High Dimensional Algebras - MIT Parallel Computing and Scientific Machine Learning (SciML)

3

Thanks, I’ve gone through the provided information, however I’m still struggling to figure out how I should solve the problem.

4

It’s supposed to use these Dual numbers to calculate derivatives via automatic differentiation. It’s expected behavior.

Do you expect the solve to go fast (based on what information)? Is your problem well conditioned? (What are the param values?)

5

What is the problem to solve? The only question you’ve asked is what are those Dual values, and those are the numeric values during automatic differentiation. You haven’t posted any issues though with it beyond not knowing what the Duals are?

6

Each time that this happens the solution process freezes.

However, now according to your replies, I realize that the freezing is probably not due to autodiff or dual numbers, as per your suggestion the output seems to be normal.

I’ll investigate further.

Thanks!

7

Freezing, what do you mean by freezing? After auto-diff you’ll have to do a linear solve, since that means it’s an implicit method. The linear solve is usually the most expensive operation. However, it shouldn’t “freeze” for a 10x10 matrix.

Can you give an MWE that displays the issue?

8

Hi Chris, while playing around with the code, I think I have found what the problem is. It’s the sign(x[1]) which multiplies
Pᵧ(α(x[2], x[3]) makes the whole calculation to hang indefinitely. If I remove the sign function, the code runs fine with or without autodiff.

The sign function is needed because it makes Pᵧ oppose the action of the driving pressure, i.e. going from contraction to expansion. However, I guess there’s some instability going on when the sign change occurs.

Is there a better way to do it?

9

Using the sign function is fine for AD, but is your f discontinuous? If so, then the DAE doesn’t necessarily have a smooth solution which would cause the solver to have issues. You would need to handle that with a continuous callback.

1 Like
10

Sorry to bother with this again! I have used a continuous callback to flip the sign of a new parameter, mysign, which then multiplies Pᵧ, however I get the following error after some time:

MethodError: no method matching setindex!(::Tuple{var"#3#6", Vector{Float64}, Vector{Float64}, Float64}, ::Float64, ::Int64)

This is the relevant part of the code:

# solve system of ODE
r0 = Vector(range(a₀, stop=b₀, length=nᵣ))
mysign = 1.0
params = (P, T, r0, mysign)
M = [1 0 0 0 0 0 0 0 0 0
     0 1 0 0 0 0 0 0 0 0
     0 0 1 0 0 0 0 0 0 0
     0 0 0 1 0 0 0 0 0 0
     0 0 0 0 1 0 0 0 0 0
     0 0 0 0 0 1 0 0 0 0
     0 0 0 0 0 0 1 0 0 0
     0 0 0 0 0 0 0 0 0 0
     0 0 0 0 0 0 0 0 1 0
     0 0 0 0 0 0 0 0 0 1]

function condition(u, t, integrator)
    u[1]
end
    
function affect!(integrator)
    integrator.p[4] = -integrator.p[4]
end    

cb = ContinuousCallback(condition, affect!)
    
f = ODEFunction(TwoPhaseFlow, mass_matrix = M)
problem = ODEProblem(f, [da₀; a₀; b₀; ρg₀; m₀; Tg₀; YբR₀; Pg₀; YբP₀; YբI₀], tspan, params)
integrator = init(problem, FBDF(autodiff=true), reltol=1e-6, maxiters=1e8, callback = cb)

tstart = tspan[1]
tstep = 0
       
function run_loop(tstep, T, tstart)
    while (tstart < tspan[2])
        tstep += 1
        T_solution = solve_Tsolid(T, (tstart, tstep, exo), P(tstart))
        R = T_solution[1]
        T = T_solution[2]
        integrator.p = params = (P, T, Vector(range(integrator.sol(tstart)[2], stop=integrator.sol(tstart)[3], length=nᵣ)), mysign)
        step!(integrator, dt, true)
        T_of_R = Interpolation(T, R)
        tstart += dt
    end
end
run_loop(0, T, tstart)

I’m not sure what the error message means in this regard.

11

Tuples are immutable. You can use Setfield.jl to work around this.

12

Thanks vettert!
Can you please elaborate on how to use it in my case?
I tried with:

function affect!(integrator)
    @set integrator.p[4] = -integrator.p[4]
end

but that resulted in an error.

13

Update:
Instead of using the sign(x) function, I used tanh(x), which effectively does the same, but does not introduce a discontinuity and has controlled steepness to simulate sign(x) as closely as needed (multiply the argument x with a suitably large number). Problem solved. No callbacks or mutation tricks needed.