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

**URL:** <https://discourse.julialang.org/t/non-numeric-dual-forwarddiff-tag-output-when-solving-dae/109589>\
**Category:** General Usage\
**Created:** [February 1, 2024, 7:04pm UTC](https://discourse.julialang.org/t/non-numeric-dual-forwarddiff-tag-output-when-solving-dae/109589 "2024-02-01T19:04:05Z")\
**Posts on this page:** 13\
**Page:** 1

<div class="post-metadata">

**Author:** ![David\_F](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/david_f/32/42706_2.png) [@David\_F](https://discourse.julialang.org/u/David_F)\
**Post date:** [February 1, 2024, 7:04pm UTC](https://discourse.julialang.org/t/non-numeric-dual-forwarddiff-tag-output-when-solving-dae/109589/1 "2024-02-01T19:04:05Z")

</div>

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:

```julia
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

```

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [February 1, 2024, 9:03pm UTC](https://discourse.julialang.org/t/non-numeric-dual-forwarddiff-tag-output-when-solving-dae/109589/2 "2024-02-01T21:03:03Z")

</div>

> [@David\_F](#):
>
> 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.

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)](https://book.sciml.ai/notes/08-Forward-Mode_Automatic_Differentiation_(AD)_via_High_Dimensional_Algebras/)

---

<div class="post-metadata">

**Author:** ![David\_F](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/david_f/32/42706_2.png) [@David\_F](https://discourse.julialang.org/u/David_F)\
**Post date:** [February 1, 2024, 9:13pm UTC](https://discourse.julialang.org/t/non-numeric-dual-forwarddiff-tag-output-when-solving-dae/109589/3 "2024-02-01T21:13:34Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![vettert](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/vettert/32/30599_2.png) [@vettert](https://discourse.julialang.org/u/vettert)\
**Post date:** [February 1, 2024, 9:21pm UTC](https://discourse.julialang.org/t/non-numeric-dual-forwarddiff-tag-output-when-solving-dae/109589/4 "2024-02-01T21:21:02Z")

</div>

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?)

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [February 1, 2024, 11:02pm UTC](https://discourse.julialang.org/t/non-numeric-dual-forwarddiff-tag-output-when-solving-dae/109589/5 "2024-02-01T23:02:01Z")

</div>

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?

---

<div class="post-metadata">

**Author:** ![David\_F](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/david_f/32/42706_2.png) [@David\_F](https://discourse.julialang.org/u/David_F)\
**Post date:** [February 2, 2024, 11:19am UTC](https://discourse.julialang.org/t/non-numeric-dual-forwarddiff-tag-output-when-solving-dae/109589/6 "2024-02-02T11:19:41Z")

</div>

> [@ChrisRackauckas](#):
>
> What is the problem to solve?

> 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!

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [February 2, 2024, 11:21am UTC](https://discourse.julialang.org/t/non-numeric-dual-forwarddiff-tag-output-when-solving-dae/109589/7 "2024-02-02T11:21:19Z")

</div>

> [@David\_F](#):
>
> 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.

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?

---

<div class="post-metadata">

**Author:** ![David\_F](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/david_f/32/42706_2.png) [@David\_F](https://discourse.julialang.org/u/David_F)\
**Post date:** [February 4, 2024, 7:18pm UTC](https://discourse.julialang.org/t/non-numeric-dual-forwarddiff-tag-output-when-solving-dae/109589/8 "2024-02-04T19:18:45Z")

</div>

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?

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [February 5, 2024, 2:45am UTC](https://discourse.julialang.org/t/non-numeric-dual-forwarddiff-tag-output-when-solving-dae/109589/9 "2024-02-05T02:45:56Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![David\_F](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/david_f/32/42706_2.png) [@David\_F](https://discourse.julialang.org/u/David_F)\
**Post date:** [February 5, 2024, 9:10am UTC](https://discourse.julialang.org/t/non-numeric-dual-forwarddiff-tag-output-when-solving-dae/109589/10 "2024-02-05T09:10:55Z")

</div>

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:

```julia
# 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.

---

<div class="post-metadata">

**Author:** ![vettert](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/vettert/32/30599_2.png) [@vettert](https://discourse.julialang.org/u/vettert)\
**Post date:** [February 5, 2024, 9:41am UTC](https://discourse.julialang.org/t/non-numeric-dual-forwarddiff-tag-output-when-solving-dae/109589/11 "2024-02-05T09:41:13Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![David\_F](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/david_f/32/42706_2.png) [@David\_F](https://discourse.julialang.org/u/David_F)\
**Post date:** [February 5, 2024, 11:49am UTC](https://discourse.julialang.org/t/non-numeric-dual-forwarddiff-tag-output-when-solving-dae/109589/12 "2024-02-05T11:49:24Z")

</div>

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

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

```

but that resulted in an error.

---

<div class="post-metadata">

**Author:** ![David\_F](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/david_f/32/42706_2.png) [@David\_F](https://discourse.julialang.org/u/David_F)\
**Post date:** [February 7, 2024, 6:28am UTC](https://discourse.julialang.org/t/non-numeric-dual-forwarddiff-tag-output-when-solving-dae/109589/13 "2024-02-07T06:28:41Z")

</div>

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.