Using PH functions with modelingToolkit.jl

[Faisal-Shaikh20]

Hello,
A very impressive package.
I've been trying to use the Clapeyron.PH.temperature function. The function works on its own, but produces an error when used inside a modelingToolkit non=linear problem.
Here is the code to reproduce. It was orignally a model of heat exchanger, but I've tried to cut it to the minimum to reproduce the error.

using ModelingToolkit
using NonlinearSolve
using Clapeyron

const model2 = PR("air")
T_ph(p,h,z) = Clapeyron.PH.temperature(model2,p,h,z,phase=:unknown)
@register_symbolic T_ph(p,h,z::AbstractArray)

T_ph(1e5,6200,[1]) # Works fine

z = [1]
p = 1e5
T = 350
k = 10
h1 = 6200
T1 = T_ph(p, h1, z)
vars = @variables T2 h2 Q

eqs = [h2 ~ h1 + Q,
   T2 ~ T_ph(p, h2, z),
   Q ~ (T - (T1+T2)/2 ) * k,
   ]
sys_ = NonlinearSystem(eqs,vars,[];name=:test)
sys = structural_simplify(sys_)
u0 = [Q => 500]
prob = NonlinearProblem(sys,u0)
sol = solve(prob)

And the error I get:

ERROR: MethodError: no method matching Float64(::ForwardDiff.Dual{ForwardDiff.Tag{DifferentiationInterface.FixTail{…}, Float64}, Float64, 1})
The type `Float64` exists, but no method is defined for this combination of argument types when trying to construct it.

Closest candidates are:
  (::Type{T})(::Real, ::RoundingMode) where T<:AbstractFloat
   @ Base rounding.jl:265
  (::Type{T})(::T) where T<:Number
   @ Core boot.jl:900
  Float64(::IrrationalConstants.Invπ)
   @ IrrationalConstants C:\Users\faisa\.julia\packages\IrrationalConstants\lWTip\src\macro.jl:131

