Solving a nonlinear equation using NonlinearSolve.jl with AutoFiniteDiff() yields errors about dual number

Specific Domains Modelling & Simulations
,
1

I am trying to solve a nonlinear equation, which is constructed using the thermophysical properties of fluids (via CoolProp package).

An MWE is shown below.

using NonlinearSolve, CoolProp

function myfun_P_MWE(p, (T, init_m))
    r_cont_ads = 10e-3
    r_ads_heat = 2.5e-3
    L = 0.5
    row_GGHS = 3
    row_cont = 3
    row_ads = 10
    delta_r_ads = (r_cont_ads - r_ads_heat) / row_ads

    rho_fre = zeros(typeof(p), row_ads)
    for i in eachindex(rho_fre)
        # Obtain the density of nitrogen according to temperature and pressure
        rho_fre[i] = PropsSI("D", "T", T[row_GGHS+row_cont+i], "P", p, "nitrogen")
    end

    # Construct the nonlinear equation
    f = 0.0
    for i in eachindex(rho_fre)
        f = f + rho_fre[i] * ((r_cont_ads - (i - 1) * delta_r_ads)^2 - (r_cont_ads - i * delta_r_ads)^2) / 2 * L
    end
    return f - init_m
end


# Solving the equation
T = collect(range(200.0, 500.0, 19))
p = 1e5
init_m = 2e-4

# myfun_P_MWE(p,(T, init_m))
p_new = solve(
    NonlinearProblem(myfun_P_MWE, p, (T, init_m))
    # RobustMultiNewton(; autodiff=AutoFiniteDiff())
).u

In this case, I am trying to solve the pressure (p) corresponding to a given temperature array (T) and mass (init_m). The error message is as follows.

ERROR: MethodError: no method matching Float64(::ForwardDiff.Dual{ForwardDiff.Tag{SciMLBase.JacobianWrapper{…}, Float64}, Float64, 1})

Closest candidates are:
  (::Type{T})(::Real, ::RoundingMode) where T<:AbstractFloat
   @ Base rounding.jl:207
  (::Type{T})(::T) where T<:Number
   @ Core boot.jl:792
  Float64(::IrrationalConstants.Sqrt2)
   @ IrrationalConstants C:\Users\qlx\.julia\packages\IrrationalConstants\vp5v4\src\macro.jl:112
  ...

Stacktrace:
  [1] convert(::Type{Float64}, x::ForwardDiff.Dual{ForwardDiff.Tag{SciMLBase.JacobianWrapper{…}, Float64}, Float64, 1})
    @ Base .\number.jl:7
  [2] cconvert(T::Type, x::ForwardDiff.Dual{ForwardDiff.Tag{SciMLBase.JacobianWrapper{…}, Float64}, Float64, 1})
    @ Base .\essentials.jl:543
  [3] PropsSI(output::String, name1::String, value1::Float64, name2::String, value2::ForwardDiff.Dual{…}, fluid::String)
    @ CoolProp C:\Users\qlx\.julia\packages\CoolProp\RDEcq\src\CoolProp.jl:81
  [4] myfun_P_MWE(p::ForwardDiff.Dual{ForwardDiff.Tag{…}, Float64, 1}, ::Tuple{Vector{…}, Float64})
    @ Main c:\Users\qlx\Desktop\RB4_AP4_multi_stage_Julia_CoolProp\test2.jl:15
  [5] NonlinearFunction
    @ C:\Users\qlx\.julia\packages\SciMLBase\NjslX\src\scimlfunctions.jl:2185 [inlined]
  [6] JacobianWrapper
    @ C:\Users\qlx\.julia\packages\SciMLBase\NjslX\src\function_wrappers.jl:103 [inlined]
  [7] __value_derivative
    @ C:\Users\qlx\.julia\packages\NonlinearSolve\R8LAH\src\internal\jacobian.jl:186 [inlined]
  [8] JacobianCache
    @ C:\Users\qlx\.julia\packages\NonlinearSolve\R8LAH\src\internal\jacobian.jl:118 [inlined]
  [9] JacobianCache
    @ C:\Users\qlx\.julia\packages\NonlinearSolve\R8LAH\src\internal\jacobian.jl:106 [inlined]
 [10] __step!(cache::NonlinearSolve.GeneralizedFirstOrderAlgorithmCache{…}; recompute_jacobian::Nothing, kwargs::@Kwargs{})
    @ NonlinearSolve C:\Users\qlx\.julia\packages\NonlinearSolve\R8LAH\src\core\generalized_first_order.jl:208
 [11] __step!
    @ C:\Users\qlx\.julia\packages\NonlinearSolve\R8LAH\src\core\generalized_first_order.jl:204 [inlined]
 [12] #step!#210
    @ C:\Users\qlx\.julia\packages\NonlinearSolve\R8LAH\src\core\generic.jl:55 [inlined]
 [13] step!
    @ C:\Users\qlx\.julia\packages\NonlinearSolve\R8LAH\src\core\generic.jl:50 [inlined]
 [14] solve!(cache::NonlinearSolve.GeneralizedFirstOrderAlgorithmCache{…})
    @ NonlinearSolve C:\Users\qlx\.julia\packages\NonlinearSolve\R8LAH\src\core\generic.jl:13
 [15] #__solve#209
    @ C:\Users\qlx\.julia\packages\NonlinearSolve\R8LAH\src\core\generic.jl:4 [inlined]
 [16] __solve
    @ C:\Users\qlx\.julia\packages\NonlinearSolve\R8LAH\src\core\generic.jl:1 [inlined]
 [17] macro expansion
    @ C:\Users\qlx\.julia\packages\NonlinearSolve\R8LAH\src\default.jl:278 [inlined]
 [18] __solve(::NonlinearProblem{…}, ::NonlinearSolvePolyAlgorithm{…}; alias_u0::Bool, verbose::Bool, kwargs::@Kwargs{…})
    @ NonlinearSolve C:\Users\qlx\.julia\packages\NonlinearSolve\R8LAH\src\default.jl:248
 [19] __solve
    @ C:\Users\qlx\.julia\packages\NonlinearSolve\R8LAH\src\default.jl:248 [inlined]
 [20] #__solve#329
    @ C:\Users\qlx\.julia\packages\NonlinearSolve\R8LAH\src\default.jl:511 [inlined]
 [21] __solve
    @ C:\Users\qlx\.julia\packages\NonlinearSolve\R8LAH\src\default.jl:508 [inlined]
 [22] #__solve#72
    @ C:\Users\qlx\.julia\packages\DiffEqBase\O8cUq\src\solve.jl:1394 [inlined]
 [23] __solve
    @ C:\Users\qlx\.julia\packages\DiffEqBase\O8cUq\src\solve.jl:1386 [inlined]
 [24] #solve_call#44
    @ C:\Users\qlx\.julia\packages\DiffEqBase\O8cUq\src\solve.jl:612 [inlined]
 [25] solve_call
    @ C:\Users\qlx\.julia\packages\DiffEqBase\O8cUq\src\solve.jl:569 [inlined]
 [26] #solve_up#53
    @ C:\Users\qlx\.julia\packages\DiffEqBase\O8cUq\src\solve.jl:1072 [inlined]
 [27] solve_up
    @ C:\Users\qlx\.julia\packages\DiffEqBase\O8cUq\src\solve.jl:1066 [inlined]
 [28] #solve#52
    @ C:\Users\qlx\.julia\packages\DiffEqBase\O8cUq\src\solve.jl:1060 [inlined]
 [29] solve(::NonlinearProblem{…})
    @ DiffEqBase C:\Users\qlx\.julia\packages\DiffEqBase\O8cUq\src\solve.jl:1050
 [30] top-level scope
    @ c:\Users\qlx\Desktop\RB4_AP4_multi_stage_Julia_CoolProp\test2.jl:33
Some type information was truncated. Use `show(err)` to see complete types.

It looks like the error is because PropsSI() cannot accept a dual number as an argument. Following the FAQ: The solver tried to set a Dual Number in my Vector of Floats. How do I fix that? in the NonlinearSolve.jl documentation, I have initialized the rho_fre with the same element type as p, but that doesn’t help.

I am not familiar with dual number, by google what it is, I think it is used when automatically differentiating the objective function. Therefore, I tried another solution in the FAQ, which is to specify the autodiff to be AutoFiniteDiff(). The revised code is as follows, but it yields the same error message, complaining about the dual number usage in PropsSI().

using NonlinearSolve, CoolProp