Stacktrace:
  [1] convert(::Type{Float64}, x::ForwardDiff.Dual{ForwardDiff.Tag{DifferentiationInterface.FixTail{…}, Float64}, Float64, 1})
    @ Base .\number.jl:7
  [2] update_state
    @ C:\Users\faisa\.julia\packages\Roots\E1WQf\src\DerivativeFree\secant.jl:38 [inlined]
  [3] update_state
    @ C:\Users\faisa\.julia\packages\Roots\E1WQf\src\DerivativeFree\secant.jl:29 [inlined]
  [4] solve!(𝐙::Roots.ZeroProblemIterator{…}; verbose::Bool)
    @ Roots C:\Users\faisa\.julia\packages\Roots\E1WQf\src\hybrid.jl:54
  [5] solve!
    @ C:\Users\faisa\.julia\packages\Roots\E1WQf\src\hybrid.jl:30 [inlined]
  [6] #solve#47
    @ C:\Users\faisa\.julia\packages\Roots\E1WQf\src\find_zero.jl:492 [inlined]
  [7] solve (repeats 2 times)
    @ C:\Users\faisa\.julia\packages\Roots\E1WQf\src\find_zero.jl:484 [inlined]
  [8] __Tproperty(model::PR{…}, p::Float64, prop::ForwardDiff.Dual{…}, z::StaticArraysCore.SVector{…}, property::typeof(enthalpy), rootsolver::Roots.Order0, phase::Symbol, abstol::Float64, reltol::Float64, threaded::Bool, T0::Float64)
    @ Clapeyron C:\Users\faisa\.julia\packages\Clapeyron\2fEkA\src\methods\property_solvers\Tproperty.jl:262
  [9] Tproperty_pure(model::PR{…}, p::Float64, prop::ForwardDiff.Dual{…}, z::StaticArraysCore.SVector{…}, property::typeof(enthalpy), rootsolver::Roots.Order0, phase::Symbol, abstol::Float64, reltol::Float64, verbose::Bool, threaded::Bool, T0::Nothing)
    @ Clapeyron C:\Users\faisa\.julia\packages\Clapeyron\2fEkA\src\methods\property_solvers\Tproperty.jl:221
 [10] _Tproperty(model::PR{…}, p::Float64, prop::ForwardDiff.Dual{…}, z::StaticArraysCore.SVector{…}, property::typeof(enthalpy); rootsolver::Roots.Order0, phase::Symbol, abstol::Float64, reltol::Float64, T0::Nothing, verbose::Bool, threaded::Bool)
    @ Clapeyron C:\Users\faisa\.julia\packages\Clapeyron\2fEkA\src\methods\property_solvers\Tproperty.jl:85
 [11] _Tproperty
    @ C:\Users\faisa\.julia\packages\Clapeyron\2fEkA\src\methods\property_solvers\Tproperty.jl:64 [inlined]
 [12] px_flash_pure(model::PR{…}, p::Float64, x::ForwardDiff.Dual{…}, z::Vector{…}, spec::typeof(enthalpy), T0::Nothing)
    @ Clapeyron C:\Users\faisa\.julia\packages\Clapeyron\2fEkA\src\methods\property_solvers\multicomponent\flash\general_flash.jl:798
 [13] ph_flash(model::PR{…}, p::Float64, h::ForwardDiff.Dual{…}, z::Vector{…}, method::GeneralizedXYFlash{…})
    @ Clapeyron C:\Users\faisa\.julia\packages\Clapeyron\2fEkA\src\methods\property_solvers\multicomponent\flash\PH.jl:41
 [14] #ph_flash#508
    @ C:\Users\faisa\.julia\packages\Clapeyron\2fEkA\src\methods\property_solvers\multicomponent\flash\PH.jl:22 [inlined]
 [15] ph_flash
    @ C:\Users\faisa\.julia\packages\Clapeyron\2fEkA\src\methods\property_solvers\multicomponent\flash\PH.jl:20 [inlined]
 [16] PH_property(model::PR{…}, p::Float64, h::ForwardDiff.Dual{…}, z::Vector{…}, f::typeof(Clapeyron.temperature), phase::Symbol, T0::Nothing, threaded::Bool)
    @ Clapeyron C:\Users\faisa\.julia\packages\Clapeyron\2fEkA\src\methods\PH.jl:23
 [17] #temperature#4
    @ C:\Users\faisa\.julia\packages\Clapeyron\2fEkA\src\methods\PH.jl:48 [inlined]
 [18] temperature
    @ C:\Users\faisa\.julia\packages\Clapeyron\2fEkA\src\methods\PH.jl:47 [inlined]
 [19] T_ph
    @ .\REPL[45]:1 [inlined]
 [20] macro expansion
    @ C:\Users\faisa\.julia\packages\SymbolicUtils\6fncq\src\code.jl:391 [inlined]
 [21] macro expansion
    @ C:\Users\faisa\.julia\packages\Symbolics\Xv70J\src\build_function.jl:348 [inlined]
 [22] macro expansion
    @ C:\Users\faisa\.julia\packages\RuntimeGeneratedFunctions\M9ZX8\src\RuntimeGeneratedFunctions.jl:163 [inlined]
 [23] macro expansion
    @ .\none:0 [inlined]
 [24] generated_callfunc
    @ .\none:0 [inlined]
 [25] (::RuntimeGeneratedFunctions.RuntimeGeneratedFunction{…})(::Vector{…}, ::Vector{…})
    @ RuntimeGeneratedFunctions C:\Users\faisa\.julia\packages\RuntimeGeneratedFunctions\M9ZX8\src\RuntimeGeneratedFunctions.jl:150
 [26] f
    @ C:\Users\faisa\.julia\packages\ModelingToolkit\K8zNC\src\systems\nonlinear\nonlinearsystem.jl:357 [inlined]
 [27] NonlinearFunction
    @ C:\Users\faisa\.julia\packages\SciMLBase\tWwhl\src\scimlfunctions.jl:2469 [inlined]
 [28] FixTail
    @ C:\Users\faisa\.julia\packages\DifferentiationInterface\TFo7D\src\utils\context.jl:7 [inlined]
 [29] vector_mode_dual_eval!
    @ C:\Users\faisa\.julia\packages\ForwardDiff\UBbGT\src\apiutils.jl:31 [inlined]
 [30] vector_mode_jacobian(f!::DifferentiationInterface.FixTail{…}, y::Vector{…}, x::Vector{…}, cfg::ForwardDiff.JacobianConfig{…})
    @ ForwardDiff C:\Users\faisa\.julia\packages\ForwardDiff\UBbGT\src\jacobian.jl:139
 [31] jacobian
    @ C:\Users\faisa\.julia\packages\ForwardDiff\UBbGT\src\jacobian.jl:40 [inlined]
 [32] jacobian(f!::NonlinearFunction{…}, y::Vector{…}, prep::DifferentiationInterfaceForwardDiffExt.ForwardDiffTwoArgJacobianPrep{…}, backend::AutoForwardDiff{…}, x::Vector{…}, contexts::DifferentiationInterface.Constant{…})
    @ DifferentiationInterfaceForwardDiffExt C:\Users\faisa\.julia\packages\DifferentiationInterface\TFo7D\ext\DifferentiationInterfaceForwardDiffExt\twoarg.jl:433
 [33] construct_jacobian_cache(prob::NonlinearProblem{…}, alg::GeneralizedFirstOrderAlgorithm{…}, f::NonlinearFunction{…}, fu::Vector{…}, u::Vector{…}, p::MTKParameters{…}; stats::SciMLBase.NLStats, autodiff::AutoForwardDiff{…}, vjp_autodiff::AutoFiniteDiff{…}, jvp_autodiff::AutoForwardDiff{…}, linsolve::Nothing)
    @ NonlinearSolveBase C:\Users\faisa\.julia\packages\NonlinearSolveBase\Kek5u\src\jacobian.jl:76
 [34] construct_jacobian_cache
    @ C:\Users\faisa\.julia\packages\NonlinearSolveBase\Kek5u\src\jacobian.jl:34 [inlined]
 [35] __init(::NonlinearProblem{…}, ::GeneralizedFirstOrderAlgorithm{…}; stats::SciMLBase.NLStats, alias_u0::Bool, maxiters::Int64, abstol::Nothing, reltol::Nothing, maxtime::Nothing, termination_condition::Nothing, internalnorm::Function, linsolve_kwargs::@NamedTuple{}, initializealg::NonlinearSolveBase.NonlinearSolveDefaultInit, kwargs::@Kwargs{…})
    @ NonlinearSolveFirstOrder C:\Users\faisa\.julia\packages\NonlinearSolveFirstOrder\3kzAL\src\solve.jl:161
 [36] __init
    @ C:\Users\faisa\.julia\packages\NonlinearSolveFirstOrder\3kzAL\src\solve.jl:121 [inlined]
 [37] __solve(::NonlinearProblem{…}, ::GeneralizedFirstOrderAlgorithm{…}; kwargs::@Kwargs{…})
    @ NonlinearSolveBase C:\Users\faisa\.julia\packages\NonlinearSolveBase\Kek5u\src\solve.jl:5
 [38] __solve
    @ C:\Users\faisa\.julia\packages\NonlinearSolveBase\Kek5u\src\solve.jl:1 [inlined]
 [39] macro expansion
    @ C:\Users\faisa\.julia\packages\NonlinearSolveBase\Kek5u\src\solve.jl:173 [inlined]
 [40] __generated_polysolve(::NonlinearProblem{…}, ::NonlinearSolvePolyAlgorithm{…}; stats::SciMLBase.NLStats, alias_u0::Bool, verbose::Bool, initializealg::NonlinearSolveBase.NonlinearSolveDefaultInit, kwargs::@Kwargs{…})
    @ NonlinearSolveBase C:\Users\faisa\.julia\packages\NonlinearSolveBase\Kek5u\src\solve.jl:130
 [41] __generated_polysolve
    @ C:\Users\faisa\.julia\packages\NonlinearSolveBase\Kek5u\src\solve.jl:130 [inlined]
 [42] #__solve#154
    @ C:\Users\faisa\.julia\packages\NonlinearSolveBase\Kek5u\src\solve.jl:127 [inlined]
 [43] __solve
    @ C:\Users\faisa\.julia\packages\NonlinearSolveBase\Kek5u\src\solve.jl:124 [inlined]
 [44] #__solve#18
    @ C:\Users\faisa\.julia\packages\NonlinearSolve\IChU2\src\default.jl:27 [inlined]
 [45] __solve
    @ C:\Users\faisa\.julia\packages\NonlinearSolve\IChU2\src\default.jl:24 [inlined]
 [46] #__solve#72
    @ C:\Users\faisa\.julia\packages\DiffEqBase\R2Vjs\src\solve.jl:1438 [inlined]
 [47] __solve
    @ C:\Users\faisa\.julia\packages\DiffEqBase\R2Vjs\src\solve.jl:1429 [inlined]
 [48] #solve_call#44
    @ C:\Users\faisa\.julia\packages\DiffEqBase\R2Vjs\src\solve.jl:634 [inlined]
 [49] solve_call
    @ C:\Users\faisa\.julia\packages\DiffEqBase\R2Vjs\src\solve.jl:591 [inlined]
 [50] #solve_up#53
    @ C:\Users\faisa\.julia\packages\DiffEqBase\R2Vjs\src\solve.jl:1107 [inlined]
 [51] solve_up
    @ C:\Users\faisa\.julia\packages\DiffEqBase\R2Vjs\src\solve.jl:1101 [inlined]
 [52] #solve#52
    @ C:\Users\faisa\.julia\packages\DiffEqBase\R2Vjs\src\solve.jl:1095 [inlined]
 [53] solve(::NonlinearProblem{…})
    @ DiffEqBase C:\Users\faisa\.julia\packages\DiffEqBase\R2Vjs\src\solve.jl:1085
Some type information was truncated. Use `show(err)` to see complete types.

I didn't get any similar errors using the PT functions, but in the flowsheet models I am making it seems more natural to keep track of p and h, and calculate other properties when needed.

[longemen3000]

Hello,

This probably requires symbolically registering some functions (Symbolics.@register_symbolic). You can try registering that specific function for now (Clapeyron.PH.temperature)

[Faisal-Shaikh20]

Hi, thanks for replying so quickly. I had registered the function, it doesn't seem to be that. If I don't register the function I get a more understandable error message about using Num in the wrong context.
This seems to be something about ForwardDiff interacting with the Roots.jl package. So maybe it is an issue for them. But I don't know how I can produce a reasonable mwe since this only seems to occur when called within NonLinearSolve.
I can probably work around it for now by just sticking to using the PT functions.

[longemen3000]

The support for using AD on volume was recently added, so if tou find anything, please let me know.

Coming back to the original problem, there is probably a propagation problem on our part that can be fixed. Ideally, the same thing that we did with volume can be done here (do the iterative calculations with the primal values and then calculate the derivatives using the primal result).

Im probably taking this problem to test if AD can be propagated though the stack of functions and iterative procedures that is PH.temperature