function myfun_P_MWE(p, (T, init_m))
    r_cont_ads = 10e-3
    r_ads_heat = 2.5e-3
    L = 0.5
    row_GGHS = 3
    row_cont = 3
    row_ads = 10
    delta_r_ads = (r_cont_ads - r_ads_heat) / row_ads

    rho_fre = zeros(typeof(p), row_ads)
    for i in eachindex(rho_fre)
        # Obtain the density of nitrogen according to temperature and pressure
        rho_fre[i] = PropsSI("D", "T", T[row_GGHS+row_cont+i], "P", p, "nitrogen")
    end

    # Construct the nonlinear equation
    f = 0.0
    for i in eachindex(rho_fre)
        f = f + rho_fre[i] * ((r_cont_ads - (i - 1) * delta_r_ads)^2 - (r_cont_ads - i * delta_r_ads)^2) / 2 * L
    end
    return f - init_m
end


# Solving the equation
T = collect(range(200.0, 500.0, 19))
p = 1e5
init_m = 2e-4

# myfun_P_MWE(p,(T, init_m))
p_new = solve(
    NonlinearProblem(myfun_P_MWE, p, (T, init_m)),
    RobustMultiNewton(; autodiff=AutoFiniteDiff())
).u

ERROR: MethodError: no method matching Float64(::ForwardDiff.Dual{ForwardDiff.Tag{SciMLBase.JacobianWrapper{…}, Float64}, Float64, 1})

Closest candidates are:
  (::Type{T})(::Real, ::RoundingMode) where T<:AbstractFloat
   @ Base rounding.jl:207
  (::Type{T})(::T) where T<:Number
   @ Core boot.jl:792
  Float64(::IrrationalConstants.Sqrt2)
   @ IrrationalConstants C:\Users\qlx\.julia\packages\IrrationalConstants\vp5v4\src\macro.jl:112
  ...

Stacktrace:
  [1] convert(::Type{Float64}, x::ForwardDiff.Dual{ForwardDiff.Tag{SciMLBase.JacobianWrapper{…}, Float64}, Float64, 1})
    @ Base .\number.jl:7
  [2] cconvert(T::Type, x::ForwardDiff.Dual{ForwardDiff.Tag{SciMLBase.JacobianWrapper{…}, Float64}, Float64, 1})
    @ Base .\essentials.jl:543
  [3] PropsSI(output::String, name1::String, value1::Float64, name2::String, value2::ForwardDiff.Dual{…}, fluid::String)
    @ CoolProp C:\Users\qlx\.julia\packages\CoolProp\RDEcq\src\CoolProp.jl:81
  [4] myfun_P_MWE(p::ForwardDiff.Dual{ForwardDiff.Tag{…}, Float64, 1}, ::Tuple{Vector{…}, Float64})
    @ Main c:\Users\qlx\Desktop\RB4_AP4_multi_stage_Julia_CoolProp\test2.jl:15
  [5] NonlinearFunction
    @ C:\Users\qlx\.julia\packages\SciMLBase\NjslX\src\scimlfunctions.jl:2185 [inlined]
  [6] JacobianWrapper
    @ C:\Users\qlx\.julia\packages\SciMLBase\NjslX\src\function_wrappers.jl:103 [inlined]
  [7] __value_derivative
    @ C:\Users\qlx\.julia\packages\NonlinearSolve\R8LAH\src\internal\jacobian.jl:186 [inlined]
  [8] JacobianCache
    @ C:\Users\qlx\.julia\packages\NonlinearSolve\R8LAH\src\internal\jacobian.jl:118 [inlined]
  [9] JacobianCache
    @ C:\Users\qlx\.julia\packages\NonlinearSolve\R8LAH\src\internal\jacobian.jl:106 [inlined]
 [10] __step!(cache::NonlinearSolve.GeneralizedFirstOrderAlgorithmCache{…}; recompute_jacobian::Nothing, kwargs::@Kwargs{})
    @ NonlinearSolve C:\Users\qlx\.julia\packages\NonlinearSolve\R8LAH\src\core\generalized_first_order.jl:208
 [11] __step!
    @ C:\Users\qlx\.julia\packages\NonlinearSolve\R8LAH\src\core\generalized_first_order.jl:204 [inlined]
 [12] #step!#210
    @ C:\Users\qlx\.julia\packages\NonlinearSolve\R8LAH\src\core\generic.jl:55 [inlined]
 [13] step!
    @ C:\Users\qlx\.julia\packages\NonlinearSolve\R8LAH\src\core\generic.jl:50 [inlined]
 [14] solve!(cache::NonlinearSolve.GeneralizedFirstOrderAlgorithmCache{…})
    @ NonlinearSolve C:\Users\qlx\.julia\packages\NonlinearSolve\R8LAH\src\core\generic.jl:13
 [15] #__solve#209
    @ C:\Users\qlx\.julia\packages\NonlinearSolve\R8LAH\src\core\generic.jl:4 [inlined]
 [16] __solve
    @ C:\Users\qlx\.julia\packages\NonlinearSolve\R8LAH\src\core\generic.jl:1 [inlined]
 [17] macro expansion
    @ C:\Users\qlx\.julia\packages\NonlinearSolve\R8LAH\src\default.jl:278 [inlined]
 [18] __solve(::NonlinearProblem{…}, ::NonlinearSolvePolyAlgorithm{…}; alias_u0::Bool, verbose::Bool, kwargs::@Kwargs{})
    @ NonlinearSolve C:\Users\qlx\.julia\packages\NonlinearSolve\R8LAH\src\default.jl:248
 [19] __solve
    @ C:\Users\qlx\.julia\packages\NonlinearSolve\R8LAH\src\default.jl:248 [inlined]
 [20] #solve_call#44
    @ C:\Users\qlx\.julia\packages\DiffEqBase\O8cUq\src\solve.jl:612 [inlined]
 [21] solve_call
    @ C:\Users\qlx\.julia\packages\DiffEqBase\O8cUq\src\solve.jl:569 [inlined]
 [22] #solve_up#53
    @ C:\Users\qlx\.julia\packages\DiffEqBase\O8cUq\src\solve.jl:1072 [inlined]
 [23] solve_up
    @ C:\Users\qlx\.julia\packages\DiffEqBase\O8cUq\src\solve.jl:1066 [inlined]
 [24] #solve#52
    @ C:\Users\qlx\.julia\packages\DiffEqBase\O8cUq\src\solve.jl:1060 [inlined]
 [25] solve(prob::NonlinearProblem{…}, args::NonlinearSolvePolyAlgorithm{…})
    @ DiffEqBase C:\Users\qlx\.julia\packages\DiffEqBase\O8cUq\src\solve.jl:1050
 [26] top-level scope
    @ c:\Users\qlx\Desktop\RB4_AP4_multi_stage_Julia_CoolProp\test2.jl:33
Some type information was truncated. Use `show(err)` to see complete types.

My questions are:

  1. Should the PropsSI() function in the CoolProp package be modified to accept a dual number as an argument?

  2. Why does the AutoFiniteDiff() still yield the error about dual number? I expect this error only happens when automatic differentiation is used.

I would greatly appreciate any guidance or workarounds to successfully solve this equation.

Here is my Julia and package version info:

Julia Version 1.10.2
Commit bd47eca2c8 (2024-03-01 10:14 UTC)
Build Info:
  Official https://julialang.org/ release
Platform Info:
  OS: Windows (x86_64-w64-mingw32)
  CPU: 12 × Intel(R) Core(TM) i5-10400F CPU @ 2.90GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-15.0.7 (ORCJIT, skylake)
Threads: 12 default, 0 interactive, 6 GC (on 12 virtual cores)
Environment:
  JULIA_EDITOR = code

  [6e4b80f9] BenchmarkTools v1.5.0
  [e084ae63] CoolProp v0.2.0
  [0c46a032] DifferentialEquations v7.13.0
  [6a86dc24] FiniteDiff v2.23.0
  [f6369f11] ForwardDiff v0.10.36
  [86223c79] Graphs v1.10.0
  [10e44e05] MATLAB v0.8.4
⌃ [961ee093] ModelingToolkit v9.9.0
  [8913a72c] NonlinearSolve v3.9.1
  [7f7a1694] Optimization v3.24.3
  [4e6fcdb7] OptimizationNLopt v0.2.0
⌃ [91a5bcdd] Plots v1.40.3
Info Packages marked with ⌃ have new versions available and may be upgradable.
2

This looks like a bug. Report it to NonlinearSolve.jl

3

Thanks for your reply. I have reported a github issue to NonlinearSolve.jl. By the way, I found a similar issue TrustRegion(radius_update_scheme=RadiusUpdateSchemes.Bastin) doesn’t respect autodiff setting being reported, but using TrustRegion() still does not resolve my problem.

1 Like